Flightradar24 - እንዴት ነው የሚሰራው? ክፍል 2፣ ADS-B ፕሮቶኮል

ሰላም ሀብር። ምናልባትም ቢያንስ አንድ ጊዜ ከዘመዶቻቸው ወይም ከጓደኞቻቸው ጋር በአውሮፕላን የተገናኙ ወይም ያዩ ሰዎች የ Flightradar24 አገልግሎትን ተጠቅመዋል። ይህ የአውሮፕላኑን አቀማመጥ በእውነተኛ ጊዜ ለመከታተል በጣም ምቹ መንገድ ነው.

В የመጀመሪያው ክፍል። የእንደዚህ አይነት የመስመር ላይ አገልግሎት አሠራር መርህ ተገልጿል. አሁን የበለጠ እንሄዳለን እና ከአውሮፕላኑ ወደ መቀበያ ጣቢያ የሚተላለፈውን እና የሚቀበለውን መረጃ ለማወቅ እና ፒቲንን በመጠቀም እራሳችንን ዲኮድ እናደርጋለን።

История

በግልጽ ለማየት እንደሚቻለው የአውሮፕላን መረጃ ተጠቃሚዎች በስማርት ስልኮቻቸው ላይ እንዲያዩት እየተጋራ አይደለም። ስርዓቱ ADS-B (ራስ-ሰር ጥገኛ ክትትል-ብሮድካስት) ተብሎ ይጠራል, እና ስለ አውሮፕላኑ መረጃን ወደ መቆጣጠሪያ ማእከል በራስ-ሰር ለማስተላለፍ ጥቅም ላይ ይውላል - መለያው ፣ መጋጠሚያዎች ፣ አቅጣጫ ፣ ፍጥነት ፣ ከፍታ እና ሌሎች መረጃዎች ይተላለፋሉ። ከዚህ ቀደም እንደዚህ አይነት ስርዓቶች ከመምጣቱ በፊት ላኪው በራዳር ላይ አንድ ነጥብ ብቻ ማየት ይችላል. ብዙ አውሮፕላኖች ሲኖሩ ይህ በቂ አልነበረም።

በቴክኒካዊ ሁኔታ ኤ.ዲ.ኤስ-ቢ የአውሮፕላን አስተላላፊን ያካትታል ፣ መረጃን በየጊዜው በ 1090 ሜኸር ድግግሞሽ (ሌሎች ሁነታዎች አሉ ፣ ግን ለእኛ በጣም አስደሳች አይደሉም ፣ ምክንያቱም መጋጠሚያዎቹ የሚተላለፉት እዚህ ብቻ ነው)። በእርግጥ ከማስተላለፊያው በተጨማሪ በአውሮፕላን ማረፊያው ውስጥ የሆነ ቦታ ተቀባይ አለ, ነገር ግን ለእኛ እንደ ተጠቃሚዎች, የራሳችን መቀበያ ፍላጎት ነው.

በነገራችን ላይ፣ ለማነፃፀር፣ የመጀመሪያው እንዲህ ዓይነቱ ስርዓት፣ ለመደበኛ ተጠቃሚዎች የተነደፈው Airnav Radarbox፣ በ2007 ታየ እና ወደ 900 ዶላር ገደማ ያስወጣል፣ እና በግምት። 250$ በዓመት ለኔትወርክ አገልግሎቶች የደንበኝነት ምዝገባ ክፍያ።

የእነዚያ የመጀመሪያዎቹ የሩሲያ ባለቤቶች ግምገማዎች በመድረኩ ላይ ሊነበቡ ይችላሉ። ሬዲዮ ስካነር. አሁን የ RTL-SDR ተቀባዮች በብዛት ይገኛሉ ፣ ተመሳሳይ መሳሪያ በ 30 ዶላር ሊገጣጠም ይችላል ፣ የበለጠ ስለዚህ ጉዳይ በ ውስጥ ነበር ። የመጀመሪያው ክፍል።. ወደ ፕሮቶኮሉ ራሱ እንሸጋገራለን - እንዴት እንደሚሰራ እንይ.

የምልክት መቀበያ

በመጀመሪያ, ምልክቱ መቅዳት አለበት. ሙሉ ሲግናል የሚቆይበት ጊዜ 120 ማይክሮ ሰከንድ ብቻ ነው ስለዚህ ክፍሎቹን በምቾት ለመበተን ቢያንስ 5 ሜኸር የናሙና መጠን ያለው ኤስዲአር ተቀባይ ተፈላጊ ነው።

