Salam Habr. Yəqin ki, qohumları və ya dostları ilə təyyarədə görüşən və ya yola salan hər kəs pulsuz Flightradar24 xidmətindən istifadə edib. Bu, real vaxt rejimində təyyarənin mövqeyini izləmək üçün çox rahat üsuldur.
В Belə bir onlayn xidmətin iş prinsipi təsvir edilmişdir. İndi biz davam edəcəyik və təyyarədən qəbuledici stansiyaya hansı məlumatların göndərildiyini və qəbul edildiyini anlayacağıq və Python istifadə edərək özümüz deşifrə edəcəyik.
Hekayə
Aydındır ki, təyyarə məlumatları istifadəçilərin smartfonlarında görməsi üçün ötürülmür. Sistem ADS-B (Avtomatik asılı müşahidə—yayım) adlanır və hava gəmisi haqqında məlumatı idarəetmə mərkəzinə avtomatik ötürmək üçün istifadə olunur - onun identifikatoru, koordinatları, istiqaməti, sürəti, hündürlüyü və digər məlumatlar ötürülür. Əvvəllər, belə sistemlərin meydana çıxmasından əvvəl dispetçer yalnız radarda bir nöqtə görə bilirdi. Həddindən artıq çox təyyarə olanda bu artıq kifayət etmirdi.
Texniki cəhətdən, ADS-B, vaxtaşırı 1090 MHz kifayət qədər yüksək tezlikdə məlumat paketlərini göndərən bir təyyarədəki ötürücüdən ibarətdir (başqa rejimlər var, lakin biz onlarla o qədər də maraqlanmırıq, çünki koordinatlar yalnız burada ötürülür). Təbii ki, ötürücüdən əlavə, hava limanında hardasa qəbuledici də var, amma istifadəçilər olaraq bizim üçün öz qəbuledicimiz maraqlıdır.
Yeri gəlmişkən, müqayisə üçün qeyd edək ki, adi istifadəçilər üçün nəzərdə tutulmuş ilk belə sistem Airnav Radarbox 2007-ci ildə yaranıb və təxminən 900 dollara başa gəlir, şəbəkə xidmətlərinə abunə isə ildə daha 250 dollara başa gəlir.
Həmin ilk rus sahiblərinin rəylərini forumda oxuya bilərsiniz . İndi RTL-SDR qəbulediciləri geniş şəkildə əlçatan olub, oxşar cihazı 30 dollara yığmaq olar; bu barədə ətraflı . Gəlin protokolun özünə keçək - onun necə işlədiyini görək.
Siqnalların qəbulu
Əvvəlcə siqnal qeyd edilməlidir. Bütün siqnalın müddəti cəmi 120 mikrosaniyədir, ona görə də onun komponentlərini rahat şəkildə sökmək üçün ən azı 5 MHz seçmə tezliyi olan SDR qəbuledicisi arzu edilir.
Qeydiyyatdan sonra biz 5000000 nümunə/san seçmə sürəti ilə WAV faylı alırıq; belə qeydin 30 saniyəsi təxminən 500MB "çəkisi". Media pleyeri ilə onu dinləmək, əlbəttə ki, faydasızdır - faylda səs yoxdur, ancaq birbaşa rəqəmsal radio siqnalı var - Software Defined Radio məhz belə işləyir.
Python istifadə edərək faylı açıb emal edəcəyik. Öz başına təcrübə etmək istəyənlər nümunə qeydini yükləyə bilərlər .
Gəlin faylı yükləyək və içində nə olduğunu görək.
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()
Nəticə: fon səs-küyünə qarşı açıq-aşkar “nəbzlər” görürük.
Hər bir "nəbz" bir siqnaldır, qrafikdəki qətnaməni artırsanız, quruluşu aydın görünür.
Gördüyünüz kimi, şəkil yuxarıdakı təsvirdə verilənə kifayət qədər uyğundur. Məlumatları emal etməyə başlaya bilərsiniz.
Deşifrə
Əvvəlcə bir az axın əldə etməlisiniz. Siqnalın özü Mançester kodlaşdırmasından istifadə edərək kodlanır:
Nibblesdəki səviyyə fərqindən həqiqi "0" və "1" almaq asandır.
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"
Siqnalın özünün quruluşu belədir:
Sahələrə daha ətraflı baxaq.
DF (Downlink Format, 5 bits) - mesajın növünü müəyyən edir. Bir neçə növ var:
()
Bizi yalnız DF17 növü maraqlandırır, çünki... Təyyarənin koordinatlarını ehtiva edən budur.
ICAO (24 bit) - təyyarənin beynəlxalq unikal kodu. Təyyarəni kodu ilə yoxlaya bilərsiniz (təəssüf ki, müəllif məlumat bazasını yeniləməyi dayandırdı, lakin hələ də aktualdır). Məsələn, 3c5ee2 kodu üçün bizdə aşağıdakı məlumatlar var:
Redaktə: in ICAO kodunun təsviri daha ətraflı verilmişdir, maraqlananlara onu oxumağı tövsiyə edirəm.
MƏLUMAT (56 və ya 112 bit) - deşifrə edəcəyimiz faktiki məlumatlar. Məlumatın ilk 5 biti sahədir Kod yazın, saxlanılan məlumatların alt növünü ehtiva edir (DF ilə qarışdırılmamalıdır). Bu növlərdən bir neçəsi var:
()
Paketlərin bir neçə nümunəsinə baxaq.
Təyyarə identifikasiyası
İkili formada nümunə:
00100 011 000101 010111 000111 110111 110001 111000
Məlumat sahələri:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, hər C1-C8 simvolunda sətirdəki indekslərə uyğun kodlar var:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Simli deşifrə etməklə, təyyarə kodunu əldə etmək asandır: 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('#', ''))
Hava mövqeyi
Əgər ad sadədirsə, deməli koordinatlar daha mürəkkəbdir. Onlar 2, cüt və tək kadrlar şəklində ötürülür. Sahə kodu TC = 01011b = 11.
Cüt və tək paketlərə misal:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Koordinatların hesablanması özü olduqca çətin bir düstura görə baş verir:
()
Mən GIS mütəxəssisi deyiləm, ona görə də onun haradan gəldiyini bilmirəm. Kim bilir, şərhlərdə yazın.
Hündürlük daha sadə hesab olunur - xüsusi bitdən asılı olaraq, ya 25, ya da 100 futun qatı kimi təmsil oluna bilər.
Hava sürəti
TC=19 ilə paket. Burada maraqlı olan odur ki, sürət ya dəqiq, yerə nisbətən (Yer Sürəti) ola bilər, ya da təyyarə sensoru (Airspeed) ilə ölçülən havada ola bilər. Bir çox müxtəlif sahələr də ötürülür:
()
Nəticə
Gördüyünüz kimi, standart yalnız peşəkarlar üçün deyil, həm də adi istifadəçilər üçün faydalı olduqda, ADS-B texnologiyası maraqlı bir simbioza çevrildi. Ancaq əlbəttə ki, bunda əsas rolu rəqəmsal SDR qəbuledicilərinin daha ucuz texnologiyası oynadı ki, bu da cihaza "quruşa görə" gigahertzdən yuxarı tezlikli siqnalları sözün əsl mənasında qəbul etməyə imkan verir.
Standartın özündə, əlbəttə ki, daha çox şey var. Maraqlananlar PDF-ə səhifədən baxa bilərlər və ya yuxarıda qeyd olunan birinə baş çəkin .
Yuxarıda göstərilənlərin hamısının çoxları üçün faydalı olacağı ehtimalı azdır, amma ən azı bunun necə işlədiyinə dair ümumi fikir, ümid edirəm ki, qalır.
Yeri gəlmişkən, Python-da hazır dekoder artıq mövcuddur, onu öyrənə bilərsiniz . SDR qəbuledicilərinin sahibləri isə hazır ADS-B dekoderini yığıb işə sala bilərlər , bu barədə daha ətraflı bəhs edildi .
Məqalədə təsvir edilən təhlilçinin mənbə kodu kəsilmənin altında verilmişdir. Bu, istehsal olduğunu iddia etməyən bir sınaq nümunəsidir, lakin bəzi şeylər orada işləyir və yuxarıda qeyd olunan faylı təhlil etmək üçün istifadə edilə bilər.
Mənbə kodu (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()
Ümid edirəm kimsə maraqlandı, diqqətinizə görə təşəkkür edirəm.
Mənbə: www.habr.com