Hello Habr. Wataƙila duk wanda ya taɓa saduwa ko ganin dangi ko abokai a cikin jirgin sama ya yi amfani da sabis ɗin Flightradar24 kyauta. Wannan hanya ce mai dacewa don bin diddigin matsayin jirgin a ainihin lokacin.
В An bayyana ka'idar aiki na irin wannan sabis na kan layi. Yanzu za mu ci gaba da gano bayanan da ake aikawa da karɓa daga jirgin zuwa tashar karɓa da kuma yanke shi da kanmu ta hanyar amfani da Python.
История
Babu shakka, ba a watsa bayanan jirgin sama don masu amfani su gani a wayoyinsu na zamani. Ana kiran tsarin ADS-B (Automatic dogaro sa ido — watsa shirye-shirye), kuma ana amfani dashi don watsa bayanai ta atomatik game da jirgin zuwa cibiyar sarrafawa - mai ganowa, daidaitawa, jagora, saurin, tsayi, da sauran bayanai ana watsa su. A baya can, kafin zuwan irin waɗannan tsarin, mai aikawa zai iya ganin aya kawai akan radar. Wannan bai isa ba lokacin da jirage suka yi yawa.
A fasaha, ADS-B ya ƙunshi mai watsawa a kan jirgin sama wanda ke aika fakitin bayanai lokaci-lokaci a daidaitaccen mitar 1090 MHz (akwai wasu hanyoyin, amma ba mu da sha'awar su sosai, tunda ana watsa haɗin kai kawai a nan). Tabbas, ban da mai watsawa, akwai kuma mai karɓa a wani wuri a filin jirgin sama, amma a gare mu, a matsayin masu amfani, mai karɓar namu yana da ban sha'awa.
Af, don kwatanta, na farko irin wannan tsarin, Airnav Radarbox, wanda aka ƙera don talakawa masu amfani, ya bayyana a cikin 2007, kuma ya kashe kimanin $900; biyan kuɗin da sabis na cibiyar sadarwa ya biya wani $250 a shekara.
Reviews na wadanda na farko Rasha masu za a iya karanta a kan forum . Yanzu da masu karɓar RTL-SDR suka zama ko'ina, ana iya haɗa irin wannan na'urar akan $ 30; ƙarin game da wannan yana cikin . Bari mu matsa zuwa ga yarjejeniya kanta - bari mu ga yadda take aiki.
Karɓar sigina
Da farko, ana buƙatar rikodin siginar. Gabaɗayan siginar yana da tsawon daƙiƙa 120 kawai, don haka don ƙwace kayan aikinta cikin nutsuwa, mai karɓar SDR tare da mitar samfur na aƙalla 5 MHz yana da kyawawa.
Bayan yin rikodi, muna karɓar fayil ɗin WAV tare da ƙimar ƙima na 5000000 samfurori / sec; 30 seconds na irin wannan rikodin "nauyi" game da 500MB. Sauraron sa tare da na'urar watsa labarai, ba shakka, ba shi da amfani - fayil ɗin ba ya ƙunshi sauti, amma siginar rediyo kai tsaye da aka ƙirƙira - wannan shine yadda Software Defined Rediyo ke aiki.
Za mu buɗe da sarrafa fayil ɗin ta amfani da Python. Wadanda suke son yin gwaji da kansu suna iya zazzage rikodin misali .
Bari mu zazzage fayil ɗin mu ga abin da ke ciki.
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()
Sakamako: muna ganin “bugu-bugu” a fili a kan hayaniyar baya.
Kowane "pulse" sigina ne, tsarin da yake bayyane a fili idan kun ƙara ƙuduri akan jadawali.
Kamar yadda kuke gani, hoton ya yi daidai da abin da aka bayar a bayanin da ke sama. Kuna iya fara sarrafa bayanan.
Ƙaddamarwa
Da farko, kuna buƙatar samun rafi kaɗan. Siginar da kanta an rufa masa asiri ta amfani da rufaffiyar Manchester:
Daga bambancin matakin a cikin nibbles yana da sauƙi don samun ainihin "0" da "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"
Tsarin siginar kanta shine kamar haka:
Bari mu dubi filayen daki-daki.
DF (Format Downlink, 5 bits) - yana ƙayyade nau'in saƙon. Akwai iri da yawa:
()
Muna sha'awar nau'in DF17 kawai, saboda ... Wannan shi ne ya ƙunshi haɗin gwiwar jirgin.
ICAO (24 ragowa) - lambar musamman ta duniya na jirgin sama. Kuna iya duba jirgin ta lambar sa (Abin takaici, marubucin ya daina sabunta bayanan bayanan, amma har yanzu yana da dacewa). Misali, don lambar 3c5ee2 muna da bayanan masu zuwa:
Tace: in An ba da bayanin lambar ICAO daki-daki; Ina ba da shawarar masu sha'awar su karanta shi.
DATA (56 ko 112 bits) - ainihin bayanan da za mu yanke. Na farko 5 na bayanai sune filin Nau'in lamba, yana ɗauke da nau'in nau'in bayanan da ake adanawa (kada a ruɗe da DF). Akwai kadan daga cikin ire-iren wadannan:
()
Bari mu kalli wasu misalan fakiti.
Ganewar jirgin sama
Misali a sigar binary:
00100 011 000101 010111 000111 110111 110001 111000
Filin bayanai:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, kowane hali C1-C8 ya ƙunshi lambobin da suka dace da fihirisa a cikin layi:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ###################################################################################################################################
Ta hanyar yanke kirtani, yana da sauƙi don samun lambar jirgin sama: 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('#', ''))
Matsayin iska
Idan sunan yana da sauƙi, to, haɗin gwiwar sun fi rikitarwa. Ana watsa su a cikin nau'i na 2, ko da kuma firam masu ban sha'awa. Lambar filin TC = 01011b = 11.
Misali na fakiti masu ma'ana da ban mamaki:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Lissafin haɗin kai da kansa yana faruwa ne bisa ga wata dabara mai banƙyama:
()
Ni ba masanin GIS ba ne, don haka ban san inda ya fito ba. Wanene ya sani, rubuta a cikin sharhi.
Ana la'akari da tsayi mafi sauƙi - dangane da ƙayyadaddun ƙayyadaddun, ana iya wakilta shi azaman ko dai nau'in ƙafa 25 ko 100.
Gudun Jirgin Sama
Kunshin tare da TC=19. Abu mai ban sha'awa a nan shi ne, gudun zai iya zama daidai ko dai, dangane da kasa (Ground Speed), ko kuma ta iska, wanda aka auna shi ta hanyar firikwensin jirgi (Airspeed). Ana kuma watsa fage daban-daban da yawa:
()
ƙarshe
Kamar yadda kake gani, fasahar ADS-B ta zama alama mai ban sha'awa, lokacin da ma'auni yana da amfani ba kawai ga masu sana'a ba, har ma ga masu amfani na yau da kullum. Amma ba shakka, muhimmiyar rawa a cikin wannan fasaha ce mai rahusa na masu karɓar SDR na dijital, wanda ke ba da damar na'urar don karɓar sigina a zahiri tare da mitoci sama da gigahertz "don pennies."
A cikin ma'auni kanta, ba shakka, akwai ƙari mai yawa. Masu sha'awar suna iya duba PDF akan shafin ko ziyarci wanda aka riga aka ambata a sama .
Yana da wuya cewa duk abin da ke sama zai zama da amfani ga mutane da yawa, amma aƙalla ra'ayin gaba ɗaya na yadda yake aiki, ina fata, ya kasance.
Af, shirye-shiryen decoder a Python ya riga ya wanzu, zaku iya yin nazarinsa . Kuma masu karɓar SDR na iya haɗawa da ƙaddamar da na'urar ADS-B da aka ƙera , an tattauna wannan dalla-dalla a ciki .
An ba da lambar tushe na parser da aka kwatanta a cikin labarin a ƙasa da yanke. Wannan misali ne na gwaji wanda ba ya yin kamar ana samarwa, amma wasu abubuwa suna aiki a ciki, kuma ana iya amfani da shi don tantance fayil ɗin da aka rubuta a sama.
Lambar tushe (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()
Ina fata wani yana sha'awar, na gode da kulawar ku.
source: www.habr.com