GRAIL Level 1 Data
Algorithm Theoretical Basis Document
Nate Harvey
Eugene Fahnestock
Daniel Kahan
Alexander Konopliv
Gerhard Kruizinga
Kamal Oudrhiri
Meegyeong Paik
DahNing Yuan
Sami Asmar
Mike Watkins
April 1, 2014
JPL Document D75862
III. Time Coordination Packets
GRAIL USO Frequency Determination
Orbit Determination and Propagation
II. Other (NonGrail) Documents
APPENDIX A. GRAIL Timing Flow Diagram
APPENDIX B. Level 1A for GRAILA (INCLUDES ORBIT DETERMINATION AND PROPAGATION INPUT PREP)
APPENDIX C. Level 1A à Level 1B for GRAIL A and B
APPENDIX D. KBR1A à KBR1C and KBR CAL GRAIL A
APPENDIX E. Orbit Determination & Propagation
Gravity Recovery and Interior Laboratory (GRAIL) Level 1 data processing involves:
1) Telemetry preprocessing
2) Clocktiming correlation
3) Biased dual oneway intersatellite range computation
4) Orbit determination
This document outlines Level 1 algorithms for telemetry preprocessing, clocktiming correlation, intersatellite range computation, and orbit determination. For Level 1 software flow diagrams, see Appendices B, C, and D.
The GRAIL mission processes data from two spacecraft, GRAILA and GRAILB (formally called Ebb and Flow), flying in formation orbits around the moon, with a payload that precisely measures the distance between GRAILA and GRAILB, from which we solve for a spherical harmonic expansion of the lunar gravity field. Each GRAIL spacecraft consists of:
1) Rectangular bus
2) Fixed (after deployment) solar panels
3) Titanium diaphragm fuel tank
4) Ultra stable oscillator (USO), which drives onboard science clocks
5) Attitude control system (ACS), consisting of:
a. Four reaction wheels to change attitude
b. Inertial Measurement Unit (IMU) to measure the rate components of angular rotation
c. Star Tracker to measure the absolute attitude
d. Sun Sensor
e. Four pairs of torque thrusters
6) Main engine thruster
7) Two low gain antennas (LGA) for Sband communication with the Deep Space Network (DSN)
8) A Radio Science Beacon (RSB) with two antennas, which transmits oneway Xband carrier
9) Sband interspacecraft Time Transfer System (TTS)
10) Kaband carrier phase tracking interspacecraft receiver/transmitter
Our description of Level 1 processing runs through:
1) Absolute Timing – Correlate three clocks onboard each GRAIL satellite (LGRS, BTC, RTC) with two universal timing standards (TDB, UTC), based on PulsePerSecond (PPS) signals that sync internal clocks and Time Correlation (TC) packets linking onboard and universal time.
2) Telemetry PreProcessing – Create Level 0 data from GRAIL packets transmitted to the ground. Channelized queries provide access to most data types, but science and engineering packets contain Blackjack format data and time correlation packets follow a special format.
3) USO Frequency Determination – Deep Space Network (DSN) stations measure GRAIL UltraStableOscillator (USO) frequencies based on Xband signals.
4) Relative Timing – An Sband Time Transfer System (TTS) measures GRAILAGRAILB relative timing. For greater accuracy, we apply carrier phase smoothing to TTS measurements, a standard algorithm described in a separate section.
5) Clock Model – Assimilate absolute timing, USO frequency, and relative timing measurements in a noncausal Kalman filter, which produces an estimate of GRAILA LGRS and GRAILB LGRS deviations from TDB.
6) KBR Processing – A threestep process flags phase breaks (KBR_debreak), converts data from LGRS to TDB (KBR_order), and produces biased Dual OneWay Range (DOWR) measurements (KBR_compress). Separate sections describe timetag conversion by Lagrangian data resampling and CRN filters to remove highfrequency DOWR noise.
Note that most steps in our Level 1 processing depend directly on the accuracy of:
7) Orbital Solutions – Assimilate oneway frequency data, twoway Doppler, and intersatellite rangerate data in a Kalman filter, running JPL’s Multiple Interferometric Ranging Analysis Using GIPSY Ensemble (MIRAGE) software, solving for dynamic and kinematic spacecraft parameters.
For ease of notation, this section discusses GRAILA timing in isolation. GRAILB timing runs the same algorithms in parallel.
GRAILA absolute timing requires coordination of five clocks:
1) LGRS: Lunar Gravity Ranging System clock. Clock for onboard Kaband (KBR), Xband (RSB), and Sband (TTS) instruments. Driven by an UltraStable Oscillator (USO). Set to 0 when booted. Very stable clock.
2) BTC: Base Time Clock. Onboard satellite clock. Also known as Spacecraft Clock (SCLK). Pretty poor clock, comparable to a wristwatch. More or less synced to UTC at launch time.
3) RTC: Real Time Clock. Flight software clock. Set to 0 when booted. Very unstable clock.
4) UTC: Coordinated Universal Time. Coordinate time in a geocentric reference frame.
5) TDB: Barycentric Dynamical Time. Coordinate time in a solarsystem barycentric frame.
Since the onboard satellite clock drifts, every few weeks the GRAIL spacecraft team solves for a satellite clock correction. Subsequently, the GRAIL spacecraft team produces a corrected clock, the spacecraft event time (SCET), from SCLK by linear extrapolation.
We produce seven Level 1A products correlating clocks.
Product 
From 
To 
TC11A 

