Project desert tortoise
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.
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.
