Hello Habr. Malaha qof kasta oo waligiis la kulmay ama ku arkay qaraabo ama saaxiibo diyaarad ayaa isticmaalay adeegga Flightradar24 ee bilaashka ah. Tani waa hab aad ugu habboon oo lagula socon karo booska diyaaradda waqtiga dhabta ah.
В Mabda'a hawlgalka ee adeegga internetka ee noocan oo kale ah ayaa lagu tilmaamay. Hadda waxaan sii socon doonaa oo aan ogaan doonaa xogta loo dirayo iyo laga helayo diyaaradda oo la geynayo saldhigga iyo anaga lafteena ayaa go'aamin doona annaga oo isticmaalaya Python.
История
Sida cad, xogta diyaarada looma gudbiyo isticmaalayaasha si ay ugu arkaan taleefanadooda casriga ah. Nidaamka waxaa loo yaqaan ADS-B (ilaalinta tooska ah ee ku-tiirsanaanta-baahinta), waxaana loo isticmaalaa in si toos ah loogu gudbiyo macluumaadka ku saabsan diyaaradda xarunta kontoroolka - aqoonsigeeda, isku-duwayaasha, jihada, xawaaraha, joogitaanka iyo xogta kale ayaa la kala qaadaa. Markii hore, ka hor inta aan la helin nidaamyadan, soo diruhu wuxuu arki karaa oo kaliya dhibic radar ah. Tani kuma filna markii ay jireen diyaarado aad u badan.
Farsamo ahaan, ADS-B waxay ka kooban tahay gudbiye diyaarad kaas oo si xilliyo ah u soo dira xirmooyin macluumaad ah inta jeer ee aad u sareysa ee 1090 MHz (waxaa jira habab kale, laakiin aad uma xiisaynno iyaga, maaddaama isku-duwayaasha lagu kala qaado halkan oo keliya). Dabcan, marka lagu daro gudbiyaha, waxaa sidoo kale jira aqbale meel ka mid ah garoonka diyaaradaha, laakiin annaga, isticmaalayaasha, qaataha noo gaar ah ayaa xiiso leh.
Jid ahaan, marka la barbardhigo, nidaamkii ugu horreeyay ee noocaan ah, Airnav Radarbox, oo loogu talagalay isticmaaleyaasha caadiga ah, ayaa soo muuqday 2007, waxayna ku kacday qiyaastii $ 900; ku-qoritaanka adeegyada shabakadda waxay ku kacaysaa $ 250 kale sanadkii.
Dib u eegista milkiilayaasha Ruushka ee ugu horreeya ayaa laga akhriyi karaa madasha . Hadda oo qaataha RTL-SDR ay noqdeen kuwo si weyn loo heli karo, qalab la mid ah ayaa lagu soo ururin karaa $30; wax badan oo ku saabsan tani waxay ahayd . Aan u gudubno borotokoolka laftiisa - aan aragno sida uu u shaqeeyo.
Helitaanka calaamadaha
Marka hore, signalka wuxuu u baahan yahay in la duubo. Calaamadaha oo dhami waxay soconayaan 120 microsecond oo keliya, si si raaxo leh loo kala saaro qaybaheeda, qaataha SDR oo leh inta jeer ee muunad ee ugu yaraan 5 MHz waa la jecel yahay.
Duubista ka dib, waxaan helnaa faylka WAV oo leh qiyaas muunad ah 5000000 muunado/sek; 30 ilbiriqsi ee duubista noocaas ah "miisaam" qiyaastii 500MB. Ku dhegeysiga qalabka warbaahinta, dabcan, waa wax aan faa'iido lahayn - feylku kuma jiro dhawaaq, laakiin calaamad raadiyaha tooska ah ee digitized - tani waa sida saxda ah ee Raadiyaha lagu qeexay Software.
Waxaan furi doonaa oo aan ka baaraandegi doonaa faylka anagoo adeegsanayna Python. Kuwa raba inay iskood tijaabiyaan waxay soo dejisan karaan duubista tusaalaha .
Aan soo dejino faylka oo aragno waxa ku jira.
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()
Natiijadu: waxaan aragnaa "pulses" cad oo ka soo horjeeda qaylada dambe.
Mid kasta oo "pulse" waa calaamad, qaab-dhismeedka kaas oo si cad u muuqda haddii aad kordhiso xallinta garaafka.
Sida aad arki karto, sawirku aad buu ula socdaa waxa lagu sheegay sharraxa sare. Waxaad bilaabi kartaa habaynta xogta
Dejinta
Marka hore, waxaad u baahan tahay inaad hesho qulqul yar. Calaamadda lafteeda ayaa lagu dhejiyay iyadoo la adeegsanayo codaynta Manchester:
Marka loo eego farqiga heerka ee nibbles way fududahay in la helo "0" iyo "1" dhab ah.
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"
Qaab dhismeedka ishaarada lafteedu waa sida soo socota:
Aynu si faahfaahsan u eegno beeraha.
DF (Qaabka Downlink, 5 bits) - wuxuu go'aamiyaa nooca fariinta. Waxaa jira dhowr nooc:
()
Waxaan daneyneynaa kaliya nooca DF17, sababtoo ah ... Waa tan ka kooban isku-duwayaasha diyaaradda.
ICAO (24-bits) - koodka caalamiga ah ee diyaaradda. Waxaad ku hubin kartaa diyaaradda koodkeeda (nasiib darro, qoraagu wuu joojiyay cusboonaysiinta xogta xogta, laakiin wali waa ku habboon tahay). Tusaale ahaan, code 3c5ee2 waxaanu haynaa macluumaadka soo socda:
Wax ka beddel: gudaha Sharaxaadda koodhka ICAO si faahfaahsan ayaa loo bixiyaa, waxaan ku talinayaa kuwa danaynaya inay akhriyaan.
MACLUUMAADKA (56 ama 112 bits) - xogta dhabta ah ee aan go'aamin doono. 5-bit ee ugu horreeya ee xogta waa goobta Nambarka Nooca ah, oo ka kooban nooca-hoosaadka xogta la kaydinayo (aan lagu khaldin DF). Waxaa jira dhowr nooc oo kuwan ah:
()
Aan eegno dhowr tusaale oo xirmo ah.
Aqoonsiga diyaaradda
Tusaale ahaan qaabka laba-geesoodka ah:
00100 011 000101 010111 000111 110111 110001 111000
Goobaha xogta:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, xaraf kasta C1-C8 waxa uu ka kooban yahay koodh u dhigma tusmooyinka xariiqa:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ###############################################0123456789
Marka la go'aamiyo xadhigga, way fududahay in la helo lambarka diyaaradda: 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('#', ''))
Booska hawada
Haddii magacu fudud yahay, markaa isku-duwayaasha ayaa aad u dhib badan. Waxa lagu kala qaadaa qaab 2 ah, xattaa oo aan caadi ahayn. Koodhka goobta TC = 01011b = 11.
Tusaalaha xirmooyinka xitaa iyo kuwa aan fiicneyn:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Xisaabinta isku-duwayaasha lafteedu waxay u dhacdaa iyadoo loo eegayo qaacido aad u dhib badan:
()
Anigu ma ihi khabiir GIS, marka ma garanayo meesha ay ka timid. Yaa garanaya, ku qor faallooyinka.
Dhererka waxaa loo arkaa mid fudud - iyadoo ku xiran qaniinyada gaarka ah, waxaa loo soo bandhigi karaa mid ka mid ah dhowr 25 ama 100 cagood.
Xawaaraha Hawada
Xidhmada leh TC=19. Arrinta xiisaha leh ee halkan ku jirta ayaa ah in xawaaruhu uu noqon karo mid sax ah, marka loo eego dhulka (Xawaaraha dhulka), ama hawada, oo lagu cabbiro dareeraha diyaaradda (Airspeed). Qaybo badan oo kala duwan ayaa sidoo kale la kala qaadaa:
()
gunaanad
Sida aad arki karto, tignoolajiyada ADS-B waxay noqotay symbiosis xiiso leh, marka heerku faa'iido u leeyahay maaha oo kaliya xirfadlayaasha, laakiin sidoo kale isticmaalayaasha caadiga ah. Laakiin dabcan, doorka muhiimka ah ee tan waxaa ciyaaray tikniyoolajiyadda jaban ee qaataha SDR dhijitaalka ah, taas oo u oggolaanaysa aaladda in ay si dhab ah u hesho calaamado leh xad dhaaf ka sarreeya gigahertz “sannad lacageedka.”
Halbeegga laftiisa, dabcan, wax badan ayaa jira. Kuwa danaynaya waxay arki karaan PDF-ka bogga ama booqo kan aan kor ku soo sheegnay .
Uma badna in dhammaan kuwa kor ku xusan ay faa'iido u yeelan doonaan qaar badan, laakiin ugu yaraan fikradda guud ee sida ay u shaqeyso, waxaan rajeynayaa, inay sii jirto.
Jid ahaan, decoder diyaar ah oo Python ku jira ayaa horay u jiray, waad baran kartaa . Iyo mulkiilayaasha qaata SDR-ku way soo ururin karaan oo bilaabi karaan qalabeeyaha ADS-B diyaarsan , arrintan ayaa si faahfaahsan looga hadlay .
Koodhka isha ee falanqaynta lagu sharraxay maqaalka ayaa lagu bixiyaa hoos goynta. Tani waa tusaale tijaabo ah oo aan iska dhigin wax soo saar, laakiin waxyaabaha qaar ayaa ka dhex shaqeeya, waxaana loo isticmaali karaa in lagu kala saaro faylka kor lagu duubay.
Koodhka isha (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()
Waxaan rajeynayaa in qof uu xiiseynayay, waad ku mahadsan tahay dareenkaaga.
Source: www.habr.com