🚋 💆🏿 😡 Flightradar24 - como funciona? Parte 2, Protocolo ADS-B 🤵🏻 💓 🤽🏽

Oi Habr. Provavelmente todos que já conheceram ou acompanharam parentes ou amigos em um avião usaram o serviço gratuito Flightradar24. Essa é uma maneira muito conveniente de rastrear a posição da aeronave em tempo real.

imagem

A primeira parte descreveu o princípio de operação de um serviço on-line. Agora iremos além, descobriremos quais dados são transmitidos e recebidos da aeronave para a estação receptora e decodificá-los independentemente usando Python.

A história

Obviamente, os dados sobre a aeronave não são transmitidos para que os usuários possam vê-los em seus smartphones. O sistema é chamado ADS - B (vigilância dependente automática - transmissão) e é usado para transmitir automaticamente informações da aeronave a um centro de controle - seu identificador, coordenadas, direção, velocidade, altitude e outros dados são transmitidos. Anteriormente, antes do advento de tais sistemas, o despachante via apenas um ponto no radar. Isso não foi suficiente quando havia muitos aviões.

Tecnicamente, o ADS-B consiste em um transmissor em uma aeronave que envia periodicamente pacotes com informações a uma frequência bastante alta de 1090 MHz (existem outros modos, mas eles não são tão interessantes para nós, porque as coordenadas são transmitidas apenas aqui). É claro que, além do transmissor, também há um receptor em algum lugar do aeroporto, mas para nós, assim como para os usuários, nosso próprio receptor é interessante.

A propósito, para comparação, o primeiro sistema desse tipo, o Airnav Radarbox, projetado para usuários comuns, apareceu em 2007 e custou cerca de US $ 900, cerca de US $ 250 por ano, valendo uma assinatura dos serviços de rede.

Você pode ler as resenhas desses primeiros proprietários russos no fórum do radioscanner . Agora que os receptores RTL-SDR estão amplamente disponíveis, um dispositivo semelhante pode ser montado por US $ 30, mais sobre isso na primeira parte . Vamos prosseguir para o próprio protocolo - vamos ver como ele funciona.

Recebendo sinais

Para começar, o sinal precisa ser gravado. Todo o sinal tem uma duração de apenas 120 microssegundos, portanto, para desmontar confortavelmente seus componentes, é desejável um receptor SDR com uma frequência de amostragem de pelo menos 5 MHz.

imagem

Após a gravação, obtemos um arquivo WAV com uma frequência de amostragem de 5.000.000 amostras / s, 30 segundos dessa gravação “pesam” cerca de 500 MB. Ouvi-lo com um media player é, obviamente, inútil - o arquivo não contém som, mas um sinal de rádio digitalizado diretamente - é assim que o Rádio Definido por Software funciona.

Vamos abrir e processar o arquivo usando Python. Aqueles que desejam experimentar por conta própria podem fazer o download da amostra de gravação no link .

Baixe o arquivo e veja 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 “pulsos” óbvios no contexto do ruído.

Cada "impulso" é um sinal cuja estrutura é claramente visível se você aumentar a resolução no gráfico.

Como você pode ver, a imagem é consistente com o descrito na descrição acima. Você pode começar a processar os dados.

Decodificação

Primeiro você precisa obter um fluxo de bits. O sinal em si é codificado usando a codificação manchester:

A partir da diferença de níveis em mordidelas, é 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 próprio sinal é a seguinte:

Vamos considerar os campos com mais detalhes.

DF (formato de downlink, 5 bits) - define o tipo de mensagem. Existem vários tipos deles:

( origem da tabela )

Estamos interessados apenas no tipo DF17, como contém as coordenadas da aeronave.

ICAO (24 bits) é um código internacional exclusivo de aeronave. Você pode verificar a aeronave por seu código no site (infelizmente, o autor parou de atualizar o banco de dados, mas ainda é relevante). Por exemplo, para o código 3c5ee2, temos as seguintes informações:

Editar: no comentário do artigo, a descrição do código da ICAO é fornecida com mais detalhes; eu recomendo que você se familiarize com os interessados.

DADOS (56 ou 112 bits) - na verdade, os dados que iremos decodificar. Os primeiros 5 bits de dados são o campo Código de tipo que contém o subtipo dos dados armazenados (não deve ser confundido com o DF). Existem alguns desses tipos:

( origem da tabela )

Vejamos alguns exemplos de pacotes.

Identificação da aeronave

Exemplo 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 caractere C1-C8 contém códigos correspondentes aos índices na sequência:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ ##### _ ################# 0123456789 ######

Após decodificar a linha, é fácil obter o código do avião: 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('#', ''))

Posição no ar

Se o nome for simples, as coordenadas serão mais complicadas. Eles são transmitidos na forma de 2x, quadros pares e ímpares. Código do campo TC = 01011b = 11.

Exemplo de pacotes pares e ímpares:

 01011 000 000101110110 00 10111000111001000 10000110101111001 01011 000 000110010000 01 10010011110000110 10000011110001000

O cálculo das coordenadas ocorre de acordo com uma fórmula bastante inteligente:

( fonte )

Como não sou especialista em GIS, não sei de onde vem. Quem sabe, escreva nos comentários.

A altitude é considerada mais fácil - dependendo de um bit específico, ela pode aparecer como um múltiplo de 25 ou 100 pés.

Velocidade no ar

Pacote com TC = 19. O interessante aqui é que a velocidade pode ser precisa, em relação ao solo (velocidade no solo) ou no ar, conforme medido por um sensor de avião (velocidade do ar). Muitos campos diferentes também são transmitidos:

( fonte )

Conclusão

Como você pode ver, a tecnologia ADS-B se tornou uma simbiose interessante quando um padrão é útil não apenas para profissionais, mas também para usuários comuns. Mas é claro que o papel principal disso foi desempenhado pelo barateamento da tecnologia dos receptores digitais SDR, que permite que o dispositivo receba sinais com uma frequência acima de gigahertz literalmente "por um centavo".

No próprio padrão, é claro, muito mais que tudo. Os interessados podem ver o PDF na página da ICAO ou visitar o site já mencionado acima.

É improvável que muitos dos itens acima sejam úteis, mas pelo menos a idéia geral de como isso funciona, espero, permanece.

A propósito, já existe um decodificador pronto em Python, que pode ser estudado aqui . E os proprietários dos receptores SDR podem montar e executar o decodificador ADS-B finalizado a partir da página , mais sobre isso na primeira parte .

O código-fonte do analisador descrito no artigo é fornecido sob o corte. Este é um exemplo de teste que não pretende ser produção, mas algo funciona nele, e você pode 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: %fs" % (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ém esteja interessado, obrigado por sua atenção.

Flightradar24 - como funciona? Parte 2, Protocolo ADS-B

A história

Recebendo sinais

Decodificação

Conclusão

More articles: