Halo Habr. Mungkin semua orang yang pernah bertemu atau mengantar saudara atau teman di pesawat pernah menggunakan layanan gratis Flightradar24. Ini adalah cara yang sangat mudah untuk melacak posisi pesawat secara real time.
В Prinsip pengoperasian layanan online tersebut telah dijelaskan. Sekarang kita akan melanjutkan dan mencari tahu data apa yang dikirim dan diterima dari pesawat ke stasiun penerima dan mendekodekannya sendiri menggunakan Python.
Cerita
Jelas sekali, data pesawat tidak dikirimkan untuk dilihat pengguna di ponsel cerdas mereka. Sistem ini disebut ADS-B (Pengawasan ketergantungan otomatis—siaran), dan digunakan untuk mengirimkan informasi tentang pesawat secara otomatis ke pusat kendali - pengidentifikasi, koordinat, arah, kecepatan, ketinggian, dan data lainnya dikirimkan. Sebelumnya, sebelum munculnya sistem seperti itu, petugas operator hanya dapat melihat satu titik di radar. Ini tidak lagi cukup karena jumlah pesawat terlalu banyak.
Secara teknis, ADS-B terdiri dari pemancar di pesawat yang secara berkala mengirimkan paket informasi pada frekuensi yang cukup tinggi yaitu 1090 MHz (ada mode lain, tetapi kami tidak begitu tertarik padanya, karena koordinat hanya dikirimkan di sini). Tentunya selain transmitter, ada juga receiver di suatu tempat di bandara, namun bagi kami sebagai pengguna, receiver kami sendiri yang menarik.
Sebagai perbandingan, sistem pertama, Airnav Radarbox, dirancang untuk pengguna biasa, muncul pada tahun 2007, dan berharga sekitar $900; berlangganan layanan jaringan dikenakan biaya tambahan $250 per tahun.
Ulasan pemilik Rusia pertama dapat dibaca di forum . Kini setelah receiver RTL-SDR tersedia secara luas, perangkat serupa dapat dirakit seharga $30; informasi lebih lanjut mengenai hal ini dapat dilihat di . Mari beralih ke protokol itu sendiri - mari kita lihat cara kerjanya.
Menerima sinyal
Pertama, sinyal perlu direkam. Seluruh sinyal memiliki durasi hanya 120 mikrodetik, jadi untuk membongkar komponennya dengan nyaman, diperlukan penerima SDR dengan frekuensi sampling minimal 5 MHz.
Setelah perekaman, kami menerima file WAV dengan kecepatan pengambilan sampel 5000000 sampel/detik; rekaman 30 detik tersebut “berbobot” sekitar 500 MB. Mendengarkannya dengan pemutar media, tentu saja, tidak ada gunanya - file tersebut tidak berisi suara, tetapi sinyal radio yang langsung didigitalkan - begitulah cara kerja Software Defined Radio.
Kami akan membuka dan memproses file menggunakan Python. Bagi yang ingin bereksperimen sendiri dapat mendownload contoh rekamannya .
Ayo unduh filenya dan lihat isinya.
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()
Hasilnya: kita melihat “denyut” yang jelas terhadap kebisingan latar belakang.
Setiap “denyut” adalah sinyal, yang strukturnya terlihat jelas jika Anda meningkatkan resolusi pada grafik.
Seperti yang Anda lihat, gambarannya cukup sesuai dengan apa yang diberikan pada uraian di atas. Anda dapat mulai memproses data.
Penguraian kode
Pertama, Anda perlu mendapatkan sedikit aliran. Sinyal itu sendiri dikodekan menggunakan pengkodean Manchester:
Dari perbedaan level camilan, mudah untuk mendapatkan “0” dan “1” yang sebenarnya.
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 sinyalnya sendiri adalah sebagai berikut:
Mari kita lihat bidangnya lebih detail.
DF (Format Downlink, 5 bit) - menentukan jenis pesan. Ada beberapa jenis:
()
Kami hanya tertarik pada tipe DF17, karena... Di sinilah letak koordinat pesawat.
ICAO (24 bit) - kode unik internasional pesawat. Anda dapat memeriksa pesawat dengan kodenya (sayangnya penulis sudah berhenti mengupdate database, namun masih relevan). Misalnya, untuk kode 3c5ee2 kami memiliki informasi berikut:
Sunting: masuk Deskripsi kode ICAO diberikan lebih detail, saya anjurkan bagi yang berminat membacanya.
DATA (56 atau 112 bit) - data aktual yang akan kami dekode. 5 bit data pertama adalah bidang Jenis Kode, berisi subtipe data yang disimpan (jangan bingung dengan DF). Ada beberapa jenis ini:
()
Mari kita lihat beberapa contoh paket.
Identifikasi pesawat
Contoh dalam bentuk biner:
00100 011 000101 010111 000111 110111 110001 111000
Bidang data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, setiap karakter C1-C8 berisi kode yang sesuai dengan indeks pada baris:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Dengan memecahkan kode string tersebut, mudah untuk mendapatkan 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 mengudara
Kalau namanya sederhana, maka koordinatnya lebih rumit. Mereka ditransmisikan dalam bentuk 2 frame genap dan ganjil. Kode bidang TC = 01011b = 11.
Contoh paket genap dan ganjil:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Perhitungan koordinatnya sendiri dilakukan menurut rumus yang agak rumit:
()
Saya bukan ahli GIS, jadi tidak tahu dari mana asalnya. Siapa tahu, tulis di komentar.
Ketinggian dianggap lebih sederhana - tergantung pada bit tertentu, dapat direpresentasikan sebagai kelipatan 25 atau 100 kaki.
Kecepatan Lintas Udara
Paket dengan TC=19. Hal yang menarik di sini adalah kecepatannya bisa akurat, relatif terhadap darat (Ground Speed), atau di udara, diukur dengan sensor pesawat (Airspeed). Banyak bidang berbeda juga ditransmisikan:
()
Kesimpulan
Seperti yang Anda lihat, teknologi ADS-B telah menjadi simbiosis yang menarik, ketika sebuah standar bermanfaat tidak hanya bagi para profesional, tetapi juga bagi pengguna biasa. Namun tentu saja, peran penting dalam hal ini dimainkan oleh teknologi penerima SDR digital yang lebih murah, yang memungkinkan perangkat menerima sinyal dengan frekuensi di atas gigahertz “dengan harga yang sangat murah.”
Dalam standarnya sendiri tentunya masih banyak lagi. Bagi yang berminat dapat melihat PDF di halaman tersebut atau kunjungi yang sudah disebutkan di atas .
Tidak mungkin semua hal di atas akan bermanfaat bagi banyak orang, tetapi setidaknya gambaran umum tentang cara kerjanya, saya harap, tetap ada.
Omong-omong, decoder siap pakai dengan Python sudah ada, Anda bisa mempelajarinya . Dan pemilik receiver SDR dapat merakit dan meluncurkan decoder ADS-B yang sudah jadi , ini dibahas lebih detail di .
Kode sumber parser yang dijelaskan dalam artikel diberikan di bawah potongan. Ini adalah contoh pengujian yang tidak berpura-pura menjadi produksi, tetapi ada beberapa hal yang berfungsi di dalamnya, dan dapat digunakan untuk mengurai file yang direkam di atas.
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()
Saya harap ada yang tertarik, terima kasih atas perhatian Anda.
Sumber: www.habr.com