License: CC BY 4.0
arXiv:2307.01724v2 [astro-ph.IM] 09 Apr 2024

Calibration of the in-orbit center-of-mass of TaiJi-1

Xiaotong Wei weixiaotong181@mails.ucas.ac.cn University of Chinese Academy of Sciences (UCAS), Beijing 100049, China International Centre for Theoretical Physics Asia-Pacific, UCAS, Beijing 100190, China Taiji Laboratory for Gravitational Wave Universe (Beijing/Hangzhou), UCAS, Beijing 100190, China    Li Huang University of Chinese Academy of Sciences (UCAS), Beijing 100049, China International Centre for Theoretical Physics Asia-Pacific, UCAS, Beijing 100190, China Taiji Laboratory for Gravitational Wave Universe (Beijing/Hangzhou), UCAS, Beijing 100190, China    Tingyang Shen University of Chinese Academy of Sciences (UCAS), Beijing 100049, China International Centre for Theoretical Physics Asia-Pacific, UCAS, Beijing 100190, China Taiji Laboratory for Gravitational Wave Universe (Beijing/Hangzhou), UCAS, Beijing 100190, China    Zhiming Cai Innovation Academy for Microsatellites, Chinese Academy of Sciences, Shanghai 201304, China    Jibo He jibo.he@ucas.ac.cn University of Chinese Academy of Sciences (UCAS), Beijing 100049, China International Centre for Theoretical Physics Asia-Pacific, UCAS, Beijing 100190, China Taiji Laboratory for Gravitational Wave Universe (Beijing/Hangzhou), UCAS, Beijing 100190, China Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
(April 9, 2024)
Abstract

Taiji program is a space mission aiming to detect gravitational waves in the low frequency band. Taiji-1 is the first technology demonstration satellite of the Taiji Program in Space, with the gravitational reference sensor (GRS) serving as one of its key scientific payloads. For accurate accelerometer measurements, the test-mass center of the GRS must be positioned precisely at the center of gravity of the satellite to avoid measurement disturbances caused by angular acceleration and gradient. Due to installation and measurement errors, fuel consumption during in-flight phase, and other factors, the offset between the test-mass center and the center-of-mass (COM) of the satellite can be significant, degrading the measurement accuracy of the accelerometer. Therefore, the offset needs to be estimated and controlled within the required range by the center-of-mass adjustment mechanism during the satellite’s lifetime. In this paper, we present a novel method, the Extended Kalman Filter combined with Rauch-Tung-Striebel Smoother, to estimate the offset, while utilizing the chi-square test to eliminate outliers. Additionally, the nonlinear Least Squares estimation algorithm is employed as a crosscheck to estimate the offset of COM. The two methods are shown to give consistent results, with the offset estimated to be dx𝑑𝑥absentdx\approxitalic_d italic_x ≈--0.190.190.190.19 mm, dy0.64𝑑𝑦0.64dy\approx 0.64italic_d italic_y ≈ 0.64 mm, and dz𝑑𝑧absentdz\approxitalic_d italic_z ≈--0.820.820.820.82 mm. The results indicate a significant improvement on the noise level of GRS after the COM calibration, which will be of great help for the future Taiji program.

preprint: APS/123-QED

I Introduction

In 2016, the LIGO collaboration made a groundbreaking discovery by detecting gravitational waves (GW) [1], thereby providing a direct verification of the predictions made by Albert Einstein in his general theory of relativity a century ago [2]. This discovery has had a profound impact on basic scientific research worldwide.

Space-based gravitational wave detection represents an intriguing frontier for future studies of the gravitational universe, as it can extend the reach of gravitational wave astronomy beyond that of the ground-based detectors, thereby allowing for a wider range of gravitational radiation sources to be observed [3], which will provide invaluable information to deepen our understanding of the evolution of early universe and the nature of gravity. Several space-borne gravitational-wave observatories have been proposed, such as LISA [4, 5, 6], DECIGO [7], ASTROD [8], Taiji [3, 9], and Tianqin [10].

The Taiji program, initiated by the Chinese Academy of Sciences (CAS), is a space mission that aims to detect gravitational waves (GW) in the frequency band between 0.1 mHz and 1.0 Hz, which is important in the fields of astronomy and cosmology [9, 11, 12]. The Taiji program proposes to detect GW signals using the Michelson laser interferometer principle, where each end of the interferometer contains a test mass (TM) serving as a reference body. This reference body is required to be free from spurious accelerations relative to its local inertial frame, and any spurious accelerations will affect the detection of tidal deformations caused by gravitational waves.

To facilitate the development of technology for the Taiji program, a three-step road map has been proposed [13, 14]. As the first step, a technology demonstrator satellite, Taiji-1 [14, 15], was launched on August 31, 2019. One of the key technologies validated by Taiji-1 is the Gravity Reference Sensor (GRS), which serves as an accelerometer and consists of sensors and electronic components [16]. The sensor comprises an electrode housing and a TM surrounded by the sensing electrode, as shown in Fig. 1.

GRS has three axes, including one non-sensitive axis and two sensitive axes, the directions of +XX+\textrm{X}+ X, +ZZ+\textrm{Z}+ Z, and +YY+\textrm{Y}+ Y in Fig. 1 correspond to the non-sensitive axis (radial direction), the first sensitive axis (flight direction), and the second sensitive axis, respectively. The sensor utilizes capacitive sensing technology to measure the disturbance acceleration of Taiji-1. The resulting data is sent to the drag-free controller, which instructs the thruster to apply force to compensate for the disturbing force. As the reference body of the future Taiji program interferometer, GRS needs to effectively mitigate the impact of all non-gravitational accelerations so that to reach the desired sensitivity of Taiji. Moreover, the accuracy of GRS is crucial for the drag-free control system, as the GRS readout results are used as inputs for issuing commands to control the spacecraft. GRS is susceptible to a range of noise sources, including Brownian noise arising from the surrounding air near TM, charge-induced noise due to TM charge accumulation, readout noise from voltage signals, temperature gradients, circuit noise, magnetic noise, self-gravity noise.