ከተቀረጸ በኋላ የ WAV ፋይልን በናሙና 5000000 ናሙና / ሰከንድ እናገኛለን ፣ እንደዚህ ያለ ቀረጻ 30 ሴኮንድ 500 ሜባ ያህል “ይመዝናል ። እርግጥ ነው, በሚዲያ ማጫወቻ ማዳመጥ ምንም ፋይዳ የለውም - ፋይሉ ድምጽ የለውም, ነገር ግን በቀጥታ ዲጂታይዝ የተደረገ የሬዲዮ ምልክት - በዚህ መንገድ ነው የሶፍትዌር ዲፋይንድ ራዲዮ የሚሰራው.

ፒቲንን በመጠቀም ፋይሉን ከፍተን እናስኬዳለን። በራሳቸው መሞከር የሚፈልጉ የምሳሌ ቀረጻ ማውረድ ይችላሉ። ማያያዣ.

ፋይሉን እናውርድ እና ውስጥ ያለውን እንይ።

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()

ውጤት፡ በድምፅ ዳራ ላይ ግልጽ የሆኑ "ግፊቶችን" እናያለን።

እያንዳንዱ "ግፊት" ምልክት ነው, በግራፉ ላይ ያለውን ጥራት ከጨመሩ አወቃቀሩ በግልጽ ይታያል.

እንደሚመለከቱት, ስዕሉ ከላይ ባለው መግለጫ ውስጥ ከተሰጠው ጋር በጣም የሚስማማ ነው. ውሂቡን ማካሄድ መጀመር ይችላሉ።

መፍታት

በመጀመሪያ, ትንሽ ዥረት ማግኘት ያስፈልግዎታል. ምልክቱ ራሱ የማንቸስተር ኢንኮዲንግ በመጠቀም ነው የተመሰጠረው፡-

ከኒብል ደረጃ ልዩነት እውነተኛ "0" እና "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"

የምልክቱ አወቃቀር ራሱ እንደሚከተለው ነው-

መስኮችን በበለጠ ዝርዝር እንመልከታቸው.

DF (Downlink ቅርጸት፣ 5 ቢት) — የመልእክቱን አይነት ይገልጻል። በርካታ ዓይነቶች አሉ:

(የጠረጴዛ ምንጭ)

እኛ በ DF17 ዓይነት ላይ ብቻ ፍላጎት አለን, ምክንያቱም የአውሮፕላኑን መጋጠሚያዎች ይዟል.

ICAO (24 ቢት) - የአውሮፕላኑ ዓለም አቀፍ ልዩ ኮድ። አውሮፕላኑን በኮዱ ማረጋገጥ ይችላሉ። የመስመር ላይ (እንደ አለመታደል ሆኖ ደራሲው የውሂብ ጎታውን ማዘመን አቁሟል, ግን አሁንም ጠቃሚ ነው). ለምሳሌ፣ ለኮዱ 3c5ee2 የሚከተለው መረጃ አለን።

አርትዕ፡ ውስጥ በጽሑፉ ላይ አስተያየቶች የ ICAO ኮድ መግለጫ በበለጠ ዝርዝር ተሰጥቷል, ፍላጎት ያላቸው እንዲያነቡት እመክራለሁ.

የውሂብ (56 ወይም 112 ቢት) - የምንፈታው ትክክለኛ ውሂብ። የመጀመሪያዎቹ 5 ቢት ዳታዎች መስክ ናቸው። ኮድ ይተይቡየተከማቸ ውሂብ ንዑስ ዓይነት የያዘ (ከዲኤፍ ጋር መምታታት የለበትም)። ከእነዚህ ዓይነቶች ውስጥ በጣም ጥቂት ናቸው-

(የጠረጴዛ ምንጭ)

የጥቅል ምሳሌዎችን እንመልከት።

የአውሮፕላን መለያ

ሁለትዮሽ ምሳሌ፡-

00100 011 000101 010111 000111 110111 110001 111000

የውሂብ መስኮች:

+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+

TC = 00100b = 4፣ እያንዳንዱ ቁምፊ C1-C8 በሕብረቁምፊው ውስጥ ካሉ ኢንዴክሶች ጋር የሚዛመዱ ኮዶችን ይይዛል።
#ABCDEFGHIJKLMNOPQRSTUVWXYZ############################0123456789#####

