Xin chào Habr. Chắc hẳn ai từng gặp hoặc tiễn người thân, bạn bè trên máy bay đều đã sử dụng dịch vụ Flightradar24 miễn phí. Đây là một cách rất thuận tiện để theo dõi vị trí của máy bay trong thời gian thực.
В Nguyên tắc hoạt động của một dịch vụ trực tuyến như vậy đã được mô tả. Bây giờ chúng ta sẽ tiếp tục tìm hiểu dữ liệu nào đang được gửi và nhận từ máy bay đến trạm nhận và tự giải mã dữ liệu đó bằng Python.
Câu chuyện
Rõ ràng, dữ liệu máy bay không được truyền đi để người dùng xem trên điện thoại thông minh của họ. Hệ thống này được gọi là ADS-B (Giám sát phụ thuộc tự động—phát sóng) và được sử dụng để tự động truyền thông tin về máy bay đến trung tâm điều khiển - mã nhận dạng, tọa độ, hướng, tốc độ, độ cao và các dữ liệu khác của nó được truyền đi. Trước đây, trước khi có những hệ thống như vậy, người điều phối chỉ có thể nhìn thấy một điểm trên radar. Điều này không còn đủ khi có quá nhiều máy bay.
Về mặt kỹ thuật, ADS-B bao gồm một máy phát trên máy bay định kỳ gửi các gói thông tin ở tần số khá cao 1090 MHz (có các chế độ khác, nhưng chúng tôi không quan tâm lắm đến chúng vì tọa độ chỉ được truyền ở đây). Tất nhiên, ngoài máy phát, còn có một máy thu ở đâu đó tại sân bay, nhưng đối với chúng tôi, với tư cách là người dùng, máy thu của chính chúng tôi thật thú vị.
Nhân tiện, để so sánh, hệ thống đầu tiên như vậy, Airnav Radarbox, được thiết kế cho người dùng thông thường, xuất hiện vào năm 2007 và có giá khoảng 900 USD; thuê bao dịch vụ mạng tốn thêm 250 USD một năm.
Đánh giá của những chủ sở hữu Nga đầu tiên có thể được đọc trên diễn đàn . Giờ đây, máy thu RTL-SDR đã trở nên phổ biến rộng rãi, một thiết bị tương tự có thể được lắp ráp với giá 30 USD; thông tin thêm về điều này có trong . Hãy chuyển sang giao thức - hãy xem nó hoạt động như thế nào.
Nhận tín hiệu
Đầu tiên, tín hiệu cần được ghi lại. Toàn bộ tín hiệu có thời lượng chỉ 120 micro giây, do đó, để tháo rời các thành phần của nó một cách thoải mái, cần có một máy thu SDR có tần số lấy mẫu ít nhất là 5 MHz.
Sau khi ghi, chúng tôi nhận được một tệp WAV với tốc độ lấy mẫu là 5000000 mẫu/giây, 30 giây ghi như vậy “nặng” khoảng 500 MB. Tất nhiên, nghe nó bằng trình phát đa phương tiện là vô ích - tệp không chứa âm thanh mà là tín hiệu vô tuyến được số hóa trực tiếp - đây chính xác là cách hoạt động của Radio Xác định Phần mềm.
Chúng tôi sẽ mở và xử lý tệp bằng Python. Những ai muốn tự mình thử nghiệm có thể tải xuống bản ghi ví dụ .
Hãy tải tập tin xuống và xem những gì bên trong.
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()
Kết quả: chúng ta thấy các “xung” rõ ràng so với tiếng ồn xung quanh.
Mỗi “xung” là một tín hiệu, cấu trúc của tín hiệu này có thể thấy rõ nếu bạn tăng độ phân giải trên biểu đồ.
Như bạn có thể thấy, hình ảnh khá phù hợp với những gì được đưa ra trong mô tả ở trên. Bạn có thể bắt đầu xử lý dữ liệu.
Giải mã
Đầu tiên, bạn cần có được một luồng bit. Bản thân tín hiệu được mã hóa bằng mã hóa Manchester:
Từ sự chênh lệch cấp độ trong các lần nhấm nháp, bạn có thể dễ dàng nhận được “0” và “1” thực sự.
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"
Cấu trúc của tín hiệu như sau:
Chúng ta hãy nhìn vào các lĩnh vực chi tiết hơn.
DF (Định dạng đường xuống, 5 bit) - xác định loại tin nhắn. Có một số loại:
()
Chúng tôi chỉ quan tâm đến loại DF17, bởi vì... Đây là nơi chứa tọa độ của máy bay.
ICAO (24 bit) - mã duy nhất quốc tế của máy bay. Bạn có thể kiểm tra máy bay bằng mã của nó (rất tiếc là tác giả đã ngừng cập nhật cơ sở dữ liệu nhưng vẫn còn phù hợp). Ví dụ: đối với mã 3c5ee2 chúng ta có thông tin sau:
Chỉnh sửa: trong Mô tả mã ICAO được đưa ra chi tiết hơn; tôi khuyên những ai quan tâm nên đọc nó.
DỮ LIỆU (56 hoặc 112 bit) - dữ liệu thực tế mà chúng tôi sẽ giải mã. 5 bit dữ liệu đầu tiên là trường Mã loại, chứa kiểu con của dữ liệu đang được lưu trữ (đừng nhầm với DF). Có khá nhiều loại như thế này:
()
Hãy xem xét một vài ví dụ về các gói.
Nhận dạng máy bay
Ví dụ ở dạng nhị phân:
00100 011 000101 010111 000111 110111 110001 111000
Các trường dữ liệu:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, mỗi ký tự C1-C8 chứa mã tương ứng với các chỉ số trong dòng:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_################0123456789######
Giải mã chuỗi sẽ dễ dàng lấy được mã máy bay: 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('#', ''))
Vị trí trên không
Nếu tên đơn giản thì tọa độ sẽ phức tạp hơn. Chúng được truyền dưới dạng 2 khung chẵn và lẻ. Mã trường TC = 01011b = 11.
Ví dụ về gói chẵn và gói lẻ:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Việc tính toán tọa độ diễn ra theo một công thức khá phức tạp:
()
Tôi không phải là chuyên gia về GIS nên tôi không biết nó đến từ đâu. Ai biết thì viết bình luận nhé.
Chiều cao được coi là đơn giản hơn - tùy thuộc vào bit cụ thể, nó có thể được biểu thị dưới dạng bội số của 25 hoặc 100 feet.
Vận tốc trên không
Gói hàng có TC=19. Điều thú vị ở đây là tốc độ có thể chính xác, so với mặt đất (Tốc độ mặt đất) hoặc trên không, được đo bằng cảm biến máy bay (Airspeed). Nhiều lĩnh vực khác nhau cũng được truyền tải:
()
Kết luận
Như bạn có thể thấy, công nghệ ADS-B đã trở thành một sự cộng sinh thú vị khi một tiêu chuẩn không chỉ hữu ích cho các chuyên gia mà còn cho cả người dùng thông thường. Nhưng tất nhiên, một vai trò quan trọng trong việc này là do công nghệ rẻ hơn của máy thu SDR kỹ thuật số, cho phép thiết bị nhận tín hiệu theo đúng nghĩa đen với tần số trên gigahertz “chỉ với một xu”.
Tất nhiên, trong bản thân tiêu chuẩn còn có nhiều hơn thế nữa. Những người quan tâm có thể xem bản PDF trên trang hoặc truy cập trang đã được đề cập ở trên .
Không chắc rằng tất cả những điều trên sẽ hữu ích với nhiều người, nhưng ít nhất ý tưởng chung về cách thức hoạt động của nó, tôi hy vọng, vẫn còn.
Nhân tiện, đã có sẵn bộ giải mã bằng Python, bạn có thể nghiên cứu nó . Và chủ sở hữu máy thu SDR có thể lắp ráp và khởi chạy bộ giải mã ADS-B làm sẵn , điều này đã được thảo luận chi tiết hơn trong .
Mã nguồn của trình phân tích cú pháp được mô tả trong bài viết được đưa ra bên dưới phần cắt. Đây là một ví dụ thử nghiệm không giả vờ là sản phẩm nhưng có một số tính năng hoạt động trong đó và nó có thể được sử dụng để phân tích cú pháp tệp được ghi ở trên.
Mã nguồn (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()
Tôi hy vọng ai đó quan tâm, cảm ơn sự quan tâm của bạn.
Nguồn: www.habr.com