Flightradar24 - bagaimana ia berfungsi? Bahagian 2, protokol ADS-B

Hello Habr. Mungkin semua orang yang pernah bertemu atau berjumpa saudara atau rakan di dalam kapal terbang telah menggunakan perkhidmatan Flightradar24 percuma. Ini adalah cara yang sangat mudah untuk mengesan kedudukan pesawat dalam masa nyata.

В bahagian pertama Prinsip operasi perkhidmatan dalam talian sedemikian telah diterangkan. Kami kini akan meneruskan dan memikirkan data yang dihantar dan diterima daripada pesawat ke stesen penerima dan menyahkodnya sendiri menggunakan Python.

Kisah

Jelas sekali, data pesawat tidak dihantar untuk dilihat oleh pengguna pada telefon pintar mereka. Sistem ini dipanggil ADS-B (Penyiaran bergantung automatik—siaran), dan digunakan untuk menghantar maklumat secara automatik tentang pesawat ke pusat kawalan - pengecam, koordinat, arah, kelajuan, ketinggian dan data lain dihantar. Sebelum ini, sebelum kemunculan sistem sedemikian, penghantar hanya dapat melihat satu titik pada radar. Ini tidak lagi mencukupi apabila terdapat terlalu banyak pesawat.

Secara teknikal, ADS-B terdiri daripada pemancar pada pesawat yang secara berkala menghantar paket maklumat pada frekuensi yang agak tinggi 1090 MHz (terdapat mod lain, tetapi kami tidak begitu berminat dengannya, kerana koordinat dihantar hanya di sini). Sudah tentu, sebagai tambahan kepada pemancar, terdapat juga penerima di suatu tempat di lapangan terbang, tetapi bagi kami, sebagai pengguna, penerima kami sendiri adalah menarik.

Ngomong-ngomong, sebagai perbandingan, sistem pertama seperti itu, Airnav Radarbox, direka untuk pengguna biasa, muncul pada tahun 2007, dan berharga kira-kira $900; langganan perkhidmatan rangkaian berharga $250 setahun lagi.

Ulasan pemilik Rusia pertama tersebut boleh dibaca di forum pengimbas radio. Kini setelah penerima RTL-SDR telah tersedia secara meluas, peranti serupa boleh dipasang dengan harga $30; lebih lanjut mengenai perkara ini ada dalam bahagian pertama. Mari kita beralih kepada protokol itu sendiri - mari lihat cara ia berfungsi.

Menerima isyarat

Pertama, isyarat perlu dirakam. Keseluruhan isyarat mempunyai tempoh hanya 120 mikrosaat, jadi untuk membuka komponennya dengan selesa, penerima SDR dengan frekuensi pensampelan sekurang-kurangnya 5 MHz adalah wajar.

Selepas rakaman, kami menerima fail WAV dengan kadar pensampelan 5000000 sampel/saat; 30 saat rakaman sedemikian "berberat" kira-kira 500MB. Mendengarnya dengan pemain media, tentu saja, tidak berguna - fail itu tidak mengandungi bunyi, tetapi isyarat radio yang didigitalkan secara langsung - ini adalah cara Radio Defined Perisian berfungsi.

Kami akan membuka dan memproses fail menggunakan Python. Mereka yang ingin mencuba sendiri boleh memuat turun contoh rakaman по ссылке.

Mari muat turun fail 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()

Keputusan: kita melihat "nadi" yang jelas terhadap bunyi latar belakang.

Setiap "nadi" ialah isyarat, strukturnya jelas kelihatan jika anda meningkatkan resolusi pada graf.

Seperti yang anda lihat, gambar itu agak konsisten dengan apa yang diberikan dalam penerangan di atas. Anda boleh mula memproses data.

Penyahkodan

Pertama, anda perlu mendapatkan sedikit aliran. Isyarat itu sendiri dikodkan menggunakan pengekodan Manchester:

Daripada perbezaan tahap dalam nibbles, adalah mudah untuk mendapatkan "0" dan "1" sebenar.

    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 isyarat itu sendiri adalah seperti berikut:

