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

Tere, Habr. Ilmselt on igaüks, kes on vähemalt korra saatnud või saatnud sugulasi või sõpru lennuki peale, kasutanud tasuta teenust Flightradar24. See on väga mugav viis lennuki asukoha jälgimiseks reaalajas.

Uues esimeses osas oli kirjeldatud sellise veebiteenuse tööpõhimõtet. Nüüd liigume edasi ja uurime, milliseid andmeid edastatakse ja võetakse vastu õhusõiduki ja vastuvõtja jaama vahel ning dekodeerime neid iseseisvalt Pythonis.

Ajalugu

On iseenesestmõistetav, et lennukite andmeid ei edastata selleks, et kasutajad saaksid neid oma nutitelefonides näha. Süsteem nimetatakse ADS–B (Automatic dependent surveillance—broadcast) ja see on mõeldud õhusõiduki teabe automatiseeritud edastamiseks lennujuhtimiskeskusele — edastatakse selle identifikaator, koordinaadid, suund, kiirus, kõrgus ja muud andmed. Enne selliste süsteemide tekkimist nägid lennujuhtimistöötajad radaril ainult punkti. Seda ei olnud piisavalt, kui lennukeid hakkas liiga palju olema.

Tehniliselt koosneb ADS-B õhusõiduki saatjast, mis perioodiliselt saadab teavet sisaldavaid pakette piisavalt kõrgel sagedusel 1090 MHz (on ka teisi režiime, kuid need ei huvitavaid nii palju, kuna koordinaate edastatakse ainult siin). Loomulikult on lisaks saatjale ka vastuvõtja kuskil lennujaamas, kuid meie kasutajatena huvitab meid meie enda vastuvõtja.

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

Esimeste Venemaa omanike tagasisidet saab lugeda foorumist radioscanner. Praegu, kui RTL-SDR vastuvõtjad on laialdaselt kergesti kättesaadavad, saab sarnase seadme kokku panna 30 dollari eest, sellest oli rohkem juttu esimeses osas. Liigume nüüd tegelikult protokolli juurde — vaatame, kuidas see töötab.

Signaalide vastuvõtt

Alustuseks tuleb signaal salvestada. Kogu signaali kestus on vaid 120 mikrosekundit, seetõttu on soovitatav kasutada SDR-vastuvõtjat, mille proovivõtmisfrekvents on vähemalt 5 MHz.

Pärast salvestamist saame WAV-faili, mille näidustussagedus on 5000000 näidist/s, ja 30 sekundi salvestus „kaalub” umbes 500MB. Selle kuulamine meediapleieriga on muidugi mõttetu — fail sisaldab mitte heli, vaid otseselt digitaliseeritud raadiosignaali — just nii töötab Software Defined Radio.

Avame ja töötleme faili Pythoniga. Kes soovib iseseisvalt katsetada, saavad alla laadida salvestuse näite. linki pidi.

Laadime faili ja vaatame, mis selle 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, et taustamüra seas on selged „impulsid”.

Iga „impulss” on signaal, mille struktuuri on hästi näha, kui skaalat joonisel suurendada.

Nagu näha, on pilt täiesti kooskõlas eespool toodud kirjeldustega. Saame alustada andmete töötlemist.

Dekodeerimine

Esialgu peame saama bitivoog. Iga signaal on kodeeritud Manchesteri kodeerimisega:

Poolebaitide tasemeerinevusest on lihtne saada tõelisi „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älju lähemalt.

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

(allika tabel)

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

ICAO (24 bitti) — rahvusvaheline ainulaadne õhusõiduki kood. Saate lennukit tema koodi järgi kontrollida. veebisaidil (kahjuks on autor andmebaasi uuendamise lõpetanud, kuid see on siiski veel aktuaalne). Näiteks koodi 3c5ee2 puhul on meil järgmine teave:

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

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

(allika tabel)

Vaadakem mõningaid pakettide näiteid.

Lennukite tuvastamine

Näide binaarvormingus:

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, ei ole keeruline saada lennuki koodi: 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 identifitseerimine:", code_str.replace('#', ''))

Õhusolek

Kui nimi on lihtne, siis koordinaatidega on keerulisem. Need edastatakse kahes, paarist ja paaritust raamist. Koode väli TC = 01011b = 11.

Näide paaritud ja paaris pakettidest:

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

Koordinaatide arvutamine toimub üsna keerulise valemi järgi:

(allikas)

Ma ei ole GIS-i ekspert, seega ei tea, kust see tuleneb. Kes on kursis, kirjutage kommentaarides.

Kõrgus arvutatakse lihtsamalt — sõltuvalt teatud bitist, võib see olla kas 25 või 100 jalga.

Õhusolekukiirus

Pakett TC=19. Siin on huvitav see, et kiirus võib olla kas täpne, maapinna suhtes (Ground Speed), või õhus, lennuki anduriga mõõdetud (Airspeed). Samuti edastatakse palju erinevaid välju:

(allikas)

Kokkuvõte

Nagu näha, on ADS-B tehnoloogia huvitav sümbioos, kus mingi standard on kasulik mitte ainult professionaalidele, vaid ka tavalistele kasutajatele. Kuid loomulikult mängis võtmerolli tehnoloogia digitaalsete SDR-vastuvõtjate odavnemine, mis võimaldab seadmel sõna otseses mõttes "mugavusest" vastu võtta gigaheertsise sagedusega signaale.

Kindlasti sisaldab standard palju enamat. Soovijad saavad vaadata PDF-faili lehelt ICAO või külastada juba varem mainitud veebisaidile.

Pigem ei ole see, mida enamus vajab, kuid vähemalt lootuses, et põhikontseptsioon sellest, kuidas see töötab, jääb.

Muide, valmis dekooder Pythonis on juba olemas, seda saab uurida siin. Ja SDR-vastuvõtjate omanikud saavad koguda ja käivitada valmis ADS-B dekooderi lehelt, sellest on räägitud esimeses osas.

Artiklis kirjeldatud parsi algkood on allpool. See on katse näide, mis ei pretendeeri produktsioonile, kuid mõningaid asju töötab ja selle abil saab eelnevalt salvestatud faili parsida.
Algkood (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 see huvitav, aitäh tähelepanu eest.

Allikas: habr.com