Summary

I. Introduction

Gravity Recovery and Interior Laboratory (GRAIL) Level 1 data processing involves:

1) Telemetry pre-processing

2) Clock-timing correlation

3) Biased dual one-way inter-satellite range computation

4) Orbit determination

This document outlines Level 1 algorithms for telemetry pre-processing, clock-timing correlation, inter-satellite range computation, and orbit determination. For Level 1 software flow diagrams, see Appendices B, C, and D.

II. GRAIL Overview

The GRAIL mission processes data from two spacecraft, GRAIL-A and GRAIL-B (formally called Ebb and Flow), flying in formation orbits around the moon, with a payload that precisely measures the distance between GRAIL-A and GRAIL-B, 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 S-band communication with the Deep Space Network (DSN)

8) A Radio Science Beacon (RSB) with two antennas, which transmits one-way X-band carrier

9) S-band inter-spacecraft Time Transfer System (TTS)

10) Ka-band carrier phase tracking inter-spacecraft receiver/transmitter

III. Level 1 Processing

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 Pulse-Per-Second (PPS) signals that sync internal clocks and Time Correlation (TC) packets linking onboard and universal time.

2) Telemetry Pre-Processing – 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 Ultra-Stable-Oscillator (USO) frequencies based on X-band signals.

4) Relative Timing – An S-band Time Transfer System (TTS) measures GRAIL-A-GRAIL-B 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 non-causal Kalman filter, which produces an estimate of GRAIL-A LGRS and GRAIL-B LGRS deviations from TDB.

6) KBR Processing – A three-step process flags phase breaks (KBR_debreak), converts data from LGRS to TDB (KBR_order), and produces biased Dual One-Way Range (DOWR) measurements (KBR_compress). Separate sections describe time-tag conversion by Lagrangian data resampling and CRN filters to remove high-frequency DOWR noise.

Note that most steps in our Level 1 processing depend directly on the accuracy of:

7) Orbital Solutions – Assimilate one-way frequency data, two-way Doppler, and inter-satellite range-rate data in a Kalman filter, running JPL’s Multiple Interferometric Ranging Analysis Using GIPSY Ensemble (MIRAGE) software, solving for dynamic and kinematic spacecraft parameters.

GRAIL Absolute Timing

For ease of notation, this section discusses GRAIL-A timing in isolation. GRAIL-B timing runs the same algorithms in parallel.

I. Clocks

GRAIL-A absolute timing requires coordination of five clocks:

1) LGRS: Lunar Gravity Ranging System clock. Clock for on-board Ka-band (KBR), X-band (RSB), and S-band (TTS) instruments. Driven by an Ultra-Stable Oscillator (USO). Set to 0 when booted. Very stable clock.

2) BTC: Base Time Clock. On-board 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 solar-system 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.

II. Products

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 on-board clocks for relativistic time dilation. We define proper time as the time which a perfect clock on-board GRAIL-A would measure, including dilation. LGRS, BTC, and RTC reflect (GRAIL-A proper time + clock drift).

III. Time Coordination Packets

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 pulse-per-second (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. GRAIL-A 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, GRAIL-A 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 on-board 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 bit-rate, but can be more than a second.

c. Light-travel time from GRAIL-A 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 RTC-UTC.

IV. Time Coordination

Since GRAIL-A does not create our various time observables synchronously, we interpolate to form pseudo-synchronous 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 bit-rate dependent empirical correction C.

Date rate (kbps)	C GRAIL-A (s)	C GRAIL-B (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 beginning-of-packet 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 DSS-24. We compute TDB to UTC at DSS-24 because our DSS-24 based TTSDTE measurements (described below) provide our only precise measure of absolute time.

LGRS®TDB:

Combining information from TC41A, TC51A, and TC61A, adjusted for light-travel time (LTM1A), CLK1A coordinates LGRS and TDB, but does not account for hardware/software/transmitter delays or GRAIL-A relativistic time dilation.

For star tracker data, we need the conversion:

BTC (SCLK)®TDB: Combine information for TC31A, TC51A, and TC61A, adjusted for light-travel 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 S-band range code Time Transfer System (TTS) communications between GRAIL-A and GRAIL-B, 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 Direct-to-Earth (TTSDTE) observations in TDE00 products, and calibrate CLK1A latency biases based on TTSDTE.

V. Bias Time

For the primary mission and the extended mission, we chose constant offsets, “Bias time”, which approximate the UTC-LGRS 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

GRAIL Telemetry

I. Introduction

Four sources produce Level 0 data:

1) Deep Space Network (DSN)

2) Radio Science Receivers (RSR)

