Hallo ЄабŃ. Wahrscheinlich hat jeder mal Verwandte oder Freunde zum Flugzeug begleitet und den kostenlosen Service Flightradar24 genutzt. Dies ist eine sehr praktische Möglichkeit, den Standort eines Flugzeugs in Echtzeit zu verfolgen.
Im Es wurde das Prinzip der Funktionsweise eines solchen Online-Services beschrieben. Jetzt gehen wir weiter und klÀren, welche Daten von der Luftfahrzeug zu der Bodenstation gesendet und empfangen werden, und dekodieren sie mit Python selbst.
Geschichte
Offensichtlich werden die Daten ĂŒber Flugzeuge nicht ĂŒbermittelt, damit die Nutzer sie auf ihren Smartphones sehen können. Das System heiĂt ADSâB (Automatic Dependent Surveillance-Broadcast) und dient zur automatischen Ăbermittlung von Informationen ĂŒber das Luftfahrzeug an das Kontrollzentrum â es werden seine Identifikationsnummer, Koordinaten, Richtung, Geschwindigkeit, Höhe und weitere Daten ĂŒbermittelt. FrĂŒher, vor dem Aufkommen solcher Systeme, konnte der Dispatcher nur einen Punkt auf dem Radar sehen. Das wurde nicht mehr genug, als es zu viele Flugzeuge gab.
Technisch besteht ADS-B aus einem Transmitter im Flugzeug, der regelmĂ€Ăig Pakete mit Informationen auf einer ausreichend hohen Frequenz von 1090 MHz sendet (es gibt auch andere Modi, aber die interessieren uns nicht so sehr, da die Koordinaten nur hier ĂŒbermittelt werden). NatĂŒrlich gibt es neben dem Transmitter auch einen EmpfĂ€nger irgendwo am Flughafen, aber fĂŒr uns als Nutzer ist unser eigener EmpfĂ€nger von Interesse.
Ăbrigens, zum Vergleich: Das erste solche System, Airnav Radarbox, das fĂŒr gewöhnliche Nutzer gedacht war, wurde 2007 eingefĂŒhrt und kostete etwa 900$, und die AbonnementgebĂŒhr fĂŒr die Netzwerkanwendungen betrug noch einmal etwa 250$ pro Jahr.
Die RĂŒckmeldungen der ersten russischen Besitzer können im Forum gelesen werden . Jetzt, da RTL-SDR-EmpfĂ€nger massenhaft verfĂŒgbar geworden sind, kann ein Ă€hnliches GerĂ€t fĂŒr 30$ zusammengestellt werden, mehr darĂŒber wurde in geschrieben. Wir widmen uns nun dem Protokoll â sehen uns an, wie es funktioniert.
Empfang von Signalen
ZunĂ€chst muss das Signal aufgezeichnet werden. Das gesamte Signal hat eine Dauer von nur 120 Mikrosekunden, daher ist ein SDR-EmpfĂ€nger mit einer Abtastrate von mindestens 5 MHz wĂŒnschenswert, um die Komponenten komfortabel auseinanderzuhalten.
Nach der Aufnahme erhalten wir eine WAV-Datei mit einer Abtastrate von 5000000 Samples/Sekunde, 30 Sekunden dieser Aufnahme wiegen etwa 500 MB. Sie mit einem Media Player abzuspielen, macht verstĂ€ndlicherweise keinen Sinn â die Datei enthĂ€lt keinen Ton, sondern ein direkt digitalisiertes Radiosignal â genau so funktioniert Software Defined Radio.
Wir werden die Datei mit Python öffnen und bearbeiten. Wer selbst experimentieren möchte, kann ein Aufnahmebeispiel herunterladen. .
Laden wir die Datei hoch und sehen wir, was sich darin befindet.
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()
Das Ergebnis: Wir sehen deutliche âImpulseâ vor dem Hintergrund von Rauschen.
Jeder âImpulsâ ist das Signal, dessen Struktur gut sichtbar ist, wenn man die Auflösung im Diagramm erhöht.
Wie zu sehen ist, entspricht das Bild gut dem, was oben beschrieben wurde. Wir können mit der Datenverarbeitung beginnen.
Dekodierung
Zuerst mĂŒssen wir den Bitstrom erhalten. Das Signal ist mit Manchester-Codierung codiert:
Aus den Unterschiede der Pegel in den Halbbytes lassen sich leicht reale â0â und â1â ableiten.
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"
Die Struktur des Signals sieht wie folgt aus:
Betrachten wir die Felder genauer.
DF (Downlink Format, 5 Bit) â definiert den Nachrichtentyp. Es gibt mehrere Typen:
()
Uns interessiert nur der Typ DF17, da dieser die Koordinaten des Luftfahrzeugs enthÀlt.
ICAO (24 Bit) â internationaler eindeutiger Code des Luftfahrzeugs. Man kann das Flugzeug anhand seines Codes ĂŒberprĂŒfen (leider hat der Autor die Datenbank nicht aktualisiert, aber sie ist noch aktuell). Zum Beispiel haben wir fĂŒr den Code 3c5ee2 die folgenden Informationen:
Korrektur: in eine detailliertere Beschreibung des ICAO-Codes zu finden, empfehle ich den Interessierten, sich damit auseinanderzusetzen.
DATA (56 oder 112 Bit) â die eigentlichen Daten, die wir dekodieren werden. Die ersten 5 Bit der Daten sind das Feld Typcode, das den Untertyp der gespeicherten Daten enthĂ€lt (nicht mit DF verwechseln). Es gibt viele solcher Typen:
()
Lassen Sie uns einige Beispiele fĂŒr Pakete durchgehen.
Flugzeugidentifikation
Beispiel in binÀrer Form:
00100 011 000101 010111 000111 110111 110001 111000
Datenfelder:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, jedes Zeichen C1-C8 enthÀlt Codes, die den Indizes in der Zeichenkette entsprechen:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Durch das Dekodieren der Zeichenkette lÀsst sich leicht der Flugzeugcode erhalten: 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("Flugzeugidentifikation:", code_str.replace('#', ''))
Luftposition
WĂ€hrend der Name einfach ist, sind die Koordinaten etwas komplizierter. Sie werden in Form von 2, geraden und ungeraden Frames ĂŒbertragen. Der Code fĂŒr das Feld TC = 01011b = 11.
Beispiel fĂŒr gerade und ungerade Pakete:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Die Berechnung der Koordinaten erfolgt nach einer ziemlich raffinierten Formel:
()
Ich bin kein GIS-Spezialist, daher weiĂ ich nicht, woher es stammt. Wer mehr weiĂ, kann gerne in die Kommentare schreiben.
Die Höhenmessung ist einfacher â abhĂ€ngig von einem bestimmten Bit kann sie entweder als Vielfaches von 25 oder 100 FuĂ dargestellt werden.
Luftgeschwindigkeit
Paket mit TC=19. Interessant ist, dass die Geschwindigkeit sowohl genau, relativ zur Erde (Ground Speed), als auch die Luftgeschwindigkeit ist, die vom Sensor des Flugzeugs gemessen wird (Airspeed). Es werden auch viele verschiedene Felder ĂŒbertragen:
()
Fazit
Wie man sehen kann, ist die Technologie ADS-B eine interessante Symbiose, bei der ein Standard nicht nur Profis, sondern auch normalen Nutzern zugutekommt. Aber selbstverstĂ€ndlich spielte die Senkung der Kosten fĂŒr digitale SDR-EmpfĂ€nger eine SchlĂŒsselrolle, die es ermöglichen, Signale mit Frequenzen ĂŒber einem Gigahertz buchstĂ€blich fĂŒr pennies zu empfangen.
Im Standard selbst gibt es natĂŒrlich viel mehr. Interessierte können die PDF auf der Seite einsehen oder die bereits oben erwĂ€hnte .
Es ist unwahrscheinlich, dass das oben Geschriebene fĂŒr viele nĂŒtzlich ist, aber zumindest die allgemeine Idee, wie es funktioniert, hoffe ich, ist erhalten geblieben.
Ăbrigens, ein funktionsfĂ€higer Decoder in Python existiert bereits, er kann genauer studiert werden . Und Besitzer von SDR-EmpfĂ€ngern können einen fertigen ADS-B-Decoder zusammenbauen und starten , ausfĂŒhrlicher darĂŒber wurde in .
Der Quellcode des in dem Artikel beschriebenen Parsers ist unter dem Artikel zu finden. Dies ist ein Testbeispiel, das nicht fĂŒr die Produktion gedacht ist, aber einige Dinge funktionieren darin, und man kann damit die oben aufgezeichnete Datei parsen.
Quellcode (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()
Ich hoffe, es war fĂŒr jemanden interessant, danke fĂŒr die Aufmerksamkeit.
Quelle: habr.com
