A Systematic Approach of Inertial Sensor Calibrations for Gravity Recovery Satellites and Its Application to Taiji-1 Mission

1. Introduction

Space-borne high-precision inertial sensors (IS) or accelerometers (ACC) can provide us references of the inertial or free-falling motions in gravitational field in space, and therefore play the key role in satellite missions related to (Newtonian or relativistic) gravitational field measurements, especially for gravity recovery missions and gravitational wave detections in space. Considering the different implementation technologies, IS that based on the electrostatic suspension and servo control technologies is still of the most precise and reliable inertial payload in present days. Electrostatic IS when working in the ACC mode can measure precisely the non-gravitational forces that exert on the satellites, and had already served the series of gravity recovery missions since the beginning of this century, including the CHAMP [1], GRACE/GFO [2,3,4,5], GOCE [6] and also Taiji-1 [7,8] missions. Aided by the technology of drag-free controls, electrostatic inertial sensors could reach the unprecedented ultra-precision level ($\sim {10}^{-15}m/s^2/H{z}^{1/2}@3mHz$) of inertial motions in space, which was successfully demonstrated by the LISA PathFinder mission [9]. Such ultra-precision IS will be the key payloads of future Gravitational Wave (GW) antennas in space, including LISA [10], Taiji [11] as well as TianQin [12].

For the high-low (such as the CHAMP and Taiji-1 missions) or low-low (like the GRACE and GFO missions) satellite-to-satellite tracking gravity missions, the details of Earth geopotentials are encoded in the orbital motions or relative motions of satellites. To map out precisely the global gravity field, the high-precision and in-orbit measurements of the non-gravitational forces that perturb the satellite orbits are required. For the present-day gravity missions including CHAMP, GRACE/GFO, this is achieved by on-board electrostatic IS systems that working in the ACC mode. The designs and the working principles for these IS systems are basically the same. The IS contains a Test Mass (TM) suspended inside an electrode cage as the reference of inertial motions, and a front end electronics (FEE) unit to read out and adjust the relative motions between the TM and the cage. In the ACC mode, the compensation or actuation voltages that pushing and maintaining the TM back to its nominal position will give rise to the precision measurements of the non-gravitational perturbations on the satellite (see detailed explanations in subsection 2.2). Such ACC data is then included in the modelling of the satellite orbits and the fittings of the Earth geopotentials. For GOCE and concepts of the Next Generation Gravity Mission (NGGM), the electrostatic IS could also work in the darg-free mode, that the non-gravitational perturbations on the satellite could be measured by the FEE of the IS but compensated by pushing, with $\mu$-N thrusters, the satellite to follow the inertial motions of the TM. For both the ACC and drag-free modes, the controlled dynamics of the TM relative to the cage is determined by the combined action of the non-gravitational forces on the spacecraft, weak disturbances on the TM, and also the compensation forces from the control loop. Therefore, to accurately interpret and make use of the IS data in gravity inversions, the parameters involved in defining the characteristics of the device need to be carefully measured and calibrated. These include, generally, the scale factors of each axis that transform the control voltages imposed on the TM into the non-gravitational forces exerted by the satellites; the bias voltages or accelerations in the readouts that come from the environmental DC forces and imbalances of the FEE; and also the offset vector from the mass center of the satellite which gives out to the confusing inertial accelerations in the readouts.

In 2022, the Taiji-1 was in its final extended phase and operated in the high-low satellite-to-satellite tracking mode, and during this year radical experiments for Taiji-1 were performed including especially the monthly global gravity field recovery experiment [8] and the re-calibrations of the key measurement system (laser interferometer and IS) with satellite maneuvers after its two years operation. Taiji-1 is the first technology demonstration satellite of the “Taiji Program in Space”, which was approved by the Chinese Academy of Science (CAS) in AUG 2018. The Taiji-1 IS system has similar designs as those onboard the gravity missions CHAMP [13], GRACE/GFO [14], GOCE [15], etc, and can work in both the ACC mode and the drag-free control mode. The “Taiji Program in Space”, initiated to expert demonstration in 2008 and officially released by the CAS in 2016, is for China’s space-borne GW observatory, namely the Taiji mission [11,16,17,18]. The successful operation of Taiji-1 in 2019 and 2020 had demonstrated and confirmed the designed performances of the scientific payloads and the satellite platform, and verified the most important individual technologies of China’s space GW antenna and possible gravity missions, including high precision laser interferometers, drag-free control system, $\mu$-N thrusters, the ultra-stable and clean platform, and especially the electrostatic inertial sensor [7,19].

Based on the classical works [20,21,22,23,24,25,26,27,28,29] and especially the valuable experiences of the calibration experiments of the GRACE and GFO missions [20,21,27,28,29], a new and systematic method for the complete operating parameters calibration of electrostatic IS system for gravity missions, including the scale factors, the acceleration biases, the offset of the Center of Mass (COM), was suggested by us and approved by the Taiji science team in the end of 2021. During 2022, a set of satellite maneuvers for Taiji-1 IS calibration were conducted, including the high frequency swings (period $\sim 30s$) and rollings (period $\sim 724s$) of the satellite along certain axes. With the observational data from the IS, star trackers and precision orbit determinations during the calibration phase, the IS operating parameters are determined with high accuracies to study their possible variations and drifts during the two years operations. Compared with the original values that determined by ground-based and in-orbit experiments [19,30], variations of these parameters can be identified, which might be caused by the mechanical disturbances during the launch, changes of the mass center of the spacecraft due to consumption of cold gas, and the aging of the electronics unit in the past three years. Based on such re-calibrated parameters, we revisit the pre-processing of the IS data, and the resolution of one of the sensitive axes of Taiji-1 IS (the z-axis) is found to be slightly improved in the interested frequency band. Such a systematic approach could be applied to the ACC or IS calibrations of gravity recovery missions like the GFO, the Chinese GRACE-type mission, and the future planned Next Generation Gravity Missions. Moreover, this approach could also shed some light on the in-orbit calibrations of the ultra-precision IS for future GW space antennas, since the principle and technology have inheritance between these two generations of the electrostatic IS payloads.

