Flightradar24 - 어떻게 작동하나요? 2부, ADS-B 프로토콜

안녕하세요 하브르입니다. 아마도 비행기에서 친척이나 친구를 만나거나 배웅해 본 적이 있는 모든 사람은 무료 Flightradar24 서비스를 사용해 본 적이 있을 것입니다. 이는 실시간으로 항공기의 위치를 ​​추적하는 매우 편리한 방법입니다.

В 첫 번째 부분 이러한 온라인 서비스의 운영 원리를 설명했습니다. 이제 항공기에서 수신국으로 어떤 데이터가 송수신되는지 파악하고 Python을 사용하여 직접 디코딩하겠습니다.

이야기

당연히 항공기 데이터는 사용자가 스마트폰에서 볼 수 있도록 전송되지 않습니다. 이 시스템은 ADS-B(자동 종속 감시 - 방송)라고 하며 항공기에 대한 정보를 관제 센터에 자동으로 전송하는 데 사용됩니다. 즉, 식별자, 좌표, 방향, 속도, 고도 및 기타 데이터가 전송됩니다. 이전에는 이러한 시스템이 출현하기 전에는 파견자가 레이더의 한 지점만 볼 수 있었습니다. 비행기가 너무 많으면 더 이상 충분하지 않았습니다.

기술적으로 ADS-B는 1090MHz의 상당히 높은 주파수에서 정보 패킷을 주기적으로 전송하는 항공기의 송신기로 구성됩니다(다른 모드도 있지만 좌표는 여기에서만 전송되기 때문에 그다지 관심이 없습니다). 물론 송신기 ​​외에 공항 어딘가에 수신기도 있지만 사용자인 우리에게는 자체 수신기가 흥미롭습니다.

그런데 비교를 위해 일반 사용자를 위해 설계된 최초의 시스템인 Airnav Radarbox가 2007년에 출시되었으며 가격은 약 900달러이며, 네트워크 서비스 가입 비용은 연간 250달러입니다.

최초의 러시아 소유자에 대한 리뷰는 포럼에서 읽을 수 있습니다. 라디오 스캐너. 이제 RTL-SDR 수신기가 널리 보급되었으므로 유사한 장치를 30달러에 조립할 수 있습니다. 이에 대한 자세한 내용은 첫 번째 부분. 프로토콜 자체로 넘어가서 어떻게 작동하는지 살펴보겠습니다.

신호 수신

먼저 신호를 녹음해야 합니다. 전체 신호의 지속 시간은 120마이크로초에 불과하므로 해당 구성 요소를 편안하게 분해하려면 샘플링 주파수가 5MHz 이상인 SDR 수신기가 바람직합니다.

녹음 후에 우리는 샘플링 속도가 5000000 샘플/초인 WAV 파일을 수신합니다. 이러한 녹음의 30초 분량은 약 500MB입니다. 물론 미디어 플레이어로 듣는 것은 쓸모가 없습니다. 파일에는 사운드가 포함되어 있지 않지만 직접 디지털화된 라디오 신호가 포함되어 있습니다. 이것이 바로 소프트웨어 정의 라디오가 작동하는 방식입니다.

Python을 사용하여 파일을 열고 처리하겠습니다. 스스로 실험하고 싶은 분은 예제 녹음을 다운로드할 수 있습니다. 링크.

파일을 다운로드해서 내용을 살펴보겠습니다.

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()

결과: 배경 소음에 비해 분명한 "펄스"가 보입니다.

각 "펄스"는 신호이며, 그래프의 해상도를 높이면 구조가 명확하게 표시됩니다.

보시다시피, 그림은 위의 설명에 제공된 내용과 상당히 일치합니다. 데이터 처리를 시작할 수 있습니다.

디코딩

먼저 비트 스트림을 얻어야 합니다. 신호 자체는 맨체스터 인코딩을 사용하여 인코딩됩니다.

니블의 레벨 차이를 통해 실제 "0"과 "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 += "0"
        elif avg1 > avg2:
            bits_str += "1"

신호 자체의 구조는 다음과 같습니다.

필드를 좀 더 자세히 살펴보겠습니다.