A precise positioning of the test-mass of the accelerometer at the center-of-mass (COM) of the Taiji-1 satellite is crucial to suppress non-gravitational and perturbed accelerations such as angular motion-related accelerations and accelerations due to gravity gradients [17, 18]. The COM of the accelerometer is adjusted to be at the COM of the satellite before launching. However, during the satellite’s operation, the consumption of propellant causes the COM of the satellite to change relative to the satellite frame, leading to a shift of the accelerometer COM relative to the satellite COM over time. Therefore, it is crucial to regularly measure the deviation of the COM position of the two during the entire life cycle of the satellite, and use the COM adjustment mechanism to perform in-orbit adjustments so that the deviation is within a certain range [19, 20, 21].

In Section II, the principle of the COM calibration is presented. The accelerometer measurement model is described in Section III. We discuss the use of the Extended Kalman filter model and Rauch-Tung-Striebel Smoother for COM calibration in Section IV, while the detection and removal of outliers principle is explained in Section V. The performance of the COM calibration is evaluated in Section VI, and the results and conclusions are summarized in Section VII.

II Principle of COM calibration during in-orbit

As shown in Fig. 1, the primary components of the high-precision electrostatic levitation accelerometer include a free fall TM, a set of capacitive electrode plate that surrounds the TM (together forming a sensitive probe), and a peripheral capacitive sensing and electrostatic feedback control circuit. This circuit enables the detection of position and attitude changes between the TM and the electrodes, as well as the measurement of acceleration through its feedback voltage. There are six sensing, control, and feedback circuits that use the same principle to measure three translational accelerations and three angular accelerations of the TM concurrently.

An offset of the COM of the electrostatic accelerometer from the COM of the satellite can cause measurement disturbances, primarily attributed to the angular motion of the satellite. Therefore, an estimate of the offset can be achieved using appropriate algorithms based on the relationship between the angular motion and the measurement disturbance.

Refer to caption
Figure 1: The geometric structure of the core mechanical assembly of the GRS [22].

The on-orbit center-of-mass calibration experiment is conducted using the Taiji-1 attitude control system, which consists of the attitude sensors and the actuators, as shown in Fig. 2. The additude sensors includes star sensors, three-axis magnetometers, and gyroscopes. The actuators include cold-gas micro-thrusters, Hall-effect micro-thrusters, magnetic torquers, and momentum wheel. The X-axis magnetic torquer is mounted on the +ZZ+\textrm{Z}+ Z side panel of the satellite, the Y-axis magnetic torquer is mounted on the bottom panel, and the Z-axis magnetic torquer is mounted on the top panel. The direction of the positive magnetic moment aligns with the satellite’s three axes, as shown in the left figure of Fig. 2.

In the calibration experiment, periodic torques, generated by the magnetic torquers, is applied to the satellite, inducing a cyclic oscillation by the Earth’s magnetic field, which will cause disturbance to the accelerometer measurements and ensure its magnitude is sufficiently large to disregard other disturbance effect, such as solar pressure torque and aerodynamic torque. Then, the offset can be modeled from readout of accelerometer and star tracker. Moreover, the accuracy of GRS is crucial for the drag-free control system, as the GRS readout results are used as inputs for issuing commands to control the spacecraft.

The TM is made of titanium alloy TC4, with a magnetic susceptibility of 3.2×1063.2superscript1063.2\times 10^{-6}3.2 × 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT, therefore the effects of the coupling between the TM and the magnetic torquers is much smaller than the present GRS accuracy and can be neglected.

Refer to caption
Figure 2: The layout of the Taiji-1 attitude control system (MTQ stands for magnetic torquer).

III Accelerometer measurement Model

The accelerometer measurement output model of the relative acceleration of TM and the electrode cage for Taiji-1 is presented as follows [15, 17, 18],

aout=d¨+ω˙×d+2ω×d˙+ω×(ω×d)+ag+angsubscript𝑎out¨𝑑˙𝜔𝑑2𝜔˙𝑑𝜔𝜔𝑑subscript𝑎gsubscript𝑎nga_{\textrm{out}}=\ddot{d}+\dot{\omega}\times d+2\omega\times\dot{d}+\omega% \times(\omega\times d)+a_{\textrm{g}}+a_{\textrm{ng}}italic_a start_POSTSUBSCRIPT out end_POSTSUBSCRIPT = over¨ start_ARG italic_d end_ARG + over˙ start_ARG italic_ω end_ARG × italic_d + 2 italic_ω × over˙ start_ARG italic_d end_ARG + italic_ω × ( italic_ω × italic_d ) + italic_a start_POSTSUBSCRIPT g end_POSTSUBSCRIPT + italic_a start_POSTSUBSCRIPT ng end_POSTSUBSCRIPT (1)

where aoutsubscript𝑎outa_{\textrm{out}}italic_a start_POSTSUBSCRIPT out end_POSTSUBSCRIPT represents the theoretical measurements of the accelerometer, d𝑑ditalic_d represent the COM offset between the accelerometer and the satellite, d˙˙𝑑\dot{d}over˙ start_ARG italic_d end_ARG and d¨¨𝑑\ddot{d}over¨ start_ARG italic_d end_ARG denote the first and second partial derivation of d𝑑ditalic_d relative to time. ω𝜔\omegaitalic_ω and ω˙˙𝜔\dot{\omega}over˙ start_ARG italic_ω end_ARG are the angular velocity and the angular acceleration respectively. agsubscript𝑎ga_{\textrm{g}}italic_a start_POSTSUBSCRIPT g end_POSTSUBSCRIPT is the acceleration due to the gravitational gradient, the accelerometer measurement model ignores the perturbations due to the solid Earth tides, ocean tides, rotational deformations, the planets including the Sun and Moon, and general relativity which can be neglected as for the calibration span is short enough. Furthermore, angsubscript𝑎nga_{\textrm{ng}}italic_a start_POSTSUBSCRIPT ng end_POSTSUBSCRIPT represents the non-gravitational acceleration on the satellite, such as atmospheric drag, solar radiation pressure, and the Earth radiation pressure.

During the offset calibration period, the deviation of the COM of the TM from the COM of the satellite can be approximated as a constant offset, given the short measurement time. Therefore, aoutsubscript𝑎outa_{\textrm{out}}italic_a start_POSTSUBSCRIPT out end_POSTSUBSCRIPT can be expressed as,

aout=ω˙×d+ω×(ω×d)+ag+ang.subscript𝑎out˙𝜔𝑑𝜔𝜔𝑑subscript𝑎gsubscript𝑎nga_{\textrm{out}}=\dot{\omega}\times d+\omega\times(\omega\times d)+a_{\textrm{% g}}+a_{\textrm{ng}}.italic_a start_POSTSUBSCRIPT out end_POSTSUBSCRIPT = over˙ start_ARG italic_ω end_ARG × italic_d + italic_ω × ( italic_ω × italic_d ) + italic_a start_POSTSUBSCRIPT g end_POSTSUBSCRIPT + italic_a start_POSTSUBSCRIPT ng end_POSTSUBSCRIPT . (2)

Considering the relatively smooth change of acceleration caused by non-conservative forces and gravity gradient, it can be approximated as a linear change within a limited time span. Therefore, the calibration interval is designed to be several minutes.

The attitude control system of the Taiji-1 satellite is equipped with a star tracker and a gyroscope to respectively measure the attitude angle of the satellite relative to the inertial coordinate system and the angular velocity ω𝜔\omegaitalic_ω. Additionally, the angular accelerations ω˙˙𝜔\dot{\omega}over˙ start_ARG italic_ω end_ARG are obtained by a second-order polynomial fitting method.

Assuming negligible scale factors and misalignment errors, after substituting the ω𝜔\omegaitalic_ω and ω˙˙𝜔\dot{\omega}over˙ start_ARG italic_ω end_ARG, the accelerometer measurement model can be expressed as follows,

Aout=A~d+αt+β+An,subscript𝐴out~𝐴𝑑𝛼𝑡𝛽subscript𝐴nA_{\textrm{out}}=\widetilde{A}d+\alpha t+\beta+A_{\textrm{n}},italic_A start_POSTSUBSCRIPT out end_POSTSUBSCRIPT = over~ start_ARG italic_A end_ARG italic_d + italic_α italic_t + italic_β + italic_A start_POSTSUBSCRIPT n end_POSTSUBSCRIPT , (3)

where A~~𝐴\tilde{A}over~ start_ARG italic_A end_ARG can be expressed as

A~=[ωy2ωz2ωxωyω˙zωxωz+ω˙yωxωy+ω˙zωx2ωz2ωzωyω˙xωxωzω˙yωzωy+ω˙xωy2ωx2].~𝐴delimited-[]superscriptsubscript𝜔y2superscriptsubscript𝜔z2subscript𝜔xsubscript𝜔ysubscript˙𝜔zsubscript𝜔xsubscript𝜔zsubscript˙𝜔ysubscript𝜔xsubscript𝜔ysubscript˙𝜔zsuperscriptsubscript𝜔x2superscriptsubscript𝜔z2subscript𝜔zsubscript𝜔ysubscript˙𝜔xsubscript𝜔xsubscript𝜔zsubscript˙𝜔ysubscript𝜔zsubscript𝜔ysubscript˙𝜔xsuperscriptsubscript𝜔y2superscriptsubscript𝜔x2\tilde{A}=\left[\begin{array}[]{ccc}-\omega_{\textrm{y}}^{2}-\omega_{\textrm{z% }}^{2}&\omega_{\textrm{x}}\omega_{\textrm{y}}-\dot{\omega}_{\textrm{z}}&\omega% _{\textrm{x}}\omega_{\textrm{z}}+\dot{\omega}_{\textrm{y}}\\ \omega_{\textrm{x}}\omega_{\textrm{y}}+\dot{\omega}_{\textrm{z}}&-\omega_{% \textrm{x}}^{2}-\omega_{\textrm{z}}^{2}&\omega_{\textrm{z}}\omega_{\textrm{y}}% -\dot{\omega}_{\textrm{x}}\\ \omega_{\textrm{x}}\omega_{\textrm{z}}-\dot{\omega}_{\textrm{y}}&\omega_{% \textrm{z}}\omega_{\textrm{y}}+\dot{\omega}_{\textrm{x}}&-\omega_{\textrm{y}}^% {2}-\omega_{\textrm{x}}^{2}\end{array}\right].over~ start_ARG italic_A end_ARG = [ start_ARRAY start_ROW start_CELL - italic_ω start_POSTSUBSCRIPT y end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ω start_POSTSUBSCRIPT z end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_CELL start_CELL italic_ω start_POSTSUBSCRIPT x end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT y end_POSTSUBSCRIPT - over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT z end_POSTSUBSCRIPT end_CELL start_CELL italic_ω start_POSTSUBSCRIPT x end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT z end_POSTSUBSCRIPT + over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT y end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT x end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT y end_POSTSUBSCRIPT + over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT z end_POSTSUBSCRIPT end_CELL start_CELL - italic_ω start_POSTSUBSCRIPT x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ω start_POSTSUBSCRIPT z end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_CELL start_CELL italic_ω start_POSTSUBSCRIPT z end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT y end_POSTSUBSCRIPT - over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT x end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT x end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT z end_POSTSUBSCRIPT - over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT y end_POSTSUBSCRIPT end_CELL start_CELL italic_ω start_POSTSUBSCRIPT z end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT y end_POSTSUBSCRIPT + over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT x end_POSTSUBSCRIPT end_CELL start_CELL - italic_ω start_POSTSUBSCRIPT y end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ω start_POSTSUBSCRIPT x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_CELL end_ROW end_ARRAY ] . (4)

