Hola Habr. Probablemente, todos los que alguna vez han recibido o despedido a familiares o amigos en el aeropuerto han utilizado el servicio gratuito Flightradar24. Es una forma muy conveniente de rastrear la ubicación de un avión en tiempo real.
En se describió el principio de funcionamiento de este tipo de servicio en línea. Ahora vamos a avanzar y averiguar qué datos se transmiten y reciben del avión a la estación receptora, y los decodificaremos nosotros mismos con Python.
Historia
Es obvio que los datos sobre los aviones no se transmiten para que los usuarios los vean en sus teléfonos inteligentes. El sistema se llama ADS–B (Broadcast de vigilancia dependiente automática) y sirve para la transmisión automática de información sobre la aeronave al centro de control de tráfico aéreo—su identificador, coordenadas, dirección, velocidad, altitud y otros datos se transmiten. Antes, sin este tipo de sistemas, el controlador solo podía ver un punto en el radar. Esto se volvió insuficiente cuando hubo demasiados aviones.
Técnicamente, el ADS-B se compone de un transmisor en la aeronave que envía periódicamente paquetes con información en una frecuencia bastante alta de 1090 MHz (hay otros modos, pero no son tan interesantes para nosotros, ya que las coordenadas solo se transmiten aquí). Por supuesto, además del transmisor, hay un receptor en algún lugar del aeropuerto, pero para nosotros, como usuarios, nos interesa nuestro propio receptor.
Por cierto, para comparar, el primer sistema de este tipo, Airnav Radarbox, diseñado para usuarios comunes, apareció en 2007 y costaba alrededor de 900$, además de que la suscripción a los servicios de red costaba otros 250$ al año.
Las opiniones de esos primeros propietarios rusos se pueden leer en el foro . Ahora, cuando los receptores RTL-SDR están ampliamente disponibles, se puede ensamblar un dispositivo similar por 30$, más detalles sobre esto se encontraron en . Pasemos, entonces, al protocolo: veamos cómo funciona.
Recepción de señales
Para empezar, es necesario grabar la señal. La duración total de la señal es de solo 120 microsegundos, así que, para descomponer sus componentes cómodamente, se recomienda un receptor SDR con una frecuencia de muestreo de al menos 5 MHz.
Después de grabar, obtenemos un archivo WAV con una frecuencia de muestreo de 5000000 muestras/seg, 30 segundos de esta grabación "pesa" alrededor de 500 MB. Escucharlo con un reproductor multimedia es, por supuesto, inútil: el archivo no contiene sonido, sino la señal de radio digitalizada directamente; así es como funciona la Radio Definida por Software.
Abriremos y procesaremos el archivo usando Python. Quienes deseen experimentar por su cuenta, pueden descargar un ejemplo de grabación. .
Carguemos el archivo y veamos qué hay 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 claros "impulsos" en medio del ruido.
Cada "impulso" es la propia señal, cuya estructura se puede observar bien si aumentamos la resolución en el gráfico.
Como se puede ver, la imagen corresponde bastante a lo que se describe arriba. Se puede proceder al procesamiento de los datos.
Decodificación
Para comenzar, necesitamos obtener el flujo de bits. La señal misma está codificada utilizando el manchester encoding:
A partir de la diferencia de niveles en los medio bytes, es fácil obtener los reales "0" y "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 += "1"
La estructura de la señal misma es la siguiente:
Veamos los campos con más detalle.
DF (Downlink Format, 5 bits) — define el tipo de mensaje. Hay varios tipos:
()
Solo nos interesa el tipo DF17, ya que es el que contiene las coordenadas de la aeronave.
ICAO (24 bits) — código único internacional de la aeronave. Puedes verificar la aeronave por su código (desafortunadamente, el autor dejó de actualizar la base, pero todavía es actual). Por ejemplo, para el código 3c5ee2 tenemos la siguiente información:
Corrección: en la descripción del código ICAO se proporciona con más detalle; recomiendo a los interesados que se familiaricen con ello.
DATA (56 o 112 bits) — los propios datos, que decodificaremos. Los primeros 5 bits de datos son el campo Type Code, que contiene el subtipo de datos almacenados (no confundir con DF). Hay muchos tipos:
()
Analicemos algunos ejemplos de paquetes.
Identificación de aeronaves
Ejemplo 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 símbolo C1-C8 contiene códigos correspondientes a los índices en la cadena:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Decodificando la cadena, es fácil obtener el código del avión: 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("Identificación del avión:", code_str.replace('#', ''))
Posición en vuelo
Si el nombre es sencillo, las coordenadas son más complicadas. Se transmiten en forma de tramas pares e impares. El código del campo TC = 01011b = 11.
Ejemplo de paquetes par e impar:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
El cálculo de las coordenadas se realiza mediante una fórmula bastante ingeniosa:
()
No soy un especialista en GIS, así que no sé de dónde se obtiene. Quien esté al tanto, por favor, escriba en los comentarios.
La altitud se calcula de manera más sencilla: dependiendo de un bit específico, puede representarse como múltiplo de 25 o 100 pies.
Velocidad en vuelo
Paquete con TC=19. Lo interesante aquí es que la velocidad puede ser tanto exacta, relativa al suelo (Ground Speed), como en aire, medida por un sensor del avión (Airspeed). También se transmiten muchos otros campos diferentes:
()
Conclusión
Como se puede ver, la tecnología ADS-B se ha convertido en una interesante combinación, donde algún estándar es útil no solo para profesionales, sino también para usuarios comunes. Pero, por supuesto, un papel clave en esto ha sido el abaratamiento de la tecnología de receptores SDR digitales, que permiten recibir señales a frecuencias superiores a un gigaherzio de manera económica.
En el propio estándar, por supuesto, hay mucho más. Los interesados pueden consultar el PDF en la página o visitar el ya mencionado .
Es improbable que a muchos les sirva todo lo anterior, pero al menos espero que la idea general de cómo funciona haya quedado clara.
Por cierto, ya existe un decodificador terminado en Python, se puede estudiar . Y los propietarios de receptores SDR pueden ensamblar y poner en funcionamiento un decodificador ADS-B ya preparado , se habló más sobre esto en .
El código fuente del parser, descrito en el artículo, se incluye a continuación. Es un ejemplo de prueba que no pretende ser producción, pero algo en él funciona y puede utilizarse para analizar el archivo mencionado anteriormente.
Código fuente (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 a alguien le haya resultado interesante, gracias por su atención.
Fuente: habr.com