BTC arrival of LGRS PPS signal 
TC21A 
LGRS 
BTC 
TC31A 
BTC 
RTC 
TC41A 
LGRS 
RTC 
TC51A 
RTC 
UTC 
TC61A 
UTC 
TDB 
CLK1A 
LGRS 
TDB 
TC21A – TC61A and CLK1A provide a direct relationship between clocks. TC11A, on the other hand, contains information used to create the TC21A product, but does not directly translate LGRS to BTC. For diagrams, see Appendix A.
In these Level 1A products, we do not adjust onboard clocks for relativistic time dilation. We define proper time as the time which a perfect clock onboard GRAILA would measure, including dilation. LGRS, BTC, and RTC reflect (GRAILA proper time + clock drift).
A USO drives the LGRS clock, running at a reference frequency of 4.832 MHz. In addition to running a clock, the LGRS counts USO cycles, and emits a pulsepersecond (PPS) signal after 4.832e6 cycles have passed.
Each PPS signal follows two independent routes (see Figure C in Appendix A):
1) LGRS counter of PPS signals.
a. Starting with the 43^{rd} pulse, at every tenth pulse the LGRS generates an LGRS timing message containing the pulse number, and sends this timing message to the onboard Command and Data Handling (CD&H) satellite computer, to be processed by flight software. This timing message travels from the LGRS counter to the CD&H in about 30 ms.
b. GRAILA flight software knows the RTC time of message arrival, and calculates an approximate BTC (~BTC) time. The CD&H creates an LGRS Standard Formatted Data Unit (SFDU) ([Clark, 2008], [Wilson, 2003], [Wilson, 2012]) which takes the format:
Header 
~BTC 
RTC 
Pulse number = LGRS pulse time, Other stuff 
c. The CD&H sends this LGRS SFDU to the transmitter, to be included in a telemetry frame and sent to Earth.
2) CD&H Module Interface Card (CMIC) clock cycle counter in the CD&H.
a. The CMIC counts BTC clock cycles. Each BTC cycle = 15 ms.
b. The LGRS sends a copy of each PPS to the CMIC, with a travel time on the close order of 30 ms.
c. Upon reception of a PPS, the CMIC clock cycle counter resets to 0.
d. Every 5 seconds, the CD&H creates an SFDU containing the information:
RTC 
BTC 
Number of BTC clock cycles since last PPS received 
e. The CD&H sends this SFDU to the transmitter, to be included in a telemetry frame and sent to Earth. We can access this information by a channelized query to the Ground Data System (GDS). We include the information on BTC time and BTC clock cycles since the last PPS in our TC11A product.
We expect the LGRS to reboot numerous times over the course of the GRAIL mission. Whenever the LGRS reboots, the LGRS PPS counter will reset to 0.
Additionally, every ten minutes, geometry permitting, GRAILA sends a time correlation (TC) packet with an RTC reading from the onboard computer to the onboard transmitter, to be included in a telemetry frame and sent to Earth. Reception at Earth occurs after three delays:
a. The TC packet moves from the onboard computer to the transmitter in a few centiseconds.
b. The time required to buffer the TC packet into a telemetry frame at the transmitter depends on transmission bitrate, but can be more than a second.
c. Lighttravel time from GRAILA to Earth equals about 1.3 seconds.
We can calculate (c) accurately, but (a) and (b) are likely to vary substantially from packet to packet, limiting our knowledge of RTCUTC.
Since GRAILA does not create our various time observables synchronously, we interpolate to form pseudosynchronous measurements, from which we derive our time products. We indicate interpolation with dashed lines in Figure A of Appendix A, and direct observables with solid lines.
LGRS®BTC®RTC:
1) The Level 0 product STC00 contains a list of BTC time and BTC clock cycles for every fifth LGRS PPS, which we include in TC11A.
2) STC00 also lists RTC, from which we derive a BTC/RTC relationship for TC31A.
3) Each LGRS SFDU lists ~BTC, RTC, and LGRS time.
4) Since the approximation ~BTC is fairly accurate, we can match information from TC11A and LGRS SFDUs to coordinate LGRS and BTC in TC21A.
5) We process information from TC21A and TC31A to coordinate LGRS and RTC in TC41A.
RTC®UTC:
TC51A lists RTC times from each TC packet and UTC packet arrival time at end of packet. Note that these times have not been corrected for unknown hardware and software delays or light time. We do solve for a bit rate delay:
where we read the bit rate r from the secondary Compressed Header Data Object (CHDO) ([Clark, 2008]) and the packet length L from the Threshold field of the data CHDO (defined below). The encode factor E = 1116/1279 by definition.
Since bit rate affects hardware and software delays, we also apply a bitrate dependent empirical correction C.
Date rate (kbps) 
C GRAILA (s) 
C GRAILB (s) 
1 
1.3747 
1.3717 
2 
0.6940 
0.6870 
4 
0.3477 
0.38775 
8 
0.17696 
0.17802 
16 
0.09161 
0.09270 
32 
0.048925 
0.0500 
64 
0.02759 
0.0287 
128 
0.01690 
0.01802 
We read end of packet Earth receive time from the SFDU quaternary CHDO, and compute SCLK beginningofpacket time:
where SCLK_{obs} is the SCLK time contained in our TC packet.
UTC®TDB:
TC61A contains a table of standard TDB to UTC conversions, computed by a call to the MIRAGE program ETTAIV ([Moyer, 2000]), based at DSS24. We compute TDB to UTC at DSS24 because our DSS24 based TTSDTE measurements (described below) provide our only precise measure of absolute time.
LGRS®TDB:
Combining information from TC41A, TC51A, and TC61A, adjusted for lighttravel time (LTM1A), CLK1A coordinates LGRS and TDB, but does not account for hardware/software/transmitter delays or GRAILA relativistic time dilation.
For star tracker data, we need the conversion:
BTC (SCLK)®TDB: Combine information for TC31A, TC51A, and TC61A, adjusted for lighttravel time (LTM1A), not accounting for delays or time dilation.
An additional source of timing information, not originally included in the GRAIL mission, became available. A dedicated DSN station occasionally eavesdrops on Sband range code Time Transfer System (TTS) communications between GRAILA and GRAILB, directly correlating LGRS®UTC, without the unknown latency biases of our usual telemetry data. After corrections for tropospheric delay, ionospheric delay, DSN clock drift, and DSN radio science receiver hardware delays [Esterhuizen, 2012], we save these TTS DirecttoEarth (TTSDTE) observations in TDE00 products, and calibrate CLK1A latency biases based on TTSDTE.
For the primary mission and the extended mission, we chose constant offsets, “Bias time”, which approximate the UTCLGRS bias. We tag many Level 1A products with LGRS + Bias time. For the main mission, LGRS + Bias time differs from UTC by about 20 seconds; for the extended mission, LGRS + Bias time differs from UTC by a fraction of a second.
Mission 
Bias Time (s) 
Primary 
382581795 
Extended 
398306639 
Four sources produce Level 0 data:
1) Deep Space Network (DSN)
2) Radio Science Receivers (RSR)
3) GRAIL spacecraft (GRAILA/GRAILB): Telemetry Delivery Subsystem (TDS)
4) Distributed Object Manager (DOM) – database which contains GRAIL spacecraft teamcomputed spacecraft mass and center of mass estimates. Data available after propulsive events.
The DSN provides media calibrations ([Runge, 2008]), earth orientation parameters ([Gross, 2009]), and closedloop twoway Doppler and oneway Doppler data, in TRK218 Orbit Data Files (ODF) ([Stipanuk, 2000]). RSR, located at DSN sites, produce openloop oneway Doppler in Tracking Data Message (TDM) format ([CCSDS, 2007]), which we convert to TRK218 files. A dedicated RSR at DSS24 produces TTSDTE timing data.
Data from DSN, RSR, and DOM enter the Planetary Data System. Lockheed Martin (LMA) and the Jet Propulsion Laboratory (JPL) maintain a Ground Data System (GDS), which receives SFDU telemetry packets from GRAILA and GRAILB, processes them, and stores raw packets and processed results.
This section introduces the GDS.
GRAIL transmits data to Earth in Standard Formatted Data Units (SFDUs). For a description of the SFDU format, see [Clark, 2008], [Wilson, 2003], and [Wilson, 2012]. Each SFDU consists of a header followed by a Compressed Header Data Object (CHDO) containing data.
The GDS timestamps each SFDU packet received, and additionally notes the Application Packet ID (APID) contained in the SFDU header, which identifies the type of data contained in a packet. The GDS extracts data from SFDUs, and stores extracted data in a channelized database. The GDS also stores raw SFDUs.
GDS users can query for raw packets with a particular APID within a specified timespan, or submit channelized queries for specific types of data extracted from SFDUs. Channelized queries produce a Level 0 text file with a header and column oriented, space separated data. For a description of Level 0 files, see the GRAIL “Data Product Software Interface Specification” ([Kahan DPSIS, 2012]).
For most APIDs, channelized queries suffice, so a user does not need to revisit raw data. For channelized query data formats, see [Kahan DPSIS, 2012].
The GDS does not channelize three data types:
APID 
Type 
Description 
72 
Engineering 
Log messages and LGRS health status (Kaband SNR, Sband SNR, relative time tag offset) 
73 
Science 
Science data from GRAIL mission, not channelized 
01 
TC 
Time correlation between TDS and Earth 
Engineering SFDU CHDOs contain nothing but Blackjack data, saved under the Blackjack datalink protocol, described in [Farrington, 2001]. Blackjack identifies each packet type with a fourletter packet ID, described in [Rogstad, 2012]. To read engineering SFDUs, concatenate data CHDOs from packet to packet into a single file. Decompose as a sequence of Blackjack packets, undoing the Blackjack datalink protocol by “deescaping” and applying a CRC errorchecking integrity check [Farrington, 2001].
Science SFDU data CHDOs contain:
1) Sixbyte representation of Base Clock Time (BTC). First four bytes contain a bigendian unsigned integer number of seconds, last two bytes contain a bigendian unsigned integer number of subseconds, where one subsecond = 1/65536 seconds.
2) Eightbyte representation of Real Clock Time (RTC) as a bigendian double precision number, swapping first four bytes for last four bytes.
3) A sequence of Blackjack packets.
Blackjack packet IDs include:
Packet ID 
Associated Product 
Description 
QFIT 
SBR KBR 
Sband TTS phase and range Kaband phase 
LOGM 
ILG 
LGRS log messages 
ADCP 
IHK 
LGRS voltage, temperature, current 
MEOK 
IHS 
LGRS health status 
TIME 
SNV 
Ancillary information for TTS 
PPST 
PPS 
Pulse per second time 
QFIT packets with PseudoRandom Noise (PRN) identifier artificially set = 1 or 2 contain SBR data. QFIT packets claiming PRN = 50 contain KBR data.
Time correlation SFDU data CHDOs contain, in order:
Contents 
Length (bytes) 
Threshold 
4 
Rate 
4 
Encoding 
4 
Counter 
4 
RTC upper 
4 
RTC lower 
4 
Subseconds 
2 
Pad 
2 
We only use Threshold, RTC upper, RTC lower, and subseconds. Threshold contains a fourbyte bigendian integer data length in bytes. RTC upper + RTC lower provide a bigendian long integer number of SCLK (BTC) seconds, swapping the first four bytes for the last four bytes. Subseconds contains a bigendian short integer number of SCLK subseconds, where one subsecond = 1/65536 seconds.
The USO1A product tabulates GRAIL USO frequency estimates. To produce USO1A:
1) Openloop RSR receivers at each DSN complex record real and imaginary (I and Q) components of signal from the GRAILA/B Xband Radio Science Beacon (RSB). We solve for Xband oneway frequency (F1) relative to a predicted frequency at each epoch based on spectral analysis of local Fast Fourier Transforms (FFTs) ([Paik, 2011]). Add predicted frequency and relative frequency to produce the sky frequency of the Xband signal.
2) DSN stations produce Sband closedloop twoway Doppler (F2).
3) JPL runs MIRAGE software (see the orbit determination and propagation section, below) to solve for GRAIL orbits, dividing data into sub2day arcs, processing, among other data types, F1 and F2. The orbit solution for each arc solves for “local” variables such as solar pressure scaling, empirical periodic accelerations, KBRR (intersatellite rangerate) bias and drift, time tag offsets for KBRR observables, and F1 frequency corrections which stochastically update about once per orbit. During MIRAGE orbit determination, we apply media corrections to DSN data and relativistic corrections based on a degree40 order40 spherical harmonic lunar gravity potential and a truncation of the GGM02C terrestrial gravity field solution.
4) MIRAGE produces Sband F1 frequency corrections. Multiply by 11/3 to convert to Xband, and subtract Xband frequency corrections from expected frequency. Each estimate holds for a 2hour time span; assign a timetag corresponding to the midpoint of that time span, and list results in USO1A frequency records:
Reception time (UTC) 
Xband frequency 
Since USO frequency, Xband frequency, and Kaband frequency differ by constant multiples, we can trivially compute USO and Ka frequencies.
GRAIL 
From 
To 
Frequency Multiple 
A 
USO 
Ka 
12×564 
A 
USO 
X 
8×2×(109+83385/262144) 
B 
USO 
Ka 
12×564 
B 
USO 
X 
8×2×(109+83476/262144) 
For more details, see [Duncan, 2012].
GRAIL Sband processing employs carrierphase smoothing, a standard processing technique. Suppose we have a sequence of continuously tracked phase measurements
and a sequence of cotemporal range measurements
Range measurements suffer from greater noise, while phase measurements carry an unknown overall bias pertrack. Carrierphase smoothing produces rangelike phase measurements by removing the mean difference between range and phase measurements:
An Sband Time Transfer System (TTS) correlates GRAILA and GRAILB LGRS clocks. Each satellite includes a phasemodulated Sband receiver/transmitter, driven by the LGRS UltraStable Oscillator (USO). We difference GRAILA®GRAILB and GRAILB®GRAILA carrier phase smoothed Sband pseudorange measurements to eliminate range contributions, leaving intersatellite LGRS clock offset. We list GRAILAGRAILB intersatellite clock offset by TDB time in DEL1A.
To allow carrierphase smoothing, GRAILA and GRAILB track carrier phase for TTS signals. GRAILA Sband runs at ~2.032 GHz, while GRAILB Sband runs at ~2.207 GHz. Each satellite beats the incoming TTS signal against a virtual local oscillator, at beat frequencies of ~3.608 MHz for GRAILA®GRAILB and ~2.858 MHz for GRAILB®GRAILA. For details, see [Duncan, 2012]. To preserve precision, the onboard computers count carrier phase cycles modulo 10^{10}, hence carrier phase counts wrap about once an hour, and a user must unwrap carrier phase cycle counts.
Our algorithm to produce DEL1A starts from a pair of SBR1A files containing Sband range and phase measurements, a PLT1A file listing interspacecraft light time estimates, and a CLK1B file correlating LGRS and TDB time. To create DEL1A:
1) Divide Sband data into time intervals with continuously tracked phase data – no lossoflock.
2) In each time interval, process carrierphase smoothed Sband range, occasionally adjusted for range jumps due to code cycle integer changes at GPA resets on each spacecraft. Letting c represent the speed of light, GRAILA transmission code cycle length in seconds equals:
For GRAILB, transmission code cycle length equals:
Modify GRAILB received range by GRAILA transmission code cycles as needed; modify GRAILA received range by GRAILB transmission code cycles as needed. We observe the following code cycle jumps:
Receiving Satellite 
Start Time (LGRS + Bias time) 
End Time (LGRS + Bias time) 
Code Cycles Added To Measurement 
A 
387437849 
387979259 
+1× 
B 
387979259 
390985361 
+1× 
B 
399616122 
399641504 
+1× 
3) By an iterative algorithm, described in [Kruizinga, 2012], determine the LGRS clock offset of GRAILA with respect to GRAILB. In DEL1A, we call this offset “eps_time.”
4) Since we carrierphase smooth eps_time on each time interval independently, solutions don’t match at interval end points. We calculate a series of biases to force nearcontinuity at interval end points, and apply these biases to create a continuous version of eps_time, which we save in DEL1A as “eps_drift.”
It should be noted that we see systematic deviations between Sphase and Srange
measurements, and between Kaphase and Sphase measurements. We expect these
deviations to affect DEL1A at a level below timing system requirements. These
deviations are related to temperature, and may also be related to Sband SNR peaks
which show up as GRAILB flies over the north or south pole.
From timetotime, TTS Sband data is not available. During these outages, we
fill in DEL1A based on Kaband data, releveled to match Sband.
Having calculated clock offset, we also calculate the approximate range t from GRAILA to GRAILB based on Sband data (t includes an ~130m arcdependent bias), and save t in an SBR1B file. For details, see “Timing of Science Data for the GRAIL Mission” ([Kruizinga, 2012]). Although our processing doesn’t use SBR1B directly, we do use SBR1B to look for glitches in KBR1C, our primary science product. (The level 1B KBR product is designated as ‘1C’ to distinguish it from earlier versions of KBR1B which did not contain the final four columns of information on temperature range corrections.)
Our final timing products, CLK1B and USO1B, account for absolute and relative timing, processing measurements through a noncausal Kalman filter. Inputs:
1) CLK1A_A/B. As described above, CLK1A relates LGRS to TDB time. Transmitter and receiver delays add an unknown bias to CLK1A measurements.
2) Twice every two weeks, DSS24 eavesdrops on the TTS system linking GRAILA and GRAILB, and creates TTSDTE measurements, from which we create the TDE00 product. Our clock Kalman filter processes TTSDTE measurements, and we calibrate CLK1A biases based on TTSDTE.
3) DEL1A relates GRAILA and GRAILB clocks.
4) REL1A_A/B products list relativistic time dilation effects.
5) USO1A_A/B. As described above, USO1A tabulates USO frequency estimates.
Algorithm:
1) Remove relativistic effects from CLK1A, TTSDTE, DEL1A, and USO1A.
2) Process all four measurement types through a noncausal Kalman filter, solving for bias, rate, and drift on each clock. Allow bias, rate, and drift parameters to random walk.
3) Delete outliers. In particular, we exclude some bad early CLK1A data.
4) Add white noise bias updates at clock resets. Allow drift parameters to random walk more quickly near solar events. Recalibrate CLK1A biases based primarily on TTSDTE.
5) Rerun filter.
6) Add relativistic effects back in, and tabulate estimated LGRS to TDB clock corrections in CLK1B_A/B and LGRS to TDB frequency corrections in USO1B_A/B.
Since CLK1B becomes an input to production of CLK1A, DEL1A, and USO1A, as well as our orbit estimates and relativistic corrections, LGRS clock estimation requires an iterative loop, in which we reestimate orbits, CLK1B inputs, CLK1B, orbits, CLK1B inputs, CLK1B…
Results seem consistent with microsecond level accuracy.
We frequently resample data to new time tags applying Lagrange interpolation, particularly when shifting observation time tags from one clock to another. To reduce floatingpoint error, our Lagrange interpolation algorithm processes time as an integer plus a double, by a computation borrowed from GRACE software ([Wu, 2006]). As a typical example, when processing Kaband data in KBR_order:
1) An input KBR1A data file contains Kaband observations, tagged at LGRS + Bias time, saved as:
= integer number of seconds
= integer number of microseconds
2) For each time tag, we read a CLK1B correction from LGRS + Bias time to TDB time.
3) Express each time tag in TDB as an integer number of seconds plus an integer number of microseconds plus a double precision correction:
4) Perform secondorder Lagrange interpolation to convert observations to evenly spaced TDB epochs.
In these equations, = output time in TDB, = input data times in TDB, = Lagrange interpolation coefficients, and = Kaband observations before and after resampling. When differencing times, difference integer second, integer microsecond, and double precision correction components separately.
GRAILA and GRAILB record and transmit 10 Hz Kaband measurements, which we typically process at 0.5 Hz to reduce computational time. Simple data decimation would alias signals above the Nyquist frequency (in this case, 0.25 Hz) into the Nyquist band, so we apply a CRN filter to reduce aliasing.
As described in [Thomas, 1999], a CRN filter applies the Cfold convolution of a rectangular window to a lowpass filter to create a finitetime approximation of a lowpass filter:
1) Given an input 10 Hz data sequence x, our CRN filter produces an output 0.5 Hz data sequence y. A perfect filter would provide unit response to frequencies below 0.25 Hz, and no response to frequencies above 0.25 Hz. We tune our CRN filter to minimize gain ripple (deviation from unit response at low frequencies) and aliasing (response at high frequencies, which mimics a low frequency response in a discrete filter). Let G(f) denote filter response at frequency f. For this document, we define gain ripple, normalized for unit response at twiceperrev frequency = 0.28 mHz, as:
For any integer n, at any target frequency f, with a 0.5 Hz sampling rate we cannot distinguish between a signal at f and a signal at f+0.5n. We define aliasing as the squarerootofthesumofsquares of normalized gain at high frequency signals differing from our target low frequency signal by a multiple of 0.5 Hz:
2) Specify a convolution order C and a filter length N (number of input points in convolved window). N must be an odd number. Based on our target bandwidth B = 0.25 Hz, compute the number of frequency points spanned by B within the Discrete Fourier Transform (DFT) of our finite filter:
3) Compute unnormalized frequency response for a lowpassed Cfold rectangular window convolution.
4) Construct filter as a DFT of frequency response.
5) Normalize filter coefficients F(j) for unit response at 0.28 mHz (twiceperrev frequency on GRAIL).
6) Apply filter to x and generate y.
As filter length increases, ripple and aliasing decrease. On the other hand, longer filter lengths increase computational time and extend corruption from bad data.
As convolution number increases, ripple and aliasing improve at low frequencies, and deteriorate at high frequencies.
At certain lengths and convolution numbers, gain ripple looks particularly good at low frequencies. Since the GRAIL mission places particular importance on global – i.e. low frequency – signals, we looked for CRN parameters with unusually good low frequency ripple. Complete search criteria:
1) Input data rate = 10 Hz.
2) Bandwidth = 0.25 Hz.
3) CRN convolution number 7, 9, or 11.
4) Filter length between 600 and 900 (60 to 90 seconds).
5) For frequencies below 0.15 Hz, ripple and aliasing below 1e6.
6) Reduced low frequency ripple.
We found two acceptable candidates: CRN9747 and CRN11825. We chose the shorter candidate: CRN9747.
We plot gain and ripple sidebyside with the GRACE 0.1 Hz bandwidth filter, CRN7707 (chosen under similar criteria). Gain and ripple improve by several orders of magnitude at all frequencies.
We define related CRN filters for rate:
and acceleration:
normalizing both for unit response at 0.28 mHz.
GRAIL estimates a lunar gravity field based on the relative motion of GRAILA and GRAILB. We process data from an intersatellite Kaband system to estimate crafttocraft separation. Each satellite includes a Kaband receiver/transmitter, and GRAILA tracks Ka from GRAILB while GRAILB tracks Ka from GRAILA.
GRAILA and GRAILB track carrier phase for Kaband signals. GRAILA Kaband runs at ~32.703 MHz, while GRAILB Kaband runs at ~32.704 MHz; Kafrequency differs by 670.032 KHz. Each satellite beats the incoming signal against the local signal. For details, see [Duncan, 2012]. To preserve precision, the onboard computers count carrier phase cycles modulo 10^{8}, hence carrier phase counts wrap once every 149 seconds, and a user must unwrap carrier phase cycle counts.
GRAIL Kaband processing borrows extensively from procedures developed for GRACE ([Wu, 2006]). KBR1A files contain raw Kaband carrier phase measurements. A threepart process generates biased dual oneway range in KBR1C:
1) KBR_debreak. Identify and flag phase breaks.
2) KBR_order. Resample to convert time tags from LGRS to TDB.
3) KBR_compress. Fill in short gaps. Compute biased dual oneway range. Apply CRN filter to produce range, rangerate, and rangeacceleration. Compute lighttime correction from PLT1A inputs, and apply CRN filter. Compute correction for phase center offset from center of mass from PCI1A inputs. Set flags.
The seventh column of KBR1A flags data quality:
Bit 0 = Undefined, set in KBR_debreak to identify possible phase break
Bit 1 = Phase break occurred in Ka
Bit 3 = Cycle slip detected in Ka
Bit 4 = Insane Ka polynomial coefficient
Bit 7 = Ka SNR < 450
The sixteenth column of KBR1C flags data quality:
Bit 0 = Phase break
Bit 1 = Unreliable PCI data for antenna center correction
Bit 2 = Interpolated PCI data for antenna center correction
Bit 3 = Extrapolated clock correction > 5s from fit center
Bit 4 = Extrapolated clock correction < 5s from fit center
Bit 5 = Data corrected for time tag bias of Ka phase
Bit 6 = Filled data > 5s from fit center
Bit 7 = Filled data < 5s from fit center
For more detail on some of these flags, read our descriptions of KBR_debreak, KBR_order, and KBR_compress, below.
KBR_debreak reads a KBR1A file, identifies phase breaks based on data gaps, and flags the input KBR1A file.
1) Rather than send Ka phase data, the GRAIL onboard computers compute a bestfit polynomial to 10 seconds of phase data, and GRAIL sends polynomial coefficients plus residuals to the ground. We exclude data points with anomalous Kaphase polynomial coefficients (which we call “insane” coefficients).
2) If a data gap exceeds 21 seconds, set Bit 1 = 1 to identify a phase break at the first postgap point.
3) If a data gap exists, but does not exceed 21 seconds, set Bit 0 = 1 to identify a possible phase break at the first postgap point.
KBR_order reads a KBR1A file and a CLK1B file and resamples to convert time tags from LGRS to TDB.
1) Linearly interpolate CLK1B corrections and apply to KBR1A.
2) Resample data, as described above for this particular example. Do not resample through phase breaks or possible phase breaks.
3) Save results in KBR1A format, but with TDB time tags.
KBR_compress reads a pair of KBR1A format files with TDB time tags, a pair of PCI1A files, a pair of USO1B files, and a pair of PLT1A files, and produces a KBR1C DOWR file.
1) Read GRAILA and GRAILB Kaband carrier frequencies f_{A} and f_{B} from the input USO1B files.
2) Read GRAILA and GRAILB Kaband phase measurements and from the input KBR1A format files with TDB time tags.
3) Compute biased Dual One Way Range (DOWR):
where c = speed of light. Note that a Kaband phase break introduces a new bias on DOWR. Since we process data in arcs, day boundaries also introduce a new bias.
4) For short data gaps with at least three DOWR points available on each side, fill the gap by least squares cubic interpolation through up to 100 points on each side; with fewer than three points on either side, interpolate linearly. For data gaps longer than 21 seconds, do not interpolate.
5) Read = light time from GRAILA to GRAILB, = light time from GRAILB to GRAILA, and mooncentered solar system barycentric positions from the input PLT1A files. Compute = instantaneous Euclidean intersatellite range. At time t, generate by seventhorder Lagrange interpolation. Calculate the Time of Flight Correction (TOF) to DOWR:
6) CRN filter TOF.
7) Read antenna phase center range, rangerate, and rangeacceleration corrections for both spacecraft from the input PCI1A files. These corrections have already been CRNfiltered. Add corrections to compute an overall Antenna Phase Center Correction for DOWR, with an associated rate and acceleration.
8) Write a KBR1C file with range, rangerate, rangeacceleration, temperature, temperaturerate, temperatureacceleration, and TOF and APCC corrections.
Flag KBR1C data:
1) From PCI1A, flag when from raw data for Ka boresight calibration slew
For each spacecraft, an onboard Kalman filter generates quaternions that represent spacecraft attitude in SSB coordinates (for a discussion of quaternions in a similar context, see [Wu, 2006]), processing starcamera and high rate gyro data (not transmitted to the ground). We save these quaternions, indexed by BTC time, in an SCA1A file. Before orbit determination, we convert quaternions to an interspacecraft range correction by a twopart process:
1) SCA1A2SCA1B. Resample from BTC to TDB.
2) SCA2PCI. Compute antennaphasecenteroffsetinterspacecraftrange correction.
From an input SCA1A file, SCA1A2SCA1B produces an output SCA1B file by:
1) Flip overall signs of quaternions in SCA1A as necessary to maintain continuity of coefficients.
2) Invoke the attitude interpolation algorithm “slerp” ([Shoemaker, 1985]) to fill any gaps ≤ 30 seconds.
3) Convert BTC time tags to TDB.
4) Call slerp to resample data to even TDB intervals.
From an input SCA1B file, SCA2PCI produces an output PCI1A file containing range, rangerate, and rangeacceleration corrections for antenna phase center (PC) offset from spacecraft center of mass. For simplicity, we consider GRAILA:
1) Read the GRAILA antenna PC offset vector (VKB) in the bodyfixedframe as a function of TDB time from the VKB1B product.
2) At each epoch, calculate the lineofsight vector from GRAILA to GRAILB in SSB coordinates.
3) Rotate VKB by each GRAILA SCA1B quaternion. Offset vectors now lie in SSB coordinates.
4) Calculate a range correction at each epoch by projection of VKB onto lineofsight.
5) CRN filter range corrections to produce range, rangerate, and rangeacceleration corrections.
GRAIL Level1 processing rests on orbit determination for GRAILA and GRAILB. Since orbit determination informs product derivation, and product derivation in turn improves orbit determination, we iterate between orbit determination and product derivation to produce a final solution.
JPL orbit estimation runs a software package called Multiple Interferometric Ranging Analysis Using GPS Ensemble (MIRAGE). JPL developed MIRAGE as part of the GPS Data Processing Facility (GDPF) created in support of the TOPEX/Poseidon mission, and later added capabilities to accommodate the GRACE mission. MIRAGE inherits from the JPL legacy Orbit Determination Program (ODP) ([Moyer, 2000]).
Our MIRAGE runs operate in mooncentered solar system barycentric frame, with TDB time.
MIRAGE orbit determination involves:
1) Trajectory propagation
2) Observation processing
3) Least squares filtering
MIRAGE least squares filtering follows Bierman’s Square Root Information Filter (SRIF) algorithm ([Bierman, 1977]).
This section describes orbit determination data inputs and our dynamic and kinematic models. Content borrows heavily from the GRACE Level2 standards document ([Watkins, 2012]), and a published article on GRAIL methodology ([Park, 2012]).
JPL orbit determination runs assimilate:
1) Openloop oneway Xband Radio Science Beacon (RSB) purephase transmissions (F1)
2) Closedloop twoway Sband Doppler, produced at DSN stations (F2)
3) KBR1C intersatellite rangerate (KBRR) measurements
4) A Small Forces File (SFF), that lists thrust event times, magnitudes, and durations
5) SCA1B spacecraft attitude quaternions
Other inputs to MIRAGE include Earth Orientation Parameters, planetary ephemerides, solar radiation pressure (SRP) parameters, and DSNtracking media calibrations.
We typically process data in 36hour arcs, running from 18:00 UTC to 06:00 UTC. Successive 36hour arcs overlap by 12 hours, and we gauge orbit quality based on the consistency of overlapping orbits. Large thruster events, however, force an arc break and thwart overlap evaluation.
Our orbit determination runs estimate the following parameters:
Parameter name 
Parameter description 
for each in series 
Max count 
For each of GRAILA and GRAILB 

