Habari Habr. Huenda kila mtu ambaye amewahi kukutana au kuonana na jamaa au marafiki kwenye ndege ametumia huduma ya bure ya Flightradar24. Hii ni njia rahisi sana ya kufuatilia nafasi ya ndege kwa wakati halisi.
В Kanuni ya uendeshaji wa huduma hiyo ya mtandaoni ilielezwa. Sasa tutaendelea na kufahamu ni data gani inatumwa na kupokewa kutoka kwa ndege hadi kituo cha kupokelea na kuzitatua sisi wenyewe kwa kutumia Chatu.
Hadithi
Ni wazi, data ya ndege haisambazwi kwa watumiaji kuona kwenye simu zao mahiri. Mfumo huo unaitwa ADS-B (Ufuatiliaji tegemezi wa Kiotomatiki-matangazo), na hutumiwa kusambaza kiotomatiki habari kuhusu ndege hadi kituo cha udhibiti - kitambulisho chake, viwianishi, mwelekeo, kasi, urefu na data zingine hupitishwa. Hapo awali, kabla ya ujio wa mifumo hiyo, mtumaji aliweza kuona tu uhakika kwenye rada. Hii haikutosha tena wakati kulikuwa na ndege nyingi.
Kitaalam, ADS-B ina kipeperushi kwenye ndege ambayo hutuma pakiti za habari mara kwa mara kwa masafa ya juu ya 1090 MHz (kuna njia zingine, lakini hatupendezwi nazo, kwani kuratibu hupitishwa hapa tu). Kwa kweli, pamoja na kipeperushi, pia kuna mpokeaji mahali pengine kwenye uwanja wa ndege, lakini kwetu sisi, kama watumiaji, mpokeaji wetu mwenyewe anavutia.
Kwa njia, kwa kulinganisha, mfumo wa kwanza kama huo, Airnav Radarbox, iliyoundwa kwa watumiaji wa kawaida, ulionekana mnamo 2007, na gharama ya karibu $ 900; usajili kwa huduma za mtandao uligharimu $250 nyingine kwa mwaka.
Mapitio ya wamiliki hao wa kwanza wa Kirusi yanaweza kusomwa kwenye jukwaa . Kwa kuwa vipokezi vya RTL-SDR vimepatikana kwa wingi, kifaa kama hicho kinaweza kukusanywa kwa $30; zaidi kuhusu hili . Hebu tuendelee kwenye itifaki yenyewe - hebu tuone jinsi inavyofanya kazi.
Kupokea ishara
Kwanza, ishara inahitaji kurekodi. Ishara nzima ina muda wa microseconds 120 tu, hivyo ili kutenganisha vipengele vyake kwa urahisi, mpokeaji wa SDR na mzunguko wa sampuli wa angalau 5 MHz ni kuhitajika.
Baada ya kurekodi, tunapokea faili ya WAV yenye kiwango cha sampuli ya sampuli 5000000/sekunde; sekunde 30 za rekodi kama hiyo "uzani" wa takriban MB 500. Kuisikiliza na kicheza media, kwa kweli, haina maana - faili haina sauti, lakini ishara ya redio ya dijiti moja kwa moja - hii ndio jinsi Redio Iliyofafanuliwa ya Programu inavyofanya kazi.
Tutafungua na kusindika faili kwa kutumia Python. Wale ambao wanataka kujaribu wenyewe wanaweza kupakua rekodi ya mfano .
Hebu kupakua faili na kuona nini ndani.
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()
Matokeo: tunaona "mipigo" dhahiri dhidi ya kelele ya chinichini.
Kila "pulse" ni ishara, muundo ambao unaonekana wazi ikiwa unaongeza azimio kwenye grafu.
Kama unaweza kuona, picha ni sawa kabisa na kile kilichotolewa katika maelezo hapo juu. Unaweza kuanza kuchakata data.
Kusimbua
Kwanza, unahitaji kupata mkondo kidogo. Ishara yenyewe imesimbwa kwa kutumia usimbaji wa Manchester:
Kutoka kwa tofauti ya ngazi katika nibbles ni rahisi kupata halisi "0" na "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"
Muundo wa ishara yenyewe ni kama ifuatavyo.
Hebu tuangalie mashamba kwa undani zaidi.
DF (Downlink Format, 5 bits) - huamua aina ya ujumbe. Kuna aina kadhaa:
()
Tunavutiwa na aina ya DF17 pekee, kwa sababu... Ni hii ambayo ina kuratibu za ndege.
ICAO (Biti 24) - nambari ya kipekee ya kimataifa ya ndege. Unaweza kuangalia ndege kwa nambari yake (kwa bahati mbaya, mwandishi ameacha kusasisha hifadhidata, lakini bado inafaa). Kwa mfano, kwa msimbo 3c5ee2 tunayo habari ifuatayo:
Hariri: katika Maelezo ya msimbo wa ICAO yametolewa kwa undani zaidi; ninapendekeza wale wanaopenda waisome.
DATA (Biti 56 au 112) - data halisi ambayo tutaamua. Biti 5 za kwanza za data ni uwanja Nambari ya Aina, iliyo na aina ndogo ya data inayohifadhiwa (isichanganywe na DF). Kuna baadhi ya aina hizi:
()
Hebu tuangalie mifano michache ya vifurushi.
Kitambulisho cha ndege
Mfano katika fomu ya binary:
00100 011 000101 010111 000111 110111 110001 111000
Sehemu za data:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, kila herufi C1-C8 ina misimbo inayolingana na fahirisi kwenye mstari:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_################0123456789#######
Kwa kusimbua kamba, ni rahisi kupata nambari ya ndege: 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('#', ''))
Nafasi ya hewa
Ikiwa jina ni rahisi, basi kuratibu ni ngumu zaidi. Wao hupitishwa kwa namna ya 2, muafaka hata na isiyo ya kawaida. Msimbo wa uga TC = 01011b = 11.
Mfano wa pakiti sawa na isiyo ya kawaida:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Hesabu ya kuratibu yenyewe hufanyika kulingana na fomula ngumu zaidi:
()
Mimi sio mtaalam wa GIS, kwa hivyo sijui inatoka wapi. Nani anajua, andika kwenye maoni.
Urefu unachukuliwa kuwa rahisi - kulingana na biti maalum, inaweza kuwakilishwa kama kizidishi cha futi 25 au 100.
Kasi ya Anga
Kifurushi chenye TC=19. Jambo la kufurahisha hapa ni kwamba mwendo unaweza kuwa sahihi, ukilinganisha na ardhi (Ground Speed), au airborne, inayopimwa na sensor ya ndege (Airspeed). Sehemu nyingi tofauti pia hupitishwa:
()
Hitimisho
Kama unaweza kuona, teknolojia ya ADS-B imekuwa symbiosis ya kuvutia, wakati kiwango ni muhimu sio tu kwa wataalamu, bali pia kwa watumiaji wa kawaida. Lakini kwa kweli, jukumu muhimu katika hili lilichezwa na teknolojia ya bei nafuu ya vipokeaji vya SDR ya dijiti, ambayo inaruhusu kifaa kupokea mawimbi kihalisi na masafa juu ya gigahertz "kwa senti."
Katika kiwango yenyewe, bila shaka, kuna mengi zaidi. Wale wanaovutiwa wanaweza kutazama PDF kwenye ukurasa au tembelea ile iliyotajwa hapo juu .
Haiwezekani kwamba yote yaliyo hapo juu yatakuwa na manufaa kwa wengi, lakini angalau wazo la jumla la jinsi inavyofanya kazi, natumaini, bado.
Kwa njia, decoder iliyotengenezwa tayari katika Python tayari ipo, unaweza kuisoma . Na wamiliki wa vipokezi vya SDR wanaweza kukusanyika na kuzindua avkodare iliyotengenezwa tayari ya ADS-B , hii ilijadiliwa kwa undani zaidi katika .
Nambari ya chanzo cha kichanganuzi kilichoelezewa katika kifungu kinapewa chini ya kata. Huu ni mfano wa majaribio ambao haujifanyi kuwa uzalishaji, lakini baadhi ya mambo hufanya kazi ndani yake, na inaweza kutumika kuchanganua faili iliyorekodiwa hapo juu.
Msimbo wa chanzo (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()
Natumai mtu alipendezwa, asante kwa umakini wako.
Chanzo: mapenzi.com