Mari lihat medan dengan lebih terperinci.

DF (Format pautan bawah, 5 bit) - menentukan jenis mesej. Terdapat beberapa jenis:

(sumber jadual)

Kami hanya berminat dengan jenis DF17, kerana... Ini adalah yang mengandungi koordinat pesawat.

ICAO (24 bit) - kod unik antarabangsa pesawat. Anda boleh menyemak pesawat dengan kodnya talian (malangnya, penulis telah berhenti mengemas kini pangkalan data, tetapi ia masih relevan). Sebagai contoh, untuk kod 3c5ee2 kami mempunyai maklumat berikut:

Sunting: dalam komen kepada artikel tersebut Penerangan tentang kod ICAO diberikan dengan lebih terperinci; Saya mengesyorkan agar mereka yang berminat membacanya.

DATA (56 atau 112 bit) - data sebenar yang akan kami dekod. 5 bit pertama data ialah medan Taip Kod, mengandungi subjenis data yang disimpan (jangan dikelirukan dengan DF). Terdapat beberapa jenis ini:

(sumber jadual)

Mari lihat beberapa contoh pakej.

Pengenalan pesawat

Contoh dalam bentuk binari:

00100 011 000101 010111 000111 110111 110001 111000

Medan 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 aksara C1-C8 mengandungi kod yang sepadan dengan indeks dalam baris:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######

Dengan menyahkod rentetan, mudah untuk mendapatkan kod 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('#', ''))

Kedudukan bawaan udara

Jika namanya mudah, maka koordinatnya lebih rumit. Mereka dihantar dalam bentuk bingkai 2, genap dan ganjil. Kod medan TC = 01011b = 11.

Contoh paket genap dan ganjil:

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

Pengiraan koordinat itu sendiri berlaku mengikut formula yang agak rumit:

(sumber)

Saya bukan pakar GIS, jadi saya tidak tahu dari mana asalnya. Siapa tahu, tulis dalam komen.

Ketinggian dianggap lebih mudah - bergantung pada bit tertentu, ia boleh diwakili sama ada gandaan 25 atau 100 kaki.

Halaju Bawaan Udara

Pakej dengan TC=19. Perkara yang menarik di sini ialah kelajuan boleh sama ada tepat, relatif kepada tanah (Ground Speed), atau bawaan udara, diukur oleh sensor pesawat (Airspeed). Banyak bidang yang berbeza juga dihantar:

(sumber)

Kesimpulan

Seperti yang anda lihat, teknologi ADS-B telah menjadi simbiosis yang menarik, apabila standard berguna bukan sahaja kepada profesional, tetapi juga kepada pengguna biasa. Tetapi sudah tentu, peranan penting dalam hal ini dimainkan oleh teknologi penerima SDR digital yang lebih murah, yang membolehkan peranti menerima isyarat secara literal dengan frekuensi melebihi gigahertz "untuk sen."

Dalam standard itu sendiri, sudah tentu, terdapat banyak lagi. Mereka yang berminat boleh melihat PDF di halaman ICAO atau lawati yang telah disebutkan di atas laman web.

Tidak mungkin semua perkara di atas berguna kepada ramai, tetapi sekurang-kurangnya idea umum tentang cara ia berfungsi, saya harap, kekal.

By the way, penyahkod siap pakai dalam Python sudah wujud, anda boleh mengkajinya di sini. Dan pemilik penerima SDR boleh memasang dan melancarkan penyahkod ADS-B siap pakai daripada halaman, ini telah dibincangkan dengan lebih terperinci dalam bahagian pertama.

Kod sumber penghurai yang diterangkan dalam artikel diberikan di bawah potongan. Ini ialah contoh ujian yang tidak berpura-pura sebagai pengeluaran, tetapi beberapa perkara berfungsi di dalamnya, dan ia boleh digunakan untuk menghuraikan fail yang direkodkan di atas.
Kod 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: %f s" % (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()

Saya harap ada yang berminat, terima kasih atas perhatian anda.

Sumber: www.habr.com