3) GRAIL spacecraft (GRAIL-A/GRAIL-B): Telemetry Delivery Subsystem (TDS)

4) Distributed Object Manager (DOM) – database which contains GRAIL spacecraft team-computed 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 closed-loop two-way Doppler and one-way Doppler data, in TRK-2-18 Orbit Data Files (ODF) ([Stipanuk, 2000]). RSR, located at DSN sites, produce open-loop one-way Doppler in Tracking Data Message (TDM) format ([CCSDS, 2007]), which we convert to TRK-2-18 files. A dedicated RSR at DSS-24 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 GRAIL-A and GRAIL-B, processes them, and stores raw packets and processed results.

This section introduces the GDS.

II. Ground Data System

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 time-span, 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]).

III. Types of GDS Data

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 (Ka-band SNR, S-band 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 four-letter 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 “de-escaping” and applying a CRC error-checking integrity check [Farrington, 2001].

Science SFDU data CHDOs contain:

1) Six-byte representation of Base Clock Time (BTC). First four bytes contain a big-endian unsigned integer number of seconds, last two bytes contain a big-endian unsigned integer number of subseconds, where one subsecond = 1/65536 seconds.

2) Eight-byte representation of Real Clock Time (RTC) as a big-endian 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	S-band TTS phase and range Ka-band 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 Pseudo-Random 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 sub-seconds. Threshold contains a four-byte big-endian integer data length in bytes. RTC upper + RTC lower provide a big-endian long integer number of SCLK (BTC) seconds, swapping the first four bytes for the last four bytes. Subseconds contains a big-endian short integer number of SCLK subseconds, where one subsecond = 1/65536 seconds.

GRAIL USO Frequency Determination

The USO1A product tabulates GRAIL USO frequency estimates. To produce USO1A:

1) Open-loop RSR receivers at each DSN complex record real and imaginary (I and Q) components of signal from the GRAIL-A/B X-band Radio Science Beacon (RSB). We solve for X-band one-way 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 X-band signal.

2) DSN stations produce S-band closed-loop two-way Doppler (F2).

3) JPL runs MIRAGE software (see the orbit determination and propagation section, below) to solve for GRAIL orbits, dividing data into sub-2-day 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 (inter-satellite range-rate) 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 degree-40 order-40 spherical harmonic lunar gravity potential and a truncation of the GGM02C terrestrial gravity field solution.

4) MIRAGE produces S-band F1 frequency corrections. Multiply by 11/3 to convert to X-band, and subtract X-band frequency corrections from expected frequency. Each estimate holds for a 2-hour time span; assign a time-tag corresponding to the mid-point of that time span, and list results in USO1A frequency records:

Reception time (UTC)

X-band frequency

Since USO frequency, X-band frequency, and Ka-band 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].

Carrier Phase Smoothing

GRAIL S-band processing employs carrier-phase smoothing, a standard processing technique. Suppose we have a sequence of continuously tracked phase measurements

and a sequence of co-temporal range measurements

Range measurements suffer from greater noise, while phase measurements carry an unknown overall bias per-track. Carrier-phase smoothing produces range-like phase measurements by removing the mean difference between range and phase measurements:

GRAIL Relative Timing

An S-band Time Transfer System (TTS) correlates GRAIL-A and GRAIL-B LGRS clocks. Each satellite includes a phase-modulated S-band receiver/transmitter, driven by the LGRS Ultra-Stable Oscillator (USO). We difference GRAIL-A®GRAIL-B and GRAIL-B®GRAIL-A carrier phase smoothed S-band pseudo-range measurements to eliminate range contributions, leaving inter-satellite LGRS clock offset. We list GRAIL-A-GRAIL-B inter-satellite clock offset by TDB time in DEL1A.

