Hej Habr. Nok hver, der har sagt farvel til eller mÞdt familie eller venner i lufthavnen, har brugt den gratis tjeneste Flightradar24. Det er en meget praktisk mÄde at fÞlge med i flyets position i realtid.
I beskrev princippet for, hvordan en sÄdan online-service fungerer. Nu gÄr vi videre og undersÞger, hvilke data der sendes og modtages fra flyet til modtagerstationen, og vi vil selv afkode dem med Python.
Historik
Det er tydeligt, at data om fly ikke sendes, sĂ„ brugerne kan se dem pĂ„ deres smartphones. Systemet kaldes ADSâB (Automatic dependent surveillanceâbroadcast) og tjener til automatisk at transmittere oplysninger om flyet til kontrolcentret â dets ID, koordinater, retning, hastighed, hĂžjde og andre data sendes. Tidligere, fĂžr sĂ„danne systemer kom frem, kunne kontroltĂ„rnet kun se et punkt pĂ„ radaren. Dette blev utilstrĂŠkkeligt, da antallet af fly blev for stort.
Teknisk set bestÄr ADS-B af en transmitter pÄ flyet, der periodisk sender pakker med oplysninger pÄ en tilstrÊkkeligt hÞj frekvens, 1090 MHz (der findes ogsÄ andre tilstande, men de er ikke sÄ interessante for os, da koordinater kun sendes her). Naturligvis er der udover transmitteren ogsÄ en modtager et sted i lufthavnen, men for os som brugere er vores egen modtager det interessante.
Til sammenligning blev det fÞrste sÄdanne system, Airnav Radarbox, designet til almindelige brugere, lanceret i 2007 og kostede omkring 900 dollar, mens abonnementet pÄ netvÊrkstjenester kostede yderligere omkring 250 dollar om Äret.
Anmeldelser fra de fĂžrste russiske ejere kan lĂŠses pĂ„ forummet . Nu, hvor RTL-SDR modtagere er blevet bredt tilgĂŠngelige, kan en tilsvarende enhed samles for 30 dollar, som det blev beskrevet i . Lad os gĂ„ videre til protokollen â lad os se, hvordan det fungerer.
Modtagelse af signaler
FÞrst skal signalet optages. Selve signalet varer kun 120 mikrosekunder, sÄ for komfortabelt at kunne adskille dets komponenter er en SDR-modtager med en samplingfrekvens pÄ mindst 5 MHz Þnskelig.
Efter optagelsen fĂ„r vi en WAV-fil med en samplingsfrekvens pĂ„ 5000000 prĂžver/sek, 30 sekunders optagelse "vejer" omkring 500 MB. At lytte til den med en medieafspiller er selvfĂžlgelig meningslĂžst â filen indeholder ikke lyd, men den direkte digitaliserede radiosignal â sĂ„dan fungerer Software Defined Radio.
Vi vil Äbne og behandle filen ved hjÊlp af Python. De, der Þnsker at eksperimentere selv, kan downloade et eksempel pÄ optagelsen. .
Lad os uploade filen og se, hvad der er indeni.
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()
Resultatet: vi ser tydelige "impulser" mod stĂžjbaggrunden.
Hver "impuls" er selve signalet, hvis struktur er godt synlig, hvis man zoomer ind pÄ grafikken.
Som man kan se, svarer billedet fuldstĂŠndig til det, der er nĂŠvnt ovenfor. Vi kan begynde at behandle dataene.
Dekodning
FÞrst skal vi fÄ bitstrÞmmen. Selve signalet er kodet ved hjÊlp af manchester encoding:
Fra forskellen i niveauer i halvbytes er det let at fÄ de reelle "0" og "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"
Strukturen af selve signalet ser sÄdan ud:
Lad os se nÊrmere pÄ felterne.
DF (Downlink Format, 5 bit) â definerer typen af besked. Der er flere typer:
()
Vi er kun interesseret i typen DF17, da den netop indeholder oplysninger om flyets position.
ICAO (24 bit) â en international unik kode for flyet. Du kan tjekke flyet ved hjĂŠlp af dets kode (desvĂŠrre er forfatteren stoppet med at opdatere databasen, men den er stadig relevant). For eksempel, for koden 3c5ee2 har vi fĂžlgende oplysninger:
Redigering: i er beskrivelsen af ICAO-koden givet mere detaljeret, anbefales det at se nÊrmere pÄ.
DATA (56 eller 112 bit) â selve dataene, som vi vil dekode. De fĂžrste 5 bit af dataene er feltet Type Code, der indeholder subtype af de gemte data (mĂ„ ikke forveksles med DF). Der er mange sĂ„danne typer:
()
Lad os analysere et par eksempler pÄ pakker.
Aircraft identification
Eksempel i binĂŠr form:
00100 011 000101 010111 000111 110111 110001 111000
Datafelterne:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, hver tegn C1-C8 indeholder koder, der svarer til indeksene i strengen:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Ved at afkode strengen er det ikke svÊrt at fÄ flyets kode: 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("Flyidentifikation:", code_str.replace('#', ''))
LuftbÄren position
Hvis navnet er enkelt, er der lidt mere kompleksitet med koordinaterne. De overfĂžres i form af par af lige og ulige rammer. Feltkoden TC = 01011b = 11.
Eksempel pÄ lige og ulige pakker:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Selve beregningen af koordinaterne udfĂžres efter en ret smart formel:
()
Jeg er ikke specialist i GIS, sÄ jeg ved ikke, hvor det kommer fra. Hvis nogen ved det, sÄ skriv det i kommentarfeltet.
HĂžjden beregnes simplere â afhĂŠngigt af et bestemt bit, kan den reprĂŠsenteres som enten en multipel af 25 eller 100 fod.
LuftbÄren hastighed
Pakke med TC=19. Det interessante her er, at hastigheden kan vÊre prÊcis, i forhold til jorden (Ground Speed), og luftbÄren, mÄlt af flyets sensor (Airspeed). Der overfÞres ogsÄ mange forskellige felter:
()
Konklusion
Som man kan se, er ADS-B-teknologien blevet en interessant symbiose, hvor en standard er nyttig for bÄde professionelle og almindelige brugere. Men naturligvis spillede den faldende pris pÄ digitale SDR-modtagere, der gÞr det muligt at modtage signaler med frekvenser over gigaHertz, en afgÞrende rolle.
I standarden er der naturligvis meget mere. De interesserede kan se PDF'en pÄ siden eller besÞge den tidligere nÊvnte .
Det er mÄske ikke relevant for mange, men forhÄbentlig er den overordnede idé om, hvordan det fungerer, blevet forstÄet.
Forresten findes der allerede en fĂŠrdig dekoder i Python, som man kan studere . Ejere af SDR-modtagere kan samle og kĂžre en klar ADS-B dekoder , mere om dette blev beskrevet i .
Kildekoden til parseren, der er beskrevet i artiklen, er givet nedenfor. Dette er et test eksempel, der ikke har til hensigt at vĂŠre production-ready, men der er noget funktionelt ved det, og man kan bruge det til at parse de nĂŠvnte filer.
Kildekode (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()
Jeg hÄber, at nogen fandt det interessant, tak for opmÊrksomheden.
Kilde: habr.com
