Flightradar24 - چگونه کار می کند؟ قسمت 2، پروتکل ADS-B

سلام حبر. احتمالاً همه کسانی که تا به حال اقوام یا دوستان خود را در هواپیما ملاقات کرده اند یا از آنها خارج شده اند از سرویس رایگان Flightradar24 استفاده کرده اند. این یک راه بسیار راحت برای ردیابی موقعیت هواپیما در زمان واقعی است.

В بخش اول اصل عملیاتی چنین سرویس آنلاین توضیح داده شد. اکنون ادامه می دهیم و متوجه می شویم که چه داده هایی از هواپیما به ایستگاه دریافت کننده ارسال و دریافت می شود و خودمان آن را با استفاده از پایتون رمزگشایی می کنیم.

داستان

بدیهی است که داده های هواپیما برای مشاهده کاربران در گوشی های هوشمند خود مخابره نمی شود. این سیستم ADS-B (پخش خودکار وابسته به نظارت-پخش) نامیده می شود و برای انتقال خودکار اطلاعات مربوط به هواپیما به مرکز کنترل استفاده می شود - شناسه، مختصات، جهت، سرعت، ارتفاع و سایر داده ها منتقل می شود. پیش از این، قبل از ظهور چنین سیستم‌هایی، دیسپچر فقط می‌توانست نقطه‌ای را در رادار ببیند. وقتی هواپیماهای زیادی وجود داشت این دیگر کافی نبود.

از نظر فنی، ADS-B شامل یک فرستنده در هواپیما است که به صورت دوره ای بسته های اطلاعاتی را با فرکانس نسبتاً بالای 1090 مگاهرتز ارسال می کند (حالت های دیگری نیز وجود دارد، اما ما چندان علاقه ای به آنها نداریم، زیرا مختصات فقط در اینجا منتقل می شوند). البته در فرودگاه علاوه بر فرستنده، گیرنده ای هم وجود دارد، اما برای ما کاربران، گیرنده خودمان جالب است.

به هر حال، برای مقایسه، اولین سیستم از این دست، Airnav Radarbox، که برای کاربران عادی طراحی شده بود، در سال 2007 ظاهر شد و حدود 900 دلار قیمت داشت؛ اشتراک خدمات شبکه نیز 250 دلار دیگر در سال هزینه داشت.

نظرات اولین صاحبان روسی را می توان در انجمن خواند رادیواسکنر. اکنون که گیرنده های RTL-SDR به طور گسترده در دسترس قرار گرفته اند، دستگاه مشابهی را می توان با قیمت 30 دلار مونتاژ کرد؛ اطلاعات بیشتر در مورد این موضوع در بخش اول. بیایید به خود پروتکل برویم - بیایید ببینیم چگونه کار می کند.

دریافت سیگنال ها

ابتدا باید سیگنال ضبط شود. مدت زمان کل سیگنال تنها 120 میکروثانیه است، بنابراین برای جداسازی راحت اجزای آن، یک گیرنده SDR با فرکانس نمونه برداری حداقل 5 مگاهرتز مطلوب است.

پس از ضبط، یک فایل WAV با نرخ نمونه برداری 5000000 نمونه در ثانیه دریافت می کنیم؛ 30 ثانیه از چنین ضبطی حدود 500 مگابایت وزن دارد. البته گوش دادن به آن با پخش کننده رسانه بی فایده است - فایل حاوی صدا نیست، بلکه یک سیگنال رادیویی مستقیم دیجیتالی شده است - دقیقاً اینگونه است که رادیو تعریف شده نرم افزار کار می کند.

ما فایل را با استفاده از پایتون باز و پردازش می کنیم. کسانی که می خواهند به تنهایی آزمایش کنند، می توانند نمونه ضبط شده را دانلود کنند по ссылке.

بیایید فایل را دانلود کنیم و ببینیم داخل آن چیست.

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 علاقه مند هستیم، زیرا ... این است که شامل مختصات هواپیما است.

ایکائو (24 بیت) - کد بین المللی منحصر به فرد هواپیما. می توانید هواپیما را با کد آن بررسی کنید آنلاین (متاسفانه، نویسنده به روز رسانی پایگاه داده را متوقف کرده است، اما همچنان مرتبط است). به عنوان مثال، برای کد 3c5ee2 اطلاعات زیر را داریم:

ویرایش: در نظرات در مورد مقاله توضیحات کد ایکائو با جزئیات بیشتری ارائه شده است؛ به علاقه مندان توصیه می کنم آن را مطالعه کنند.

داده ها (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. نکته جالب اینجاست که سرعت می تواند دقیق باشد، نسبت به زمین (سرعت زمین) یا هوابرد، که توسط یک سنسور هواپیما (Airspeed) اندازه گیری می شود. بسیاری از زمینه های مختلف نیز منتقل می شوند:

(منبع)

نتیجه

همانطور که می بینید، فناوری ADS-B به یک همزیستی جالب تبدیل شده است، زمانی که یک استاندارد نه تنها برای حرفه ای ها، بلکه برای کاربران عادی نیز مفید است. اما مطمئناً، فناوری ارزان‌تر گیرنده‌های دیجیتال SDR نقش کلیدی را ایفا کرد که به دستگاه اجازه می‌دهد به معنای واقعی کلمه سیگنال‌هایی را با فرکانس‌های بالاتر از گیگاهرتز دریافت کند.

البته در خود استاندارد موارد بسیار بیشتری وجود دارد. علاقه مندان می توانند فایل پی دی اف را در صفحه مشاهده کنند ایکائو یا به یکی از موارد ذکر شده در بالا مراجعه کنید سایت اینترنتی.

بعید است که همه موارد فوق برای بسیاری مفید باشد، اما حداقل ایده کلی نحوه کار، امیدوارم، باقی بماند.

به هر حال، یک رمزگشای آماده در پایتون از قبل وجود دارد، می توانید آن را مطالعه کنید اینجا. و صاحبان گیرنده های 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()

امیدوارم کسی علاقه مند بوده باشد، از توجه شما متشکرم.

منبع: www.habr.com