Kia ora Habr. Ko nga tangata katoa kua tutaki, kua kite ranei i o whanaunga, hoa ranei i runga i te waka rererangi kua whakamahi i te ratonga Flightradar24 koreutu. He huarahi tino watea tenei ki te whai i te tuunga o te waka rererangi i te waa tuturu.
В I whakaahuatia te kaupapa whakahaere o taua ratonga ipurangi. Inaianei ka anga whakamua tatou ki te whakaaro he aha nga raraunga e tukuna ana me te whiwhi mai i te waka rererangi ki te teihana tango me te wetewete ma te whakamahi i te Python.
kōrero
Ma te mohio, kaore nga raraunga rererangi e tukuna kia kitea e nga kaiwhakamahi i runga i o raatau waea atamai. Ko te punaha e kiia ana ko ADS-B (Te Tirotiro Tiaki Aunoa—Paoho), ka whakamahia ki te tuku aunoa i nga korero mo te waka rererangi ki te pokapu whakahaere - ka tukuna tona tohu, taunga, ahunga, tere, teitei me etahi atu raraunga. I mua, i mua i te taenga mai o enei punaha, ka taea e te kaikawe te kite i tetahi waahi i runga i te radar. Kua kore tenei e ranea i te mea he maha rawa nga rererangi.
Ma te hangarau, ko te ADS-B he kaiwhakawhiti i runga i te waka rererangi e tuku ana i nga paanui korero i ia wa i te waa teitei o te 1090 MHz (kei kona ano etahi atu momo, engari kaore matou i te tino aro ki a raatau, na te mea ka tukuna nga taunga ki konei anake). Ko te tikanga, i tua atu i te kaikawe, he kaiwhiwhi ano i tetahi waahi i te taunga rererangi, engari mo matou, hei kaiwhakamahi, he mea whakamere to tatou ake kaiwhiwhi.
Ma te ara, mo te whakataurite, ko te punaha tuatahi, Airnav Radarbox, i hangaia mo nga kaiwhakamahi noa, i puta mai i te tau 2007, me te utu mo te $900; ko te ohaurunga ki nga ratonga whatunga he utu $250 ia tau.
Ko nga korero mo nga rangatira tuatahi o Ruhia ka taea te panui i runga i te huinga . Inaianei kua watea whanuitia nga kaiwhiwhi RTL-SDR, ka taea te whakahiato he taputapu rite mo te $30; ko etahi atu mo tenei kei roto . Ka anga atu ki te kawa ake - kia kite tatou ka pehea te mahi.
Te whiwhi tohu
Tuatahi, me tuhi te tohu. Ko te tohu katoa he 120 moroiti noa te roa, no reira kia pai te wetewete i ona waahanga, he pai te kaikawe SDR me te auau tauira o te 5 MHz neke atu.
Whai muri i te rekoata, ka whiwhi matou i tetahi konae WAV me te reiti tauira 5000000 tauira/hekona; 30 hēkona o taua rekoata "te taumaha" tata ki te 500MB. Ko te whakarongo ki a ia me te kaitoro pāpāho, he pono, he horihori - karekau he oro i roto i te konae, engari he tohu reo irirangi whakarorohiko tika - koinei tonu te mahi a Software Defined Radio.
Ka whakatuwherahia, ka tukatukahia te konae ma te whakamahi i te Python. Ko te hunga e hiahia ana ki te whakamatau i a raatau ake ka taea te tango i tetahi tauira rekoata .
Me tiki ake te konae ka kite he aha kei roto.
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()
Hua: ka kite tatou i nga "pupuhi" e pa ana ki te haruru o muri.
Ko ia "pupuhi" he tohu, ka tino kitea te hanganga ki te whakanuia e koe te taumira ki te kauwhata.
Ka taea e koe te kite, he rite tonu te pikitia ki nga mea e homai ana i te whakaahuatanga i runga ake nei. Ka taea e koe te timata ki te tukatuka i nga raraunga.
Wetewaehere
Tuatahi, me whiwhi koe i tetahi awa iti. Ko te tohu ano ka whakawaeheretia ma te whakawaehere Manchester:
Mai i te rereketanga o nga taumata o nga nibbles he ngawari ki te tiki i te "0" me te "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"
Ko te hanganga o te tohu ano e whai ake nei:
Kia ata tirohia nga mara.
DF (Whakaahua Raro, 5 moka) - ka whakatau i te momo karere. He maha nga momo:
()
Kei te hiahia noa matou ki te momo DF17, na te mea ... Koinei te mea kei roto nga taunga o te waka rererangi.
ICAO (24 bits) - waehere ahurei o te ao o te rererangi. Ka taea e koe te tirotiro i te waka rererangi ma tana waehere (aua koa, kua mutu te whakahou a te kaituhi i te putunga raraunga, engari he mea tika tonu). Hei tauira, mo te waehere 3c5ee2 kei a matou nga korero e whai ake nei:
Whakatika: in Ko te whakamaaramatanga o te waehere ICAO he nui ake nga korero; Ka tūtohu ahau kia panuihia e te hunga e hiahia ana.
RARAUNGA (56, 112 paraka ranei) - nga raraunga tuturu ka wetewetehia e matou. Ko nga moka raraunga tuatahi e 5 ko te mara Pato Waehere, kei roto te momoroto o nga raraunga e rongoa ana (kaore e pohehe ki te DF). He iti noa enei momo:
()
Kia titiro tatou ki etahi tauira o nga kete.
Tautuhinga rererangi
Tauira i roto i te ahua rua:
00100 011 000101 010111 000111 110111 110001 111000
Ngā āpure Raraunga:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, kei ia pūāhua C1-C8 ngā waehere e hāngai ana ki ngā taupū i te rārangi:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######
Ma te wetewete i te aho, he ngawari ki te tiki i te waehere rererangi: 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('#', ''))
Te tūnga rererangi
Mena he ngawari te ingoa, he uaua ake nga taunga. Ka tukuna i roto i te ahua o nga papa 2, rite me te rereke. Waehere mara TC = 01011b = 11.
He tauira o nga paatete taurite me te kerekere:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ko te tatauranga o nga taunga ka puta mai i runga i te tauira tino uaua:
()
Ehara ahau i te tohunga GIS, no reira kare au e mohio no hea. Ko wai ka mohio, tuhia ki nga korero.
He maamaa ake te teitei - i runga i te moka motuhake, ka taea te tohu hei taurea 25, 100 putu ranei.
Te tere o te rangi
Mōkī me te TC=19. Ko te mea whakamere i konei ko te tere ka taea te tika, e pa ana ki te whenua (Te tere o te whenua), ma te hau ranei, ka inehia e te puoro rererangi (Airspeed). He maha nga waahi rereke ka tukuna ano:
()
mutunga
Ka taea e koe te kite, ko te hangarau ADS-B kua waiho hei tohu whakamere, i te wa e whai hua ana te paerewa ehara i te mea ngaio anake, engari ki nga kaiwhakamahi noa. Engari ko te tikanga, ko tetahi mahi nui i roto i tenei i takarohia e te hangarau iti o nga kaiwhiwhi SDR mamati, e taea ai e te taputapu te whiwhi tohu me nga iarere i runga ake i te gigahertz "mo nga pene."
I roto i te paerewa ake, o te akoranga, he nui atu. Ka taea e te hunga e hiahia ana te tiro i te PDF kei runga i te wharangi toro atu ranei ki te mea kua whakahuatia ake nei .
Kaore pea ka whai hua nga mea katoa o runga ake nei ki te tini, engari ko te whakaaro whanui mo te mahi, e tumanako ana ahau, ka mau tonu.
Ma te ara, kei te noho kee he kaiwhakamahara kua rite ki Python, ka taea e koe te ako . A ka taea e nga rangatira o nga kaiwhiwhi SDR te whakahiato me te whakarewa i tetahi kaiwhakamahara ADS-B kua rite , i korerohia tenei mo etahi atu korero i roto i .
Ko te waehere puna o te parser e whakaahuatia ana i roto i te tuhinga ka hoatu i raro i te tapahi. He tauira whakamatautau tenei kaore e kii he hanga, engari kei te mahi etahi mea, ka taea te whakamahi ki te tarai i te konae kua tuhia ki runga ake nei.
Waehere puna (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()
Ko te tumanako i aro mai tetahi, whakawhetai mo to aro.
Source: will.com