Revisiting GRACE Follow-On KBR Antenna Phase Center Calibration by Addressing Multipath Noise

Abstract

The Gravity Recovery and Climate Experiment Follow-On (GRFO) mission precisely measures the inter-satellite range between the centers of mass of its twin satellites to map the earth’s gravity field. The baseline ranging measurement is achieved using the K-band ranging (KBR) system, which is sensitive to satellite attitude variations caused by the offset between the satellite center of mass and the KBR antenna phase center. Accurate decoupling of the KBR range from attitude variations requires precise determination of the KBR’s antenna offset vectors (AOVs). To address this, GRFO conducted eight KBR calibration maneuvers on 17 and 28 September 2020. However, these maneuvers exaggerated the impact of microwave multipath noise, complicating AOV estimation. Existing studies have not fully mitigated this noise. This study introduces a new frequency-domain method to estimate AOVs by leveraging double-difference signals and analyzing their spectral characteristics, along with those of the KBR range during calibration maneuvers, to suppress multipath noise. Our recalibrated AOVs achieve good alignment between the KBR and laser ranging interferometer (LRI) ranging signals. We validate our recalibrated AOVs by comparing the residuals between the LRI and KBR ranging signals corrected using both recalibrated AOVs and documented AOVs. The results show that, for the majority (58.4%) of the analyzed period (from January 2020 to June 2023), the residuals corrected by the recalibrated AOVs are closer to the LRI ranging signal. These findings demonstrate the effectiveness of the proposed method in addressing multipath noise and improving the accuracy of KBR range measurements. This work provides a framework for future gravity missions requiring precise calibration of multipath effects in inter-satellite ranging systems.

1. Introduction

The Gravity Recovery and Climate Experiment Follow-On (GRFO), launched on 22 May 2018, continues the work of the Gravity Recovery and Climate Experiment (GRACE). Over its six years of operation so far, GRFO has demonstrated its effectiveness in earth science, providing critical insights into sea level change, ice melt dynamics, flood monitoring and earth’s free oscillations [1,2,3,4]. As the successor of GRACE, GRFO is equipped with the K-band ranging (KBR) system, star cameras (SCA) for attitude determination, accelerometers (ACC) for measuring the non-gravitational forces and global positioning system (GPS) receivers for precise orbit determination (POD) [5]. A unique feature of GRFO is the laser ranging interferometer (LRI), which validates the key technologies for future gravity missions and space-based gravitational wave antennas [6,7]. Among these payloads, the KBR and LRI systems simultaneously and independently measure the distance variations between the two polar-orbit low-earth-orbit (LEO) spacecrafts. These measurements, essential for recovering the earth’s gravity field, are affected by internal error sources, such as satellite attitude jitter and carrier frequency instability. The decoupling between KBR range and attitude requires the precise determination of the KBR antenna offset vectors (AOVs) on individual GRFO satellites. This study presents a new frequency-domain approach for estimating the AOVs by revisiting multipath noise, a dominant error in AOV estimation that has not been fully evaluated in previous studies [8,9,10,11].

The KBR system operates in a dual one-way ranging (DOWR) scheme, transmitting and receiving K- and Ka-band microwave between the antenna phase centers (APCs) of two satellites [12]. The APC of each satellite is offset from its center of mass (CM) by the AOV, and variations in the microwave link due to satellite attitude jitter introduce time-varying errors that couple with the earth’s gravity field. In the Level-1 processing workflow [13,14], the corrections for these errors are performed by projecting the microwave link onto the line-of-sight (LOS) baseline defined by the twin satellites’ CMs, resulting in the so-called antenna offset correction (AOC) for Level-1B KBR (KBR1B) data products. To estimate the AOVs, GRFO conducted both the pre-launch ground calibration and on-board KBR calibration maneuvers on 17 September 2020 and 28 September 2020, using maneuver strategies similar to those employed for GRACE. The ground calibration results, stored in the VKB1B dataset, provide baseline AOV values.

Previous methods for AOV estimation, such as those outlined by Huang [8,9], relied on time-domain analyses that combined maneuver data with their mirror counterparts to suppress multipath noise. While effective to some extent, these methods were limited in their ability to fully suppress multipath interference. Recent advancements in the recalibration of systematic errors in GRACE KBR data have demonstrated improvements. For instance, recalibration has been shown to reduce biases in pitch and yaw alignment angles, enhance gravity field solutions by up to 50% for zonal coefficients and improve agreement with in situ ocean bottom pressure observations, as well as estimates of Antarctic ice sheet mass variations [10]. Building on these developments, Ellmer [15] introduced a more sophisticated approach by combining total least-squares adjustment (TLSA), a priori calibrated values and a stochastic model to co-estimate the KBR APC alongside the Stokes coefficients of the gravity field. This method revealed biases in the x-coordinate estimates, which were partially reduced—though not entirely eliminated—by incorporating AOC covariance matrices. The application of TLSA further corrected these biases, achieving near-complete mitigation. However, our analysis indicates that even these advanced techniques may not entirely eliminate multipath noise, suggesting the need for further refinement in AOV estimation methodologies.

Satellite attitude variations also contribute to tilt-to-length (TTL) coupling, a significant noise in LRI ranging measurement at low frequencies. TTL coupling arises from misalignment in the satellite’s CM and LRI’s reference reflection point, thereby introducing errors in the ranging data. To achieve minimal TTL coupling, active beam pointing based on differential wave front sensing compensates spacecraft attitude fluctuations; thus, the TTL coupling in the LRI ranging measurement was proven to be lower than 150 µm/rad [5]. Furthermore, due to the laser-based measuring scheme with phase locking, the Level-1B LRI (LRI1B) measurements are free from multipath noise [16].

To revisit the effects of multipath noise on AOV estimation, we use the K- and Ka-band double-difference signals derived from GRFO’s Level-1A data products. These signals are formed by combining four one-way KBR phase measurements with the frequencies of each carrier as the scales [17]. Our analysis shows that this signal reveals the spectral characteristics of multipath noise during KBR calibration.

This study presents a new algorithm to calibrate the AOVs of individual GRFO satellites in the frequency domain using the measurements during on-board KBR calibration maneuver by analyzing and addressing multipath noise. Section 2 establishes the observation equations for KBR calibration. Section 3 and Section 4 detail the KBR calibration maneuver strategy and then analyze the spectral characteristics of AOCs and double-difference signals, which inform our frequency-domain approach. In Section 5, we propose an algorithm to estimate the optimal angular frequency of each sub-maneuver, the multipath noise’s amplitude at this frequency and the AOV of each satellite. Section 6 presents the results, introduces a model-independent method to calibrate multipath noise in the inter-satellite microwave link and compares the LRI ranging measurement with the KBR ranging measurement corrected by recalibrated AOVs from this paper and the AOVs provided by VKB1B.

2. Observational Equation

As illustrated in Figure 1, for each satellite ($i=C, D$), the APC (red point) is deviated from its CM (blue point) by the AOV ${\mathit{d}}_i$, and the KBR system measures the distance between the individual APCs. However, the inter-satellite range between the CMs ${\rho }_{COM}$ is required for gravity field recovery. Therefore, for the misalignment of the inter-satellite pointing of the two satellites, the AOVs ${\mathit{d}}_C$ and ${\mathit{d}}_D$ are projected on the LOS to correct the KBR range ${\rho }_{KBR}$ between two APCs, and this correction is the so-called AOC [15,18].

The AOC $\Delta {\rho }_{AOC}$ is given by

$$\begin{array}{c}\Delta {\rho }_{AOC}=⟨{\mathit{u}}_D,{\mathit{R}}_D{\mathit{d}}_D⟩-⟨{\mathit{u}}_C,{\mathit{R}}_C{\mathit{d}}_C⟩,\end{array}$$

(1)

where ${\mathit{r}}_C$ and ${\mathit{r}}_D$ are the positions of the twin satellites in the inertial frame; ${\mathit{R}}_C$ and ${\mathit{R}}_D$ are their orientation from the satellite reference frame (SRF) to the inertial reference frame (IRF); and ${\mathit{u}}_C$ and ${\mathit{u}}_D$ are the AOVs in individual SRFs,. $⟨ ⟩$ denotes the inner product of two vectors.

In addition to this geometric correction, the distance between both satellites’ CMs ${\rho }_{COM}$ at the same epoch differs from the measured KBR biased range ${\rho }_{KBR}$ due to the travel time of the microwave between the twin spacecrafts. This is corrected by the light-time correction (LTC), denoted as $\Delta {\rho }_{KLTC}$. Taking into account other error sources, including integer ambiguity $N$, multipath noise $\xi$ and stochastic noise $e$, the GRFO range observations are described by the following ionosphere-free combination [18]

$$\begin{array}{ll}{\rho }_{COM}& ={\displaystyle \frac{16}{7}}\left({\rho }_{KBR}^{Ka}+{\xi }_{Ka}+N_{Ka}+e_{Ka}\right)-{\displaystyle \frac{9}{7}}\left({\rho }_{KBR}^K+{\xi }_K+N_K+e_K\right)+\Delta {\rho }_{KLTC}+\Delta {\rho }_{AOC}\\ & ={\rho }_{KBR}+\Delta {\rho }_{KLTC}+\Delta {\rho }_{AOC}+\xi +N+e,\end{array}$$

(2)

where ${\rho }_{KBR}={\displaystyle \frac{16}{7}}{\rho }_{KBR}^{Ka}-{\displaystyle \frac{9}{7}}{\rho }_{KBR}^K$, $\xi ={\displaystyle \frac{16}{7}}{\xi }_{Ka}-{\displaystyle \frac{9}{7}}{\xi }_K$, $N={\displaystyle \frac{16}{7}}{N}_{Ka}-{\displaystyle \frac{9}{7}}{N}_K$, $e={\displaystyle \frac{16}{7}}{e}_{Ka}-{\displaystyle \frac{9}{7}}{e}_K$, ${\rho }_{KBR}^K$ and ${\rho }_{KBR}^{Ka}$ are the KBR biased range; ${\xi }_K$ and ${\xi }_{Ka}$ are the multipath noise; $N_K$ and $N_{Ka}$ are the integer ambiguities; and $e_K$ and $e_{Ka}$ are the stochastic noise for K- and Ka-bands.

