Sæll Habr. Sennilega hafa allir sem einhvern tíma hitt eða séð ættingja eða vini í flugvél notað ókeypis Flightradar24 þjónustuna. Þetta er mjög þægileg leið til að fylgjast með staðsetningu flugvélarinnar í rauntíma.
В Lýst var rekstrarreglu slíkrar netþjónustu. Við munum nú halda áfram og reikna út hvaða gögn eru send og móttekin frá flugvélinni til móttökustöðvarinnar og afkóða þau sjálf með Python.
Story
Augljóslega eru flugvélagögn ekki send til notenda að sjá á snjallsímum sínum. Kerfið er kallað ADS-B (Automatic dependent surveillance—broadcast), og er notað til að senda sjálfkrafa upplýsingar um flugvélina til stjórnstöðvar - auðkenni hennar, hnit, stefnu, hraði, hæð og önnur gögn eru send. Áður fyrr, áður en slík kerfi komu til sögunnar, gat sendimaðurinn aðeins séð punkt á ratsjánni. Þetta var ekki lengur nóg þegar það voru of margar flugvélar.
Tæknilega séð samanstendur ADS-B af sendi á flugvél sem sendir reglulega pakka af upplýsingum á nokkuð hári tíðni 1090 MHz (það eru aðrar stillingar, en við höfum ekki svo mikinn áhuga á þeim, þar sem hnit eru send aðeins hér). Auk sendisins er auðvitað líka móttakari einhvers staðar á flugvellinum, en fyrir okkur sem notendur er okkar eigin móttakari áhugaverður.
Við the vegur, til samanburðar, fyrsta slíka kerfið, Airnav Radarbox, hannað fyrir venjulega notendur, kom fram árið 2007 og kostaði um $900; áskrift að sérþjónustu kostaði aðra $250 á ári.
Umsagnir um þessa fyrstu rússnesku eigendur má lesa á spjallborðinu . Nú þegar RTL-SDR móttakarar eru orðnir víða fáanlegir er hægt að setja saman svipað tæki fyrir $30; meira um þetta var í . Við skulum halda áfram að samskiptareglunni sjálfri - við skulum sjá hvernig hún virkar.
Að taka á móti merkjum
Fyrst þarf að skrá merkið. Allt merkið er aðeins 120 míkrósekúndur, svo til að taka íhluti þess í sundur á þægilegan hátt er SDR móttakari með sýnatökutíðni að minnsta kosti 5 MHz æskilegt.
Eftir upptöku fáum við WAV skrá með sýnatökutíðni upp á 5000000 sýni/sek; 30 sekúndur af slíkri upptöku „vega“ um 500MB. Að hlusta á það með margmiðlunarspilara er auðvitað gagnslaust - skráin inniheldur ekki hljóð, heldur beint stafrænt útvarpsmerki - þetta er nákvæmlega hvernig Software Defined Radio virkar.
Við munum opna og vinna úr skránni með Python. Þeir sem vilja gera tilraunir á eigin spýtur geta sótt dæmi um upptöku .
Við skulum hlaða niður skránni og sjá hvað er inni.
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()
Niðurstaða: við sjáum augljósa „púlsa“ á móti bakgrunnshljóði.
Hver „púls“ er merki, uppbygging þess sést vel ef þú eykur upplausnina á línuritinu.
Eins og þú sérð er myndin alveg í samræmi við það sem gefið er upp í lýsingunni hér að ofan. Þú getur byrjað að vinna úr gögnunum.
Afkóðun
Fyrst þarftu að fá smá straum. Merkið sjálft er kóðað með Manchester kóðun:
Út frá stigsmuninum á nibblum er auðvelt að fá raunverulegt „0“ og „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"
Uppbygging merkisins sjálfs er sem hér segir:
Við skulum skoða reitina nánar.
DF (Downlink Format, 5 bitar) - ákvarðar tegund skilaboða. Það eru nokkrar gerðir:
()
Við höfum aðeins áhuga á gerð DF17, vegna þess að... Það er þetta sem inniheldur hnit flugvélarinnar.
ICAO (24 bitar) - alþjóðlegur einstakur kóði loftfarsins. Þú getur athugað flugvélina með kóða hennar (því miður er höfundur hætt að uppfæra gagnagrunninn en hann á samt við). Til dæmis, fyrir kóðann 3c5ee2 höfum við eftirfarandi upplýsingar:
Edit: í Lýsingin á ICAO kóðanum er gefin nánar, ég mæli með því að áhugasamir lesi hann.
GÖGN (56 eða 112 bita) - raunveruleg gögn sem við munum afkóða. Fyrstu 5 gagnabitarnir eru reiturinn Tegundarkóði, sem inniheldur undirtegund gagna sem verið er að geyma (ekki að rugla saman við DF). Það eru til nokkrar af þessum tegundum:
()
Við skulum skoða nokkur dæmi um pakka.
Auðkenning loftfars
Dæmi í tvöfaldri mynd:
00100 011 000101 010111 000111 110111 110001 111000
Gagnareitir:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, hver stafur C1-C8 inniheldur kóða sem samsvara vísitölum í línunni:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Með því að afkóða strenginn er auðvelt að fá flugvélakóðann: 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('#', ''))
Loftborinn staða
Ef nafnið er einfalt, þá eru hnitin flóknari. Þeir eru sendar í formi 2, sléttra og stakra ramma. Reitkóði TC = 01011b = 11.
Dæmi um jafna og staka pakka:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Útreikningur á hnitum fer sjálfur fram samkvæmt frekar erfiðri formúlu:
()
Ég er ekki GIS sérfræðingur, svo ég veit ekki hvaðan það kemur. Hver veit, skrifaðu í athugasemdir.
Hæð er talin einfaldari - allt eftir tilteknum bita er hægt að tákna hana sem annað hvort margfeldi af 25 eða 100 fetum.
Flughraði
Pakki með TC=19. Það áhugaverða hér er að hraðinn getur annað hvort verið nákvæmur, miðað við jörðu (Ground Speed), eða í lofti, mældur með flugskynjara (Airspeed). Mörg mismunandi svið eru einnig send:
()
Ályktun
Eins og þú sérð hefur ADS-B tækni orðið áhugavert sambýli, þegar staðall er gagnlegur ekki aðeins fyrir fagfólk heldur einnig venjulegum notendum. En auðvitað var lykilhlutverk í þessu gegnt af ódýrari tækni stafrænna SDR móttakara, sem gerir tækinu kleift að taka á móti merki með tíðni yfir gígahertz "fyrir smáaura."
Í staðlinum sjálfum er auðvitað miklu meira. Áhugasamir geta skoðað PDF á síðunni eða heimsækja þann sem þegar er nefndur hér að ofan .
Það er ólíklegt að allt ofangreint muni nýtast mörgum, en að minnsta kosti er almenn hugmynd um hvernig það virkar, ég vona, áfram.
Við the vegur, tilbúinn afkóðari í Python er þegar til, þú getur rannsakað það . Og eigendur SDR móttakara geta sett saman og sett af stað tilbúinn ADS-B afkóðara , var fjallað um þetta nánar í .
Frumkóði þáttarans sem lýst er í greininni er gefinn fyrir neðan klippið. Þetta er prufudæmi sem gefur sig ekki út fyrir að vera framleiðsla, en sumt virkar í því og það er hægt að nota til að flokka skrána sem skráð er hér að ofan.
Upprunakóði (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()
Ég vona að einhver hafi haft áhuga, takk fyrir athyglina.
Heimild: www.habr.com