Flightradar24 — kuidas see töötab? Osa 2, ADS-B protokoll

Tere, Habr. Ilmselt on igaühel, kes on kunagi saatnud või saatnud sugulasi või sõpru lennukisse, olnud võimalus kasutada tasuta teenust Flightradar24. See on väga mugav viis lennuki asukoha reaalajas jälgimiseks.

V esimeses osas oli kirjeldatud sellise veebiteenuse tööpõhimõtet. Nüüd liikume edasi ja selgitame välja, milliseid andmeid edastatakse ja võetakse vastu lennuki ja vastuvõtustatsiooni vahel ning dekodeerime need iseseisvalt Pythoniga.

Ajalugu

On ilmne, et lennukite andmeid edastatakse mitte selleks, et kasutajad näeksid neid oma nutitelefoni ekraanil. Süsteem nimega ADS–B (Automatic dependent surveillance—broadcast) teenib automaatset teabe edastamist lennuki kohta dispetšerkeskusesse — edastatakse selle identifikaator, koordinaadid, suund, kiirus, kõrgus ja muud andmed. Varem, enne selliste süsteemide tekkimist, nägi dispetšer ainult punkti radaril. Seda oli liiga vähe, kui lennukite arv suurenes liiga suureks.

Tehniliselt tähendab ADS-B seda, et lennuki peal on edastaja, mis perioodiliselt saadab pakette teavet piisavalt kõrgel sagedusel 1090 MHz (on ka teisi režiime, kuid need ei huvitavad meid nii väga, kuna koordinaate edastatakse ainult siin). Loomulikult on edastaja kõrval ka vastuvõtja kuskil lennujaamas, kuid meie kui kasutajate jaoks on huvitav meie enda vastuvõtja.

Muide, võrdluseks, esimene selline süsteem, Airnav Radarbox, mis oli mõeldud tavalistele kasutajatele, ilmus 2007. aastal ja maksis umbes 900 dollarit, veel umbes 250 dollarit maksis aastane tellimus võrgu teenustele.

Neid varaseid vene omanike arvustusi saab lugeda foorumist radioscanner. Praegu, kui RTL-SDR vastuvõtjad on laialdaselt kergesti kättesaadavad, saab sarnase seadme kokku panna 30 dollari eest, selle kohta oli rohkem teavet esimeses osas. Läheme nüüd protokolli juurde — vaatame, kuidas see töötab.

Signaalide vastuvõtt

Alustuseks on signaal vaja salvestada. Kogu signaal kestab vaid 120 mikrosekundit, seega on soovitav, et SDR-vastuvõtja oleks vähemalt 5 MHz näidussagedusega.

Pärast salvestamist saame WAV-faili, mille skaadi sagedus on 5000000 näidist/sekundis, 30 sekundi pikkune selline salvestus kaalub umbes 500 MB. Mängida seda meedia mängijaga on muidugi mõttetu — fail ei sisalda helisid, vaid otseselt digitaalset raadiolingid — just nii töötab Software Defined Radio.

Avame ja töötleme faili Pythoniga. Need, kes soovivad ise katsetada, võivad alla laadida näidisfaili. lingi kaudu.

Laadime faili ja vaatame, mis sees on.

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

Tulemus: näeme selgeid 'impulsse' taustamüra seas.

Iga 'impuls' ongi signaal, mille struktuur on hästi nähtav, kui suurendada graafiku eraldusvõimet.

Nagu näha, vastab pilt täiesti sellele, mis on eespool kirjeldatud. Võime andmete töötlemisega alustada.

Dekodeerimine

Alustuseks tuleb saada bitivool. Iga signaal on kodeeritud manchesteri kodeerimisega:

Taseme erinevuse põhjal on lihtne saada reaalsed '0' ja '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"

Signaali struktuur on järgmine:

Vaatame väljad lähemalt.

DF (Downlink Format, 5 bitti) — määrab sõnumi tüübi. Neid on mitu tüüpi:

(allika tabel)

Meid huvitab ainult tüüp DF17, kuna just see sisaldab õhusõiduki koordinaate.

ICAO (24 bitti) — rahvusvaheline ainulaadne õhusõiduki kood. Lennuki kontrollimiseks selle koodi järgi saab veebisaidil (kahjuks on autor andmebaasi uuendamise lõpetanud, kuid see on endiselt aktuaalne). Näiteks koodi 3c5ee2 kohta on meil järgmine teave:

Muudatus: artikli kommentaarides on ICAO koodi kirjeldus esitatud põhjalikumalt, huvilistele soovitan tutvuda.

DATA (56 või 112 bitti) — tegelikult andmed, mida me dekodeerime. Esimesed 5 bitti andmetest — väli Type Code, mis sisaldab salvestatud andmete alamtüüpi (ära segad DF-ga). Selliseid tüüpe on üsna palju:

(allika tabel)

Vaatame mitu paketti näiteid.

Lennuki identifitseerimine

Näide binaarses vormis:

00100 011 000101 010111 000111 110111 110001 111000

Andmeväljad:

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

TC = 00100b = 4, iga sümbol C1-C8 sisaldab koode, mis vastavad indeksitele reas:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######

Dekodeerides rida, on lihtne saada lennuki kood: 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("Lennuki tuvastamine:", code_str.replace('#', ''))

Õhus viibimise asukoht

Kui nimetus on lihtne, siis koordinaatidega on keerulisem. Need edastatakse kahes, paarides ja paaritutes raamistikes. Välja koode TC = 01011b = 11.

Näide paarilisest ja paaritust paketist:

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

Koordinaatide arvutamine toimub üsna nutika valemi abil:

(allikas)

Ma ei ole GIS-i spetsialist, nii et kust see saadakse, ei tea. Kes teab, kirjutage kommentaaridesse.

Kõrgus on lihtsam arvutada — sõltuvalt kindlast bitist võib see esindada kas 25 või 100 jalaga.

Õhus viibimise kiirus

Paketid TC=19. Siin on huvitav see, et kiirus võib olla kas täpne, maapinna suhtes (Ground Speed) või õhus, mida mõõdab lennuki andur (Airspeed). Edastatakse ka palju erinevaid välju:

(allikas)

Kokkuvõte

Nagu näha, on ADS-B tehnoloogia saanud huvitavaks sünteesiks, kus mõni standard on kasulik mitte ainult spetsialistidele, vaid ka tavakasutajatele. Loomulikult mängis siin üliolulist rolli digitaalsete SDR-vastuvõtjate tehnoloogia odavnemine, mis võimaldab seadmel sõna otseses mõttes „väikese raha eest” vastu võtta üle gigaheerti signaale.

Ise standardses on loomulikult palju rohkem. Huvi korral võivad soovijad vaadata PDF-i lehelt ICAO või külastada juba eelpool mainitud veebileht.

Tõenäoliselt ei ole paljudele kasu kogu eeltoodust, kuid vähemalt üldine idee sellest, kuidas see töötab, loodetavasti jäi alles.

Muide, valmis dekooder Pythonis juba eksisteerib, seda saab uurida siit. Ja SDR-vastuvõtjate omanikud saavad kokku panna ja käivitada valmis ADS-B dekooderi lehelt, sellest räägiti lähemalt esimeses osas.

Artiklis kirjeldatud parseri lähtekood on toodud allpool. See on testnäide, mis ei pretendeeri tootmisstandardile, kuid mõned asjad selles töötavad, ja salvestatud faili saab sellega parsida.
Lähtekood (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()

Loodan, et kellelegi oli huvitav, aitäh tähelepanu eest.

Allikas: habr.com