Sugeng Habr. Mbokmenawa saben wong sing wis tau ketemu utawa ndeleng sederek utawa kanca ing pesawat wis nggunakake layanan Flightradar24 gratis. Iki minangka cara sing trep banget kanggo nglacak posisi pesawat ing wektu nyata.
В Prinsip operasi layanan online kasebut diterangake. Saiki kita bakal nerusake lan ngerteni apa data sing dikirim lan ditampa saka pesawat menyang stasiun panampa lan dekoding dhewe nggunakake Python.
История
Temenan, data pesawat ora dikirim supaya pangguna bisa ndeleng ing smartphone. Sistem kasebut diarani ADS-B (Otomatis gumantung ndjogo-siaran), lan digunakake kanggo ngirim informasi kanthi otomatis babagan pesawat menyang pusat kontrol - pengenal sawijining, koordinat, arah, kacepetan, dhuwur lan data liyane ditularaké. Sadurunge, sadurunge tekane saka sistem kuwi, dispatcher mung bisa ndeleng titik ing radar. Iki ora cukup maneh nalika akeh banget pesawat.
Teknis, ADS-B kasusun saka pemancar ing pesawat sing periodik ngirim paket informasi ing frekuensi cukup dhuwur saka 1090 MHz (ana mode liyane, nanging kita ora dadi kasengsem ing wong-wong mau, amarga koordinat ditularaké mung kene). Mesthi, saliyane pemancar, ana uga panrima ing endi wae ing bandara, nanging kanggo kita, minangka pangguna, panrima kita dhewe menarik.
Miturut cara, kanggo mbandhingake, sistem pisanan kasebut, Airnav Radarbox, dirancang kanggo pangguna biasa, muncul ing taun 2007, lan regane udakara $900; langganan layanan jaringan regane $250 saben taun.
Ulasan saka pamilik Rusia sing sepisanan bisa diwaca ing forum kasebut . Saiki panrima RTL-SDR wis kasedhiya, piranti sing padha bisa dirakit kanthi rega $30; luwih akeh babagan iki ana ing . Ayo pindhah menyang protokol dhewe - ayo ndeleng cara kerjane.
Nampa sinyal
Kaping pisanan, sinyal kasebut kudu direkam. Sinyal kabeh nduweni durasi mung 120 mikrodetik, supaya bisa mbongkar komponen kanthi nyaman, panrima SDR kanthi frekuensi sampling paling sethithik 5 MHz.
Sawise ngrekam, kita nampa file WAV kanthi tingkat sampling 5000000 conto / detik; 30 detik rekaman kasebut "bobot" udakara 500MB. Ngrungokake karo pamuter media, mesthi, ora ana gunane - file ora ngemot swara, nanging sinyal radio langsung digitized - iki persis carane Software Defined Radio dianggo.
Kita bakal mbukak lan ngolah file nggunakake Python. Sing pengin nyoba dhewe bisa ngundhuh conto rekaman .
Ayo ngundhuh file lan ndeleng apa sing ana ing njero.
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()
Asil: kita ndeleng "pulsa" sing jelas nglawan gangguan latar mburi.
Saben "pulsa" minangka sinyal, struktur sing katon jelas yen sampeyan nambah resolusi ing grafik.
Kaya sing sampeyan ngerteni, gambar kasebut cukup konsisten karo apa sing diwenehake ing katrangan ing ndhuwur. Sampeyan bisa miwiti ngolah data.
Dekoding
Pisanan, sampeyan kudu njaluk stream dicokot. Sinyal kasebut dhewe dienkode nggunakake enkoding Manchester:
Saka prabédan tingkat ing nibbles iku gampang kanggo njaluk nyata "0" lan "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"
Struktur sinyal kasebut dhewe yaiku:
Ayo ndeleng lapangan kanthi luwih rinci.
DF (Format Downlink, 5 bit) - nemtokake jinis pesen. Ana sawetara jinis:
()
Kita mung kasengsem ing jinis DF17, amarga ... Iki sing ngemot koordinat pesawat.
ICAO (24 bit) - kode unik internasional pesawat. Sampeyan bisa mriksa pesawat kanthi kode (sayange, penulis wis mandheg nganyari database, nanging isih relevan). Contone, kanggo kode 3c5ee2 kita duwe informasi ing ngisor iki:
Sunting: ing Katrangan kode ICAO diwenehake kanthi luwih rinci; Aku nyaranake supaya sing kasengsem maca.
DATA (56 utawa 112 bit) - data nyata sing bakal kita decode. 5 bit pisanan data minangka lapangan Kode Ketik, ngemot subtipe data sing disimpen (ora bakal bingung karo DF). Ana sawetara jinis iki:
()
Ayo goleki sawetara conto paket.
Identifikasi pesawat
Tuladha ing wangun binar:
00100 011 000101 010111 000111 110111 110001 111000
Kolom data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, saben karakter C1-C8 ngemot kode sing cocog karo indeks ing baris:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######
Kanthi dekoding senar, gampang entuk kode pesawat: 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('#', ''))
Posisi udhara
Yen jeneng iku prasaja, banjur koordinat luwih rumit. Padha ditularaké ing wangun 2, pigura malah lan aneh. Kode lapangan TC = 01011b = 11.
Contoh paket genap dan ganjil:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Pitungan koordinat dhewe dumadi miturut rumus sing rada angel:
()
Aku dudu ahli GIS, mula aku ora ngerti asale saka ngendi. Sapa ngerti, tulis ing komentar.
Dhuwur dianggep luwih prasaja - gumantung saka bit tartamtu, bisa diwakili minangka kelipatan 25 utawa 100 kaki.
Kacepetan udhara
Paket karo TC = 19. Sing menarik ing kene yaiku kacepetan bisa akurat, relatif marang lemah (Ground Speed), utawa udhara, diukur dening sensor pesawat (Airspeed). Akeh lapangan sing beda-beda uga ditularake:
()
kesimpulan
Kaya sing sampeyan ngerteni, teknologi ADS-B wis dadi simbiosis sing menarik, nalika standar migunani ora mung kanggo para profesional, nanging uga kanggo pangguna biasa. Nanging mesthine, peran penting ing iki dimainake dening teknologi panrima SDR digital sing luwih murah, sing ngidini piranti kasebut kanthi harfiah nampa sinyal kanthi frekuensi ing ndhuwur gigahertz "kanggo dhuwit".
Ing standar dhewe, mesthi, ana luwih. Sing kasengsem bisa ndeleng PDF ing kaca kasebut utawa ngunjungi sing wis kasebut ing ndhuwur .
Ora mungkin kabeh sing kasebut ing ndhuwur bakal migunani kanggo akeh, nanging paling ora ide umum babagan cara kerjane, muga-muga tetep.
Miturut cara, decoder siap digawe ing Python wis ana, sampeyan bisa sinau . Lan pamilik panrima SDR bisa ngumpul lan miwiti dekoder ADS-B sing wis siap , iki dibahas luwih rinci ing .
Kode sumber parser sing diterangake ing artikel kasebut diwenehi ing ngisor potongan kasebut. Iki minangka conto tes sing ora nyamar dadi produksi, nanging ana sawetara perkara sing bisa digunakake, lan bisa digunakake kanggo ngurai file sing direkam ing ndhuwur.
Kode sumber (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()
Muga-muga ana sing kasengsem, matur nuwun kanggo perhatian sampeyan.
Source: www.habr.com