Pozdrav Habr. Vjerojatno je svaki onaj tko je barem jednom ispratio rodbinu ili prijatelje na avion koristio besplatnu uslugu Flightradar24. To je vrlo praktičan način praćenja položaja aviona u stvarnom vremenu.
U opisano je načelo rada takve online usluge. Sada ćemo ići dalje i saznati koje se informacije šalju i primaju od zrakoplova do prijemne stanice, te dekodirati ih samostalno pomoću Pythona.
Povijest
Očigledno je da se informacije o zrakoplovima ne šalju kako bi ih korisnici vidjeli na svojim pametnim telefonima. Sustav se naziva ADS–B (Automatic Dependent Surveillance–Broadcast) i služi za automatsko slanje informacija o zrakoplovu u kontrolni centar — šalju se njegov identifikator, koordinate, smjer, brzina, visina i drugi podaci. Prije nego što su se pojavili takvi sustavi, kontrolor je mogao vidjeti samo točku na radaru. To je postalo nedovoljno kada je zrakoplova postalo previše.
Tehnički, ADS-B se sastoji od odašiljača na zrakoplovu, koji povremeno šalje pakete s informacijama na prilično visokoj frekvenciji od 1090 MHz (ima i drugih modova, ali nas oni nisu toliko zanimaju, jer se koordinate šalju samo ovdje). Naravno, osim odašiljača, postoji i prijemnik negdje u zračnoj luci, ali nama, kao korisnicima, zanimljiv je naš vlastiti prijemnik.
Usput, za usporedbu, prvi takav sustav, Airnav Radarbox, namijenjen običnim korisnicima, pojavio se 2007. godine i koštao je oko 900$, dok je godišnja pretplata na mrežne usluge iznosila oko 250$.
Recenzije tih prvih ruskih vlasnika možete pročitati na forumu . Sada, kada su RTL-SDR prijemnici masovno postali dostupni, sličan uređaj može se sastaviti za 30$, detaljnije o tome bilo je u . Mi ćemo se sada prebaciti na protokol - pogledajte kako to funkcionira.
Prijem signala
Za početak, signal treba snimiti. Cjelokupni signal traje samo 120 mikrosekundi, stoga je poželjno imati SDR prijemnik s frekvencijom uzorkovanja od najmanje 5 MHz kako bismo ugodno razdvojili njegove komponente.
Nakon snimanja dobijamo WAV datoteku s frekvencijom uzorkovanja od 5000000 uzoraka/u sekundi, 30 sekundi takvog zapisa 'teži' oko 500MB. Slušati je putem medija plejera naravno, besmisleno je – datoteka ne sadrži zvuk, već direktno digitalizirani radio signal – ovako funkcioniše Software Defined Radio.
Otvorit ćemo i obraditi datoteku pomoću Pythona. Oni koji žele eksperimentirati sami mogu preuzeti primjer snimke. .
Učitat ćemo datoteku i pogledati što se unutra nalazi.
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()
Rezultat: vidimo jasne 'impulse' na pozadini šuma.
Svaki 'impuls' zapravo je signal, čija je struktura dobro vidljiva ako povećamo razlučivost na grafu.
Kao što se može vidjeti, slika sasvim odgovara onome što je opisano iznad. Može se krenuti s obradom podataka.
Dekodiranje
Za početak, potrebno je dobiti bitni tok. Sam signal je kodiran pomoću manchester encodinga:
Iz razlike nivoa u polubajtovima lako se mogu dobiti stvarne '0' i '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 += "1"
Struktura samog signala ima sljedeći oblik:
Pogledajmo polja detaljnije.
DF (Downlink Format, 5 bitova) – određuje tip poruke. Ima nekoliko tipova:
()
Nas zanimaju samo tip DF17, jer upravo on sadrži koordinate zrakoplova.
ICAO (24 bita) – međunarodni jedinstveni kod zrakoplova. Provjeriti avion prema njegovom kodu možete (nažalost, autor je prestao ažurirati bazu, ali ona je još uvijek aktualna). Na primjer, za kod 3c5ee2 imamo sljedeće informacije:
Ispravak: u opis koda ICAO je prikazan detaljnije, zainteresiranima preporučujem da se upoznaju.
DATA (56 ili 112 bita) – zapravo podaci koje ćemo dekodirati. Prvih 5 bita podataka – polje Type Code, koje sadrži podtip pohranjenih podataka (ne miješati s DF). Takvih tipova postoji dosta:
()
Razmotrimo nekoliko primjera paketa.
Aircraft identification
Primjer u binarnom obliku:
00100 011 000101 010111 000111 110111 110001 111000
Polja podataka:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, svaki znak C1-C8 sadrži kodove koji odgovaraju indeksima u nizu:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Dekodiranjem niza, lako je dobiti kod aviona: EWG7184
symbols = "#ABCDEFGHIJKLMNOPQRSTUVWXY####_###############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("Identifikacija aviona:", code_str.replace('#', ''))
Položaj u zraku
Ako je naziv jednostavan, s koordinatama je malo složenije. Prenose se u obliku parova, parnih i neparnih okvira. Kod polja TC = 01011b = 11.
Primjer parnog i neparnog paketa:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Samo izračunavanje koordinata odvija se prema prilično složenoj formuli:
()
Nisam stručnjak za GIS, tako da ne znam odakle dolazi. Tko zna, neka napiše u komentarima.
Visina se smatra jednostavnije — ovisno o određenom bitu, može se predstavljati ili kao višekratnik od 25 ili 100 stopa.
Brzina u zraku
Paketi s TC=19. Zanimljivo je da brzina može biti precizna, u odnosu na tlo (Ground Speed), ili zračna, mjereno senzorom aviona (Airspeed). Također se prenosi mnoštvo različitih polja:
()
Zaključak
Kao što se može vidjeti, tehnologija ADS-B postala je zanimljiv spoj, kada jedan standard koristi ne samo profesionalcima, već i običnim korisnicima. No, naravno, ključnu ulogu odigralo je sniženje cijene tehnologije digitalnih SDR-prijemnika, omogućavajući na uređaju gotovo 'za sitniš' primanje signala s frekvencijama iznad gigaherca.
U samom standardu naravno, ima mnogo više od svega. zainteresirani mogu pogledati PDF na stranici ili posjetiti već spomenutu .
Teško da će mnogo toga biti korisno, ali barem bi opća ideja kako to funkcionira, nadam se, ostala.
Usput, gotov dekoder na Pythonu već postoji, može se proučiti . A vlasnici SDR-prijemnika mogu sastaviti i pokrenuti gotov ADS-B dekoder , detaljnije o tome je objašnjeno u .
Izvorni kod parsera, opisan u članku, naveden je ispod. Ovo je testni primjer, ne pretendira na produkciju, ali nešto u njemu funkcioniše, i njime se može parsirati gore navedeni zapisani fajl.
Izvorni kod (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()
Nadam se da je nekome bilo zanimljivo, hvala na pažnji.
Izvor: habr.com
