Geekly articles weekly

Flightradar24 - bagaimana cara kerjanya? Bagian 2, Protokol ADS-B

Hai Habr. Mungkin setiap orang yang pernah bertemu atau menemani kerabat atau teman di pesawat menggunakan layanan Flightradar24 gratis. Ini adalah cara yang sangat nyaman untuk melacak posisi pesawat secara real time.

gambar

Bagian pertama menggambarkan prinsip pengoperasian layanan online semacam itu. Sekarang kita akan melangkah lebih jauh, dan mencari tahu data apa yang dikirimkan dan diterima dari pesawat ke stasiun penerima, dan mendekodekannya sendiri menggunakan Python.

Ceritanya


Jelas, data pesawat tidak ditransmisikan sehingga pengguna dapat melihatnya di smartphone mereka. Sistem ini disebut ADS-B (Automatic-surveillance-broadcast - broadcast), dan digunakan untuk secara otomatis mengirimkan informasi pesawat ke pusat kendali - pengidentifikasi, koordinat, arah, kecepatan, ketinggian, dan data lainnya dikirimkan. Sebelumnya, sebelum munculnya sistem seperti itu, operator hanya bisa melihat titik di radar. Ini tidak cukup ketika ada terlalu banyak pesawat.

Secara teknis, ADS-B terdiri dari pemancar di pesawat terbang yang secara berkala mengirimkan paket dengan informasi pada frekuensi yang cukup tinggi 1090 MHz (ada mode lain, tetapi mereka tidak begitu menarik bagi kami, karena koordinatnya hanya dikirimkan di sini). Tentu saja, selain pemancar, ada juga penerima di suatu tempat di bandara, tetapi bagi kami, seperti bagi pengguna, penerima kami sendiri menarik.

Ngomong-ngomong, untuk perbandingan, sistem seperti itu yang pertama, Airnav Radarbox, dirancang untuk pengguna biasa, muncul pada 2007, dan harganya sekitar $ 900, sekitar $ 250 setahun layak berlangganan layanan jaringan.



Anda dapat membaca ulasan pemilik Rusia pertama di forum radioscanner . Sekarang penerima RTL-SDR telah tersedia secara luas, perangkat serupa dapat dirakit seharga $ 30, lebih lanjut tentang ini di bagian pertama . Kami akan melanjutkan ke protokol itu sendiri - mari kita lihat cara kerjanya.

Menerima sinyal


Untuk memulai, sinyal perlu direkam. Seluruh sinyal memiliki durasi hanya 120 mikrodetik, oleh karena itu, untuk dapat membongkar komponennya dengan nyaman, penerima SDR dengan frekuensi pengambilan sampel minimal 5 MHz diperlukan.

gambar

Setelah merekam, kami mendapatkan file WAV dengan frekuensi sampling 5.000.000 sampel / detik, 30 detik dari rekaman "menimbang" sekitar 500 MB. Mendengarkannya dengan pemutar media tentu saja tidak berguna - file tidak mengandung suara, tetapi sinyal radio yang langsung didigitalkan - ini adalah cara kerja Radio yang Ditentukan Perangkat Lunak.

Kami akan membuka dan memproses file menggunakan Python. Mereka yang ingin bereksperimen sendiri dapat mengunduh rekaman sampel dari tautan .

Unduh file dan lihat apa yang ada di dalamnya.

from scipy.io import wavfile import matplotlib.pyplot as plt import numpy as np fs, data = wavfile.read("adsb_20190311_191728Z_1090000kHz_RF.wav") data = data.astype(float) I, Q = data[:, 0], data[:, 1] A = np.sqrt(I*I + Q*Q) plt.plot(A) plt.show()

Hasil: kita melihat "impuls" yang jelas terhadap latar belakang kebisingan.



Setiap "impuls" adalah sinyal yang strukturnya terlihat jelas jika Anda meningkatkan resolusi pada grafik.



Seperti yang Anda lihat, gambar konsisten dengan apa yang dijelaskan dalam uraian di atas. Anda dapat mulai memproses data.

Decoding


Pertama, Anda perlu mendapatkan bitstream. Sinyal itu sendiri dikodekan menggunakan manchester encoding:



Dari perbedaan level dalam camilan, mudah untuk mendapatkan "0" dan "1" yang nyata.

  bits_str = "" for p in range(8): pos = start_data + bit_len*p p1, p2 = A[pos: pos + bit_len/2], A[pos + bit_len/2: pos + bit_len] avg1, avg2 = np.average(p1), np.average(p2) if avg1 < avg2: bits_str += "0" elif avg1 > avg2: bits_str += "1"

Struktur sinyal itu sendiri adalah sebagai berikut:



Mari kita perhatikan kolom-kolomnya lebih detail.