This work expands as follows: In Sec. 2 we briefly introduce the Taiji-1 satellite and the IS payload. The requirement of the IS calibrations and the systematic method we adopted are described in Sec. 3. In Sec. 4, the observational data and the processing procedure are introduced, and the re-calibrated operating parameters are compared with the original used ones; furthermore, we re-process the representative data of the Taiji-1 IS and re-evaluate its performance. The conclusion of this work could be found in Sec. 5.

2. Inertial Sensor of Taiji-1

2.1. Taiji-1 satellite

According to the three step road map of the Taiji program [16,17], Taiji-1 is the first technology demonstration satellite. The Taiji-1 satellite weights about 180 kg, and the key measurement system contains the drag-free control system and the optical metrology system. Taiji-1 was launched to a circular dawn/dusk Sun-synchronous orbit, with the altitude about $600km$ and inclination angle $97.{67}^o$. The orbit has a stable Sun-facing angle, which can provide a constant power supply for the battery and also the stable temperature gradient for the platform. The orbit coordinates system is defined as follows, the $+X$-direction is along the flight direction, $+Z$ the radial direction and $+Y$ is defined by the right-hand rule.

The IS was installed at the mass center of the spacecraft (SC), having the nominal offset $\le 150\mu m$. The IS, the drag-free control unit, and the two types of micro-thrusters, that the Hall and radio frequency iron $\mu$N thrusters, together constitute the drag-free control system. The optical metrology system contains an optical bench, a high-precision phasemeter and two Na-YAG laser sources. The TM interferometer can provide the independent readout of the position of the TM in the x-axis of the IS and the optical bench interferometer serves as a reference. Both the interferometers reached the resolutions $\le 100pm/H{z}^{1/2}$ [7]. The ultra-stable and clean satellite platform has a highly stable thermal control system, which provides the $\sim \pm 1200mK$ thermal stability of the satellite environment, $\sim \pm 350mK$ in the middle cabin and about $\pm 2.6mK$ for the key measurement system [7,30]. The reaction wheel installed along the Y-direction of the satellite is the only movable unit, which is for reliability considerations. The attitude and orbit control system contains the star trackers, gyroscopes, sun sensor, magnetometer, magnetorquer, cold gas thrusters and the controller. The architecture of Taiji-1 is shown in Figure (1) and Figure (2).

2.2. Inertial sensor

The electrostatic IS system of Taiji-1 contains mainly the mechanical assembly, the FEE unit, and auxiliary subsystems like vacuum chambers, etc. The mechanical assembly is consisted of a 72 g and $4cm\times 4cm\times 1cm$ parallelepipedic TM of titanium alloy and an electrode cage made of ultra low expansion silica that encloses the TM. Both the TM and the cage are gold coated, and inside the cage there are six pairs of electrodes facing the TM side faces, see Figure (3) for illustration and the definitions of the measurement axes in the IS frame. The TM, serves as the inertial reference, is suspended electrostatically inside the cage. When operating, the position variations of the TM relative to the cage causes the changes of the capacitance between the TM side faces and electrodes, which induces signals through Wheatstone bridges that can be picked up by the FEE to give rise to the measurements of the TM position and attitude. Based on such data and the PID algorithm, the TM is servo-controlled to its nominal position by applying low frequency actuation voltages through the same electrode in the accelerometer mode or by pushing the spacecraft through $\mu N$-thrusters in the drag-free control mode. In the normal science operations of Taiji-1, the +z direction of the IS points to the flight direction, the +x the radial direction, and +y is defined by the right-hand rule.

The IS science and housekeeping data is archived and processed at the Taiji-1 data processing center of CAS in Beijing, and the data management and the detailed processing flows can be found in [31]. The level 1 science data product contains the actuation voltages on the six electrodes. With the calibrated scale factors and biases, the actuation voltages are transformed into linear accelerations and angular accelerations of the TM relative to the cage, and with the COM offset corrected these are then written into the level 2 science data product. The position sensor data together with the IS state monitoring data including temperatures of the core assembly and the FEE unit, biased reference voltages are all packaged as IS housekeeping data product.

The nominal precision level (or the acceleration noise level) of Taiji-1’s IS is $3\times {10}^{-9}m/s^2/H{z}^{1/2}@10mHz$, see Table. Table 1 for the key design requirements. The in-orbit performance was evaluated with the measurement of the y-axis, which is one of the sensitive axes pointing to the orbital normal direction. The amplitude spectrum density (ASD) of the acceleration measured by the y-axis was $\le 2\times {10}^{-9}m/s^2/H{z}^{1/2}$ [7,19], that fulfilled the design requirements. The IS couples to the space environment and satellite platform in a rather complicated way. To be more specific, the dynamics of the TM or the residual acceleration ($a_R^i\left(t\right),i=x,y,z$) of the TM relative to the platform can be written down as

$$\begin{array}{cc}\hfill a_R^i\left(t\right)& =a_{TM}^i\left(t\right)-a_{SC}^i\left(t\right)\hfill \\ \hfill & =a_{grav,TM}^i\left(t\right)+a_{para,TM}^i\left(t\right)-a_{grav,SC}^i\left(t\right)-a_{para,SC}^i\left(t\right)+a_{in,TM}^i\left(t\right)+a_c^i(t-{\tau }_c^i)\hfill \\ \hfill & =a_c^i(t-{\tau }_c^i)-a_{para,SC}^i\left(t\right)+G^{ij}\left(t\right)d_j+a_{para,TM}^i\left(t\right).\hfill \end{array}$$

