Salama Habr. Angamba izay rehetra nihaona na nahita havana na namana tao anaty fiaramanidina dia nampiasa ny serivisy Flightradar24 maimaim-poana. Ity dia fomba tena mety hanaraha-maso ny toeran'ny fiaramanidina amin'ny fotoana tena izy.
В Nofaritana ny fitsipiky ny fiasan'ny serivisy an-tserasera toy izany. Handeha isika izao ary hamantatra izay angon-drakitra alefa sy azo avy amin'ny fiaramanidina mankany amin'ny tobim-pandraisana ary mamadika azy io amin'ny alàlan'ny Python.
История
Mazava ho azy fa tsy ampitaina ho hitan'ny mpampiasa amin'ny findainy ny angon-drakitra fiaramanidina. Ny rafitra dia antsoina hoe ADS-B (Automatic dependent surveillance-broadcast), ary ampiasaina handefasana ny vaovao momba ny fiaramanidina ho any amin'ny foibe fanaraha-maso - ny famantarana, ny fandrindrana, ny fitarihana, ny hafainganam-pandeha, ny haavony ary ny angona hafa dia alefa. Teo aloha, talohan'ny nahatongavan'ny rafitra toy izany, ny dispatcher dia tsy nahita afa-tsy teboka iray tamin'ny radara. Tsy ampy intsony izany rehefa be loatra ny fiaramanidina.
Ara-teknika, ADS-B dia ahitana mpanelanelana amin'ny fiaramanidina izay mandefa tsindraindray fonosana ny vaovao amin'ny somary avo matetika ny 1090 MHz (misy fomba hafa, fa tsy dia liana amin'izy ireo, satria ny fandrindrana ihany no nampitaina eto). Mazava ho azy fa ankoatry ny émetteur dia misy ihany koa ny resevera any ho any amin'ny seranam-piaramanidina, fa ho antsika mpampiasa dia mahaliana ny mpandray antsika.
Raha ny marina, raha ampitahaina, ny rafitra voalohany toy izany, Airnav Radarbox, natao ho an'ny mpampiasa tsotra, dia niseho tamin'ny 2007, ary mitentina $ 900 eo ho eo; ny famandrihana amin'ny serivisy tambajotra dia mitentina $ 250 isan-taona.
Ny hevitra momba ireo tompona Rosiana voalohany dia azo vakiana ao amin'ny forum . Amin'izao fotoana izao dia lasa azo ampiasaina be dia be ny mpandray RTL-SDR, fitaovana mitovy amin'izany dia azo amboarina amin'ny $30; bebe kokoa momba izany dia ao amin'ny . Andao hiroso amin'ny protocol mihitsy - hojerentsika ny fomba fiasa.
Mandray famantarana
Voalohany, mila raketina ny famantarana. Ny famantarana manontolo dia manana faharetan'ny 120 microsegondra fotsiny, noho izany dia ilaina ny mpandray SDR miaraka amin'ny fatran'ny santionany farafahakeliny 5 MHz, mba hanesorana ireo singa ao aminy.
Aorian'ny firaketana dia mahazo rakitra WAV izahay miaraka amin'ny santionany 5000000 santionany / seg; 30 segondra amin'ny fandraketana toy izany dia "lanja" eo amin'ny 500MB. Ny fihainoana azy miaraka amin'ny mpilalao haino aman-jery, mazava ho azy, dia tsy misy ilàna azy - tsy misy feo ny rakitra, fa mari-pamantarana onjam-peo mivantana - izany no fomba fiasan'ny Software Defined Radio.
Hanokatra sy hikarakara ilay rakitra amin'ny alàlan'ny Python izahay. Ireo izay te-hanandrana samirery dia afaka misintona rakitra ohatra .
Aleo alaina ny rakitra dia hojerentsika izay ao anatiny.
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()
Vokany: mahita “pulse” miharihary manoloana ny tabataba ambadika.
Ny "pulse" tsirairay dia famantarana, ny firafitra izay hita mazava tsara raha ampitomboinao ny fanapahan-kevitra eo amin'ny grafika.
Araka ny hitanao dia mifanaraka tsara amin'ny voalaza etsy ambony ny sary. Afaka manomboka manodina ny angona ianao.
Decoding
Voalohany dia mila maka rano kely ianao. Ny mari-pamantarana mihitsy dia voakodia amin'ny alàlan'ny famandrihana Manchester:
Avy amin'ny fahasamihafan'ny haavon'ny nibbles dia mora ny mahazo ny tena "0" sy "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"
Ny firafitry ny famantarana dia toy izao manaraka izao:
Andeha hojerentsika amin'ny antsipiriany ny saha.
DF (Format Downlink, 5 bit) - mamaritra ny karazana hafatra. Misy karazany maromaro:
()
Karazana DF17 ihany no mahaliana anay, satria... Io no ahitana ny coordinates ny fiaramanidina.
ICAO (24 bits) - kaody tokana iraisam-pirenena momba ny fiaramanidina. Azonao atao ny manamarina ny fiaramanidina amin'ny alàlan'ny kaody (indrisy fa najanon'ny mpanoratra ny fanavaozana ny angon-drakitra, fa mbola ilaina ihany). Ohatra, ho an'ny code 3c5ee2 dia manana ireto fampahalalana manaraka ireto izahay:
Ahitsio : in Ny famaritana ny kaody ICAO dia omena amin'ny antsipiriany bebe kokoa; Manoro hevitra aho mba hamaky azy ireo izay liana.
NY FANAZAVANA (56 na 112 bits) - ny tena angon-drakitra izay ho decode. Ny data 5 bit voalohany dia ny saha Karazana kaody, misy ny sobika amin'ny angona voatahiry (tsy afangaro amin'ny DF). Misy vitsivitsy amin'ireto karazany ireto:
()
Andeha isika hijery ohatra vitsivitsy momba ny fonosana.
Famantarana fiaramanidina
Ohatra amin'ny endrika binary:
00100 011 000101 010111 000111 110111 110001 111000
Data saha:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, ny tarehintsoratra C1-C8 tsirairay dia misy kaody mifanaraka amin'ny indices ao amin'ny tsipika:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######
Amin'ny alàlan'ny famongorana ny tady dia mora ny mahazo ny kaody fiaramanidina: 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('#', ''))
Toetran'ny rivotra
Raha tsotra ny anarana dia sarotra kokoa ny coordinates. Izy ireo dia ampitaina amin'ny endrika 2, mitovy sy hafahafa. Kaody saha TC = 01011b = 11.
Ohatra amin'ny fonosana mitovy sy hafahafa:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ny fikajiana ny koordinate dia mitranga araka ny fomba fiasa somary sarotra:
()
Tsy manam-pahaizana momba ny GIS aho, ka tsy fantatro hoe avy aiza izy io. Iza no mahalala, soraty ao amin'ny fanehoan-kevitra.
Ny haavony dia heverina ho tsotra kokoa - miankina amin'ny bitika manokana, dia azo aseho amin'ny endrika maromaro 25 na 100 metatra.
Velocity an'habakabaka
Package miaraka amin'ny TC=19. Ny mahaliana eto dia ny hafainganam-pandeha dia mety ho marina, raha oharina amin'ny tany (Ground Speed), na an'habakabaka, refesina amin'ny alàlan'ny sensor fiaramanidina (Airspeed). Ampitaina ihany koa ny sehatra maro samihafa:
()
famaranana
Araka ny hitanao, ny teknolojia ADS-B dia lasa symbiose mahaliana, raha ny fenitra dia mahasoa tsy ho an'ny matihanina ihany, fa koa ho an'ny mpampiasa tsotra. Saingy mazava ho azy, anjara lehibe amin'izany no nilalao ny teknolojia mora kokoa amin'ny mpandray SDR nomerika, izay ahafahan'ilay fitaovana mandray ara-bakiteny ireo famantarana miaraka amin'ny frequences mihoatra ny gigahertz "ho an'ny pennies."
Ao amin'ny fenitra mihitsy, mazava ho azy, misy zavatra maro hafa. Ireo liana dia afaka mijery ny PDF amin'ny pejy na tsidiho ilay efa voalaza etsy ambony .
Tsy azo inoana fa hahasoa ny maro ireo voalaza etsy ambony ireo, fa farafaharatsiny ny hevitra ankapobeny momba ny fomba fiasa, manantena aho fa mijanona.
Raha ny marina, efa misy ny decoder efa vita amin'ny Python, azonao atao ny mianatra azy io . Ary ny tompon'ny mpandray SDR dia afaka mivory sy manangana décoder ADS-B efa vita , noresahina tamin'ny antsipiriany bebe kokoa izany tao amin'ny .
Ny kaody loharanon'ny parser voalaza ao amin'ny lahatsoratra dia omena eo ambanin'ny fanapahana. Ity dia ohatra andrana izay tsy miseho ho famokarana, fa misy zavatra miasa ao, ary azo ampiasaina handinihana ny rakitra voarakitra etsy ambony.
Kaody loharano (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()
Manantena aho fa misy olona liana, misaotra amin'ny fiheveranao.
Source: www.habr.com