Geekly articles weekly

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

Vou adivinhar e dizer que todos cujos amigos ou familiares já viajaram de avião usaram o Flightradar24 - um serviço gratuito e conveniente para rastrear voos em tempo real.

imagem

Na primeira parte, as idéias básicas de operação foram descritas. Agora vamos mais longe e descobrir que dados estão transmitindo e recebendo exatamente entre a aeronave e uma estação terrestre. Também decodificaremos esses dados usando Python.

História


Deveria ser óbvio que os dados dos aviões não se destinam apenas aos usuários que desejam vê-los em seus smartphones. Este sistema é chamado ADS-B (vigilância dependente automática - transmissão) e foi projetado para a transmissão automática dos dados de informações da aeronave para o centro de controle - diferentes parâmetros, como coordenadas, velocidade, rumo, altitude e outros, estão enviando. Anteriormente, o expedidor podia ver apenas um ponto na tela do radar. E definitivamente não se tornou suficiente quando o número de aviões aumentou drasticamente.

Tecnicamente, o ADS-B consiste em um transmissor dentro da aeronave, que envia periodicamente quadros de dados de informações a uma frequência relativamente 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 somente aqui). ) Claro, também há um receptor em algum lugar do aeroporto, mas para nós, assim como para os usuários, nosso próprio receptor é mais interessante.

Apenas para comparação, o primeiro sistema projetado para usuários comuns, o Airnav Radarbox, apareceu em 2007 e custou cerca de US $ 900, e cerca de US $ 250 por ano custaram uma assinatura de seus serviços de rede (a principal idéia desse sistema é coletar e compartilhar dados de muitos receptores , um receptor independente é relativamente inútil).



Agora, quando os receptores RTL-SDR se tornam muito mais disponíveis, um dispositivo semelhante pode ser fabricado por US $ 30. Ele pode ser encontrado na primeira parte do artigo , e iremos além e descreveremos o próprio protocolo - vamos ver como ele funciona.

Recebendo um sinal


Primeiro, precisamos gravar uma amostra de sinal. Todo o sinal tem apenas 120 microssegundos de duração e, para ver seus detalhes em uma "resolução" boa, é melhor ter rádio SDR com taxa de amostragem de pelo menos 5 MHz.

imagem

Após a gravação, estamos obtendo um arquivo WAV com uma taxa de amostragem de 5.000.000 de amostras / s, 30 segundos desse registro tem aproximadamente 500 MB de tamanho. Obviamente, é inútil ouvi-lo com seu media player favorito - o arquivo não contém um som, mas um sinal de rádio diretamente digitalizado - é exatamente assim que o Rádio Definido por Software funciona.

Podemos abrir e processar esse arquivo com Python. Aqueles que desejam fazer o experimento por conta própria, podem baixar uma amostra deste link .

Vamos abrir um arquivo e ver o que está 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 alguns 'impulsos' sobre o ruído.



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



Como podemos ver, a imagem é totalmente consistente com a descrição acima. Agora podemos processar esses dados.

Decodificação


Primeiro, precisamos ter um pouco de fluxo. O sinal em si é codificado com uma codificação manchester:



A partir das diferenças de meia mordida, podemos facilmente 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 em si é assim:



Vamos ver os campos com mais detalhes.

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


( origem da página )

Estamos interessados ​​apenas no tipo DF17, porque apenas este contém as coordenadas da aeronave.

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



DATA (56 ou 112 bits) - são os dados em si, que decodificaremos. Os primeiros 5 bits de dados são o campo Código de tipo , que contém o subtipo dos dados que estão sendo armazenados (não devem ser misturados ao campo DF ). Existem muitos desses tipos:


( origem da tabela )

Vejamos alguns exemplos.

Identificação da aeronave

Um 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, e cada símbolo C1-C8 contém códigos, que devem corresponder aos índices nesta cadeia de caracteres:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ ##### _ ################# 0123456789 ######

Após decodificar, é fácil obter o nome 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('#', ''))

Posição no ar

A decodificação de nomes era simples, mas as coordenadas são mais complicadas. Eles são transmitidos na forma de 2 quadros, pares e ímpares. O código do campo TC = 01011b = 11.



Exemplo de quadros de dados pares e ímpares:

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

O cálculo das coordenadas em si usa uma fórmula um pouco complicada:


( fonte )

Não sou especialista em GIS, portanto não sei de onde vem. Quem o conhece melhor, por favor escreva nos comentários.

O cálculo da altitude é mais simples - dependendo de um bit específico, ele pode ser representado como um múltiplo de 25 ou 100 pés.

Velocidade no ar

Dataframe com TC = 19. O interessante aqui é que a velocidade pode ser relativa ao solo (mais precisa, denominada Velocidade do Solo) e Airspeed, medida pelo sensor de ar da aeronave (pode ser menos precisa por causa do vento). Muitos outros campos diferentes também são transmitidos:


( fonte )

Conclusão


Como podemos ver, a tecnologia ADS-B se tornou uma simbiose interessante, quando um padrão é amplamente utilizável não apenas para profissionais, mas também para usuários comuns. Mas, definitivamente, o papel principal disso foi o barateamento da tecnologia dos receptores digitais SDR, que permite receber sinais com uma frequência maior que gigahertz em um dispositivo muito barato.

O próprio padrão, é claro, possui muito mais dados. Os interessados ​​podem visualizar o PDF na página ICAO ou visitar o site mode-s.org já mencionado acima. É improvável que este artigo seja realmente usado pela maioria dos leitores, mas espero que, pelo menos a ideia geral de como ele funcione, seja agora mais claro.

A propósito, o decodificador ADS-B Python já existe, pode ser investigado no github . Os proprietários dos receptores SDR também podem criar e executar o decodificador ADS-B pronto para uso a partir desta página e (repetirei novamente) alguns detalhes que também apresentamos na primeira parte deste artigo.

O código fonte do analisador, descrito acima, é fornecido sob o spoiler. Este é apenas um exemplo de teste que não pretende a qualidade da produção, mas em geral funciona e pode analisar o arquivo WAV , 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 tenha sido útil, obrigado pela leitura.

Source: https://habr.com/ru/post/pt447078/

More articles:


All Articles