Hallo Habr. Waarschijnlijk heeft iedereen die ooit familieleden of vrienden in een vliegtuig heeft ontmoet of gezien, gebruik gemaakt van de gratis dienst Flightradar24. Dit is een zeer handige manier om de positie van het vliegtuig in realtime te volgen.
В Het werkingsprincipe van een dergelijke online dienst werd beschreven. We gaan nu uitzoeken welke gegevens er van het vliegtuig naar het ontvangststation worden verzonden en ontvangen en deze zelf decoderen met Python.
Verhaal
Het is duidelijk dat vliegtuiggegevens niet worden verzonden zodat gebruikers deze op hun smartphones kunnen bekijken. Het systeem heet ADS-B (Automatic depend surveillance-broadcast) en wordt gebruikt om automatisch informatie over het vliegtuig naar het controlecentrum te verzenden - de identificatie, coördinaten, richting, snelheid, hoogte en andere gegevens worden verzonden. Voorheen, vóór de komst van dergelijke systemen, kon de coördinator alleen een punt op de radar zien. Dit was niet langer voldoende toen er te veel vliegtuigen waren.
Technisch gezien bestaat ADS-B uit een zender in een vliegtuig die periodiek informatiepakketten verzendt met een vrij hoge frequentie van 1090 MHz (er zijn andere modi, maar daar zijn we niet zo in geïnteresseerd, omdat coördinaten alleen hier worden verzonden). Uiteraard staat er naast de zender ook ergens op de luchthaven een ontvanger, maar voor ons als gebruikers is een eigen ontvanger interessant.
Ter vergelijking: het eerste dergelijke systeem, Airnav Radarbox, ontworpen voor gewone gebruikers, verscheen in 2007 en kostte ongeveer $900; een abonnement op netwerkdiensten kostte nog eens $250 per jaar.
Recensies van die eerste Russische eigenaren zijn te lezen op het forum . Nu RTL-SDR-ontvangers algemeen verkrijgbaar zijn, kan een soortgelijk apparaat voor $ 30 in elkaar worden gezet; hierover stond meer in . Laten we verder gaan met het protocol zelf - laten we kijken hoe het werkt.
Signalen ontvangen
Eerst moet het signaal worden opgenomen. Het gehele signaal heeft een duur van slechts 120 microseconden, dus om de componenten comfortabel te demonteren is een SDR-ontvanger met een bemonsteringsfrequentie van minimaal 5 MHz wenselijk.
Na de opname ontvangen we een WAV-bestand met een bemonsteringssnelheid van 5000000 samples/sec; 30 seconden van zo'n opname "weegt" ongeveer 500 MB. Luisteren met een mediaspeler heeft uiteraard geen zin - het bestand bevat geen geluid, maar een direct gedigitaliseerd radiosignaal - dit is precies hoe Software Defined Radio werkt.
We zullen het bestand openen en verwerken met Python. Wie zelf wil experimenteren, kan een voorbeeldopname downloaden .
Laten we het bestand downloaden en kijken wat erin zit.
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()
Resultaat: we zien duidelijke “pulsen” tegen het achtergrondgeluid.
Elke “puls” is een signaal waarvan de structuur duidelijk zichtbaar is als je de resolutie in de grafiek verhoogt.
Zoals u kunt zien, komt het beeld redelijk overeen met wat in de bovenstaande beschrijving wordt gegeven. U kunt beginnen met het verwerken van de gegevens.
decoderen
Eerst moet je een beetje stroom krijgen. Het signaal zelf is gecodeerd met behulp van Manchester-codering:
Uit het niveauverschil in hapjes is het gemakkelijk om echte “0” en “1” te krijgen.
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"
De structuur van het signaal zelf is als volgt:
Laten we de velden in meer detail bekijken.
DF (Downlink-formaat, 5 bits) - bepaalt het type bericht. Er zijn verschillende soorten:
()
Wij zijn alleen geïnteresseerd in type DF17, omdat... Hierin staan de coördinaten van het vliegtuig.
ICAO (24 bits) - internationale unieke code van het vliegtuig. Je kunt het vliegtuig controleren aan de hand van de code (helaas is de auteur gestopt met het updaten van de database, maar deze is nog steeds relevant). Voor code 3c5ee2 hebben we bijvoorbeeld de volgende informatie:
Bewerken: binnen De beschrijving van de ICAO-code wordt gedetailleerder gegeven; ik raad geïnteresseerden aan deze te lezen.
GEGEVENS (56 of 112 bits) - de feitelijke gegevens die we zullen decoderen. De eerste 5 bits aan gegevens vormen het veld Type Code, met daarin het subtype van de gegevens die worden opgeslagen (niet te verwarren met DF). Er zijn nogal wat van deze soorten:
()
Laten we een paar voorbeelden van pakketten bekijken.
Identificatie van vliegtuigen
Voorbeeld in binaire vorm:
00100 011 000101 010111 000111 110111 110001 111000
Data velden:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, elk teken C1-C8 bevat codes die overeenkomen met indices in de regel:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_#############0123456789######
Door de string te decoderen, is het eenvoudig om de vliegtuigcode te verkrijgen: 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('#', ''))
Positie in de lucht
Als de naam eenvoudig is, zijn de coördinaten ingewikkelder. Ze worden verzonden in de vorm van 2, even en oneven frames. Veldcode TC = 01011b = 11.
Voorbeeld van even en oneven pakketten:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
De berekening van de coördinaten zelf gebeurt volgens een nogal lastige formule:
()
Ik ben geen GIS-expert, dus ik weet niet waar het vandaan komt. Wie weet, schrijf het in de reacties.
Hoogte wordt als eenvoudiger beschouwd: afhankelijk van het specifieke bit kan deze worden weergegeven als een veelvoud van 25 of 100 meter.
Snelheid in de lucht
Pakket met TC=19. Het interessante hier is dat de snelheid nauwkeurig kan zijn, ten opzichte van de grond (Ground Speed), of in de lucht, gemeten door een vliegtuigsensor (Airspeed). Er worden ook veel verschillende velden verzonden:
()
Conclusie
Zoals u kunt zien, is de ADS-B-technologie een interessante symbiose geworden, waarbij een standaard niet alleen nuttig is voor professionals, maar ook voor gewone gebruikers. Maar uiteraard speelde de goedkopere technologie van digitale SDR-ontvangers hierin een sleutelrol, waardoor het apparaat letterlijk signalen kan ontvangen met frequenties boven een gigahertz “voor centen”.
In de standaard zelf zit natuurlijk nog veel meer. Geïnteresseerden kunnen de PDF op de pagina bekijken of bezoek de hierboven al genoemde .
Het is onwaarschijnlijk dat al het bovenstaande voor velen nuttig zal zijn, maar ik hoop dat het algemene idee van hoe het werkt in ieder geval blijft bestaan.
Er bestaat trouwens al een kant-en-klare decoder in Python, je kunt deze bestuderen . En eigenaren van SDR-ontvangers kunnen een kant-en-klare ADS-B-decoder samenstellen en lanceren , dit werd in meer detail besproken in .
De broncode van de parser die in het artikel wordt beschreven, wordt onder de tekst weergegeven. Dit is een testvoorbeeld dat niet pretendeert productie te zijn, maar sommige dingen werken erin, en het kan worden gebruikt om het hierboven opgenomen bestand te parseren.
Broncode (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()
Ik hoop dat iemand geïnteresseerd was, bedankt voor uw aandacht.
Bron: www.habr.com