X 
X, Y, Z components of position in mooncentered solar system barycentric frame 
n/a (at initial epoch) 
1 
Y 
1 

Z 
1 

DX 
X, Y, Z components of velocity w.r.t. mooncentered solar system barycentric frame 
n/a (at initial epoch) 
1 
DY 
1 

DZ 
1 

SOLCOF 
SRP scale factor in direction of suntospacecraft vector 
whole arc 
1 
GX 
SRP scaling factors in two directions orthogonal to suntospacecraft vector 
1 

GY 
1 

D1T 
spacecraft to Earth F1 bias 
whole arc for first two MIRAGE passes; 6813 s = 1.8925 hr, 6630 s = 1.8417 hr during PM & XM, respectively, for later passes 
50 
D2T 
spacecraft to Earth F1 drift 
50 

D3T 
spacecraft to Earth F1 drift rate 
50 

PR 
empirical accel., radial, constant (bias) 
6813 s = 1.8925 hr, 6630 s = 1.8417 hr during PM & XM, respectively (i.e. approximately matching 1 orbital period) 
24 
CR1 
empirical accel., radial, amp. of 
24 

SR1 
empirical accel., radial, amp. of 
24 

CR2 
empirical accel., radial, amp. of 
24 

SR2 
empirical accel., radial, amp. of 
24 

