Molo Habr. Mhlawumbi wonke umntu owakhe wadibana okanye wabona izalamane okanye abahlobo kwinqwelomoya uye wasebenzisa inkonzo yasimahla yeFlightradar24. Le yindlela elula kakhulu yokulandelela indawo yenqwelomoya ngexesha lokwenyani.
В Umgaqo wokusebenza wenkonzo enjalo ye-intanethi uchazwe. Ngoku siza kuqhubela phambili kwaye sifumanise ukuba yeyiphi idatha ethunyelwayo kwaye ifunyenwe kwinqwelomoya ukuya kwisikhululo sokufumana kwaye siyichaze ngokwethu sisebenzisa iPython.
История
Ngokucacileyo, idatha yenqwelomoya ayithunyelwa ukuba abasebenzisi bayibone kwii-smartphones zabo. Inkqubo ibizwa ngokuba yi-ADS-B (i-Automatic dependent surveillance-broadcast), kwaye isetyenziselwa ukuthumela ngokuzenzekelayo ulwazi malunga nenqwelomoya kwiziko lolawulo - isihlonzi sayo, ulungelelwaniso, ulwalathiso, isantya, ukuphakama kunye nezinye iinkcukacha zihanjiswa. Ngaphambili, ngaphambi kokufika kweenkqubo ezinjalo, i-dispatcher inokubona kuphela inqaku kwi-radar. Oku kwakunganelanga xa kwakukho iinqwelo-moya ezininzi kakhulu.
Ngokobuchwephesha, i-ADS-B iqulethe i-transmitter kwinqwelomoya ehlala ithumela iipakethi ngolwazi kwi-frequency ephezulu ye-1090 MHz (kukho ezinye iindlela, kodwa asinamdla kuzo, kuba ulungelelwaniso luhanjiswa apha kuphela). Ewe kunjalo, ukongeza kwi-transmitter, kukwakho umamkeli kwindawo ethile kwisikhululo seenqwelomoya, kodwa kuthi, njengabasebenzisi, owamkeli wethu unomdla.
Ngendlela, ukuthelekisa, inkqubo yokuqala enjalo, i-Airnav Radarbox, eyenzelwe abasebenzisi abaqhelekileyo, yavela kwi-2007, kwaye ixabisa malunga ne-$ 900; ukubhaliselwa kwiinkonzo zenethiwekhi kubiza enye i-$ 250 ngonyaka.
Ukuphononongwa kwabo banini bokuqala baseRashiya banokufundwa kwiforum . Ngoku ekubeni abamkeli be-RTL-SDR sele befumaneka ngokubanzi, isixhobo esifanayo sinokudityaniswa nge-30 yeedola; ngakumbi malunga noku bekungaphakathi. . Masiqhubele phambili kwiprotocol ngokwayo - masibone ukuba isebenza njani.
Ukufumana imiqondiso
Okokuqala, isignali kufuneka irekhodwe. Umqondiso wonke unobude be-microseconds ezili-120 kuphela, ngoko ke ukudibanisa ngokukhululekileyo amacandelo ayo, umamkeli we-SDR kunye nesampulu yesampulu ubuncinane ye-5 MHz iyanqweneleka.
Emva kokurekhoda, sifumana ifayile ye-WAV kunye nesantya sesampuli ye-5000000 isampuli / isekhondi; Imizuzwana ye-30 yokurekhoda okunjalo "ubunzima" malunga ne-500MB. Ukuyiphulaphula kunye nomdlali wemidiya, ngokuqinisekileyo, akunamsebenzi - ifayile ayinayo isandi, kodwa isignali yerediyo edijithisiweyo ngokuthe ngqo - yile ndlela kanye i-Software Defined Radio isebenza ngayo.
Siza kuvula kwaye siqhube ifayile sisebenzisa iPython. Abo bafuna ukuzama ngokwabo banokukhuphela umzekelo wokurekhoda .
Masikhuphele ifayile kwaye sibone okungaphakathi.
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()
Isiphumo: sibona “ii-pulses” ezicacileyo ngokubhekiselele kwingxolo yangasemva.
"I-pulse" nganye isignali, isakhiwo esibonakala ngokucacileyo ukuba ukwandisa isisombululo kwigrafu.
Njengoko ubona, umfanekiso uhambelana ngokupheleleyo noko kunikwe kule nkcazo ingentla. Ungaqalisa ukucubungula idatha.
Ukuguqula iikhowudi
Okokuqala, kufuneka ufumane umsinga omncinci. Umqondiso ngokwawo ufakwe ngekhowudi usebenzisa i-encoding yaseManchester:
Ukusuka kumahluko wenqanaba kwii-nibbles kulula ukufumana i-real "0" kunye ne "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"
Ubume bomqondiso ngokwawo buhamba ngolu hlobo lulandelayo:
Makhe sijonge amasimi ngokubanzi.
DF (Downlink Format, 5 bits) - imisela uhlobo lomyalezo. Kukho iintlobo ezininzi:
()
Sinomdla kuphela kuhlobo lwe-DF17, kuba... Yile nto equlethe ulungelelwaniso lwenqwelomoya.
ICAO (iibhithi ezingama-24) - ikhowudi eyodwa yamazwe ngamazwe yenqwelomoya. Unokuyijonga indiza ngekhowudi yayo (Ngelishwa, umbhali uyekile ukuhlaziya isiseko sedatha, kodwa sisasebenza). Umzekelo, kwikhowudi 3c5ee2 sinolwazi olulandelayo:
Hlela: kwi Inkcazo yekhowudi ye-ICAO inikwe ngakumbi; Ndincoma ukuba abo banomdla bayifunde.
ULWAZI (amasuntswana angama-56 okanye angama-112) - eyona datha siza kuyicacisa. I-5 bits yokuqala yedatha yintsimi Uhlobo lweKhowudi, equlathe udidi olusezantsi lwedatha egciniweyo (akufunekanga ibhidaniswe neDF). Kukho uninzi lwezi ntlobo:
()
Makhe sijonge imizekelo embalwa yeepakethe.
Ukuchongwa kwenqwelomoya
Umzekelo kwimo yokubini:
00100 011 000101 010111 000111 110111 110001 111000
Iindawo zedatha:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
I-TC = 00100b = 4, uphawu ngalunye C1-C8 luqulethe iikhowudi ezihambelana ne-indices kumgca:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_#################0123456789#######
Ngokucacisa umtya, kulula ukufumana ikhowudi yenqwelomoya: 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('#', ''))
Indawo yasemoyeni
Ukuba igama lilula, ke ulungelelwaniso lunzima ngakumbi. Zihanjiswa ngendlela ye-2, kunye nezakhelo ezingaqhelekanga. Ikhowudi yendawo TC = 01011b = 11.
Umzekelo weepakethi ezilinganayo kunye nezingaqhelekanga:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ubalo lolungelelwaniso ngokwalo lwenzeka ngokwefomula enenkohliso:
()
Andiyoncutshe yeGIS, ngoko ke andazi ukuba ivela phi. Ngubani owaziyo, bhala kwizimvo.
Ubude buthathwa bulula - ngokuxhomekeke kwisuntswana elithile, bunokumelwa njengophindaphindo lweenyawo ezingama-25 okanye ezili-100.
Isantya sasemoyeni
Ipakethe ene-TC=19. Into enomdla apha kukuba isantya sinokuchaneka, ngokubhekiselele kumhlaba (iSpeed Speed), okanye i-airborne, ilinganiswa yinzwa yendiza (i-Airspeed). Iinkalo ezininzi ezahlukeneyo nazo ziyahanjiswa:
()
isiphelo
Njengoko ubona, iteknoloji ye-ADS-B iye yaba yi-symbiosis enomdla, xa umgangatho uluncedo kungekuphela nje kwiingcali, kodwa nakubasebenzisi abaqhelekileyo. Kodwa ewe, indima ephambili koku idlalwe yitekhnoloji engabizi kakhulu yabamkeli be-SDR yedijithali, evumela isixhobo ukuba sifumane ngokoqobo imiqondiso enemiza engaphezulu kwegigahertz “yepeni.”
Kumgangatho ngokwawo, kunjalo, kuninzi ngakumbi. Abo banomdla banokujonga iPDF kwiphepha okanye undwendwele le sele ikhankanyiwe ngasentla .
Akunakwenzeka ukuba konke oku kungasentla kuya kuba luncedo kwabaninzi, kodwa ubuncinci umbono oqhelekileyo wendlela esebenza ngayo, ndiyathemba, uhlala uhleli.
Ngendlela, i-decoder esele ilungile kwiPython sele ikhona, unokuyifunda . Kwaye abanini abafumana i-SDR banokuhlanganisana kwaye baqalise i-decoder ye-ADS-B esele yenziwe , oku kwaxoxwa ngokubanzi kwi .
Ikhowudi yomthombo we-parser echazwe kwinqaku inikwe ngaphantsi kwe-cut. Lo ngumzekelo wokuvavanya ongenzi ngathi imveliso, kodwa ezinye izinto zisebenza kuyo, kwaye ingasetyenziselwa ukucazulula ifayile erekhodiweyo ngasentla.
Ikhowudi yomthombo (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()
Ndiyathemba ukuba kukho umntu onomdla, enkosi ngengqalelo yakho.
umthombo: www.habr.com