Project desert tortoise

From NebarnixWiki
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

Doppler Model

The model is complex enough to need its own page... [1]

Time

DCS/2 (NOAA-15 and NOAA-18) appear to have less precise hardware than ADCS (NOAA-19). No correction factor seems needed for ADCS, but one is needed for DCS/2 spacecraft For example:

TXID:FB3AD000 SPID:NOAA-19 ADCS
ID	Local Time	ADCS Time	Freq       	LTdelta	STdelta
196	111.7084	10060.62	-89901.75000	50.0525	**50.00
414	161.7609	10110.62	-93257.56250	49.9911	**50.00
611	211.7520	10160.62	-97425.96875	49.9558	**50.00
838	261.7078	10210.62	-101951.8438	50.531	**50.00
1030	312.2388	10260.62	-106160.9063	99.4706	**100.00
1402	411.7094	10360.62	-112104.6875	50.0451	**50.00
1582	461.7545	10410.62	-113888.6563	50.0914	**50.00
1755	511.8459	10460.62	-115123.4375		
TXID:FB3AD000 SPID:NOAA-18 DCS/2
Raw, and with Correction Factor of 0.95494
ID	LTime	        Stime	SDT	CSTime	        SCDT     	LDT
117	169.1584	5066.25	*104.71	4837.964775	99.99177	99.9318
391	269.0902	5170.96	*52.36	4937.956542	50.00066	50.1962
586	319.2864	5223.32	*52.36	4987.957201	50.00066	49.8038
757	369.0902	5275.68	*104.7	5037.957859	99.98222	99.9971
1110	469.0873	5380.38	*52.37	5137.940077	50.01021	50.0827
1276	519.1700	5432.75	*52.36	5187.950285	50.00066	50.0039
1427	569.1739	5485.11		5237.950943		
TXID:FB3AD000 SPID:NOAA-15 DCS/2
Raw, and with Correction Factor of 0.95494
ID	LTime   	Stime   	SDT	CSTime     	SCDT       	LDT
418	281.9768	12293.59	52.35	11739.64083	49.991109	50.0000
576	331.9768	12345.94	52.37	11789.63194	50.0102078	50.0000
728	381.9768	12398.31	52.36	11839.64215	50.0006584	50.0000
886	431.9768	12450.67	52.33	11889.64281	49.9720102	50.0356
1025	482.0124	12503.00	52.36	11939.61482	50.0006584	49.9644
1163	531.9768	12555.36		11989.61548		

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. Some good data I have comes from a Baltic sea area transmitter (fishing ship? sea float?) with the ID 0x24000530. By analyzing the bytes, we can figure out how to combine these three bytes to produce the frequency as received by the spacecraft.

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, and the time and frequency data. There is a function to fine tune the orbital period and and satellite radius to minimize the error directly at the frequency inflection point by providing a fudge factor for each (which is different for every solution set). This script will attempt to find the best position and spit out the latitude and longitude. 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

Result 1 Nov 2015

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 errors of the curves (but where are the errors coming from? Ionisphere? spacecraft calibration?). Both models solve for frequency offsets of the transmitter.
DopplerCurveFit.jpg DopplerModelResultsVisualization2.jpg DopplerModelResultsVisualization1.jpg DopplerModelResults.jpg DopplerModelResults3vs4.jpg

Updated results 11 Nov 2015

I learned that geocentric latitude (what the model and ephemris uses) is NOT the same as geodetic latitude, which is what google maps uses (which is because the earth is an oblate spheroid). The (simplest) formula to convert is (it gets more complex as you start getting super accurate, like the WGS84 datum that your GPS uses)

GeodLat = atand(tand(GeoCentricLat)/(1-0.00666667))

When you correct for the difference it shifts everything north by about 20 arcminutes.....

Phoenix Transmitter.png

Result 2 Feb 2016

First discovery of actual location of transmitter occured while driving along Chincoteague road by the NOAA facility on the NASA Wallops base. The decoded bits show an ID 01D81000, which was found in in the embedded telemetry DCS stream. I used this to fine-tune my model, also this is the first test to use WGS95 elliptical earth model.

WFF DCS5.jpg WFF DCS3.jpg WFF DCS4.jpg WFF DCS2.jpg WFF DCS1.jpg

List of Discovered Transmitters

68961013 - NASA Goddard, Maryland, or in the area.
01B90000 - NASA Goldstone, around (35.335948, -116.806606)
FB3AD000 - Telonics, Mesa, AZ (33.384355, -111.811219)
01D81000 - NASA Wallops NESDIS (37.94648, -75.46059)