The corresponding axis of the spacecraft can be represented by ωisubscript𝜔i\omega_{\textrm{i}}italic_ω start_POSTSUBSCRIPT i end_POSTSUBSCRIPT and ω˙isubscript˙𝜔i\dot{\omega}_{\textrm{i}}over˙ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT i end_POSTSUBSCRIPT (i=x,y,zixyz\textrm{i}=\textrm{x},\textrm{y},\textrm{z}i = x , y , z) for angular acceleration and angular velocity, respectively. The linear slope is represented by α𝛼\alphaitalic_α and the constant bias by β𝛽\betaitalic_β. In this study, we removed the linear effect by detrending Aoutsubscript𝐴outA_{\textrm{out}}italic_A start_POSTSUBSCRIPT out end_POSTSUBSCRIPT. Therefore, the model of Aoutsubscript𝐴outA_{\textrm{out}}italic_A start_POSTSUBSCRIPT out end_POSTSUBSCRIPT utilized in this study is,

Aout=A~d+Ansubscript𝐴out~𝐴𝑑subscript𝐴nA_{\rm out}=\widetilde{A}d+A_{\textrm{n}}italic_A start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT = over~ start_ARG italic_A end_ARG italic_d + italic_A start_POSTSUBSCRIPT n end_POSTSUBSCRIPT (5)

where Ansubscript𝐴nA_{\textrm{n}}italic_A start_POSTSUBSCRIPT n end_POSTSUBSCRIPT is the measurement noise.

IV Extended Kalman filter model and Rauch-Tung-Striebel Smoother for COM calibration

Kalman filter is a high-efficiency recursive filter [23]. The filtering theory proposed by Kalman is only applicable to linear systems. An Extended Kalman Filter (EKF) was proposed in Ref. [24, 25], and can be applied to the nonlinear field.

The accelerometer measurement model used in this study is given by Eq. 5, and one can define the state variable as the COM offset and use the Kalman filter to estimate the offset. Here, the state vector is denoted as X=[dx,dy,dz]𝑋subscript𝑑xsubscript𝑑ysubscript𝑑zX=[d_{\textrm{x}},d_{\textrm{y}},d_{\textrm{z}}]italic_X = [ italic_d start_POSTSUBSCRIPT x end_POSTSUBSCRIPT , italic_d start_POSTSUBSCRIPT y end_POSTSUBSCRIPT , italic_d start_POSTSUBSCRIPT z end_POSTSUBSCRIPT ] with di˙=0˙subscript𝑑i0\dot{d_{\textrm{i}}}=0over˙ start_ARG italic_d start_POSTSUBSCRIPT i end_POSTSUBSCRIPT end_ARG = 0 for i=x,y,zixyz\textrm{i}=\textrm{x},\textrm{y},\textrm{z}i = x , y , z. The state equation, as derived in Ref. [20], is,

X^k=Φk,k1Xk1.subscript^𝑋𝑘subscriptΦ𝑘𝑘1subscript𝑋𝑘1\hat{X}_{k}=\Phi_{k,k-1}X_{k-1}.over^ start_ARG italic_X end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = roman_Φ start_POSTSUBSCRIPT italic_k , italic_k - 1 end_POSTSUBSCRIPT italic_X start_POSTSUBSCRIPT italic_k - 1 end_POSTSUBSCRIPT . (6)

Here, k𝑘kitalic_k signify the step of filter, Φk,k1subscriptΦ𝑘𝑘1\Phi_{k,k-1}roman_Φ start_POSTSUBSCRIPT italic_k , italic_k - 1 end_POSTSUBSCRIPT represents the state transition matrix from step k1𝑘1k-1italic_k - 1 to step k𝑘kitalic_k. As no manipulation is performed on the state vector, the state transition matrix can be written as Φk,k1=IsubscriptΦ𝑘𝑘1𝐼\Phi_{k,k-1}=Iroman_Φ start_POSTSUBSCRIPT italic_k , italic_k - 1 end_POSTSUBSCRIPT = italic_I, where I𝐼Iitalic_I is the identity matrix.

The output of the accelerometer can be defined as the observation equation,

Zout,𝑘=HkXk+Vk,subscript𝑍out,𝑘subscript𝐻𝑘subscript𝑋𝑘subscript𝑉𝑘Z_{\textrm{out,{k}}}=H_{k}{X}_{k}+V_{k},italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT = italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_V start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , (7)

where Hk=A~ksubscript𝐻𝑘subscript~𝐴𝑘H_{k}=\tilde{A}_{k}italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = over~ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT, and Vksubscript𝑉𝑘V_{k}italic_V start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT is the discrete measurement noise that satisfies,

E{Vk}=0,Cov{Vk}=Rk.formulae-sequence𝐸subscript𝑉𝑘0𝐶𝑜𝑣subscript𝑉𝑘subscript𝑅𝑘E\left\{V_{k}\right\}=0,Cov\left\{V_{k}\right\}=R_{k}.italic_E { italic_V start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT } = 0 , italic_C italic_o italic_v { italic_V start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT } = italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT . (8)

Here, Rksubscript𝑅𝑘R_{k}italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT denotes the variance matrix of the measurement noise, assuming that all Vksubscript𝑉𝑘V_{k}italic_V start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT are independent, unbiased, and possess finite variance which implies that Rksubscript𝑅𝑘R_{k}italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT is a diagonal matrix.

The predicted observation equation is given by,

Z^out,𝑘=HkX^k.subscript^𝑍out,𝑘subscript𝐻𝑘subscript^𝑋𝑘\hat{Z}_{\textrm{out,{k}}}=H_{k}\hat{X}_{k}.over^ start_ARG italic_Z end_ARG start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT = italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT over^ start_ARG italic_X end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT . (9)

The estimated covariance matrix is given by,

P^k=𝚽k,k1Pk1𝚽k,k1T+Qk1,subscript^𝑃𝑘subscript𝚽𝑘𝑘1subscript𝑃𝑘1superscriptsubscript𝚽𝑘𝑘1Tsubscript𝑄𝑘1\hat{P}_{k}=\bm{\Phi}_{k,k-1}P_{k-1}\bm{\Phi}_{k,k-1}^{\mathrm{T}}+Q_{k-1},over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = bold_Φ start_POSTSUBSCRIPT italic_k , italic_k - 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_k - 1 end_POSTSUBSCRIPT bold_Φ start_POSTSUBSCRIPT italic_k , italic_k - 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT + italic_Q start_POSTSUBSCRIPT italic_k - 1 end_POSTSUBSCRIPT , (10)

