Project desert tortoise

From NebarnixWiki
Revision as of 06:19, 26 February 2016 by NebarnixWikiSysop (Talk | contribs) (Added page from cut-pasted parts of the demodulator page)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Problem Statement

I found an embedded data stream in the NOAA DSB broadcasts that contains doppler data from ground-based transmitters. These messages contain short bits of data (up to 64 bytes) and the time vs frequency data required to geolocate them to a km or so (in theory). My friends at the hackerspace kept making jokes about using to find a desert tortoise with a huge transmitter glued to its back, and thus the project name was born...

Work

Orbital Determination

In order to correlate the doppler shifts to geo-locations we need to know the position of the spacecraft at each moment. Matlab does not have readily available tools to do this, so I wrote a quick python script using pyephem which one can query directly from matlab or interface with directly. When given a UTC time and date, the script will spit out instantaneous lat/lon of the satellite as well as orbital period, inclination, and altitude/radius. We will use all of these parameters to initialize and optimize our doppler model.

#!/usr/bin/env python

import time
import datetime
import ephem
import optparse
import sys

parser = optparse.OptionParser()
parser.add_option("-t", "--timestamp", dest="timestamp",
                  help="Date and time to look for stuffs",                  
                  metavar="TIMESTAMP")
parser.add_option("-l", "--lattitude", action="store_const", dest="mode", const=1,
                  help="ONLY return the lattitude at datetime",                  
                  )
parser.add_option("-o", "--longitude", action="store_const", dest="mode", const=2,
                  help="ONLY return the longitude at datetime",                  
                  )
parser.add_option("-r", "--radius", action="store_const", dest="mode", const=3,
                  help="ONLY return the radius of the spacecraft at the datetime",                  
                  )
parser.add_option("-a", "--altitude", action="store_const", dest="mode", const=4,
                  help="ONLY return the altitude of the spacecraft above sealevel at the datetime",                  
                  )
parser.add_option("-p", "--period", action="store_const", dest="mode", const=5,
                  help="ONLY return the orbital period of the spacecraft",                  
                  )
parser.add_option("-i", "--inclination", action="store_const", dest="mode", const=6,
                  help="ONLY return the inclination of the spacecraft",                  
                  )
options,args = parser.parse_args()
if options.timestamp:
    # Try parsing the date argument
    try:
        timestamp = datetime.datetime.strptime(options.timestamp, "%Y-%m-%d %H:%M:%S.%f")
    except:
        print "Error parsing date input:",sys.exc_info()
        sys.exit(1)

#example TLE name ./TLE/noaa2016-02-14.txt
filename = "noaa" + datetime.datetime.strftime(timestamp, "%Y-%m-%d") + ".txt"
try:
   fo = open("/home/pi/pyephem/TLE/" + filename, "r")
except:
   print "TLE file not available for this date :(<br>"
   print filename
   sys.exit(1)
   
line = fo.readline()
while line != "":
   if "NOAA 18" in line:
      line1=line
      line2=fo.readline()
      line3=fo.readline()
      break
   line = fo.readline()
   
fo.close

NOAA = ephem.readtle(line1, line2, line3)
NOAA.compute(timestamp)

pi = 3.14159265359;

if   options.mode == 1:
   print('%f' % ((180.0/pi)*NOAA.sublat))
elif options.mode == 2:
   print('%f' % ((180.0/pi)*NOAA.sublong))
elif options.mode == 3:
   print('%f' % (6378 + NOAA.elevation/1000))
elif options.mode == 4:
   print('%f' % (NOAA.elevation/1000))
elif options.mode == 5:
   print('%f' % (86164.0 / float(line3.split()[7])))
elif options.mode == 6:
   print('%f' % (float(line3.split()[2])))
else:
   print('%f %f %f' % ((180.0/pi)*NOAA.sublat, (180.0/pi)*NOAA.sublong, NOAA.elevation/1000))

Doppler Data

Doppler shifts produce sigmoid functions, so we should see a very specific shape with an inflection point where the satellite passes overhead and relative velocity reaches zero, then reverses. The last three bytes hold the key. The most data I have so far comes from a Baltic sea area transmitter (fishing ship? sea float?) with the ID 0x24000530 (in hex). There is probably enough data here to determine rough location, but more data to characterize the receiver range (units of frequency shift?) would be nice. Also, all doppler calcs will show two identical points for every solution -- multiple passes must be captured in order to eliminate the incorrect solution.

Three separate bytes, doesn't look right
Fixed point offset binary -- now we're talking! Doppler ahoj!


Latest Doppler Shift Model for Geolocating ARGOS Transmitters

A matlab script that interfaces with an updated version of the above python script has been created that does pretty much all of the hard work for you. You must provide Spacecraft ID (15, 18, 19) T0 time and date of the time axis (xdata vs ydata provided) and there is a small calibration factor to the orbital period that you can use to run a second pass if required after the first pass. This script will attempt to find the best position and spit out the latitude and longitude but please be aware that there are always TWO possible solutions, and they become more and more fuzzy the fewer the data points you feed it. Multiple satellite passes are required to discern which location is the correct one!

http://wiki.nebarnix.com/wiki/File:Doppler3DRotationInclination3FitandWebephem.m

Results are shown below, the plotted points represent the solutions from two passes 20 minutes apart (NOAA15 and NOAA18). The distantly spaced models uses the 3 DOF model which assume a doppler scaling of 1.0 (just trust the raw frequency data) whereas the more closely spaced data uses the 4 DOF model that solves for slight calibration errors of the spacecraft. Both models solve for frequency offsets of the transmitter.

DopplerCurveFit.jpg DopplerModelResultsVisualization2.jpg DopplerModelResultsVisualization1.jpg DopplerModelResults.jpg DopplerModelResults3vs4.jpg