Equations (1) and (2) imply the relationship between the AOVs and the other geometrical quantities given by the GRFO data product. For a time span $\{t_1,t_2,\dots , t_m\}$, we rewrite Equations (1) and (2) to obtain the observation equation

$$\begin{array}{c}\mathit{b}=\mathit{A}\mathit{x}\end{array}$$

(3)

where

$$\mathit{A}=\left[\begin{array}{cc}\mathit{u}\left(t_1\right){\mathit{R}}_D\left(t_1\right)& -\mathit{u}\left(t_1\right){\mathit{R}}_C\left(t_1\right)\\ \mathit{u}\left(t_2\right){\mathit{R}}_D\left(t_2\right)& -\mathit{u}\left(t_2\right){\mathit{R}}_C\left(t_2\right)\\ ⋮& ⋮\\ \mathit{u}\left(t_m\right){\mathit{R}}_D\left(t_m\right)& -\mathit{u}\left(t_m\right){\mathit{R}}_C\left(t_m\right)\end{array}\right]$$

(4)

$$\mathit{x}=\left[\begin{array}{c}{\mathit{d}}_C\\ {\mathit{d}}_D\end{array}\right]={\left[d_{C_x},d_{C_y},d_{C_z},d_{D_x},d_{D_y},d_{D_z}\right]}^{⊺},$$

(5)

$$\mathit{b}=\left[\begin{array}{c}{\rho }_{COM}\left(t_1\right)-{\rho }_{KBR}\left(t_1\right)-\Delta {\rho }_{KLTC}\left(t_1\right)-\xi \left(t_1\right)-N-e\left(t_1\right)\\ {\rho }_{COM}\left(t_2\right)-{\rho }_{KBR}\left(t_2\right)-\Delta {\rho }_{KLTC}\left(t_2\right)-\xi \left(t_2\right)-N-e\left(t_2\right)\\ ⋮\\ {\rho }_{COM}\left(t_m\right)-{\rho }_{KBR}\left(t_m\right)-\Delta {\rho }_{KLTC}\left(t_m\right)-\xi \left(t_m\right)-N-e\left(t_m\right)\end{array}\right]$$

(6)

From Equation (3), the AOVs in individual SRFs can be estimated by solving the linear system. For clarity, the data sources of the quantities in $\mathit{A}$ and $\mathit{b}$ are listed in

Table 1. Data sources.

Quantity	Data Product	Description
$\mathit{u}$	GNI1B	The unit LOS vector
${\rho }_{COM}$	GNI1B	The inter-satellite range between two CMs
${\rho }_{KBR}$	KBR1B	The KBR biased range between two APCs
${\rho }_{KLTC}$	KBR1B	The light-time correction for ${\rho }_{KBR}$
${\mathit{R}}_C$	SCA1B	The rotation matrix from IRF to SRF of GRACE-C
${\mathit{R}}_D$	SCA1B	The rotation matrix from IRF to SRF of GRACE-D

.

3. KBR Calibration Parameter

Although the APCs of the twin satellites are precalibrated on the ground before launch, the on-orbit calibration is necessary due to the inevitable changes in APC positions in the space environment. Based on the analysis in [8], the periodical multipath noise during the calibration maneuvers is not negligible (see Section 4.2 for the mechanism). There are two ways to suppress the multipath error. The first is to build an appropriate multipath noise model, and the second is to reduce the magnitude of multipath error by mirror maneuvers. The GRFO calibration design team adopted the latter option.

In 2020, the on-board KBR calibration maneuver was conducted, consisting of eight sub-maneuvers.

Table 2. KBR calibration parameters.

Index (k)	Starting Time (UTC)	GRACE ID	Direction	${\overline{\mathit{\eta }}}_{\mathit{k}}$ (°)	${\overline{\mathsf{\Theta }}}_{\mathit{k}}$ (°)	${\mathit{T}}_{\mathit{k}}$ (s)	Duration (s)
$1^{\prime }$	2020-09-17T05:37:00	C	−pitch	−2	1	250	3750
1	2020-09-17T08:46:00	C	+pitch	+2	1	250	3750
$2^{\prime }$	2020-09-17T11:55:00	C	−yaw	−2	1	250	3750
2	2020-09-17T15:04:00	C	+yaw	+2	1	250	3750
$3^{\prime }$	2020-09-28T05:05:00	D	−pitch	−2	1	250	3750
3	2020-09-28T08:15:00	D	+pitch	+2	1	250	3750
$4^{\prime }$	2020-09-28T11:55:00	D	−yaw	−2	1	250	3750
4	2020-09-28T14:35:00	D	+yaw	+2	1	250	3750

provides the details of the GRFO’s KBR calibration maneuver schedule. For simplicity, the eight sub-maneuvers are indexed by $k$, with the mirror sub-maneuvers denoted using a prime symbol. Each $k$-th sub-maneuver involves the following two steps:

-

Initial Rotation: Rotate one GRFO satellite along the yaw/pitch axis relative to the nominal attitude by an angle bias ${\overline{\eta }}_k$.

-

Swing Maneuver: Swing this satellite along the same axis, and its attitude, ${\overline{\theta }}_k$, varies according to Equation (7).

$$\begin{array}{c}{\overline{\theta }}_k={\overline{\eta }}_k+{\overline{\mathsf{\Theta }}}_k\sin \left({\displaystyle \frac{2\pi }{T_k}}t\right)={\overline{\eta }}_k+{\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi f_{k}t\right)={\overline{\eta }}_k+{\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}\left({\overline{\omega }}_{k}t\right),\end{array}$$

(7)

where ${\overline{\eta }}_k$ is the angle bias; ${\overline{\mathsf{\Theta }}}_k$ is the maneuver amplitude; $T_k$ is the maneuver period; $f_k$ is the maneuver frequency; and ${\overline{\omega }}_k$ is the maneuver angular frequency. Figure 2 shows the attitudes relative to the LOS vector on 17 September 2020 and 28 September 2020.

During the KBR calibration maneuvers, the satellite under calibration adjusted its attitude using attitude control thrusters. In addition to the required angular accelerations, the thruster firings would also inevitably cause linear accelerations. Both the KBR biased range and the inter-satellite range between two CMs measure the range variation caused by the linear accelerations. The range variations in ${\rho }_{KBR}$ and ${\rho }_{COM}$ cancel each other in the vector $\mathit{b}$ in Equation (3). Thus, the thruster-induced range variation does not affect the KBR calibration.

4. Spectral Characteristics

Based on the calibration maneuver strategy above, for the $k$-th sub-maneuver, the satellite attitude follows a sinusoidal variation. In the frequency domain, the attitude’s energy is entirely localized at the frequencies $\pm {\overline{\omega }}_k$. After the geometry projection represented in Equation (1), the attitude’s spectral characteristics are transmitted into the AOCs, providing useful information in the estimation of the AOVs. Meanwhile, as outlined in [8], multipath noise remains the primary error source in KBR calibration. In the following sections, we propose a novel method to evaluate and suppress the multipath effect by introducing the double-difference signal of the KBR phase signals.

4.1. Spectral Characteristics of AOCs

During the calibration maneuver, the misalignment of the inter-satellite pointing results in a rotation transformation from the SRF to the IRF. The rotation matrix for the $i$-th satellite is expressed as

$$\begin{array}{c}{\mathit{R}}_i={\mathit{R}}_i^{LOSF-IRF}{\mathit{R}}_i^{SRF-LOSF},\end{array}$$

(8)

where ${\mathit{R}}_i^{LOSF-IRF}$ represents the rotation matrix from the line-of-sight frame (LOSF) of the $i$-th satellite to the IRF, and ${\mathit{R}}_i^{SRF-LOSF}$ represents the rotation matrix from the SRF to the LOSF of the $i$-th satellite. Since the nominal attitude of GRFO satellites deviates from the LOS only by a few mrad, we neglect the attitude jitter for simplicity. Therefore, when the $i$-th satellite is under the $k$-th sub-maneuver, $({\mathit{R}}_i^{SRF-LOSF}{)}^k$ is derived from Equation (7) directly