(1)

In the first line, $a_{TM}^i\left(t\right)$ and $a_{SC}^i\left(t\right)$ denote the acceleration of the TM and spacecraft with respect to the local inertial frame. In the second equal, the accelerations could be expanded as the followings. $a_{grav,SC}^i\left(t\right)$ and $a_{grav,TM}^i\left(t\right)$ are the gravitational accelerations of the SC and TM respectively. $a_{para,SC}^i\left(t\right)$ and $a_{para,TM}^i\left(t\right)$ are of parasitic accelerations, where $a_{para,SC}^i\left(t\right)$ mainly comes from non-gravitational forces from solar radiation, air drag, and Earth albedo acting on the satellite and also the mechanical disturbances from the satellite platform, and $a_{para,TM}^i\left(t\right)$ are from noise forces acting directly on the TM including actuation noises, spring-like couplings, radiometric effect, magnetic couplings, thermo-noises from gold wire attached to the TM and etc. $a_{in,TM}^i\left(t\right)$ is the inertial acceleration that comes from the relative attitude variation between the TM and the satellite. $a_c^i(t-{\tau }_c^i)$ is the compensation acceleration to keep the TM to its nominal position, and ${\tau }_c$ is the delay time of the control loop. In the ACC mode $a_c^i\left(t\right)$ is the electrostatic actuation force acting on the TM, while in the drag-free control mode $a_c^i\left(t\right)$ is the compensation force by the $\mu N$-thrusters acting on the satellite. At last, in the third line, the gravitational gradients and the inertia acceleration can be summarized into the term that is proportional to the COM offset $d_i$

$$G^{ij}\left(t\right)d_j=T^{ij}\left(t\right)d_j+{\omega }^{ik}\left(t\right){\omega }_k^j\left(t\right)d_j+{\dot{\omega }}^{ij}\left(t\right)d_j,$$

(2)

where $T^{ij}$ denotes the components of the gravitational tidal matrix and ${\omega }^{ij}$ the angular velocity matrices of the satellite relative to the local inertial frame.

In this work, we consider the ACC mode, the compensation force $a_c^i\left(t\right)$ is read out in terms of the actuation voltages and is of the science data of the IS system.

$$\begin{array}{ccc}\hfill a_c^i\left(t\right)& =& a_{para,SC}^i\left(t\right)-a_{para,TM}^i\left(t\right)-G^{ij}\left(t\right)d_j+a_R^i\left(t\right),\hfill \end{array}$$

(3)

where the actuation acceleration reads

$$\begin{array}{ccc}\hfill a_c^i\left(t\right)& =& b^i+k_1^{i\alpha }{V}_{\alpha }\left(t\right)+k_2^{i\alpha \beta }{V}_{\alpha }\left(t\right)V_{\beta }\left(t\right).\hfill \end{array}$$

(4)

Here $V_{\alpha },(\alpha =x1,x2,x3,y,z1,z2)$ are actuation voltages on each electrode, $b^i$ denotes the acceleration bias, $k^{i\alpha }$ the linear scale factors that transform the voltages to accelerations, and $k^{i\alpha \beta }$ the quadratic factors. In the normal science mode of the IS, the TM is controlled to tightly follow the motions of the electrode cage or the satellite platform, that their relative motions are $\le {10}^{2}pm/H{z}^{1/2}$ in the sensitive band. This means that the residual acceleration term $a_R^i\left(t\right)$ in Eq. (3) can be ignored. The other term that can be ignored in practical use is the quadratic term $k_2^{i\alpha \beta }{V}_{\alpha }\left(t\right)V_{\beta }\left(t\right)$ in Eq. (4). Therefore, with the COM offset term $G^{ij}\left(t\right)d_j$ been corrected, and the scale factors, bias determined, the voltages data or the actuation accelerations data can give rise to the measurements of the non-gravitational forces exerted by the satellite,

$$a_{para,SC}^i\left(t\right)=k_1^{i\alpha }{V}_{\alpha }\left(t\right)+G^{ij}\left(t\right)d_j+a_{para,TM}^i\left(t\right)+b^i.$$

(5)

The parasitic acceleration noises $a_{para,TM}^i\left(t\right)$ acting on the TM determines the noise floor of the IS system.

3. Principle of IS calibration

As discussed in the previous section, to correctly interpret and make use of the IS data, one needs to carefully determine the operating parameters of the IS device. Even though some of the relevant parameters were calibrated with ground based experiments before launch or in-orbit experiments in the commissioning phase, large disturbances during the launch, consumption of working material, agings of the electric units and etc. may still cause the changes of the characteristic of the IS device. Therefore, for Taiji-1’s IS system and missions that carry the similar electrostatic IS payloads, it is necessary to calibrate the basic set of operating parameters, including the scale factors $k^{i\alpha }$, linear bias $b^i$, and the COM offset $d_i$, with the in-orbit data and regularly within the mission lifetime.

In the followings, we discuss the calibration principles of this set of parameters and the related satellite maneuver strategies, that are adopted for Taiji-1’s calibration. The key considerations here are to try to complete the IS calibrations with less satellite maneuvers and shorter calibration time durations, and also try to reduce the possible risks as much as possible.

3.1. principle of scale factors and COM offset calibrations

