Nyob zoo Habr. Tej zaum txhua tus neeg uas tau ntsib lossis pom tawm ntawm cov txheeb ze lossis phooj ywg ntawm lub dav hlau tau siv qhov kev pabcuam Flightradar24 dawb. Qhov no yog ib txoj hauv kev yooj yim heev los taug qab lub dav hlau txoj haujlwm hauv lub sijhawm.
В Lub hauv paus ntsiab lus ntawm kev ua haujlwm ntawm qhov kev pabcuam online no tau piav qhia. Tam sim no peb yuav mus tom ntej thiab txheeb xyuas seb cov ntaub ntawv twg raug xa thiab tau txais los ntawm lub dav hlau mus rau qhov chaw txais thiab txiav txim siab nws tus kheej siv Python.
История
Obviously, cov ntaub ntawv dav hlau tsis xa mus rau cov neeg siv pom ntawm lawv cov smartphones. Lub kaw lus hu ua ADS-B (Tsis siv neeg nyob soj ntsuam - tshaj tawm), thiab yog siv los xa cov ntaub ntawv hais txog lub dav hlau mus rau qhov chaw tswj hwm - nws tus cim, kev tswj hwm, kev taw qhia, ceev, qhov siab thiab lwm yam ntaub ntawv raug xa mus. Yav dhau los, ua ntej qhov tshwm sim ntawm cov tshuab zoo li no, tus neeg xa khoom tuaj yeem pom qhov taw qhia ntawm radar. Qhov no tsis txaus thaum muaj dav hlau ntau dhau lawm.
Technically, ADS-B muaj lub tshuab xa tawm ntawm lub dav hlau uas xa cov ntaub ntawv ib ntus ntawm qhov zaus ntawm 1090 MHz (muaj lwm hom, tab sis peb tsis txaus siab rau lawv, txij li kev sib koom tes tau kis ntawm no xwb). Ntawm chav kawm, ntxiv rau lub transmitter, kuj muaj ib tug receiver qhov chaw ntawm lub tshav dav hlau, tab sis rau peb, raws li cov neeg siv, peb tus kheej receiver yog nthuav.
Los ntawm txoj kev, rau kev sib piv, thawj qhov system, Airnav Radarbox, tsim rau cov neeg siv zoo tib yam, tau tshwm sim hauv xyoo 2007, thiab raug nqi txog $ 900; kev tso npe rau cov kev pabcuam network raug nqi ntxiv $ 250 ib xyoos.
Kev tshuaj xyuas ntawm thawj tus tswv Lavxias tuaj yeem nyeem hauv lub rooj sab laj . Tam sim no uas RTL-SDR receivers tau dhau los ua dav, cov cuab yeej zoo sib xws tuaj yeem sib sau rau $ 30; ntau ntxiv txog qhov no yog nyob rau hauv . Cia peb txav mus rau qhov kev cai nws tus kheej - cia saib nws ua haujlwm li cas.
Txais cov cim
Ua ntej, lub teeb liab yuav tsum tau sau tseg. Tag nrho cov teeb liab muaj lub sijhawm tsuas yog 120 microseconds, yog li kom yooj yim tshem tawm nws cov khoom, ib tus neeg txais SDR nrog qhov piv txwv ntawm tsawg kawg 5 MHz yog qhov xav tau.
Tom qab kaw, peb tau txais cov ntaub ntawv WAV nrog tus nqi piv txwv ntawm 5000000 qauv / sec; 30 vib nas this ntawm cov ntaub ntawv no "hnyav" txog 500MB. Mloog nws nrog cov neeg siv xov xwm, tau kawg, tsis muaj txiaj ntsig - cov ntaub ntawv tsis muaj suab, tab sis lub xov tooj cua ncaj qha digitized - qhov no yog raws nraim li Software Defined Radio ua haujlwm.
Peb yuav qhib thiab ua cov ntaub ntawv siv Python. Cov neeg uas xav sim lawv tus kheej tuaj yeem rub tawm ib qho piv txwv kaw .
Wb download tau cov ntaub ntawv thiab saib dab tsi nyob rau hauv.
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()
Qhov tshwm sim: peb pom pom tseeb "pulses" tiv thaiv lub suab nrov tom qab.
Txhua "pulse" yog lub teeb liab, cov qauv uas pom meej yog tias koj nce qhov kev daws teeb meem ntawm daim duab.
Raws li koj tuaj yeem pom, daim duab zoo ib yam nrog cov lus piav qhia saum toj no. Koj tuaj yeem pib ua cov ntaub ntawv.
Kev txiav txim siab
Ua ntej, koj yuav tsum tau txais me ntsis kwj. Lub teeb liab nws tus kheej yog encoded siv Manchester encoding:
Los ntawm qib sib txawv hauv nibbles nws yooj yim kom tau txais tiag tiag "0" thiab "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"
Cov qauv ntawm lub teeb liab nws tus kheej yog raws li nram no:
Cia peb saib cov teb kom ntxaws ntxiv.
DF (Downlink Format, 5 bits) - txiav txim siab hom lus. Muaj ntau ntau hom:
()
Peb tsuas yog txaus siab rau hom DF17, vim tias ... Nws yog qhov no uas muaj cov coordinates ntawm lub aircraft.
ICAO (24 khoom) - thoob ntiaj teb cim cim ntawm lub dav hlau. Koj tuaj yeem tshawb xyuas lub dav hlau los ntawm nws cov cai (Hmoov tsis zoo, tus sau tau tso tseg tsis hloov kho cov ntaub ntawv, tab sis nws tseem cuam tshuam). Piv txwv li, rau code 3c5ee2 peb muaj cov ntaub ntawv hauv qab no:
Khob: in Cov lus piav qhia ntawm ICAO code yog muab rau hauv kev nthuav dav ntxiv; Kuv xav kom cov neeg nyiam nyeem nws.
NTAUB NTAWV (56 lossis 112 khoom) - cov ntaub ntawv tiag tiag uas peb yuav txiav txim siab. Thawj 5 cov ntaub ntawv yog cov teb Hom Code, uas muaj cov subtype ntawm cov ntaub ntawv khaws cia (tsis txhob tsis meej pem nrog DF). Muaj ob peb yam ntawm cov no:
()
Cia peb saib ob peb yam piv txwv ntawm pob.
Kev txheeb xyuas lub dav hlau
Piv txwv nyob rau hauv binary daim ntawv:
00100 011 000101 010111 000111 110111 110001 111000
Cov ntaub ntawv teb:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, txhua tus cim C1-C8 muaj cov lej sib raug rau cov lej hauv kab:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##############################################
Los ntawm kev txiav txim siab txoj hlua, nws yog qhov yooj yim kom tau txais lub dav hlau code: 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('#', ''))
Huab cua txoj hauj lwm
Yog hais tias lub npe yog yooj yim, ces lub coordinates yuav nyuaj. Lawv kis tau nyob rau hauv daim ntawv ntawm 2, txawm thiab khib thav ntawv. Field code TC = 01011b = 11.
Piv txwv ntawm cov pob ntawv sib npaug thiab khib:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Kev suav ntawm kev tswj hwm nws tus kheej tshwm sim raws li cov qauv yooj yim:
()
Kuv tsis yog tus kws tshaj lij GIS, yog li kuv tsis paub qhov twg los. Leej twg paub, sau rau hauv cov lus.
Qhov siab yog suav hais tias yooj yim - nyob ntawm qhov tshwj xeeb me ntsis, nws tuaj yeem sawv cev ua ntau yam ntawm 25 lossis 100 ko taw.
Huab cua Velocity
Pob nrog TC = 19. Qhov nthuav ntawm no yog qhov ceev tuaj yeem yog qhov tseeb, txheeb ze rau hauv av (Ground Speed), lossis huab cua, ntsuas los ntawm lub dav hlau sensor (Airspeed). Ntau qhov chaw sib txawv kuj kis tau:
()
xaus
Raws li koj tuaj yeem pom, ADS-B thev naus laus zis tau dhau los ua ib qho kev nthuav dav, thaum tus qauv tseem ceeb tsis yog rau cov kws tshaj lij xwb, tab sis kuj rau cov neeg siv zoo tib yam. Tab sis tau kawg, lub luag haujlwm tseem ceeb hauv qhov no tau ua los ntawm cov thev naus laus zis pheej yig dua ntawm cov neeg txais SDR digital, uas tso cai rau lub cuab yeej kom tau txais cov lus qhia nrog ntau zaus tshaj li gigahertz "rau pennies."
Nyob rau hauv tus qauv nws tus kheej, ntawm chav kawm, muaj ntau ntau. Cov neeg txaus siab tuaj yeem saib PDF ntawm nplooj ntawv los yog mus ntsib tus uas twb hais los saum no lawm .
Nws tsis zoo li tag nrho cov saum toj no yuav muaj txiaj ntsig zoo rau ntau tus, tab sis tsawg kawg yog lub tswv yim dav dav ntawm nws ua haujlwm li cas, kuv vam tias, tseem nyob.
Los ntawm txoj kev, npaj txhij-ua decoder hauv Python twb muaj lawm, koj tuaj yeem kawm nws . Thiab cov tswv ntawm SDR receivers tuaj yeem sib sau ua ke thiab tso tawm ADS-B decoder npaj txhij , qhov no tau tham ntau yam ntxiv hauv .
Lub hauv paus code ntawm tus parser piav nyob rau hauv tsab xov xwm yog muab hauv qab no qhov kev txiav. Qhov no yog ib qho piv txwv sim uas tsis ua txuj ua ntau lawm, tab sis qee yam ua haujlwm hauv nws, thiab nws tuaj yeem siv los txheeb xyuas cov ntaub ntawv sau tseg saum toj no.
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()
Kuv vam tias ib tug neeg xav tau, ua tsaug rau koj cov lus qhia.
Tau qhov twg los: www.hab.com