$$\begin{array}{c}\begin{array}{c}({\mathit{R}}_i^{SRF-LOSF}{)}^k=\{\begin{array}{c}{\mathit{R}}_y\left({\overline{\theta }}_k\right),\left(k=1,1^{\prime },3,3^{\prime }\right)\\ {\mathit{R}}_z\left({\overline{\theta }}_k\right),\left(k=2,2^{\prime },4,4^{\prime }\right)\end{array}, \end{array}\end{array}$$

(9)

where

$${\mathit{R}}_y\left(\theta \right)=\left[\begin{array}{ccc}\cos \theta & 0& -\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ 0& 1& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\theta & 0& \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right],$$

(10)

$${\mathit{R}}_z\left(\theta \right)=\left[\begin{array}{ccc}\cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right],$$

(11)

For the satellite not under maneuver, $({\mathit{R}}_j^{SRF-LOSF}{)}^k$ is approximated as a unit diagonal matrix because it maintains a nominal attitude, that is

$$\begin{array}{c}\begin{array}{c}\begin{array}{c}({\mathit{R}}_j^{SRF-LOSF}{)}^k=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\end{array}, \end{array}\end{array}$$

(12)

Additionally

$${\left(\mathit{u}{\mathit{R}}_i^{SRF-IRF}\right)}^{⊺}={\mathit{R}}_i^{LOSF-SRF}{\mathit{R}}_i^{IRF-LOSF}{\mathit{u}}^{⊺},$$

(13)

where ${\mathit{R}}_i^{IRF-LOSF}{\mathit{u}}^{⊺}$ indicates the transposed unit LOS vector in the $i$-th satellite’s LOSF, that is, ${\mathit{R}}_C^{IRF-LOSF}{\mathit{u}}^{⊺}=[-\mathrm{1,0},0{]}^{⊺}$, ${\mathit{R}}_D^{IRF-LOSF}{\mathit{u}}^{⊺}=[\mathrm{1,0},0{]}^{⊺}$.

Furthermore, the elements related to the maneuver angle ${\overline{\theta }}_k$ occupy two rows of ${\mathit{R}}_C^{SRF-LOSF}$; therefore, only two columns in the model matrix $\mathit{A}$ are sensitive to one sub-maneuver. The two sensitive rows of $\mathit{x}$ are denoted as $x_{s_{1:}}$ and $x_{s_{2:}}$, and $x_{s_{1:}}=d_{i_x}$ is the x component of AOV of the $i$-th satellite.

Substitute Equations (9) and (12) into Equation (1), and the AOCs for the $k$-th sub-maneuver are expanded as follows:

$$\begin{array}{ll}\Delta {\rho }_{AOC}^k\left(t\right)& \approx x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\theta }}_k-x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\theta }}_k+d_{j_x}\\ & =x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}({\overline{\eta }}_k+{\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t\left)\right)\\ & -x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\eta }}_k+{\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t\left)\right)+d_{j_x}\\ & =(x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k-x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k)\mathrm{c}\mathrm{o}\mathrm{s}({\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t\left)\right)\\ & -(x_{s_{1:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k+x_{s_{2:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k)\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t\left)\right)+d_{j_x},\end{array}(k=\mathrm{1,2},\mathrm{3,4})$$

(14)

Because the $\overline{\mathsf{\Theta }}\mathrm{s}\mathrm{i}\mathrm{n}\left({\overline{\omega }}_{k}t\right)$ in Equation (14) is small, we obtain

$$\begin{array}{ll}\mathrm{c}\mathrm{o}\mathrm{s}({\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t))& =1-2{\mathrm{s}\mathrm{i}\mathrm{n}}^2({\displaystyle \frac{{\overline{\mathsf{\Theta }}}_k}{2}}\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t))\approx 1-2({\displaystyle \frac{{\overline{\mathsf{\Theta }}}_k}{2}}\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t))\\ & \approx 1-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}+{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\mathrm{c}\mathrm{o}\mathrm{s}(2{\overline{\omega }}_{k}t),\end{array}$$

(15)

$$\begin{array}{c}\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}({\overline{\omega }}_{k}t\left)\right)\approx {\overline{\mathsf{\Theta }}}_k\mathrm{s}\mathrm{i}\mathrm{n}\left({\overline{\omega }}_{k}t\right).\end{array}\end{array}$$

(16)

Then, Equation (1) is approximated as follows:

$$\begin{array}{ll}\Delta {\rho }_{AOC}^k\left(t\right)& \approx d_{j_x}+\left(x_{s_{1:}}\cos {\overline{\eta }}_k-x_{s_{2:}}\sin {\overline{\eta }}_k\right)\left(1-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\right)\\ & -\overline{\mathsf{\Theta }}(x_{s_{1:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k+x_{s_{2:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k)\mathrm{s}\mathrm{i}\mathrm{n}\left({\overline{\omega }}_{k}t\right)\\ & +{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}(x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k-x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k)\mathrm{c}\mathrm{o}\mathrm{s}(2{\overline{\omega }}_{k}t),\end{array}$$

(17)

The three lines on the right-hand side of Equation (17) are the direct current component $|\Delta {\tilde{\rho }}_{AOC}^k(0\left)\right|$ and the harmonic components at the maneuver frequency $|\Delta {\tilde{\rho }}_{AOC}^k({\overline{\omega }}_k\left)\right|$ and the doubling maneuver frequency $|\Delta {\tilde{\rho }}_{AOC}^k(2{\overline{\omega }}_k\left)\right|$, respectively (Figure 3b,d).

$$\begin{array}{c}|\Delta {\tilde{\rho }}_{AOC}^k(0)|\approx d_{j_x}+\left(x_{s_{1:}}\cos {\overline{\eta }}_k-x_{s_{2:}}\sin {\overline{\eta }}_k\right)\left(1-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\right),\end{array}$$

(18)

$$\begin{array}{c}|\Delta {\tilde{\rho }}_{AOC}^k({\overline{\omega }}_k\left)\right|\approx {\overline{\mathsf{\Theta }}}_k(x_{s_{1:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k+x_{s_{2:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k),\end{array}$$

(19)

$$\begin{array}{c}|\Delta {\tilde{\rho }}_{AOC}^k(2{\overline{\omega }}_k)|\approx {\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\left(x_{s_{1:}}\cos {\overline{\eta }}_k-x_{s_{2:}}\sin {\overline{\eta }}_k\right),\end{array}$$

(20)

The same derivation can be easily adapted to obtain the approximation of its mirror sub-maneuver, i.e., $k^{\prime }$

$$\begin{array}{ll}\Delta {\rho }_{AOC}^{k^{\prime }}\left(t\right)& \approx d_{j_x}+\left(x_{s_{1:}}\cos {\overline{\eta }}_{k^{\prime }}-x_{s_{2:}}\sin {\overline{\eta }}_{k^{\prime }}\right)\left(1-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_{k^{\prime }}{)}^2}{4}}\right)\\ & -{\overline{\mathsf{\Theta }}}_{k^{\prime }}(x_{s_{1:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_{k^{\prime }}+x_{s_{2:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }})\mathrm{s}\mathrm{i}\mathrm{n}\left({\overline{\omega }}_{k^{\prime }}t\right)\\ & +{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_{k^{\prime }}{)}^2}{4}}(x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }}-x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_{k^{\prime }})\mathrm{c}\mathrm{o}\mathrm{s}(2{\overline{\omega }}_{k^{\prime }}t),\end{array}$$

(21)

In Table 2, $\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k=\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }}$ and $\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k=-\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_{k^{\prime }}$; therefore, the coefficients of $x_{s_{1:}}$ in Equations (17) and (21) are the opposites of each other. Thus, we solve $x_{s_{1:}}$ with the linear combination of Equations (17) and (21). However, $|\Delta {\tilde{\rho }}_{AOC}^k(0\left)\right|$ and $|\Delta {\tilde{\rho }}_{AOC}^{k^{\prime }}(0\left)\right|$ are coupled with integer ambiguity $N$ in the KBR biased range ${\rho }_{KBR}$, and the x component of AOV $x_{s_{1:}}$ dominates the harmonic components at the doubling maneuver frequency:

$$\begin{array}{c}{\displaystyle \frac{x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k}{x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k}}\approx {10}^{-5},\end{array}$$

(22)

Therefore, it is the harmonic components at the maneuver frequency $|\Delta {\tilde{\rho }}_{AOC}^k({\overline{\omega }}_k\left)\right|$ and $|\Delta {\tilde{\rho }}_{AOC}^{k^{\prime }}({\overline{\omega }}_{k^{\prime }}\left)\right|$ that give a hint to solve $x_{s_{1:}}$.

Meanwhile, the AOCs $\Delta {\mathit{\rho }}_{AOC}$ are not direct observations of the GRFO but are viewed as a derived quantity from ${\mathit{\rho }}_{COM}$, ${\mathit{\rho }}_{KBR}$ and $\Delta {\mathit{\rho }}_{KLTC}$, as shown in Equation (2). Therefore, in the time domain, Equation (3) indicates that $\mathit{b}$ is a version of $\Delta {\mathit{\rho }}_{AOC}$ with stochastic noise $\mathit{e}$, integer ambiguity $\mathit{N}$ and multipath noise $\mathit{\xi }$. In the frequency domain, the integer ambiguity $\mathit{N}$ exerts no effect on $\left|\tilde{b}\right({\overline{\omega }}_k\left)\right|$ and $\left|\tilde{b}\left(2{\overline{\omega }}_k\right)\right|$, and the harmonic component of stochastic noise $\mathit{e}$ at the maneuver frequency $\left|\tilde{e}\right({\overline{\omega }}_k\left)\right|<2 \mathsf{\mu }\mathrm{m}$; thus, the boresight error caused by $\left|\tilde{e}\right({\overline{\omega }}_k\left)\right|$ is less than 100 $\mathsf{\mu }$rad. Multipath noise is analyzed in the following section. Figure 3 presents the time series and ASDs of the measurement vector $\mathit{b}$. Outside of the eight shaded regions in Figure 3a,c, notable low-frequency noise is introduced by the POD measurements.

We substitute $\Delta {\mathit{\rho }}_{AOC}$ with the vector $\mathit{b}$ (Figure 3) to estimate the AOVs:

$$\begin{array}{c}|{\tilde{b}}^k({\overline{\omega }}_k)|\approx {\overline{\mathsf{\Theta }}}_k(x_{s_{1:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k+x_{s_{2:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k)+|{\tilde{\xi }}^k({\overline{\omega }}_k)|,\end{array}$$

(23)

$$\begin{array}{c}\begin{array}{c}|{\tilde{b}}^{k^{\prime }}({\overline{\omega }}_{k^{\prime }})|\approx -{\overline{\mathsf{\Theta }}}_{k^{\prime }}(x_{s_{1:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_{k^{\prime }}+x_{s_{2:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }})+|{\tilde{\xi }}^{k^{\prime }}({\overline{\omega }}_{k^{\prime }})|,\end{array}\end{array}$$

(24)

$$\begin{array}{c}|{\tilde{b}}^k(2{\overline{\omega }}_k)|\approx {\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}(x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k-x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k)+|{\tilde{\xi }}^k(2{\overline{\omega }}_k)|,\end{array}$$

(25)

$$\begin{array}{c}|{\tilde{b}}^{k^{\prime }}(2{\overline{\omega }}_{k^{\prime }})|\approx {\displaystyle \frac{({\overline{\mathsf{\Theta }}}_{k^{\prime }}{)}^2}{4}}(x_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }}-x_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_{k^{\prime }})+|{\tilde{\xi }}^{k^{\prime }}(2{\overline{\omega }}_{k^{\prime }})|.\end{array}$$

(26)

where ${\xi }_k$ and ${\xi }_{k^{\prime }}$ are multipath noise during the $k$-th sub-maneuver and its mirror sub-maneuver.

In Equations (23) and (24), there are four unknown parameters: (i) the AOV-related parameters $x_{s_{1:}}$ and $x_{s_{2:}}$ and (ii) the multipath noise $\left|{\tilde{\xi }}^k\left({\overline{\omega }}_k\right)\right|$ and $\left|{\tilde{\xi }}^{k^{\prime }}\left({\overline{\omega }}_{k^{\prime }}\right)\right|$. Once the $|{\tilde{\xi }}^k({\overline{\omega }}_k)|$ and $|{\tilde{\xi }}^{k^{\prime }}({\overline{\omega }}_{k^{\prime }})|$ are determined, $x_{s_{1:}}$ and $x_{s_{2:}}$ are straightforwardly solved. Wang treated the multipath noise during the $k$-th sub-maneuver and its mirror sub-maneuver as equal, i.e., $|{\tilde{\xi }}^k({\overline{\omega }}_k)|=|{\tilde{\xi }}^{k^{\prime }}({\overline{\omega }}_{k^{\prime }})|$ [8]. $x_{s_{2:}}$ is solved from the difference between Equations (23) and (24):

$$\begin{array}{c}{x}_{s_{2:}}={\displaystyle \frac{\left|{\tilde{b}}^k\left({\overline{\omega }}_k\right)\right|-\left|{\tilde{b}}^{k^{\prime }}\left({\overline{\omega }}_{k^{\prime }}\right)\right|}{{\overline{\mathsf{\Theta }}}_k\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k+{\overline{\mathsf{\Theta }}}_{k^{\prime }}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }}}},\end{array}$$

(27)

However, from Equations (22), (25) and (26), we obtain the following relationship:

$$\begin{array}{c}{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_{k^{\prime }}{)}^2}{4}}{x}_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }}={\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}{x}_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k\gg {\displaystyle \frac{({\overline{\mathsf{\Theta }}}_{k^{\prime }}{)}^2}{4}}{x}_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_{k^{\prime }}=-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}{x}_{s_{2:}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k,\end{array}$$

(28)

We further evaluate the harmonic components of multipath noise at $2{\overline{\omega }}_k$ by

$$\begin{array}{c}\left|{\tilde{\xi }}^k\left(2{\overline{\omega }}_k\right)\right|\approx \left|{\tilde{b}}^k\left(2{\overline{\omega }}_k\right)\right|-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}{x}_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k,\end{array}$$

