Hello Habr. Mohlomong motho e mong le e mong ea kileng a kopana kapa a bona beng ka eena kapa metsoalle ka sefofaneng o sebelisitse tšebeletso ea mahala ea Flightradar24. Ona ke mokhoa o bonolo haholo oa ho lekola boemo ba sefofane ka nako ea nnete.
В Ho ile ha hlalosoa molao-motheo oa ts'ebetso ea ts'ebeletso e joalo ea inthanete. Hona joale re tla tsoela pele ho fumana hore na ke lintlha life tse rometsoeng le tse amoheloang ho tloha sefofaneng ho ea seteisheneng se amohelang 'me re itlhalose re sebelisa Python.
История
Ho hlakile hore data ea sefofane ha e fetisoe hore basebelisi ba e bone ho li-smartphone tsa bona. Sistimi ena e bitsoa ADS-B (Automatic dependent surveillance-broadcast), 'me e sebelisoa ho fetisetsa tlhahisoleseling ka sefofane setsing sa taolo - sekhetho sa sona, likhokahano, tataiso, lebelo, bophahamo le lintlha tse ling lia fetisoa. Nakong e fetileng, pele ho qaleha ha mekhoa e joalo, motho ea tsamaisang litaba o ne a ka bona ntlha feela ho radar. Sena se ne se se se sa lekana ha ho ne ho e-na le lifofane tse ngata haholo.
Ha e le hantle, ADS-B e na le mochine o fetisang sefofaneng se romellang lipakete tsa tlhahisoleseding nako le nako ka lebelo le phahameng la 1090 MHz (ho na le mekhoa e meng, empa ha re e thahaselle hakaalo, kaha li-coordinate li fetisoa mona feela). Ehlile, ho phaella ho transmitter, ho boetse ho na le moamoheli kae-kae boema-fofane, empa ho rona, joalo ka basebelisi, moamoheli oa rona oa thahasellisa.
Ha e le hantle, ha ho bapisoa, tsamaiso ea pele e joalo, Airnav Radarbox, e etselitsoeng basebelisi ba tloaelehileng, e hlahile ka 2007, 'me e bitsa chelete e ka bang $ 900; ho ngolisa litšebeletso tsa marang-rang ho bitsa $ 250 ka selemo.
Litlhahlobo tsa beng ba pele ba Russia li ka baloa sethaleng . Kaha joale ba amohelang RTL-SDR ba se ba fumaneha hohle, sesebelisoa se ts'oanang se ka bokelloa ka $30; tse ling ka sena li ne li le teng. . Ha re feteleng ho protocol ka boeona - ha re bone hore na e sebetsa joang.
Ho amohela matshwao
Taba ea pele, lets'oao le hloka ho rekota. Letšoao lohle le na le nako ea li-microseconds tse 120 feela, kahoo ho arola likarolo tsa eona hantle, moamoheli oa SDR o nang le sampole ea bonyane 5 MHz oa lakatseha.
Ka mor'a ho rekota, re fumana faele ea WAV e nang le tekanyo ea sampole ea 5000000 sampuli / sec; Metsotsoana ea 30 ea ho rekota joalo e "bekha" e ka bang 500MB. Ho e mamela ka sebapali sa mecha ea phatlalatso, ha e le hantle, ha ho na thuso - faele ha e na molumo, empa e le letšoao la seea-le-moea ka ho toba - sena ke hantle kamoo Software Defined Radio e sebetsang kateng.
Re tla bula le ho sebetsana le faele re sebelisa Python. Ba batlang ho etsa liteko ba le bang ba ka khoasolla mohlala oa ho rekota .
Ha re khoasolle faele mme re bone se ka hare.
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()
Sephetho: re bona "li-pulse" tse hlakileng khahlano le lerata le ka morao.
"Pulse" e 'ngoe le e' ngoe ke pontšo, eo sebopeho sa eona se bonahalang ka ho hlaka haeba u eketsa qeto ho graph.
Joalokaha u ka bona, setšoantšo se lumellana hantle le se fanoeng tlhalosong e ka holimo. U ka qala ho sebetsana le data.
Decoding
Pele, u lokela ho fumana molapo o monyenyane. Letšoao ka boeona le kentsoe ho sebelisoa encoding ea Manchester:
Ho tloha phapang ea boemo ba li-nibbles ho bonolo ho fumana "0" le "1" ea sebele.
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"
Sebopeho sa pontšo ka boeona ke se latelang:
Ha re shebeng masimo ka botlalo.
DF (Downlink Format, 5 bits) - e khetholla mofuta oa molaetsa. Ho na le mefuta e 'maloa:
()
Re thahasella feela mofuta oa DF17, hobane... Ke eona ena e nang le likhokahano tsa sefofane.
ICAO (24 bits) - khoutu e ikhethang ea machaba ea sefofane. U ka hlahloba sefofane ka khoutu ea sona (ka bomalimabe, mongoli o emisitse ho ntlafatsa database, empa e ntse e le ea bohlokoa). Mohlala, bakeng sa khoutu 3c5ee2 re na le lintlha tse latelang:
Fetola: ho Tlhaloso ea khoutu ea ICAO e fanoe ka botlalo; Ke khothaletsa ba thahasellang hore ba e bale.
Data (56 kapa 112 bits) - data ea sebele eo re tla e khetholla. Li-bits tse 5 tsa pele tsa data ke tšimo Mofuta oa Khoutu, e nang le subtype ea data e bolokiloeng (e se ke ea ferekanngoa le DF). Ho na le mefuta e mengata ea mefuta ena:
()
Ha re shebeng mehlala e seng mekae ea liphutheloana.
Boitsebiso ba sefofane
Mohlala ka sebopeho sa binary:
00100 011 000101 010111 000111 110111 110001 111000
Libaka tsa data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, tlhaku ka 'ngoe C1-C8 e na le likhoutu tse tsamaellanang le li-indices moleng:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##################0123456789#####_################XNUMX
Ka ho hlophisa khoele, ho bonolo ho fumana khoutu ea sefofane: 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('#', ''))
Boemo ba moea
Haeba lebitso le le bonolo, joale li-coordinate li rarahane haholoanyane. Li fetisoa ka mokhoa oa 2, esita le liforeimi tse sa tloaelehang. Khoutu ea tšimo TC = 01011b = 11.
Mohlala oa lipakete tse tšoanang le tse sa tloaelehang:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Palo ea li-coordinate ka boeona e etsahala ho latela mokhoa o qhekellang:
()
Ha ke setsebi sa GIS, kahoo ha ke tsebe hore na e tsoa hokae. Ke mang ea tsebang, ngola litlhalosong.
Bophahamo bo nkuoa bo le bobebe - ho ipapisitsoe le sets'oants'o se ikhethileng, bo ka emeloa joalo ka palo ea maoto a 25 kapa 100.
Lebelo la Sefofane
Sephutheloana se nang le TC=19. Ntho e thahasellisang mona ke hore lebelo le ka ba le nepahetseng, le amanang le fatše (Ground Speed), kapa le fofang, le lekantsoeng ke sensor ea sefofane (Airspeed). Likarolo tse ngata tse fapaneng li boetse li fetisoa:
()
fihlela qeto e
Joalokaha u ka bona, theknoloji ea ADS-B e fetohile symbiosis e thahasellisang, ha tekanyetso e sa sebetse feela ho litsebi, empa le ho basebelisi ba tloaelehileng. Empa ehlile, karolo ea bohlokoa ho sena e bapaloa ke thekenoloji e theko e tlaase ea ba amohelang SDR ea dijithale, e lumellang sesebelisoa ho amohela matšoao ka maqhubu a holimo ho gigahertz "bakeng sa lipeni."
Ha e le hantle, litekanyetsong ka boeona, ho na le tse ling tse ngata. Ba thahasellang ba ka sheba PDF leqepheng lena kapa etela e seng e boletsoe ka holimo .
Ha ho na monyetla oa hore tsohle tse kaholimo li tla ba molemo ho ba bangata, empa bonyane mohopolo o akaretsang oa hore na o sebetsa joang, ke ts'epa, o ntse o le teng.
Ka tsela, decoder e lokiselitsoeng ho Python e se e ntse e le teng, u ka e ithuta . 'Me beng ba li-receiver tsa SDR ba ka bokana 'me ba qala mochine o hlophisitsoeng oa ADS-B , sena se ile sa tšohloa ka ho qaqileng haholoanyane ka .
Khoutu ea mohloli oa parser e hlalositsoeng sehloohong e fanoe ka tlase ho sehiloeng. Ona ke mohlala oa teko o sa iketse eka ke oa tlhahiso, empa lintho tse ling li sebetsa ho oona, 'me o ka sebelisoa ho hlalosa faele e tlalehiloeng ka holimo.
Mohloli oa khoutu (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()
Ke tšepa hore ho na le motho ea neng a thahasella, ke leboha tlhokomelo ea hau.
Source: www.habr.com