Merhaba Habr. Muhtemelen uçakta akrabaları veya arkadaşlarıyla tanışan veya onları gören herkes ücretsiz Flightradar24 hizmetini kullanmıştır. Bu, uçağın konumunu gerçek zamanlı olarak takip etmenin çok kullanışlı bir yoludur.
В Böyle bir çevrimiçi hizmetin çalışma prensibi açıklandı. Şimdi devam edeceğiz ve uçaktan alıcı istasyona hangi verilerin gönderilip alındığını bulacağız ve Python kullanarak kendimiz çözeceğiz.
Öykü
Açıkçası, uçak verileri kullanıcıların akıllı telefonlarında görmesi için aktarılmıyor. Sistem ADS-B (Otomatik bağımlı gözetim - yayın) olarak adlandırılır ve uçak hakkındaki bilgileri otomatik olarak kontrol merkezine iletmek için kullanılır - tanımlayıcı, koordinatlar, yön, hız, yükseklik ve diğer veriler iletilir. Daha önce, bu tür sistemlerin ortaya çıkmasından önce, sevk görevlisi radarda yalnızca bir noktayı görebiliyordu. Çok fazla uçak varken bu artık yeterli olmuyordu.
Teknik olarak ADS-B, periyodik olarak 1090 MHz gibi oldukça yüksek bir frekansta bilgi paketleri gönderen bir uçak üzerindeki bir vericiden oluşur (başka modlar da vardır, ancak koordinatlar yalnızca burada iletildiği için onlarla o kadar ilgilenmiyoruz). Elbette vericinin yanı sıra havaalanında bir yerde bir alıcı da var ama biz kullanıcılar için kendi alıcımız ilginç.
Bu arada, karşılaştırma yapmak gerekirse, sıradan kullanıcılar için tasarlanan bu tür ilk sistem Airnav Radarbox 2007'de ortaya çıktı ve maliyeti yaklaşık 900 dolardı; ağ hizmetlerine aboneliğin maliyeti ise yılda 250 dolardı.
Bu ilk Rus sahiplerin yorumları forumda okunabilir . Artık RTL-SDR alıcıları yaygın olarak bulunabildiğinden, benzer bir cihaz 30 dolara monte edilebiliyor; bununla ilgili daha fazla bilgi şuradaydı: . Protokolün kendisine geçelim - nasıl çalıştığını görelim.
Sinyal alma
Öncelikle sinyalin kaydedilmesi gerekiyor. Sinyalin tamamı yalnızca 120 mikrosaniyelik bir süreye sahiptir, bu nedenle bileşenlerini rahatça sökmek için örnekleme frekansı en az 5 MHz olan bir SDR alıcısı tercih edilir.
Kayıttan sonra, örnekleme hızı 5000000 örnek/sn olan bir WAV dosyası alıyoruz; böyle bir kaydın 30 saniyesi yaklaşık 500 MB "ağırlığındadır". Bunu bir medya oynatıcıyla dinlemek elbette işe yaramaz - dosya ses içermez, doğrudan sayısallaştırılmış radyo sinyali içerir - Yazılım Tanımlı Radyo tam olarak bu şekilde çalışır.
Dosyayı Python kullanarak açıp işleyeceğiz. Kendi başına deney yapmak isteyenler örnek kaydı indirebilir .
Dosyayı indirelim ve içinde ne olduğuna bakalım.
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()
Sonuç: Arka plan gürültüsüne karşı bariz "darbeler" görüyoruz.
Her "darbe", grafikteki çözünürlüğü artırırsanız yapısı açıkça görülebilen bir sinyaldir.
Gördüğünüz gibi resim yukarıdaki açıklamada verilenlerle oldukça tutarlı. Verileri işlemeye başlayabilirsiniz.
kod çözme
İlk önce biraz akış almanız gerekiyor. Sinyalin kendisi Manchester kodlaması kullanılarak kodlanır:
Yarım baytlardaki seviye farkından gerçek “0” ve “1” elde etmek kolaydır.
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"
Sinyalin yapısı aşağıdaki gibidir:
Alanlara daha detaylı bakalım.
DF (Downlink Formatı, 5 bit) - mesajın türünü belirler. Birkaç türü vardır:
()
Biz sadece DF17 tipiyle ilgileniyoruz çünkü... Uçağın koordinatlarını içeren şey budur.
ICAO (24 bit) - uçağın uluslararası benzersiz kodu. Uçağı kodundan kontrol edebilirsiniz (maalesef yazar veritabanını güncellemeyi durdurdu, ancak bu hala alakalı). Örneğin, 3c5ee2 kodu için aşağıdaki bilgilere sahibiz:
Düzenleme: içinde ICAO kodunun açıklaması daha ayrıntılı olarak verilmiştir, ilgilenenlerin okumasını tavsiye ederim.
VERİ (56 veya 112 bit) - kodunu çözeceğimiz gerçek veriler. Verinin ilk 5 biti alandır Tür kodu, depolanan verilerin alt türünü içerir (DF ile karıştırılmamalıdır). Bu türlerden epeyce var:
()
Birkaç paket örneğine bakalım.
uçak kimliği
İkili biçimde örnek:
00100 011 000101 010111 000111 110111 110001 111000
Veri alanları:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, her C1-C8 karakteri satırdaki indekslere karşılık gelen kodları içerir:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############0123456789######
Dizinin kodunu çözerek uçak kodunu almak kolaydır: 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('#', ''))
Havadaki konum
İsim basitse koordinatlar daha karmaşıktır. Çift ve tek olmak üzere 2 çerçeve şeklinde iletilir. Alan kodu TC = 01011b = 11.
Çift ve tek paketlere örnek:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Koordinatların hesaplanması oldukça karmaşık bir formüle göre gerçekleşir:
()
Ben bir GIS uzmanı değilim, bu yüzden nereden geldiğini bilmiyorum. Kim bilir yorumlara yazın.
Yükseklik daha basit kabul edilir; spesifik bit'e bağlı olarak 25 veya 100 feet'in katları olarak temsil edilebilir.
Havadaki Hız
TC=19'lu paket. Buradaki ilginç şey, hızın yere göre doğru (Yer Hızı) veya havadaki bir uçak sensörü tarafından ölçülen (Hava Hızı) doğru olabilmesidir. Birçok farklı alan da iletilir:
()
Sonuç
Gördüğünüz gibi ADS-B teknolojisi, bir standardın yalnızca profesyoneller için değil sıradan kullanıcılar için de yararlı olduğu ilginç bir simbiyoz haline geldi. Ancak elbette bunda önemli bir rol, cihazın kelimenin tam anlamıyla "bir kuruş karşılığında" gigahertz'in üzerindeki frekanslara sahip sinyalleri almasına olanak tanıyan daha ucuz dijital SDR alıcı teknolojisi tarafından oynandı.
Standardın kendisinde elbette çok daha fazlası var. İlgilenenler sayfadaki PDF'yi görüntüleyebilir veya yukarıda belirtilen yeri ziyaret edin .
Yukarıdakilerin hepsinin birçok kişi için yararlı olması pek olası değildir, ancak en azından nasıl çalıştığına dair genel fikir kalır, umarım.
Bu arada, Python'da hazır bir kod çözücü zaten var, onu inceleyebilirsiniz . Ve SDR alıcılarının sahipleri hazır bir ADS-B kod çözücüyü monte edip çalıştırabilir , bu daha ayrıntılı olarak tartışıldı .
Makalede anlatılan ayrıştırıcının kaynak kodu kesimin altında verilmiştir. Bu, üretim gibi görünmeyen bir test örneğidir ancak içinde bazı şeyler çalışır ve yukarıda kaydedilen dosyayı ayrıştırmak için kullanılabilir.
Kaynak kodu (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()
Umarım birileri ilgilenmiştir, ilginiz için teşekkürler.
Kaynak: habr.com