To allow carrier-phase smoothing, GRAIL-A and GRAIL-B track carrier phase for TTS signals. GRAIL-A S-band runs at ~2.032 GHz, while GRAIL-B S-band 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 GRAIL-A®GRAIL-B and ~2.858 MHz for GRAIL-B®GRAIL-A. For details, see [Duncan, 2012]. To preserve precision, the onboard computers count carrier phase cycles modulo 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 S-band range and phase measurements, a PLT1A file listing inter-spacecraft light time estimates, and a CLK1B file correlating LGRS and TDB time. To create DEL1A:

1) Divide S-band data into time intervals with continuously tracked phase data – no loss-of-lock.

2) In each time interval, process carrier-phase smoothed S-band 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, GRAIL-A transmission code cycle length in seconds equals:

For GRAIL-B, transmission code cycle length equals:

Modify GRAIL-B received range by GRAIL-A transmission code cycles as needed; modify GRAIL-A received range by GRAIL-B 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 GRAIL-A with respect to GRAIL-B. In DEL1A, we call this offset “eps_time.”

4) Since we carrier-phase smooth eps_time on each time interval independently, solutions don’t match at interval end points. We calculate a series of biases to force near-continuity 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 S-phase and S-range

measurements, and between Ka-phase and S-phase 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 S-band SNR peaks

which show up as GRAIL-B flies over the north or south pole.

From time-to-time, TTS S-band data is not available. During these outages, we

fill in DEL1A based on Ka-band data, re-leveled to match S-band.

Having calculated clock offset, we also calculate the approximate range t from GRAIL-A to GRAIL-B based on S-band data (t includes an ~130m arc-dependent 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.)

LGRS Clock Error Model

Our final timing products, CLK1B and USO1B, account for absolute and relative timing, processing measurements through a non-causal 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, DSS-24 eavesdrops on the TTS system linking GRAIL-A and GRAIL-B, 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 GRAIL-A and GRAIL-B 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 non-causal 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) Re-run 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 re-estimate orbits, CLK1B inputs, CLK1B, orbits, CLK1B inputs, CLK1B…

Results seem consistent with micro-second level accuracy.

Data Resampling

We frequently resample data to new time tags applying Lagrange interpolation, particularly when shifting observation time tags from one clock to another. To reduce floating-point 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 Ka-band data in KBR_order:

1) An input KBR1A data file contains Ka-band 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 second-order 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 = Ka-band observations before and after resampling. When differencing times, difference integer second, integer microsecond, and double precision correction components separately.

CRN Filter

GRAIL-A and GRAIL-B record and transmit 10 Hz Ka-band 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 C-fold convolution of a rectangular window to a low-pass filter to create a finite-time approximation of a low-pass 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 twice-per-rev 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 square-root-of-the-sum-of-squares 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 un-normalized frequency response for a low-passed C-fold 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 (twice-per-rev 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 1e-6.

6) Reduced low frequency ripple.

We found two acceptable candidates: CRN-9-747 and CRN-11-825. We chose the shorter candidate: CRN-9-747.

We plot gain and ripple side-by-side with the GRACE 0.1 Hz bandwidth filter, CRN-7-707 (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.

KBR Processing

I. Introduction

GRAIL estimates a lunar gravity field based on the relative motion of GRAIL-A and GRAIL-B. We process data from an inter-satellite Ka-band system to estimate craft-to-craft separation. Each satellite includes a Ka-band receiver/transmitter, and GRAIL-A tracks Ka from GRAIL-B while GRAIL-B tracks Ka from GRAIL-A.

GRAIL-A and GRAIL-B track carrier phase for Ka-band signals. GRAIL-A Ka-band runs at ~32.703 MHz, while GRAIL-B Ka-band runs at ~32.704 MHz; Ka-frequency 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⁸, hence carrier phase counts wrap once every 149 seconds, and a user must unwrap carrier phase cycle counts.

GRAIL Ka-band processing borrows extensively from procedures developed for GRACE ([Wu, 2006]). KBR1A files contain raw Ka-band carrier phase measurements. A three-part process generates biased dual one-way 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 one-way range. Apply CRN filter to produce range, range-rate, and range-acceleration. Compute light-time 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.

II. KBR_debreak

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 best-fit polynomial to 10 seconds of phase data, and GRAIL sends polynomial coefficients plus residuals to the ground. We exclude data points with anomalous Ka-phase 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 post-gap 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 post-gap point.

III. KBR_order

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.

IV. KBR_compress

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 GRAIL-A and GRAIL-B Ka-band carrier frequencies f_A and f_B from the input USO1B files.

2) Read GRAIL-A and GRAIL-B Ka-band 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 Ka-band 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 GRAIL-A to GRAIL-B, = light time from GRAIL-B to GRAIL-A, and moon-centered solar system barycentric positions from the input PLT1A files. Compute = instantaneous Euclidean inter-satellite range. At time t, generate by seventh-order Lagrange interpolation. Calculate the Time of Flight Correction (TOF) to DOWR:

