Bună Habr. Probabil că toți cei care au însoțit sau au luat la revedere rude sau prieteni pe aeroport au folosit serviciul gratuit Flightradar24. Este o modalitate foarte convenabilă de a urmări poziția unui avion în timp real.
În A fost descris principiul de funcționare al unui astfel de serviciu online. Acum vom merge mai departe și vom descoperi ce date sunt transmise și primite de la aeronava la stația de recepție, și le vom decoda singuri cu ajutorul Python.
Istoria
Este evident că datele despre avioane nu sunt transmise pentru ca utilizatorii să le vadă pe smartphone-urile lor. Sistemul se numește ADS–B (supraveghere automată dependentă—broadcast) și este folosit pentru a transmite automat informații despre aeronave la centrul de control — sunt transmise identificatorul, coordonatele, direcția, viteza, altitudinea și alte date. În trecut, înainte de apariția unor astfel de sisteme, dispecerul putea vedea doar un punct pe radar. Asta a devenit insuficient, când numărul de avioane a crescut prea mult.
Tehnic vorbind, ADS-B constă dintr-un transmițător pe aeronavă, care trimite periodic pachete cu informații la o frecvență destul de înaltă de 1090 MHz (există și alte moduri, dar nu ne interesează atât de mult, deoarece coordonatele sunt transmise doar aici). Evident, pe lângă transmițător, există și un receptor undeva în aeroport, dar pentru noi, ca utilizatori, receptorul care ne interesează este cel propriu.
Apropo, pentru comparație, primul sistem de acest fel, Airnav Radarbox, destinat utilizatorilor obișnuiți, a apărut în 2007 și costa aproximativ 900$, iar abonamentul pentru serviciile de rețea costa încă aproximativ 250$ pe an.
Recenziile primilor proprietari ruși pot fi citite pe forumul . Acum, când receptorii RTL-SDR au devenit larg disponibili, un dispozitiv similar poate fi asamblat cu 30$, mai multe detalii despre acest lucru au fost discutate în . Vom trece acum la protocol — să vedem cum funcționează acesta.
Recepția semnalelor
Pentru început, semnalul trebuie înregistrat. Întregul semnal are o durată de doar 120 microsecunde, deci pentru a analiza confortabil componentele sale, este recomandat un receptor SDR cu o frecvență de eșantionare de cel puțin 5MHz.
După înregistrare, obținem un fișier WAV cu o rată de probare de 5000000 de eșantioane/sec, 30 de secunde de înregistrare „cântăresc” aproximativ 500MB. Ascultarea acestuia cu un media player este, desigur, inutilă — fișierul conține nu sunet, ci semnal radio digitalizat — așa funcționează Software Defined Radio.
Vom deschide și procesa fișierul folosind Python. Cei care doresc să experimenteze pe cont propriu pot descărca exemplul de înregistrare. .
Să încărcăm fișierul și să vedem ce este în interior.
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()
Rezultatul: vedem „impulsuri” clare pe fundalul zgomotului.
Fiecare „impuls” este, de fapt, semnalul a cărui structură este bine vizibilă dacă mărim rezoluția pe grafic.
După cum se poate observa, imaginea se aliniază perfect cu ceea ce am descris mai sus. Putem începe procesarea datelor.
Decodare
În primul rând, trebuie să obținem fluxul de biți. Semnalul este codificat folosind codarea manchester:
Din diferența de niveluri în semi-biturile de date, este ușor să obțineți reale „0” și „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"
Structura semnalului are următorul format:
Să examinăm câmpurile mai în detaliu.
DF (Downlink Format, 5 biți) — determină tipul mesajului. Există mai multe tipuri:
()
Ne interesează doar tipul DF17, deoarece acesta conține coordonatele aeronavei.
ICAO (24 biți) — codul unic internațional al aeronavei. Puteți verifica avionul după codul său (din păcate, autorul a încetat să actualizeze baza, dar aceasta este încă relevantă). De exemplu, pentru codul 3c5ee2 avem următoarele informații:
Corectare: în descrierea codului ICAO este prezentată mai detaliat, celor interesați le recomand să se familiarizeze.
DATA (56 sau 112 biți) — datele propriu-zise pe care le vom decoda. Primele 5 biți de date — câmpul Type Code, conținând subtipul datelor stocate (nu trebuie confundat cu DF). Există destul de multe tipuri:
()
Să analizăm câteva exemple de pachete.
Identificarea aeronavei
Exemplu în formă binară:
00100 011 000101 010111 000111 110111 110001 111000
Câmpuri de date:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, fiecare simbol C1-C8 conține coduri corespunzătoare indicilor din șir:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Decodificând șirul, nu este greu să obții codul avionului: 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("Identificarea aeronavei:", code_str.replace('#', ''))
Poziția airborne
Dacă pentru denumire totul este simplu, coordonatele sunt mai complicate. Acestea sunt transmise sub formă de frame-uri pare și impare. Codul câmpului TC = 01011b = 11.
Exemplu de pachete pare și impare:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Însă calculul coordonatelor se face printr-o formulă destul de ingenioasă:
()
Nu sunt specialist în GIS, așa că nu știu de unde se obține. Cine cunoaște, să scrie în comentarii.
Înălțimea este calculată mai simplu — în funcție de un anumit bit, aceasta poate fi reprezentată fie ca un multiplu de 25, fie de 100 de picioare.
Viteza airborne
Pachet cu TC=19. Este interesant de spus că viteza poate fi atât exactă, relativ la pământ (Ground Speed), cât și aerială, măsurată de senzorul avionului (Airspeed). Se transmit, de asemenea, multe câmpuri diferite:
()
Concluzie
După cum se poate vedea, tehnologia ADS-B a devenit un simbioză interesantă, când un anumit standard este util nu doar profesioniștilor, ci și utilizatorilor obișnuiți. Dar, desigur, rolul cheie în aceasta a fost jucat de reducerea costurilor tehnologiei receptorilor SDR digitale, permițând pe dispozitiv să primească efectiv semnale cu frecvențe mai mari de un gigaherț.
În standard sunt, desigur, mult mai multe. Cei interesați pot viziona PDF-ul de pe pagina sau pot vizita pagina menționată mai sus .
Puțin probabil că multora le va fi util tot ce a fost scris mai sus, dar sper că cel puțin ideea generală despre cum funcționează a rămas.
Apropo, un decodor gata pe Python există deja, poate fi studiat . Iar proprietarii receptorilor SDR pot construi și rula un decodor ADS-B gata , despre care s-a vorbit mai detaliat în .
Codul sursă al parser-ului, descris în articol, este prezentat mai jos. Acesta este un exemplu de test, care nu pretinde a fi în producție, dar are unele funcții care funcționează, iar cu el se poate parsa fișierul înregistrat mai sus.
Codul sursă (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()
Sper că a fost interesant pentru cineva, mulțumesc pentru atenție.
Sursa: habr.com
