Hello Habr. Boya gbogbo eniyan ti o ti pade tabi ri awọn ibatan tabi awọn ọrẹ lori ọkọ ofurufu ti lo iṣẹ Flightradar24 ọfẹ. Eyi jẹ ọna ti o rọrun pupọ lati tọpa ipo ọkọ ofurufu ni akoko gidi.
В Ilana iṣẹ ti iru iṣẹ ori ayelujara ni a ṣe apejuwe. A yoo lọ siwaju ati ṣawari kini data ti a firanṣẹ ati gba lati inu ọkọ ofurufu si ibudo gbigba ati iyipada ti ara wa nipa lilo Python.
История
O han ni, data ọkọ ofurufu ko gbejade fun awọn olumulo lati rii lori awọn fonutologbolori wọn. Eto naa ni a pe ni ADS-B (Kakiri igbẹkẹle aifọwọyi—igbohunsafẹfẹ), ati pe o lo lati gbe alaye laifọwọyi nipa ọkọ ofurufu si ile-iṣẹ iṣakoso - idamọ rẹ, awọn ipoidojuko, itọsọna, iyara, giga ati awọn data miiran ti wa ni gbigbe. Ni iṣaaju, ṣaaju dide ti iru awọn ọna šiše, dispatcher le nikan ri aaye kan lori radar. Eleyi je ko to gun nigba ti nibẹ wà ju ọpọlọpọ awọn ofurufu.
Ni imọ-ẹrọ, ADS-B ni atagba lori ọkọ ofurufu ti o firanṣẹ awọn apo-iwe ti alaye lorekore ni igbohunsafẹfẹ giga ti 1090 MHz (awọn ipo miiran wa, ṣugbọn a ko nifẹ si wọn, nitori pe awọn ipoidojuko ti gbejade nibi nikan). Nitoribẹẹ, ni afikun si atagba, olugba tun wa ni ibikan ni papa ọkọ ofurufu, ṣugbọn fun wa, bi awọn olumulo, olugba tiwa jẹ igbadun.
Nipa ọna, fun lafiwe, iru eto akọkọ, Airnav Radarbox, ti a ṣe apẹrẹ fun awọn olumulo lasan, farahan ni ọdun 2007, o si jẹ nipa $ 900; ṣiṣe alabapin si awọn iṣẹ nẹtiwọki jẹ $ 250 miiran ni ọdun kan.
Awọn atunyẹwo ti awọn oniwun Russian akọkọ ni a le ka lori apejọ naa . Ni bayi pe awọn olugba RTL-SDR ti wa ni ibigbogbo, iru ẹrọ kan le pejọ fun $30; diẹ sii nipa eyi wa ninu . Jẹ ki a lọ si ilana funrararẹ - jẹ ki a wo bii o ṣe n ṣiṣẹ.
Gbigba awọn ifihan agbara
Ni akọkọ, ifihan agbara nilo lati gba silẹ. Gbogbo ifihan agbara ni iye akoko 120 microseconds nikan, nitorinaa lati ṣajọpọ awọn paati rẹ ni itunu, olugba SDR kan pẹlu igbohunsafẹfẹ iṣapẹẹrẹ ti o kere ju 5 MHz jẹ iwunilori.
Lẹhin igbasilẹ, a gba faili WAV kan pẹlu oṣuwọn iṣapẹẹrẹ ti awọn ayẹwo 5000000 / iṣẹju-aaya; 30 iṣẹju-aaya ti iru gbigbasilẹ “wọn” nipa 500MB. Nfeti si pẹlu ẹrọ orin media, nitorinaa, ko wulo - faili ko ni ohun ninu, ṣugbọn ifihan agbara redio ti o ni digitized taara - eyi ni deede bi Software Defined Redio ṣiṣẹ.
A yoo ṣii ati ṣe ilana faili nipa lilo Python. Awọn ti o fẹ lati ṣe idanwo lori ara wọn le ṣe igbasilẹ igbasilẹ apẹẹrẹ kan .
Jẹ ki a ṣe igbasilẹ faili naa ki o wo kini inu.
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()
Abajade: a rii “awọn iṣọn” ti o han gbangba lodi si ariwo lẹhin.
Ọkọọkan “pulse” jẹ ifihan agbara kan, eto eyiti o han gbangba ti o ba pọ si ipinnu lori iyaya naa.
Bi o ti le ri, aworan naa jẹ ibamu pẹlu ohun ti a fun ni apejuwe loke. O le bẹrẹ sisẹ data naa.
Yiyipada
Ni akọkọ, o nilo lati gba ṣiṣan diẹ. Awọn ifihan agbara funrararẹ ti ni koodu nipa lilo fifi koodu Manchester:
Lati iyatọ ipele ni awọn nibbles o rọrun lati gba gidi "0" ati "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"
Ilana ti ifihan funrararẹ jẹ bi atẹle:
Jẹ ki a wo awọn aaye ni alaye diẹ sii.
DF (Iwe kika isalẹ, 5 bits) - pinnu iru ifiranṣẹ naa. Awọn oriṣi pupọ lo wa:
()
A nifẹ nikan ni iru DF17, nitori ... O jẹ eyi ti o ni awọn ipoidojuko ti ọkọ ofurufu naa.
ICAO (24 die-die) - koodu alailẹgbẹ agbaye ti ọkọ ofurufu naa. O le ṣayẹwo ọkọ ofurufu nipasẹ koodu rẹ (Laanu, onkọwe ti dẹkun mimu imudojuiwọn data, ṣugbọn o tun wulo). Fun apẹẹrẹ, fun koodu 3c5ee2 a ni alaye wọnyi:
Ṣatunkọ: in Apejuwe koodu ICAO ni alaye diẹ sii; Mo ṣeduro pe awọn ti o nifẹ si ka.
DATA (56 tabi 112 die-die) - data gangan ti a yoo pinnu. Awọn die-die 5 akọkọ ti data jẹ aaye naa Koodu Iru, ti o ni awọn subtype ti data ti wa ni ipamọ (kii ṣe idamu pẹlu DF). Diẹ ninu awọn iru wọnyi wa:
()
Jẹ ki a wo awọn apẹẹrẹ diẹ ti awọn idii.
Idanimọ ọkọ ofurufu
Apẹẹrẹ ni fọọmu alakomeji:
00100 011 000101 010111 000111 110111 110001 111000
Awọn aaye data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, ohun kikọ kọọkan C1-C8 ni awọn koodu ti o baamu si awọn atọka ninu laini:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#######################0123456789#####
Nipa yiyipada okun, o rọrun lati gba koodu ọkọ ofurufu: 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('#', ''))
Ipo ti afẹfẹ
Ti orukọ ba rọrun, lẹhinna awọn ipoidojuko jẹ idiju diẹ sii. Wọn ti wa ni zqwq ni awọn fọọmu ti 2, ani ati odd awọn fireemu. Koodu aaye TC = 01011b = 11.
Apẹẹrẹ ti awọn apo-iwe paapaa ati aibikita:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Iṣiro ti awọn ipoidojuko funrararẹ waye ni ibamu si agbekalẹ ẹtan kuku:
()
Emi kii ṣe amoye GIS, nitorina Emi ko mọ ibiti o ti wa. Tani o mọ, kọ ninu awọn asọye.
Giga ni a ka pe o rọrun - da lori bit kan pato, o le jẹ aṣoju bi boya ọpọ ti 25 tabi 100 ẹsẹ.
Iyara ti afẹfẹ
Package pẹlu TC=19. Ohun ti o nifẹ si nibi ni pe iyara le jẹ deede, ni ibatan si ilẹ (Iyara Ilẹ), tabi ti afẹfẹ, ti iwọn nipasẹ sensọ ọkọ ofurufu (Airspeed). Ọpọlọpọ awọn aaye oriṣiriṣi tun wa ni gbigbe:
()
ipari
Bii o ti le rii, imọ-ẹrọ ADS-B ti di symbiosis ti o nifẹ, nigbati boṣewa kan wulo kii ṣe si awọn alamọja nikan, ṣugbọn si awọn olumulo lasan. Ṣugbọn nitorinaa, ipa pataki ninu eyi ni a ṣe nipasẹ imọ-ẹrọ ti o din owo ti awọn olugba SDR oni nọmba, eyiti o fun laaye ẹrọ naa lati gba awọn ifihan agbara gangan pẹlu awọn loorekoore loke gigahertz “fun awọn pennies.”
Ninu boṣewa funrararẹ, nitorinaa, pupọ wa. Awọn ti o nifẹ le wo PDF lori oju-iwe naa tabi ṣabẹwo si eyi ti a ti sọ tẹlẹ loke .
Ko ṣeeṣe pe gbogbo awọn ti o wa loke yoo wulo fun ọpọlọpọ, ṣugbọn o kere ju imọran gbogbogbo ti bii o ṣe n ṣiṣẹ, Mo nireti, wa.
Nipa ọna, decoder ti o ti ṣetan ni Python ti wa tẹlẹ, o le kawe rẹ . Ati awọn oniwun ti awọn olugba SDR le pejọ ati ṣe ifilọlẹ decoder ADS-B ti o ti ṣetan , yi ti a sísọ ni diẹ apejuwe awọn ni .
Awọn koodu orisun ti parser ti a ṣalaye ninu nkan naa ni a fun ni isalẹ gige. Eyi jẹ apẹẹrẹ idanwo ti ko ṣe dibọn pe o jẹ iṣelọpọ, ṣugbọn diẹ ninu awọn nkan ṣiṣẹ ninu rẹ, ati pe o le ṣee lo lati ṣe itupalẹ faili ti o gbasilẹ loke.
Koodu orisun (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()
Mo nireti pe ẹnikan nifẹ, o ṣeun fun akiyesi rẹ.
orisun: www.habr.com