Flightradar24 - nola funtzionatzen du? 2. zatia, ADS-B protokoloa

Kaixo Habr. Seguruenik, hegazkinean senideak edo lagunak ezagutu edo ikusi dituzten guztiek doako Flightradar24 zerbitzua erabili dute. Hau oso modu erosoa da hegazkinaren posizioa denbora errealean jarraitzeko.

В Lehenengo zatia Lineako zerbitzu horren funtzionamendu-printzipioa deskribatu zen. Aurrera egingo dugu orain hegazkinetik hartzailera zer datu bidaltzen eta jasotzen ari diren eta guk geuk deskodetuko ditugu Python erabiliz.

Story

Jakina, hegazkinen datuak ez dira transmititzen erabiltzaileek beren telefono mugikorretan ikusteko. Sistemari ADS-B (Automatic dependent surveillance—broadcast) deitzen zaio eta hegazkinari buruzko informazioa automatikoki transmititzeko erabiltzen da kontrol zentrora -bere identifikatzailea, koordenatuak, norabidea, abiadura, altitudea eta bestelako datuak transmititzen dira-. Aurretik, halako sistemak agertu baino lehen, bidaltzaileak radarrean puntu bat besterik ez zuen ikusten. Hau jada ez zen nahikoa hegazkin gehiegi zeudenean.

Teknikoki, ADS-B hegazkin bateko transmisore batek osatzen du, aldian-aldian informazio paketeak 1090 MHz-ko maiztasun nahiko altuan bidaltzen dituena (beste modu batzuk badaude, baina ez zaizkigu hain interesatzen, koordenatuak hemen bakarrik transmititzen baitira). Noski, transmisoreaz gain, aireportuan hargailu bat ere badago nonbait, baina guretzat, erabiltzaile gisa, gure hargailu propioa interesgarria da.

Bide batez, konparazio baterako, erabiltzaile arruntentzat diseinatutako Airnav Radarbox sistema hori 2007an agertu zen eta 900 dolar inguru balio zuen; sareko zerbitzuetarako harpidetzak urtean beste 250 dolar balio zuen.

Lehen errusiar jabe horien iritziak foroan irakur daitezke irrati-eskanerra. Orain, RTL-SDR hargailuak eskuragarri egon direnez, antzeko gailu bat munta daiteke 30 $-ren truke; honi buruzko informazio gehiago ikusi zen. Lehenengo zatia. Goazen protokoloari berari, ikus dezagun nola funtzionatzen duen.

Seinaleak jasotzea

Lehenik eta behin, seinalea grabatu behar da. Seinale osoak 120 mikrosegundoko iraupena du, beraz, bere osagaiak eroso desmuntatzeko, gutxienez 5 MHz-ko laginketa-maiztasuna duen SDR hargailu bat komeni da.

Grabatu ondoren, WAV fitxategi bat jasotzen dugu 5000000 lagin/seg-ko laginketa-abiadura duena; grabazio horren 30 segundoek 500 MB inguru "pisatzen dute". Multimedia erreproduzitzaile batekin entzutea, noski, ez du ezertarako balio - fitxategiak ez du soinurik, zuzenean digitalizatutako irrati-seinalea baizik - horrela funtzionatzen du Software Defined Radio.

Fitxategia ireki eta prozesatuko dugu Python erabiliz. Beren kabuz esperimentatu nahi dutenek grabazio adibide bat deskarga dezakete по ссылке.

Deskargatu dezagun fitxategia eta ikus dezagun zer dagoen barruan.

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

Emaitza: atzeko zarataren kontrako “pultsak” nabariak ikusten ditugu.

"Pultsu" bakoitza seinale bat da, eta horren egitura argi ikusten da grafikoaren bereizmena handitzen baduzu.

Ikus dezakezunez, irudia nahiko koherentea da goiko deskribapenean ematen denarekin. Datuak prozesatzen has zaitezke.

Deskodetzea

Lehenik eta behin, korronte pixka bat lortu behar duzu. Seinalea bera Manchester kodeketa erabiliz kodetzen da:

Nibbles maila-diferentziatik, erraza da benetako "0" eta "1" lortzea.

    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"