6) CRN filter TOF.

7) Read antenna phase center range, range-rate, and range-acceleration corrections for both spacecraft from the input PCI1A files. These corrections have already been CRN-filtered. 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, range-rate, range-acceleration, temperature, temperature-rate, temperature-acceleration, and TOF and APCC corrections.

Flag KBR1C data:

1) From PCI1A, flag when from raw data for Ka boresight calibration slew

Quaternions

I. Introduction

For each spacecraft, an on-board Kalman filter generates quaternions that represent spacecraft attitude in SSB coordinates (for a discussion of quaternions in a similar context, see [Wu, 2006]), processing star-camera 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 inter-spacecraft range correction by a two-part process:

1) SCA1A2SCA1B. Resample from BTC to TDB.

2) SCA2PCI. Compute antenna-phase-center-offset-inter-spacecraft-range correction.

II. SCA1A2SCA1B

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.

III. SCA2PCI

From an input SCA1B file, SCA2PCI produces an output PCI1A file containing range, range-rate, and range-acceleration corrections for antenna phase center (PC) offset from spacecraft center of mass. For simplicity, we consider GRAIL-A:

1) Read the GRAIL-A antenna PC offset vector (VKB) in the body-fixed-frame as a function of TDB time from the VKB1B product.

2) At each epoch, calculate the line-of-sight vector from GRAIL-A to GRAIL-B in SSB coordinates.

3) Rotate VKB by each GRAIL-A SCA1B quaternion. Offset vectors now lie in SSB coordinates.

4) Calculate a range correction at each epoch by projection of VKB onto line-of-sight.

5) CRN filter range corrections to produce range, range-rate, and range-acceleration corrections.

Orbit Determination and Propagation

I. Introduction

GRAIL Level-1 processing rests on orbit determination for GRAIL-A and GRAIL-B. 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 moon-centered 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 Level-2 standards document ([Watkins, 2012]), and a published article on GRAIL methodology ([Park, 2012]).

II. Data Inputs

JPL orbit determination runs assimilate:

1) Open-loop one-way X-band Radio Science Beacon (RSB) pure-phase transmissions (F1)

2) Closed-loop two-way S-band Doppler, produced at DSN stations (F2)

