Hallo Habr. Seker almal wat al familielede of vriende op 'n vliegtuig ontmoet of afgesien het, het die gratis Flightradar24-diens gebruik. Dit is 'n baie gerieflike manier om die vliegtuig se posisie in reële tyd op te spoor.
В Die bedryfsbeginsel van so 'n aanlyn diens is beskryf. Ons sal nou voortgaan en uitvind watter data van die vliegtuig na die ontvangstasie gestuur en ontvang word en dit self met Python dekodeer.
Story
Dit is duidelik dat vliegtuigdata nie oorgedra word vir gebruikers om op hul slimfone te sien nie. Die stelsel word ADS-B (Automatic dependent surveillance—broadcast) genoem, en word gebruik om inligting oor die vliegtuig outomaties na die beheersentrum oor te dra – sy identifiseerder, koördinate, rigting, spoed, hoogte en ander data word versend. Voorheen, voor die koms van sulke stelsels, kon die versender net 'n punt op die radar sien. Dit was nie meer genoeg as daar te veel vliegtuie was nie.
Tegnies bestaan ADS-B uit 'n sender op 'n vliegtuig wat periodiek pakkies inligting teen 'n redelik hoë frekwensie van 1090 MHz stuur (daar is ander modusse, maar ons is nie so geïnteresseerd in hulle nie, aangesien koördinate net hier versend word). Natuurlik, benewens die sender, is daar ook iewers op die lughawe 'n ontvanger, maar vir ons as gebruikers is ons eie ontvanger interessant.
Terloops, ter vergelyking, die eerste so 'n stelsel, Airnav Radarbox, ontwerp vir gewone gebruikers, het in 2007 verskyn en het sowat $900 gekos; 'n intekening op netwerkdienste het nog $250 per jaar gekos.
Resensies van daardie eerste Russiese eienaars kan op die forum gelees word . Noudat RTL-SDR-ontvangers wyd beskikbaar geword het, kan 'n soortgelyke toestel vir $30 saamgestel word; meer hieroor was in . Kom ons gaan aan na die protokol self – kom ons kyk hoe dit werk.
Ontvangs van seine
Eerstens moet die sein aangeteken word. Die hele sein het 'n duur van slegs 120 mikrosekondes, dus om sy komponente gemaklik uitmekaar te haal, is 'n SDR-ontvanger met 'n monsterfrekwensie van minstens 5 MHz wenslik.
Na opname ontvang ons 'n WAV-lêer met 'n steekproeftempo van 5000000 monsters/sek; 30 sekondes van so 'n opname "weeg" ongeveer 500MB. Om daarna met 'n mediaspeler te luister, is natuurlik nutteloos - die lêer bevat nie klank nie, maar 'n direk gedigitaliseerde radiosein - dit is presies hoe Software Defined Radio werk.
Ons sal die lêer oopmaak en verwerk met Python. Diegene wat op hul eie wil eksperimenteer, kan 'n voorbeeldopname aflaai .
Kom ons laai die lêer af en kyk wat binne is.
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()
Gevolg: ons sien duidelike "pulse" teen die agtergrondgeraas.
Elke "puls" is 'n sein, waarvan die struktuur duidelik sigbaar is as jy die resolusie op die grafiek verhoog.
Soos u kan sien, stem die prentjie redelik ooreen met wat in die beskrywing hierbo gegee word. U kan die data begin verwerk.
Dekodering
Eerstens moet jy 'n bietjie stroom kry. Die sein self word geënkodeer met Manchester-kodering:
Van die vlakverskil in nibbles is dit maklik om regte "0" en "1" te kry.
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"
Die struktuur van die sein self is soos volg:
Kom ons kyk na die velde in meer besonderhede.
DF (Afskakel-formaat, 5 bisse) - bepaal die tipe boodskap. Daar is verskeie tipes:
()
Ons stel net belang in tipe DF17, want... Dit is dit wat die koördinate van die vliegtuig bevat.
ICAO (24 bisse) - internasionale unieke kode van die vliegtuig. Jy kan die vliegtuig deur sy kode nagaan (ongelukkig het die skrywer opgehou om die databasis by te werk, maar dit is steeds relevant). Byvoorbeeld, vir kode 3c5ee2 het ons die volgende inligting:
Redigeer: in Die beskrywing van die ICAO-kode word in meer besonderhede gegee; Ek beveel aan dat belangstellendes dit lees.
DATA (56 of 112 bisse) - die werklike data wat ons sal dekodeer. Die eerste 5 stukkies data is die veld Tik kode, wat die subtipe bevat van die data wat gestoor word (moet nie met DF verwar word nie). Daar is 'n hele paar van hierdie tipes:
()
Kom ons kyk na 'n paar voorbeelde van pakkette.
Vliegtuig identifikasie
Voorbeeld in binêre vorm:
00100 011 000101 010111 000111 110111 110001 111000
Datavelde:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, elke karakter C1-C8 bevat kodes wat ooreenstem met indekse in die reël:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Deur die string te dekodeer, is dit maklik om die vliegtuigkode te kry: 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('#', ''))
Posisie in die lug
As die naam eenvoudig is, dan is die koördinate meer ingewikkeld. Hulle word oorgedra in die vorm van 2, ewe en onewe rame. Veldkode TC = 01011b = 11.
Voorbeeld van ewe en onewe pakkies:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Die berekening van koördinate self vind plaas volgens 'n taamlik moeilike formule:
()
Ek is nie 'n GIS-kenner nie, so ek weet nie waar dit vandaan kom nie. Wie weet, skryf in die kommentaar.
Hoogte word as eenvoudiger beskou - afhangende van die spesifieke bietjie, kan dit as 'n veelvoud van 25 of 100 voet voorgestel word.
Snelheid in die lug
Pakket met TC=19. Die interessante ding hier is dat die spoed óf akkuraat kan wees, relatief tot die grond (Ground Speed), óf in die lug, gemeet deur 'n vliegtuigsensor (Airspeed). Baie verskillende velde word ook oorgedra:
()
Gevolgtrekking
Soos u kan sien, het ADS-B-tegnologie 'n interessante simbiose geword, wanneer 'n standaard nie net vir professionele persone nuttig is nie, maar ook vir gewone gebruikers. Maar natuurlik is 'n sleutelrol hierin gespeel deur die goedkoper tegnologie van digitale SDR-ontvangers, wat die toestel in staat stel om letterlik seine te ontvang met frekwensies bo 'n gigahertz "vir pennies."
In die standaard self is daar natuurlik baie meer. Belangstellendes kan die PDF op die bladsy sien of besoek die een wat reeds hierbo genoem is .
Dit is onwaarskynlik dat al die bogenoemde vir baie nuttig sal wees, maar die algemene idee van hoe dit werk, hoop ek, bly ten minste.
Terloops, 'n klaargemaakte dekodeerder in Python bestaan reeds, jy kan dit bestudeer . En eienaars van SDR-ontvangers kan 'n klaargemaakte ADS-B-dekodeerder saamstel en bekendstel , is dit in meer besonderhede bespreek in .
Die bronkode van die ontleder wat in die artikel beskryf word, word onder die snit gegee. Dit is 'n toetsvoorbeeld wat nie voorgee om produksie te wees nie, maar sommige dinge werk daarin, en dit kan gebruik word om die lêer wat hierbo aangeteken is, te ontleed.
Bronkode (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()
Ek hoop iemand het belanggestel, dankie vir jou aandag.
Bron: will.com