Përshëndetje Habr. Ndoshta të gjithë ata që kanë takuar ose larguar ndonjëherë të afërmit ose miqtë në një avion kanë përdorur shërbimin falas Flightradar24. Kjo është një mënyrë shumë e përshtatshme për të gjurmuar pozicionin e avionit në kohë reale.
В Parimi i funksionimit të një shërbimi të tillë në internet u përshkrua. Tani do të vazhdojmë dhe do të kuptojmë se çfarë të dhënash dërgohen dhe merren nga avioni në stacionin marrës dhe do t'i deshifrojmë vetë duke përdorur Python.
Histori
Natyrisht, të dhënat e avionit nuk transmetohen që përdoruesit t'i shohin në telefonat e tyre inteligjentë. Sistemi quhet ADS-B (Automatic dependent surveillance-transmetim) dhe përdoret për të transmetuar automatikisht informacionin rreth avionit në qendrën e kontrollit - identifikuesi i tij, koordinatat, drejtimi, shpejtësia, lartësia dhe të dhëna të tjera transmetohen. Më parë, para ardhjes së sistemeve të tilla, dispeçeri mund të shihte vetëm një pikë në radar. Kjo nuk mjaftonte më kur kishte shumë avionë.
Teknikisht, ADS-B përbëhet nga një transmetues në një avion që dërgon periodikisht pako informacioni në një frekuencë mjaft të lartë prej 1090 MHz (ka mënyra të tjera, por ne nuk jemi aq të interesuar për to, pasi koordinatat transmetohen vetëm këtu). Sigurisht, përveç transmetuesit, ka edhe një marrës diku në aeroport, por për ne, si përdorues, është interesant marrësi ynë.
Nga rruga, për krahasim, sistemi i parë i tillë, Airnav Radarbox, i krijuar për përdoruesit e zakonshëm, u shfaq në vitin 2007 dhe kushtoi rreth 900 dollarë; një abonim në shërbimet e rrjetit kushtonte 250 dollarë të tjerë në vit.
Shqyrtimet e atyre pronarëve të parë rusë mund të lexohen në forum . Tani që marrësit RTL-SDR janë bërë gjerësisht të disponueshëm, një pajisje e ngjashme mund të montohet për 30 dollarë; më shumë rreth kësaj ishte në . Le të kalojmë në vetë protokollin - le të shohim se si funksionon.
Marrja e sinjaleve
Së pari, sinjali duhet të regjistrohet. I gjithë sinjali ka një kohëzgjatje prej vetëm 120 mikrosekonda, kështu që për të çmontuar me lehtësi përbërësit e tij, është i dëshirueshëm një marrës SDR me një frekuencë kampionimi prej të paktën 5 MHz.
Pas regjistrimit, ne marrim një skedar WAV me një shpejtësi kampionimi prej 5000000 mostra/sek; 30 sekonda e një regjistrimi të tillë "peshojnë" rreth 500 MB. Dëgjimi i tij me një media player, natyrisht, është i padobishëm - skedari nuk përmban zë, por një sinjal radio të dixhitalizuar drejtpërdrejt - kjo është pikërisht mënyra se si funksionon Radio Defined Software.
Ne do të hapim dhe përpunojmë skedarin duke përdorur Python. Ata që duan të eksperimentojnë vetë mund të shkarkojnë një shembull regjistrimi .
Le të shkarkojmë skedarin dhe të shohim se çfarë ka brenda.
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()
Rezultati: ne shohim "pulse" të dukshme kundrejt zhurmës së sfondit.
Çdo "puls" është një sinjal, struktura e të cilit është qartë e dukshme nëse rritni rezolucionin në grafik.
Siç mund ta shihni, fotografia është mjaft në përputhje me atë që është dhënë në përshkrimin e mësipërm. Ju mund të filloni përpunimin e të dhënave.
Dekodimi
Së pari, ju duhet të merrni pak transmetim. Vetë sinjali është i koduar duke përdorur kodimin Manchester:
Nga diferenca e nivelit në thithka është e lehtë të marrësh "0" dhe "1" të vërtetë.
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"
Struktura e vetë sinjalit është si më poshtë:
Le t'i shikojmë fushat në më shumë detaje.
DF (Formati i lidhjes zbritëse, 5 bit) - përcakton llojin e mesazhit. Ka disa lloje:
()
Na intereson vetëm tipi DF17, sepse... Është kjo që përmban koordinatat e avionit.
ICAO (24 bit) - kodi unik ndërkombëtar i avionit. Ju mund ta kontrolloni aeroplanin me kodin e tij (për fat të keq, autori ka ndaluar përditësimin e bazës së të dhënave, por ajo është ende e rëndësishme). Për shembull, për kodin 3c5ee2 kemi informacionin e mëposhtëm:
Edit: në Përshkrimi i kodit të ICAO është dhënë më në detaje, të interesuarve u rekomandoj ta lexojnë atë.
TË DHËNAT (56 ose 112 bit) - të dhënat aktuale që do të deshifrojmë. 5 bitet e para te te dhenave jane fusha Lloji i Kodit, që përmban nëntipin e të dhënave që ruhen (të mos ngatërrohet me DF). Ka mjaft nga këto lloje:
()
Le të shohim disa shembuj të paketave.
Identifikimi i avionit
Shembull në formë binare:
00100 011 000101 010111 000111 110111 110001 111000
Fushat e të dhënave:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, çdo karakter C1-C8 përmban kode që korrespondojnë me indekset në rresht:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789#######
Duke deshifruar vargun, është e lehtë të merret kodi i avionit: 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('#', ''))
Pozicioni në ajër
Nëse emri është i thjeshtë, atëherë koordinatat janë më të ndërlikuara. Ato transmetohen në formën e 2 kornizave çift dhe tek. Kodi i fushës TC = 01011b = 11.
Shembull i paketave çift dhe tek:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Llogaritja e koordinatave në vetvete ndodh sipas një formule mjaft të ndërlikuar:
()
Unë nuk jam ekspert i GIS, kështu që nuk e di nga vjen. Kush e di, shkruani në komente.
Lartësia konsiderohet më e thjeshtë - në varësi të bitit specifik, mund të përfaqësohet si shumëfish i 25 ose 100 këmbëve.
Shpejtësia ajrore
Paketa me TC=19. Gjëja interesante këtu është se shpejtësia mund të jetë ose e saktë, në raport me tokën (Ground Speed), ose ajrore, e matur nga një sensor avioni (Airspeed). Shumë fusha të ndryshme transmetohen gjithashtu:
()
Përfundim
Siç mund ta shihni, teknologjia ADS-B është bërë një simbiozë interesante, kur një standard është i dobishëm jo vetëm për profesionistët, por edhe për përdoruesit e zakonshëm. Por sigurisht, një rol kyç në këtë luajti teknologjia më e lirë e marrësve dixhital SDR, e cila lejon pajisjen të marrë fjalë për fjalë sinjale me frekuenca mbi një gigahertz "për qindarkë".
Në vetë standardin, natyrisht, ka shumë më tepër. Të interesuarit mund ta shikojnë PDF-në në faqe ose vizitoni atë të përmendur më lart .
Nuk ka gjasa që të gjitha sa më sipër do të jenë të dobishme për shumë njerëz, por të paktën ideja e përgjithshme se si funksionon, shpresoj, mbetet.
Nga rruga, një dekoder i gatshëm në Python tashmë ekziston, ju mund ta studioni atë . Dhe pronarët e marrësve SDR mund të mbledhin dhe lëshojnë një dekoder të gatshëm ADS-B , kjo u diskutua më në detaje në .
Kodi burimor i analizuesit të përshkruar në artikull është dhënë më poshtë. Ky është një shembull provë që nuk pretendon të jetë prodhim, por disa gjëra funksionojnë në të dhe mund të përdoret për të analizuar skedarin e regjistruar më sipër.
Kodi burimor (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()
Shpresoj se dikush ishte i interesuar, faleminderit për vëmendjen tuaj.
Burimi: www.habr.com