Silav Habr. Dibe ku her kesê ku di balafirê de xizm û hevalên xwe dîtiye an dîtiye, karûbarê belaş Flightradar24 bikar aniye. Ev rêgezek pir hêsan e ku meriv pozîsyona balafirê di wextê rast de bişopîne.
В Prensîba xebitandinê ya karûbarek wusa serhêl hate vegotin. Em ê naha pêşde biçin û fêhm bikin ka çi dane ji balafirê berbi stasyona wergirtinê ve têne şandin û wergirtin û bi karanîna Python bi xwe wê deşîfre dikin.
История
Eşkere ye, daneyên balafirê ji bo bikarhêneran nayên şandin ku li ser têlefonên xwe bibînin. Pergalê jê re ADS-B (Automatic dependent surveillance-weşan) tê gotin, û tê bikar anîn ku bixweber agahdariya li ser balafirê ji navenda kontrolê re bişîne - nasnav, koordînat, rêgez, lez, bilindahî û daneyên din têne şandin. Berê, berî hatina pergalên weha, belavker tenê dikaribû xalek li ser radarê bibîne. Dema ku balafir zêde bûn êdî ev têr nedikir.
Ji hêla teknîkî ve, ADS-B ji veguhezkerek li ser balafirê pêk tê ku bi periyodîk pakêtên agahdarî bi frekansek pir zêde ya 1090 MHz dişîne (modulên din jî hene, lê em ew qas bi wan re ne eleqedar in, ji ber ku koordînat tenê li vir têne şandin). Bê guman, ji bilî veguhestinê, li cîhek li balafirgehê wergirek jî heye, lê ji bo me, wekî bikarhêner, wergirê me bixwe balkêş e.
Bi awayê, ji bo berhevdanê, pergala yekem a bi vî rengî, Airnav Radarbox, ku ji bo bikarhênerên asayî hatî sêwirandin, di sala 2007-an de xuya bû, û bi qasî 900 $ lêçû; abonetiya karûbarên torê salê 250 $ din lêçû.
Nirxên wan xwediyên yekem ên rûsî dikarin li ser forumê werin xwendin . Naha ku wergirên RTL-SDR bi berfirehî peyda bûne, amûrek wusa dikare bi 30 $ were berhev kirin; bêtir li ser vê yekê di . Ka em biçin ser protokolê bixwe - em bibînin ka ew çawa dixebite.
Sinyalan distînin
Pêşîn, pêdivî ye ku sînyala were tomar kirin. Tevahiya sînyalê bi tenê 120 mîkroçirkeyan dirêj dike, ji ber vê yekê ji bo ku bi rihetî pêkhateyên wê ji hev veqetînin, wergirê SDR bi frekansa nimûneyê ya herî kêm 5 MHz tê xwestin.
Piştî tomarkirinê, em pelek WAV bi rêjeya nimûneyê 5000000 nimûne / sankî werdigirin; 30 saniyeyên weha tomar "giraniya" bi qasî 500 MB. Guhdariya wê bi lîstikvanek medyayê re, bê guman, bêkêr e - pel deng nagire, lê îşaretek radyoyê rasterast dîjîtalkirî - bi vî rengî Radyoya Defined Software bi vî rengî dixebite.
Em ê pelê bi karanîna Python vekin û pêvajoyê bikin. Kesên ku dixwazin bi tena serê xwe ceribandinê bikin dikarin tomarek nimûne dakêşin .
Ka em pelê dakêşin û bibînin ka hundurê çi ye.
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()
Encam: em li hember dengê paşerojê "pişk"ên eşkere dibînin.
Her "puls" îşaretek e, heke hûn çareseriya li ser grafîkê zêde bikin avahiya wê bi zelalî xuya dibe.
Wekî ku hûn dikarin bibînin, wêne bi ya ku di danasîna li jor de hatî dayîn bi tevahî hevaheng e. Hûn dikarin dest bi hilberandina daneyan bikin.
Deşîfrekirin
Pêşîn, hûn hewce ne ku hinekî tîrêjê bigirin. Nîşan bixwe bi karanîna kodkirina Manchesterê tê kod kirin:
Ji cûdahiya astê di nibbles de hêsan e ku meriv "0" û "1" ya rastîn bistîne.
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"
Struktura sînyalê bixwe wiha ye:
Ka em bi hûrgulî li zeviyan binêrin.
DF (Forma dakêşanê, 5 bit) - celebê peyamê diyar dike. Gelek celeb hene:
()
Em tenê bi tîpa DF17 re eleqedar dibin, ji ber ku ... Ev e ku koordînatên balafirê dihewîne.
ICAO (24 bit) - koda bêhempa ya navneteweyî ya balafirê. Hûn dikarin balafirê bi koda wê kontrol bikin (mixabin, nivîskar nûvekirina databasê rawestandiye, lê ew hîn jî têkildar e). Mînakî, ji bo koda 3c5ee2 me agahdariya jêrîn heye:
Biguherîne: in Danasîna koda ICAO bi hûrgulî tête dayîn; Ez pêşniyar dikim ku kesên eleqedar wê bixwînin.
JIMARE (56 an 112 bit) - Daneyên rastîn ên ku em ê deşîfre bikin. Daneyên yekem 5 bit zevî ne Koda Tîpa, binecureya daneyên ku têne hilanîn vedihewîne (ku bi DF re neyê tevlihev kirin). Gelek ji van celeb hene:
()
Ka em li çend nimûneyên pakêtan binêrin.
Nasnameya balafirê
Mînak di forma binary de:
00100 011 000101 010111 000111 110111 110001 111000
Qadên daneyê:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, her karaktera C1-C8 kodên ku bi nîşaneyên di rêzê de têkildar in vedihewîne:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789#######
Bi deşîfrekirina rêzê, hêsan e ku hûn koda balafirê bistînin: 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('#', ''))
Helwesta hewayê
Ger nav sade be, wê hingê hevrêzî tevlihevtir in. Ew di forma çarçoveyên 2, heta û xerîb de têne veguhestin. Koda zeviyê TC = 01011b = 11.
Mînaka pakêtên jin û mêr:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Hesabkirina koordînatan bixwe li gorî formulek pir dijwar pêk tê:
()
Ez ne pisporê GIS im, ji ber vê yekê ez nizanim ew ji ku tê. Kî dizane, di şîroveyan de binivîse.
Bilindahî hêsantir tê hesibandin - li gorî bitek taybetî, ew dikare wekî pirjimarek 25 an 100 lingan were destnîşan kirin.
Leza hewayê
Pakêta bi TC=19. Tiştê balkêş li vir ew e ku lez dikare an rast be, li gorî erdê (Zemîn Leza), an jî hewa be, ku ji hêla senzorek balafirê ve were pîvandin (Airspeed). Gelek qadên cûda jî têne şandin:
()
encamê
Wekî ku hûn dikarin bibînin, teknolojiya ADS-B bûye sembiozek balkêş, dema ku standardek ne tenê ji pisporan re, lê ji bo bikarhênerên asayî jî bikêr e. Lê bê guman, di vê yekê de rolek sereke ji hêla teknolojiya erzantir a wergirên SDR-ya dîjîtal ve hate lîstin, ku destûrê dide cîhazê ku bi rastî îşaretan bi frekansên li jor gigahertz "ji bo peran" werbigire.
Di standard bixwe de, bê guman, pir zêde heye. Kesên ku dixwazin dikarin li ser rûpelê PDF-ê bibînin an serdana yê ku berê li jor hatî behs kirin .
Ne mimkûn e ku hemî ya jorîn ji gelek kesan re bikêr be, lê bi kêmanî ramana gelemperî ya çawa dixebite, ez hêvî dikim, dimîne.
Bi awayê, di Python de dekoderek amade jixwe heye, hûn dikarin wê bixwînin . Û xwedan wergirên SDR dikarin dekoderek ADS-B ya amade bicivînin û bidin destpêkirin , li ser vê yekê bi berfirehî hate nîqaş kirin .
Koda çavkaniyê ya parserê ku di gotarê de hatî destnîşan kirin li jêr qutkirî tê dayîn. Ev mînakek ceribandinê ye ku wekî hilberandinê naxuye, lê hin tişt tê de dixebitin, û ew dikare were bikar anîn da ku pelê ku li jor hatî tomarkirin pars bike.
Koda çavkaniyê (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()
Ez hêvî dikim ku kesek eleqedar bû, spas ji bo baldariya we.
Source: www.habr.com