For electrostatic IS systems with parallelepipedic TMs such as the cases for Taiji-1, GRACE/GRACE-FO and etc, the scale factors appeared in Eq. (4) can be divided into two sets, that the linear scale factors $[k^x,k^y,k^z]$ and angular scale factors $[{\beta }^x,{\beta }^y,{\beta }^z]$, which transform the actuation voltages imposed on the electrodes into the corresponding compensation linear accelerations $a_c^i$ and angular accelerations ${\dot{\omega }}_c^i$ respectively. For Taiji-1, given the geometrical and mechanical parameters of the TM and the electrodes, the nominal values of the two sets of scale factors can be derived,

$$\begin{array}{cc}\hfill k^x& =\frac{2{ϵ}_0{S}_x{V}_P}{M{D}_x^2},k^y=\frac{2{ϵ}_0{S}_y{V}_P}{M{D}_y^2},k^z=\frac{2{ϵ}_0{S}_z{V}_P}{M{D}_z^2},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\beta }^x& =\frac{{ϵ}_0{S}_z{l}_{\varphi }{V}_P}{J_{\varphi }{D}_z^2},{\beta }^y=\frac{{ϵ}_0{S}_x{l}_{\theta }{V}_P}{J_{\theta }{D}_x^2},{\beta }^z=\frac{{ϵ}_0{S}_x{l}_{\psi }{V}_P}{J_{\psi }{D}_x^2}.\hfill \end{array}$$

(7)

Here M stands for the mass of the TM, $J_{\varphi }$, $J_{\theta }$, $J_{\psi }$ denote the mass moment of the TM along the x, y and z axes. $S_i$ is the total area of electrode surface of the ith axis, $D_i$ the nominal distance between the TM surface and the electrode, and $l_i$ denotes the force arm of the electrodes pair that control the $ith$ rotation degree of freedom, see again Figure 3. ${ϵ}_0$ stands for the vacuum permittivity, and $V_P$ the preload biased voltage. The transformation relations between the actuation voltages and compensation accelerations are shown in Table. 2. According to the designs of the Taiji-1’s IS system, we have the following useful relations in calibrating the linear scale factors and the COM offset,

$$\begin{array}{cc}\hfill {\beta }^x& =\frac{M{l}_{\varphi }}{2J_{\varphi }}{k}^z,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\beta }^y& =\frac{M{l}_{\theta }}{2J_{\theta }}{k}^x,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\beta }^z& =\frac{M{l}_{\psi }}{2J_{\psi }}{k}^x.\hfill \end{array}$$

(10)

These relations remained unchanged during the mission lifetime since they involve only the geometrical and mechanical properties of the TM and the electrode cage. The high machining accuracy ($\delta l/l\sim {10}^{-4}$) of the TM and cage structures ensures that the relations between the scale factors are sufficiently accurate. In this case, that $l_{\theta }=l_{\psi }$ and $J_{\theta }=J_{\psi }$, we have ${\beta }_y={\beta }_z$. Another important property is that, during the normal science operation of the IS in its ACC mode, the TM is controlled to tightly follow the motions of the electrode cage or the spacecraft. For Taiji-1, the position fluctuations of the TM relative to each electrode surface is $\le {10}^{2}pm/H{z}^{1/2}$ in the sensitive band. This means that the rotations of the TM and the spacecraft could be treated as precisely synchronized, that one has

$$\begin{array}{c}\hfill {\overrightarrow{\omega }}_{TM}={\overrightarrow{\omega }}_{SC},{\overrightarrow{\dot{\omega }}}_{TM}={\overrightarrow{\dot{\omega }}}_{SC}.\end{array}$$

(11)

Here, ${\overrightarrow{\omega }}_{TM}$, ${\overrightarrow{\dot{\omega }}}_{TM}$ and ${\overrightarrow{\omega }}_{SC}$, ${\overrightarrow{\dot{\omega }}}_{SC}$ denote the angular velocities and angular accelerations of the TM and the spacecraft respectively. Therefore, despite the offset between the installation orientations of the IS system and the star tracks, the measured angular velocities and accelerations of the spacecraft and the TM are interchangeable.

The rotations or attitude variations ${\overrightarrow{\omega }}_{SC}$ and ${\overrightarrow{\dot{\omega }}}_{SC}$ of the spacecraft could be independently measured by on-board star tracks. This motivates us to make use of such attitude measurements to calibrate the scale factors, which is different from the former methods based on precision orbit determination (POD) data [21,32]. One could swing the spacecraft periodically along certain axis with a rather large angular accelerations and with relative higher frequency compare to the signals band of air drags and Solar radiations, that could be clearly identified and precisely measured by the IS system. With the inputs of the angular accelerations derived by the star track measurements and the actuation voltages readout by the front end electric unit of the IS system, one can fit the angular scale factors ${\beta }^i$ based on the equations in Table. Table 2 with the least squares estimation or the Kalman filter algorithms. According to the relations between angular and linear scale factors in Eq. (8) - (), the linear scale factors can be further determined. For Taiji-1’ IS, controls along the y-axis is independent of other degrees of freedom, its actuation voltage does not involve with any rotation controls of the TM. Therefore, the linear scale factor $k^y$ for Taiji-1’ IS system can not be calibrated with this method and is left for blank in this work. Please see Figure 4 for the illustration of this calibration method.