ሕብረቁምፊውን ከገለበጥን በኋላ የአውሮፕላኑን ኮድ ለማግኘት ቀላል ነው፡ 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('#', ''))

የአየር ወለድ አቀማመጥ

ሁሉም ነገር በስሙ ቀላል ከሆነ ፣ ከዚያ ከመጋጠሚያዎች ጋር የበለጠ የተወሳሰበ ነው። እንደ 2x፣ እንግዳ እና አልፎ ተርፎም ፍሬሞች ይተላለፋሉ። የመስክ ኮድ TC = 01011b = 11.

የእኩል እና ያልተለመዱ እሽጎች ምሳሌ፡-

01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000

የመጋጠሚያዎቹ ስሌት ራሱ የሚከናወነው በጣም አስቸጋሪ በሆነ ቀመር መሠረት ነው-

(ምንጩ)

እኔ የጂአይኤስ ስፔሻሊስት አይደለሁም፣ ስለዚህ ከየት እንደመጣ አላውቅም። ማን ያውቃል, በአስተያየቶቹ ውስጥ ይፃፉ.

ከፍታ እንደ ቀላል ይቆጠራል - እንደ ልዩ ቢት ፣ እንደ 25 ወይም 100 ጫማ ብዜት ሊወከል ይችላል።

የአየር ወለድ ፍጥነት

ጥቅል ከTC=19 ጋር። እዚህ ላይ የሚያስደንቀው ነገር ፍጥነቱ ከመሬት አንጻር (የመሬት ፍጥነት) እና አየር በአውሮፕላኑ ዳሳሽ (Airspeed) የሚለካው ሁለቱም ትክክለኛ ሊሆኑ ይችላሉ። ብዙ የተለያዩ መስኮችም አልፈዋል፡-

(ምንጩ)

መደምደሚያ

እንደሚመለከቱት ፣ የ ADS-B ቴክኖሎጂ አስደሳች ሲምባዮሲስ ሆኗል ፣ አንድ ደረጃ ለባለሙያዎች ብቻ ሳይሆን ለተለመዱ ተጠቃሚዎችም ጠቃሚ ነው። ግን በእርግጥ በዚህ ውስጥ ቁልፍ ሚና የተጫወተው በዲጂታል ኤስዲአር ተቀባዮች ርካሽ ቴክኖሎጂ ነው ፣ ይህም መሣሪያው በጥሬው “ለአንድ ሳንቲም” ከአንድ ጊጋኸርትዝ በላይ ድግግሞሽ ያላቸውን ምልክቶች እንዲቀበል ያስችለዋል።

በደረጃው በራሱ, በእርግጥ, ብዙ ተጨማሪ. ፍላጎት ያላቸው በገጹ ላይ ያለውን ፒዲኤፍ ማየት ይችላሉ። ICAO ወይም ከላይ የተጠቀሱትን ይጎብኙ ድር ጣቢያ.

ከላይ ያሉት ሁሉም ነገሮች ለብዙዎች ጠቃሚ ይሆናሉ ተብሎ የማይታሰብ ነገር ነው ፣ ግን ቢያንስ እንዴት እንደሚሰራ አጠቃላይ ሀሳብ ፣ ተስፋ አደርጋለሁ ፣ ይቀራል።

በነገራችን ላይ, ዝግጁ የሆነ የ Python ዲኮደር አስቀድሞ አለ, ሊያጠኑት ይችላሉ እዚህ. እና የኤስዲአር መቀበያዎች ባለቤቶች ዝግጁ የሆነ ADS-B ዲኮደርን ሰብስበው ማስኬድ ይችላሉ። ከገጽ፣ ይህ በ ውስጥ የበለጠ በዝርዝር ተብራርቷል የመጀመሪያው ክፍል።.

በአንቀጹ ውስጥ የተገለጸው የትንታኑ ምንጭ ኮድ በቆራጩ ስር ተሰጥቷል። ይህ እንደ ምርት የማይመስል የሙከራ ምሳሌ ነው, ነገር ግን አንድ ነገር በውስጡ ይሠራል, እና ከላይ የተቀዳውን ፋይል ሊተነተኑ ይችላሉ.
ምንጭ ኮድ (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()

አንድ ሰው ፍላጎት እንደነበረው ተስፋ አደርጋለሁ፣ ስለ እርስዎ ትኩረት እናመሰግናለን።

ምንጭ: hab.com