Halo Habr. Is dòcha gu bheil a h-uile duine a choinnich no a chunnaic a-riamh air càirdean no caraidean air plèana air an t-seirbheis Flightradar24 an-asgaidh a chleachdadh. Is e dòigh air leth goireasach a tha seo airson sùil a chumail air suidheachadh an itealain ann an àm fìor.
В Chaidh iomradh a thoirt air prionnsapal obrachaidh seirbheis air-loidhne mar sin. Thèid sinn air adhart a-nis gus faighinn a-mach dè an dàta a thathas a’ cur agus a’ faighinn bhon itealan chun stèisean faighinn agus ga dhì-chòdachadh sinn fhìn a’ cleachdadh Python.
История
Gu follaiseach, chan eil dàta itealain air a ghluasad airson luchd-cleachdaidh fhaicinn air na fònaichean sgairteil aca. Canar ADS-B ris an t-siostam (sgrùdadh eisimeileach fèin-ghluasadach - craoladh), agus tha e air a chleachdadh gus fiosrachadh mun itealan a chuir chun ionad smachd gu fèin-ghluasadach - tha an aithnichear, co-chomharran, stiùireadh, astar, àirde agus dàta eile air a ghluasad. Roimhe sin, mus tàinig na siostaman sin, cha robh an neach-tagraidh a 'faicinn ach puing air an radar. Cha robh seo gu leòr tuilleadh nuair a bha cus phlèanaichean ann.
Gu teicnigeach, tha ADS-B a ’toirt a-steach inneal-sgaoilidh air itealan a bhios bho àm gu àm a’ cur pacaidean fiosrachaidh aig tricead gu math àrd de 1090 MHz (tha modhan eile ann, ach chan eil uimhir de dh ’ùidh againn annta, leis gu bheil co-chomharran air an gluasad a-mhàin an seo). Gu dearbh, a bharrachd air an inneal-sgaoilidh, tha inneal-glacaidh an àiteigin aig a 'phort-adhair cuideachd, ach dhuinne, mar luchd-cleachdaidh, tha an cuidhteas againn fhèin inntinneach.
Co-dhiù, airson coimeas a dhèanamh, nochd a’ chiad shiostam den leithid, Airnav Radarbox, a chaidh a dhealbhadh airson luchd-cleachdaidh àbhaisteach, ann an 2007, agus chosg e timcheall air $900; cosgaidh ballrachd gu seirbheisean lìonra $250 eile sa bhliadhna.
Faodar lèirmheasan air na ciad luchd-seilbh Ruiseanach sin a leughadh air an fhòram . A-nis gu bheil luchd-glacaidh RTL-SDR rim faighinn gu farsaing, faodar inneal coltach ris a chruinneachadh airson $ 30; bha barrachd mu dheidhinn seo ann an . Gluaisidh sinn air adhart chun phròtacal fhèin - chì sinn mar a tha e ag obair.
A 'faighinn comharran
An toiseach, feumar an comharra a chlàradh. Chan eil fad a’ chomharra gu lèir ach 120 microseconds, agus mar sin gus na co-phàirtean aige a thoirt air falbh gu comhfhurtail, tha inneal-glacaidh SDR le tricead samplachaidh co-dhiù 5 MHz ion-mhiannaichte.
Às deidh dhuinn clàradh, gheibh sinn faidhle WAV le ìre samplachaidh de 5000000 sampall / diog; 30 diogan den leithid de chlàradh “cuideam” timcheall air 500MB. Tha a bhith ag èisteachd ris le cluicheadair meadhanan, gu dearbh, gun fheum - chan eil fuaim anns an fhaidhle, ach comharra rèidio didseatach gu dìreach - seo dìreach mar a tha Bathar-bog Mìnichte Radio ag obair.
Fosglaidh agus giullachd sinn am faidhle a’ cleachdadh Python. Faodaidh an fheadhainn a tha airson feuchainn leotha fhèin clàradh eisimpleir a luchdachadh sìos .
Leig leinn am faidhle a luchdachadh sìos agus faicinn dè a tha a-staigh.
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()
Toradh: chì sinn “pulses” follaiseach an aghaidh fuaim a’ chùl.
Tha gach “pulse” na chomharradh, agus tha an structar ri fhaicinn gu soilleir ma mheudaicheas tu an rùn air a’ ghraf.
Mar a chì thu, tha an dealbh gu math co-chòrdail ris na tha air a thoirt seachad san tuairisgeul gu h-àrd. Faodaidh tu tòiseachadh air an dàta a ghiullachd.
Dì-chòdachadh
An toiseach, feumaidh tu beagan sruth fhaighinn. Tha an comharra fhèin air a chòdachadh a’ cleachdadh còdachadh Manchester:
Bhon eadar-dhealachadh ìre ann an nibbles tha e furasta faighinn fìor “0” agus “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"
Tha structar a 'chomharra fhèin mar a leanas:
Bheir sinn sùil nas mionaidiche air na raointean.
DF (Cruth Downlink, 5 pìosan) - a 'dearbhadh an seòrsa teachdaireachd. Tha grunn sheòrsaichean ann:
()
Chan eil ùidh againn ach ann an seòrsa DF17, oir ... Is e seo a tha a’ toirt a-steach co-chomharran an itealain.
ICAO (24 pìosan) - còd sònraichte eadar-nàiseanta an itealain. Faodaidh tu sgrùdadh a dhèanamh air an itealan leis a’ chòd aige (gu mì-fhortanach, tha an t-ùghdar air stad a chur air ùrachadh an stòr-dàta, ach tha e fhathast buntainneach). Mar eisimpleir, airson còd 3c5ee2 tha am fiosrachadh a leanas againn:
Deasaich: anns Tha an tuairisgeul air còd ICAO air a thoirt seachad nas mionaidiche; Tha mi a’ moladh gun leugh an fheadhainn le ùidh e.
DÀTA (56 no 112 pìosan) - an fhìor dàta a dhì-chòdaicheas sinn. Is e a’ chiad 5 pìosan dàta an raon Còd seòrsa, anns a bheil an subtype den dàta air a stòradh (gun a bhith air a mheasgadh le DF). Tha grunn de na seòrsaichean sin ann:
()
Bheir sinn sùil air beagan eisimpleirean de phasganan.
Aithneachadh itealain
Eisimpleir ann an cruth binary:
00100 011 000101 010111 000111 110111 110001 111000
Raointean dàta:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, tha còdan anns gach caractar C1-C8 a fhreagras air clàran-amais san loidhne:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ####_##################
Le bhith a’ còdachadh an t-sreang, tha e furasta còd an itealain fhaighinn: 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('#', ''))
Suidheachadh adhair
Ma tha an t-ainm sìmplidh, tha na co-chomharran nas iom-fhillte. Tha iad air an gluasad ann an cruth 2, eadhon agus frèaman neònach. Còd achaidh TC = 01011b = 11.
Eisimpleir de phasgan cothromach is neo-àbhaisteach:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Tha àireamhachadh nan co-chomharran fhèin a’ tachairt a rèir foirmle caran duilich:
()
Chan e eòlaiche GIS a th’ annam, agus mar sin chan eil fhios agam cò às a tha e. Cò aig tha fios, sgrìobh anns na beachdan.
Thathas den bheachd gu bheil àirde nas sìmplidh - a rèir an ìre shònraichte, faodar a riochdachadh mar iomadachadh de 25 no 100 troigh.
Luas an adhair
Pacaid le TC=19. Is e an rud inntinneach an seo gum faod an astar a bhith neo-mhearachdach, an coimeas ris an talamh (Ground Speed), no èadhar, air a thomhas le sensor itealain (Airspeed). Tha grunn raointean eadar-dhealaichte air an sgaoileadh cuideachd:
()
co-dhùnadh
Mar a chì thu, tha teicneòlas ADS-B air a thighinn gu bhith na symbiosis inntinneach, nuair a tha inbhe feumail chan ann a-mhàin do phroifeiseantaich, ach cuideachd do luchd-cleachdaidh àbhaisteach. Ach gu dearbh, bha prìomh àite ann an seo air a chluich leis an teicneòlas as saoire de ghlacadairean didseatach SDR, a leigeas leis an inneal comharran fhaighinn gu litearra le tricead os cionn gigahertz “airson sgillinn.”
Anns an ìre fhèin, gu dearbh, tha tòrr a bharrachd ann. Faodaidh an fheadhainn le ùidh am PDF fhaicinn air an duilleag no tadhal air an fhear a chaidh ainmeachadh gu h-àrd .
Chan eil e coltach gum bi a h-uile rud gu h-àrd feumail do mhòran, ach co-dhiù tha am beachd coitcheann air mar a tha e ag obair, tha mi an dòchas, fhathast ann.
Co-dhiù, tha decoder deiseil ann am Python ann mu thràth, faodaidh tu sgrùdadh a dhèanamh air . Agus faodaidh sealbhadairean luchd-glacaidh SDR decoder deiseil ADS-B a chruinneachadh agus a chuir air bhog , chaidh seo a dheasbad nas mionaidiche ann an .
Tha còd stòr a’ pharser a tha air a mhìneachadh san artaigil air a thoirt seachad fon ghearradh. Is e seo eisimpleir deuchainn nach eil a ’leigeil a-mach gur e cinneasachadh a th’ ann, ach tha cuid de rudan ag obair ann, agus faodar a chleachdadh gus am faidhle a chaidh a chlàradh gu h-àrd a pharsadh.
Còd stòr (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()
Tha mi an dòchas gun robh ùidh aig cuideigin, tapadh leat airson d’ aire.
Source: www.habr.com