Seinalearen egitura bera honako hau da:

Ikus ditzagun eremuak zehatzago.

DF (Downlink Format, 5 bit) - mezu mota zehazten du. Hainbat mota daude:

(taularen iturria)

DF17 mota bakarrik interesatzen zaigu, zeren... Hau da hegazkinaren koordenatuak biltzen dituena.

ICAO (24 bit) - hegazkinaren nazioarteko kode bakarra. Hegazkina bere kodearen arabera egiaztatu dezakezu online (zoritxarrez, egileak datu-basea eguneratzeari utzi dio, baina oraindik garrantzitsua da). Adibidez, 3c5ee2 kodearako informazio hau dugu:

Editatu: in artikuluari iruzkinak ICAO kodearen deskribapena zehatzago ematen da; interesatuei irakur dezatela gomendatzen diet.

DATUAK (56 edo 112 bit) - deskodetuko ditugun benetako datuak. Datuen lehenengo 5 bitak eremua dira Mota kodea, gordetzen ari diren datuen azpimota jasotzen duena (ez nahastu DFrekin). Mota hauetako dezente daude:

(taularen iturria)

Ikus ditzagun paketeen adibide batzuk.

Hegazkinaren identifikazioa

Adibidea forma bitarrean:

00100 011 000101 010111 000111 110111 110001 111000

Datu-eremuak:

+------+------+------+------+------+------+------+------+------+------+
| 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 karaktere bakoitzak lerroko indizeei dagozkien kodeak ditu:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######

Katea deskodetuz, erraza da hegazkinaren kodea eskuratzea: 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('#', ''))

Aireko posizioa

Izena sinplea bada, orduan koordenatuak konplikatuagoak dira. 2 fotograma, bikoiti eta bakoitien moduan transmititzen dira. Eremu-kodea TC = 01011b = 11.

Pakete bikoitien eta bakoitien adibidea:

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

Koordenatuen kalkulua bera formula nahiko zailaren arabera gertatzen da:

(iturri)

Ez naiz GIS aditua, beraz, ez dakit nondik datorren. Nork daki, idatzi iruzkinetan.

Altuera sinpleagotzat jotzen da - bit zehatzaren arabera, 25 edo 100 oineko multiplo gisa irudika daiteke.

Aireko abiadura

TC=19 duen paketea. Gauza interesgarria da hemen abiadura zehatza izan daitekeela, lurrearekiko (Ground Speed) edo airekoa, hegazkinaren sentsore batek neurtuta (Airspeed). Eremu ezberdin asko ere transmititzen dira:

(iturri)

Ondorioa

Ikus dezakezunez, ADS-B teknologia sinbiosi interesgarri bihurtu da, estandar bat profesionalentzat ez ezik, erabiltzaile arruntentzat ere erabilgarria denean. Baina, jakina, horretan funtsezko eginkizuna SDR hargailu digitalen teknologia merkeagoak izan zuen, gailuak literalki gigahertz batetik gorako maiztasunak dituzten seinaleak jasotzeko aukera ematen baitu "zentimoengatik".

Estandarrean bertan, noski, askoz gehiago dago. Interesa dutenek orrialdean PDFa ikus dezakete ICAO edo lehen aipatutakoa bisitatu web.

Nekez da aurreko guztia askorentzat erabilgarria izango denik, baina gutxienez nola funtzionatzen duen ideia orokorra geratzen da, espero dut.

Bide batez, Python-en prest egindako deskodetzaile bat badago jadanik, aztertu dezakezu Hemen. Eta SDR hargailuen jabeek prest egindako ADS-B deskodetzaile bat muntatu eta abiarazi dezakete orrialdetik, honetan eztabaidatu zen zehatzago Lehenengo zatia.

Artikuluan deskribatutako analizatzailearen iturburu kodea ebakiaren azpian ematen da. Hau proba-adibide bat da, produkzioa denik itxuratzen ez duena, baina gauza batzuk bertan funtzionatzen dute, eta goian grabatutako fitxategia analizatzeko erabil daiteke.
Iturburu kodea (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()

Espero dut norbait interesatuta egotea, eskerrik asko zure arretagatik.

Iturria: www.habr.com