Saluton Habr. Verŝajne ĉiuj, kiuj iam renkontis aŭ vidis parencojn aŭ amikojn en aviadilo, uzis la senpagan servon Flightradar24. Ĉi tio estas tre oportuna maniero spuri la pozicion de la aviadilo en reala tempo.
В La funkcia principo de tia reta servo estis priskribita. Ni nun daŭrigos kaj eltrovos kiajn datumojn oni sendas kaj ricevas de la aviadilo al la riceva stacio kaj mem malkodos ĝin uzante Python.
История
Evidente, aviadilaj datumoj ne estas transdonitaj por ke uzantoj vidu sur siaj inteligentaj telefonoj. La sistemo nomiĝas ADS-B (Aŭtomata dependa gvatado—elsendo), kaj estas uzata por aŭtomate transdoni informojn pri la aviadilo al la kontrolcentro - ĝia identigilo, koordinatoj, direkto, rapideco, alteco kaj aliaj datumoj estas transdonitaj. Antaŭe, antaŭ la apero de tiaj sistemoj, la sendanto povis vidi nur punkton sur la radaro. Ĉi tio ne plu sufiĉis, kiam estis tro da aviadiloj.
Teknike, ADS-B konsistas el dissendilo sur aviadilo, kiu periode sendas pakaĵojn da informoj kun sufiĉe alta ofteco de 1090 MHz (ekzistas aliaj reĝimoj, sed ni ne tiom interesiĝas pri ili, ĉar koordinatoj estas transdonitaj nur ĉi tie). Kompreneble, krom la dissendilo, estas ankaŭ ricevilo ie en la flughaveno, sed por ni, kiel uzantoj, nia propra ricevilo estas interesa.
Cetere, por komparo, la unua tia sistemo, Airnav Radarbox, desegnita por ordinaraj uzantoj, aperis en 2007, kaj kostis ĉirkaŭ 900 USD; abono al retservoj kostis pliajn 250 USD jare.
Recenzoj pri tiuj unuaj rusaj posedantoj legeblas sur la forumo . Nun kiam RTL-SDR-riceviloj fariĝis vaste haveblaj, simila aparato povas esti kunvenita kontraŭ $ 30; pli pri tio estis en . Ni transiru al la protokolo mem - ni vidu kiel ĝi funkcias.
Ricevante signalojn
Unue, la signalo devas esti registrita. La tuta signalo havas daŭron de nur 120 mikrosekundoj, do por komforte malmunti ĝiajn komponantojn, SDR-ricevilo kun specimena ofteco de almenaŭ 5 MHz estas dezirinda.
Post registrado, ni ricevas WAV-dosieron kun specimena rapideco de 5000000 specimenoj/se; 30 sekundoj de tia registrado "pezas" ĉirkaŭ 500MB. Aŭskulti ĝin per plurmedia ludilo, kompreneble, estas senutila - la dosiero ne enhavas sonon, sed rekte ciferecigitan radiosignalon - ĝuste tiel funkcias Software Defined Radio.
Ni malfermos kaj prilaboros la dosieron per Python. Tiuj, kiuj volas eksperimenti memstare, povas elŝuti ekzemplon de registrado .
Ni elŝutu la dosieron kaj vidu kio estas ene.
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()
Rezulto: ni vidas evidentajn "pulsojn" kontraŭ la fona bruo.
Ĉiu "pulso" estas signalo, kies strukturo estas klare videbla se vi pliigas la rezolucion sur la grafikaĵo.
Kiel vi povas vidi, la bildo estas sufiĉe kongrua kun tio, kio estas donita en la supra priskribo. Vi povas komenci prilabori la datumojn.
Malkodado
Unue, vi devas ricevi iom da fluo. La signalo mem estas ĉifrita uzante Manĉestrokodigon:
De la niveldiferenco en mordetoj estas facile akiri verajn "0" kaj "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"
La strukturo de la signalo mem estas kiel sekvas:
Ni rigardu la kampojn pli detale.
DF (Downlink Format, 5 bitoj) - determinas la specon de mesaĝo. Estas pluraj tipoj:
()
Ni interesiĝas nur pri tipo DF17, ĉar... Ĝuste ĉi tio enhavas la koordinatojn de la aviadilo.
ICAO (24 bitoj) - internacia unika kodo de la aviadilo. Vi povas kontroli la aviadilon per ĝia kodo (bedaŭrinde, la aŭtoro ĉesis ĝisdatigi la datumbazon, sed ĝi ankoraŭ estas grava). Ekzemple, por kodo 3c5ee2 ni havas la jenajn informojn:
Redakto: en La priskribo de la ICAO-kodo estas donita pli detale; mi rekomendas, ke la interesatoj legu ĝin.
DATUMOJ (56 aŭ 112 bitoj) - la realaj datumoj, kiujn ni malkodos. La unuaj 5 bitoj da datumoj estas la kampo Tajpu Kodon, enhavante la subspecon de la datenoj stokitaj (malsama al DF). Estas sufiĉe multaj el ĉi tiuj tipoj:
()
Ni rigardu kelkajn ekzemplojn de pakoj.
Identigo de aviadiloj
Ekzemplo en binara formo:
00100 011 000101 010111 000111 110111 110001 111000
Datumaj kampoj:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, ĉiu signo C1-C8 enhavas kodojn respondajn al indeksoj en la linio:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Malkodante la ŝnuron, estas facile akiri la aviadilkodon: 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('#', ''))
Aera pozicio
Se la nomo estas simpla, tiam la koordinatoj estas pli komplikaj. Ili estas transdonitaj en formo de 2, paraj kaj neparaj kadroj. Kampokodo TC = 01011b = 11.
Ekzemplo de paraj kaj neparaj pakaĵoj:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
La kalkulo de koordinatoj mem okazas laŭ sufiĉe malfacila formulo:
()
Mi ne estas spertulo pri GIS, do mi ne scias de kie ĝi venas. Kiu scias, skribu en la komentoj.
Alteco estas konsiderita pli simpla - depende de la specifa bito, ĝi povas esti reprezentita kiel aŭ oblo de 25 aŭ 100 futoj.
Aera Rapideco
Pako kun TC=19. La interesa afero ĉi tie estas, ke la rapideco povas esti aŭ preciza, relative al la grundo (Ground Speed), aŭ aera, mezurita per aviadila sensilo (Airspeed). Multaj malsamaj kampoj ankaŭ estas elsenditaj:
()
konkludo
Kiel vi povas vidi, ADS-B-teknologio fariĝis interesa simbiozo, kiam normo utilas ne nur por profesiuloj, sed ankaŭ por ordinaraj uzantoj. Sed kompreneble, ŝlosila rolo en ĉi tio ludis la pli malmultekosta teknologio de ciferecaj SDR-riceviloj, kiu permesas al la aparato laŭvorte ricevi signalojn kun frekvencoj super gigaherco "por pencoj".
En la normo mem, kompreneble, estas multe pli. Interesitoj povas vidi la PDF sur la paĝo aŭ vizitu tiun jam supre menciitan .
Estas neverŝajne, ke ĉio ĉi-supra estos utila al multaj, sed almenaŭ la ĝenerala ideo pri kiel ĝi funkcias, mi esperas, restas.
Cetere, preta malĉifrilo en Python jam ekzistas, vi povas studi ĝin . Kaj posedantoj de SDR-riceviloj povas kunveni kaj lanĉi pretan ADS-B-malĉifrilon , ĉi tio estis diskutita pli detale en .
La fontkodo de la analizilo priskribita en la artikolo estas donita sub la tranĉo. Ĉi tio estas prova ekzemplo, kiu ne ŝajnigas esti produktado, sed iuj aferoj funkcias en ĝi, kaj ĝi povas esti uzata por analizi la dosieron registritan supre.
Fontkodo (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()
Mi esperas, ke iu interesiĝis, dankon pro via atento.
fonto: www.habr.com