Hola Habr. Probablemente todos los que alguna vez han conocido o despedido a familiares o amigos en un avión han utilizado el servicio gratuito Flightradar24. Esta es una forma muy cómoda de seguir la posición del avión en tiempo real.
В Se describió el principio de funcionamiento de dicho servicio en línea. Ahora continuaremos y descubriremos qué datos se envían y reciben desde la aeronave a la estación receptora y los decodificaremos nosotros mismos usando Python.
historia
Obviamente, los datos de los aviones no se transmiten para que los usuarios los vean en sus teléfonos inteligentes. El sistema se llama ADS-B (vigilancia dependiente automática: transmisión) y se utiliza para transmitir automáticamente información sobre la aeronave al centro de control: se transmiten su identificador, coordenadas, dirección, velocidad, altitud y otros datos. Anteriormente, antes de la aparición de este tipo de sistemas, el despachador sólo podía ver un punto en el radar. Esto ya no era suficiente cuando había demasiados aviones.
Técnicamente, ADS-B consiste en un transmisor en un avión que envía periódicamente paquetes de información a una frecuencia bastante alta de 1090 MHz (hay otros modos, pero no nos interesan tanto, ya que las coordenadas solo se transmiten aquí). Por supuesto, además del transmisor, en algún lugar del aeropuerto también hay un receptor, pero para nosotros, como usuarios, nuestro propio receptor es interesante.
Por cierto, a modo de comparación, el primer sistema de este tipo, Airnav Radarbox, diseñado para usuarios comunes, apareció en 2007 y costaba alrededor de 900 dólares; una suscripción a los servicios de red costaba otros 250 dólares al año.
Las reseñas de aquellos primeros propietarios rusos se pueden leer en el foro. . Ahora que los receptores RTL-SDR están ampliamente disponibles, se puede montar un dispositivo similar por 30 dólares; más información sobre esto se encuentra en . Pasemos al protocolo en sí: veamos cómo funciona.
Recibir señales
Primero, es necesario grabar la señal. La señal completa tiene una duración de sólo 120 microsegundos, por lo que para desmontar cómodamente sus componentes es deseable un receptor SDR con una frecuencia de muestreo de al menos 5 MHz.
Después de la grabación, recibimos un archivo WAV con una velocidad de muestreo de 5000000 de muestras/segundo; 30 segundos de dicha grabación “pesan” alrededor de 500 MB. Escucharlo con un reproductor multimedia, por supuesto, es inútil: el archivo no contiene sonido, sino una señal de radio directamente digitalizada; así es exactamente como funciona Software Defined Radio.
Abriremos y procesaremos el archivo usando Python. Aquellos que quieran experimentar por su cuenta pueden descargar una grabación de ejemplo. .
Descarguemos 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 "pulsos" obvios contra el ruido de fondo.
Cada "pulso" es una señal, cuya estructura es claramente visible si aumenta la resolución del gráfico.
Como puede ver, la imagen coincide bastante con lo que se da en la descripción anterior. Puede comenzar a procesar los datos.
Descodificación
Primero, necesitas obtener un flujo de bits. La señal en sí está codificada mediante codificación Manchester:
A partir de la diferencia de nivel en los mordiscos es fácil obtener un “0” y un “1” reales.
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"
La estructura de la señal en sí es la siguiente:
Veamos los campos con más detalle.
DF (Formato de enlace descendente, 5 bits): determina el tipo de mensaje. Hay varios tipos:
()
Sólo nos interesa el tipo DF17, porque... Es este el que contiene las coordenadas del avión.
OACI (24 bits): código único internacional de la aeronave. Puedes consultar el avión por su código. (Desafortunadamente, el autor dejó de actualizar la base de datos, pero sigue siendo relevante). Por ejemplo, para el código 3c5ee2 tenemos la siguiente información:
Editar: en La descripción del código OACI se da con más detalle, recomiendo a los interesados leerlo.
DATOS (56 o 112 bits): los datos reales que decodificaremos. Los primeros 5 bits de datos son el campo. Código tipo, que contiene el subtipo de los datos que se almacenan (no confundir con DF). Hay bastantes de estos tipos:
()
Veamos algunos ejemplos de paquetes.
Identificación de aeronaves
Ejemplo en forma binaria:
00100 011 000101 010111 000111 110111 110001 111000
Campos de información:
+------+------+------+------+------+------+------+------+------+------+
| 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 contiene códigos correspondientes a índices en la línea:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Al decodificar la cadena, es fácil obtener el código de la 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 en el aire
Si el nombre es simple, las coordenadas son más complicadas. Se transmiten en forma de 2 fotogramas, pares e impares. Código de campo TC = 01011b = 11.
Ejemplo de paquetes pares e impares:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
El cálculo de las coordenadas en sí se realiza según una fórmula bastante complicada:
()
No soy un experto en SIG, así que no sé de dónde viene. Quién sabe, escribe en los comentarios.
La altura se considera más simple: dependiendo de la broca específica, se puede representar como un múltiplo de 25 o 100 pies.
Velocidad en el aire
Paquete con TC=19. Lo interesante aquí es que la velocidad puede ser precisa, en relación con el suelo (Ground Speed), o en el aire, medida por un sensor de avión (Airspeed). También se transmiten muchos campos diferentes:
()
Conclusión
Como puede ver, la tecnología ADS-B se ha convertido en una simbiosis interesante, cuando el estándar es útil no solo para los profesionales, sino también para los usuarios comunes. Pero, por supuesto, un papel clave lo desempeñó la tecnología más barata de los receptores SDR digitales, que permite que el dispositivo reciba literalmente señales con frecuencias superiores a un gigahercio "por unos centavos".
En la propia norma, por supuesto, hay mucho más. Los interesados pueden ver el PDF en la página o visitar el ya mencionado anteriormente .
Es poco probable que todo lo anterior sea útil para muchos, pero espero que al menos quede una idea general de cómo funciona.
Por cierto, ya existe un decodificador listo para usar en Python, puedes estudiarlo . Y los propietarios de receptores SDR pueden ensamblar y ejecutar un decodificador ADS-B listo para usar , esto se discutió con más detalle en .
El código fuente del analizador descrito en el artículo se encuentra debajo del corte. Este es un ejemplo de prueba que no pretende ser producción, pero algunas cosas funcionan en él y se puede usar para analizar el archivo registrado 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 interesado, gracias por su atención.
Fuente: habr.com