Hei Habr. Luultavasti jokainen, joka on koskaan tavannut tai nähnyt sukulaisia tai ystäviä lentokoneessa, on käyttänyt ilmaista Flightradar24-palvelua. Tämä on erittäin kätevä tapa seurata lentokoneen sijaintia reaaliajassa.
В Tällaisen verkkopalvelun toimintaperiaate kuvattiin. Nyt mennään eteenpäin ja selvitetään, mitä dataa lähetetään ja vastaanotetaan lentokoneesta vastaanottoasemalle ja puretaan ne itse Pythonilla.
Tarina
Ilmeisesti lentokonetietoja ei välitetä käyttäjien nähtäväksi älypuhelimissaan. Järjestelmää kutsutaan nimellä ADS-B (Automatic dependent valve-broadcast), ja sen avulla lähetetään automaattisesti tietoja lentokoneesta lennonjohtokeskukseen - sen tunniste, koordinaatit, suunta, nopeus, korkeus ja muut tiedot lähetetään. Aikaisemmin, ennen tällaisten järjestelmien tuloa, lähettäjä saattoi nähdä vain pisteen tutkalla. Tämä ei enää riittänyt, kun lentokoneita oli liikaa.
Teknisesti ADS-B koostuu lentokoneen lähettimestä, joka lähettää ajoittain tietopaketteja melko suurella taajuudella 1090 MHz (on muitakin tiloja, mutta emme ole niin kiinnostuneita niistä, koska koordinaatit lähetetään vain täällä). Tietysti jossain lentokentällä on lähettimen lisäksi vastaanotin, mutta meille käyttäjille oma vastaanotin on mielenkiintoinen.
Muuten, vertailun vuoksi ensimmäinen tällainen tavallisille käyttäjille suunniteltu järjestelmä Airnav Radarbox ilmestyi vuonna 2007 ja maksoi noin 900 dollaria, verkkopalveluiden tilaus maksoi vielä 250 dollaria vuodessa.
Näiden ensimmäisten venäläisten omistajien arvostelut voidaan lukea foorumilla . Nyt kun RTL-SDR-vastaanottimia on tullut laajalti saataville, samanlainen laite voidaan koota 30 dollarilla; lisää tästä oli . Siirrytään itse protokollaan - katsotaan kuinka se toimii.
Signaalien vastaanotto
Ensinnäkin signaali on tallennettava. Koko signaalin kesto on vain 120 mikrosekuntia, joten sen komponenttien purkamiseksi mukavasti on toivottava SDR-vastaanotin, jonka näytteenottotaajuus on vähintään 5 MHz.
Tallennuksen jälkeen saamme WAV-tiedoston, jonka näytteenottotaajuus on 5000000 30 500 näytettä/s; XNUMX sekuntia tällaista tallennetta "painoa" noin XNUMX Mt. Sen kuunteleminen mediasoittimella on tietysti turhaa - tiedosto ei sisällä ääntä, vaan suoraan digitoitua radiosignaalia - juuri näin Software Defined Radio toimii.
Avaamme ja käsittelemme tiedoston Pythonilla. Ne, jotka haluavat kokeilla itse, voivat ladata esimerkkitallenteen .
Ladataan tiedosto ja katsotaan mitä sen sisällä on.
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()
Tulos: näemme ilmeisiä "pulsseja" taustamelua vasten.
Jokainen "pulssi" on signaali, jonka rakenne näkyy selvästi, jos lisäät kaavion resoluutiota.
Kuten näette, kuva on melko yhdenmukainen yllä olevan kuvauksen kanssa. Voit aloittaa tietojen käsittelyn.
Dekoodaus
Ensinnäkin sinun täytyy saada vähän virtaa. Itse signaali on koodattu Manchesterin koodauksella:
Nappuloiden tasoerosta on helppo saada oikeat "0" ja "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"
Itse signaalin rakenne on seuraava:
Katsotaanpa kenttiä tarkemmin.
DF (Downlink Format, 5 bits) - määrittää viestin tyypin. On olemassa useita tyyppejä:
()
Olemme kiinnostuneita vain tyypistä DF17, koska... Se sisältää lentokoneen koordinaatit.
ICAO (24 bittiä) - lentokoneen kansainvälinen yksilöllinen koodi. Voit tarkistaa koneen sen koodin perusteella (Valitettavasti kirjoittaja on lopettanut tietokannan päivittämisen, mutta se on edelleen ajankohtainen). Esimerkiksi koodia 3c5ee2 varten meillä on seuraavat tiedot:
Edit: sisään ICAO-koodin kuvaus on annettu tarkemmin, suosittelen kiinnostuneiden lukemaan sen.
DONNEES (56 tai 112 bittiä) - todelliset tiedot, jotka puretaan. Ensimmäiset 5 databittiä ovat kenttä Tyyppi koodi, joka sisältää tallennettavien tietojen alatyypin (ei pidä sekoittaa DF:ään). Näitä tyyppejä on useita:
()
Katsotaanpa muutamia esimerkkejä paketeista.
Lentokoneen tunnistus
Esimerkki binäärimuodossa:
00100 011 000101 010111 000111 110111 110001 111000
Tietokentät:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, jokainen merkki C1-C8 sisältää koodit, jotka vastaavat rivin indeksejä:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Dekoodaamalla merkkijono on helppo saada lentokonekoodi: 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('#', ''))
Asento ilmassa
Jos nimi on yksinkertainen, koordinaatit ovat monimutkaisempia. Ne lähetetään 2, parillisen ja parittoman kehyksen muodossa. Kenttäkoodi TC = 01011b = 11.
Esimerkki parillisista ja parittomista paketeista:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Itse koordinaattien laskenta tapahtuu melko hankalan kaavan mukaan:
()
En ole GIS-asiantuntija, joten en tiedä mistä se tulee. Kuka tietää, kirjoita kommentteihin.
Korkeutta pidetään yksinkertaisempana - tietystä bitistä riippuen se voidaan esittää joko 25 tai 100 jalan kerrannaisena.
Ilman nopeus
Paketti, jossa TC=19. Mielenkiintoista tässä on, että nopeus voi olla joko tarkka suhteessa maahan (Ground Speed) tai ilmassa, mitattuna lentokoneen anturilla (Airspeed). Myös monia erilaisia kenttiä lähetetään:
()
Johtopäätös
Kuten näette, ADS-B-teknologiasta on tullut mielenkiintoinen symbioosi, kun standardi on hyödyllinen paitsi ammattilaisille, myös tavallisille käyttäjille. Mutta tietysti avainrooli tässä oli digitaalisten SDR-vastaanottimien halvemmalla tekniikalla, jonka avulla laite voi kirjaimellisesti vastaanottaa signaaleja gigahertsin taajuuksilla "penneillä".
Itse standardissa on tietysti paljon muutakin. Kiinnostuneet voivat katsoa PDF-tiedoston sivulta tai vieraile jo edellä mainitussa .
On epätodennäköistä, että kaikki edellä mainitut ovat hyödyllisiä monille, mutta ainakin yleinen käsitys siitä, miten se toimii, toivottavasti säilyy.
Muuten, valmis dekooderi Pythonissa on jo olemassa, voit tutkia sitä . Ja SDR-vastaanottimien omistajat voivat koota ja käynnistää valmiin ADS-B-dekooderin , tästä keskusteltiin tarkemmin julkaisussa .
Artikkelissa kuvatun jäsentimen lähdekoodi on annettu leikkauksen alla. Tämä on testiesimerkki, joka ei teeskentele tuotantoa, mutta jotkut asiat toimivat siinä, ja sen avulla voidaan jäsentää yllä tallennettua tiedostoa.
Lähdekoodi (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()
Toivottavasti joku oli kiinnostunut, kiitos huomiosta.
Lähde: will.com