Hello Habr. Tingali ang tanan nga nakahimamat o nakakita sa mga paryente o higala sa usa ka eroplano naggamit sa libre nga serbisyo sa Flightradar24. Kini usa ka kombenyente nga paagi aron masubay ang posisyon sa eroplano sa tinuud nga oras.
В Ang prinsipyo sa operasyon sa ingon nga serbisyo sa online gihulagway. Magpadayon kami karon ug mahibal-an kung unsang mga datos ang gipadala ug nadawat gikan sa eroplano hangtod sa estasyon sa pagdawat ug gi-decode kini sa among kaugalingon gamit ang Python.
История
Dayag nga, ang data sa ayroplano wala ipadala aron makita sa mga tiggamit sa ilang mga smartphone. Ang sistema gitawag ug ADS-B (Automatic dependent surveillance—broadcast), ug gigamit sa awtomatikong pagpasa sa impormasyon bahin sa ayroplano ngadto sa control center - ang identifier, coordinates, direksyon, speed, altitude ug uban pang datos niini gipasa. Kaniadto, sa wala pa ang pag-abut sa ingon nga mga sistema, ang dispatser makakita lamang sa usa ka punto sa radar. Dili na kini igo kung daghan kaayo ang mga eroplano.
Sa teknikal nga paagi, ang ADS-B naglangkob sa usa ka transmitter sa usa ka ayroplano nga matag karon ug unya nagpadala sa mga pakete sa kasayuran sa medyo taas nga frequency sa 1090 MHz (adunay ubang mga mode, apan dili kami interesado niini, tungod kay ang mga koordinasyon gipasa lamang dinhi). Siyempre, dugang sa transmitter, adunay usa usab ka receiver sa usa ka dapit sa airport, apan alang kanamo, isip mga tiggamit, ang among kaugalingong receiver makapaikag.
Pinaagi sa dalan, alang sa pagtandi, ang una nga ingon nga sistema, ang Airnav Radarbox, nga gidisenyo alang sa ordinaryong mga tiggamit, nagpakita kaniadtong 2007, ug nagkantidad mga $900; ang usa ka suskrisyon sa mga serbisyo sa network nagkantidad ug $250 sa usa ka tuig.
Ang mga pagsusi sa una nga mga tag-iya sa Russia mabasa sa forum . Karon nga ang RTL-SDR receiver kay kaylap na nga magamit, ang susamang device mahimong matigom sa $30; dugang pa bahin niini anaa sa . Mopadayon kita sa protocol mismo - tan-awon naton kung giunsa kini molihok.
Pagdawat ug signal
Una, kinahanglan nga irekord ang signal. Ang tibuuk nga signal adunay gidugayon nga 120 microseconds, aron komportable nga madisassemble ang mga sangkap niini, usa ka tigdawat sa SDR nga adunay frequency sa sampling nga labing menos 5 MHz ang gitinguha.
Pagkahuman sa pagrekord, nakadawat kami usa ka WAV file nga adunay sampling rate nga 5000000 nga mga sample / sec; 30 segundos sa ingon nga pagrekord "gitimbang" mga 500MB. Ang pagpamati niini gamit ang usa ka media player, siyempre, walay kapuslanan - ang file wala'y sulod nga tingog, apan usa ka direkta nga digitized nga signal sa radyo - mao gayud kini kung giunsa ang Software Defined Radio nagtrabaho.
Among ablihan ug iproseso ang file gamit ang Python. Kadtong gusto nga mag-eksperimento sa ilang kaugalingon maka-download sa usa ka pananglitan nga pagrekord .
Atong i-download ang file ug tan-awon kung unsa ang naa sa sulod.
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()
Resulta: nakita namon ang klaro nga "pulso" batok sa kasaba sa background.
Ang matag "pulso" usa ka signal, ang istruktura niini klaro nga makita kung imong dugangan ang resolusyon sa graph.
Sama sa imong nakita, ang litrato nahiuyon sa gihatag sa paghulagway sa ibabaw. Mahimo nimong sugdan ang pagproseso sa datos.
Pag-decode
Una, kinahanglan nimo nga makakuha usa ka gamay nga sapa. Ang signal mismo gi-encode gamit ang Manchester encoding:
Gikan sa kalainan sa lebel sa mga nibbles dali nga makuha ang tinuod nga "0" ug "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"
Ang istruktura sa signal mismo mao ang mosunod:
Atong tan-awon ang mga uma sa mas detalyado.
DF (Downlink Format, 5 bits) - nagtino sa matang sa mensahe. Adunay pipila ka mga matang:
()
Interesado lang kami sa tipo nga DF17, tungod kay... Kini mao ang naglangkob sa mga coordinate sa eroplano.
ICAO (24 bits) - internasyonal nga talagsaon nga code sa eroplano. Mahimo nimong susihon ang eroplano pinaagi sa code niini (kasubo, ang tagsulat mihunong sa pag-update sa database, apan kini may kalabutan gihapon). Pananglitan, alang sa code 3c5ee2 kami adunay mosunod nga impormasyon:
Edit: sa Ang deskripsyon sa ICAO code gihatag sa mas detalyado; Girekomenda nako nga basahon kini sa mga interesado.
data (56 o 112 bits) - ang aktuwal nga datos nga atong i-decode. Ang unang 5 ka tipik sa datos mao ang field Type Code, nga adunay sulod nga subtype sa datos nga gitipigan (dili malibog sa DF). Adunay pipila niini nga mga matang:
()
Atong tan-awon ang pipila ka mga pananglitan sa mga pakete.
Pag-ila sa ayroplano
Pananglitan sa binary nga porma:
00100 011 000101 010111 000111 110111
Mga natad sa datos:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, matag karakter C1-C8 adunay mga code nga katumbas sa mga indeks sa linya:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######
Pinaagi sa pag-decode sa pisi, dali nga makuha ang code sa eroplano: 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('#', ''))
Posisyon sa hangin
Kung ang ngalan yano, nan ang mga koordinasyon mas komplikado. Gipasa sila sa porma sa 2, parehas ug katingad-an nga mga bayanan. Field code TC = 01011b = 11.
Pananglitan sa bisan ug katingad-an nga mga pakete:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Ang pagkalkula sa mga koordinasyon mismo mahitabo sumala sa usa ka medyo malisud nga pormula:
()
Dili ko eksperto sa GIS, mao nga wala ko kahibalo kung diin kini gikan. Kinsa ang nahibal-an, isulat sa mga komento.
Ang gitas-on gikonsiderar nga mas simple - depende sa espesipikong bit, kini mahimong irepresentar nga usa ka multiple nga 25 o 100 ka mga tiil.
Airborne Velocity
Pakete nga adunay TC=19. Ang makaiikag nga butang dinhi mao nga ang katulin mahimong tukma, relatibo sa yuta (Ground Speed), o airborne, gisukod sa sensor sa ayroplano (Airspeed). Daghang lainlaing natad ang gipasa usab:
()
konklusyon
Sama sa imong makita, ang teknolohiya sa ADS-B nahimong usa ka makapaikag nga symbiosis, kung ang usa ka sumbanan mapuslanon dili lamang sa mga propesyonal, kondili usab sa mga ordinaryong tiggamit. Apan siyempre, usa ka mahinungdanong papel niini ang gidula sa mas barato nga teknolohiya sa digital SDR receivers, nga nagtugot sa device nga literal nga makadawat og mga signal nga adunay mga frequency nga labaw sa gigahertz "alang sa mga pennies."
Sa estandard mismo, siyempre, adunay daghan pa. Kadtong mga interesado makatan-aw sa PDF sa panid o bisitaha ang nahisgotan na sa ibabaw .
Dili tingali nga ang tanan sa ibabaw mahimong mapuslanon sa kadaghanan, apan labing menos ang kinatibuk-ang ideya kung giunsa kini molihok, gilauman ko, nagpabilin.
Pinaagi sa dalan, adunay usa ka andam nga decoder sa Python, mahimo nimo kini tun-an . Ug ang mga tag-iya sa mga tigdawat sa SDR mahimong magtigum ug maglansad sa usa ka andam nga ADS-B decoder , kini gihisgutan sa mas detalyado sa .
Ang source code sa parser nga gihulagway sa artikulo gihatag ubos sa cut. Kini usa ka panig-ingnan sa pagsulay nga wala magpakaaron-ingnon nga produksiyon, apan adunay pipila ka mga butang nga nagtrabaho niini, ug kini magamit aron ma-parse ang file nga natala sa ibabaw.
Source code (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()
Nanghinaut ko nga adunay interesado, salamat sa imong pagtagad.
Source: www.habr.com