(29)

$$\begin{array}{c}\left|{\tilde{\xi }}^{k^{\prime }}\left(2{\overline{\omega }}_{k^{\prime }}\right)\right|\approx \left|{\tilde{b}}^{k^{\prime }}\left(2{\overline{\omega }}_{k^{\prime }}\right)\right|-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_{k^{\prime }}{)}^2}{4}}{x}_{s_{1:}}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_{k^{\prime }},\end{array}$$

(30)

In the next section, we will discuss the relationship between $\left|{\tilde{\xi }}^k\left({\overline{\omega }}_k\right)\right|$ and $\left|{\tilde{\xi }}^k\left(2{\overline{\omega }}_k\right)\right|$ by leveraging the double-difference signal. Using this relationship, we integrate Equations (23) and (24) with two additional equations (the relationship between $\left|{\tilde{\xi }}^k\left({\overline{\omega }}_k\right)\right|$ and $\left|{\tilde{\xi }}^k\left(2{\overline{\omega }}_k\right)\right|$ and the relationship between $\left|{\tilde{\xi }}^{k^{\prime }}\left({\overline{\omega }}_{k^{\prime }}\right)\right|$ and $\left|{\tilde{\xi }}^{k^\prime }\left(2{\overline{\omega }}_{k^{\prime }}\right)\right|$) to solve the four aforementioned unknown parameters.

4.2. Spectral Characteristics of Double-Difference Signal

The microwave wavefront received by one GRFO satellite is approximately planar and perpendicular to the LOS vector. During each sub-maneuver, the multipath noise of the satellite under calibration exhibits spectral characteristics similar to those of the AOC ($\mathit{b}$). In contrast, the satellite not under calibration does not exhibit these spectral characteristics.

This alignment of spectral characteristics indicates that multipath noise is one of the predominant contributors to errors in estimating the AOVs, as it shares the same period as $\mathit{b}$. To mitigate the impact of multipath noise on AOV estimation, we introduce the double-difference combination of KBR phase signals. This approach effectively reduces the influence of multipath noise, enabling a more accurate estimation of AOVs.

One-way phase measurements of K- and Ka-bands at any time $t$ contain gravity signal proportional to the corresponding K/Ka-band carrier frequency. This enables us to eliminate the gravity field signal, carrier frequency noise, high-frequency noise of the ionosphere and time-tag noise by combining four one-way phase measurements with K- and Ka-bands carrier frequencies as the scaling factors. Thus, this combination can be used to evaluate the stochastic noise and residual multipath noise. We call this linear combination the double-difference signal [17]:

$$\begin{array}{ll}\sigma \left(t\right)& ={\displaystyle \frac{c}{f_C^K+f_D^K}}\left[\left({\varphi }_K^{D-C}\left(t\right)-{\displaystyle \frac{3}{4}}{\varphi }_{Ka}^{D-C}\left(t\right)\right)-\left({\varphi }_K^{C-D}\left(t\right)-{\displaystyle \frac{3}{4}}{\varphi }_{Ka}^{C-D}\left(t\right)\right)\right]\\ & =\left({\rho }_K^{D-C}\right(t)-{\rho }_{Ka}^{D-C}(t\left)\right)-\left({\rho }_K^{C-D}\left(t\right)-{\rho }_{Ka}^{C-D}\left(t\right)\right),\end{array}$$

(31)

where $f_C^K$ and $f_D^K$ represent the K-band carrier frequency for each GRFO satellite; ${\varphi }_K^{C-D}\left(t\right)$ and ${\varphi }_{Ka}^{C-D}\left(t\right)$ indicate separately the K- and Ka-band carrier phase signal transmitted by GRACE-D and received by GRACE-C; ${\varphi }_K^{D-C}\left(t\right)$ and ${\varphi }_{Ka}^{D-C}\left(t\right)$ indicate separately the K- and Ka-band carrier phase signal transmitted by GRACE-C and received by GRACE-D; ${\rho }_K^{C-D}\left(t\right)$, ${\rho }_{Ka}^{C-D}\left(t\right)$, ${\rho }_K^{D-C}\left(t\right)$ and ${\rho }_{Ka}^{D-C}\left(t\right)$ are the corresponding inter-satellite ranges. Figure 4 illustrates the double-difference signals during each sub-maneuver in the time and frequency domain. The harmonic components at the maneuver frequency and its doubling are notable for all sub-maneuvers, except the positive pitch pone for GRACE-C.

During the $k$-th sub-maneuver, the double-difference signal is expanded as follows:

$$\begin{array}{ll}{\sigma }^k\left(t\right)& ={\xi }_K^k\left(t\right)-{\xi }_{Ka}^k\left(t\right)+e_{{\sigma }^k}+N_{{\sigma }^k}\\ & +2{\overline{\mathsf{\Theta }}}_k\sin \left({\displaystyle \frac{{\overline{\omega }}_k}{2}}\tau \right)\left[\left(x_{s_{1:}}^K-x_{s_{1:}}^{Ka}\right)\sin {\overline{\eta }}_k+\left(x_{s_{2:}}^K-x_{s_{2:}}^{Ka}\right)\cos {\overline{\eta }}_k\right]\cdot \cos \left({\overline{\omega }}_k\left(t+{\displaystyle \frac{\tau }{2}}\right)\right)\\ & +{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{2}}\sin \left({\overline{\omega }}_k\tau \right)\left[\left(x_{s_{1:}}^K-x_{s_{1:}}^{Ka}\right)\cos {\overline{\eta }}_k-\left(x_{s_{2:}}^K-x_{s_{2:}}^{Ka}\right)\sin {\overline{\eta }}_k\right]\cdot \sin \left(2{\overline{\omega }}_k\left(t+{\displaystyle \frac{\tau }{2}}\right)\right)\\ & \approx {\xi }_K^k\left(t\right)-{\xi }_{Ka}^k\left(t\right)+e_{{\sigma }^k}+N_{{\sigma }^k},\end{array}$$

(32)

where ${\xi }_K^k$ and ${\xi }_{Ka}^k$ represent the K- and Ka-band multipath noise during the $k$-th sub-maneuver; $e_{{\sigma }^k}$ represents the stochastic noise and residual ionosphere noise; $N_{{\sigma }^k}$ represents the residual integer ambiguity; and $\tau$ represents the flight time of carrier. For GRFO, $\tau \approx 1\times {10}^{-3}$ second, then ${\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{2}}\sin \left({\overline{\omega }}_k\tau \right)\approx 3\times {10}^{-9}$ and $2{\overline{\mathsf{\Theta }}}_k\sin \left({\displaystyle \frac{{\overline{\omega }}_k}{2}}\tau \right)\approx 4\times {10}^{-7}$. These values indicate that the fifth and sixth terms in Equation (32) are negligibly small and can be disregarded in the analysis. Consequently, the double-difference signal reflects the spectral characteristics of multipath noise:

$$|{\tilde{\sigma }}^k({\overline{\omega }}_k)|=\left|\left|{\tilde{\xi }}_K^k\left({\overline{\omega }}_k\right)\right|-\left|{\tilde{\xi }}_{Ka}^k\left({\overline{\omega }}_k\right)\right|\right|,$$

(33)

$$|{\tilde{\sigma }}^k(2{\overline{\omega }}_k)|=\left|\left|{\tilde{\xi }}_K^k\left(2{\overline{\omega }}_k\right)\right|-\left|{\tilde{\xi }}_{Ka}^k\left(2{\overline{\omega }}_k\right)\right|\right|,$$

(34)

The coefficients of K- and Ka-band multipath noise in the double-difference signal ($[1,-1]$) differ from those in the GRFO range observation ($\left[-{\displaystyle \frac{9}{7}},{\displaystyle \frac{16}{7}}\right]$), indicating that the multipath noise effects in the double-difference signal and GRFO range are linearly independent from each other. Meanwhile, during each sub-maneuver, the K- and Ka-band carriers share identical reflection points, and they exhibit consistent reflection angles and attenuation factors. Therefore, we make an assumption as follows:

$$\begin{array}{c}\begin{array}{c}\begin{array}{c}{\displaystyle \frac{\left|{\tilde{\xi }}_K^k\left(2{\overline{\omega }}_k\right)\right|}{\left|{\tilde{\xi }}_K^k\left({\overline{\omega }}_k\right)\right|}}={\displaystyle \frac{\left|{\tilde{\xi }}_{Ka}^k\left(2{\overline{\omega }}_k\right)\right|}{\left|{\tilde{\xi }}_{Ka}^k\left({\overline{\omega }}_k\right)\right|}}={\displaystyle \frac{\left|\left|{\tilde{\xi }}_K^k\left(2{\overline{\omega }}_k\right)\right|-\left|{\tilde{\xi }}_{Ka}^k\left(2{\overline{\omega }}_k\right)\right|\right|}{\left|\left|{\tilde{\xi }}_K^k\left({\overline{\omega }}_k\right)\right|-\left|{\tilde{\xi }}_{Ka}^k\left({\overline{\omega }}_k\right)\right|\right|}}={\alpha }_k\end{array}\end{array},\end{array}$$

(35)

For the $k$-th sub-maneuver, the ratio of harmonic components at ${\overline{\omega }}_k$ and $2{\overline{\omega }}_k$ demonstrates a consistent ratio across both K-band and Ka-band carriers. This ratio is identical to that observed in the double-difference signal, denoted as ${\alpha }_k$.

Note that when GRACE-D is under calibration, i.e., $k=3,3^{\prime },4,4^{\prime }$, the signs of ${\xi }^k$ in the double-difference signal are opposite to those in Equation (26) [17].

5. Methodology

We propose a new algorithm to estimate the harmonic components of multipath noise at the maneuver frequency and the AOVs $\mathit{x}$ together, as outlined in Figure 5.

According to [15], the AOVs of GRACE are affected by solar illumination, and spectral analysis reveals that the period of AOV change corresponds to a full revolution of the longitude of the ascending node for the GRACE orbital plane. Due to the resemblance of the MWIs of GRACE and GRFO, we assume that the AOVs of the GRFO remain constant during KBR calibration.

Taking the calibration from 17 September 2020 as an example, we reformulate Equations (23)–(26) and (33)–(35) to a linear system of equations as follows, and we replace $\left|{\tilde{\xi }}_K^k\right(2{\overline{\omega }}_k\left)\right|$ and $\left|{\tilde{\xi }}_{Ka}^k\right(2{\overline{\omega }}_k\left)\right|$ with ${\alpha }_k\left|{\tilde{\xi }}_K^k\right({\overline{\omega }}_k\left)\right|$ and ${\alpha }_k\left|{\tilde{\xi }}_{Ka}^k\right({\overline{\omega }}_k\left)\right|$.

$$\begin{array}{c}[\mathsf{\Psi },\mathsf{\Pi }]\mathit{\gamma }=\mathit{\beta },\end{array}$$

(36)

where $\mathsf{\Psi }$ is the design matrix related to the AOV ${\mathit{d}}_C$, and $\mathsf{\Pi }$ is the design matrix related to multipath noise.

$$\mathsf{\Psi }=\left[\begin{array}{ccc}{\overline{\mathsf{\Theta }}}_1\sin {\overline{\eta }}_1& 0& {\overline{\mathsf{\Theta }}}_1\cos {\overline{\eta }}_1\\ {\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_1\right)}^2}{4}}\cos {\overline{\eta }}_1& 0& -{\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_1\right)}^2}{4}}\sin {\overline{\eta }}_1\\ 0& 0& 0\\ -{\overline{\mathsf{\Theta }}}_{1^{\prime }}\sin {\overline{\eta }}_{1^{\prime }}& 0& -{\overline{\mathsf{\Theta }}}_{1^{\prime }}\cos {\overline{\eta }}_{1^{\prime }}\\ {\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_{1^{\prime }}\right)}^2}{4}}\cos {\overline{\eta }}_{1^{\prime }}& 0& -{\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_1\right)}^2}{4}}\sin {\overline{\eta }}_1\\ 0& 0& 0\\ {\overline{\mathsf{\Theta }}}_2\sin {\overline{\eta }}_2& {\overline{\mathsf{\Theta }}}_2\cos {\overline{\eta }}_2& 0\\ {\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_2\right)}^2}{4}}\cos {\overline{\eta }}_2& -{\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_2\right)}^2}{4}}\sin {\overline{\eta }}_2& 0\\ 0& 0& 0\\ -{\overline{\mathsf{\Theta }}}_{2^{\prime }}\sin {\overline{\eta }}_{2^{\prime }}& -{\overline{\mathsf{\Theta }}}_{2^{\prime }}\cos {\overline{\eta }}_{\prime }& 0\\ {\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_{2^{\prime }}\right)}^2}{4}}\cos {\overline{\eta }}_{2^{\prime }}& -{\displaystyle \frac{{\left({\overline{\mathsf{\Theta }}}_{2^{\prime }}\right)}^2}{4}}\sin {\overline{\eta }}_{2^{\prime }}& 0\\ 0& 0& 0\end{array}\right],$$

(37)

$$\mathit{\gamma }=[d_{C_x},d_{C_y},d_{C_z},\left|{\tilde{\xi }}_{Ka}^1\left({\overline{\omega }}_1\right)\right|,\left|{\tilde{\xi }}_K^1\left({\overline{\omega }}_1\right)\right|,\left|{\tilde{\xi }}_{Ka}^{1^{\prime }}\left({\overline{\omega }}_{1^{\prime }}\right)\right|,\left|{\tilde{\xi }}_K^{1^{\prime }}\left({\overline{\omega }}_{1^{\prime }}\right)\right|,\left|{\tilde{\xi }}_{Ka}^2\left({\overline{\omega }}_2\right)\right|,\left|{\tilde{\xi }}_K^2\left({\overline{\omega }}_2\right)\right|,\left|{\tilde{\xi }}_{Ka}^{2^{\prime }}\left({\overline{\omega }}_{2^{\prime }}\right)\right|,\left|{\tilde{\xi }}_K^{2^{\prime }}\left({\overline{\omega }}_{2^{\prime }}\right)\right|{]}^{⊺},$$

(38)

$$\mathsf{\Pi }=\left[\begin{array}{cccccccc}{\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}& 0& 0& 0& 0& 0& 0\\ {\displaystyle \frac{16}{7}}{\alpha }_1& -{\displaystyle \frac{9}{7}}{\alpha }_1& 0& 0& 0& 0& 0& 0\\ 1& -1& 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{16}{7}}{\alpha }_{1^{\prime }}& -{\displaystyle \frac{9}{7}}{\alpha }_{1^{\prime }}& 0& 0& 0& 0\\ 0& 0& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{16}{7}}{\alpha }_2& -{\displaystyle \frac{9}{7}}{\alpha }_2& 0& 0\\ 0& 0& 0& 0& 1& -1& 0& 0\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{16}{7}}{\alpha }_{2^{\prime }}& -{\displaystyle \frac{9}{7}}{\alpha }_{2^{\prime }}\\ 0& 0& 0& 0& 0& 0& 1& -1\end{array}\right],$$

(39)

$$\mathit{\beta }=[\left|{\tilde{b}}^1\left({\overline{\omega }}_1\right)\right|, \left|{\tilde{b}}^1\left(2{\overline{\omega }}_1\right)\right|,\left|{\tilde{\sigma }}^1\left({\overline{\omega }}_1\right)\right|,\left|{\tilde{b}}^{1^{\prime }}\left({\overline{\omega }}_{1^{\prime }}\right)\right|,\left|{\tilde{b}}^{1^{\prime }}\left(2{\overline{\omega }}_{1^{\prime }}\right)\right|,\left|{\tilde{\sigma }}^{1^{\prime }}\left({\overline{\omega }}_{1^{\prime }}\right)\right|,\left|{\tilde{b}}^2\left({\overline{\omega }}_2\right)\right|,\left|{\tilde{b}}^2\left(2{\overline{\omega }}_2\right)\right|,\left|{\tilde{\sigma }}^2\left({\overline{\omega }}_2\right)\right|,\left|{\tilde{b}}^{2^{\prime }}\left({\overline{\omega }}_{2^{\prime }}\right)\right|,\left|{\tilde{b}}^{2^{\prime }}\left(2{\overline{\omega }}_{2^{\prime }}\right)\right|,\left|{\tilde{\sigma }}^{2^{\prime }}\left({\overline{\omega }}_{2^{\prime }}\right)\right|{]}^{⊺},$$

(40)

However, several questions need to be addressed:

• Replace elements in $\mathsf{\Psi }$ with measurements: solving ${\mathit{d}}_C$ using Equation (36) requires precise knowledge of $\mathsf{\Psi }$, but it is impossible for the actual angle bias ${\overline{\eta }}_k$ and the actual maneuver amplitude ${\overline{\mathsf{\Theta }}}_k$ to maintain their corresponding nominal values. We denote the two sensitive columns of ${\mathit{A}}^k$ as ${\mathit{A}}_{:s_1}^k$ and ${\mathit{A}}_{:s_2}^k$. The non-zero elements in $\mathsf{\Psi }$ are the theoretical values of the harmonic components at ${\overline{\omega }}_k$ and $2{\overline{\omega }}_k$ of ${\mathit{A}}_{:s_1}^k$ and ${\mathit{A}}_{:s_2}^k$, that is, ideally

  $$\begin{array}{c}\left|{\tilde{A}}_{:s_1}^k\left({\overline{\omega }}_k\right)\right|={\overline{\mathsf{\Theta }}}_k\sin {\overline{\eta }}_k,\end{array}$$

  (41)

  $$\left|{\tilde{A}}_{:s_2}^k\left({\overline{\omega }}_k\right)\right|={\overline{\mathsf{\Theta }}}_k\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\eta }}_k$$

  (42)

  $$\begin{array}{c}\left|{\tilde{A}}_{:s_1}^k\left(2{\overline{\omega }}_k\right)\right|={\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\cos {\overline{\eta }}_k,\end{array}$$

  (43)

  $$\begin{array}{c}\end{array}-\left|{\tilde{A}}_{:s_2}^k\left(2{\overline{\omega }}_k\right)\right|=-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k$$

  (44)

  $$\begin{array}{c}\end{array}-\left|{\tilde{A}}_{:s_2}^k\left(2{\overline{\omega }}_k\right)\right|=-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\eta }}_k.$$

  (45)

  Therefore, we substitute $\left\{\left|{\tilde{A}}_{:s_1}^k\left({\overline{\omega }}_k\right)\right|,\left|{\tilde{A}}_{:s_2}^k\left({\overline{\omega }}_k\right)\right|,\left|{\tilde{A}}_{:s_1}^k\left(2{\overline{\omega }}_k\right)\right|,-\left|{\tilde{A}}_{:s_2}^k\left(2{\overline{\omega }}_k\right)\right|\right\}$ for $\left\{{\overline{\mathsf{\Theta }}}_k\sin {\overline{\eta }}_k,{\overline{\mathsf{\Theta }}}_k\cos {\overline{\eta }}_k,{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\cos {\overline{\eta }}_k,-{\displaystyle \frac{({\overline{\mathsf{\Theta }}}_k{)}^2}{4}}\sin {\overline{\eta }}_k\right\}$.
• Due to the complicated attitude control strategy, it is difficult to evaluate whether the actual maneuver period $T_k$ maintains the nominal value ${\overline{T}}_k$ from the time domain. Direct discrete Fourier transform (DFT) may cause spectral leakage when ${\overline{T}}_k\ne T_k$, thus affecting $\mathsf{\Psi }$. As presented in [19], we estimate harmonic components using the least-squares (LS) method by finding an optimal $T_k^{\star }$.
• The low-frequency noise in POD measurements is another main noise source when estimating the harmonic components of ${\mathit{b}}^k$ (Figure 3). Wang presented a method by fitting the low-frequency noise in ${\mathit{b}}^k$ using a fourth-order polynomial function [8]. For simplicity, we use the second-order difference of ${\mathit{b}}^k$ to eliminate low-frequency noise, which is equivalent to fitting low-frequency noise using a second-order polynomial function. The second-order difference with signal $\mathit{y}$ is denoted as $\ddot{\mathit{y}}$, i.e., $\ddot{y}\left(t_i\right)=y\left(t_{i+2}\right)-2y\left(t_{i+1}\right)+y\left(t_i\right),i=\mathrm{0,1},\dots ,n-2$, and its Fourier transform is denoted as $\tilde{\ddot{y}}$.

We find the optimal ${\omega }_k^{\star }$ to minimize the $L_{\mathrm{2,1}}$ norm of the following error matrix $\mathit{W}$ using a grid search.

$$\begin{array}{c}\end{array}\begin{array}{c}\underset{{\omega }_k^{\star }\in S}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n} }f({\omega }_k^{\star })=\{{\omega }_k^{\star }\in S:f({\omega }_k^{\star })\le f\left({\omega }_k\right) \mathrm{for} \mathrm{all}{ \omega }_k\in S\},\end{array}$$

(46)

where

$$\begin{array}{c}\end{array}\begin{array}{c}\begin{array}{c}f\left({\omega }_k\right)=\Vert \mathit{W}{\Vert }_{\mathrm{2,1}},\end{array}\end{array}$$

(47)

$$\begin{array}{c}\mathit{W}=\mathit{M}-\mathit{H}\times \mathit{Y},\end{array}$$

(48)

$$\mathit{M}=[{\ddot{\mathit{b}}}^k,{\ddot{\mathit{A}}}_{:s_1}^k,{\ddot{\mathit{A}}}_{:s_2}^k,{\ddot{\mathit{\sigma }}}^k],$$

(49)

$$\mathit{H}=\left[\begin{array}{cc}{e}^{i{\omega }_k{t}_0}& e^{i2{\omega }_k{t}_0}\\ e^{i{\omega }_k{t}_1}& e^{i2{\omega }_k{t}_1}\\ ⋮& ⋮\\ e^{i{\omega }_k{t}_{n-2}}& e^{i2{\omega }_k{t}_{n-2}}\end{array}\right],$$

(50)

$$\mathit{Y}=\left[\begin{array}{c}\left|{\tilde{\ddot{b}}}^k\left({\omega }_k\right)\right|,\\ \left|{\tilde{\ddot{b}}}^k\left(2{\omega }_k\right)\right|\end{array}\begin{array}{c}\left|{\tilde{\ddot{A}}}_{:s_1}^k\left({\omega }_k\right)\right|,\\ ,\left|{\tilde{\ddot{A}}}_{:s_1}^k\left(2{\omega }_k\right)\right|,\end{array}\begin{array}{c}\left|{\tilde{\ddot{A}}}_{:s_2}^k\left({\omega }_k\right)\right|,\\ \left|{\tilde{\ddot{A}}}_{:s_2}^k\left(2{\omega }_k\right)\right|,\end{array}\begin{array}{c}\left|{\tilde{\ddot{\sigma }}}^k\left({\omega }_k\right)\right|\\ \left|{\tilde{\ddot{\sigma }}}^k\left(2{\omega }_k\right)\right|\end{array}\right],$$

(51)

where $\left|{\tilde{\ddot{A}}}_{:s_1}^k\left(\omega \right)\right|$, $\left|{\tilde{\ddot{A}}}_{:s_2}^k\left(\omega \right)\right|$, $\left|{\tilde{\ddot{b}}}^k\left(\omega \right)\right|$ and $\left|{\tilde{\ddot{\sigma }}}^k\left(\omega \right)\right|$ represent the harmonic components of ${\ddot{\mathit{b}}}^k$, ${\ddot{\mathit{A}}}_{:s_1}^k$, ${\ddot{\mathit{A}}}_{:s_2}^k$ and ${\ddot{\mathit{\sigma }}}^k$ at $\omega$; $\mathit{H}$ is the inverse discrete Fourier transform (IDFT) functional matrix; $\mathit{M}$ is stacked by ${\ddot{\mathit{b}}}^k$, ${\ddot{\mathit{A}}}_{:s_1}^k$, ${\ddot{\mathit{A}}}_{:s_2}^k$ and ${\ddot{\mathit{\sigma }}}^k$.

Denote the optimal $\mathit{Y}$ as ${\mathit{Y}}^{\star }$, and the $i$-th column of ${\mathit{Y}}^{\star }$ is obtained by

$$\begin{array}{c}{\mathit{Y}}_{:i}^{\star }=\left[\right({\mathit{H}}^{\star }{)}^{⊺}{\mathit{H}}^{\star }{]}^{-1}({\mathit{H}}^{\star }{)}^{⊺}{\mathit{M}}_{:i},(i=\mathrm{1,2},\mathrm{3,4})\end{array}$$

(52)

Here, using Equations (52) and (53), we reformulate Equation (36) into its second-order difference as follows:

$$\begin{array}{c}\left[\ddot{\mathsf{\Psi }},\ddot{\mathsf{\Pi }}\right]\ddot{\mathit{\gamma }}=\ddot{\mathit{\beta }},\end{array}$$

(53)

where

$$\ddot{\mathsf{\Psi }}=\left[\begin{array}{ccc}\left|{\tilde{\ddot{A}}}_{:1}^1\left({\omega }_1^{\star }\right)\right|& 0& \left|{\tilde{\ddot{A}}}_{:3}^1\left({\omega }_1^{\star }\right)\right|\\ \left|{\tilde{\ddot{A}}}_{:1}^1\left(2{\omega }_1^{\star }\right)\right|& 0& -\left|{\tilde{\ddot{A}}}_{:3}^1\left(2{\omega }_1^{\star }\right)\right|\\ 0& 0& 0\\ \left|{\tilde{\ddot{A}}}_{:1}^{1^{\prime }}\left({\omega }_{1^{\prime }}^{\star }\right)\right|& 0& -\left|{\tilde{\ddot{A}}}_{:3}^{1^{\prime }}\left({\omega }_{1^{\prime }}^{\star }\right)\right|\\ \left|{\tilde{\ddot{A}}}_{:1}^{1^{\prime }}\left(2{\omega }_{1^{\prime }}^{\star }\right)\right|& 0& \left|{\tilde{\ddot{A}}}_{:3}^{1^{\prime }}\left(2{\omega }_{1^{\prime }}^{\star }\right)\right|\\ 0& 0& 0\\ \left|{\tilde{\ddot{A}}}_{:1}^2\left({\omega }_2^{\star }\right)\right|& \left|{\tilde{\ddot{A}}}_{:2}^2\left({\omega }_2^{\star }\right)\right|& 0\\ \left|{\tilde{\ddot{A}}}_{:1}^2\left(2{\omega }_2^{\star }\right)\right|& -\left|{\tilde{\ddot{A}}}_{:2}^2\left(2{\omega }_2^{\star }\right)\right|& 0\\ 0& 0& 0\\ \left|{\tilde{\ddot{A}}}_{:1}^{2^{\prime }}\left({\omega }_{2^{\prime }}^{\star }\right)\right|& -\left|{\tilde{\ddot{A}}}_{:2}^{2^{\prime }}\left({\omega }_{2^{\prime }}^{\star }\right)\right|& 0\\ \left|{\tilde{\ddot{A}}}_{:1}^{2^{\prime }}\left(2{\omega }_{2^{\prime }}^{\star }\right)\right|& \left|{\tilde{\ddot{A}}}_{:2}^{2^{\prime }}\left(2{\omega }_{2^{\prime }}^{\star }\right)\right|& 0\\ 0& 0& 0\end{array}\right],$$