3) KBR1C inter-satellite range-rate (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 DSN-tracking media calibrations.

We typically process data in 36-hour arcs, running from 18:00 UTC to 06:00 UTC. Successive 36-hour 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.

III. Estimated Parameters

Our orbit determination runs estimate the following parameters:

Parameter name	Parameter description	for each in series	Max count
For each of GRAIL-A and GRAIL-B
X	X, Y, Z components of position in moon-centered solar system barycentric frame	n/a (at initial epoch)	1
Y			1
Z			1
DX	X, Y, Z components of velocity w.r.t. moon-centered solar system barycentric frame	n/a (at initial epoch)	1
DY			1
DZ			1
SOLCOF	SRP scale factor in direction of sun-to-spacecraft vector	whole arc	1
GX	SRP scaling factors in two directions orthogonal to sun-to-spacecraft 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. moon-centered 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 GRAIL-A and GRAIL-B together (inter-spacecraft 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	varies (nominally not used = 0 s)	1
Total for inter-spacecraft tracking			203
OVERALL TOTAL			1265

Note that most of these parameters constrain the effect of (“soak up”) un-modeled non-gravitational 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 GRAIL-B (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.

List of Acronyms

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 self-convolutions of a rectangular window

CSP Control statement processor

DFT Discrete Fourier Transform

DLP Data link protocol

DOM Distributed object manager

DOWR Dual one-way range

DPSIS Data Product Software Interface Specification

DSN Deep space network

DSS Deep space station

ETTAIV Ephemeris time – international atomic time (vector)

F1 One-way doppler

F2 Two-way 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

GRAIL-A GRAIL satellite A

GRAIL-B GRAIL satellite B

IMU Inertial measurement unit

JPL Jet Propulsion Laboratory

KBR Ka-band range

KBRR Ka-band range rate

LGA Low-gain 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 Pulse-per-second signal

PRN Pseudo-random 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 Signal-to-noise ratio

SRIF Square-root 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 direct-to-earth

USO Ultra stable oscillator

UTC Coordinated universal time

Bibliography

I. GRAIL Project Documents

[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, Dah-Ning 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 D-61501, 2012.

[Fahnestock, 2012] Eugene G. Fahnestock, Ryan S. Park, Dah-Ning 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 13-16, 2012, paper AIAA 2012-4428.

[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, March-April 2012.

[Rogstad, 2012] M. Girard, T. Rogstad, “GRAIL Telemetry Dictionary (TD).” JPL-D-71987, Revision E, 2012.

[Turyshev frames, 2012] Slava G. Turyshev, Olivier L. Minazzoli, and Viktor T. Toth, “Accelerating Relativistic Reference Frames in Minkowski Space-Time.” 2012.

[Turyshev grail, 2012] Slava G. Turyshev, Viktor T. Toth, and Mikhail V. Sazhin, “General Relativistic Observables of the GRAIL Mission.” 2012.

II. Other (Non-Grail) Documents

[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 20546-0001, USA.

[Clark, 2008] Tracy Clark, “0161-Telecomm: Telemetry Standard Formatted Data Unit (SFDU) Interface.” JPL-D-16765, Revision A, 2008.

[Esterhuizen, 2012] S. Esterhuizen, "Moon-to-Earth: Eavesdropping on the GRAIL inter-spacecraft time-transfer link using a large antenna and a software receiver," ION GNSS, ION, 2012.

[Farrington, 2001] Allen H. Farrington, “BlackJack Data Link Protocol.” JPL-D-20675, 2001.

[Gross, 2009] Richard Gross, “TRK-2-21: DSN Tracking System Earth Orientation Parameter Data Interface.” JPL D-16765, 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 Open-loop Radio Science Investigations,” Proceedings of the IEEE, May 2011, Vol 99, Issue 5, pp. 881-888.

[Runge, 2008] Thomas Runge, “TRK-2-23: Media Calibration Interface.” JPL-D-16765, 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. 245-254.

[Stipanuk, 2000] J. Stipanuk, “TRK-2-18: Tracking System Interfaces, Orbit Data File Interface.” JPL-D-16765, Change 3, 2000.

[Thomas, 1999] J. B. Thomas, An Analysis of Gravity Field Estimation Based on Intersatellite Dual-1-Way Biased Ranging. JPL Publication 98-15, 1999.

[Watkins, 2012] Michael M. Watkins and Dah-Ning Yuan, GRACE JPL Level-2 Processing Standards Document for Level-2 Product Release 05. GRACE 327-744 (v 5.0) March 17, 2012.

[Wilson, 2003] E. Wilson, “0171-Telecomm-NJPL: JPL Created SFDU Structures.” JPL-D-16765, Revision A, 2003.

[Wilson, 2012] E. Wilson, “0172-Telecomm-CHDO: DSN Created CHDO Structures.” JPL-D-16765, Revision E, 2012.

[Wu, 2006] Sien-Chong Wu, Gerard Kruizinga, and Willy Bertiger, Algorithm Theoretical Basis for GRACE Level-1B Data Processing V1.2. GRACE 327-741 (JPL D-27672), 2006.

APPENDIX A. GRAIL Timing Flow Diagram

APPENDIX B. Level 1A for GRAIL-A (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