Olá Habr. Provavelmente, todos que já se despediram ou receberam amigos ou familiares no aeroporto utilizaram o serviço gratuito Flightradar24. É uma maneira bastante conveniente de acompanhar a localização do avião em tempo real.
Em O princípio de funcionamento desse tipo de serviço online foi descrito. Agora vamos mais longe e descobrir quais dados são transmitidos e recebidos da aeronave pela estação receptora, e como decodificá-los nós mesmos com o Python.
Histórico
É evidente que as informações sobre os aviões não são transmitidas para que os usuários as vejam em seus smartphones. O sistema é chamado de ADS–B (Automatic Dependent Surveillance–Broadcast) e serve para a transmissão automática de informações sobre a aeronave ao centro de controle — seu identificador, coordenadas, direção, velocidade, altura e outros dados são enviados. Anteriormente, antes do surgimento de tais sistemas, o despachante conseguia ver apenas um ponto no radar. Isso tornou-se insuficiente quando o número de aeronaves aumentou consideravelmente.
Tecnicamente, o ADS-B consiste em um transmissor na aeronave que periodicamente envia pacotes de informação em uma frequência bastante alta de 1090 MHz (existem outros modos, mas não são tão interessantes para nós, pois as coordenadas são transmitidas apenas aqui). Obviamente, além do transmissor, há um receptor em algum lugar no aeroporto, mas para nós, como usuários, o receptor que nos interessa é o nosso próprio.
A propósito, para comparação, o primeiro sistema desse tipo, o Airnav Radarbox, voltado para usuários comuns, surgiu em 2007 e custava cerca de 900 dólares; a assinatura para os serviços de rede custava mais 250 dólares por ano.
As avaliações dos primeiros proprietários russos podem ser lidas no fórum . Agora, quando os receptores RTL-SDR se tornaram amplamente disponíveis, um dispositivo semelhante pode ser montado por 30 dólares; mais detalhes sobre isso foram discutidos em . Vamos agora ao protocolo — vamos ver como isso funciona.
Recebendo sinais
Para começar, o sinal precisa ser gravado. Todo o sinal tem uma duração de apenas 120 microssegundos, portanto, para analisar suas componentes confortavelmente, é desejável ter um receptor SDR com frequência de amostragem de pelo menos 5 MHz.
Após a gravação, obtemos um arquivo WAV com uma taxa de amostragem de 5000000 samples/segundo; 30 segundos dessa gravação 'pesam' cerca de 500 MB. Ouvir isso em um reprodutor de mídia é, evidentemente, inútil — o arquivo contém não som, mas sim o sinal de rádio digitalizado — é assim que funciona o Software Defined Radio.
Vamos abrir e processar o arquivo usando Python. Aqueles que desejam experimentar por conta própria podem baixar um exemplo da gravação. .
Faremos o upload do arquivo e veremos o que há 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 'impulsos' claros em meio ao ruído.
Cada 'impulso' é o próprio sinal, cuja estrutura é bem visível ao ampliar a resolução no gráfico.
Como podemos ver, a imagem corresponde bem ao que foi descrito acima. Podemos começar a processar os dados.
Decodificação
Para começar, precisamos obter o fluxo de bits. O próprio sinal foi codificado usando codificação Manchester:
A partir da diferença de níveis em meio-byte, é fácil obter os reais '0' e '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"
A estrutura do próprio sinal é a seguinte:
Vamos analisar os campos com mais detalhes.
DF (Downlink Format, 5 bits) — define o tipo de mensagem. Existem vários tipos:
()
Estamos interessados apenas no tipo DF17, pois é este que contém as coordenadas da aeronave.
ICAO (24 bits) — código único internacional da aeronave. Você pode verificar a aeronave pelo seu código (infelizmente, o autor parou de atualizar a base, mas ela ainda é atual). Por exemplo, para o código 3c5ee2 temos as seguintes informações:
Correção: em a descrição do código ICAO é apresentada de forma mais detalhada. Recomendo que os interessados consultem.
DATA (56 ou 112 bits) — são os dados que iremos decodificar. Os primeiros 5 bits dos dados — campo Type Code, contém o subtipo dos dados armazenados (não confundir com DF). Existem muitos tipos:
()
Vamos analisar alguns exemplos de pacotes.
Identificação da aeronave
Exemplo em formato binário:
00100 011 000101 010111 000111 110111 110001 111000
Campos de dados:
+------+------+------+------+------+------+------+------+------+------+
| 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 contém códigos correspondentes aos índices na string:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######
Decodificando a string, não é difí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("Identificação da Aeronave:", code_str.replace('#', ''))
Posição em voo
Se o nome é simples, as coordenadas são mais complicadas. Elas são transmitidas em dois quadros, pares e ímpares. O código do campo TC = 01011b = 11.
Exemplo de pacotes par e ímpar:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
O cálculo das coordenadas é baseado em uma fórmula bastante complexa:
()
Eu não sou especialista em GIS, então não sei de onde isso é derivado. Quem souber, por favor, comente.
A altura é calculada de forma mais simples — dependendo de um bit específico, pode ser representada como múltipla de 25 ou 100 pés.
Velocidade em voo
Pacote com TC=19. O interessante aqui é que a velocidade pode ser tanto a exata, em relação ao solo (Ground Speed), quanto a aérea, medida pelo sensor da aeronave (Airspeed). Além disso, muitos campos diferentes são transmitidos:
()
Conclusão
Como pode ser visto, a tecnologia ADS-B se tornou um interessante simpósio, onde algum padrão é útil não apenas para os profissionais, mas também para os usuários comuns. Mas, claro, o fator chave nisso foi a redução de custos na tecnologia de receptores SDR digitais, permitindo que um dispositivo receba sinais a frequências acima de gigahercios por um preço irrisório.
No próprio padrão, há claro muito mais. Aqueles que desejarem podem conferir o PDF na página ou visitar o já mencionado .
Dificilmente tudo isso será útil para muitos, mas pelo menos a ideia geral de como isso funciona, espero, permaneceu.
A propósito, um decodificador pronto em Python já existe, e pode ser estudado . E os proprietários de receptores SDR podem montar e executar um decodificador de ADS-B pronto , detalhes sobre isso foram discutidos em .
O código-fonte do parser, descrito no artigo, está disponível abaixo. Este é um exemplo de teste, não pretende ser produção, mas algumas partes funcionam, e podem ser usadas para analisar o arquivo gravado acima.
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 tenha sido interessante para alguém, obrigado pela atenção.
Fonte: habr.com
