Përshëndetje Habr. Me siguri, secili që ndonjëherë ka takuar ose ka shoqëruar të afërm apo miq në avion, ka përdorur shërbimin falas Flightradar24. Ky është një mënyrë shumë e përshtatshme për të ndjekur pozitat e avionëve në kohë reale.
Në u përshkrua principi i funksionimit të një shërbimi të tillë online. Tani do të shkojmë më tej dhe do të zbulojmë se cilat të dhëna dërgohen dhe priten nga aeronavigacioni në stacionin e pranimit, dhe do t'i dekodojmë ato vetë duke përdorur Python.
Historia
ĂshtĂ« e qartĂ« se tĂ« dhĂ«nat mbi avionĂ«t nuk dĂ«rgohen pĂ«r tĂ« lejuar pĂ«rdoruesit t'i shohin ato nĂ« telefonat e tyre inteligjentĂ«. Sistemi quhet ADSâB (Automatic dependent surveillanceâbroadcast) dhe shĂ«rben pĂ«r transmetimin automatik tĂ« informacionit rreth aeroplanit nĂ« qendrĂ«n e kontrollit tĂ« trafikut ajror â dĂ«rgohen identifikuesi, 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 numri i avionĂ«ve u rrit.
Në mënyrë teknike, ADS-B përbëhet nga një transmetues në aeroplan që periodikisht dërgon paketa me informacion në një frekuencë mjaft të lartë prej 1090 MHz (ka edhe mënyra të tjera, por ato nuk janë aq interesante për ne, pasi koordinatat dërgohen vetëm këtu). Sigurisht, përveç transmetuesit, ka edhe një marrës diku në aeroport, por për ne, si përdorues, është e rëndësishme marrësi ynë i vet.
Për krahasim, sistemi i parë i tillë, Airnav Radarbox, i projektuar për përdoruesit e zakonshëm, u shfaq në vitin 2007 dhe kishte një çmim prej rreth 900$, ndërsa abonimi për shërbime online kushtonte rreth 250$ në vit.
Mund tĂ« lexoni mendimet e atyre pronarĂ«ve tĂ« parĂ« rusĂ« nĂ« forumin . Tani, kur marrĂ«sit RTL-SDR janĂ« bĂ«rĂ« masivisht tĂ« disponueshĂ«m, njĂ« pajisje e ngjashme mund tĂ« ndĂ«rtohet pĂ«r 30$, mĂ« shumĂ« rreth kĂ«tij procesi Ă«shtĂ« diskutuar nĂ« . Tani do tĂ« kalojmĂ« nĂ« protokollin â tĂ« shohim si funksionon ai.
Pranimi i sinjaleve
Për të filluar, sinjali duhet regjistruar. I gjithë sinjali ka një gjatësi vetëm 120 mikrosekonda, kështu që për të analizuar ngjyrat e tij, është e dëshirueshme një marrës SDR me frekuencë të numërimit prej të paktën 5MHz.
Pas regjistrimit, ne marrim njĂ« skedar WAV me frekuencĂ« numĂ«rimi prej 5000000 samples/sec, 30 sekonda regjistrimi tĂ« tillĂ« peshojnĂ« rreth 500MB. TĂ« dĂ«gjosh atĂ« me njĂ« lojtar media, sigurisht, Ă«shtĂ« e padobishme â skedari pĂ«rmban jo muzikĂ«, por sinjalin radio tĂ« digitalizuar â kjo Ă«shtĂ« mĂ«nyra se si funksionon Software Defined Radio.
Ne do ta hapim dhe përpunojmë skedarin me ndihmën e Python. Ata që dëshirojnë të eksperimentojnë vetë mund të shkarkojnë një shembuj regjistrimi .
Do të ngarkojmë skedarin dhe do 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 "impulsin" e qartë mbi sfondin e zhurmës.
Ădo "impuls" Ă«shtĂ« nĂ« tĂ« vĂ«rtetĂ« sinjali, struktura e tĂ« cilit duket qartĂ« nĂ«se e rritim zgjidhjen nĂ« grafik.
Siç mund të shihet, figura është plotësisht në përputhje me përshkrimin e mësipërm. Mund të fillojmë përpunimin e të dhënave.
Dekodimi
Për të filluar, duhet të marrim rrjedhën e bitëve. Sinjali vetë është koduar me anë të manchester encoding:
Nga ndryshimi i niveleve në gjysmabytetesh është lehtë të merret "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 këtë formë:
Të shohim fushat më në detaje.
DF (Downlink Format, 5 bit) â pĂ«rcakton llojin e mesazhit. Ka disa lloje:
()
Na intereson vetëm lloji DF17, pasi ai përmban koordinatat e aeroplanit.
ICAO (24 bit) â kodi unik ndĂ«rkombĂ«tar i aeroplanit. Mund ta kontrollosh aeroplanin sipas kodit tĂ« tij (pĂ«r fat tĂ« keq, autori ka ndaluar pĂ«rditĂ«simin e bazĂ«s, por ajo akoma Ă«shtĂ« aktuale). PĂ«r shembull, pĂ«r kodin 3c5ee2 kemi informacionin e mĂ«poshtĂ«m:
Korrigjim: në përshkrimi i kodit ICAO është dhënë më në detaje, ata që janë të interesuar rekomandohen që të informohen.
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 Type Code, qĂ« pĂ«rmban naltĂ«sinĂ« e tĂ« dhĂ«nave tĂ« ruhen (mos e ngatĂ«rro me DF). Ka shumĂ« lloje tĂ« tillĂ«:
()
Le të analizojmë disa shembuj paketa.
Identifikimi i aviacionit
Shembulli në formën binar:
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 kodet që korrespondohen me indekset në rresht:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Duke dekoduar rreshtin, nuk është e zakonshme të marrësh kodin e aeroplanit: 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 Aviacionit:", code_str.replace('#', ''))
Pozita ajrore
Nëse emri është i qartë, koordinatat janë më të ndërlikuara. Ato transmetohen në formën e dy framave, çift dhe tek. Kodi i fushës TC = 01011b = 11.
Shembulli i paketave çift dhe tek:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Llogaritja e vetë koordinatave bëhet përmes një formule mjaft të sofistikuar:
()
Nuk jam specialist në GIS, prandaj nuk e di se si nxirret. Kush e di, le të shkruajë në komentet.
LartĂ«sia llogaritet mĂ« thjesht â nĂ« varĂ«si tĂ« njĂ« bita tĂ« caktuar, ajo mund tĂ« paraqitet si e shumĂ«fishueshme me 25 ose 100 kĂ«mbĂ«.
Shpejtësia ajrore
Paketa me TC=19. E veçanta është se shpejtësia mund të jetë e saktë, në lidhje me tokën (Shpejtësia në Tokë), ashtu si edhe ajo ajrore, e cila matet nga sensori i avionit (Shpejtësia Ajrore). Po ashtu, transmetohet një sërë fushash të ndryshme:
()
PĂ«rfundimi
Siç mund të shihet, teknologjia ADS-B është bërë një simbiozë interesante, ku një standard shërben jo vetëm profesionistëve, por edhe përdoruesve të zakonshëm. Sigurisht, roli kyç në këtë e ka luajtur ulja e kostove të teknologjisë së marrësve SDR digjitalë, të cilat lejojnë që në pajisje të pranohet sinjale me frekuenca mbi gigaherc.
Në standardin e vet ka shumë më tepër. Ata që janë të interesuar mund të shikojnë PDF-në në faqen ose të vizitojnë tashmë të përmendurin më sipër .
Nuk besoj se do t'i nevojitet shumëkujt gjithçka e mësipërme, por të paktën shpresa është se ideja e përgjithshme se si funksionon ka mbetur.
Apo, ekziston një dekoder i gatshëm në Python, të cilin mund ta studiojnë . Ndërsa pronarët e marrësve SDR mund të mbledhin dhe nisin një dekoder të gatshëm ADS-B , për më shumë, është treguar në .
Kodi burimor i parser-it, i përshkruar në artikull, është dhënë poshtë. Ky është një shembull testimi, që nuk pretendon të jetë prodhim, por disa gjëra funksionojnë në të, dhe me të mund të parse një skedë të 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 dikujt i ka pëlqyer, faleminderit për vëmendjen.
Burimi: habr.com