PT 
empirical accel., transverse, constant (bias) 
24 

CT1 
empirical accel., transverse, amp. of 
24 

ST1 
empirical accel., transverse, amp. of 
24 

CT2 
empirical accel., transverse, amp. of 
24 

ST2 
empirical accel., transverse, amp. of 
24 

P 
empirical accel., normal, constant (bias) 
24 

CN1 
empirical accel., normal, amp. of 
24 

SN1 
empirical accel., normal, amp. of 
24 

CN2 
empirical accel., normal, amp. of 
24 

SN2 
empirical accel., normal, amp. of 
24 

IDLX0 
X, Y, Z components of small force (thruster) event accelerations, w.r.t. mooncentered solar system barycentric frame 
n/a (with each event as detailed in SFF files) 
4 
IDLY0 
4 

IDLZ0 
4 

Total for each spacecraft 
531 

For both GRAILA and GRAILB together (interspacecraft tracking) 

RRB 
KBRR bias 
whole arc 
50 
RRD 
KBRR drift 
50 

RRCP1 
KBRR periodic term, amp. of 
50 

RRSP1 
KBRR periodic term, amp. of 
50 

TIMETAG_252_1 
KBRR time tag bias until 1^{st} relative clock offset resetting event (such as GPA reboot) 
varies (nominally = whole arc) 
1 
TIMETAG_252_2 
ditto, between 1^{st} and 2^{nd} such events 
varies (nominally not used = 0 s) 
1 
TIMETAG_252_3 
ditto, after 2^{nd} such event 
1 

