Ola Habr. Probablemente, todos os que coñeceron ou viran a familiares ou amigos nun avión utilizaron o servizo gratuíto Flightradar24. Esta é unha forma moi cómoda de rastrexar a posición da aeronave en tempo real.
В Describiuse o principio de funcionamento deste servizo en liña. Agora imos adiante e descubrir que datos se envían e reciben desde a aeronave á estación receptora e decodificaremos nós mesmos usando Python.
Historia
Obviamente, os datos da aeronave non se transmiten para que os usuarios os vexan nos seus teléfonos intelixentes. O sistema chámase ADS-B (Automatic dependent surveillance—broadcast) e úsase para transmitir automaticamente información sobre a aeronave ao centro de control: transmítense o seu identificador, coordenadas, dirección, velocidade, altitude e outros datos. Anteriormente, antes da aparición deste tipo de sistemas, o despachador só podía ver un punto no radar. Isto xa non era suficiente cando había demasiados avións.
Tecnicamente, ADS-B consiste nun transmisor nunha aeronave que envía periodicamente paquetes de información a unha frecuencia bastante alta de 1090 MHz (hai outros modos, pero non nos interesan tanto, xa que as coordenadas só se transmiten aquí). Por suposto, ademais do transmisor, tamén hai un receptor nalgún lugar do aeroporto, pero para nós, como usuarios, é interesante o noso propio receptor.
Por certo, a modo de comparación, o primeiro sistema deste tipo, Airnav Radarbox, deseñado para usuarios comúns, apareceu en 2007 e custaba uns 900 dólares; unha subscrición aos servizos de rede custaba outros 250 dólares ao ano.
No foro pódense ler as opinións dos primeiros propietarios rusos . Agora que os receptores RTL-SDR están amplamente dispoñibles, pódese montar un dispositivo similar por $ 30; máis sobre isto foi en . Pasemos ao protocolo en si, vexamos como funciona.
Recepción de sinais
En primeiro lugar, o sinal debe ser gravado. Todo o sinal ten unha duración de só 120 microsegundos, polo que para desmontar comodamente os seus compoñentes, é desexable un receptor SDR cunha frecuencia de mostraxe de polo menos 5 MHz.
Despois da gravación, recibimos un ficheiro WAV cunha frecuencia de mostraxe de 5000000 de mostras/seg; 30 segundos desa gravación "pesan" uns 500 MB. Escoitalo cun reprodutor multimedia, por suposto, é inútil: o ficheiro non contén son, senón un sinal de radio dixitalizado directamente, así é exactamente como funciona Software Defined Radio.
Abriremos e procesaremos o ficheiro usando Python. Os que queiran experimentar pola súa conta poden descargar un exemplo de gravación .
Descarguemos o ficheiro e vexamos o que hai dentro.
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()
Resultado: vemos "pulsos" obvios contra o ruído de fondo.
Cada "pulso" é un sinal, cuxa estrutura é claramente visible se aumenta a resolución no gráfico.
Como podes ver, a imaxe é bastante consistente co que se dá na descrición anterior. Podes comezar a procesar os datos.
Decodificación
En primeiro lugar, cómpre obter un pouco de fluxo. O sinal en si está codificado usando a codificación Manchester:
A partir da diferenza de nivel dos bocados é fácil obter "0" e "1" reais.
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"
A estrutura do sinal é a seguinte:
Vexamos os campos con máis detalle.
DF (Formato de enlace descendente, 5 bits): determina o tipo de mensaxe. Hai varios tipos:
()
Só nos interesa o tipo DF17, porque... É este o que contén as coordenadas da aeronave.
ICAO (24 bits) - código único internacional da aeronave. Podes comprobar o avión polo seu código (Desafortunadamente, o autor deixou de actualizar a base de datos, pero segue sendo relevante). Por exemplo, para o código 3c5ee2 temos a seguinte información:
Edición: en A descrición do código ICAO dáse con máis detalle; recoméndolle aos interesados que o lean.
DATOS (56 ou 112 bits): os datos reais que descodificaremos. Os primeiros 5 bits de datos son o campo Código de tipo, que contén o subtipo dos datos que se almacenan (non se debe confundir con DF). Hai bastantes destes tipos:
()
Vexamos algúns exemplos de paquetes.
Identificación da aeronave
Exemplo en forma binaria:
00100 011 000101 010111 000111 110111 110001 111000
Campos de datos:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, cada carácter C1-C8 contén códigos correspondentes aos índices da liña:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Ao decodificar a cadea, é fácil obter o código da aeronave: 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('#', ''))
Posición aérea
Se o nome é sinxelo, as coordenadas son máis complicadas. Transmítense en forma de 2 fotogramas pares e impares. Código de campo TC = 01011b = 11.
Exemplo de paquetes pares e impares:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
O cálculo das coordenadas en si ocorre segundo unha fórmula bastante complicada:
()
Non son un experto en GIS, polo que non sei de onde vén. Quen sabe, escribe nos comentarios.
A altura considérase máis sinxela: dependendo do bit específico, pódese representar como un múltiplo de 25 ou 100 pés.
Velocidade aerotransportada
Paquete con TC=19. O interesante aquí é que a velocidade pode ser precisa, relativa ao chan (Ground Speed) ou aerotransportada, medida por un sensor de avión (Airspeed). Tamén se transmiten moitos campos diferentes:
()
Conclusión
Como podes ver, a tecnoloxía ADS-B converteuse nunha simbiose interesante, cando un estándar é útil non só para profesionais, senón tamén para usuarios comúns. Pero, por suposto, un papel fundamental neste xogo foi a tecnoloxía máis barata dos receptores SDR dixitais, que permite que o dispositivo reciba literalmente sinais con frecuencias superiores a un gigahercio "por centavos".
Na propia norma, por suposto, hai moito máis. Os interesados poden ver o PDF na páxina ou visita o xa mencionado anteriormente .
É improbable que todo o anterior sexa útil para moitos, pero polo menos a idea xeral de como funciona, espero, permanece.
Por certo, xa existe un decodificador listo en Python, podes estudalo . E os propietarios de receptores SDR poden montar e lanzar un descodificador ADS-B preparado , isto foi discutido con máis detalle en .
O código fonte do analizador descrito no artigo aparece debaixo do corte. Este é un exemplo de proba que non pretende ser produción, pero algunhas cousas funcionan nel e pódese usar para analizar o ficheiro gravado anteriormente.
Código fonte (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()
Espero que alguén estivese interesado, grazas pola vosa atención.
Fonte: www.habr.com