Салом Хабр. Эҳтимол ҳар касе, ки дар ҳавопаймо хешовандон ё дӯстонашро вохӯрдааст ё гусел кардааст, аз хидмати ройгони Flightradar24 истифода кардааст. Ин роҳи хеле қулайи пайгирии мавқеъи ҳавопаймо дар вақти воқеӣ аст.
В Принсипи кори чунин хидмати онлайн тавсиф шудааст. Мо ҳоло пеш меравем ва мефаҳмем, ки кадом маълумот аз ҳавопаймо ба истгоҳи қабул фиристода мешавад ва қабул карда мешавад ва худамон бо истифода аз Python рамзкушоӣ мекунем.
ҳикояи
Аён аст, ки маълумот дар бораи ҳавопаймо интиқол дода намешавад, то корбарон дар смартфони худ бубинанд. Система ADS-B (Automatic dependent deveillance — broadcast) номида мешавад ва барои ба таври худкор интиқол додани маълумот дар бораи ҳавопаймо ба маркази идоракунӣ истифода мешавад - идентификатор, координатҳо, самт, суръат, баландӣ ва дигар маълумот интиқол дода мешавад. Пештар, то пайдо шудани ин гуна системахо диспетчер танхо нуктаеро дар радар дида метавонист. Вақте ки ҳавопаймоҳо аз ҳад зиёд буданд, ин дигар кофӣ набуд.
Аз ҷиҳати техникӣ, ADS-B аз интиқолдиҳанда дар ҳавопаймо иборат аст, ки давра ба давра бастаҳои иттилоотро бо басомади хеле баланди 1090 МГс мефиристад (режимҳои дигар вуҷуд доранд, аммо мо ба онҳо он қадар таваҷҷӯҳ зоҳир намекунем, зеро координатҳо танҳо дар ин ҷо интиқол дода мешаванд). Албатта, ба гайр аз передатчик дар ягон чо дар фурудгох приёмник хам мавчуд аст, аммо барои мо, истифодабарандагон, приёмники худи мо чолиб аст.
Дар омади гап, барои муқоиса, аввалин чунин система, Airnav Radarbox, ки барои корбарони оддӣ пешбинӣ шудааст, соли 2007 пайдо шуд ва тақрибан 900 доллар арзиш дошт; обуна ба хидматҳои шабакавӣ дар як сол боз 250 доллар арзиш дорад.
Баррасиҳои ин соҳибони аввалини русро дар форум хондан мумкин аст . Ҳоло, ки қабулкунакҳои RTL-SDR ба таври васеъ дастрас шуданд, дастгоҳи шабеҳро метавон бо нархи 30 доллар ҷамъ кард; бештар дар бораи ин дар . Биёед ба худи протокол гузарем - биёед бубинем, ки он чӣ гуна кор мекунад.
Қабули сигналҳо
Аввалан, сигнал бояд сабт карда шавад. Давомнокии тамоми сигнал ҳамагӣ 120 микросония дорад, аз ин рӯ барои ба осонӣ ҷудо кардани ҷузъҳои он, қабулкунандаи SDR бо басомади интихоб на камтар аз 5 МГс матлуб аст.
Пас аз сабт, мо файли WAV-ро бо суръати интихоби 5000000 намуна / сония мегирем; 30 сонияи чунин сабт тақрибан 500 МБ "вазнаш мекунад". Гӯш кардани он бо плеери медиа, албатта, бефоида аст - файл дорои садо нест, балки сигнали радиои мустақими рақамӣ - маҳз ҳамин тавр кор мекунад Software Defined Radio.
Мо файлро бо истифода аз Python мекушоем ва коркард мекунем. Онҳое, ки мехоҳанд мустақилона озмоиш кунанд, метавонанд сабти мисолро зеркашӣ кунанд .
Биёед файлро зеркашӣ кунем ва бубинем, ки дар дохили он чӣ ҳаст.
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()
Натиҷа: мо дар муқобили садои пасзамина "импулсҳо"-и аёнро мебинем.
Ҳар як «набз» сигналест, ки сохтори он дар сурати зиёд кардани қарори графикӣ равшан намоён аст.
Тавре ки шумо мебинед, тасвир бо он чизе, ки дар тавсифи боло оварда шудааст, комилан мувофиқ аст. Шумо метавонед коркарди маълумотро оғоз кунед.
Рамзи рамзгузорӣ
Аввалан, шумо бояд каме ҷараён гиред. Худи сигнал бо истифода аз рамзгузории Манчестер рамзгузорӣ шудааст:
Аз фарқияти сатҳ дар нибблҳо ба осонӣ ба даст овардани "0" ва "1" воқеӣ аст.
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 (Формати поён, 5 бит) - намуди паёмро муайян мекунад. Якчанд намудҳо вуҷуд доранд:
()
Мо танҳо ба навъи DF17 таваҷҷӯҳ дорем, зеро... Ин аст, ки координатҳои ҳавопайморо дар бар мегирад.
ICAO (24 бит) - рамзи байналмилалии беназири ҳавопаймо. Шумо метавонед ҳавопайморо аз рӯи рамзи он тафтиш кунед (мутаассифона, муаллиф навсозии базаи маълумотро қатъ кардааст, аммо он ҳанӯз ҳам муҳим аст). Масалан, барои рамзи 3c5ee2 мо маълумоти зерин дорем:
Таҳрир: дар Тавсифи рамзи ICAO муфассалтар оварда шудааст; Ман тавсия медиҳам, ки ба хоҳишмандон онро мутолиа кунанд.
МАЪЛУМОТ (56 ё 112 бит) - маълумоти воқеие, ки мо рамзкушоӣ мекунем. 5 битҳои аввали маълумот майдон мебошанд Рамзро нависед, дорои зернамуди маълумоти захирашаванда (бо DF омехта нашавад). Шумораи зиёди ин намудҳо мавҷуданд:
()
Биёед якчанд мисоли бастаҳоро дида бароем.
Шиносоии ҳавопаймо
Намуна дар шакли дуӣ:
00100 011 000101 010111 000111 110111 110001 111000
Майдонҳои маълумот:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, ҳар як аломати C1-C8 дорои рамзҳои мувофиқи индексҳои сатр мебошад:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_#############0123456789######
Бо рамзкушоӣ кардани сатр, гирифтани рамзи ҳавопаймо осон аст: 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('#', ''))
Мавқеи ҳавоӣ
Агар ном оддӣ бошад, пас координатҳо мураккабтаранд. Онҳо дар шакли чаҳорчӯбаи 2, ҷуфт ва тоқ интиқол дода мешаванд. Рамзи майдон TC = 01011b = 11.
Намунаи бастаҳои ҷуфт ва тоқ:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ҳисобкунии координатҳо аз рӯи формулаи хеле душвор сурат мегирад:
()
Ман коршиноси GIS нестам, бинобар ин намедонам, ки он аз куҷо меояд. Кӣ медонад, дар коментария нависед.
Баландӣ соддатар ҳисобида мешавад - вобаста ба битҳои мушаххас, он метавонад ҳамчун чанд 25 ё 100 фут ифода карда шавад.
Суръати ҳавоӣ
Баста бо TC = 19. Чизи ҷолиб дар ин ҷо он аст, ки суръат метавонад дақиқ бошад, нисбат ба замин (Суръати замин) ё ҳавоӣ, ки тавассути сенсори ҳавопаймо (Airspeed) чен карда мешавад. Бисёр соҳаҳои гуногун низ интиқол дода мешаванд:
()
хулоса
Тавре ки шумо мебинед, технологияи ADS-B ба симбиози ҷолиб табдил ёфтааст, вақте ки стандарт на танҳо барои мутахассисон, балки барои корбарони оддӣ низ муфид аст. Аммо, албатта, дар ин технологияи арзонтари қабулкунакҳои рақамии SDR нақши калидӣ бозид, ки ба дастгоҳ имкон медиҳад, ки сигналҳоро бо басомадҳои аз гигагерц боло “барои пенниҳо” қабул кунад.
Дар худи стандарт, албатта, хеле зиёд аст. Хоҳишмандон метавонанд дар саҳифаи PDF дидан кунанд ё ба яке аз дар боло зикршуда боздид кунед .
Гумон аст, ки ҳамаи гуфтаҳои дар боло овардашуда барои бисёриҳо муфид бошанд, аммо ҳадди аққал идеяи умумӣ дар бораи он, ки чӣ тавр кор мекунад, ман умедворам, боқӣ мемонад.
Дар омади гап, декодери тайёр дар Python аллакай вуҷуд дорад, шумо метавонед онро омӯзед . Ва соҳибони қабулкунакҳои SDR метавонанд декодери омодаи ADS-B-ро ҷамъ ва ба кор андозанд , дар ин бора муфассалтар сухан рафт .
Рамзи сарчашмаи таҳлилгари дар мақола тавсифшуда дар зер буриш дода мешавад. Ин як мисоли санҷишӣ аст, ки худро истеҳсолот вонамуд намекунад, аммо баъзе чизҳо дар он кор мекунанд ва онро барои таҳлили файли дар боло сабтшуда истифода бурдан мумкин аст.
Рамзи манбаъ (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()
Ман умедворам, ки касе манфиатдор буд, ташаккур барои таваҷҷӯҳи шумо.
Манбаъ: will.com