Sawubona Habr. Cishe wonke umuntu owake wahlangana noma wabona izihlobo noma abangani endizeni usebenzise isevisi yamahhala ye-Flightradar24. Lena indlela elula kakhulu yokulandelela indawo yendiza ngesikhathi sangempela.
В Umgomo wokusebenza wesevisi enjalo eku-inthanethi wachazwa. Manje sizoqhubeka sithole ukuthi iyiphi idatha ethunyelwayo netholwayo esuka endizeni iye esiteshini sokwamukela futhi sizihlukanisele yona sisebenzisa iPython.
История
Ngokusobala, idatha yendiza ayidluliselwa ukuze abasebenzisi bayibone kuma-smartphone abo. Lolu hlelo lubizwa nge-ADS-B (Automatic dependent surveillance—broadcast), futhi lusetshenziselwa ukudlulisa ngokuzenzakalelayo ulwazi olumayelana nendiza esikhungweni sokulawula - isihlonzi sayo, izixhumanisi, isiqondiso sayo, isivinini, ukuphakama kanye neminye imininingwane iyadluliselwa. Ngaphambilini, ngaphambi kokufika kwezinhlelo ezinjalo, i-dispatcher yayingabona iphuzu kuphela ku-radar. Lokhu kwakungasanele lapho kunezindiza eziningi kakhulu.
Ngobuchwepheshe, i-ADS-B iqukethe i-transmitter endizeni ethumela amaphakethe olwazi ngezikhathi ezithile ngemvamisa ephezulu kakhulu ye-1090 MHz (zikhona ezinye izindlela, kepha asinandaba nazo kangako, ngoba izixhumanisi zithunyelwa lapha kuphela). Vele, ngaphezu komthumeli, kukhona nomamukeli endaweni ethile esikhumulweni sezindiza, kodwa kithina, njengabasebenzisi, umamukeli wethu siqu uyathakazelisa.
Ngendlela, ukuqhathanisa, uhlelo lokuqala olunjalo, i-Airnav Radarbox, eyenzelwe abasebenzisi abajwayelekile, yavela ngo-2007, futhi ibiza cishe u-$ 900; ukubhaliswa kwezinsizakalo zenethiwekhi kubiza enye i-$ 250 ngonyaka.
Ukubuyekezwa kwalabo abanikazi bokuqala baseRussia kungafundwa esithangamini . Manje njengoba izamukeli ze-RTL-SDR sezitholakala kabanzi, idivayisi efanayo ingahlanganiswa ngo-$30; okuningi mayelana nalokhu bekungaphakathi. . Asiqhubekele kuphrothokholi ngokwayo - ake sibone ukuthi isebenza kanjani.
Ukwamukela amasignali
Okokuqala, isignali idinga ukurekhodwa. Isignali yonke inobude besikhathi sama-microseconds angu-120 kuphela, ngakho-ke ukuze uhlukanise kahle izingxenye zayo, isamukeli se-SDR esinemvamisa yesampula okungenani engu-5 MHz siyathandeka.
Ngemva kokurekhoda, sithola ifayela le-WAV elinenani lesampula lamasampuli angu-5000000/isekhondi; imizuzwana engu-30 yokurekhoda okunjalo “isisindo” cishe esingu-500MB. Ukuyilalela ngesidlali semidiya, vele, akusizi - ifayela alinawo umsindo, kodwa isignali yomsakazo ekhonjiswe ngqo - yiyona ndlela I-Software Defined Radio isebenza ngayo.
Sizovula futhi sicubungule ifayela sisebenzisa iPython. Labo abafuna ukuzivivinya ngokwabo bangalanda isibonelo sokurekhoda .
Ake silande ifayela futhi sibone ukuthi yini ngaphakathi.
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()
Umphumela: sibona “ama-pulse” asobala ngomsindo wangemuva.
I-"pulse" ngayinye iyisiginali, isakhiwo esibonakala ngokucacile uma ukhulisa isinqumo kugrafu.
Njengoba ubona, isithombe siyahambisana nalokho okushiwo encazelweni engenhla. Ungaqala ukucubungula idatha.
Ukuqopha
Okokuqala, udinga ukuthola ukusakaza kancane. Isiginali ngokwayo ibhalwe ngekhodi kusetshenziswa umbhalo we-Manchester:
Kusukela kumehluko wezinga kuma-nibbles kulula ukuthola u-“0” no-“1” wangempela.
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"
Isakhiwo sesignali ngokwayo simi kanje:
Ake sibheke izinkambu ngokuningiliziwe.
DF (Ifomethi ye-Downlink, amabhithi angu-5) - inquma uhlobo lomlayezo. Kunezinhlobo eziningana:
()
Sithanda kuphela uhlobo lwe-DF17, ngoba... Yilokhu okuqukethe izixhumanisi zendiza.
I-ICAO (24 bits) - ikhodi yamazwe ngamazwe eyingqayizivele yendiza. Ungahlola indiza ngekhodi yayo (ngeshwa, umbhali uyekile ukubuyekeza i-database, kodwa isabalulekile). Isibonelo, kukhodi 3c5ee2 sinolwazi olulandelayo:
Hlela: ku Incazelo yekhodi ye-ICAO inikezwe imininingwane eyengeziwe; ngincoma ukuthi labo abanentshisekelo bayifunde.
IDATHA (amabhithi angu-56 noma angu-112) - idatha yangempela esizoyinquma. Izingcezu zokuqala ezi-5 zedatha ziyinkambu Uhlobo Lekhodi, equkethe uhlobo oluncane lwedatha egcinwayo (akumele kudidaniswe ne-DF). Kukhona ezimbalwa zalezi zinhlobo:
()
Ake sibheke izibonelo ezimbalwa zamaphakheji.
Ukuhlonza indiza
Isibonelo ngefomu kanambambili:
00100 011 000101 010111 000111 110111 110001 111000
Izinkambu 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, uhlamvu ngalunye C1-C8 luqukethe amakhodi ahambisana nezinkomba emgqeni:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_#################0123456789#######
Ngokuqopha intambo, kulula ukuthola ikhodi yendiza: 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 esemoyeni
Uma igama lilula, khona-ke izixhumanisi ziyinkimbinkimbi kakhulu. Asakazwa ngendlela yamafreyimu angu-2, alinganayo futhi angajwayelekile. Ikhodi yenkundla TC = 01011b = 11.
Isibonelo samaphakethe alinganayo nayinqaba:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ukubalwa kwezixhumanisi ngokwazo kwenzeka ngokuya ngefomula ekhohlisayo:
()
Angiyena uchwepheshe we-GIS, ngakho-ke angazi ukuthi ivelaphi. Kwazi bani, bhala kumazwana.
Ubude bubhekwa bulula - kuye ngokuthi ibhithi ethile, bungamelwa njengokuphindaphinda kwamafidi angama-25 noma ayi-100.
I-Airborne Velocity
Iphakheji ene-TC=19. Into ethokozisayo lapha ukuthi ijubane lingaba elinembile, uma liqhathaniswa nomhlabathi (Ijubane Lomhlaba), noma libe semoyeni, likalwa ngenzwa yendiza (Airspeed). Izinkambu eziningi ezahlukene nazo zisakazwa:
()
isiphetho
Njengoba ubona, ubuchwepheshe be-ADS-B buye baba yi-symbiosis ethakazelisayo, lapho izinga liwusizo hhayi kuphela kochwepheshe, kodwa nakubasebenzisi abavamile. Kodwa-ke, indima ebalulekile kulokhu yadlalwa ubuchwepheshe obushibhile bezamukeli zedijithali ze-SDR, okuvumela idivayisi ukuthi yamukele ngokoqobo amasiginali ngamaza angaphezulu kwegigahertz “ngamapeni.”
Endinganisweni ngokwayo, yiqiniso, kukhona okuningi. Labo abanentshisekelo bangabuka i-PDF ekhasini noma vakashela lesi esishiwo ngenhla .
Akunakwenzeka ukuthi konke okungenhla kuzoba usizo kwabaningi, kodwa okungenani umqondo ojwayelekile wokuthi usebenza kanjani, ngethemba, usekhona.
Ngendlela, i-decoder esenziwe ngomumo ePython isivele ikhona, ungayifunda . Futhi abanikazi babamukeli be-SDR bangahlangana futhi baqalise isiqophi se-ADS-B esenziwe ngomumo , lokhu kwaxoxwa ngakho kabanzi ku .
Ikhodi yomthombo yomhlahleli echazwe esihlokweni inikezwe ngezansi kokusikwa. Lesi isibonelo sokuhlola esingazenzisi ukukhiqiza, kodwa ezinye izinto zisebenza kuso, futhi singasetshenziswa ukuhlaziya ifayela elirekhodwe ngenhla.
Ikhodi 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()
Ngethemba ukuthi othile ube nentshisekelo, ngiyabonga ngokunaka kwakho.
Source: www.habr.com