where Q𝑄Qitalic_Q denotes the system noise variance matrix, which is assumed to be zero in this study.

The Kalman gain is

Kk=P^kHkT(HkP^kHkT+Rk)1,subscript𝐾𝑘subscript^𝑃𝑘superscriptsubscript𝐻𝑘Tsuperscriptsubscript𝐻𝑘subscript^𝑃𝑘superscriptsubscript𝐻𝑘Tsubscript𝑅𝑘1K_{k}=\hat{P}_{k}H_{k}^{\mathrm{T}}\left(H_{k}\hat{P}_{k}H_{k}^{\textrm{T}}+R_% {k}\right)^{\mathrm{-1}},italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT ( italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT T end_POSTSUPERSCRIPT + italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (11)

and the state update value after Kalman filter is,

Xk=X^k+Kk(Zout,𝑘Z^out,𝑘).subscript𝑋𝑘subscript^𝑋𝑘subscript𝐾𝑘subscript𝑍out,𝑘subscript^𝑍out,𝑘X_{k}=\hat{X}_{k}+K_{k}\left(Z_{\textrm{out,{k}}}-\hat{Z}_{\textrm{out,{k}}}% \right).italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = over^ start_ARG italic_X end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ( italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT - over^ start_ARG italic_Z end_ARG start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT ) . (12)

Here, X^ksubscript^𝑋𝑘\hat{X}_{k}over^ start_ARG italic_X end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT represents the prior estimate value, and Xksubscript𝑋𝑘X_{k}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT is the posterior estimate.

The update of the error covariance matrix is,

Pk=(IKkHk)P^k(IKkHk)T+KkRkKkT.subscript𝑃𝑘𝐼subscript𝐾𝑘subscript𝐻𝑘subscript^𝑃𝑘superscript𝐼subscript𝐾𝑘subscript𝐻𝑘Tsubscript𝐾𝑘subscript𝑅𝑘superscriptsubscript𝐾𝑘TP_{k}=\left(I-K_{k}H_{k}\right)\hat{P}_{k}\left(I-K_{k}H_{k}\right)^{\mathrm{T% }}+K_{k}R_{k}K_{k}^{\mathrm{T}}.italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = ( italic_I - italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ( italic_I - italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT + italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT . (13)

After the Kalman filter is applied, the difference between the predicted and observed measurements is known as the filtered residual. The filtered residual can be used to evaluate the accuracy of the filter. The covariance matrix of the filtered residual provides information on the uncertainty of the estimation. Therefore, it is important to analyze both the filtered residuals rksubscript𝑟𝑘r_{k}italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT and the covariance matrix of the filtered residuals RkKsuperscriptsubscript𝑅𝑘KR_{k}^{\textrm{K}}italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT K end_POSTSUPERSCRIPT, as follows,

rk=Zout,𝑘HkXk,subscript𝑟𝑘subscript𝑍out,𝑘subscript𝐻𝑘subscript𝑋𝑘r_{k}=Z_{\textrm{out,{k}}}-H_{k}X_{k},italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT - italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , (14)
RkK=RkHkPkHkT.superscriptsubscript𝑅𝑘Ksubscript𝑅𝑘subscript𝐻𝑘subscript𝑃𝑘superscriptsubscript𝐻𝑘TR_{k}^{\textrm{K}}=R_{k}-H_{k}P_{k}H_{k}^{\mathrm{T}}.italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT K end_POSTSUPERSCRIPT = italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT - italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT . (15)

The characteristic property of the Kalman filter is that the filtered residual vectors are uncorrelated, and even independent in the case of Gaussian distribution.

After applying the aforementioned Kalman Filter, the state estimator of the dynamic system can be further refined using the Rauch-Tung-Striebel (RTS) smoother, as described in Ref. [26]. The smoothing equations are given by the following formulas,

X^k+1=Φk+1,kXk,subscript^𝑋𝑘1subscriptΦ𝑘1𝑘subscript𝑋𝑘\hat{X}_{k+1}=\Phi_{k+1,k}X_{k},over^ start_ARG italic_X end_ARG start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT = roman_Φ start_POSTSUBSCRIPT italic_k + 1 , italic_k end_POSTSUBSCRIPT italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , (16)
P^k+1=Φk+1,kPkΦk+1,kT+Qk,subscript^𝑃𝑘1subscriptΦ𝑘1𝑘subscript𝑃𝑘superscriptsubscriptΦ𝑘1𝑘Tsubscript𝑄𝑘\hat{P}_{k+1}=\Phi_{k+1,k}P_{k}\Phi_{k+1,k}^{\mathrm{T}}+Q_{k},over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT = roman_Φ start_POSTSUBSCRIPT italic_k + 1 , italic_k end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT roman_Φ start_POSTSUBSCRIPT italic_k + 1 , italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT + italic_Q start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , (17)
Gk=PkΦk+1,kTP^k+11,subscript𝐺𝑘subscript𝑃𝑘superscriptsubscriptΦ𝑘1𝑘Tsuperscriptsubscript^𝑃𝑘11G_{k}=P_{k}\Phi_{k+1,k}^{\mathrm{T}}\hat{P}_{k+1}^{\mathrm{-1}},italic_G start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT roman_Φ start_POSTSUBSCRIPT italic_k + 1 , italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (18)
XkS=Xk+Gk(Xk+1SX^k+1),superscriptsubscript𝑋𝑘Ssubscript𝑋𝑘subscript𝐺𝑘superscriptsubscript𝑋𝑘1Ssubscript^𝑋𝑘1X_{k}^{\textrm{S}}=X_{k}+G_{k}\left(X_{k+1}^{\textrm{S}}-\hat{X}_{k+1}\right),italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT = italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_G start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ( italic_X start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT - over^ start_ARG italic_X end_ARG start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ) , (19)
PkS=Pk+Gk(Pk+1SP^k+1)GkT,superscriptsubscript𝑃𝑘Ssubscript𝑃𝑘subscript𝐺𝑘superscriptsubscript𝑃𝑘1Ssubscript^𝑃𝑘1superscriptsubscript𝐺𝑘TP_{k}^{\textrm{S}}=P_{k}+G_{k}\left(P_{k+1}^{\textrm{S}}-\hat{P}_{k+1}\right)G% _{k}^{\mathrm{T}},italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT = italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_G start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ( italic_P start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT - over^ start_ARG italic_P end_ARG start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ) italic_G start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT , (20)
rkS=Zout,𝑘HkXkS,superscriptsubscript𝑟𝑘Ssubscript𝑍out,𝑘subscript𝐻𝑘superscriptsubscript𝑋𝑘Sr_{k}^{\textrm{S}}=Z_{\textrm{out,{k}}}-H_{k}X_{k}^{\textrm{S}},italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT = italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT - italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT , (21)
RkS=RkHkPkSHkT.superscriptsubscript𝑅𝑘Ssubscript𝑅𝑘subscript𝐻𝑘superscriptsubscript𝑃𝑘Ssuperscriptsubscript𝐻𝑘TR_{k}^{\textrm{S}}=R_{k}-H_{k}P_{k}^{\textrm{S}}H_{k}^{\mathrm{T}}.italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT = italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT - italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT . (22)

The RTS smoother provides smoothed estimates of the state mean and state covariance at time step k𝑘kitalic_k, denoted as XkSsuperscriptsubscript𝑋𝑘SX_{k}^{\textrm{S}}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT and PkSsuperscriptsubscript𝑃𝑘SP_{k}^{\textrm{S}}italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT, respectively. The smoother gain on time step k𝑘kitalic_k, denoted as Gksubscript𝐺𝑘G_{k}italic_G start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT, corrects the RTS smoother estimate. The recursion is initialized at the last time step T of Kalman filter with XTS=XTsuperscriptsubscript𝑋TSsubscript𝑋TX_{\textrm{T}}^{\textrm{S}}=X_{\textrm{T}}italic_X start_POSTSUBSCRIPT T end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT = italic_X start_POSTSUBSCRIPT T end_POSTSUBSCRIPT and PTS=PTsuperscriptsubscript𝑃TSsubscript𝑃TP_{\textrm{T}}^{\textrm{S}}=P_{\textrm{T}}italic_P start_POSTSUBSCRIPT T end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT = italic_P start_POSTSUBSCRIPT T end_POSTSUBSCRIPT. The smoothed residuals are denoted as rkSsuperscriptsubscript𝑟𝑘Sr_{k}^{\textrm{S}}italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT, and the covariance matrix of smoothed residuals is denoted as RkSsuperscriptsubscript𝑅𝑘SR_{k}^{\textrm{S}}italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT.

To obtain a more accurate result, a combined method based on Kalman filter and the RTS Smoother (KF-RTS) and outliers removal is proposed. The data is filtered with the Kalman filter, and smoothed by the RTS smoother, which takes into account all the available valid data points. A chi-square confidence test is performed during the smoothing process for assessing the data quality and removing outliers. This KF-RTS process is iterated until there are no outliers.

A possible drawback of the filter algorithm is the fact that one needs an initial value of the state vector together with its covariance matrix. This can be obtained by fitting a small number of measurements at the start of the track by a conventional least-squares fit, but this is not an elegant solution. The other possibility is to start with an arbitrary state vector and an infinite covariance matrix, i.e. a large multiple of the identity matrix. This is completely in the spirit of the filtering approach, but may lead to numerical instabilities in the computation of the gain matrix, since the infinities have to cancel in order to give a finite gain matrix. This may be difficult on a computer with a short word length.

In this article the initial value of the state vector is set to zero and its covariance is chosen to be 0.001, which is large enough. The measurement noise is calculated by selecting a segment of stationary data from the corresponding data.

Here, we use the nonlinear least squares (NLLS) method [27] for parameter estimation, as a crosscheck of the KF-RTS method. As is well-known, the key challenge for NLLS is to find the value of θ^^𝜃\hat{\theta}over^ start_ARG italic_θ end_ARG that minimizes the function F(θ)𝐹𝜃F(\theta)italic_F ( italic_θ )

F(θ)12i=1m(fi(θ))2𝐹𝜃12superscriptsubscript𝑖1𝑚superscriptsubscript𝑓𝑖𝜃2F(\mathbf{\theta})\equiv\frac{1}{2}\sum_{i=1}^{m}\left(f_{i}(\mathbf{\theta})% \right)^{\mathrm{2}}italic_F ( italic_θ ) ≡ divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_θ ) ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT (23)

where fi(θ)yiModel(θ,input)subscript𝑓𝑖𝜃subscript𝑦𝑖Model𝜃inputf_{i}(\theta)\equiv y_{i}-\mathrm{Model}(\theta,\mathrm{input})italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_θ ) ≡ italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - roman_Model ( italic_θ , roman_input ). The Levenberg-Marquardt (LM) algorithm [28, 29], which is an efficient method to minimize F(θ)𝐹𝜃F(\theta)italic_F ( italic_θ ), is used to find the optimal parameters in this article.

V Detection and removal of outliers

During the Kalman filtering and RTS smoothing, each data point was utilized to obtain an optimal estimate of the state value, and this process also allows for an assessment of the quality of each data point.

The chi-square value [30] is commonly used for assessing the quality of data and detecting outliers, which can be caused by spacecraft maneuvering, non-Gaussian noise, electronic noise, or other coupled noise, and deviate significantly from the normal sequence of measurements.

The residuals of the global fit can be utilized to identify measurements with large residuals as potential outliers. Kalman filters and smoother enable the exploitation of complete information locally, to determine the validity of a measurement with high probability. The smoothed residual chi-square value can serve as a useful decision criterion for the data quality of the measurement,

χ2=(rkS)T(RkS)1rkS.superscript𝜒2superscriptsuperscriptsubscript𝑟𝑘STsuperscriptsuperscriptsubscript𝑅𝑘S1superscriptsubscript𝑟𝑘S\chi^{2}=(r_{k}^{\textrm{S}})^{\mathrm{T}}(R_{k}^{\textrm{S}})^{\mathrm{-1}}r_% {k}^{\textrm{S}}.italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = ( italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT ( italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT . (24)

It is demonstrated that the test on smoothed residual chi-square consistently outperforms the test on filtered residual chi-square [31]. Thus, searching for possible outliers should be carried out during smoothing, as it allows the utilization of complete parameter information.

During the smoothing process, the measurement point Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT can be removed from the smoothed estimate XkSsuperscriptsubscript𝑋𝑘SX_{k}^{\textrm{S}}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT to obtain XkS*superscriptsuperscriptsubscript𝑋𝑘S{X_{k}^{\textrm{S}}}^{*}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT, which represents the optimal estimate of the system state at step k𝑘kitalic_k using all data information except Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT. This optimal estimate can be utilized for the detection and removal of outliers. To remove Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT from the estimate XkSsuperscriptsubscript𝑋𝑘SX_{k}^{\textrm{S}}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT, an inverse Kalman filter can be applied with the covariance matrix of Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT taken as negative. This step of the filter is described in Ref. [31], and the smoothed estimate of Xksubscript𝑋𝑘X_{k}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT without using Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT, XkS*superscriptsuperscriptsubscript𝑋𝑘S{X_{k}^{\textrm{S}}}^{*}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT, can be calculated as,

XkS*=XkS+KkS*(Aout,𝑘HkXkS),superscriptsuperscriptsubscript𝑋𝑘Ssuperscriptsubscript𝑋𝑘Ssuperscriptsuperscriptsubscript𝐾𝑘Ssubscript𝐴out,𝑘subscript𝐻𝑘superscriptsubscript𝑋𝑘S{X_{k}^{\textrm{S}}}^{*}={X}_{k}^{\textrm{S}}+{K_{k}^{\textrm{S}}}^{*}\left(A_% {\textrm{out,{k}}}-{H}_{k}{X}_{k}^{\textrm{S}}\right),italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT = italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT + italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ( italic_A start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT - italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT ) , (25)

in which

KkS*=PkSHkT(HkPkSHkTRk)1,superscriptsuperscriptsubscript𝐾𝑘Ssuperscriptsubscript𝑃𝑘Ssuperscriptsubscript𝐻𝑘Tsuperscriptsubscript𝐻𝑘superscriptsubscript𝑃𝑘Ssuperscriptsubscript𝐻𝑘Tsubscript𝑅𝑘1{K_{k}^{\textrm{S}}}^{*}={P}_{k}^{\textrm{S}}{H_{k}}^{\mathrm{T}}\left({H}_{k}% {P}_{k}^{\textrm{S}}{H_{k}}^{\mathrm{T}}-{R}_{k}\right)^{\mathrm{-1}},italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT = italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT ( italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT - italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (26)

and

PkS*=(IKkS*Hk)PkS.superscriptsuperscriptsubscript𝑃𝑘S𝐼superscriptsuperscriptsubscript𝐾𝑘Ssubscript𝐻𝑘superscriptsubscript𝑃𝑘S{P_{k}^{\textrm{S}}}^{*}=\left({I}-{K_{k}^{\textrm{S}}}^{*}{H}_{k}\right){P}_{% k}^{\textrm{S}}.italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT = ( italic_I - italic_K start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT italic_H start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT . (27)

If Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT is a valid measurement and the covariance matrix of its Gaussian readout error is known, the quantity χ2superscript𝜒2\chi^{2}italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT follows a chi-square distribution with Nzsubscript𝑁zN_{\textrm{z}}italic_N start_POSTSUBSCRIPT z end_POSTSUBSCRIPT degrees of freedom, where Nzsubscript𝑁zN_{\textrm{z}}italic_N start_POSTSUBSCRIPT z end_POSTSUBSCRIPT is the dimension of Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT.

The measurement can be identified as an outlier if the value of χ2superscript𝜒2\chi^{2}italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT exceeds a certain threshold c𝑐citalic_c. This threshold is chosen as the (1γ)1𝛾(1-\gamma)( 1 - italic_γ ) quantile of the corresponding χ2superscript𝜒2\chi^{2}italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT distribution, where γ𝛾\gammaitalic_γ represents the probability of rejecting a valid measurement and is chosen to be γ=0.001𝛾0.001\gamma=0.001italic_γ = 0.001 in this paper. Other choices of γ𝛾\gammaitalic_γ value, 0.005, 0.01, 0.05 and 0.1 have been tried, and the effects on the estimated COM offset is found to be negligible.

The measurement Zout,𝑘subscript𝑍out,𝑘Z_{\textrm{out,{k}}}italic_Z start_POSTSUBSCRIPT out, italic_k end_POSTSUBSCRIPT can be removed permanently from the list as an outlier, and the RTS smoother can continue with XkS*superscriptsuperscriptsubscript𝑋𝑘S{X_{k}^{\textrm{S}}}^{*}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT and PkS*superscriptsuperscriptsubscript𝑃𝑘S{P_{k}^{\textrm{S}}}^{*}italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT instead of XkSsuperscriptsubscript𝑋𝑘S{X_{k}^{\textrm{S}}}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT and PkSsuperscriptsubscript𝑃𝑘S{P_{k}^{\textrm{S}}}italic_P start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT for updating the estimates XjSsuperscriptsubscript𝑋𝑗SX_{j}^{\textrm{S}}italic_X start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT start_POSTSUPERSCRIPT S end_POSTSUPERSCRIPT when j>k𝑗𝑘j>kitalic_j > italic_k. To remove all outliers, this Kalman filter and RTS smoother must be recomputed without the outliers and iterated until convergence is achieved.

VI Results of COM calibration

The COM bias estimation in this paper is performed using the two algorithms discussed earlier. Figure 3 presents a comparison between the linear acceleration results obtained from the NLLS estimation and the original data, while Fig. 4 illustrates the comparison between the linear acceleration results obtained from the extended Kalman filter algorithm and the original data. Both methods exhibit excellent agreement with the original data. The results of COM calibration using the KF-RTS algorithm before removing outliers are shown in Fig. 5, with the shaded area indicating the range of one standard deviation.

Figure 6 illustrates the comparison between the fit results obtained using the NLLS algorithm and the experimental data with outliers removed. The final round results of the offset obtained using the KF-RTS algorithm are presented in Fig. 7, with the shaded area indicating the range of one standard deviation. Figure 8 displays the comparison between the final round results obtained using the KF-RTS algorithm and the experimental data.

Refer to caption
Figure 3: The calibration experiment data and the NLLS fit results (with outliers).
Refer to caption
Figure 4: The calibration experiment data, together with the first round of outliers detected through the χ2superscript𝜒2\chi^{2}italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT test and KF-RTS smoother results obtained before removing the outliers.
Refer to caption
Figure 5: The first KF-RTS smoother result of dxsubscript𝑑xd_{\textrm{x}}italic_d start_POSTSUBSCRIPT x end_POSTSUBSCRIPT dysubscript𝑑yd_{\textrm{y}}italic_d start_POSTSUBSCRIPT y end_POSTSUBSCRIPT dzsubscript𝑑zd_{\textrm{z}}italic_d start_POSTSUBSCRIPT z end_POSTSUBSCRIPT.
Refer to caption
Figure 6: The calibration experiment data and the NLLS fit results obtained after removing outliers.
Refer to caption
Figure 7: The final KF-RTS smoother results of dxsubscript𝑑xd_{\textrm{x}}italic_d start_POSTSUBSCRIPT x end_POSTSUBSCRIPT dysubscript𝑑yd_{\textrm{y}}italic_d start_POSTSUBSCRIPT y end_POSTSUBSCRIPT dzsubscript𝑑zd_{\textrm{z}}italic_d start_POSTSUBSCRIPT z end_POSTSUBSCRIPT.
Refer to caption
Figure 8: The final round outliers detected through χ2superscript𝜒2\chi^{2}italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT test and KF-RTS smoother result and calibration experiment data.

Table 1 presents the estimation results for the COM offset obtained using the two methods. Despite the use of different estimation algorithms, the results obtained from both methods are highly consistent, thus increasing the confidence in the accuracy of the results. It is worth mentioning that, the value of χ2/nofsuperscript𝜒2nof\chi^{2}/\mathrm{nof}italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / roman_nof as the goodness of a fit is 2.33 for the first round and is 0.88 for the final round, where nofnof\mathrm{nof}roman_nof represents the number of degrees of freedom.

Table 1: The results of estimation of the offset of COM.
Axis X (μμ\upmuroman_μm) Y (μμ\upmuroman_μm) Z (μμ\upmuroman_μm)
NLLS with outliers --175 ±plus-or-minus\pm± 41 610 ±plus-or-minus\pm± 101 --777 ±plus-or-minus\pm± 31
First KF-RTS --173 ±plus-or-minus\pm± 27 606 ±plus-or-minus\pm± 66 --777 ±plus-or-minus\pm± 20
NLLS without outliers --191 ±plus-or-minus\pm± 27 643 ±plus-or-minus\pm± 66 --818 ±plus-or-minus\pm± 20
Final KF-RTS --189 ±plus-or-minus\pm± 27 638 ±plus-or-minus\pm± 68 --818 ±plus-or-minus\pm± 20

Figure 9 and Figure 10 depict the comparison of the amplitude spectral density (ASD) between the calibration experimental data and flight data before and after the COM calibration. Figure 9 shows that the peak observed in the experimental data is caused by a modulation signal, and its influence is highly suppressed after the COM calibration, which is consistent with our expectations. Figure 10 demonstrates a significant reduction in the acceleration noise level of GRS readings in the frequency range of 0.001–0.1 Hz after the COM calibration.

Refer to caption
Figure 9: The ASD of the calibration experiment data before and after the COM calibration.
Refer to caption
Figure 10: The ASD of the flight data before and after the COM calibration.

VII Discussion and Conclusions

Detecting and reducing the deviation between the COM of the inspection load and the COM of the satellite is crucial for high-precision accelerometers, as it can improve their accuracy. It is also a significant step in the development of space-based gravitational wave observatory to achieve their scientific objectives. Furthermore, the calibration can help obtain a high-precision gravity field, which is valuable for conducting geoscience research with greater accuracy.

In this study, the offset of the COM between the inspection load and the satellite is estimated using the KF-RTS smoother. Outliers are detected using the chi-square test, and the inverse Kalman filter is applied to remove them. The LM algorithm, as a cross check, is used to find the optimal offset parameters for the NLLS method. The results obtained with both methods are in very good agreement, and the offset of COM between the inspection load and the satellite is estimated with an accuracy of 𝒪(10μm)𝒪10μm\mathcal{O}(10\,\upmu\mathrm{m})caligraphic_O ( 10 roman_μ roman_m ). After obtaining the COM offset, one can reduce it using the COM adjustment mechanism, and the effects of the COM offset can also be suppressed in the data-processing. The COM calibration is crucial for improving the accuracy of the accelerometer which will directly impact the detection sensitivity of the final space-based gravitational wave observatory.

For the space-borne gravitational-wave observatories, such as Taiji-3, there are three satellites and each satellite is equipped with two TMs, which means that the COM of TMs does not coincide with that of the satellite. However, as a reference body for the satellite, it would cause residual acceleration if the TMs are away from their nominal positions. Therefore, it is necessary to periodically estimate or monitor the deviation of the COM of TM from a fixed point and make adjustments when necessary. The same calibration principle and methods as discussed in this paper can be used. The Taiji-3 satellite will be equipped with higher-precision star-sensitive instruments, which will enable more accurate results.

Acknowledgements

This work is partially supported by the Strategic Priority Research Program of the Chinese Academy of Sciences Grant No. XDA15020700 and XDA15021100, and by the Fundamental Research Funds for the Central Universities. We acknowledge support from the National Space Science Data Center, National Science and Technology Infrastructure of China.

References