For the COM offset calibrations, one notices that, according to Eq. (2) and (3), the periodic swing of the spacecraft will also couple to the COM offset and produce periodic linear accelerations along the axes that are perpendicular to the rotation axis due to inertial effects, see again Figure 4. According to Table 2, one can then use the common mode readouts of the actuation voltages of each axis instead of the differential mode used in the scale factor calibrations, together with the spacecraft attitude data by the star tracks or the IS readouts itself to fit and calibrate the COM offset vector. Such method had been carefully studied and employed by the GRACE and GRACE-FO team [2,20,33,34]. While, one notices that a possible interference may come from the gravity gradient signals, since the spacecraft attitude variations would also produce periodic projections of the local gravity tidal force with the same frequency. This, on the contrary, forces us to choose satellite maneuvers with small magnitude of attitude variations. In fact, for Taiji-1 and GRACE-type missions, the magnitudes of gravity gradients $\sim {10}^{-6}/s^2$. Therefore, according to Eq. (2), for COM offset ~${10}^{-3}m$, attitude variations $\delta \omega \sim {10}^{-3}rad$ will give rise to interference signals ~${10}^{-12}m/s^2$, which could be safely ignored. This again drives us to swing the spacecraft with high frequencies to produce larger angular variations or accelerations. If the above considerations are satisfied, the remained interference from the gravity gradients can be treated as a linear term due to the orbit evolutions.

For clarity, based on Eq. (3) one can re-write the observation equations for the COM calibration as

$${\overrightarrow{a}}_c\left(t\right)=\mathbf{A}\left(t\right)\overrightarrow{d}+{\overrightarrow{a}}_{I}t+{\overrightarrow{b}}_I,$$

(12)

where

$$\mathbf{A}=\left(\begin{array}{ccc}-{\omega }_y^2-{\omega }_z^2& {\omega }_x{\omega }_y-{\dot{\omega }}_z& {\omega }_x{\omega }_z+{\dot{\omega }}_y\\ {\omega }_x{\omega }_y+{\dot{\omega }}_z& -{\omega }_x^2-{\omega }_z^2& {\omega }_y{\omega }_z-{\dot{\omega }}_x\\ {\omega }_x{\omega }_z+{\dot{\omega }}_y& {\omega }_y{\omega }_z+{\dot{\omega }}_x& -{\omega }_z^2-{\omega }_x^2\end{array}\right).$$

(13)

Here, ${\overrightarrow{a}}_{I}t+{\overrightarrow{b}}_I$ are the linear term from the non-gravitational accelerations acting on the spacecraft and the gravity tidal accelerations coupled to the TM. Given the swing maneuvers discussed above, the COM offset vector $\overrightarrow{d}$ could then be fitted out.

To summarize, we suggest to calibrate the IS scale factors and COM offset with one round swing maneuver of the Taiji-1 satellite. To enhance the signal-to-noise ratio (SNR) and reduce the possible interferences, the swing maneuver should be of high frequency compared with the signal band of non-gravitational forces, and the swing amplitude should be small to reduce the interference signals from gravitational tidal forces. Also, the time span of the maneuver should be short to make sure that the linearity of the tidal force model remains accurate enough. At last but not least, the attitude maneuvers should not be driven by thrusters, since beside the disturbances caused by propulsion the misalignment of the thrusters could produce large interference signals in linear accelerations. With these considerations, the maneuvers conducted by Taiji-1 were swings of the satellite driven by the magnetic torquers along certain axis with period about $25\sim 30s$ and total time spans $<300s$. To enlarge the angular acceleration, we operated the magnetic torquers at their full powers, and the wave-trains of the satellite angular velocity were triangular waves with magnitudes about $1\times {10}^{-4}rad/s$. The swing maneuvers were conducted on 18 MAY 2022, see Figure 6 and Figure 9 for illustrations, and the data processing and fittings are discussed in the next section.

3.2. principle of IS bias calibration

From the physical point of view, the intrinsic bias $b^i$ in the actuation acceleration measurements in Eq. (5) mainly comes from the asymmetry of the electrodes on the opposite sides of the same axis and the imperfection in FEE unit. The imbalance of mass distributions surrounding the IS system, couplings between the TM and residual magnetic field and etc., may also contribute to the intrinsic bias accelerations. Therefore, the intrinsic bias along each axis is stable, and its changes could be ignored in short time measurements. On the contrary, the projections of the DC or very low frequency non-gravitational forces along each axis change not only with the orbit positions but also with the attitude of the satellite.

Generally, the long term energy loss due to orbital decays based on the POD data and the work done by the drag forces evaluated by the IS data need to be balanced, which gives us a method to find out the intrinsic biases. While, such calibration method requires rather long term and continuous observations, and also the precision data of Earth geopotentials as inputs. For related missions, to avoid these technical difficulties and to make use of the IS data in time, we suggest here to roll the satellite to give rise to a quick calibration of the intrinsic biases with only the in-orbit measurements as inputs.

According to Eq. (3) and (4), for the rolling maneuver, we re-write the actuation acceleration measurements as

$$\begin{array}{ccc}\hfill a_c^i\left(t\right)& =& b^i+k_1^{i\alpha }{V}_{\alpha }\left(t\right)\hfill \\ & =& \sum _{J=1}^3{a}_{para,SC}^J\left(t\right)cos({\Theta }^{iJ}\left(t\right)+{\Theta }_0^{iJ})-G^{ij}\left(t\right)d_j+a_{para,TM}^i,\hfill \\ \hfill \Rightarrow k_1^{i\alpha }{V}_{\alpha }\left(t\right)& =& -b^i+\sum _{J=1}^3{a}_{para,SC}^J\left(t\right)cos({\Theta }^{iJ}\left(t\right)+{\Theta }_0^{iJ})-G^{ij}\left(t\right)d_j+a_{para,TM}^i.\hfill \end{array}$$

(14)

Here $a_{para,SC}^J$ with $(J=X,Y,Z)$ denotes the components of the non-gravitational accelerations exerted by the satellite in the orbit coordinates system. ${\Theta }^{iJ}$ is the angle between the $ith$ axis of the IS system and the $Jth$ axis of the orbit coordinates system. This rolling modulation will separate the DC and low frequency non-gravitational forces from the intrinsic biases of the IS in the linear acceleration measurements, and could be subtracted or averaged out from the data to suppress their effects on bias estimations.

