Flightradar24 — как работи? Част 2, протокол ADS-B

Здравейте, Хабр. Вероятно всеки, който е изпращал или посрещал роднини или приятели на летище, е използвал безплатния сервис Flightradar24. Това е много удобен начин за проследяване на местоположението на самолета в реално време.

В първата част беше описан принципът на работа на такъв онлайн сервис. Сега ще продължим и ще разберем какви данни се предават и приемат от въздухоплавателното средство до приемната станция, и ще ги декодираме сами с помощта на Python.

История

Очевидно е, че данните за самолетите не се предават, за да могат потребителите да ги виждат на смартфоните си. Системата се нарича ADS–B (Automatic dependent surveillance—broadcast) и служи за автоматично предаване на информация за въздухоплавателното средство в диспетчерския център — предават се идентификатор, координати, посока, скорост, височина и други данни. Преди появата на такива системи, диспетчерът можеше да вижда само точка на радара. Това стана недостатъчно, когато самолетите станаха твърде много.

Технически, ADS-B се състои от предавател на въздухоплавателното средство, който периодично изпраща пакети с информация на достатъчно висока честота от 1090 МГц (има и други режими, но те не ни интересуват, тъй като координатите се предават само тук). Разбира се, освен предавателя, има и приемник някъде на летището, но за нас, като потребители, интересен е нашият собствен приемник.

За сравнение, първата такава система, Airnav Radarbox, предназначена за обикновени потребители, се появи през 2007 година и струваше около 900$, а подписката за мрежови услуги беше около 250$ на година.

Отзиви на първите руски собственици можете да прочетете на форума radioscanner. Сега, когато масово се появиха 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 на страницата 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()

Надявам се, че на някого му беше интересно, благодаря за вниманието.

Източник: habr.com