DF (Downlink Format, 5 bit) - mendefinisikan jenis pesan. Ada beberapa jenis:


( sumber tabel )

Kami hanya tertarik pada tipe DF17, seperti ini berisi koordinat pesawat.

ICAO (24 bit) adalah kode pesawat unik internasional. Anda dapat memeriksa pesawat dengan kodenya di situs web (sayangnya, penulis telah berhenti memperbarui database, tetapi masih relevan). Misalnya, untuk kode 3c5ee2 kami memiliki informasi berikut:



Sunting: dalam komentar pada artikel, deskripsi kode ICAO diberikan secara lebih rinci, saya sarankan Anda membacanya untuk mereka yang tertarik.

DATA (56 atau 112 bit) - sebenarnya data yang akan kami dekode. 5 bit pertama data adalah bidang Kode Jenis yang berisi subtipe dari data yang disimpan (jangan dikacaukan dengan DF). Ada beberapa jenis ini:


( sumber tabel )

Mari kita lihat beberapa paket contoh.

Identifikasi pesawat

Contoh biner:

00100 011 000101 010111 000111 110111 110001 111000

Bidang Data:

 +------+------+------+------+------+------+------+------+------+------+ | TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 | +------+------+------+------+------+------+------+------+------+------+

TC = 00100b = 4, setiap karakter C1-C8 berisi kode yang sesuai dengan indeks dalam string:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ ##### _ ################ 0123456789 ######

Setelah memecahkan kode baris, mudah untuk mendapatkan kode pesawat: EWG7184

 symbols = "#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######" code_str = "" for p in range(8): c = int(bits_str[8 + 6*p:8 + 6*(p + 1)], 2) code_str += symbols[c] print("Aircraft Identification:", code_str.replace('#', ''))

Posisi di udara

Jika namanya sederhana, maka koordinatnya lebih rumit. Mereka ditransmisikan dalam bentuk frame 2x, genap dan ganjil. Kode bidang TC = 01011b = 11.



Contoh paket genap dan ganjil:

 01011 000 000101110110 00 10111000111001000 10000110101111001 01011 000 000110010000 01 10010011110000110 10000011110001000

Penghitungan koordinat dilakukan berdasarkan rumus yang agak pintar:


( sumber )

Saya bukan spesialis GIS, jadi saya tidak tahu dari mana asalnya. Siapa yang tahu, tulis di komentar.

Ketinggian dianggap lebih mudah - tergantung pada bit tertentu, itu dapat muncul sebagai kelipatan 25 atau 100 kaki.

Kecepatan udara

Paket dengan TC = 19. Yang menarik di sini adalah bahwa kecepatannya bisa akurat, relatif terhadap ground (Kecepatan Tanah), atau udara, yang diukur dengan sensor pesawat (Kecepatan Udara). Banyak bidang yang berbeda juga ditransmisikan:


( sumber )

Kesimpulan


Seperti yang Anda lihat, teknologi ADS-B telah menjadi simbiosis yang menarik ketika standar berguna tidak hanya untuk para profesional, tetapi juga untuk pengguna biasa. Tapi tentu saja, peran utama dalam hal ini dimainkan oleh semakin murahnya teknologi penerima SDR digital, yang memungkinkan perangkat untuk menerima sinyal dengan frekuensi yang lebih tinggi daripada gigahertz pada perangkat yang secara harfiah “untuk satu sen”.

Dalam standar itu sendiri, tentu saja, lebih dari segalanya. Mereka yang tertarik dapat melihat PDF di halaman ICAO atau mengunjungi situs yang telah disebutkan di atas.

Rasanya tidak mungkin banyak hal di atas berguna, tetapi setidaknya gagasan umum tentang cara kerjanya, saya harap, tetap ada.

Omong-omong, dekoder siap pakai dengan Python sudah ada, dapat dipelajari di sini . Dan pemilik penerima SDR dapat merakit dan menjalankan dekoder ADS-B yang telah selesai dari halaman , lebih lanjut tentang ini di bagian pertama .

