Flightradar24 - kumaha jalanna? Bagian 2, protokol ADS-B

Salam Habr. Panginten sadayana anu kantos tepang atanapi ningali baraya atanapi réréncangan dina pesawat parantos nganggo jasa Flightradar24 gratis. Ieu mangrupikeun cara anu saé pikeun ngalacak posisi pesawat sacara real waktos.

В bagian kahiji Prinsip operasi layanan online sapertos dijelaskeun. Urang ayeuna bakal teraskeun sareng terangkeun data naon anu dikirim sareng ditampi tina pesawat ka stasiun panampi sareng decoding sorangan nganggo Python.

dongeng

Jelas, data pesawat teu dikirimkeun pikeun pamaké pikeun nempo dina smartphone maranéhanana. Sistim nu disebut ADS-B (Otomatis gumantung panjagaan-siaran), sarta dipaké pikeun otomatis ngirimkeun informasi ngeunaan pesawat ka puseur kontrol - identifier na, koordinat, arah, speed, luhurna jeung data sejenna dikirimkeun. Saméméhna, saméméh mecenghulna sistem misalna, dispatcher nu ukur bisa ningali titik dina radar nu. Ieu henteu cekap deui nalika seueur teuing pesawat.

Téhnisna, ADS-B diwangun ku pamancar dina pesawat anu périodik ngirimkeun pakét inpormasi dina frékuénsi anu cukup luhur 1090 MHz (aya modus séjén, tapi urang teu jadi kabetot dina aranjeunna, saprak koordinat dikirimkeun ngan di dieu). Tangtosna, salian pamancar, aya ogé panarima dimana waé di bandara, tapi pikeun urang, salaku pangguna, panarima urang sorangan anu pikaresepeun.

Ngomong-ngomong, pikeun babandingan, sistem anu munggaran, Airnav Radarbox, dirancang pikeun pangguna biasa, muncul dina taun 2007, sareng hargana sakitar $ 900; langganan jasa jaringan ngarugikeun $ 250 deui sataun.

Ulasan ngeunaan pamilik Rusia anu munggaran tiasa dibaca dina forum radioscanner. Ayeuna yén panarima RTL-SDR parantos sayogi, alat anu sami tiasa dirakit pikeun $ 30; langkung seueur ngeunaan ieu dina bagian kahiji. Hayu urang ngalih ka protokol sorangan - hayu urang tingali kumaha gawéna.

Narima sinyal

Kahiji, sinyal perlu dirékam. Sakabéh sinyal boga durasi ngan 120 microseconds, jadi mun comfortably ngabongkar komponén na, panarima SDR kalawan frékuénsi sampling sahanteuna 5 MHz desirable.

Saatos ngarékam, kami nampi file WAV kalayan laju sampling 5000000 conto / detik; 30 detik rekaman sapertos "beurat" sakitar 500MB. Ngadangukeunana sareng pamuter média, tangtosna, henteu aya gunana - filena henteu ngandung sora, tapi sinyal radio langsung didigitalkeun - ieu persis kumaha jalanna Software Defined Radio.

Kami bakal muka sareng ngolah file nganggo Python. Anu hoyong ékspérimén nyalira tiasa ngaunduh conto rekaman link.

Hayu urang unduh file sareng tingali naon anu aya di jerona.

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()

Hasilna: urang ningali atra "pulsa" ngalawan noise latar.

Unggal "pulsa" mangrupikeun sinyal, strukturna jelas katingali upami anjeun ningkatkeun résolusi dina grafik.

Sakumaha anjeun tiasa tingali, gambar éta rada konsisten sareng anu dijelaskeun dina pedaran di luhur. Anjeun tiasa ngamimitian ngolah data.

Decoding

Kahiji, anjeun kudu meunang saeutik stream. Sinyal sorangan disandikeun nganggo encoding Manchester:

Tina bédana tingkat dina nibbles gampang pikeun meunangkeun nyata "0" jeung "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 sorangan nyaéta kieu:

Hayu urang nempo widang dina leuwih jéntré.

