Hola Habr. Probablement tothom que hagi conegut o vist familiars o amics en un avió hagi utilitzat el servei gratuït Flightradar24. Aquesta és una manera molt convenient de fer un seguiment de la posició de l'avió en temps real.
В Es va descriure el principi de funcionament d'aquest servei en línia. Ara avançarem i esbrinarem quines dades s'envien i reben des de l'avió a l'estació receptora i les descodifiquem nosaltres mateixos amb Python.
Història
Òbviament, les dades de l'aeronau no es transmeten perquè els usuaris les vegin als seus telèfons intel·ligents. El sistema s'anomena ADS-B (Automatic dependent surveillance—broadcast) i s'utilitza per transmetre automàticament informació sobre l'aeronau al centre de control: es transmeten el seu identificador, coordenades, direcció, velocitat, altitud i altres dades. Anteriorment, abans de l'arribada d'aquests sistemes, el despatxador només podia veure un punt al radar. Això ja no era suficient quan hi havia massa avions.
Tècnicament, l'ADS-B consisteix en un transmissor en una aeronau que envia periòdicament paquets d'informació a una freqüència força alta de 1090 MHz (hi ha altres modes, però no ens interessen tant, ja que les coordenades només es transmeten aquí). Per descomptat, a més de l'emissor, també hi ha un receptor en algun lloc de l'aeroport, però per a nosaltres, com a usuaris, és interessant el nostre propi receptor.
Per cert, per comparar, el primer sistema d'aquest tipus, Airnav Radarbox, dissenyat per a usuaris normals, va aparèixer el 2007 i va costar uns 900 dòlars; una subscripció als serveis de xarxa va costar 250 dòlars més a l'any.
Les ressenyes dels primers propietaris russos es poden llegir al fòrum . Ara que els receptors RTL-SDR s'han fet àmpliament disponibles, es pot muntar un dispositiu similar per 30 dòlars; més informació sobre això es troba a . Passem al protocol en si, veurem com funciona.
Recepció de senyals
En primer lloc, cal gravar el senyal. Tot el senyal té una durada de només 120 microsegons, de manera que per desmuntar còmodament els seus components, és desitjable un receptor SDR amb una freqüència de mostreig d'almenys 5 MHz.
Després de la gravació, rebem un fitxer WAV amb una freqüència de mostreig de 5000000 de mostres/s; 30 segons d'aquesta gravació "pesen" uns 500 MB. Escoltar-lo amb un reproductor multimèdia, per descomptat, no serveix per a res: el fitxer no conté so, sinó un senyal de ràdio digitalitzat directament; així és exactament com funciona la ràdio definida per programari.
Obrirem i processarem el fitxer amb Python. Aquells que vulguin experimentar pel seu compte poden descarregar-se un exemple de gravació .
Baixem el fitxer i veiem què hi ha dins.
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()
Resultat: veiem "polsos" evidents contra el soroll de fons.
Cada "pols" és un senyal, l'estructura del qual és clarament visible si augmenta la resolució del gràfic.
Com podeu veure, la imatge és força coherent amb el que es dóna a la descripció anterior. Podeu començar a processar les dades.
Descodificació
En primer lloc, cal obtenir una mica de flux. El senyal en si es codifica mitjançant la codificació Manchester:
A partir de la diferència de nivell en els mordis, és fàcil obtenir "0" i "1" reals.
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"
L'estructura del senyal en si és la següent:
Vegem els camps amb més detall.
DF (Format d'enllaç descendent, 5 bits): determina el tipus de missatge. Hi ha diversos tipus:
()
Només ens interessa el tipus DF17, perquè... És aquest el que conté les coordenades de l'avió.
OACI (24 bits) - codi únic internacional de l'aeronau. Podeu comprovar l'avió pel seu codi (malauradament, l'autor ha deixat d'actualitzar la base de dades, però encara és rellevant). Per exemple, per al codi 3c5ee2 tenim la informació següent:
Edita: en La descripció del codi OACI es dóna amb més detall; recomano que els interessats el llegiu.
DATA (56 o 112 bits): les dades reals que descodificarem. Els primers 5 bits de dades són el camp Tipus de codi, que conté el subtipus de les dades que s'emmagatzemen (no s'ha de confondre amb DF). Hi ha bastants d'aquests tipus:
()
Vegem alguns exemples de paquets.
Identificació de l'aeronau
Exemple en forma binària:
00100 011 000101 010111 000111 110111 110001 111000
Camps de dades:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, cada caràcter C1-C8 conté codis corresponents als índexs de la línia:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
En descodificar la cadena, és fàcil obtenir el codi de l'avió: 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('#', ''))
Posició aerotransportada
Si el nom és senzill, les coordenades són més complicades. Es transmeten en forma de 2 fotogrames, parells i senars. Codi de camp TC = 01011b = 11.
Exemple de paquets parells i senars:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
El càlcul de les coordenades en si es produeix segons una fórmula força complicada:
()
No sóc un expert en GIS, així que no sé d'on ve. Qui sap, escriu als comentaris.
L'alçada es considera més simple: depenent del bit específic, es pot representar com un múltiple de 25 o 100 peus.
Velocitat aerotransportada
Paquet amb TC=19. El més interessant aquí és que la velocitat pot ser precisa, en relació amb el terra (Ground Speed) o en l'aire, mesurada per un sensor d'avió (Airspeed). També es transmeten molts camps diferents:
()
Conclusió
Com podeu veure, la tecnologia ADS-B s'ha convertit en una simbiosi interessant, quan un estàndard és útil no només per als professionals, sinó també per als usuaris normals. Però, per descomptat, un paper clau en això va tenir la tecnologia més barata dels receptors digitals SDR, que permet que el dispositiu rebi literalment senyals amb freqüències per sobre d'un gigahertz "per cèntims".
A l'estàndard en si, és clar, hi ha molt més. Els interessats poden consultar el PDF a la pàgina o visiteu el ja esmentat anteriorment .
És poc probable que tot l'anterior sigui útil per a molts, però almenys la idea general de com funciona, espero, es manté.
Per cert, ja existeix un descodificador fet a Python, podeu estudiar-lo . I els propietaris de receptors SDR poden muntar i llançar un descodificador ADS-B preparat , això es va parlar amb més detall a .
El codi font de l'analitzador descrit a l'article es mostra a sota del tall. Aquest és un exemple de prova que no pretén ser producció, però algunes coses hi funcionen i es pot utilitzar per analitzar el fitxer gravat anteriorment.
Codi font (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 que algú estigui interessat, gràcies per la vostra atenció.
Font: www.habr.com