Moien Habr. Wahrscheinlech jiddereen deen jeemools Famill oder Frënn op engem Fliger begéint oder gesinn huet, huet de gratis Flightradar24 Service benotzt. Dëst ass e ganz praktesche Wee fir d'Positioun vum Fliger an Echtzäit ze verfolgen.
В De Betribsprinzip vun esou engem Online-Service gouf beschriwwen. Mir wäerten elo virgoen an erauszefannen wéi eng Daten aus dem Fliger op d'Empfangstatioun geschéckt a kritt ginn an se selwer dekodéieren mat Python.
Geschicht
Natierlech ginn Fligerdaten net iwwerdroe fir d'Benotzer op hire Smartphones ze gesinn. De System gëtt ADS-B (Automatic dependent surveillance-broadcast) genannt, a gëtt benotzt fir automatesch Informatioun iwwer de Fliger un de Kontrollzentrum ze vermëttelen - säin Identifizéierer, Koordinaten, Richtung, Geschwindegkeet, Héicht an aner Daten ginn iwwerdroen. Virdrun, virun der Entdeckung vun esou Systemer, konnt de Sender nëmmen e Punkt um Radar gesinn. Dat goung net méi duer wann et ze vill Fligeren waren.
Technesch besteet ADS-B aus engem Sender op engem Fliger deen periodesch Informatiounspäck mat enger zimlech héijer Frequenz vun 1090 MHz schéckt (et ginn aner Modi, awer mir sinn net esou interesséiert, well d'Koordinaten nëmmen hei iwwerdroe ginn). Natierlech gëtt et nieft dem Sender och iergendwou um Fluchhafen en Empfänger, awer fir eis als Benotzer ass eisen eegenen Empfänger interessant.
Iwwregens, zum Verglach, den éischten esou System, Airnav Radarbox, entworf fir gewéinlech Benotzer, erschéngt am Joer 2007, a kascht ongeféier $ 900; en Abonnement op Netzwierkservicer kascht weider $ 250 d'Joer.
Rezensiounen vun deenen éischte russesche Besëtzer kënnen um Forum gelies ginn . Elo datt RTL-SDR Empfänger wäit verfügbar sinn, kann en ähnlechen Apparat fir $ 30 zesummegesat ginn; méi iwwer dëst war an . Loosst eis op de Protokoll selwer goen - loosst eis kucken wéi et funktionnéiert.
Empfang vun Signaler
Als éischt muss d'Signal opgeholl ginn. De ganze Signal huet eng Dauer vun nëmmen 120 Mikrosekonnen, also fir seng Komponenten bequem ze demontéieren, ass en SDR Empfänger mat enger Samplingsfrequenz vu mindestens 5 MHz wënschenswäert.
No der Opnam kréie mir eng WAV-Datei mat engem Samplingsrate vu 5000000 Echantillon / sec; 30 Sekonnen vun esou Opnam "weien" ongeféier 500MB. Mat engem Media Player nolauschteren ass natierlech nëtzlos - de Fichier enthält keen Toun, mee en direkt digitaliséierte Radiosignal - genee esou funktionnéiert Software Defined Radio.
Mir wäerten d'Datei opmaachen a veraarbecht mat Python. Déi, déi eleng experimentéiere wëllen, kënnen eng Beispillopnam eroflueden .
Loosst eis d'Datei eroflueden a kucken wat dobannen ass.
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()
Resultat: mir gesinn offensichtlech "Puls" géint den Hannergrondgeräusch.
All "Puls" ass e Signal, d'Struktur vun deem ass kloer ze gesinn wann Dir d'Resolutioun op der Grafik erhéicht.
Wéi Dir kënnt gesinn, ass d'Bild ganz konsequent mat deem wat an der Beschreiwung hei uewen uginn ass. Dir kënnt d'Veraarbechtung vun den Donnéeën ufänken.
Decodéieren
Als éischt musst Dir e bësse Stream kréien. D'Signal selwer gëtt mat Manchester Kodéierung kodéiert:
Vum Niveaudifferenz an Nibbles ass et einfach richteg "0" an "1" ze kréien.
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"
D'Struktur vum Signal selwer ass wéi follegt:
Loosst eis d'Felder méi detailléiert kucken.
DF (Downlink Format, 5 Bits) - bestëmmt den Typ vu Message. Et gi verschidden Zorte:
()
Mir sinn nëmme fir den Typ DF17 interesséiert, well ... Et ass dëst deen d'Koordinate vum Fliger enthält.
ICAO (24 Bits) - international eenzegaarteg Code vum Fliger. Dir kënnt de Fliger duerch säi Code kontrolléieren (leider huet den Auteur opgehalen d'Datebank ze aktualiséieren, awer se ass nach ëmmer relevant). Zum Beispill, fir Code 3c5ee2 hu mir déi folgend Informatioun:
Edit: an D'Beschreiwung vum ICAO Code gëtt méi detailléiert uginn; Ech recommandéieren déi interesséiert et ze liesen.
DATA (56 oder 112 Bits) - déi tatsächlech Donnéeën déi mir decodéieren. Déi éischt 5 Bits vun Daten sinn d'Feld Typ Code, mat der Ënnertyp vun den Daten déi gespäichert ginn (net mat DF ze verwiessele). Et gi zimmlech e puer vun dësen Typen:
()
Loosst eis e puer Beispiller vu Packagen kucken.
Fliger Identifikatioun
Beispill a binärer Form:
00100 011 000101 010111 000111 110111 110001 111000
Daten Felder:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, all Charakter C1-C8 enthält Coden entspriechend Indizes an der Linn:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Andeems Dir de String dekodéiert, ass et einfach de Fligercode ze kréien: 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('#', ''))
Airborne Positioun
Wann den Numm einfach ass, da sinn d'Koordinate méi komplizéiert. Si ginn a Form vun 2, souguer an komesch Rummen iwwerdroen. Feldcode TC = 01011b = 11.
Beispill vu souguer an komesch Pakete:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
D'Berechnung vun de Koordinate selwer geschitt no enger zimlech komplizéierter Formel:
()
Ech sinn keen GIS Expert, also weess ech net wou et hierkënnt. Wien weess, schreift an de Kommentaren.
Héicht gëtt als méi einfach ugesinn - ofhängeg vum spezifesche Bit kann se entweder als Multiple vu 25 oder 100 Féiss duergestallt ginn.
Airborne Geschwindegkeet
Package mat TC=19. Déi interessant Saach hei ass datt d'Geschwindegkeet entweder präzis ka sinn, relativ zum Buedem (Ground Speed), oder Airborne, gemooss vun engem Fligersensor (Airspeed). Vill verschidde Felder ginn och iwwerdroen:
()
Konklusioun
Wéi Dir kënnt gesinn, ass d'ADS-B Technologie eng interessant Symbios ginn, wann e Standard net nëmme fir Professionnelen nëtzlech ass, awer och fir gewéinlech Benotzer. Awer natierlech, eng Schlësselroll an dësem gouf vun der méi bëlleger Technologie vun digitale SDR Empfänger gespillt, wat den Apparat erlaabt wuertwiertlech Signaler mat Frequenzen iwwer e Gigahertz "fir Pennys" ze kréien.
Am Standard selwer gëtt et natierlech vill méi. Déi interesséiert kënnen den PDF op der Säit kucken oder besicht déi schonn uewen ernimmt .
Et ass onwahrscheinlech datt all déi uewe fir vill nëtzlech sinn, awer op d'mannst déi allgemeng Iddi wéi et funktionnéiert, hoffen ech, bleift.
Iwwregens, e fäerdege Decoder am Python existéiert schonn, Dir kënnt et studéieren . A Besëtzer vun SDR Empfänger kënnen e fäerdegen ADS-B Decoder montéieren a starten , Dëst gouf méi am Detail diskutéiert an .
De Quellcode vum Parser, deen am Artikel beschriwwe gëtt, gëtt ënner dem Schnëtt uginn. Dëst ass en Testbeispiel dat sech net als Produktioun mécht, awer e puer Saache funktionnéieren dran, an et kann benotzt ginn fir d'Datei hei uewen opgeholl ze analyséieren.
Source Code (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()
Ech hoffen een interesséiert war, merci fir Är Opmierksamkeet.
Source: will.com