DF (Format Downlink, 5 bit) - nangtukeun jenis pesen. Aya sababaraha jinis:

(sumber méja)

Kami ngan ukur resep kana jinis DF17, sabab ... Ieu anu ngandung koordinat pesawat.

ICAO (24 bit) - kode unik internasional pesawat. Anjeun tiasa pariksa pesawat ku kode na online (hanjakalna, panulis geus dieureunkeun ngamutahirkeun database, tapi masih relevan). Salaku conto, pikeun kode 3c5ee2 kami gaduh inpormasi ieu:

Edit: di komentar kana artikel Katerangan ngeunaan kodeu ICAO dijelaskeun sacara langkung rinci; Abdi nyarankeun yén anu resep maca éta.

DATA (56 atanapi 112 bit) - data saleresna anu bakal kami decode. 5 bit data munggaran nyaéta lapangan Kodeu Ketik, ngandung subtipe data anu disimpen (henteu aya patalina sareng DF). Aya sababaraha jenis ieu:

(sumber méja)

Hayu urang tingali sababaraha conto bungkusan.

Idéntifikasi pesawat

Contona dina wangun binér:

00100 011 000101 010111 000111 110111 110001 111000

Widang data:

+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+

TC = 00100b = 4, unggal karakter C1-C8 ngandung kodeu pakait jeung indéks dina garis:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######

Ku decoding string, éta gampang pikeun meunangkeun 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 hawa

Upami namina sederhana, maka koordinatna langkung rumit. Éta dikirimkeun dina bentuk 2, pigura genap sareng ganjil. Kode widang TC = 01011b = 11.

Conto pakét genap sareng ganjil:

01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000

Itungan koordinat sorangan lumangsung nurutkeun rumus rada tricky:

(sumber)

Abdi sanés ahli GIS, janten kuring henteu terang timana asalna. Anu terang, tulis dina komentar.

Jangkungna dianggap leuwih basajan - gumantung kana bit husus, éta bisa digambarkeun boh mangrupa sababaraha 25 atawa 100 suku.

Laju Airborne

Pakét jeung TC=19. Hal anu pikaresepeun di dieu nyaéta lajuna tiasa akurat, relatif ka taneuh (Ground Speed), atanapi hawa udara, diukur ku sénsor pesawat (Airspeed). Loba widang béda ogé dikirimkeun:

(sumber)

kacindekan

Sakumaha anjeun tiasa tingali, téknologi ADS-B parantos janten simbiosis anu pikaresepeun, nalika standar mangpaat henteu ngan ukur pikeun profésional, tapi ogé pikeun pangguna biasa. Tapi tangtosna, peran konci dina ieu dimaénkeun ku téknologi langkung mirah tina panarima SDR digital, anu ngamungkinkeun alat pikeun sacara harfiah nampi sinyal kalayan frékuénsi luhur gigahertz "pikeun pennies".

Dina standar sorangan, tangtosna, aya leuwih. Anu kabetot tiasa ningali PDF dina halaman éta ICAO atawa nganjang ka nu geus disebutkeun di luhur website.

Henteu sigana yén sadaya di luhur bakal mangpaat pikeun seueur, tapi sahenteuna ideu umum kumaha éta jalanna, kuring ngarepkeun, tetep.

Ku jalan kitu, a decoder siap-dijieun di Python geus aya, anjeun tiasa diajar eta di dieu. Sareng anu gaduh panarima SDR tiasa ngumpul sareng ngaluncurkeun dekoder ADS-B anu siap-siap ti kaca, ieu dibahas dina leuwih jéntré dina bagian kahiji.

Kodeu sumber tina parser dijelaskeun dina artikel dirumuskeun di handap cut. Ieu mangrupikeun conto tés anu henteu nyamar janten produksi, tapi aya sababaraha hal anu tiasa dianggo, sareng éta tiasa dianggo pikeun ngémutan file anu kacatet di luhur.
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()

Kuring miharep batur éta kabetot, hatur nuhun kana perhatian Anjeun.

sumber: www.habr.com