(54)

$$\begin{array}{l}\ddot{\mathit{\gamma }}=\left[\begin{array}{c}{\mathit{d}}_C\\ {\ddot{\mathit{\xi }}}_{\gamma }\end{array}\right]\\ ={\left[\begin{array}{c}{d}_{C_x},d_{C_y},d_{C_z},\left|{\widetilde{\ddot{\mathsf{\xi }}}}_{Ka}^1\left({\omega }_1^{\star }\right)\right|,\left|{\widetilde{\ddot{\mathsf{\xi }}}}_K^1\left({\omega }_1^{\star }\right)\right|,\left|{\widetilde{\ddot{\mathsf{\xi }}}}_{Ka}^{1^{\text{'}}}\left({\omega }_{1^{\text{'}}}^{\star }\right)\right|,\left|{\widetilde{\ddot{\mathsf{\xi }}}}_K^{1^{\text{'}}}\left({\omega }_{1^{\text{'}}}^{\star }\right)\right|,\\ \left|{\widetilde{\ddot{\mathsf{\xi }}}}_{Ka}^2\left({\omega }_2^{\star }\right)\right|,\left|{\widetilde{\ddot{\mathsf{\xi }}}}_K^2\left({\omega }_2^{\star }\right)\right|,\left|{\widetilde{\ddot{\mathsf{\xi }}}}_{Ka}^{2^{\text{'}}}\left({\omega }_{2^{\text{'}}}^{\star }\right)\right|,\left|{\widetilde{\ddot{\mathsf{\xi }}}}_K^{2^{\text{'}}}\left({\omega }_{2^{\text{'}}}^{\star }\right)\right|\end{array}\right]}^{⊺}\end{array}$$

(55)

$$\ddot{\mathsf{\Pi }}=\left[\begin{array}{cccccccc}{\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}& 0& 0& 0& 0& 0& 0\\ {\displaystyle \frac{16}{7}}{\ddot{\mathsf{\alpha }}}_1^{\star }& -{\displaystyle \frac{9}{7}}{\ddot{\mathsf{\alpha }}}_1^{\star }& 0& 0& 0& 0& 0& 0\\ 1& -1& 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{16}{7}}{\ddot{\mathsf{\alpha }}}_{1^{\prime }}^{\star }& -{\displaystyle \frac{9}{7}}{\ddot{\mathsf{\alpha }}}_{1^{\prime }}^{\star }& 0& 0& 0& 0\\ 0& 0& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{16}{7}}{\ddot{\mathsf{\alpha }}}_2^{\star }& -{\displaystyle \frac{9}{7}}{\ddot{\mathsf{\alpha }}}_2^{\star }& 0& 0\\ 0& 0& 0& 0& 1& -1& 0& 0\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{16}{7}}& -{\displaystyle \frac{9}{7}}\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{16}{7}}{\ddot{\mathsf{\alpha }}}_{2^{\prime }}^{\star }& -{\displaystyle \frac{9}{7}}{\ddot{\mathsf{\alpha }}}_{2^{\prime }}^{\star }\\ 0& 0& 0& 0& 0& 0& 1& -1\end{array}\right],$$

(56)

$$\ddot{\mathit{\beta }}={\left[\begin{array}{c}\left|{\tilde{\ddot{b}}}^1\left({\omega }_1^{\star }\right)\right|,\left|{\tilde{\ddot{b}}}^1\left(2{\omega }_1^{\star }\right)\right|,\left|{\tilde{\ddot{\sigma }}}^1\left({\omega }_1^{\star }\right)\right|,\left|{\tilde{\ddot{b}}}^{1^{\prime }}\left({\omega }_{1^{\prime }}^{\star }\right)\right|,\left|{\tilde{\ddot{b}}}^{1^{\prime }}\left(2{\omega }_{1^{\prime }}^{\star }\right)\right|,\left|{\tilde{\ddot{\sigma }}}^{1^{\prime }}\left({\omega }_{1^{\prime }}^{\star }\right)\right|,\\ \left|{\tilde{\ddot{b}}}^2\left({\omega }_2^{\star }\right)\right|,\left|{\tilde{\ddot{b}}}^2\left(2{\omega }_2^{\star }\right)\right|,\left|{\tilde{\ddot{\sigma }}}^2\left({\omega }_2^{\star }\right)\right|,\left|{\tilde{\ddot{b}}}^{2^{\prime }}\left({\omega }_{2^{\prime }}^{\star }\right)\right|,\left|{\tilde{\ddot{b}}}^{2^{\prime }}\left(2{\omega }_{2^{\prime }}^{\star }\right)\right|,\left|{\tilde{\ddot{\sigma }}}^{2^{\prime }}\left({\omega }_{2^{\prime }}^{\star }\right)\right|\end{array}\right]}^{⊺},$$

(57)

where ${\ddot{\alpha }}_i^{\star }=\left|{\ddot{\sigma }}^k\right(2{\omega }_k^{\star }\left)\right|/\left|{\ddot{\sigma }}^k\right({\omega }_k^{\star }\left)\right|$, $\left|{\ddot{\xi }}_K^k\right(\omega \left)\right|$ and $\left|{\ddot{\xi }}_{Ka}^k\right(\omega \left)\right|$ are the harmonic components of ${\ddot{\xi }}_K^k$ and ${\ddot{\xi }}_{Ka}^k$ at $\omega$, respectively. Due to the non-negativity of the elements of ${\ddot{\mathit{\xi }}}_{\gamma }$, the following constraints are applied:

$$\begin{array}{c}\left|{\tilde{\ddot{\xi }}}_K^k\left({\omega }_k^{\star }\right)\right|>0,\left|{\tilde{\ddot{\xi }}}_{Ka}^k\left({\omega }_k^{\star }\right)\right|>0,\end{array}$$

(58)

By solving Equations (53) and (58), we obtain the multipath noise at the maneuver frequency and the AOVs $\mathit{x}$. Likewise, the AOV of GRACE-D ${\mathit{d}}_D$ can be solved using a similar method based on measurements from 28 September 2020.

6. Results and Discussion

6.1. Optimal Angular Frequency of Each Sub-Maneuver and Multipath Noise

Figure 6 is an m-by-n mesh of pixels where $m=8$ is the number of sub-maneuvers, and $n=4001$ is the length of the grid set $S$, displaying the values of $f\left({\omega }_k\right)$ for each sub-maneuver on the grid set $S$. Each $f\left({\omega }_k\right)$ specifies the color for one pixel in Figure 6. The black downward-pointing triangle markers indicate the optimal maneuver period $T^{\star }$ of each sub-maneuver, and the corresponding angular frequency is computed as ${\omega }_k^{\star }=2\pi /T^{\star }$.

After obtaining the optimal ${\omega }_k^{\star }$, we form Equation (53) after computing ${\mathit{Y}}_{:i}^{\star }$. Due to differentiation of the Fourier transform, the amplitudes of multipath noise at ${\omega }_k^{\star }$ are approximately computed by

$$\begin{array}{c}\begin{array}{c}\begin{array}{cc}\left|{\xi }_K^k\right({\omega }_k^{\star }\left)\right|& ={\displaystyle \frac{\left|{\ddot{\xi }}_K^k\left({\omega }_k^{\star }\right)\right|}{(5{\omega }_k^{\star }{)}^2}},\end{array}\end{array}\end{array}$$

(59)

$$\begin{array}{c}\begin{array}{cc}\left|{\xi }_{Ka}^k\right({\omega }_k^{\star }\left)\right|& ={\displaystyle \frac{\left|{\ddot{\xi }}_{Ka}^k\left({\omega }_k^{\star }\right)\right|}{(5{\omega }_k^{\star }{)}^2}},\end{array}\end{array}$$

(60)

where the constant 5 stands for the sampling interval in seconds.

In the multipath noise model presented in [18], the adopted multipath noise parameter is $3 \mathsf{\mu }\mathrm{m}/\mathrm{mrad}$. Given that the amplitude of the attitude angle at ${\omega }_k^{\star }$ is approximately 17 mrad, the corresponding amplitude multipath noise $\left|\xi \right(\omega \left)\right|$ is 51,000 nm. The results in

Table 3. Multipath noise for each sub-maneuver.

nm	GRACE-C (17 September 2020)				GRACE-D (28 September 2020)
	$1^{\prime }$	$1$	$2^{\prime }$	2	$3^{\prime }$	$3$	$4^{\prime }$	4
$\left|{\ddot{\xi }}_K^k\left({\omega }_k^{\star }\right)\right|$	5367.3	-	7831.1	7098.0	6124.3	10,397.5	2276.3	6468.1
$\left|{\ddot{\xi }}_{Ka}^k\left({\omega }_k^{\star }\right)\right|$	809.6	-	457.6	568.0	7570.8	8230.8	5153.2	6118.3

