Hello Habr. Marahil lahat ng nakakilala o nakakita ng mga kamag-anak o kaibigan sa isang eroplano ay gumamit ng libreng serbisyo ng Flightradar24. Ito ay isang napaka-maginhawang paraan upang subaybayan ang posisyon ng sasakyang panghimpapawid sa real time.
В inilarawan ang prinsipyo ng pagpapatakbo ng naturang online na serbisyo. Ngayon ay lalakad pa tayo at alamin kung anong data ang ipinadala at natanggap mula sa sasakyang panghimpapawid patungo sa istasyon ng pagtanggap, at i-decode ito mismo gamit ang Python.
Kuwento
Malinaw, hindi ibinabahagi ang data ng sasakyang panghimpapawid para makita ng mga user sa kanilang mga smartphone. Ang system ay tinatawag na ADS-B (Automatic dependent surveillance—broadcast), at ginagamit upang awtomatikong magpadala ng impormasyon tungkol sa sasakyang panghimpapawid sa control center - ang pagkakakilanlan nito, mga coordinate, direksyon, bilis, altitude at iba pang data ay ipinapadala. Noong nakaraan, bago ang pagdating ng naturang mga sistema, ang dispatcher ay nakikita lamang ng isang tuldok sa radar. Hindi ito sapat noong napakaraming eroplano.
Sa teknikal, ang ADS-B ay binubuo ng isang transmiter ng sasakyang panghimpapawid na pana-panahong nagpapadala ng mga packet na may impormasyon sa medyo mataas na dalas ng 1090 MHz (mayroong iba pang mga mode, ngunit hindi sila kawili-wili sa amin, dahil ang mga coordinate ay ipinadala lamang dito). Siyempre, bilang karagdagan sa transmitter, mayroon ding receiver sa isang lugar sa paliparan, ngunit para sa amin, bilang mga gumagamit, ang aming sariling receiver ay interesado.
Sa pamamagitan ng paraan, para sa paghahambing, ang unang naturang sistema, ang Airnav Radarbox, na idinisenyo para sa mga ordinaryong gumagamit, ay lumitaw noong 2007 at nagkakahalaga ng halos $ 900, at ang isang subscription sa mga serbisyo ng network ay nagkakahalaga ng halos $ 250 bawat taon.
Ang mga pagsusuri sa mga unang may-ari ng Russia ay mababasa sa forum . Ngayon na ang mga RTL-SDR receiver ay naging malawakang magagamit, ang isang katulad na aparato ay maaaring tipunin para sa $30, higit pa tungkol dito ay nasa . Magpapatuloy tayo sa protocol mismo - tingnan natin kung paano ito gumagana.
Pagtanggap ng signal
Una, ang signal ay dapat na naitala. Ang buong signal ay may tagal lamang na 120 microseconds, kaya upang kumportableng i-disassemble ang mga bahagi nito, isang SDR receiver na may sampling rate na hindi bababa sa 5 MHz ay kanais-nais.
Pagkatapos ng pag-record, nakakakuha kami ng WAV file na may sampling rate na 5000000 samples / sec, 30 segundo ng naturang recording ay "timbang" ng humigit-kumulang 500MB. Siyempre, walang silbi na pakinggan ito gamit ang isang media player - ang file ay hindi naglalaman ng tunog, ngunit isang direktang digitized na signal ng radyo - ito ay kung paano gumagana ang Software Defined Radio.
Bubuksan at ipoproseso namin ang file gamit ang Python. Ang mga gustong mag-eksperimento nang mag-isa ay maaaring mag-download ng halimbawang pag-record .
I-download natin ang file at tingnan kung ano ang nasa loob.
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()
Resulta: nakikita namin ang malinaw na "mga impulses" laban sa background ng ingay.
Ang bawat "impulse" ay isang senyales, ang istraktura nito ay malinaw na nakikita kung tataasan mo ang resolution sa graph.
Tulad ng nakikita mo, ang larawan ay medyo pare-pareho sa kung ano ang ibinigay sa paglalarawan sa itaas. Maaari mong simulan ang pagproseso ng data.
Pagde-decode
Una, kailangan mong makakuha ng kaunting stream. Ang signal mismo ay naka-encode gamit ang manchester encoding:
Madaling makakuha ng tunay na "0" at "1" mula sa pagkakaiba ng antas sa mga nibbles.
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"
Ang istraktura ng signal mismo ay ang mga sumusunod:
Tingnan natin ang mga patlang nang mas detalyado.
DF (Format ng Downlink, 5 bits) — tinutukoy ang uri ng mensahe. Mayroong ilang mga uri:
()
Kami ay interesado lamang sa uri ng DF17, dahil naglalaman ito ng mga coordinate ng sasakyang panghimpapawid.
ICAO (24 bits) — ang internasyonal na natatanging code ng sasakyang panghimpapawid. Maaari mong suriin ang sasakyang panghimpapawid sa pamamagitan ng code nito (Sa kasamaang palad, ang may-akda ay tumigil sa pag-update ng database, ngunit ito ay may kaugnayan pa rin). Halimbawa, para sa code 3c5ee2 mayroon kaming sumusunod na impormasyon:
Edit: sa ang paglalarawan ng ICAO code ay ibinigay nang mas detalyado, inirerekomenda ko na basahin ito ng mga interesado.
DATA (56 o 112 bits) - ang aktwal na data na aming ide-decode. Ang unang 5 bits ng data ay ang field Uri ng CodeA na naglalaman ng subtype ng nakaimbak na data (hindi malito sa DF). Mayroong ilan sa mga ganitong uri:
()
Tingnan natin ang ilang halimbawa ng mga pakete.
pagkakakilanlan ng sasakyang panghimpapawid
Binary na halimbawa:
00100 011 000101 010111 000111 110111 110001 111000
Mga field ng data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, ang bawat character na C1-C8 ay naglalaman ng mga code na tumutugma sa mga indeks sa string:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######
Kapag na-decode ang string, madaling makuha ang aircraft code: 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('#', ''))
posisyon sa hangin
Kung ang lahat ay simple sa pangalan, pagkatapos ay sa mga coordinate ito ay mas kumplikado. Ang mga ito ay ipinadala bilang 2x, kakaiba at kahit na mga frame. Field code TC = 01011b = 11.
Isang halimbawa ng pantay at kakaibang packet:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ang pagkalkula ng mga coordinate mismo ay nagaganap ayon sa isang medyo nakakalito na formula:
()
Hindi ako isang GIS specialist, kaya hindi ko alam kung saan ito nanggaling. Sino ang nakakaalam, isulat sa mga komento.
Itinuturing na mas simple ang taas - depende sa partikular na bit, maaari itong katawanin bilang isang multiple na 25 o 100 talampakan.
Airborne Velocity
Packet na may TC=19. Ang kagiliw-giliw na bagay dito ay ang bilis ay maaaring parehong tumpak, na may kaugnayan sa lupa (Ground Speed), at hangin, na sinusukat ng sensor ng sasakyang panghimpapawid (Airspeed). Maraming iba't ibang larangan ang naipasa din:
()
Konklusyon
Tulad ng nakikita mo, ang teknolohiya ng ADS-B ay naging isang kawili-wiling simbiyos, kapag ang isang pamantayan ay kapaki-pakinabang hindi lamang para sa mga propesyonal, kundi pati na rin para sa mga ordinaryong gumagamit. Ngunit siyempre, ang pangunahing papel dito ay nilalaro ng mas murang teknolohiya ng mga digital SDR receiver, na nagpapahintulot sa device na literal na "para sa isang sentimos" na makatanggap ng mga signal na may dalas sa itaas ng gigahertz.
Sa pamantayan mismo, siyempre, higit pa. Maaaring tingnan ng mga interesado ang PDF sa page o bisitahin ang nabanggit na sa itaas .
Ito ay malamang na ang lahat ng nasa itaas ay magiging kapaki-pakinabang sa marami, ngunit hindi bababa sa pangkalahatang ideya kung paano ito gumagana, umaasa ako, ay nananatili.
Sa pamamagitan ng paraan, mayroon nang isang handa na Python decoder, maaari mo itong pag-aralan . At ang mga may-ari ng mga SDR receiver ay maaaring mag-assemble at magpatakbo ng isang handa na ADS-B decoder , ito ay tinalakay nang mas detalyado sa .
Ang source code ng parser na inilarawan sa artikulo ay ibinigay sa ilalim ng cut. Isa itong halimbawa ng pagsubok na hindi nagpapanggap na produksyon, ngunit may gumagana dito, at maaari nilang i-parse ang file na naitala sa itaas.
Source Code (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()
Sana may interesado, salamat sa iyong pansin.
Pinagmulan: www.habr.com