Total for interspacecraft tracking 
203 

OVERALL TOTAL 
1265 
Note that most of these parameters constrain the effect of (“soak up”) unmodeled nongravitational accelerations and measurement biases. JPL chose these parameters based on experience with GRACE, the Mars Reconnaissance Orbiter, Magellan, and other missions.
IV. Momentum Wheels
Attitude control on each spacecraft depends on four momentum wheels. During orbit determination, when wheel 3 on GRAILB (and to a lesser extent, at least one other wheel) rotates at ~628 radians/s, we frequently see outliers.
WRS1A and WRS1B products list wheel speeds in radians/s.
ACS Attitude control system
APCC Antenna phase center correction
APID Application packet ID
ATBD Algorithm Theoretical Basis Document
BTC Base time clock
CD&H Command and data handling satellite computer
CHDO Compressed header data object
CMIC CD&H module interface card
CRC Cyclic redundancy check
CRN N recursive selfconvolutions of a rectangular window
CSP Control statement processor
DFT Discrete Fourier Transform
DLP Data link protocol
DOM Distributed object manager
DOWR Dual oneway range
DPSIS Data Product Software Interface Specification
DSN Deep space network
DSS Deep space station
ETTAIV Ephemeris time – international atomic time (vector)
F1 Oneway doppler
F2 Twoway doppler
FFT Fast Fourier Transform
GDPF GPS data processing facility
GDS Ground data system
GPA Gravity recovery processor assembly
GPS Global Positioning System
GRACE Gravity Recovery and Climate Experiment
GRAIL Gravity Recovery and Interior Laboratory
GRAILA GRAIL satellite A
GRAILB GRAIL satellite B
IMU Inertial measurement unit
JPL Jet Propulsion Laboratory
KBR Kaband range
KBRR Kaband range rate
LGA Lowgain antenna
LMA Lockheed Martin
LGRS Lunar gravity ranging system
MIRAGE Multiple Interferometric Ranging Analysis Using GPS Ensemble
MWA Microwave assembly
ODF Orbit data file
ODP Orbit Determination Program
PPS Pulsepersecond signal
PRN Pseudorandom noise
RSB Radio science beacon
RSBA Radio science beacon assembly
RSR Radio science receiver
RTC Real time clock
SCET Spacecraft event time
SCLK Spacecraft clock
SDS Science data system
SFDU Standard formatted data unit
SNR Signaltonoise ratio
SRIF Squareroot information filter
SRP Solar radiation pressure
TC Time correlation packet
TDB Barycentric dynamical time
TDM Tracking data message
TDS Telemetry delivery subsystem
TOF Time of flight correction
TTA Time transfer assembly
TTS Time transfer system
TTSDTE Time transfer system directtoearth
USO Ultra stable oscillator
UTC Coordinated universal time
[Asmar, 2012] Sami W. Asmar, Alexander S. Konopliv, S. Ryan Park, Michael M. Watkins, Gerhard Kruizinga, Maria T. Zuber, David E. Smith, James G. Williams, Meegyeong Paik, DahNing Yuan, Eugene Fahnestock, Dmitry Strekalov, Nate Harvey, Daniel S. Kahan, and Wenwen Lu, “The Scientific Measurement System of the Gravity Recovery and Interior Laboratory (GRAIL) Mission.” Submitted to Space Science Reviews, 2012.
[Duncan, 2012] Courtney Duncan, “Once Upon a Fortnight." JPL D61501, 2012.
[Fahnestock, 2012] Eugene G. Fahnestock, Ryan S. Park, DahNing Yuan, and Alex S. Konopliv, “Spacecraft Thermal and Optical Modeling Impacts on Estimation of the GRAIL Lunar Gravity Field.” Proc. AIAA/AAS Astrodynamics Specialist Conference, Minneapolis, Minnesota, August 1316, 2012, paper AIAA 20124428.
[Kahan handbook, 2012] Daniel Kahan and Nate Harvey, “GRAIL Level 1 Data Product User Handbook.” 2012.
[Kahan DPSIS, 2012] Daniel Kahan, “GRAIL Data Product Software Interface Specification.” 2012.
[Kruizinga, 2012] Gerard L. H. Kruizinga and William I. Bertiger, “Timing of Science Data for the GRAIL Mission.” 2012.
[Park, 2012] Ryan S. Park, Sami W. Asmar, Eugene G. Fahnestock, Alex S. Konopliv, WenWen Lu, and Mike M. Watkins, “Gravity Recovery and Interior Laboratory Simulation of Static and Temporal Gravity Field.” Journal of Spacecraft and Rockets, Vol. 49, No. 2, MarchApril 2012.
[Rogstad, 2012] M. Girard, T. Rogstad, “GRAIL Telemetry Dictionary (TD).” JPLD71987, Revision E, 2012.
[Turyshev frames, 2012] Slava G. Turyshev, Olivier L. Minazzoli, and Viktor T. Toth, “Accelerating Relativistic Reference Frames in Minkowski SpaceTime.” 2012.
[Turyshev grail, 2012] Slava G. Turyshev, Viktor T. Toth, and Mikhail V. Sazhin, “General Relativistic Observables of the GRAIL Mission.” 2012.
[Bierman, 1977] G. Bierman, Factorization Methods for Discrete Sequential Estimation, Vol. 128, Academic Press, New York, 1977, pp. 13–32,135–161.
[CCSDS, 2007] Tracking Data Message: Recommended Standard, 2007. Published by CCSDS Secretariat, Space Communications and Navigation Office, 7L70, Space Operations Mission Directorate, NASA Headquarters, Washington DC 205460001, USA.
[Clark, 2008] Tracy Clark, “0161Telecomm: Telemetry Standard Formatted Data Unit (SFDU) Interface.” JPLD16765, Revision A, 2008.
[Esterhuizen, 2012] S. Esterhuizen, "MoontoEarth: Eavesdropping on the GRAIL interspacecraft timetransfer link using a large antenna and a software receiver," ION GNSS, ION, 2012.
[Farrington, 2001] Allen H. Farrington, “BlackJack Data Link Protocol.” JPLD20675, 2001.
[Gross, 2009] Richard Gross, “TRK221: DSN Tracking System Earth Orientation Parameter Data Interface.” JPL D16765, Revision B, 2009.
[Moyer, 2000] Theodore Moyer, Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation. Monograph 2, Deep Space Communications and Navigation Series, JPL, 2000.
[Paik, 2011] Meegyeong Paik, Sami W. Asmar, “Detecting High Dynamics Signals from Openloop Radio Science Investigations,” Proceedings of the IEEE, May 2011, Vol 99, Issue 5, pp. 881888.
[Runge, 2008] Thomas Runge, “TRK223: Media Calibration Interface.” JPLD16765, Revision C, 2008.
[Shoemaker, 1985] Ken Shoemaker, “Animating Rotation with Quaternion Curves.” SIGGRAPH ’85 Proceedings of the 12^{th} Annual Conference on Computer Graphics and Interactive Techniques, pp. 245254.
[Stipanuk, 2000] J. Stipanuk, “TRK218: Tracking System Interfaces, Orbit Data File Interface.” JPLD16765, Change 3, 2000.
[Thomas, 1999] J. B. Thomas, An Analysis of Gravity Field Estimation Based on Intersatellite Dual1Way Biased Ranging. JPL Publication 9815, 1999.
[Watkins, 2012] Michael M. Watkins and DahNing Yuan, GRACE JPL Level2 Processing Standards Document for Level2 Product Release 05. GRACE 327744 (v 5.0) March 17, 2012.
[Wilson, 2003] E. Wilson, “0171TelecommNJPL: JPL Created SFDU Structures.” JPLD16765, Revision A, 2003.
[Wilson, 2012] E. Wilson, “0172TelecommCHDO: DSN Created CHDO Structures.” JPLD16765, Revision E, 2012.
[Wu, 2006] SienChong Wu, Gerard
Kruizinga, and Willy Bertiger, Algorithm Theoretical Basis for GRACE
Level1B Data Processing V1.2. GRACE 327741 (JPL D27672), 2006.
APPENDIX A. GRAIL Timing Flow Diagram
APPENDIX B. Level 1A for GRAILA (INCLUDES ORBIT DETERMINATION AND PROPAGATION INPUT PREP)
APPENDIX C. Level 1A à Level 1B for GRAIL A and B
APPENDIX D. KBR1A à KBR1C and KBR CAL GRAIL A