Hej Habr. Sandsynligvis har alle, der nogensinde har mødt eller set slægtninge eller venner af på et fly, brugt den gratis Flightradar24-tjeneste. Dette er en meget bekvem måde at spore flyets position i realtid.
В Funktionsprincippet for en sådan onlinetjeneste blev beskrevet. Vi vil nu gå videre og finde ud af, hvilke data der sendes og modtages fra flyet til modtagestationen og selv afkode det ved hjælp af Python.
Story
Det er klart, at flydata ikke transmitteres, så brugerne kan se dem på deres smartphones. Systemet kaldes ADS-B (Automatic dependent surveillance—broadcast), og bruges til automatisk at sende information om flyet til kontrolcentret - dets identifikator, koordinater, retning, hastighed, højde og andre data transmitteres. Tidligere, før fremkomsten af sådanne systemer, kunne afsenderen kun se et punkt på radaren. Dette var ikke længere nok, når der var for mange fly.
Teknisk set består ADS-B af en sender på et fly, der med jævne mellemrum sender informationspakker med en ret høj frekvens på 1090 MHz (der er andre tilstande, men vi er ikke så interesserede i dem, da koordinater kun transmitteres her). Udover senderen er der selvfølgelig også en modtager et sted i lufthavnen, men for os som brugere er vores egen modtager interessant.
Til sammenligning dukkede det første sådan system, Airnav Radarbox, designet til almindelige brugere, op i 2007 og kostede omkring $900; et abonnement på netværkstjenester kostede yderligere $250 om året.
Anmeldelser af de første russiske ejere kan læses på forummet . Nu hvor RTL-SDR-modtagere er blevet bredt tilgængelige, kan en lignende enhed samles for $30; mere om dette var i . Lad os gå videre til selve protokollen – lad os se, hvordan den virker.
Modtagelse af signaler
Først skal signalet optages. Hele signalet har en varighed på kun 120 mikrosekunder, så for komfortabelt at adskille dets komponenter er en SDR-modtager med en samplingsfrekvens på mindst 5 MHz ønskelig.
Efter optagelse modtager vi en WAV-fil med en samplingshastighed på 5000000 samples/sek; 30 sekunder af en sådan optagelse "vejer" omkring 500 MB. At lytte til det med en medieafspiller er selvfølgelig nytteløst – filen indeholder ikke lyd, men et direkte digitaliseret radiosignal – det er præcis sådan Software Defined Radio fungerer.
Vi åbner og behandler filen ved hjælp af Python. De, der ønsker at eksperimentere på egen hånd, kan downloade en eksempeloptagelse .
Lad os downloade filen og se, hvad der er indeni.
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: vi ser tydelige "impulser" mod baggrundsstøjen.
Hver "impuls" er et signal, hvis struktur er tydeligt synlig, hvis du øger opløsningen på grafen.
Som du kan se, stemmer billedet ganske overens med det, der er givet i beskrivelsen ovenfor. Du kan begynde at behandle dataene.
Afkodning
Først skal du få lidt stream. Selve signalet er kodet ved hjælp af Manchester-kodning:
Fra niveauforskellen i nibbles er det nemt at få ægte "0" og "1".
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"
Strukturen af selve signalet er som følger:
Lad os se på felterne mere detaljeret.
DF (Downlink Format, 5 bit) - bestemmer typen af besked. Der er flere typer:
()
Vi er kun interesserede i type DF17, fordi... Det er denne, der indeholder flyets koordinater.
ICAO (24 bit) - international unik kode for flyet. Du kan tjekke flyet ved dets kode (forfatteren er desværre stoppet med at opdatere databasen, men den er stadig relevant). For eksempel, for kode 3c5ee2 har vi følgende oplysninger:
Edit: i Beskrivelsen af ICAO-koden er givet mere detaljeret; jeg anbefaler, at interesserede læser den.
DATA (56 eller 112 bit) - de faktiske data, som vi vil afkode. De første 5 databit er feltet Type kode, der indeholder undertypen af de data, der lagres (ikke at forveksle med DF). Der er en del af disse typer:
()
Lad os se på et par eksempler på pakker.
Identifikation af fly
Eksempel i binær form:
00100 011 000101 010111 000111 110111 110001 111000
Datafelter:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, hvert tegn C1-C8 indeholder koder svarende til indekser i linjen:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Ved at afkode strengen er det nemt at få flykoden: 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('#', ''))
Luftbåren position
Hvis navnet er enkelt, så er koordinaterne mere komplicerede. De transmitteres i form af 2, lige og ulige rammer. Feltkode TC = 01011b = 11.
Eksempel på lige og ulige pakker:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Selve koordinatberegningen sker efter en ret vanskelig formel:
()
Jeg er ikke GIS-ekspert, så jeg ved ikke, hvor det kommer fra. Hvem ved, skriv i kommentarerne.
Højde anses for at være enklere - afhængigt af den specifikke bit kan den repræsenteres som enten et multiplum af 25 eller 100 fod.
Luftbåren hastighed
Pakke med TC=19. Det interessante her er, at hastigheden enten kan være nøjagtig, i forhold til jorden (Ground Speed), eller luftbåren, målt af en flysensor (Airspeed). Mange forskellige felter overføres også:
()
Konklusion
Som du kan se, er ADS-B-teknologi blevet en interessant symbiose, når en standard er nyttig ikke kun for fagfolk, men også for almindelige brugere. Men selvfølgelig blev en nøglerolle i dette spillet af den billigere teknologi af digitale SDR-modtagere, som tillader enheden bogstaveligt talt at modtage signaler med frekvenser over en gigahertz "for øre."
I selve standarden er der selvfølgelig meget mere. Interesserede kan se PDF'en på siden eller besøg den, der allerede er nævnt ovenfor .
Det er usandsynligt, at alt ovenstående vil være nyttigt for mange, men i det mindste forbliver den generelle idé om, hvordan det fungerer, forhåbentlig.
Forresten findes der allerede en færdig dekoder i Python, du kan studere den . Og ejere af SDR-modtagere kan samle og lancere en færdiglavet ADS-B-dekoder , dette blev diskuteret mere detaljeret i .
Kildekoden for parseren, der er beskrevet i artiklen, er givet under udskæringen. Dette er et testeksempel, der ikke foregiver at være produktion, men nogle ting fungerer i det, og det kan bruges til at parse filen optaget ovenfor.
Kildekode (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()
Jeg håber nogen var interesseret, tak for din opmærksomhed.
Kilde: www.habr.com