-
地球重力场和海洋环流探测器(gravity field and steady-state ocean circulation explorer, GOCE)于2009年3月17日发射成功，其采用高低卫星跟踪卫星测量技术和重力梯度测量技术相结合的模式，目的是在100 km的空间分辨率下确定精度为1~2 cm的全球大地水准面和精度为1 m Gal的全球重力异常[1]。为实现这一目标，GOCE卫星搭载了高精度重力梯度仪(electrostatic gravity gradiometry, EGG)、星象仪(star tracker, STR)以及高精度全球导航卫星系统(global navigation satellite system, GNSS)星载接收机。高精度的重力梯度测量数据预处理可为恢复高精度静态地球重力场模型提供数据质量保障[2]，既是卫星重力梯度测量数据处理及应用中的重要环节，也是实现GOCE卫星预期科学目标的关键任务之一[3]。在数据预处理过程中，梯度仪部分主要是对原始数据进行计算得到标称数据，并对标称数据进行校正[4]，星象仪部分则主要是对姿态数据进行校正计算。星象仪观测的姿态数据经校正后可计算出精确的卫星姿态角速度，综合梯度仪的观测数据，可解算得到GOCE卫星level 1b(L1b)中高精度的重力梯度观测数据[5]。
GOCE卫星重力梯度level 1a-level 1b(L1aL1b)数据预处理主要包括观测数据的标定和异常数据的探测及剔除，最大限度地利用所有工程改正和标定的信息，消除测量系统本身和空间观测环境的非理想性对数据的干扰。在这一过程中，通过对GOCE卫星有效载荷的原始数据进行工程化处理以获得用于重力场恢复的科学观测数据[6]。在梯度仪数据处理校正中，Siemes等[7]通过将梯度仪校准模型扩展一个二次因子来实现，该因子的校正可大大减少主要在地球磁极附近发生的扰动，并且可以通过结合星象仪姿态数据对加速度计进行校正，进一步提高交轨方向上GOCE卫星重力梯度的精度[8-9]。后续精确计算梯度过程中，Stummer等[10]提出了一种确定角速度的新方法，基于GOCE卫星搭载的星象仪和梯度仪特性，将两个荷载观测的角速度通过维纳滤波，选择最佳相对权重在频谱域中进行组合，其解算得到的重力梯度观测值在测量带宽内显示出较为明显的改善。精确的星象仪姿态数据是重力场恢复所需的关键观测数据之一。在星象仪的数据处理上，Siemes等[11]对GOCE星象仪数据的处理上加入了温度校正，该方法以温度的线性函数来估算星象仪与梯度仪之间的相对姿态偏移，在由星象仪旋转到梯度仪坐标系的偏角中，其与星象仪温度变化相关，约为7″~9″的小偏移。文献[12-13]在对重力恢复与气候实验卫星(gravity recovery and climate experiment, GRACE)的星象仪数据处理中，使用了基于权重矩阵的组合方法，实现了对星象仪3个轴向上的角速度恢复精度提高，可在整个功率谱密度范围内实现显著改善。
目前对单个星象仪姿态数据计算得到的角速度含有较大噪声，当其从星象仪坐标系(star sensor reference frame, SSRF)转换到梯度仪坐标系(gradiometer reference frame, GRF)时，低灵敏度轴(less sensitive, LS)的角速度分量会将噪声传播到超灵敏度轴(ultra sensitive, US)的角速度分量中，进而将角速度噪声传播到后续重力梯度观测数据中。本文在传统单星象仪姿态数据处理基础上，系统分析星象仪噪声特性，基于星象仪噪声水平构建噪声分布加权矩阵，对两个或多个星象仪的姿态四元数联合求解最佳姿态四元数，并分析评估了所提方法的精度与可靠性。
-
GOCE卫星重力梯度L1b处理流程中[14-15]控制电压与加速度计内部检测质量的加速度成比例[16-17]。每个加速计将控制电压转换为线性加速度，并将其旋转到一个共同的参考系，即梯度仪参考系[18]。图 1显示了GRF内梯度计6个加速度计A1~A6的分布。
图 1 GRF内梯度计的6个加速度计位置
Figure 1. Positions of Six Accelerometers of Gradiometer in GRF
根据卫星重力梯度测量原理[19]，梯度仪中加速度计测得加速度值为：
$$ {\mathit{\boldsymbol{a}}_\mathit{\boldsymbol{i}}} = - \left( {\mathit{\boldsymbol{V}} - {\mathit{\boldsymbol{\omega }}^2} - \mathit{\boldsymbol{\dot \omega }}} \right){\mathit{\boldsymbol{r}}_\mathit{\boldsymbol{i}}} + \mathit{\boldsymbol{d}} $$ (1) 式中，下标i∈{1, 2…6}分别代表图 1中加速度计A1~A6; ai为第i个加速度计的真实加速度；V表示重力梯度张量；ω2为卫星角速度平方；$\mathit{\boldsymbol{\dot \omega }}$为卫星的角加速度；ri表示第i个加速度计与卫星质心之间的距离矢量；d为非保守力加速度[20]。重力梯度张量V除含有重力梯度外，还包括由于作用在卫星上的离心力而引起的加速度，因此要从差模加速度中减去离心加速度以分离出重力梯度[21]，计算式为：
$$ V_{x x}=-2 \frac{a_{d, 14, x}}{L_{x}}-\omega_{y}^{2}-\omega_{z}^{2} $$ (2) $$ V_{y y}=-2 \frac{a_{d, 25, y}}{L_{y}}-\omega_{x}^{2}-\omega_{z}^{2} $$ (3) $$ V_{z z}=-2 \frac{a_{d, 36, z}}{L_{z}}-\omega_{x}^{2}-\omega_{y}^{2} $$ (4) $$ V_{x y}=-\frac{a_{d, 14, y}}{L_{x}}-\frac{a_{d, 25, x}}{L_{y}}+\omega_{x} \omega_{y} $$ (5) $$ V_{x z}=-\frac{a_{d, 14, z}}{L_{x}}-\frac{a_{d, 36, x}}{L_{z}}+\omega_{x} \omega_{z} $$ (6) $$ V_{y z}=-\frac{a_{d, 25, z}}{L_{z}}-\frac{a_{d, 36, y}}{L_{y}}+\omega_{y} \omega_{z} $$ (7) 式中，V为梯度仪参考系中的重力梯度，下标x、y、z为梯度仪参考系内的轴；ad, ij, x、ad, ij, y、ad, ij, z分别为x、y、z轴的差模加速度对，数字组合ij∈{14, 25, 36}代表两两组合的加速度计对；Lx、Ly、Lz为两两加速度计对在x、y、z轴方向上的距离；ωx、ωy、ωz为梯度仪在绕x、y、z轴向上的角速度。精确求算ωx、ωy、ωz是通过星象仪四元数联合重构的角速度来实现，并合并星象仪和梯度仪角速度数据[22]。由于星象仪姿态数据相对于梯度仪数据独立，在后续的计算过程中对星象仪姿态数据进行联合校正后，单独计算分析了联合后的星象仪数据对重力梯度的影响。
进行星象仪姿态数据联合前，需要对3个星象仪STR1、STR2、STR3进行分析，星象仪确定视轴相对于惯性坐标系(inertial reference frame, IRF)的指向，即观测数据为惯性坐标系下的姿态四元数。如图 2所示，梯度仪坐标系由上标GRF表示，并与卫星主轴对齐，其中xGRF、yGRF和zGRF轴分别指向飞行方向、垂直于轨道平面方向和地心方向。3个星象仪的坐标系由对应的SSRFi(i=1, 2, 3)表示，在星象仪安置过程中，每个SSRFi的z轴已与视轴对齐，x轴和y轴均在焦平面内[23]。3个星象仪都指向稍微远离地心(zGRF负轴)的位置，STR1和STR2的视轴约在yGRFozGRF平面内，STR3的视轴略微向飞行方向倾斜。图 2中显示出两两星象仪STR1/2、STR1/3和STR2/3视轴之间的夹角约为41°、40°和55°，该3个角度为星象仪中姿态数据由SSRFi转换到GRF中的重要参数。
图 2 星象仪参考系相对于梯度仪参考系的方位
Figure 2. Position of Star Sensor Reference System Relative to Gradiometer Reference System
-
基于Romans [24]提出的卫星姿态组合方法，结合GOCE卫星中星象仪数据特点，将其应用于多星象仪姿态数据联合。该方法将噪声在星象仪坐标系各个轴分布特点建模为联合计算的加权矩阵，并通过最小二乘法结合两个星象仪的姿态四元数得到最佳四元数q* [25]。
将星象仪测量的四元数进行建模，计算式为：
$$ \boldsymbol{q}^{\mathrm{m}\ \mathrm{IRF} \rightarrow \operatorname{SSRF} x}=\boldsymbol{q}^{\text {true } \mathrm{IRF} \rightarrow \mathrm{SSRF} x} \boldsymbol{q}^{\text {noise } \mathrm{STR} x} $$ (8) 式中，qm IRF→SSRFx为测得的IRF到SSRFx的姿态四元数，即单个星象仪STRx中的测量数据，x∈{1, 2, 3}为不同星象仪的标识；qtrue IRF→SSRFx为IRF到SSRFx的真实姿态四元数；qnoise STRx为STRx测得姿态四元数与真实姿态四元数间的噪声。
qnoise STRx具体形式为：
$$ \boldsymbol{q}^{\text {noise STR}x}=\left[\begin{array}{c} 1 \\ 0.5 e_{1}^{\mathrm{STR}x} \\ 0.5 e_{2}^{\mathrm{STR} x} \\ 0.5 e_{3}^{\mathrm{STR} x} \end{array}\right] $$ (9) 式中，e1STRx、e2STRx、e3STRx为STRx中引入的小角度偏转来代替噪声；0.5为引入的比例系数，则以下条件成立：
$$ \boldsymbol{q}^{\mathrm{GRF} \rightarrow \mathrm{SSRF} x} \boldsymbol{q}^{\text {true } \mathrm{SSRF} x \rightarrow \mathrm{IRF}} \boldsymbol{q}^{\text {true } \mathrm{IRF} \rightarrow \mathrm{SSRF} y} \boldsymbol{q}^{\mathrm{SSRF} y \rightarrow \mathrm{GRF}}=\left[\begin{array}{l} 1 \\ 0 \\ 0 \\ 0 \end{array}\right] $$ (10) 式中，qGRF→SSRFx为GRF到SSRFx的姿态四元数；qtrue SSRFx→IRF为SSRFx到IRF的真实姿态四元数；qtrue IRF→SSRFy为IRF到SSRFy的真实姿态四元数，y∈{1, 2, 3}为不同星象仪的标识；qSSRFy→GRF为SSRFy到GRF的姿态四元数。
对于两个星象仪STRx、STRy的姿态四元数测量值，存在较小的相对误差：
$$ \begin{array}{c} \boldsymbol{q}^{\mathrm{GRF} \rightarrow \mathrm{SSRF} x} \boldsymbol{q}^{\mathrm{m}\ \mathrm{SSRF} x \rightarrow \mathrm{IRF}} \boldsymbol{q}^{\mathrm{m}\ \mathrm{IRF} \rightarrow \mathrm{SSRF} y} \boldsymbol{q}^{\mathrm{SSRF} y \rightarrow \mathrm{GRF}}=\\ \left[\begin{array}{c} 1 \\ 0.5 d_{1}^{\mathrm{STR} x y} \\ 0.5 d_{2}^{\mathrm{STR}xy} \\ 0.5 d_{3}^{\mathrm{STR} x y} \end{array}\right] \end{array} $$ (11) 式中，qm SSRFx→IRF为测得IRF到SSRFx姿态四元数的逆；qm IRF→SSRFy为测得的IRF到SSRFy的姿态四元数，即单个星象仪STRy中的测量数据。
STRx与STRy测量数据之间的相对误差dSTRxy与星象仪的噪声有关：
$$ \begin{gathered} \boldsymbol{d}^{\mathrm{STR} x y}=\left[\begin{array}{l} d_{1}^{\mathrm{STR} x y} \\ d_{2}^{\mathrm{STR} x y} \\ d_{3}^{ \mathrm{STR} x y} \end{array}\right]=\boldsymbol{R}^{\mathrm{SSRF} y \rightarrow \mathrm{GRF}}\left[\begin{array}{l} e_{1}^{\mathrm{STR} y} \\ e_{2}^{\mathrm{STRy}} \\ e_{3}^{\mathrm{STR} y} \end{array}\right]- \\ \boldsymbol{R}^{\mathrm{SSRF} x \rightarrow \mathrm{GRF}}\left[\begin{array}{c} e_{1}^{\mathrm{STR} x} \\ e_{2}^{\mathrm{STR} x} \\ e_{3}^{\mathrm{STR} x} \end{array}\right] \end{gathered} $$ (12) 式中，RSSRFx→GRF、RSSRFy→GRF分别为SSRFx、SS-RFy到GRF的旋转矩阵；e1STRy、e2STRy、e3STRy为STRy中引入的小角度偏转来代替噪声。
通过最小化加权平方和来获得最佳四元数，计算式为：
$$ S = \sum {{{\left( {{\mathit{\boldsymbol{e}}^{{\rm{STR}}x}}} \right)}^{\rm{T}}}} {\mathit{\boldsymbol{p}}^{{\rm{STR}}x}}{\mathit{\boldsymbol{e}}^{{\rm{STR}}x}} $$ (13) 式中，S为最小化加权平方和的约束函数；eSTRx为星象仪STRx的噪声：
$$ \boldsymbol{e}^{\mathrm{STR} x}=\left[\begin{array}{lll} e_{1}^{\mathrm{STR} x} & e_{2}^{\mathrm{STR} x} & e_{3}^{\mathrm{STR} x} \end{array}\right]^{\mathrm{T}} $$ (14) pSTRx为STRx各轴噪声的权重矩阵：
$$ {\mathit{\boldsymbol{p}}^{{\rm{STR}}x}} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ 0&0&{0.01} \end{array}} \right] $$ (15) 式中，0.01反映了星象仪视轴上的姿态精度比焦平面上其他轴的姿态精度低10倍。式(13)改写为：
$$ S = \sum {{{\left( {{{\mathit{\boldsymbol{\tilde e}}}^{{\rm{STR}}x}}} \right)}^{\rm{T}}}} {\mathit{\boldsymbol{\tilde p}}^{{\rm{STR}}x}}{\mathit{\boldsymbol{\tilde e}}^{{\rm{STR}}x}} $$ (16) 式中，$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$为转化为梯度仪坐标系下的星象仪噪声；$\tilde{\boldsymbol{p}}^{\mathrm{STR} x}$为通过方差传播计算得到的权重矩阵。计算式分别为：
$$ {{\mathit{\boldsymbol{\tilde e}}}^{{\rm{STR}}x}} = {\mathit{\boldsymbol{R}}^{{\rm{SSRF}}x \to {\rm{GRF}}}}{\mathit{\boldsymbol{e}}^{{\rm{STR}}x}} $$ (17) $$ {{\mathit{\boldsymbol{\tilde p}}}^{{\rm{STR}}x}} = {\mathit{\boldsymbol{R}}^{{\rm{SSRF}}x \to {\rm{GRF}}}}{\mathit{\boldsymbol{p}}^{{\rm{STR}}x}}{\left( {{\mathit{\boldsymbol{R}}^{{\rm{SSRF}}x \to {\rm{GRF}}}}} \right)^{\rm{T}}} $$ (18) 重新排列式(12)可得：
$$ {{\mathit{\boldsymbol{\tilde e}}}^{{\rm{STR}}y}} = {\mathit{\boldsymbol{d}}^{{\rm{STR}}xy}} + {{\mathit{\boldsymbol{\tilde e}}}^{{\rm{STR}}x}} $$ (19) 因为$\tilde{\boldsymbol{e}}^{\mathrm{STR} y}$可用$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$表达，所以选择$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$ 作为最小二乘法的参数，根据参数$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$ 表示S:
$$ \begin{gathered} S=\sum\limits_{x}\left(\tilde{\boldsymbol{e}}^{\mathrm{STR} x}\right)^{\mathrm{T}} \tilde{\boldsymbol{p}}^{\mathrm{STR} x} \tilde{\boldsymbol{e}}^{\mathrm{STR} x}= \\ \left(\tilde{\boldsymbol{e}}^{\mathrm{STR} x}\right)^{\mathrm{T}} \tilde{\boldsymbol{p}}^{\mathrm{STR} x} \tilde{\boldsymbol{e}}^{\mathrm{STR} x}+\sum\limits_{x \neq y}\left(\tilde{\boldsymbol{e}}^{\mathrm{STR} y}\right)^{\mathrm{T}} \tilde{\boldsymbol{p}}^{\mathrm{STR}y} \tilde{\boldsymbol{e}}^{\mathrm{STR} y}= \\ \left(\tilde{\boldsymbol{e}}^{\mathrm{STR} x}\right)^{\mathrm{T}} \tilde{\boldsymbol{p}}^{\mathrm{STR} x} \tilde{\boldsymbol{e}}^{\mathrm{STR} x}+\sum\limits_{x \neq y}\left(\boldsymbol{d}^{\mathrm{STR} x y}+\right. \\ \left.\tilde{\boldsymbol{e}}^{\mathrm{STR} x}\right)^{\mathrm{T}} \tilde{\boldsymbol{p}}^{\mathrm{STR} y}\left(\boldsymbol{d}^{\mathrm{STR} x y}+\tilde{\boldsymbol{e}}^{\mathrm{STR} x}\right) \end{gathered} $$ (20) 通过设置S关于$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$ 阶导数为0，求解$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$ :
$$ \tilde{\boldsymbol{e}}^{\mathrm{STR} x}=-\left(\sum\limits_{y} \tilde{\boldsymbol{p}}^{\mathrm{STR} y}\right)^{-1}\left(\sum\limits_{y \neq x} \tilde{\boldsymbol{p}}^{\mathrm{STR}y} \boldsymbol{d}^{\mathrm{STR} x y}\right) $$ (21) 将解算得到的$\tilde{\boldsymbol{e}}^{\mathrm{STR} x}$ 代入式(10)中，可计算得到最佳四元数q*:
$$ \begin{aligned} \boldsymbol{q}^{*}=&\boldsymbol{q}^{\mathrm{opt}\ \mathrm{IRF} \rightarrow \mathrm{SSRF} x}=\boldsymbol{q}^{\mathrm{m}\ \mathrm{IRF} \rightarrow \mathrm{SSRF} x}\left(\boldsymbol{q}^{\text {noise } \mathrm{STR} x}\right)^{*}=\\ &\boldsymbol{q}^{\mathrm{m} \operatorname{IRF} \rightarrow \mathrm{SSRF} x}\left[\begin{array}{c} 1 \\ -0.5 \boldsymbol{R}^{\mathrm{GRF} \rightarrow \mathrm{SSRF} x} \tilde{\boldsymbol{e}}^{\mathrm{STR} x} \end{array}\right] \end{aligned} $$ (22) -
本文选取欧洲空间局发布的2013-10-16星象仪观测数据文件(STR_VC2/3_1b)，在测试计算中，首先对STR_VC2_1b和STR_VC3_1b中四元数数据进行连续处理，剔除含有无效数据(valid data flag, Val_Flag)和视场中含有异物(big and bright object, BBO)的四元数数据，然后分别对单星象仪及多星象仪四元数进行联合校准计算，最后通过四元数方程计算各个轴的角速度w[26]:
$$ \mathit{\boldsymbol{\dot q}} = \frac{1}{2}\mathit{\boldsymbol{qw}} $$ (23) 表 1列出了本文星象仪姿态数据联合所需的数据文件，读取星象仪数据文件时选择有效(Val_Flag=1)且在视场中无异物(BBO=0)的姿态数据[27]。
表 1 星象仪姿态四元数校准和组合所需的GOCE卫星数据
Table 1. GOCE Data Required for Attitude Quaternion Calibration and Combination of Star Tracker
数据名称 数据描述 STR_VC2_1b 星象仪标识符(STR_ID)，姿态四元数(Q1, Q2, Q3, Q4)，有效标志(Val_Flag)，星象仪视场中异物(BBO) STR_VC3_1b 星象仪标识符(STR_ID)，姿态四元数(Q1, Q2, Q3, Q4)，有效标志(Val_Flag)，剩余星象仪视场中异物(BBO) AUX_EGG_DB 从星象仪参考系旋转到与卫星主轴对齐的梯度仪参考系 -
由于选择的数据含有STR1和STR2有效姿态数据，在联合校准之前，先对单个星象仪进行计算，分析各个轴的角速度的噪声精度分布情况，以便于后续在坐标系转化过程中对比噪声的传播。在SSRF1下将STR1姿态数据由式(23)解算得到相应的角速度，图 3(a)为x、y、z轴上角速度分量的时间序列图，通过对比其随时间平滑起伏波动程度可以分析各轴角速度噪声的分布特征。其中wx、wy呈现出噪声低而wz呈现出噪声高的特征，且wz的噪声明显高于wx、wy。图 3(b)为x、y、z轴上角速度平方根的功率谱密度(power spectral density, PSD) [28]，在低于1×10-2 Hz的低频区域，3个轴的角速度wx、wy、wz精度相似；约在1×10-2 Hz处，wz噪声开始高于wx、wy的噪声，其精度降低；约在1×10-1 Hz处，wz相对于wx、wy噪声高出约一个量级，并且噪声随着频率的增加显示占据主导地位。分析认为，这与GOCE卫星搭载的星象仪三轴的指向直接相关。z轴的角速度分量垂直于星象仪的焦平面并指向视轴方向，x轴和y轴位于STR的焦平面中。星象仪在飞行过程中的视轴为低灵敏轴，其相对精度为100 μrad；其余两个位于焦平面中的轴为超灵敏轴，其相对精度为1 μrad。因此，对绕SSRF轴角速度的精度而言，wz相较于wx、wy精度更低，其噪声水平更高。
图 3 SSRF1下STR1三个轴向角速度分量的时间序列和平方根的功率谱密度
Figure 3. Time Series and Square Root PSD of Three Axial Angular Velocity Components of STR1 in SSRF1
-
在GOCE卫星重力梯度数据处理过程中，为了星象仪数据与梯度仪数据相关联，通过STR1对应的旋转矩阵将角速度由SSRF1转换为GRF，图 4为GRF下单个星象仪角速度的时间序列与平方根功率谱密度。
图 4 GRF下STR1三个轴向角速度分量的时间序列和平方根的功率谱密度
Figure 4. Time Series and Square Root PSD of Three Axial Angular Velocity Components of STR1 in GRF
由图 3(a)和图 4(a)对比可知，GRF下wy随时间起伏波动变大，呈现出噪声高、精度降低的特征，且与wx相比噪声更大，精度更差。将图 4(b)与图 3(b)作对比分析，wy与wz噪声呈现出相似性；约在1×10-2 Hz处，wy的噪声开始变大、精度降低；约在1×10-1 Hz处，wy、wz相对于wx噪声高出约一个量级；而wx则几乎无变化，其精度不受坐标系转化的影响。这是由于在将STR1的角速度从SSRF1转换为GRF时，低灵敏轴的噪声会传播到超灵敏轴中。由SSRF1到GRF的转换主要是围绕xSSRF1轴进行的，这导致SSRF1下精度较差的wz将噪声传播到GRF下的wy和wz。
-
在GRF下将STR1与STR2联合计算各个轴的角速度标准差，与单个星象仪STR1计算结果的对比见表 2。表 3则是联合前后计算Allan方差得到的系数。从表 2中可以看出，STR1与STR2联合后，三轴角速度的总精度提高了1/3, x轴精度提升几乎不变，主要是y轴与z轴的精度的提升，其中y轴标准差提升约4.6个量级。表 3的Allan方差系数对比分析与表 2有类似的结果，x轴系数几乎不变，而y轴与z轴的各个系数均有降低。这一结论也体现在图 5中，即STR1与STR2联合前后角速度分量wx、wy、wz在时间序列上的对比。
表 2 GRF下STR1和STR2联合前后角速度标准差/(rad·s-1)
Table 2. Standard Deviation of Angular Velocities of STR1 and STR1/2 in GRF/(rad·s-1)
星象仪 三轴总标准差/10-9 x轴标准差/10-5 y轴标准差/10-5 z轴标准差/10-5 STR1 8.755 35 6.144 83 7.094 67 7.802 54 STR1/2 5.914 04 6.126 57 1.554 77 6.755 63 表 3 STR1与STR1/2 Allan方差系数
Table 3. Allan Variance Coefficients of STR1 and STR1/2
星象仪 轴向 量化噪声/(°) 角度随机游走/(°·h-0.5) 零偏不稳定性/(°·h-1) 角速度游走/(°·h-1·h-0.5) 角速度斜坡/(°·h-1·h-1) STR1 x 0.164 141 0.004 925 0.159 934 0.061 761 0.010 308 y 0.089 338 0.003 716 0.042 224 0.038 977 0.010 336 z 0.185 804 0.004 940 0.173 201 0.071 214 0.007 549 STR1/2 x 0.164 067 0.004 932 0.159 969 0.061 828 0.010 281 y 0.028 512 0.000 368 0.017 341 0.012 722 0.003 061 z 0.150 297 0.004 377 0.145 182 0.065 010 0.004 319 
图 5 STR1与STR2联合前后3个轴角速度时间序列
Figure 5. Time Series of Three Axial Angular Velocity Components Between STR1 and STR1/2
由图 5可以看出，双星象仪STR1与STR2联合后计算得到的角速度效果明显更好，与单星象仪解算的角速度相比具有噪声更低、精度更好的优势；联合后的wy、wz噪声已得到有效降低，而wx的改善很小。分析认为，这是由于单星和双星联合都是围绕SSRF的x轴的旋转，而SSRF的x轴为具有更高精度的视轴，STR1与STR2联合后计算得到的角速度精度均有明显的提升。
图 6为STR1与STR2联合后3个轴向角速度分量的时间序列和平方根PSD，比较图 4(b)与图 6(b)可以看出，与单个STR1相比，STR1与STR2联合后显示出y轴、z轴的角速度分量wy(GRF)、wz(GRF)的总体噪声水平明显下降；在1×10-2 Hz处约下降了1个量级，其精度达到10-5 rad/s量级；且wz与wx的噪声水平呈现出相似性；在3×10-3~3×10-2 Hz处，wy比wx的噪声水平更低、精度更高；在3×10-2 Hz以上，wy则与wx、wz的噪声水平呈现出相似性。对图 6(a)与图 6(b)分析可得，两星联合产生的角速度不会受到单星的视轴组件测量精度较低的影响，其有效地抑制了由于SSRF-GRF坐标系变换导致精度较低的角速度分量传播到其他分量。
图 6 GRF下STR1/2联合后三个轴向角速度分量的时间序列和平方根的功率谱密度
Figure 6. Time Series and Square Root PSD of Three Axial Angular Velocity Components of STR1/2 in GRF
-
基于重力梯度测量原理进一步对多星象仪联合后的姿态数据进行计算得到梯度，以分析该联合方法对重力梯度的影响[29]。为了分析STR1与STR2联合后对重力梯度的影响，暂不考虑梯度仪角速度的影响，分别对单星象仪和多星象仪联合后的姿态数据进行梯度解算，通过比较重力梯度迹的平方根PSD来评价该联合方法计算得到的梯度精度[30]，如图 7所示。
图 7 STR1与STR1/2重力梯度迹平方根的功率谱密度
Figure 7. Square Root PSD of Gravity Gradient Trace of STR1 and STR1/2
由图 7可以看出，在中低频下，STR1单星计算的梯度在1 cpr(约为1.85×10-4 Hz)的谐波具有最大的噪声PSD, STR1与STR2联合后相对于STR1单星计算的梯度整个频域都有显著改善，特别是在3~200 m Hz内有着较为明显的改善，大约低了0.5个量级，这表明多星联合法的性能明显优于单星方法，计算出的重力梯度具有更低的噪声和更高的精度，在后续解算重力场时也可以进一步分析，提高重力场的精度[31]。
-
本文提出了一种星象仪姿态数据的联合算法，系统分析了星象仪噪声特性，通过构建反映STR测量的各个轴噪声分布的加权矩阵，进行多星象仪姿态数据联合，并计算出最佳姿态四元数。结果表明，多星象仪联合产生的角速度不会受到单星象仪中视轴测量精度较低的影响，能够有效抑制由于SSRF-GRF坐标系变换而导致精度较低的角速度分量对其他分量影响。
通过单星象仪与多星象仪联合法计算的角速度分别计算重力梯度，比较了两组重力梯度迹平方根的PSD，分析认为两星联合后的数据可改善3~200 m Hz内重力梯度精度。GOCE卫星实际在轨运行搭载了3个星象仪，保证了随时都可提供至少2个星象仪的有效数据用于组合解算姿态数据。基于加权矩阵的组合方法的优点在于，将2个或者所有的3个星象仪进行组合可以有效地提高数据精度，而星象仪姿态数据精度的改善将直接反映在重力场中。因此，本文提出的多星象仪联合解算方法提供了一个较好的姿态解决方案，可以应用于今后中国重力卫星任务的预研和实施。同时，下一步重力仪梯度数据预处理工作还将综合考虑温度对星象仪姿态数据的影响，磁力计数据对星象仪姿态数据的校正，以及星象仪星图识别精化算法，这些预处理工作都将有效提高重力梯度科学数据产品的精度。
Reconstruction Method of Satellite Gravity Gradient Measurement Angular Velocity by Combining Star Tracker Quaternion
-
摘要: 在卫星重力梯度测量的level 1b (L1b)数据处理过程中，为了计算得到高精度的重力梯度，需精确测定卫星姿态角速度。在星象仪姿态数据坐标系转换中，单星象仪中的低精度角速度分量会将噪声传递到其他角速度分量，进而影响整个重力梯度的计算精度。基于星象仪噪声特性，构建星象仪各轴的噪声分布加权矩阵，对两个或多个星象仪的姿态四元数联合求解最佳姿态四元数，为后续角速度恢复提供精确姿态控制。计算结果表明，多星象仪联合解算得到的卫星姿态角速度在整个功率谱密度上实现了显著改善。与单星象仪计算方法相比，多星象仪联合得到的最佳姿态角速度wy、wz在10~100 mHz内精度约提升1个量级，噪声水平约达到10-6 rad/s量级，能够有效抑制低精度角速度分量在坐标系转换中导致的噪声传播。此外，基于多星象仪联合方法计算得到的重力梯度迹值的功率谱密度，相比单星象仪方法也有较为显著的改善。Abstract:
Objectives In order to calculate the high precision gravity gradient in level 1b (L1b) data processing of gravity field and steady-state ocean circulation explorer(GOCE)satellite, it is necessary to accurately measure the satellite attitude angular velocity. In the coordinate system transformation of the attitude data of astrology, the low precision angular velocity components of a single astrology instrument will transmit the noise to other angular velocity components, which will affect the calculation accuracy of the whole gravity gradient. Methods In this paper, based on the noise characteristics of the astrological instrument, the weighted matrix of the noise distribution on each axis of the astrological instrument is constructed, and the optimal attitude quaternion is solved jointly for the attitude quaternion of two or more astrological instruments, so as to provide accurate attitude control for subsequent angular velocity recovery. Results Compared with the single astro-meter calculation method, the accuracy of the optimal attitude angular velocity components wy and wz obtained by the combination of multiple astro-meters can reach about 1 order of magnitude higher within 10-100 mHz, and the noise level is up to 10-6 rad/s, which can effectively suppress the noise propagation caused by the low precision angular velocity components in coordinate system transformation. Conclusions The results show that the satellite attitude angular velocity calculated by the combination of multiple astrological instruments can be significantly improved in the whole power spectral density, the power spectral density of gravity gradient trace values calculated by the combined method of multiple astrologers is also significantly improved compared with the single astrologer method. -
图 2 星象仪参考系相对于梯度仪参考系的方位
Figure 2. Position of Star Sensor Reference System Relative to Gradiometer Reference System
图 3 SSRF1下STR1三个轴向角速度分量的时间序列和平方根的功率谱密度
Figure 3. Time Series and Square Root PSD of Three Axial Angular Velocity Components of STR1 in SSRF1
图 4 GRF下STR1三个轴向角速度分量的时间序列和平方根的功率谱密度
Figure 4. Time Series and Square Root PSD of Three Axial Angular Velocity Components of STR1 in GRF
图 5 STR1与STR2联合前后3个轴角速度时间序列
Figure 5. Time Series of Three Axial Angular Velocity Components Between STR1 and STR1/2
图 6 GRF下STR1/2联合后三个轴向角速度分量的时间序列和平方根的功率谱密度
Figure 6. Time Series and Square Root PSD of Three Axial Angular Velocity Components of STR1/2 in GRF
图 7 STR1与STR1/2重力梯度迹平方根的功率谱密度
Figure 7. Square Root PSD of Gravity Gradient Trace of STR1 and STR1/2
表 1 星象仪姿态四元数校准和组合所需的GOCE卫星数据
Table 1. GOCE Data Required for Attitude Quaternion Calibration and Combination of Star Tracker
数据名称 数据描述 STR_VC2_1b 星象仪标识符(STR_ID)，姿态四元数(Q1, Q2, Q3, Q4)，有效标志(Val_Flag)，星象仪视场中异物(BBO) STR_VC3_1b 星象仪标识符(STR_ID)，姿态四元数(Q1, Q2, Q3, Q4)，有效标志(Val_Flag)，剩余星象仪视场中异物(BBO) AUX_EGG_DB 从星象仪参考系旋转到与卫星主轴对齐的梯度仪参考系 
下载: 导出CSV
表 2 GRF下STR1和STR2联合前后角速度标准差/(rad·s-1)
Table 2. Standard Deviation of Angular Velocities of STR1 and STR1/2 in GRF/(rad·s-1)
星象仪 三轴总标准差/10-9 x轴标准差/10-5 y轴标准差/10-5 z轴标准差/10-5 STR1 8.755 35 6.144 83 7.094 67 7.802 54 STR1/2 5.914 04 6.126 57 1.554 77 6.755 63
下载: 导出CSV
表 3 STR1与STR1/2 Allan方差系数
Table 3. Allan Variance Coefficients of STR1 and STR1/2
星象仪 轴向 量化噪声/(°) 角度随机游走/(°·h-0.5) 零偏不稳定性/(°·h-1) 角速度游走/(°·h-1·h-0.5) 角速度斜坡/(°·h-1·h-1) STR1 x 0.164 141 0.004 925 0.159 934 0.061 761 0.010 308 y 0.089 338 0.003 716 0.042 224 0.038 977 0.010 336 z 0.185 804 0.004 940 0.173 201 0.071 214 0.007 549 STR1/2 x 0.164 067 0.004 932 0.159 969 0.061 828 0.010 281 y 0.028 512 0.000 368 0.017 341 0.012 722 0.003 061 z 0.150 297 0.004 377 0.145 182 0.065 010 0.004 319
下载: 导出CSV
-
[1] 程鹏飞, 文汉江, 刘焕玲, 等. 卫星大地测量学的研究现状及发展趋势[J]. 武汉大学学报·信息科学版, 2019, 44(1): 48-54 doi: 10.13203/j.whugis20180356 Cheng Pengfei, Wen Hanjiang, Liu Huanling, et al. Research Situation and Future Development of Satellite Geodesy [J]. Geomatics and Information Science of Wuhan University, 2019, 44(1): 48-54 doi: 10.13203/j.whugis20180356 [2] Baur O, Bock H, Höck E, et al. Comparison of GOCE-GPS Gravity Fields Derived by Different Approaches[J]. Journal of Geodesy, 2014, 88(10): 959-973 doi: 10.1007/s00190-014-0736-6 [3] 吴云龙. GOCE卫星重力梯度测量数据的预处理研究[D]. 武汉: 武汉大学, 2010 Wu Yunlong. Study on Pre-Processing of GOCE Satellite Gravity Gradiometry Data[D]. Wuhan: Wuhan University, 2010 [4] 周泽兵, 白彦峥, 祝竺, 等. 卫星重力测量中加速度计在轨参数校准方法研究[J]. 中国空间科学技术, 2009, 29(6): 74-80 doi: 10.3321/j.issn:1000-758X.2009.06.012 Zhou Zebing, Bai Yanzheng, Zhu Zhu, et al. In-Orbit Calibration Methods of Accelerometer Parameters on Satellite-Borne Gravimetry [J]. Chinese Space Science and Technology, 2009, 29(6): 74-80 doi: 10.3321/j.issn:1000-758X.2009.06.012 [5] Rummel R, Yi W Y, Stummer C. GOCE Gravitational Gradiometry[J]. Journal of Geodesy, 2011, 85(11): 777-790 doi: 10.1007/s00190-011-0500-0 [6] Frommknecht B, Lamarre D, Meloni M, et al. GOCE Level 1b Data Processing[J]. Journal of Geodesy, 2011, 85(11): 759-775 doi: 10.1007/s00190-011-0497-4 [7] Siemes C, Rexer M, Schlicht A, et al. GOCE Gradiometer Data Calibration[J]. Journal of Geodesy, 2019, 93(9): 1 603-1 630 doi: 10.1007/s00190-019-01271-9 [8] Siemes C, Haagmans R, Kern M, et al. Monitoring GOCE Gradiometer Calibration Parameters Using Accelerometer and Star Sensor Data: Methodology and First Results[J]. Journal of Geodesy, 2012, 86(8): 629-645 doi: 10.1007/s00190-012-0545-8 [9] Rispens S, Bouman J. Calibrating the GOCE Accelerations with Star Sensor Data and a Global Gravity Field Model[J]. Journal of Geodesy, 2009, 83(8): 737-749 doi: 10.1007/s00190-008-0290-1 [10] Stummer C, Fecher T, Pail R. Alternative Method for Angular Rate Determination Within the GOCE Gradiometer Processing[J]. Journal of Geodesy, 2011, 85(9): 585-596 doi: 10.1007/s00190-011-0461-3 [11] Siemes C, Rexer M, Haagmans R. GOCE Star Tracker Attitude Quaternion Calibration and Combination[J]. Advances in Space Research, 2019, 63(3): 1 133-1 146 http://www.sciencedirect.com/science/article/pii/S0273117718307993 [12] Bandikova T, Flury J. Improvement of the GRACE Star Camera Data Based on the Revision of the Com bination Method[J]. Advances in Space Research, 2014, 54(9): 1 818-1 827 doi: 10.1016/j.asr.2014.07.004 [13] Wu S C, Kruizinga G, Bertiger W. Algorithm Theoretical Basis Document for GRACE Level-1B Data Processing v1. 2[R]. Los Angeles, USA: Jet Propulsion Laboratory, 2006 [14] Cesare C, Catastini G. Gradiometer On-Orbit Calibration Procedure Analysis[R]. Paris, France: European Space Agency, 2008 [15] Cesare C, Sechi G, Catastini G. Gradiometer Ground Processing Algorithms Documentation[R]. Paris, France: European Space Agency, 2008 [16] 杨宇飞, 杨元喜, 徐君毅, 等. 低轨卫星对导航卫星星座轨道测定的增强作用[J]. 武汉大学学报·信息科学版, 2020, 45(1): 46-52 doi: 10.13203/j.whugis20180215 Yang Yufei, Yang Yuanxi, Xu Junyi, et al. Navigation Satellites Orbit Determination with the Enhancement of Low Earth Orbit Satellites[J]. Geomatics and Information Science of Wuhan University, 2020, 45(1): 46-52 doi: 10.13203/j.whugis20180215 [17] 吴汤婷, 徐新禹, 魏辉, 等. 一种计算GOCE卫星加速度的实用算法[J]. 测绘地理信息, 2017, 42(2): 8-11 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG201702003.htm Wu Tangting, Xu Xinyu, Wei Hui, et al. An Applied Algorithm to Calculate the Acceleration of GOCE Satellite[J]. Journal of Geomatics, 2017, 42(2): 8-11 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG201702003.htm [18] Gruber T, Rummel R, Abrikosov O, et al. GOCE Level 2 Product Data Handbook[R]. Novotel, Luxembourg: The GOCE High Level Processing Facility, 2010 [19] Rummel R. Satellite Gradiometry[M]// Mathematical and Numerical Techniques in Physical Geodesy. Berlin, Heidelberg: Springer-Verlag, 1986: 317-363 [20] 邹贤才, 李建成, 衷路萍, 等. 动力法校准GRACE星载加速度计[J]. 武汉大学学报·信息科学版, 2015, 40(3): 357-360 http://ch.whu.edu.cn/article/id/3210 Zou Xiancai, Li Jiancheng, Zhong Luping, et al. Calibration of the Accelerometers Onboard GRACE with the Dynamic Method[J]. Geomatics and Information Science of Wuhan University, 2015, 40(3): 357-360 http://ch.whu.edu.cn/article/id/3210 [21] 陈俊勇. 地球重力场、重力、重力梯度在三维直角坐标系中的表达式[J]. 武汉大学学报·信息科学版, 2004, 29(5): 377-379 doi: 10.13203/j.whugis2004.05.001 Chen Junyong. Expressions of Earth Gravity Field, Gravity and Gravity Gradient in 3D Cartesian Coordinates System[J]. Geomatics and Information Science of Wuhan University, 2004, 29(5): 377-379 doi: 10.13203/j.whugis2004.05.001 [22] Rispens S, Bouman J. Calibrating the GOCE Accelerations with Star Sensor Data and a Global Gravity Field Model[J]. Journal of Geodesy, 2009, 83(8): 737-749 doi: 10.1007/s00190-008-0290-1 [23] Frommknecht B. Integrated Sensor Analysis of the GRACE Mission[R]. München, Germany: Verlag der Bayerischen Akademieder Wissenschaften, 2008 [24] Romans L. Optimal Combination of Quaternions From Multiple Star Cameras[R]. Los Angeles, USA: Jet Propulsion Laboratory, 2003 [25] Siemes C. GOCE Level 1b Gravity Gradient Processing Algorithms[R]. Paris, France: European Space Agency, 2018 [26] 章仁为. 卫星轨道姿态动力学与控制[M]. 北京: 北京航空航天大学出版社, 1998 Zhang Renwei. Dynamics and Control of Satellite Orbital Attitude[M]. Beijing: Beijing University of Aeronautics and Astronautics Press, 1998 [27] SERCO/DATAMAT Consortium. GOCE L1b Products User Handbook[R]. Paris, France: European Space Agency, 2006 [28] Welch P. The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging over Short, Modified Periodograms[J]. IEEE Transactions on Audio and Eletroacoustics, 1967, 15(2): 70-73 doi: 10.1109/TAU.1967.1161901 [29] 周晓青. 探测地球重力场的卫星重力梯度指标研究与分析[D]. 武汉: 武汉大学, 2016 Zhou Xiaoqing. Research and Analysis of Satellite Gravity Gradiometry Indicator for Exploring the Earth's Gravity Field[D]. Wuhan: Wuhan University, 2016 [30] Rummel R, Yi W Y, Stummer C. GOCE Gravitational Gradiometry[J]. Journal of Geodesy, 2011, 85(11): 777-790 doi: 10.1007/s00190-011-0500-0 [31] Wan X Y, Yu J H. Derivation of the Radial Gradient of the Gravity Based on Non-Full Tensor Satellite Gravity Gradients[J]. Journal of Geodynamics, 2013, 66: 59-64 doi: 10.1016/j.jog.2013.02.005 -
点击查看大图
计量
- 文章访问数: 596
- HTML全文浏览量: 189
- PDF下载量: 63
- 被引次数: 0
