Flightradar24 - бул кантип иштейт? 2-бөлүк, ADS-B протоколу

Салам Хабр. Туугандарын же досторун учакта жолуктурган же узаткандардын баары бекер Flightradar24 кызматын колдонушкандыр. Бул реалдуу убакытта учактын абалын байкоо үчүн абдан ыңгайлуу жолу болуп саналат.

В биринчи бөлүгү Мындай онлайн кызматтын иштөө принциби сүрөттөлгөн. Эми биз учактан кабыл алуучу станцияга кандай маалыматтар жөнөтүлүп жана кабыл алынып жатканын аныктап, Python аркылуу өзүбүз чечебиз.

История

Албетте, учак маалыматтары колдонуучулардын смартфондорунда көрүшү үчүн берилбейт. Система ADS-B (Automatic dependent surveillance — уктуруу) деп аталып, учак тууралуу маалыматты башкаруу борборуна автоматтык түрдө берүү үчүн колдонулат – анын идентификатору, координаттары, багыты, ылдамдыгы, бийиктиги жана башка маалыматтар берилет. Мурда мындай системалар пайда болгонго чейин диспетчер радардан бир гана чекит көрө алчу. Самолёттор өтө көп болгон учурда бул жетишсиз болгон.

Техникалык жактан алганда, ADS-B 1090 МГц жетишерлик жогорку жыштыктагы маалымат пакеттерин мезгил-мезгили менен жөнөтүүчү учактагы өткөргүчтөн турат (башка режимдер бар, бирок биз аларды анча кызыктырбайбыз, анткени координаттар бул жерде гана берилет). Албетте, передатчиктен тышкары аэропорттун кайсы бир жеринде кабыл алгыч да бар, бирок биз үчүн, колдонуучулар катары, өзүбүздүн кабыл алгычыбыз кызыктуу.

Баса, салыштыруу үчүн, жөнөкөй колдонуучулар үчүн иштелип чыккан Airnav Radarbox биринчи системасы 2007-жылы пайда болгон жана болжол менен 900 долларды түзгөн, тармактык кызматтарга жазылуу жылына дагы 250 долларды түзөт.

Ошол биринчи орус ээлеринин сын-пикирлери форумда окуса болот радиосканер. Эми RTL-SDR ресиверлери кеңири жеткиликтүү болуп калгандыктан, окшош аппаратты 30 долларга чогултса болот; бул тууралуу көбүрөөк маалымат биринчи бөлүгү. Келгиле, протоколдун өзүнө өтөлү – анын кантип иштээрин карап көрөлү.

Сигналдарды кабыл алуу

Биринчиден, сигнал жазуу керек. Бүтүндөй сигналдын узактыгы 120 микросекундду түзөт, андыктан анын компоненттерин ыңгайлуу демонтаждоо үчүн, кеминде 5 МГц үлгү алуу жыштыгы менен SDR кабылдагыч керек.

Жаздыргандан кийин биз 5000000 30 500 үлгү/сек үлгү алуу ылдамдыгы менен WAV файлын алабыз; Мындай жазуунун XNUMX секунду болжол менен XNUMX МБ "салмагы". Аны медиа ойноткуч менен угуу, албетте, пайдасыз - файлда үн жок, бирок түздөн-түз санариптештирилген радио сигналы - Software Defined Radio так ушундай иштейт.

Биз файлды ачып, Python аркылуу иштетебиз. Өз алдынча эксперимент жасоону каалагандар мисал жазууну жүктөп алса болот байланыш.

Келгиле, файлды жүктөп алып, ичинде эмне бар экенин көрөлү.

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()

Натыйжа: фон ызы-чуусуна каршы ачык "импульстарды" көрөбүз.

Ар бир "импульс" - бул сигнал, анын түзүмү, эгерде сиз графикте резолюцияны жогорулатсаңыз, ачык көрүнүп турат.

Көрүнүп тургандай, сүрөт жогорудагы сыпаттамада айтылгандарга дал келет. Сиз маалыматтарды иштеп баштаса болот.

Декоддоо

Биринчиден, сиз бир аз агым алуу керек. Сигнал өзү Манчестер коддоосу менен коддолгон:

