Helo Habr. Mae'n debyg bod pawb sydd erioed wedi cyfarfod neu weld oddi ar berthnasau neu ffrindiau ar awyren wedi defnyddio'r gwasanaeth Flightradar24 rhad ac am ddim. Mae hon yn ffordd gyfleus iawn i olrhain lleoliad yr awyren mewn amser real.
В Disgrifiwyd egwyddor gweithredu gwasanaeth ar-lein o'r fath. Byddwn nawr yn mynd ymlaen i ddarganfod pa ddata sy'n cael ei anfon a'i dderbyn gan yr awyren i'r orsaf dderbyn a'i ddadgodio ein hunain gan ddefnyddio Python.
Stori
Yn amlwg, nid yw data awyrennau yn cael ei drosglwyddo i ddefnyddwyr ei weld ar eu ffonau smart. Gelwir y system yn ADS-B (gwyliadwriaeth ddibynnol awtomatig - darlledu), ac fe'i defnyddir i drosglwyddo gwybodaeth am yr awyren yn awtomatig i'r ganolfan reoli - mae ei dynodwr, cyfesurynnau, cyfeiriad, cyflymder, uchder a data arall yn cael eu trosglwyddo. Yn flaenorol, cyn dyfodiad systemau o'r fath, dim ond pwynt ar y radar y gallai'r anfonwr ei weld. Nid oedd hyn yn ddigon bellach pan oedd gormod o awyrennau.
Yn dechnegol, mae ADS-B yn cynnwys trosglwyddydd ar awyren sy'n anfon pecynnau gwybodaeth o bryd i'w gilydd ar amlder eithaf uchel o 1090 MHz (mae yna foddau eraill, ond nid oes gennym gymaint o ddiddordeb ynddynt, gan mai dim ond yma y trosglwyddir cyfesurynnau). Wrth gwrs, yn ychwanegol at y trosglwyddydd, mae yna hefyd dderbynnydd rhywle yn y maes awyr, ond i ni, fel defnyddwyr, mae ein derbynnydd ein hunain yn ddiddorol.
Gyda llaw, er mwyn cymharu, ymddangosodd y system gyntaf o'r fath, Airnav Radarbox, a ddyluniwyd ar gyfer defnyddwyr cyffredin, yn 2007, a chostiodd tua $900; mae tanysgrifiad i wasanaethau rhwydwaith yn costio $250 y flwyddyn arall.
Gellir darllen adolygiadau o'r perchnogion Rwsia cyntaf hynny ar y fforwm . Nawr bod derbynyddion RTL-SDR ar gael yn eang, gellir cydosod dyfais debyg am $30; roedd mwy am hyn yn . Gadewch i ni symud ymlaen at y protocol ei hun - gadewch i ni weld sut mae'n gweithio.
Derbyn signalau
Yn gyntaf, mae angen recordio'r signal. Dim ond 120 microseconds yw hyd y signal cyfan, felly i ddadosod ei gydrannau'n gyfforddus, mae derbynnydd SDR gydag amledd samplu o 5 MHz o leiaf yn ddymunol.
Ar ôl recordio, rydym yn derbyn ffeil WAV gyda chyfradd samplu o 5000000 o samplau yr eiliad; mae 30 eiliad o recordiad o'r fath yn “pwyso” tua 500MB. Mae gwrando arno gyda chwaraewr cyfryngau, wrth gwrs, yn ddiwerth - nid yw'r ffeil yn cynnwys sain, ond signal radio wedi'i ddigideiddio'n uniongyrchol - dyma'n union sut mae Radio Diffiniedig gan Feddalwedd yn gweithio.
Byddwn yn agor ac yn prosesu'r ffeil gan ddefnyddio Python. Gall y rhai sydd am arbrofi ar eu pen eu hunain lawrlwytho recordiad enghreifftiol .
Gadewch i ni lawrlwytho'r ffeil a gweld beth sydd y tu mewn.
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()
Canlyniad: gwelwn “guriadau” amlwg yn erbyn y sŵn cefndir.
Mae pob “pwls” yn signal, y mae ei strwythur i'w weld yn glir os cynyddwch y cydraniad ar y graff.
Fel y gwelwch, mae'r llun yn eithaf cyson â'r hyn a roddir yn y disgrifiad uchod. Gallwch ddechrau prosesu'r data.
Dadgodio
Yn gyntaf, mae angen i chi gael ychydig o ffrwd. Mae'r signal ei hun wedi'i amgodio gan ddefnyddio amgodio Manceinion:
O'r gwahaniaeth lefel mewn pigiadau mae'n hawdd cael “0” ac “1” go iawn.
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"
Mae strwythur y signal ei hun fel a ganlyn:
Gadewch i ni edrych ar y meysydd yn fwy manwl.
DF (Fformat Downlink, 5 did) - sy'n pennu'r math o neges. Mae yna sawl math:
()
Dim ond math DF17 sydd gennym ni ddiddordeb, oherwydd... Dyma sy'n cynnwys cyfesurynnau'r awyren.
ICAO (24 did) - cod unigryw rhyngwladol yr awyren. Gallwch wirio'r awyren yn ôl ei god (yn anffodus, mae'r awdur wedi rhoi'r gorau i ddiweddaru'r gronfa ddata, ond mae'n dal yn berthnasol). Er enghraifft, ar gyfer cod 3c5ee2 mae gennym y wybodaeth ganlynol:
Golygu: yn Rhoddir y disgrifiad o god ICAO yn fanylach; argymhellaf fod y rhai sydd â diddordeb yn ei ddarllen.
DATA (56 neu 112 did) - y data gwirioneddol y byddwn yn ei ddadgodio. Y 5 did cyntaf o ddata yw'r maes Cod Math, yn cynnwys is-fath y data sy'n cael ei storio (na ddylid ei gymysgu â DF). Mae yna ychydig iawn o'r mathau hyn:
()
Gadewch i ni edrych ar rai enghreifftiau o becynnau.
Adnabod awyrennau
Enghraifft ar ffurf ddeuaidd:
00100 011 000101 010111 000111 110111 110001 111000
Meysydd data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, mae pob nod C1-C8 yn cynnwys codau sy'n cyfateb i fynegeion yn y llinell:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ####_##################
Trwy ddadgodio'r llinyn, mae'n hawdd cael cod yr awyren: 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('#', ''))
Safle yn yr awyr
Os yw'r enw'n syml, yna mae'r cyfesurynnau'n fwy cymhleth. Fe'u trosglwyddir ar ffurf fframiau 2, eilrif ac od. Cod maes TC = 01011b = 11.
Enghraifft o becynnau eilrif ac od:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Mae cyfrifo cyfesurynnau ei hun yn digwydd yn ôl fformiwla eithaf anodd:
()
Dydw i ddim yn arbenigwr GIS, felly nid wyf yn gwybod o ble mae'n dod. Pwy a wyr, ysgrifennwch y sylwadau.
Mae uchder yn cael ei ystyried yn symlach - yn dibynnu ar y darn penodol, gellir ei gynrychioli naill ai fel lluosrif o 25 neu 100 troedfedd.
Cyflymder Awyr
Pecyn gyda TC=19. Y peth diddorol yma yw y gall y cyflymder fod naill ai'n gywir, o'i gymharu â'r ddaear (Ground Speed), neu yn yr awyr, wedi'i fesur gan synhwyrydd awyrennau (Airspeed). Mae llawer o feysydd gwahanol hefyd yn cael eu trosglwyddo:
()
Casgliad
Fel y gwelwch, mae technoleg ADS-B wedi dod yn symbiosis diddorol, pan fo safon yn ddefnyddiol nid yn unig i weithwyr proffesiynol, ond hefyd i ddefnyddwyr cyffredin. Ond wrth gwrs, chwaraewyd rhan allweddol yn hyn gan dechnoleg rhatach derbynwyr SDR digidol, sy'n caniatáu i'r ddyfais dderbyn signalau yn llythrennol ag amleddau uwchlaw gigahertz “am geiniogau.”
Yn y safon ei hun, wrth gwrs, mae llawer mwy. Gall y rhai sydd â diddordeb weld y PDF ar y dudalen neu ymwelwch â'r un a grybwyllwyd eisoes uchod .
Mae’n annhebygol y bydd pob un o’r uchod yn ddefnyddiol i lawer, ond o leiaf mae’r syniad cyffredinol o sut mae’n gweithio, rwy’n gobeithio, yn parhau.
Gyda llaw, mae datgodiwr parod yn Python eisoes yn bodoli, gallwch chi ei astudio . A gall perchnogion derbynwyr SDR ymgynnull a lansio datgodiwr ADS-B parod , trafodwyd hyn yn fanylach yn .
Rhoddir cod ffynhonnell y parser a ddisgrifir yn yr erthygl o dan y toriad. Mae hon yn enghraifft brawf nad yw'n esgus bod yn gynhyrchiad, ond mae rhai pethau'n gweithio ynddo, a gellir ei ddefnyddio i ddosrannu'r ffeil a gofnodwyd uchod.
Cod ffynhonnell (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()
Gobeithio bod gan rywun ddiddordeb, diolch am eich sylw.
Ffynhonnell: hab.com