Hallo Habr. Waarschijnlijk heeft iedereen die ooit familie of vrienden op het vliegtuig heeft uitgezwaaid, gebruikgemaakt van de gratis service Flightradar24. Dit is een zeer handige manier om de positie van een vliegtuig in real-time te volgen.
In werd het principe van zo'n online service beschreven. Nu gaan we verder en ontdekken we welke gegevens van het vliegtuig naar het ontvangstation worden verzonden en ontvangen, en decoderen we deze zelf met behulp van Python.
Geschiedenis
Het is duidelijk dat de gegevens over vliegtuigen niet worden verzonden zodat gebruikers ze op hun smartphones kunnen zien. Het systeem heet ADS-B (Automatic Dependent Surveillance—Broadcast) en dient voor de automatische verzending van informatie over het vliegtuig naar het verkeersleidingscentrum — informatie over het identificatienummer, coördinaten, richting, snelheid, hoogte en andere gegevens worden verzonden. Vroeger, voordat zulke systemen bestonden, kon de verkeersleider slechts een punt op de radar zien. Dit was niet voldoende toen er te veel vliegtuigen waren.
Technisch gezien bestaat ADS-B uit een zender op het vliegtuig die periodiek gegevenspakketten verzendt op een hoge frequentie van 1090 MHz (er zijn ook andere modi, maar die zijn voor onze doeleinden minder interessant, omdat de coördinaten alleen hier worden verzonden). Uiteraard is er, naast de zender, ook een ontvanger ergens op de luchthaven, maar voor ons als gebruikers is onze eigen ontvanger interessant.
Ter vergelijking, het eerste dergelijke systeem, Airnav Radarbox, dat gericht is op gewone gebruikers, verscheen in 2007 en kostte ongeveer 900$, daarnaast kostte het abonnement op de netwerksystemen nog eens ongeveer 250$ per jaar.
De beoordelingen van die eerste Russische eigenaren zijn te lezen op het forum . Tegenwoordig, nu RTL-SDR ontvangers massaal beschikbaar zijn, kan een vergelijkbaar apparaat worden samengesteld voor 30$, hierover was meer te lezen in . Laten we nu verdergaan naar het protocol — laten we kijken hoe het werkt.
Ontvangst van signalen
Om te beginnen moet het signaal worden vastgelegd. Het hele signaal heeft een duur van slechts 120 microseconden, dus om de componenten comfortabel te kunnen ontleden, is een SDR-ontvanger met een samplefrequentie van ten minste 5MHz gewenst.
Na het opnemen ontvangen we een WAV-bestand met een samplefrequentie van 5000000 monsters/s, 30 seconden opname 'weegt' ongeveer 500MB. Het is vanzelfsprekend nutteloos om dit met een mediaspeler te beluisteren — het bestand bevat geen geluid, maar de gedigitaliseerde radiosignalen — zo werkt Software Defined Radio.
We openen en verwerken het bestand met Python. Degenen die zelf willen experimenteren, kunnen een voorbeeld van de opname downloaden. .
Laten we het bestand laden 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 een achtergrond van ruis.
Elke 'puls' is het signaal, waarvan de structuur goed zichtbaar is bij het vergroten van de resolutie op de grafiek.
Zoals te zien is, komt de afbeelding goed overeen met de beschrijving hierboven. We kunnen beginnen met de verwerking van de gegevens.
Decoderen
Om te beginnen moeten we de bitsstroom krijgen. Het signaal is gecodeerd met behulp van Manchester-codering:
Uit de niveauverschillen in de halve bytes kunnen gemakkelijk de echte '0' en '1' worden verkregen.
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"
De structuur van het signaal zelf ziet er als volgt uit:
Laten we de velden nader bekijken.
DF (Downlink Format, 5 bits) — bepaalt het type bericht. Er zijn verschillende types:
()
Wij zijn alleen geïnteresseerd in type DF17, omdat dit de coördinaten van het vliegtuig bevat.
ICAO (24 bits) — internationaal unieke code van het vliegtuig. U kunt het vliegtuig controleren aan de hand van zijn code (helaas heeft de auteur de database niet meer bijgewerkt, maar deze is nog steeds actueel). Bijvoorbeeld, voor de code 3c5ee2 hebben we de volgende informatie:
Correctie: in is de beschrijving van de ICAO-code uitgebreider, geïnteresseerden raad ik aan om zich te verdiepen.
DATA (56 of 112 bits) — de eigenlijke gegevens die we gaan decoderen. De eerste 5 bits van de gegevens zijn het veld Type Code, dat het subtype van de opgeslagen gegevens bevat (niet te verwarren met DF). Er zijn er behoorlijk veel van:
()
Laten we een aantal voorbeelden van pakketten bekijken.
Vliegtuigidentificatie
Voorbeeld in binaire vorm:
00100 011 000101 010111 000111 110111 110001 111000
Gegevensvelden:
+------+------+------+------+------+------+------+------+------+------+
| 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 de indexen in de string:
#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("Vliegtuig Identificatie:", code_str.replace('#', ''))
Vliegende positie
Als het om de naam gaat, is het eenvoudig, maar de coördinaten zijn ingewikkelder. Ze worden verzonden in paren van even en oneven frames. Het veldcode TC = 01011b = 11.
Voorbeeld van een even en een oneven pakket:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
De berekening van de coördinaten gebeurt volgens een vrij ingewikkelde formule:
()
Ik ben geen specialist in GIS, dus ik weet niet waar het vandaan komt. Wie het weet, laat het me weten in de reacties.
De hoogte wordt eenvoudiger berekend — afhankelijk van een bepaalde bit kan deze worden weergegeven als een veelvoud van 25 of 100 voet.
Vliegende snelheid
Pakket met TC=19. Interessant is dat de snelheid zowel precies kan zijn, ten opzichte van de grond (Ground Speed), als luchtvaart, gemeten door de sensor van het vliegtuig (Airspeed). Er worden ook veel verschillende velden verzonden:
()
Conclusie
Zoals te zien is, is de ADS-B-technologie een interessante symbiose, waarbij een bepaalde standaard niet alleen nuttig is voor professionals, maar ook voor gewone gebruikers. Uiteraard speelde de prijsverlaging van digitale SDR-ontvangers een sleutelrol, waardoor men signals met frequenties boven de gigahertz letterlijk 'voor een prikkie' op een apparaat kon ontvangen.
In de standaard zelf is er natuurlijk veel meer. Geïnteresseerden kunnen de PDF op de pagina bekijken of de eerder genoemde pagina bezoeken .
Waarschijnlijk zal niet veel van wat hierboven is geschreven nuttig zijn voor velen, maar hopelijk is het algemene idee achter hoe dit werkt duidelijk.
Trouwens, er bestaat al een kant-en-klare decoder in Python, die je kunt bestuderen . En eigenaren van SDR-ontvangers kunnen een kant-en-klaar ADS-B-decodeerapparaat samenstellen en uitvoeren , meer hierover werd besproken in .
De broncode van de parser, zoals beschreven in het artikel, is hieronder te vinden. Dit is een testvoorbeeld dat niet bedoeld is voor productie, maar er werkt wel iets en je kunt het eerder genoemde bestand ermee 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 dit interessant vond, bedankt voor je aandacht.
Bron: habr.com