align with the expectations in [8], confirming that the multipath effect is smaller than the assumption by Kim [18]. Furthermore, the multipath noise of GRACE-C is much smaller than that of GRACE-D. This difference is likely due to the KBR antenna horn of GRACE-D being more affected by the space environment than that of GRACE-C, as the x-axis of the SRF of GRACE-D aligns with its flight direction, whereas the x-axis of the SRF of GRACE-C is oriented opposite to its direction of flight.

Furthermore, the estimation of the amplitude at ${\omega }_k^{\star }$ of multipath noise offers a potential approach to calibrating the multipath effect. Based on the above analysis, the amplitude of multipath noise at any arbitrary frequency can be determined by maneuvering the attitude angle of one satellite in a sinusoidal pattern at the desired frequency. This method enables the estimation of multipath noise amplitudes at various frequency points, constrained only by the satellite’s maneuvering capabilities.

For future missions requiring multipath noise calibration, if the estimated frequency points are sufficiently dense, the multipath noise can be modeled as a fitted curve derived from the amplitudes at each frequency. Additionally, the observed differences in multipath noise amplitudes between a sub-maneuver and its corresponding mirror sub-maneuver suggest that multipath noise is influenced by the attitude orientation of the calibrated satellite. By interpolating the calibrated multipath noise amplitudes across different directions, a comprehensive estimation of multipath noise can be achieved.

6.2. KBR Antenna Offset Vectors

Table 4. Antenna offset vectors for individual GRFO satellite.

mm	GRACE-C			GRACE-D
	x	y	z	x	y	z
Recalibrated	1458.2992	−0.073	−0.526	1445.1798	0.770	−0.247
VKB1B	1444.3985	−0.017	0.448	1444.4575	0.054	0.230

shows our AOC estimates of the individual GRFO and AOCs provided by the VKB1B data product [20]. The difference in the x component of the AOC for GRACE-C is approximately 1 cm, while for GRACE-D, the corresponding difference is about 1 mm. The differences in the y and z components are less than 1 mm. Additionally, the boresight difference of two AOCs of GRACE-C is 0.648 mrad, and the boresight difference of GRACE-D is 0.458 mrad.

As presented in [21], the noise of the ranging signal measured by the LRI reaches 10 ${\displaystyle \frac{\mathrm{n}\mathrm{m}}{\sqrt{\mathrm{H}\mathrm{z}}}}$ at 40 mHz and 300 ${\displaystyle \frac{\mathrm{p}\mathrm{m}}{\sqrt{\mathrm{H}\mathrm{z}}}}$ at 1 Hz, significantly below the requirement for noise in the KBR ranging signal. Additionally, the antenna offset correction for the KBR ranging signal is intended to compensate for the effect on ranging signal due to attitude jitter. In contrast, according to [21], spacecraft attitude fluctuations are compensated by active beam pointing based on differential wavefront sensing, that is, making the LRI1B ranging signal insensitive to the attitude of the twin satellites.

To validate our recalibrated AOVs, we use the LRI ranging signal as a reference or “ruler”. The residuals between the LRI ranging signal and the KBR ranging signal corrected by two different sets of AOVs represent the validity of the recalibrated AOVs.

We denote the AOCs derived from VKB1B-provided AOVs as $\Delta {\rho }_{AOC}^{VKB1B}$ and those derived from the recalibrated AOVs as $\Delta {\rho }_{AOC}^{\mathrm{r}\mathrm{e}}$. The residuals can be represented as follows:

$$\begin{array}{c}{s}^{VKB1B}={\rho }_{LRI}+\Delta {\rho }_{LLTC}-{\rho }_{KBR}-\Delta {\rho }_{KLTC}-\Delta {\rho }_{AOC}^{VKB1B}.\end{array}$$

(61)

$$\begin{array}{c}{s}^{re}={\rho }_{LRI}+\Delta {\rho }_{LLTC}-{\rho }_{KBR}-\Delta {\rho }_{KLTC}-\Delta {\rho }_{AOC}^{re}.\end{array}$$

(62)

where $s^{VKB1B}$ is the residual between the LRI range and the VKB1B-corrected KBR range, while $s^{re}$ denotes the residual between the LRI range and the KBR range corrected by our results. Figure 7 juxtaposes the residual of the LRI ranging signal and the KBR ranging signal (1) corrected by the recalibrated AOVs ($s^{re}$), (2) corrected by AOVs from VKB1B ($s^{VKB1B}$) and (3) uncorrected on the CM calibration epochs. During these dates, the variations in cone angles increase, leading to uncorrected residuals exceeding the KBR noise requirement. Notably, across all dates, the residuals corrected using our recalibrated AOVs are smaller than those corrected using the VKB1B-provided AOVs. In the mHz frequency band, where time-variable gravity field signals are concentrated, our results show significantly smaller residuals, highlighting the improved accuracy of our recalibrated AOVs in this critical range. In contrast, at low frequencies and in high-frequency regions dominated by noise, the residuals from both methods are of similar magnitude. This demonstrates that our approach provides enhanced performance in the frequency band most relevant for gravity field recovery while maintaining comparable results in the other regions.

To further validate the recalibrated AOVs, we use the ratio of the average powers of these two residuals within a milli-hertz band. The ratio ${\mathsf{\Omega }}_{f_1}^{f_2}$ is denoted as

$$\begin{array}{c}{\mathsf{\Omega }}_{f_1}^{f_2}={\displaystyle \frac{(Q^{re}{)}_{f_1}^{f_2}}{(Q^{VKB1B}{)}_{f_1}^{f_2}}}={\displaystyle \frac{{\int }_{f_1}^{f_2}{S}^{re}\left(f^{\prime }\right)d{f}^{\prime }}{{\int }_{f_1}^{f_2}{S}^{VKB1B}\left(f^{\prime }\right)d{f}^{\prime }}}.\end{array}$$

(63)

where the $f_1=1$ mHz and $f_2$ = 10 mHz, $(Q^{re}{)}_{f_1}^{f_2}$ and $(Q^{VKB1B}{)}_{f_1}^{f_2}$ represent the band powers from $f_1$ to $f_2$ of $s^{re}$ and $s^{VKB1B}$, respectively, and $S^{re}$ and $S^{VKB1B}$ represent the power spectral densities of $s^{re}$ and $s^{VKB1B}$, respectively.

Figure 8 shows the ratio in the milli-Hertz band from 1 January 2020 to 30 June 2023. The blue area indicates $(Q^{re}{)}_{1\mathrm{m}\mathrm{Hz}}^{10\mathrm{m}\mathrm{Hz}}>(Q^{VKB1B}{)}_{1\mathrm{m}\mathrm{Hz}}^{10\mathrm{m}\mathrm{Hz}}$, that is, the KBR range corrected by the AOVs from VKB1B is closer to the LRI ranging signal, while the yellow area indicates $(Q^{re}{)}_{1\mathrm{m}\mathrm{Hz}}^{10\mathrm{m}\mathrm{Hz}}<(Q^{VKB1B}{)}_{1\mathrm{m}\mathrm{Hz}}^{10\mathrm{m}\mathrm{Hz}}$, that is, the KBR ranging signal corrected by the recalibrated AOVs is closer to the LRI ranging signal. Among the 1276 days of data analyzed in this paper, 243 days (19.0%) resulted in a value greater than 1 (blue area); 745 days (58.4%) resulted in a value less than 1 (yellow area); and 251 days (19.6%) were missing (red cross). The results demonstrate the effectiveness of the recalibrated AOVs, as the number of days with ${\mathsf{\Omega }}_{1\mathrm{m}\mathrm{H}\mathrm{z}}^{10\mathrm{m}\mathrm{H}\mathrm{z}}<1$ (yellow area) is more than three times the number of days with ${\mathsf{\Omega }}_{1\mathrm{m}\mathrm{H}\mathrm{z}}^{10\mathrm{m}\mathrm{H}\mathrm{z}}>1$ (blue area). The average value and variance of $(Q^{re}{)}_{1\mathrm{m}\mathrm{Hz}}^{10\mathrm{m}\mathrm{Hz}}$ are $3.736\times {10}^{-4}$ m and $1.728\times {10}^{-4}$ m, separately, while those of $(Q^{VKB1B}{)}_{1\mathrm{m}\mathrm{Hz}}^{10\mathrm{m}\mathrm{Hz}}$ are $4.635\times {10}^{-4}$ m and $2.144\times {10}^{-4}$ m. This indicates that, for the majority of the analyzed period, the KBR ranging signal corrected by the recalibrated AOVs aligns more closely with the LRI ranging signal compared to corrections using the VKB1B-provided AOVs.

7. Conclusions

In this study, we revisited the calibration of the APCs of GRFO’s KBR system, focusing on addressing the impact of multipath noise. By employing a frequency-domain method based on the spectral characteristics of AOCs and the double-difference signals, we successfully evaluated and suppressed multipath noise using double-difference signals and recalibrated the AOVs.

Validation against the LRI ranging signal demonstrates the effectiveness of the recalibrated AOVs. The residuals between the LRI and KBR ranging signals corrected with the recalibrated AOVs consistently outperform those corrected with VKB1B-provided AOVs, especially in the milli-Hertz band. Over the 1276 days analyzed, the power ratio ${\mathsf{\Omega }}_{1\mathrm{m}\mathrm{H}\mathrm{z}}^{10\mathrm{m}\mathrm{H}\mathrm{z}}$ indicates superior performance of the recalibrated AOVs on 58.4% of the days, more than three times the proportion for VKB1B-provided AOVs. This improvement highlights the utility of our approach in mitigating multipath noise and enhancing the accuracy of inter-satellite range measurements.

Our findings contribute to a robust methodology for future gravity missions, emphasizing the importance of on-orbit calibration and the role of spectral analysis in addressing multipath noise. This work not only enhances the utility of GRFO data but also establishes a foundation for similar calibration strategies in future satellite missions.