DF (다운링크 형식, 5비트) - 메시지 유형을 결정합니다. 여러 가지 유형이 있습니다:

(테이블 소스)

우리는 DF17 유형에만 관심이 있습니다. 왜냐하면... 이것은 항공기의 좌표를 포함하는 것입니다.

ICAO (24비트) - 항공기의 국제 고유 코드입니다. 코드로 비행기를 확인할 수 있습니다 온라인으로 (불행히도 저자는 데이터베이스 업데이트를 중단했지만 여전히 관련이 있습니다). 예를 들어 코드 3c5ee2의 경우 다음 정보가 있습니다.

편집 : 에서 기사에 대한 의견 ICAO 코드에 대한 설명이 더 자세히 나와 있으니 관심 있는 분들은 꼭 읽어보시길 권합니다.

데이터 (56 또는 112비트) - 디코딩할 실제 데이터입니다. 데이터의 처음 5비트는 필드입니다. 유형 코드, 저장되는 데이터의 하위 유형을 포함합니다(DF와 혼동하지 마세요). 다음과 같은 유형이 많이 있습니다.

(테이블 소스)

패키지의 몇 가지 예를 살펴보겠습니다.

항공기 식별

바이너리 형식의 예:

00100 011 000101 010111 000111 110111 110001 111000

데이터 필드:

+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+

TC = 00100b = 4, 각 문자 C1-C8에는 줄의 인덱스에 해당하는 코드가 포함됩니다.
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######

문자열을 해독하면 항공기 코드 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('#', ''))

공수 위치

이름이 단순하면 좌표가 더 복잡해집니다. 2개의 짝수 및 홀수 프레임 형태로 전송됩니다. 필드 코드 TC = 01011b = 11.

짝수 및 홀수 패킷의 예:

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

좌표 자체의 계산은 다소 까다로운 공식에 따라 발생합니다.

(출처)

나는 GIS 전문가가 아니기 때문에 그것이 어디서 나오는지 모르겠습니다. 아시는 분 댓글에 적어주세요.

높이는 더 간단한 것으로 간주됩니다. 특정 비트에 따라 25피트 또는 100피트의 배수로 표시될 수 있습니다.

공수 속도

TC=19인 패키지. 여기서 흥미로운 점은 속도가 지면을 기준으로 정확할 수도 있고(Ground Speed) 항공기 센서로 측정된 공중 속도(Airspeed)일 수도 있다는 것입니다. 다양한 필드도 전송됩니다.

(출처)

결론

보시다시피 ADS-B 기술은 표준이 전문가뿐만 아니라 일반 사용자에게도 유용할 때 흥미로운 공생 관계가 되었습니다. 그러나 물론 여기서 중요한 역할은 저렴한 디지털 SDR 수신기 기술에 의해 수행되었습니다. 이를 통해 장치는 문자 그대로 "XNUMX센트에" 기가헤르츠 이상의 주파수로 신호를 수신할 수 있습니다.

물론 표준 자체에는 훨씬 더 많은 것이 있습니다. 관심 있는 사람은 해당 페이지에서 PDF를 볼 수 있습니다. ICAO 또는 위에서 이미 언급한 곳을 방문하세요. 웹 사이트.

위의 모든 내용이 많은 사람들에게 유용할 것 같지는 않지만 적어도 작동 방식에 대한 일반적인 아이디어는 남아 있기를 바랍니다.

그건 그렇고, 기성 Python 디코더가 이미 존재하므로 연구할 수 있습니다. 여기에. SDR 수신기 소유자는 기성품 ADS-B 디코더를 조립하고 출시할 수 있습니다. 페이지에서, 이는 에서 더 자세히 논의되었습니다. 첫 번째 부분.

기사에 설명된 파서의 소스 코드는 컷 아래에 나와 있습니다. 이것은 프로덕션을 가장하지 않은 테스트 예제이지만 그 안에서 몇 가지 작업이 수행되며 위에 기록된 파일을 구문 분석하는 데 사용할 수 있습니다.
소스 코드(파이썬)

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()

누군가가 관심을 가졌기를 바랍니다. 관심을 가져주셔서 감사합니다.

출처 : habr.com