Salom Xabr. Qarindoshlari yoki do'stlarini samolyotda uchratgan yoki kutib olgan har bir kishi bepul Flightradar24 xizmatidan foydalangandir. Bu real vaqt rejimida samolyotning holatini kuzatishning juda qulay usuli.
В Bunday onlayn xizmatning ishlash printsipi tasvirlangan. Endi biz davom etamiz va samolyotdan qabul qiluvchi stantsiyaga qanday ma'lumotlar yuborilayotganini va qabul qilinishini aniqlaymiz va Python yordamida o'zimiz dekodlaymiz.
История
Shubhasiz, samolyot ma'lumotlari foydalanuvchilarning smartfonlarida ko'rishlari uchun uzatilmaydi. Tizim ADS-B (Avtomatik qaram kuzatuv — translyatsiya) deb nomlanadi va havo kemasi haqidagi maʼlumotlarni avtomatik ravishda boshqaruv markaziga uzatish uchun ishlatiladi – uning identifikatori, koordinatalari, yoʻnalishi, tezligi, balandligi va boshqa maʼlumotlar uzatiladi. Ilgari, bunday tizimlar paydo bo'lishidan oldin, dispetcher faqat radardagi nuqtani ko'rishi mumkin edi. Samolyotlar juda ko'p bo'lganda, bu etarli emas edi.
Texnik jihatdan, ADS-B vaqti-vaqti bilan 1090 MGts chastotada ma'lumot paketlarini yuboradigan samolyotdagi uzatgichdan iborat (boshqa rejimlar mavjud, ammo biz ularni unchalik qiziqtirmaymiz, chunki koordinatalar faqat shu erda uzatiladi). Albatta, transmitterdan tashqari, aeroportning biron bir joyida qabul qiluvchi ham mavjud, ammo biz, foydalanuvchilar sifatida, o'zimizning qabul qiluvchimiz qiziq.
Aytgancha, taqqoslash uchun, oddiy foydalanuvchilar uchun mo'ljallangan birinchi Airnav Radarbox tizimi 2007 yilda paydo bo'lgan va taxminan 900 dollar turadi, tarmoq xizmatlariga obuna yiliga yana 250 dollar turadi.
Ushbu birinchi rus egalarining sharhlarini forumda o'qish mumkin . Endi RTL-SDR qabul qiluvchilar keng tarqalgan bo'lib, shunga o'xshash qurilmani 30 dollarga yig'ish mumkin; bu haqda batafsil ma'lumot . Keling, protokolning o'ziga o'taylik - uning qanday ishlashini ko'rib chiqaylik.
Signallarni qabul qilish
Birinchidan, signalni yozib olish kerak. Butun signalning davomiyligi atigi 120 mikrosekundni tashkil qiladi, shuning uchun uning tarkibiy qismlarini qulay tarzda qismlarga ajratish uchun kamida 5 MGts namuna olish chastotasiga ega SDR qabul qiluvchisi maqsadga muvofiqdir.
Yozib olgandan so'ng, biz sekundiga 5000000 30 500 namuna olish tezligiga ega WAV faylini olamiz; bunday yozuvning XNUMX soniyasi taxminan XNUMX MB og'irlikda. Media pleer bilan uni tinglash, albatta, foydasiz - faylda ovoz yo'q, lekin to'g'ridan-to'g'ri raqamlangan radio signal - Software Defined Radio aynan shunday ishlaydi.
Python yordamida faylni ochamiz va qayta ishlaymiz. O'z-o'zidan tajriba o'tkazmoqchi bo'lganlar namunali yozuvni yuklab olishlari mumkin .
Keling, faylni yuklab olamiz va ichida nima borligini ko'ramiz.
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()
Natija: biz fon shovqiniga qarshi aniq "impulslarni" ko'ramiz.
Har bir "puls" - bu signal bo'lib, uning tuzilishi, agar siz grafikdagi piksellar sonini oshirsangiz, aniq ko'rinadi.
Ko'rib turganingizdek, rasm yuqoridagi tavsifda keltirilgan narsalarga juda mos keladi. Siz ma'lumotlarni qayta ishlashni boshlashingiz mumkin.
Dekodlash
Birinchidan, siz bir oz oqim olishingiz kerak. Signalning o'zi Manchester kodlash yordamida kodlangan:
Nibbles darajasidagi farqdan haqiqiy "0" va "1" ni olish oson.
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"
Signalning tuzilishi quyidagicha:
Keling, maydonlarni batafsil ko'rib chiqaylik.
DF (Downlink Format, 5 bit) - xabar turini aniqlaydi. Bir nechta turlari mavjud:
()
Bizni faqat DF17 turi qiziqtiradi, chunki... Bu samolyotning koordinatalarini o'z ichiga oladi.
ICAO (24 bit) - samolyotning xalqaro noyob kodi. Samolyotni uning kodi bo'yicha tekshirishingiz mumkin (afsuski, muallif ma'lumotlar bazasini yangilashni to'xtatdi, lekin u hali ham dolzarbdir). Masalan, 3c5ee2 kodi uchun bizda quyidagi ma'lumotlar mavjud:
Tahrirlash: in ICAO kodining tavsifi batafsilroq berilgan, men qiziquvchilarga uni o'qishni tavsiya qilaman.
DATA (56 yoki 112 bit) - biz dekodlaydigan haqiqiy ma'lumotlar. Ma'lumotlarning dastlabki 5 biti maydondir Kodni kiriting, saqlanadigan ma'lumotlarning pastki turini o'z ichiga oladi (DF bilan adashtirmaslik kerak). Ushbu turlarning bir nechtasi mavjud:
()
Keling, paketlarning bir nechta misollarini ko'rib chiqaylik.
Samolyotni identifikatsiya qilish
Ikkilik shakldagi misol:
00100 011 000101 010111 000111 110111 110001 111000
Ma'lumotlar maydonlari:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, har bir belgi C1-C8 qatordagi indekslarga mos keladigan kodlarni o'z ichiga oladi:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Satrni dekodlash orqali samolyot kodini olish oson: 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('#', ''))
Havo holati
Agar nom oddiy bo'lsa, unda koordinatalar murakkabroq. Ular 2, juft va toq kadrlar shaklida uzatiladi. Maydon kodi TC = 01011b = 11.
Juft va toq paketlarga misol:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Koordinatalarni hisoblash juda murakkab formula bo'yicha amalga oshiriladi:
()
Men GIS mutaxassisi emasman, shuning uchun u qayerdan kelganini bilmayman. Kim biladi, izohlarda yozing.
Balandlik oddiyroq deb hisoblanadi - ma'lum bir bitga qarab, u 25 yoki 100 futning ko'paytmasi sifatida ifodalanishi mumkin.
Havo tezligi
TC = 19 bilan paket. Bu erda qiziq narsa shundaki, tezlik aniq bo'lishi mumkin, erga nisbatan (Ground Speed) yoki havoda, samolyot sensori (Airspeed) tomonidan o'lchanadi. Ko'p turli sohalar ham uzatiladi:
()
xulosa
Ko'rib turganingizdek, standart nafaqat professionallar, balki oddiy foydalanuvchilar uchun ham foydali bo'lgan ADS-B texnologiyasi qiziqarli simbiozga aylandi. Lekin, albatta, bunda asosiy rolni raqamli SDR qabul qiluvchilarining arzonroq texnologiyasi o'ynadi, bu qurilmaga tom ma'noda gigagertsdan yuqori chastotali signallarni "tinga" olish imkonini beradi.
Standartning o'zida, albatta, yana ko'p narsalar mavjud. Hohlovchilar PDF-ni sahifada ko'rishlari mumkin yoki yuqorida aytib o'tilganiga tashrif buyuring .
Yuqoridagilarning barchasi ko'pchilik uchun foydali bo'lishi dargumon, lekin hech bo'lmaganda uning qanday ishlashi haqida umumiy fikr, umid qilamanki, saqlanib qoladi.
Aytgancha, Python-da tayyor dekoder allaqachon mavjud, siz uni o'rganishingiz mumkin . Va SDR qabul qiluvchilar egalari tayyor ADS-B dekoderini yig'ish va ishga tushirishlari mumkin , bu haqida batafsilroq muhokama qilingan .
Maqolada tasvirlangan parserning manba kodi kesim ostida berilgan. Bu ishlab chiqarishga o'xshamaydigan sinov namunasi, lekin unda ba'zi narsalar ishlaydi va undan yuqorida yozilgan faylni tahlil qilish uchun foydalanish mumkin.
Manba kodi (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()
Umid qilamanki, kimdir qiziqdi, e'tiboringiz uchun rahmat.
Manba: www.habr.com