Здравейте, Хабр. Вероятно всеки, който е изпращал или посрещал роднини или приятели на летище, е използвал безплатния сервис Flightradar24. Това е много удобен начин за проследяване на местоположението на самолета в реално време.
В беше описан принципът на работа на такъв онлайн сервис. Сега ще продължим и ще разберем какви данни се предават и приемат от въздухоплавателното средство до приемната станция, и ще ги декодираме сами с помощта на Python.
История
Очевидно е, че данните за самолетите не се предават, за да могат потребителите да ги виждат на смартфоните си. Системата се нарича ADS–B (Automatic dependent surveillance—broadcast) и служи за автоматично предаване на информация за въздухоплавателното средство в диспетчерския център — предават се идентификатор, координати, посока, скорост, височина и други данни. Преди появата на такива системи, диспетчерът можеше да вижда само точка на радара. Това стана недостатъчно, когато самолетите станаха твърде много.
Технически, ADS-B се състои от предавател на въздухоплавателното средство, който периодично изпраща пакети с информация на достатъчно висока честота от 1090 МГц (има и други режими, но те не ни интересуват, тъй като координатите се предават само тук). Разбира се, освен предавателя, има и приемник някъде на летището, но за нас, като потребители, интересен е нашият собствен приемник.
За сравнение, първата такава система, Airnav Radarbox, предназначена за обикновени потребители, се появи през 2007 година и струваше около 900$, а подписката за мрежови услуги беше около 250$ на година.
Отзиви на първите руски собственици можете да прочетете на форума . Сега, когато масово се появиха RTL-SDR приемници, аналогично устройство може да се събере за 30$, по-подробно за това беше в . Ние обаче ще преминем директно към протокола — нека видим как работи.
Прием на сигнали
Първо, сигналът трябва да бъде записан. Целият сигнал има продължителност само 120 микросекунди, затова за комфортно разясняване на компонентите му, желателен е SDR приемник с честота на дискретизация не по-малка от 5 MHz.
След записване получаваме WAV файл с честота на дискретизация 5000000 семпла/сек, 30 секунди такова записване «тежи» около 500Мб. Слушането му с медия плейър е разбира се безполезно — файлът съдържа не звук, а директно оцифрован радиосигнал — точно по този начин работи 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()
Резултат: виждаме ясни «импулси» на фона на шума.
Всеки «импулс» е именно сигнал, структурата на който е добре видима, ако увеличим разрешението на графиката.
Както можем да видим, картинката напълно съответства на описаното по-горе. Можем да преминем към обработка на данните.
Декодиране
Първо, трябва да получим битовия поток. Самият сигнал е закодиран с помощта на manchester encoding:
От разликата на нивата в полубайтите лесно можем да получим реални «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 += "1"
Структура на самия сигнал изглежда по следния начин:
Нека разгледаме полетата по-подробно.
DF (Downlink Format, 5 бита) — определя типа на съобщението. Има няколко типа:
()
Интересува ни само тип DF17, тъй като именно той съдържа координатите на въздушното съдно.
ICAO (24 бита) — международен уникален код на въздушното съдно. Можете да проверите самолета по неговия код на (за съжаление, авторът е спрял да актуализира базата, но тя все още е актуална). Например, за код 3c5ee2 имаме следната информация:
Корекция: в описанието на кода ICAO е приведено по-подробно, интересуващите се препоръчвам да се запознаят.
ДАНИ (56 или 112 бита) — самите данни, които ще декодираме. Първите 5 бита от данните — поле Type Code, съдържащо подтип на съхраняваните данни (да не се бърка с 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("Идентификация на самолета:", code_str.replace('#', ''))
Въздушна позиция
Ако с наименованието всичко е просто, то с координатите е малко по-сложно. Те се предават под формата на четни и нечетни кадри. Кодът на полето 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 на страницата или да посетят вече споменатия по-горе .
Малко вероятно е много хора да се нуждаят от всичко написано досега, но поне основната идея за това как работи, надявам се, остана.
Всъщност, готов декодер на 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()
Надявам се, че на някого му беше интересно, благодаря за вниманието.
Източник: habr.com
