Bonjou Habr. Pwobableman tout moun ki te janm rankontre oswa wè fanmi oswa zanmi nan yon avyon te itilize sèvis gratis Flightradar24 la. Sa a se yon fason trè pratik yo swiv pozisyon avyon an nan tan reyèl.
В Prensip fonksyònman yon sèvis online konsa te dekri. Nou pral kounye a ale pi devan epi konnen ki done yo te voye ak resevwa soti nan avyon an nan estasyon an k ap resevwa epi dekode li tèt nou lè l sèvi avèk Python.
Istwa
Li evidan, done avyon yo pa transmèt pou itilizatè yo wè sou smartphone yo. Sistèm nan rele ADS-B (Otomatik depandan siveyans-difizyon), epi yo itilize otomatikman transmèt enfòmasyon sou avyon an nan sant kontwòl la - idantifyan li yo, kowòdone, direksyon, vitès, altitid ak lòt done yo transmèt. Précédemment, anvan avènement de sistèm sa yo, dispatcher a te kapab sèlman wè yon pwen sou rada a. Sa pa t ase ankò lè te gen twòp avyon.
Teknikman, ADS-B konsiste de yon transmetè sou yon avyon ki detanzantan voye pake enfòmasyon nan yon frekans jistis segondè nan 1090 MHz (gen lòt mòd, men nou pa tèlman enterese nan yo, depi kowòdone yo transmèt sèlman isit la). Natirèlman, nan adisyon a transmetè a, gen tou yon reseptè yon kote nan ayewopò an, men pou nou, kòm itilizatè, reseptè pwòp nou an enteresan.
By wout la, pou konparezon, premye sistèm sa yo, Airnav Radarbox, ki fèt pou itilizatè òdinè, te parèt an 2007, epi li te koute apeprè $ 900; yon abònman nan sèvis rezo koute yon lòt $ 250 yon ane.
Revizyon sou premye pwopriyetè Ris sa yo ka li sou fowòm nan . Kounye a ke reseptè RTL-SDR yo te vin lajman disponib, yon aparèy menm jan an ka rasanble pou $ 30; plis enfòmasyon sou sa a te nan . Ann ale nan pwotokòl la li menm - ann wè ki jan li fonksyone.
Resevwa siyal yo
Premyèman, siyal la bezwen anrejistre. Siyal la tout antye gen yon dire sèlman 120 mikrosgond, kidonk demonte alèz eleman li yo, yon reseptè SDR ak yon frekans echantiyon nan omwen 5 MHz se dezirab.
Apre anrejistreman, nou resevwa yon dosye WAV ak yon pousantaj echantiyon 5000000 echantiyon/sec; 30 segonn nan anrejistreman sa yo "peze" sou 500MB. Koute li ak yon jwè medya, nan kou, se initil - dosye a pa gen son, men yon siyal radyo dirèkteman nimerik - sa a se egzakteman ki jan Software Defined Radio travay.
Nou pral louvri epi trete dosye a lè l sèvi avèk Python. Moun ki vle fè eksperyans poukont yo ka telechaje yon egzanp anrejistreman .
Ann telechaje fichye a epi gade sa ki anndan an.
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()
Rezilta: nou wè evidan "puls" kont bri nan background.
Chak "batman kè" se yon siyal, estrikti nan ki vizib klèman si ou ogmante rezolisyon an sou graf la.
Kòm ou ka wè, foto a se byen konsistan ak sa yo bay nan deskripsyon an pi wo a. Ou ka kòmanse trete done yo.
Dekodaj
Premyèman, ou bezwen jwenn yon ti kouran. Siyal la tèt li kode lè l sèvi avèk kodaj Manchester:
Soti nan diferans nan nivo nan nibbles li fasil jwenn reyèl "0" ak "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"
Estrikti siyal la li menm se jan sa a:
Ann gade nan jaden yo an plis detay.
DF (Downlink Fòma, 5 Bits) - detèmine kalite mesaj la. Gen plizyè kalite:
()
Nou enterese sèlman nan kalite DF17, paske... Li se sa a ki gen kowòdone yo nan avyon an.
ICAO (24 bits) - kòd entènasyonal inik nan avyon an. Ou ka tcheke avyon an pa kòd li yo (malerezman, otè a te sispann mete ajou baz done a, men li toujou enpòtan). Pou egzanp, pou kòd 3c5ee2 nou gen enfòmasyon sa yo:
Edit: nan Yo bay deskripsyon kòd ICAO a an plis detay; mwen rekòmande pou moun ki enterese li li.
DONE (56 oswa 112 bit) - done aktyèl ke nou pral dekode. Premye 5 bits done yo se jaden an Kalite Kòd, ki gen soutip done yo estoke (pa dwe konfonn ak DF). Gen kèk nan kalite sa yo:
()
Ann gade nan kèk egzanp pakè.
Idantifikasyon avyon
Egzanp nan fòm binè:
00100 011 000101 010111 000111 110111 110001
Jaden done:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, chak karaktè C1-C8 gen kòd ki koresponn ak endis nan liy lan:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Lè yo dekode fisèl la, li fasil pou jwenn kòd avyon an: 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('#', ''))
Pozisyon Airborne
Si non an senp, kowòdone yo pi konplike. Yo transmèt nan fòm lan nan 2, menm ak enpè ankadreman. Kòd jaden TC = 01011b = 11.
Egzanp pake par ak enpè:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Kalkil la nan kowòdone tèt li fèt dapre yon fòmil olye difisil:
()
Mwen pa yon ekspè GIS, kidonk mwen pa konnen ki kote li soti. Ki moun ki konnen, ekri nan kòmantè yo.
Wotè konsidere kòm pi senp - depann sou ti jan an espesifik, li ka reprezante kòm swa yon miltip nan 25 oswa 100 pye.
Vitès Airborne
Pake ak TC=19. Bagay la enteresan isit la se ke vitès la ka swa egzat, parapò ak tè a (Ground Speed), oswa lè, mezire pa yon Capteur avyon (Airspeed). Anpil jaden diferan yo transmèt tou:
()
Konklizyon
Kòm ou ka wè, teknoloji ADS-B te vin tounen yon senbyotik enteresan, lè yon estanda itil pa sèlman pwofesyonèl, men tou, itilizatè òdinè. Men, nan kou, yon wòl kle nan sa a te jwe pa teknoloji a pi bon mache nan reseptè SDR dijital, ki pèmèt aparèy la literalman resevwa siyal ak frekans ki pi wo a yon gigahertz "pou peni."
Nan estanda nan tèt li, nan kou, gen anpil plis. Moun ki enterese ka gade PDF la sou paj la oswa vizite youn ki deja mansyone pi wo a .
Li fasil pou tout sa ki anwo yo pral itil anpil moun, men omwen lide jeneral sou fason li fonksyone, mwen espere, rete.
By wout la, yon dekodeur pare-fè nan Python deja egziste, ou ka etidye li . Ak pwopriyetè reseptè SDR yo ka rasanble ak lanse yon dekodeur ADS-B pare , sa a te diskite an plis detay nan .
Kòd sous analizeur ki dekri nan atik la bay anba a koupe a. Sa a se yon egzanp tès ki pa pretann yo dwe pwodiksyon, men gen kèk bagay ki travay nan li, epi li ka itilize analize dosye a anrejistre pi wo a.
Kòd sous (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()
Mwen espere ke yon moun te enterese, mèsi pou atansyon ou.
Sous: www.habr.com