Сайн уу Хабр. Онгоцонд хамаатан садан, найз нөхөдтэйгээ уулзаж, үдэж өгсөн хүн бүр Flightradar24 үнэгүй үйлчилгээг ашигласан байх. Энэ нь онгоцны байрлалыг бодит цаг хугацаанд хянах маш тохиромжтой арга юм.
В Ийм онлайн үйлчилгээний үйл ажиллагааны зарчмыг тайлбарлав. Одоо бид нисэх онгоцноос хүлээн авах станц руу ямар өгөгдөл илгээж, хүлээн авч байгааг олж мэдээд Python ашиглан өөрсдөө тайлах болно.
түүх
Ухаалаг гар утсан дээрээ нисэх онгоцны мэдээллийг хэрэглэгчдэд дамжуулахгүй нь ойлгомжтой. Уг системийг ADS-B (Automatic dependent surveillance—broadcast) гэж нэрлэдэг бөгөөд агаарын хөлгийн талаарх мэдээллийг удирдлагын төв рүү автоматаар дамжуулахад ашигладаг - түүний танигч, координат, чиглэл, хурд, өндөр болон бусад өгөгдлийг дамжуулдаг. Өмнө нь ийм систем гарч ирэхээс өмнө диспетчер зөвхөн радар дээрх цэгийг харж чаддаг байсан. Хэт олон онгоц байхад энэ нь хангалтгүй болсон.
Техникийн хувьд ADS-B нь 1090 МГц-ийн нэлээд өндөр давтамжтайгаар мэдээллийн багцыг үе үе илгээдэг онгоц дээрх дамжуулагчаас бүрддэг (бусад горимууд байдаг, гэхдээ координатыг зөвхөн энд дамжуулдаг тул бид тэднийг тийм ч их сонирхдоггүй). Мэдээжийн хэрэг, дамжуулагчаас гадна онгоцны буудлын хаа нэгтээ хүлээн авагч байдаг, гэхдээ хэрэглэгчдийн хувьд бидний хувьд өөрсдийн хүлээн авагч нь сонирхолтой байдаг.
Дашрамд хэлэхэд, жирийн хэрэглэгчдэд зориулагдсан анхны ийм систем болох Airnav Radarbox нь 2007 онд гарч ирсэн бөгөөд ойролцоогоор 900 долларын үнэтэй байсан бөгөөд сүлжээний үйлчилгээнд бүртгүүлэх нь жилд 250 долларын үнэтэй байдаг.
Оросын анхны эздийн талаархи тоймыг форум дээрээс уншиж болно . Одоо RTL-SDR хүлээн авагч өргөн боломжтой болсон тул ижил төстэй төхөөрөмжийг 30 доллараар угсарч болно; энэ талаар дэлгэрэнгүйг эндээс харж болно. . Протокол руугаа явцгаая - энэ нь хэрхэн ажилладагийг харцгаая.
Дохио хүлээн авах
Эхлээд дохиог бүртгэх шаардлагатай. Бүх дохио нь ердөө 120 микросекунд үргэлжлэх тул түүний бүрэлдэхүүн хэсгүүдийг эвтэйхэн задлахын тулд дор хаяж 5 МГц давтамжтай SDR хүлээн авагчийг ашиглах нь зүйтэй.
Бичлэг хийсний дараа бид 5000000 дээж/сек-ийн түүвэрлэлтийн хурдтай WAV файлыг хүлээн авдаг; ийм бичлэгийн 30 секунд нь 500 МБ орчим "жинтэй". Медиа тоглуулагчаар үүнийг сонсох нь мэдээж хэрэг ашиггүй - файл нь дуу чимээ агуулаагүй, харин шууд дижитал радио дохиог агуулдаг - Програм хангамжийн тодорхойлогдсон радио яг ийм байдлаар ажилладаг.
Бид 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 (Downlink Format, 5 бит) - мессежийн төрлийг тодорхойлно. Хэд хэдэн төрөл байдаг:
()
Бид зөвхөн DF17 төрлийг сонирхдог, учир нь... Энэ нь онгоцны координатыг агуулдаг.
ICAO (24 бит) - онгоцны олон улсын өвөрмөц код. Та онгоцыг кодоор нь шалгаж болно (Харамсалтай нь зохиогч мэдээллийн баазыг шинэчлэхээ больсон боловч энэ нь хамааралтай хэвээр байна). Жишээлбэл, 3c5ee2 кодын хувьд бидэнд дараах мэдээлэл байна.
Засварлах: in 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 хүлээн авагчийн хямд технологи гол үүрэг гүйцэтгэсэн бөгөөд энэ нь төхөөрөмжид гигагерцээс дээш давтамжтай дохиог шууд утгаараа хүлээн авах боломжийг олгодог.
Стандартад мэдээж илүү олон зүйл бий. Сонирхсон хүмүүс PDF хуудаснаас үзэх боломжтой эсвэл дээр дурдсан нэгэнд зочилно уу .
Дээр дурдсан бүх зүйл олон хүнд хэрэг болох магадлал багатай ч наад зах нь энэ нь хэрхэн ажилладаг талаар ерөнхий санаа хэвээр байна гэж найдаж байна.
Дашрамд хэлэхэд, Python дээр бэлэн декодлогч аль хэдийн байгаа тул та үүнийг судалж болно . Мөн SDR хүлээн авагчийн эзэд бэлэн ADS-B декодерыг угсарч ажиллуулж болно , энэ талаар илүү дэлгэрэнгүй авч үзсэн .
Өгүүлэлд тайлбарласан задлан шинжлэлийн эх кодыг зүсэлтийн доор өгөв. Энэ бол туршилтын жишээ бөгөөд үйлдвэрлэл гэж дүр эсгэдэггүй, гэхдээ дотор нь зарим зүйл ажилладаг бөгөөд дээр дурдсан файлыг задлан шинжлэхэд ашиглаж болно.
Эх код (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()
Хэн нэгэн сонирхож байсан гэж найдаж байна, анхаарал тавьсанд баярлалаа.
Эх сурвалж: www.habr.com