Нибблдердин деңгээли айырмасынан чыныгы "0" жана "1" алуу оңой.

    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 (Downlink Format, 5 бит) - билдирүүнүн түрүн аныктайт. бир нече түрлөрү бар:

(стол булагы)

Бизди DF17 түрү гана кызыктырат, анткени... Бул учактын координаттарын камтыйт.

ICAO (24 бит) - учактын эл аралык уникалдуу коду. Сиз анын коду боюнча учакты текшере аласыз сайтта (тилекке каршы, автор маалымат базасын жаңылоону токтотту, бирок ал актуалдуу бойдон калууда). Мисалы, 3c5ee2 коду үчүн бизде төмөнкү маалыматтар бар:

Түзөтүү: in макалага комментарийлер ИКАО кодунун сүрөттөлүшү кененирээк берилген, кызыккандарга аны окуп чыгууну сунуштайм.

DATA (56 же 112 бит) - биз чечмелей турган чыныгы маалыматтар. Маалыматтын биринчи 5 бит талаа болуп саналат Кодду териңиз, сакталып жаткан маалыматтардын түрүн камтыган (DF менен чаташтырбоо керек). Бул түрлөрүнүн бир нечеси бар:

(стол булагы)

Келгиле, пакеттердин бир нече мисалдарын карап көрөлү.

Учактын идентификациясы

Бинардык формадагы мисал:

00100 011 000101 010111 000111 110111 110001 111000

Маалымат талаалары:

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

TC = 00100b = 4, ар бир символ C1-C8 саптагы индекстерге тиешелүү коддорду камтыйт:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######

Сапты чечмелөө менен учактын кодун алуу оңой: 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('#', ''))

Аба позициясы

Эгер аты жөнөкөй болсо, координаттар татаалыраак. Алар 2, жуп жана так кадр түрүндө берилет. Талаа коду TC = 01011b = 11.

Жуп жана так пакеттердин мисалы:

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

Координаталарды эсептөө абдан татаал формула боюнча ишке ашат:

(булак)

Мен ГИС боюнча адис эмесмин, андыктан анын кайдан келгенин билбейм. Ким билет комментарийге жазгыла.

Бийиктиги жөнөкөй болуп эсептелет - белгилүү бир бит жараша, ал 25 же 100 фут эсе катары көрсөтүлүшү мүмкүн.

Абадагы ылдамдык

TC = 19 менен пакет. Бул жерде кызыктуу нерсе, ылдамдык так, жерге салыштырмалуу болушу мүмкүн (Ground Speed), же абада, учак сенсору (Airspeed) менен өлчөнөт. Көптөгөн ар кандай талаалар да өткөрүлүп берилет:

(булак)

жыйынтыктоо

Көрүнүп тургандай, ADS-B технологиясы кызыктуу симбиоз болуп калды, качан стандарт адистер үчүн гана эмес, жөнөкөй колдонуучулар үчүн да пайдалуу. Бирок, албетте, бул жерде негизги ролду санариптик SDR кабыл алгычтарынын арзаныраак технологиясы ойногон, ал аппаратка гигагерцтен жогору жыштыктагы сигналдарды "тийинчерээк" кабыл алууга мүмкүндүк берет.

Стандарттын өзүндө, албетте, дагы көп нерсе бар. Каалоочулар PDF баракчасынан көрө алышат ICAO же жогоруда айтылган бирине барыңыз сайты.

Жогоруда айтылгандардын баары көптөр үчүн пайдалуу болушу күмөн, бирок, жок эле дегенде, анын кандайча иштээри жөнүндө жалпы түшүнүк, мен үмүттөнөм.

Айтмакчы, Pythonдо даяр декодер бар, сиз аны изилдей аласыз бул жерде. Ал эми SDR кабыл алгычтарынын ээлери даяр ADS-B декодерин чогултуп, ишке киргизе алышат бетинен, бул тууралуу кененирээк сөз болду биринчи бөлүгү.

Макалада сүрөттөлгөн талдоочунун баштапкы коду кесиптин астында берилген. Бул өндүрүш болуп көрүнбөгөн сыноо мисалы, бирок анда кээ бир нерселер иштейт жана аны жогоруда жазылган файлды талдоо үчүн колдонсо болот.
Булак коду (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()

Кимдир бирөө кызыктырды деп үмүттөнөм, көңүлүңүзгө рахмат.

Source: www.habr.com