Ciao Habr. Probabilmente chiunque abbia mai accompagnato o salutato parenti o amici all'aeroporto ha utilizzato il servizio gratuito Flightradar24. È un modo molto pratico per monitorare la posizione degli aerei 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 ricevente, e li decodificheremo autonomamente utilizzando Python.
Storia
È ovvio che i dati sugli aerei non vengono trasmessi affinché gli utenti possano vederli sui loro smartphone. Il sistema si chiama ADS–B (Automatic dependent surveillance—broadcast) e serve per la trasmissione automatica delle informazioni sull'aeromobile al centro di controllo — vengono trasmessi il suo identificativo, le coordinate, la direzione, la velocità, l'altezza e altri dati. In passato, prima dell'introduzione di tali sistemi, i controllori potevano vedere solo un punto sul radar. Questo è diventato insufficiente quando il numero di aerei è aumentato considerevolmente.
Tecnicamente, l'ADS-B è composto da un trasmettitore sull'aeromobile, che invia periodicamente pacchetti di informazioni a una frequenza sufficientemente alta di 1090 MHz (ci sono altre modalità, ma non ci interessano molto, dato che le coordinate vengono trasmesse solo qui). Ovviamente, oltre al trasmettitore, esiste anche un ricevitore da qualche parte in aeroporto, ma per noi, come utenti, è interessante il nostro proprio ricevitore.
A proposito, per fare un confronto, il primo sistema di questo tipo, Airnav Radarbox, progettato per utenti normali, è apparso nel 2007 e costava circa 900$, con una sottoscrizione ai servizi di rete che costava circa 250$ all'anno.
Le recensioni di quei primi proprietari russi possono essere lette sul forum . Oggi, con la disponibilità di ricevitori RTL-SDR, è possibile assemblare un dispositivo simile per 30$, di più su questo è stato trattato in . Passiamo così al protocollo — vediamo come funziona.
Ricezione dei segnali
Per iniziare, è necessario registrare il segnale. L'intero segnale ha una durata di soli 120 microsecondi, quindi per analizzare comodamente i suoi componenti è consigliabile 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/secondo; 30 secondi di questa registrazione occupano circa 500 MB. Ascoltarlo con un media player è ovviamente inutile: il file non contiene suono, bensì il segnale radio digitalizzato — è così che funziona la Software Defined Radio.
Apriremo e elaboreremo il file utilizzando Python. Coloro che desiderano sperimentare da soli 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 chiari «impulsi» sullo sfondo del rumore.
Ogni «impulso» è il segnale, la cui struttura è ben visibile se si aumenta la risoluzione nel grafico.
Come possiamo vedere, l'immagine corrisponde perfettamente a quanto descritto sopra. Possiamo procedere all'elaborazione dei dati.
Decodifica
Per iniziare, dobbiamo ottenere il flusso di bit. Il segnale stesso è codificato con manchester encoding:
Dalla differenza di livelli nei mezzi byte, è facile ottenere i reali «0» e «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"
La struttura del segnale è la seguente:
Esaminiamo i campi più in dettaglio.
DF (Downlink Format, 5 bit) — determina il tipo di messaggio. Ci sono diversi tipi:
()
Ci interessa solo il tipo DF17, poiché contiene le coordinate dell'aeromobile.
ICAO (24 bit) — codice unico internazionale dell'aeromobile. Puoi 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 è fornita in modo più dettagliato, consiglio a chi è interessato di informarsi.
DATA (56 o 112 bit) — i dati propri che decodificheremo. I primi 5 bit dei dati sono un campo Type Code, che contiene il sottotipo dei dati memorizzati (non confondere con DF). Ci sono davvero molti di questi tipi:
()
Esaminiamo alcuni esempi di pacchetti.
Identificazione dell'aeromobile
Esempio in formato binario:
00100 011 000101 010111 000111 110111 110001 111000
Campi dei 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 che corrispondono agli indici nella stringa:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Decodificando la stringa, si ottiene facilmente il codice dell'aeromobile: 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 Aeromobile:", code_str.replace('#', ''))
Posizione in volo
Se per il nome è tutto semplice, per le coordinate è più complicato. Vengono trasmesse sotto forma di pacchetti pari e dispari. Codice 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 tramite una formula piuttosto complessa:
()
Non sono un esperto GIS, quindi non so da dove venga estratto. Chiunque ne sappia, scriva nei commenti.
L'altezza è considerata in modo più semplice: a seconda di un certo bit, può essere rappresentata come multipla di 25 o 100 piedi.
Velocità in volo
Pacchetto con TC=19. Quello che è interessante qui è che la velocità può essere sia precisa, rispetto al suolo (Ground Speed), sia aerea, misurata dal sensore dell'aeromobile (Airspeed). Inoltre, vengono trasmessi diversi campi.
()
Conclusione
Come si può notare, la tecnologia ADS-B è diventata una sintesi interessante, dove uno standard è utile non solo ai professionisti, ma anche agli utenti comuni. Naturalmente, un ruolo cruciale è stato svolto dalla riduzione dei costi della tecnologia dei ricevitori SDR digitali, che consentono di ricevere segnali a frequenze superiori a un gigahertz su un dispositivo, letteralmente "a poco prezzo".
Nel proprio standard, ovviamente, c'è molto di più. Gli interessati possono consultare il PDF sulla pagina o visitare il già menzionato .
Difficilmente a molti interesserà quanto scritto sopra, ma almeno spero che sia rimasta chiara l'idea generale di come funzioni.
A proposito, esiste già un decoder pronto in Python, che possono esaminare . E i possessori di ricevitori SDR possono raccogliere e avviare un decoder ADS-B già pronto , di questo si è parlato più dettagliatamente in .
Il codice sorgente del parser, descritto nell'articolo, è fornito sotto. Questo è un esempio di test, che non pretende di essere in produzione, ma funziona in parte 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 sia stato interessante per qualcuno, grazie per l'attenzione.
Fonte: habr.com