For practical use, this method benefits when the maneuver time span for each estimation was short, that for short orbital arcs the non-gravitational forces could be treated as varies linearly with time,

$$a_{para,SC}^J\left(t\right)=a^{J}t+a_0^J,$$

(15)

see Figure 17(b) in subSection 3.2 for illustrations. The input data sets include the angles ${\Theta }^{iJ}$, which can be derived by the POD data from GPS or Beidou system and the satellite attitude data from the star trackers, the actuation voltages that are readout by the FEE unit of the IS, and also the scale factors and COM offset been calibrated. The periodic terms $a_{para,SC}^J\left(t\right)cos({\Theta }^{iJ}\left(t\right)+{\Theta }_0^{iJ})$ on the right-hand side of Eq. (14) can be fitted and subtracted from the IS actuation accelerations. Generally, with the in-orbit mass center adjustment for the satellite platform, the COM offset term $G^{ij}\left(t\right)d_j$ could then be ignored in the data fittings. IF not, the $G^{ij}\left(t\right)d_j$ term could also be modeled with the above input data and subtracted from the IS readouts. Then, the biases can be estimated based on the above observation equation (14). For Taiji-1, to fulfill the requirement discussed above, the rolling period of the satellite was about $724s$, and to test the effectiveness of this method and also to accumulate data segments with better qualities, the entire time span of the rolling maneuver was $1.6\times {10}^{6}s$, see Figure 17 for illustrations.

To conclude this section, the complete calibration process of the scale factors, COM offset and IS biases is summarized in Figure 5.

4. Calibration results

In 2022 the final extended phase of Taiji-1, still many experiments were planned and performed. The swing maneuvers for the IS scale factors and COM offset calibrations were conducted for several trials with slightly different frequencies to obtain a more accurate triangular wave trains of the SC/TM angular velocities or square wave trains of the angular accelerations. In this work, the scale factors and COM offset calibrations are based on the data of the swing maneuver conducted on 18 MAY 2022. After the encounter of the satellite with the Earth’s shadow, the calibration experiments were continued in AUG 2022. The rolling maneuver for IS bias calibration was conducted from AUG 2022, lasted for about $1.6\times {10}^{6}s$ to accumulate enough data. The detailed data processing procedures and fitting algorithms are expended in the following subsections.

4.1. Swing maneuver and data preprocessing

As discussed in the subsection 3.1, for scale factors and COM offset calibrations, the related data products are the IS actuation voltages readouts, and the satellite attitude data from star trackers. The POD data from the Beidou or GPS system is also required to determine the position of the satellite and the local orbital coordinates system. Then, the satellite attitude variations in the local orbital coordinates system is derived first. According to the calibration principles, the Fourier components of the swing frequency and its harmonics in the data are used to fit the corresponding parameters. To reduce possible interferences from the high frequency noises (mainly comes from the FEE unit) and long time drifts, the actuation voltages and attitude data are detrended, and then smoothed with low path filters. The low-path filter used here is the CRN filter (a classical digital filter characterized by an N-th order self Convolution of Rectangular time-domain window function) with the cut-off frequency 0.1 Hz. One sees Figure (6) - (9) for the illustrations of the attitude and IS voltages data.

Figure 6. The angular velocity of the satellite measured by star tracker. The blue lines denote the data before detrend and the red lines denote the data after detrend.

Figure 6. The angular velocity of the satellite measured by star tracker. The blue lines denote the data before detrend and the red lines denote the data after detrend.

From these figures, one sees that the projections of the swings along the z-axis (IS frame) is rather obscure. This is acceptable as it is mentioned in the previous section that for Taiji-1’s IS the scale factor $k^y$ can not be calibrated by this method, and for the COM offset calibrations only swings along two different orthogonal axes are needed because of the redundancy. The outliers or spikes, for example the $V_{x1}$ data, are kept in the data, since no evident correspondence between such possible “data anomalies” with any instruments or payload events are found.

4.2. Results of scale factor calibrations

After the preparations of the IS voltages and satellite attitude data sets, base on the method discussed in subsection 3.1, the angular scale factors of Taiji-1’s IS ([${\beta }^x,{\beta }^y,{\beta }^z$] with ${\beta }^y={\beta }^z$) can be estimated by means of the following algorithms of Kalman filters

$$\begin{array}{cc}\hfill x_{n+1}& ={\Phi }_n{x}_n+{\Gamma }_n{u}_n+{{\rm Y}}_n{w}_n,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\widehat{y}}_n& =H_n{x}_n+{ϵ}_n.\hfill \end{array}$$

(17)

The above equations are the standard form of the discrete-time linear Kalman filter, and $\Phi$ is the state transition matrix, $\Gamma$ the gain of the input u in the prediction/propagation equation, and ${\rm Y}$ the gain of the noise w of the dynamics of the estimator x. The number n represents the sample or the step number and the “$\widehat{}$” labels the estimation values. Here, the estimator x, input y and model matrix H are defined as

$$\begin{array}{cc}\hfill x& ={[{\beta }^x,{\beta }^y]}^T,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill y& ={[{\dot{\omega }}_x,{\dot{\omega }}_y]}^T,\hfill \end{array}$$

(19)

$$H=\left(\begin{array}{cc}{V}_{z1}-V_{z2}& 0\\ 0& V_{x2}-V_{x3}.\end{array}\right).$$

(20)

For ${\beta }^x$ and ${\beta }^y$ are constant parameters without dynamic feature, we have

$$\begin{array}{cc}\hfill {\Phi }_k& =1,\hfill \\ \hfill {\Gamma }_k& ={{\rm Y}}_k=0.\hfill \end{array}$$