Kode sumber parser yang dijelaskan dalam artikel diberikan di bawah potongan. Ini adalah contoh pengujian yang tidak berpura-pura menjadi produksi, tetapi sesuatu berfungsi di dalamnya, dan Anda dapat mengurai file yang direkam di atas.
Kode Sumber (Python)
 from __future__ import print_function from scipy.io import wavfile from scipy import signal import matplotlib.pyplot as plt import numpy as np import math import sys def parse_message(data, start, bit_len): max_len = bit_len*128 A = data[start:start + max_len] A = signal.resample(A, 10*max_len) bits = np.zeros(10*max_len) bit_len *= 10 start_data = bit_len*8 # Parse first 8 bits bits_str = "" for p in range(8): pos = start_data + bit_len*p p1, p2 = A[pos: pos + bit_len/2], A[pos + bit_len/2: pos + bit_len] avg1, avg2 = np.average(p1), np.average(p2) if avg1 < avg2: bits_str += "0" elif avg1 > avg2: bits_str += "1" df = int(bits_str[0:5], 2) # Aircraft address (db - https://junzis.com/adb/?q=3b1c5c ) bits_str = "" for p in range(8, 32): pos = start_data + bit_len * p p1, p2 = A[pos: pos + bit_len / 2], A[pos + bit_len / 2: pos + bit_len] avg1, avg2 = np.average(p1), np.average(p2) if avg1 < avg2: bits_str += "0" elif avg1 > avg2: bits_str += "1" # print "Aircraft address:", bits_str, hex(int(bits_str, 2)) address = hex(int(bits_str, 2)) # Filter specific aircraft (optional) # if address != "0x3c5ee2": # return if df == 16 or df == 17 or df == 18 or df == 19 or df == 20 or df == 21: # print "Pos:", start, "DF:", msg_type # Data (56bit) bits_str = "" for p in range(32, 88): pos = start_data + bit_len*p p1, p2 = A[pos: pos + bit_len/2], A[pos + bit_len/2: pos + bit_len] avg1, avg2 = np.average(p1), np.average(p2) if avg1 < avg2: bits_str += "0" # bits[pos + bit_len / 2] = 50 elif avg1 > avg2: bits_str += "1" # http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html # print "Data:" # print bits_str[:8], bits_str[8:20], bits_str[20:22], bits_str[22:22+17], bits_str[39:39+17] # Type Code: tc, ec = int(bits_str[:5], 2), int(bits_str[5:8], 2) # print("DF:", df, "TC:", tc) # 1 - 4 Aircraft identification # 5 - 8 Surface position # 9 - 18 Airborne position (w/ Baro Altitude) # 19 Airborne velocities if tc >= 1 and tc <= 4: # and (df == 17 or df == 18): print("Aircraft address:", address) print("Data:") print(bits_str[:8], bits_str[8:14], bits_str[14:20], bits_str[20:26], bits_str[26:32], bits_str[32:38], bits_str[38:44]) symbols = "#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######" code_str = "" for p in range(8): c = int(bits_str[8 + 6*p:8 + 6*(p + 1)], 2) code_str += symbols[c] print("Aircraft Identification:", code_str.replace('#', '')) print() if tc == 11: print("Aircraft address:", address) print("Data: (11)") print(bits_str[:8], bits_str[8:20], bits_str[20:22], bits_str[22:22+17], bits_str[39:39+17]) # Bit 22 contains the F flag which indicates which CPR format is used (odd or even) # First frame has F flag = 0 so is even and the second frame has F flag = 1 so odd # f = bits_str[21:22] # print("F:", int(f, 2)) # Altitude alt1b = bits_str[8:20] if alt1b[-5] == '1': bits = alt1b[:-5] + alt1b[-4:] n = int(bits, 2) alt_ft = n*25 - 1000 print("Alt (ft)", alt_ft) # lat_dec = int(bits_str[22:22+17], 2) # lon_dec = int(bits_str[39:39+17], 2) # print("Lat/Lon:", lat_dec, lon_dec) # http://airmetar.main.jp/radio/ADS-B%20Decoding%20Guide.pdf print() if tc == 19: print("Aircraft address:", address) print("Data:") # print(bits_str) print(bits_str[:5], bits_str[5:8], bits_str[8:10], bits_str[10:13], bits_str[13] ,bits_str[14:24], bits_str[24], bits_str[25:35], bits_str[35:36], bits_str[36:65]) subtype = int(bits_str[5:8], 2) # https://mode-s.org/decode/adsb/airborne-velocity.html spd, hdg, rocd = -1, -1, -1 if subtype == 1 or subtype == 2: print("Velocity Subtype 1: Ground speed") v_ew_sign = int(bits_str[13], 2) v_ew = int(bits_str[14:24], 2) - 1 # east-west velocity v_ns_sign = int(bits_str[24], 2) v_ns = int(bits_str[25:35], 2) - 1 # north-south velocity v_we = -1*v_ew if v_ew_sign else v_ew v_sn = -1*v_ns if v_ns_sign else v_ns spd = math.sqrt(v_sn*v_sn + v_we*v_we) # unit in kts hdg = math.atan2(v_we, v_sn) hdg = math.degrees(hdg) # convert to degrees hdg = hdg if hdg >= 0 else hdg + 360 # no negative val if subtype == 3: print("Subtype Subtype 3: Airspeed") hdg = int(bits_str[14:24], 2)/1024.0*360.0 spd = int(bits_str[25:35], 2) vr_sign = int(bits_str[36], 2) vr = int(bits_str[36:45], 2) rocd = -1*vr if vr_sign else vr # rate of climb/descend print("Speed (kts):", spd, "Rate:", rocd, "Heading:", hdg) print() # print() def calc_coordinates(): def _cprN(lat, is_odd): nl = _cprNL(lat) - is_odd return nl if nl > 1 else 1 def _cprNL(lat): try: nz = 15 a = 1 - math.cos(math.pi / (2 * nz)) b = math.cos(math.pi / 180.0 * abs(lat)) ** 2 nl = 2 * math.pi / (math.acos(1 - a/b)) return int(math.floor(nl)) except: # happens when latitude is +/-90 degree return 1 def floor_(x): return int(math.floor(x)) lat1b, lon1b, alt1b = "10111000111010011", "10000110111111000", "000101111001" lat2b, lon2b, alt2b = "10010011101011100", "10000011000011011", "000101110111" lat1, lon1, alt1 = int(lat1b, 2), int(lon1b, 2), int(alt1b, 2) lat2, lon2, alt2 = int(lat2b, 2), int(lon2b, 2), int(alt2b, 2) # 131072 is 2^17, since CPR lat and lon are 17 bits each cprlat_even, cprlon_even = lat1/131072.0, lon1/131072.0 cprlat_odd, cprlon_odd = lat2/131072.0, lon2/131072.0 print(cprlat_even, cprlon_even) j = floor_(59*cprlat_even - 60*cprlat_odd) print(j) air_d_lat_even = 360.0 / 60 air_d_lat_odd = 360.0 / 59 # Lat lat_even = float(air_d_lat_even * (j % 60 + cprlat_even)) lat_odd = float(air_d_lat_odd * (j % 59 + cprlat_odd)) if lat_even >= 270: lat_even = lat_even - 360 if lat_odd >= 270: lat_odd = lat_odd - 360 # Lon ni = _cprN(lat_even, 0) m = floor_(cprlon_even * (_cprNL(lat_even)-1) - cprlon_odd * _cprNL(lat_even) + 0.5) lon = (360.0 / ni) * (m % ni + cprlon_even) print("Lat", lat_even, "Lon", lon) # Altitude # Q-bit (bit 48) indicates whether the altitude is encoded in multiples of 25 or 100 ft (0: 100 ft, 1: 25 ft) # The value can represent altitudes from -1000 to +50175 ft. if alt1b[-5] == '1': bits = alt1b[:-5] + alt1b[-4:] n = int(bits, 2) alt_ft = n*25 - 1000 print("Alt (ft)", alt_ft) fs, data = wavfile.read("adsb_20190311_191728Z_1090000kHz_RF.wav") T = 1/fs print("Sample rate %f MS/s" % (fs / 1e6)) print("Cnt samples %d" % len(data)) print("Duration: %fs" % (T * len(data))) data = data.astype(float) cnt = data.shape[0] # Processing only part on file (faster): # cnt = 10000000 # data = data[:cnt] print("Processing I/Q...") I, Q = data[:, 0], data[:, 1] A = np.sqrt(I*I + Q*Q) bits = np.zeros(cnt) # To see scope without any processing, uncomment # plt.plot(A) # plt.show() # sys.exit(0) print("Extracting signals...") pos = 0 avg = 200 msg_start = 0 # Find beginning of each signal while pos < cnt - 16*1024: # P1 - message start while pos < cnt - 16*1024: if A[pos] < avg and A[pos+1] > avg and pos - msg_start > 1000: msg_start = pos bits[pos] = 100 pos += 4 break pos += 1 start1, start2, start3, start4 = msg_start, 0, 0, 0 # P2 while pos < cnt - 16*1024: if A[pos] < avg and A[pos+1] > avg: start2 = pos bits[pos] = 90 pos += 1 break pos += 1 # P3 while pos < cnt - 16*1024: if A[pos] < avg and A[pos+1] > avg: start3 = pos bits[pos] = 80 pos += 1 break pos += 1 # P4 while pos < cnt - 16*1024: if A[pos] < avg and A[pos+1] > avg: start4 = pos bits[pos] = 70 pos += 1 break pos += 1 sig_diff = start4 - start1 if 20 < sig_diff < 25: bits[msg_start] = 500 bit_len = int((start4 - start1) / 4.5) # print(pos, start1, start4, ' - ', bit_len) # start = start1 + 8*bit_len parse_message(A, msg_start, bit_len) pos += 450 # For debugging: check signal start # plt.plot(A) # plt.plot(bits) # plt.show()


Semoga ada yang tertarik, terima kasih atas perhatiannya.

Source: https://habr.com/ru/post/id443498/

More articles:


All Articles