Ciao Habr. Probabilmente chiunque abbia accompagnato o salutato parenti o amici all'aeroporto ha utilizzato il servizio gratuito Flightradar24. È un modo molto comodo per monitorare la posizione di un aereo in tempo reale.
In è stato descritto il principio di funzionamento di un tale servizio online. Ora andremo oltre e scopriremo quali dati vengono trasmessi e ricevuti dall'aeromobile alla stazione di ricezione, e li decodificheremo da soli utilizzando Python.
Storia
È ovvio che i dati sugli aerei non vengono trasmessi affinché gli utenti li vedano sui loro smartphone. Il sistema si chiama ADS–B (Automatic dependent surveillance—broadcast) e serve per la trasmissione automatica di informazioni sull'aeromobile al centro di controllo — vengono trasmessi il suo identificatore, coordinate, direzione, velocità, altitudine e altri dati. In precedenza, prima dell'arrivo di tali sistemi, il controllore poteva vedere solo un punto sul radar. Questo è diventato insufficiente quando gli aerei sono aumentati drasticamente.
Tecnicamente, ADS-B è composto da un trasmettitore sull'aeromobile, che invia periodicamente pacchetti con informazioni a una frequenza abbastanza alta di 1090 MHz (ci sono altri modalità, ma non sono di nostro interesse, poiché le coordinate vengono trasmesse solo qui). Naturalmente, oltre al trasmettitore, c'è un ricevitore da qualche parte in aeroporto, ma per noi, come utenti, è interessante il nostro ricevitore personale.
A proposito, a titolo di confronto, il primo di questi sistemi, Airnav Radarbox, pensato per utenti normali, è apparso nel 2007 e costava circa 900$, con un abbonamento ai servizi di rete che costava ulteriori 250$ all'anno.
Le recensioni dei primi proprietari russi possono essere lette nel forum . Ora, quando i ricevitori RTL-SDR sono diventati ampiamente disponibili, un dispositivo simile può essere assemblato per 30$, di più su questo si è parlato in . Passiamo quindi al protocollo — vediamo come funziona.
Ricezione dei segnali
Per cominciare, è necessario registrare il segnale. L'intero segnale ha una durata di soli 120 microsecondi, quindi per analizzare comodamente i suoi componenti, è desiderabile un ricevitore SDR con una frequenza di campionamento di almeno 5 MHz.
Dopo la registrazione riceviamo un file WAV con una frequenza di campionamento di 5000000 campioni/sec, 30 secondi di tale registrazione "pesano" circa 500MB. Ascoltarlo con un lettore multimediale è ovviamente inutile: il file contiene non suoni, ma il segnale radio digitalizzato direttamente: è così che funziona la Software Defined Radio.
Apriremo e elaboreremo il file utilizzando Python. Coloro che desiderano sperimentare autonomamente possono scaricare un esempio di registrazione. .
Carichiamo il file e vediamo cosa c'è dentro.
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()
Risultato: vediamo evidenti "impulsi" sullo sfondo del rumore.
Ogni "impulso" è effettivamente il segnale, la cui struttura è ben visibile se aumentiamo la risoluzione del grafico.
Come possiamo vedere, l'immagine corrisponde perfettamente a quanto descritto sopra. Possiamo procedere con l'elaborazione dei dati.
Decodifica
Per cominciare, è necessario ottenere il flusso di bit. Il segnale stesso è codificato utilizzando la codifica manchester:
Dalla differenza dei livelli nei semi-byte è facile ottenere i "0" e "1" reali.
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"
La struttura del segnale stesso ha il seguente aspetto:
Esaminiamo i campi più nel dettaglio.
DF (Downlink Format, 5 bit) — definisce il tipo di messaggio. Ci sono diversi tipi:
()
Ci interessa solo il tipo DF17, poiché è l'unico che contiene le coordinate dell'aeromobile.
ICAO (24 bit) — codice unico internazionale dell'aeromobile. È possibile controllare l'aereo tramite il suo codice (sfortunatamente, l'autore ha smesso di aggiornare il database, ma è ancora attuale). Ad esempio, per il codice 3c5ee2 abbiamo le seguenti informazioni:
Correzione: in la descrizione del codice ICAO è riportata in modo più dettagliato, consiglio di dare un'occhiata a chi è interessato.
DATA (56 o 112 bit) — i veri dati che noi decodificheremo. I primi 5 bit di dati sono il campo Type Code, che contiene il sotto-tipo dei dati memorizzati (non confondere con DF). Ci sono molti tipi:
()
Esaminiamo alcuni esempi di pacchetti.
Identificazione dell'aeromobile
Esempio in forma binaria:
00100 011 000101 010111 000111 110111 110001 111000
Campi dati:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, ogni simbolo C1-C8 contiene codici corrispondenti agli indici nella stringa:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Decodificando la stringa, è semplice ottenere il codice dell'aereo: 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("Identificazione dell'aereo:", code_str.replace('#', ''))
Posizione in volo
Se con il nome è tutto chiaro, con le coordinate è più complicato. Esse vengono trasmesse sotto forma di 2 pacchetti, pari e dispari. Il codice del campo TC = 01011b = 11.
Esempio di pacchetti pari e dispari:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Il calcolo delle coordinate avviene secondo una formula piuttosto complessa:
()
Non sono un esperto di GIS, quindi non so da dove sia estratto. Chi è informato, può scrivere nei commenti.
L'altezza è più semplice da calcolare: a seconda di un bit specifico, può essere espressa come un multiplo di 25 o 100 piedi.
Velocità in volo
Pacchetto con TC=19. Qui è interessante notare che la velocità può essere sia quella precisa, rispetto al suolo (Ground Speed), sia quella aerea, misurata dal sensore dell'aereo (Airspeed). Vengono trasmessi anche molti altri campi diversi:
()
Conclusione
Come si può vedere, la tecnologia ADS-B è diventata un affascinante connubio, dove uno standard si rivela utile non solo ai professionisti, ma anche agli utenti comuni. Ma, naturalmente, un ruolo chiave in questo lo ha avuto la riduzione dei costi della tecnologia dei ricevitori SDR digitali, che permettono di ricevere segnali a frequenze superiori al gigahertz con dispositivi letteralmente a 'pochi spiccioli'.
Nel protocollo ovviamente c'è molto di più. Gli interessati possono consultare il PDF sulla pagina o visitare il già menzionato .
Probabilmente a molti non interesserà tutto ciò che è stato scritto, ma almeno spero che l'idea generale di come funzioni sia rimasta.
A proposito, esiste già un decodificatore pronto in Python, che si può studiare . E i proprietari di ricevitori SDR possono assemblare e avviare un decodificatore ADS-B pronto , di questo è stato parlato nel .
Il codice sorgente del parser, descritto nell'articolo, è fornito sotto. Questo è un esempio di test, non destinato alla produzione, ma funziona e può essere utilizzato per analizzare il file registrato sopra.
Codice sorgente (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()
Spero che a qualcuno sia interessato, grazie per l'attenzione.
Fonte: habr.com