Figure 7. (a) Satellite angular velocities measured by star tracker where Red lines denote the raw readouts and blue lines denote the data filtered by CRN filter. (b) Satellite angular acceleration measured by star tracker. (c) The filtered satellite angular acceleration filtered by CRN filter.

Figure 7. (a) Satellite angular velocities measured by star tracker where Red lines denote the raw readouts and blue lines denote the data filtered by CRN filter. (b) Satellite angular acceleration measured by star tracker. (c) The filtered satellite angular acceleration filtered by CRN filter.

Figure 8. IS voltage readouts alone x-axis during the swing maneuver.

Figure 8. IS voltage readouts alone x-axis during the swing maneuver.

Figure 9. IS voltage readouts alone z-axis during the swing maneuver.

Figure 9. IS voltage readouts alone z-axis during the swing maneuver.

In this case, to minimize the measurement noise $ϵ$, the Kalman filter equations can be re-written as

$$\begin{array}{cc}\hfill {\widehat{x}}_{n+1}& =x_n+K_n({\widehat{y}}_n-H_n{\widehat{x}}_n),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill K_n& =P_n{H_n}^T{[\left(H_n{P}_n{H_n}^T\right)+R_n]}^{-1},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill P_{n+1}& =[I-K_n{H}_n]P_n,\hfill \end{array}$$

(23)

where the error covariance matrix P and variance matrix R of the measurement noise are defined as

$$\begin{array}{cc}\hfill P_n& =E\left\{{\widehat{x}}_n{\widehat{x}}_n^T\right\},\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill R_n& =E\left\{{\widehat{ϵ}}_n{\widehat{ϵ}}_n^T\right\}.\hfill \end{array}$$

(25)

With this algorithm and given the prepared data from the swing maneuver, the angular scale factors ${\beta }^x$ and ${\beta }^y$ convergent rather fast to their estimated values, see Figure 10 for illustrations.

The comparisons between the angular accelerations from the star tracks and the modeled actuation voltages with calibrated scale factors are shown in Figure 11, and their residuals are shown in Figure 12 and Figure 13.

At last, with the relations Eq.(8) and (), the complete scale factors calibrated with this method together with their estimation errors are listed in Table. Table 3.

4.3. Results of COM calibration

With the re-calibrated scale factors, the actuation voltage signals of the IS are transformed to the linear compensation acceleration $a_c^i$ and angular compensation acceleration ${\dot{\omega }}_c^i$. Compared with the star track data, the attitude variation signals of the TM (or the satellite) measured by the IS have a better SNR. Therefore, based on the discussion in the previous section, that during the swing maneuver the IS was in its normal ACC mode and the TM followed tightly the rotations of the satellite platform, we will use the attitude variation data ${\dot{\omega }}_{TM}$ measured by the IS instead of ${\dot{\omega }}_{SC}$ derived by the star tracks data in the following COM offset calibrations.

The estimation algorithm based on Kalman filters is the same as that for the scale factors estimations in the previous subsection. The COM offset vector $\overrightarrow{d}$ is an unknown constant vector without dynamical features, therefore the equations of parameter estimations can be re-written in the form of Eq. (21) - (24). Here, for the case of COM offset calibration, the estimator x, input y and model H are

$$\begin{array}{cc}\hfill x& ={[d_x,d_y,d_z]}^T,\hfill \\ \hfill y& ={[a^x,a^y,a^z]}^T,\hfill \\ \hfill H& =\mathbf{A}({\omega }_i,{\dot{\omega }}_i).\hfill \end{array}$$

The convergences of the COM offset components are shown in Figure 14. The comparisons between the measured linear acceleration (after pre-processing) and the modeled inertial accelerations with the re-calibrated values $d^i$ are shown in Figure 15. The re-calibrated COM offset vector is suggested to the Taiji-1 science team, and in practical use the inertial accelerations are to be modeled given the attitude data and COM offset values, and be subtracted from the linear acceleration data (COM correction). The ASD curves of the TM’s linear accelerations before and after the COM corrections during the swing maneuver can be found in Figure 16, and with the estimated values $d_i$ the peaks of the inertial accelerations are successfully identified and removed.

Finally, the calibrated COM offset values and their estimation errors are listed in Table 4. After nearly three years operation, the magnitudes of Taiji-1 COM offset are still within the order of ${10}^2\mu m$, while its long-term changes are also evident with respect to its nominal values.

4.4. Results of bias calibration

Based on the principle of bias calibration discussed in Sec. Section 3.2, the Taiji-1 satellite performed a long time and uniform rolling maneuver alone y-axis (IS frame) in AUG 2022, which lasted for about $1.6\times {10}^{6}s$ ($\sim 19days$). The rolling periods were short, $\sim 724s$, that a data segment about $2000\sim 3000s$ long would contain several rolling periods and the linear approximations of the non-gravitational forces in Eq. (15) could be applied. A segment of the time series data of the satellite attitude evolution $\theta \left(t\right)$ alone y-axis is shown in Figure 17(a), and the linear accelerations modulated by the rollings in the x-axis and z-axis are shown in Figure 17(b). Their ASD curves can be found in Figure 18.

Figure 17. (a). Attitude evolution of the satellite alone y-axis during the rolling maneuver. (b). The TM linear accelerations along the x and z axes during the rolling maneuver.

Figure 17. (a). Attitude evolution of the satellite alone y-axis during the rolling maneuver. (b). The TM linear accelerations along the x and z axes during the rolling maneuver.

One can see that in Figure 17(b) the amplitudes evolutions of the oscillating signals in the TM linear acceleration measurements agree well with the linear model from Eq. (15),

