Përshëndetje Habr. Ndoshta të gjithë ata që kanë takuar ose shoqëruar të afërm ose miq në avion e kanë përdorur shërbimin falas Flightradar24. Kjo është një mënyrë shumë e përshtatshme për të ndjekur pozitat e avionit në kohë reale.
Në u përshkrua principi i funksionimit të një shërbimi të tillë online. Tani do të shkojmë më tutje dhe do të zbulojmë cilat të dhëna dërgohen dhe pranohen nga aeroplani në stacionin e marrjes, dhe do t'i dekodojmë ato vetë me Python.
Historia
Është e qartë se të dhënat rreth avionëve nuk dërgohen për t'u parë nga përdoruesit në telefonat e tyre të mençur. Sistemi quhet ADS–B (Automatic dependent surveillance—broadcast), dhe shërben për transmetimin automatike të informacionit rreth aeroplanit në qendrën e kontrollit — dërgohen identifikatori, koordinatat, drejtimi, shpejtësia, lartësia dhe të dhëna të tjera. Më parë, para se të shfaqeshin këto sisteme, dispeçeri mund të shihte vetëm një pikë në radar. Kjo u bë e pamjaftueshme kur avionët u bënë shumë më të shumtë.
Teknikisht, ADS-B përbëhet nga një transmetues në aeroplan, i cili dërgon rregullisht paketa me informacion në një frekuencë mjaft të lartë 1090 MHz (ka edhe mënyra të tjera, por ato nuk na interesojnë shumë, pasi koordinatat dërgohen vetëm këtu). Natyrisht, përveç transmetuesit, ka edhe një marrës diku në aeroport, por për ne, si përdorues, marrësi ynë është i rëndësishëm.
Përveç kësaj, për krahasim, sistemi i parë i tillë, Airnav Radarbox, i destinuar për përdorues të zakonshëm, u shfaq në vitin 2007, dhe kushtonte rreth 900$, ndërsa abone vitore për shërbimet në rrjet kushtonte rreth 250$.
Reviews of those first Russian owners can be read on the forum . Tani, kur marrësit RTL-SDR janë bërë masivisht të disponueshëm, një pajisje e ngjashme mund të ndihet me 30$, më shumë rreth këtij është folur në . Ne do të kalojmë tani te protokolli — do të shohim si funksionon.
Pranimi i sinjaleve
Për fillim, sinjali duhet të regjistrohet. I gjithë sinjali ka një gjatësi prej vetëm 120 mikrosekondash, kështu që për të analizuar rehat komponentët e tij, është e preferueshme një marrëveshje SDR me frekuencë mostrimi të paktën 5 MHz.
Pas regjistrimit, ne marrim një skedar WAV me frekuencë mostrimi 5000000 mostra/seck, 30 sekonda të tillë regjistrimi peshojnë rreth 500Mb. Të dëgjohet me një mediaplayer sigurisht që është e pavlerë - skedari përmban jo tingull, por sinjalin e digitalizuar të radios - kështu funksionon Software Defined Radio.
Ne do të hapim dhe përpunojmë skedarin me ndihmën e Python. Ata që duan të eksperimentojnë vetë, mund të shkarkojnë shembullin e regjistrimit. .
Le të ngarkojmë skedarin dhe të shohim ç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 "impulse" të qarta në sfondin e zhurmës.
Çdo "impuls" është vetë sinjali, struktura e të cilit duket mirë nëse rritim zgjidhjen në grafik.
Siç mund të shihet, imazhi përputhet plotësisht me përshkrimin e dhënë më lart. Mund të fillojmë përpunimin e të dhënave.
Dekodimi
Për fillim, duhen marrë rrjedhat e bitëve. Sinjali vetë është koduar me encoding manchester:
Nga ndryshimi i niveleve në polubaite është lehtë të marrim "0" dhe "1" reale.
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 e vetë sinjalit ka pamjen e mëposhtme:
Le të shqyrtojmë fushat më në detaje.
DF (Downlink Format, 5 bit) - përcakton tipin e mesazhit. Ka disa tipe:
()
Na intereson vetëm tipi DF17, pasi ai përmban koordinatat e aeroplanit.
ICAO (24 bit) - kodi unik ndërkombëtar i aeroplanit. Mund ta kontrolloni aeroplanin sipas kodit të tij (për fat të keq, autori ndaloi përditësimin e databazës, por ajo është ende relevante). Për shembull, për kodin 3c5ee2 kemi informacionin e mëposhtëm:
Korrigjim: në përshkrimi i kodit ICAO është përcaktuar më në detaje, ata që janë të interesuar rekomandoj ta lexojnë.
DATA (56 ose 112 bit) - në të vërtetë të dhënat që ne do të dekodojmë. 5 bitët e parë të të dhënave janë fusha Tipi i Kodit, duke përmbajtur nën-tipin e të dhënave të ruajtura (mos e ngatërroni me DF). Të tilla janë mjaft tipa:
()
Le të analizojmë disa shembuj të paketave.
Identifikimi i aeroplanit
Shembulli në formën 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 simbol C1-C8 përmban kode që korrespondon me indekset në varg:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Duke dekoduar vargun, nuk është e vështirë të merrni kodin e 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("Identifikimi i Avionit:", code_str.replace('#', ''))
Pozita ajrore
Nëse me emrin gjithçka është e qartë, me koordinatat është më e komplikuar. Ato transmetohen në formën e dy kornizave, çift dhe tek. Kodi i fushës TC = 01011b = 11.
Shembulli i paketave çifte dhe tek:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Vetë llogaritja e koordinatave ndodh sipas një formule mjaft të ndërlikuar:
()
Nuk jam specialist për GIS, kështu që nuk e di nga vjen, kush është i informuar, shkruani në komentet.
Shkalla llogaritet më thjeshtë — varësisht nga një bit i caktuar, ajo mund të paraqitet ose si shumëfish i 25, ose 100 këmbësh.
Shpejtësia ajrore
Paketë me TC=19. Interesante këtu është se shpejtësia mund të jetë e saktë, në lidhje me tokën (Shpejtësia e Tokës), ose ajrore, e cila matet nga ndjesi i avionit (Shpejtësia Ajrore). Po ashtu, transmetohen shumë fushat e ndryshme:
()
Përfundim
Siç mund të shihet, teknologjia ADS-B është bërë një simbiozë interesante, kur ndonjë standard i përdor përveç profesionistëve, edhe përdorues të zakonshëm. Por, sigurisht, roli kryesor këtu e luajti ulja e kostos së teknologjisë së marrësve SDR digjitalë, të cilët lejojnë marrjen e sinjaleve me frekuenca më të larta se gigaherc.
Në standard, sigurisht, ka shumë më tepër. Të interesuarit mund të shohin PDF-në në faqen ose të vizitojnë faqen e përmendur më lart .
Më vjen keq, ndoshta nuk do t'i duhet shumicës së njerëzve gjithçka e shkruar më sipër, por të paktën shpresa është që ideja e përgjithshme e asaj se si funksionon të ketë mbetur.
Për më tepër, një dekoder i gatshëm në Python tashmë ekziston, mund ta studioni . Ndërsa pronarët e marrësve SDR mund të mbledhin dhe shkarkojnë një dekoder të gatshëm ADS-B , më shumë për këtë është folur në .
Kodi burimor i parserit, i përshkruar në artikull, është sjellë poshtë. Ky është një shembull testues, nuk pretendon për prodhim, por diçka në të funksionon, dhe mund të përdoret për të parse skedarin e lartpërmendur.
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 t'i ketë interesuar dikë, faleminderit për vëmendjen.
Burimi: habr.com