$$\sum _{J=1}^3{a}_{para,SC}^J\left(t\right)cos({\Theta }^{iJ}\left(t\right)+{\Theta }_0^{iJ})=\sum _{J=1}^3[a^{J}tcos({\Theta }^{iJ}\left(t\right)+{\Theta }_0^{iJ})+a_0^{J}cos({\Theta }^{iJ}\left(t\right)+{\Theta }_0^{iJ})].$$

Other long-term trends due to orbital evolutions, see the linear acceleration in x-axis in Fig Figure 17(b) as an example, can be fitted out with quadratic or cubic polynomial fitting methods.

The estimation procedure for the intrinsic biases of Taiji-1’s IS can be summarized as the following five steps: 1) preprocessing including data segmentation and quality check, 2) fitting out trends in the TM linear accelerations due to orbital evolutions, 3) fitting out the oscillating non-gravitational signals and the COM offset coupled signals, 4) estimations of the biases, 5) statistical analysis of the estimations for each data segment. In preprocessing, the entire data set of the rolling maneuver is divided into short data segments with only a few rolling periods. The data quality of each segment were checked, and we choose only the data segments that have a more uniform rolling rate. About 300 data segments with different lengths are used in the bias calibrations, and the final estimations, see Table 5, are obtained from the statistical analysis of the calibration results for each segment.

4.5. IS performance evaluation

With all the necessary operating parameters been calibrated in this work, the performance of the Taiji-1’s IS system should be re-estimated with the updated parameters. In previous works [7,19], the performance or resolution level was obtained based on the y-axis measurements, since the projections of the non-gravitational forces, such as air drags and etc., in the orbital normal direction are small in comparing with those in the flight direction. But, as discussed in the previous sections, the scale factors $k^y$ can not be calibrated with our new method, therefore the performance of the other sensitive axis, that the z-axis along the flight direction, is re-estimated and discussed here. The original estimation of the noise along this axis can be found in Figure 16(a) from [19] and in Figure 2(d) from [7], which is above the $3\times {10}^{-9}m/s^2/H{z}^{1/2}$ level from 1 mHz to 1 Hz and reach about $\sim {10}^{-8}m/s^2/H{z}^{1/2}$ in the mHz band. Since the z-axis is more noisy, the accurate estimation of the best performance level of the IS system will not be trivial, which is beyond the scope of this paper and will be left for future works.

Here, a data set of $V_{z1},V_{z2}$ in 2022 with good quality is selected for the performance evaluation. With the updated scale factors and biases, and after the COM corrections, the ASD curves of the TM linear accelerations $a_c^z$, which can be viewed as an estimation of the noise floor $a_{para,TM}^z$ in the z-axis measurement, is found in Figure 19. It is found that, compared with the previous results, the new noise floor or resolution level is improved and now below the $3\times {10}^{-9}m/s^2/H{z}^{1/2}$ level in the sensitive band, and even below the ${10}^{-9}m/s^2/H{z}^{1/2}$ level in the mHz band. This conclusion requires more analysis and further investigations, and is beyond the scope of this work.

5. Conclusion

In this work, for gravity recovery missions, we suggest a systematic approach to calibrate the most important operation parameters for space-borne electrostatic IS system with two sets of satellite maneuvers. The operating parameters considered include the scale factors and COM offset vector that can be calibrated through one swing maneuver, and the intrinsic bias accelerations that can be calibrated through rolling maneuvers. The time spans required for these satellite maneuvers depend on the accuracies of the specific IS payloads and the requirements on the precision of the parameter calibrations. The working principle of the electrostatic IS and the estimation principles for its operating parameters are discussed, and the corresponding observation equations derived.

This systematic approach was applied to the re-calibrations of Taiji-1’s IS system to study the drift and variations of the operating parameters after its two years operation. The two satellite maneuvers were conducted in MAY 2022 and AUG 2022 respectively. For Taiji-1, the short-time swing maneuver about 200 s is sufficient for scale factors and COM offset calibrations, and a data segment of a few thousand seconds long from the rolling maneuver is enough for a preliminary calibration of the biases. The real rolling maneuver lasted for about 19 days, and this is for the accumulations of data segments with better qualities and for the improvement of the accuracy of the bias estimations. The scale factors, COM offset and the intrinsic bias accelerations of Taiji-1’s IS are precisely calibrated with the in-orbit data. The linear scale factor $k^y$ of the y-axis can not be calibrated by our new method, and its updated value is left for blank in this work. The complete set of re-calibrated parameters are suggested to the Taiji-1 science team, and are archived in the Taiji-1 data center of CAS in Beijing for future use in data processing.

One of the main objectives of Taiji-1 was to evaluate the performance of the IS payload and test the related technologies. Therefore, with the IS operating parameters being updated, a re-estimation of the performance of the IS and the comparison with the former results should be carried out. As mentioned, since the new scale factor for the best sensitive y-axis is not available, we re-processed the data of the z-axis and find that, compared with former results, the performance or noise floor is improved with the updated parameters. While a complete evaluation of the best performance of Taiji-1’s IS system is not trivial, and this is beyond the scope of this paper and needs more in-depth investigations and analysis in future works. On the other hand, in the extended phase in 2022, the Taiji-1 satellite operated in the high-low satellite-to-satellite tracking mode and could provide us monthly data of global gravity field. The re-calibration of the IS could improve the accuracies of the measurements of non-gravitational forces that disturb the orbit motions of the satellite, and the Taiji-1’s global gravity model could be further updated.

At last but not least, this systematic approach could offer high reference value for the ACC or IS calibrations of gravity recovery missions like the Chinese GRACE-type mission and the future planned Next Generation Gravity Missions. Also, it could shed some light on the in-orbit calibrations of the ultra-precision IS for future GW space antennas, since the principle and technology inheritance between these two generations of the electrostatic IS payloads.