吴云龙， 郭泽华， 肖云， 马林

1 地震大地测量重点实验室, 中国地震局地震研究所, 武汉 430071 2 防灾科技学院, 河北三河 065201 3 引力与固体潮国家野外科学观测研究站， 武汉 430071 4 西安测绘研究所， 西安 710054 5 航天东方红卫星有限公司， 北京 100094

0 引言

全球重力场和海洋环流探测器(Gravity field and steady-state Ocean Circulation Explorer, GOCE)由欧洲空间局(European Space Agency，ESA)于2009年3月17日成功发射，在轨运行4年后于2013年11月11日结束卫星任务.GOCE卫星主要任务目标是在100 km的空间分辨率下确定精度为1～2 cm的全球大地水准面和精度为1 mGal的全球重力异常(ESA, 1999).为实现这一极具难度的科学目标，GOCE卫星搭载了高精度静电重力梯度仪(Electrostatic Gravity Gradiometry，EGG)、高精度全球导航卫星系统(Global Navigation Satellite System, GNSS)星载接收机和星敏感器(Star Sensor，STR)等多个关键载荷(Cesare, 2008; Rummel et al., 2011).

高精度GOCE重力梯度观测数据在地球科学的多个领域应用广泛.国际地球重力场模型中心(International Centre for Global Earth Models, ICGEM)发布的全球最优重力场模型产品系列中，相当一部分产品来源于卫星重力梯度观测数据(Drewes et al.，2016).国内外研究同行将GOCE重力梯度观测数据及其模型产品在大地测量学、固体地球物理学、海洋学以及冰川学等领域开展了持续深入的应用研究(Sun and Okubo, 2004; Tapley et al., 2004; 孙文科, 2008; 钟敏等, 2009; Knudsen et al., 2011; Gruber et al., 2012; Becker et al., 2014; Fuchs et al., 2014; Hirt, 2014; Bouman et al., 2014).但是，所有这些地学领域应用都取决于能否获取高精度重力梯度测量观测数据.此外，为更好发挥重力梯度仪的观测性能，GOCE卫星设计运行在极低的轨道高度(～259 km)，以最大效应感知重力梯度信息，从而对卫星轨道高度维持和高精度姿态控制都提出了更高要求(Floberghagen et al., 2011; Pail et al., 2011).GOCE卫星官方数据处理部门将构建高精度重力梯度数据作为L1级数据处理流程的核心环节，主要为将静电重力梯度仪得到的加速度数据、星敏感器得到的姿态数据，联合构建高精度重力梯度分量观测数据(ESA, 2009).因此，高精度的重力梯度数据L1级预处理可为恢复高精度静态地球重力场模型提供数据质量保障，是卫星重力梯度测量数据处理及应用中的重要环节，也是实现GOCE卫星预期科学目标的关键任务之一(吴云龙，2010).

国际上，欧洲空间局GOCE卫星任务的高级数据处理部门(High Level Processing Facility，HPF)负责L1级和L2级官方数据产品的处理和发布(Frommknecht et al., 2011).其分布在欧洲多个研究机构的10个研究小组围绕整个数据处理流程进行了系统深入的研究，主要包括角速度重建、梯度观测数据构建、科学数据预处理、卫星轨道精密定轨、重力场模型恢复(快速产品和科学产品)以及数据极空白改善等(Siemes et al., 2012; Wan and Yu, 2013; Baur et al., 2014; Visser et al., 2014, 2016; Visser, 2017; Siemes, 2018a, b; Siemes et al., 2019a, b).近年来，我国围绕发展民用重力梯度测量卫星也开展了相关预研究工作，仿真设计了新型重力测量卫星模式，模拟分析了不同重力测量卫星系统配置条件、关键载荷技术指标以及噪声水平下的重力场模型反演能力(冉将军等，2015; 徐新禹等，2018)；基于地面实测重力数据给出了重力梯度测量卫星极空白的弥补改善策略(Lu et al.，2020)；系统研究了重力梯度测量卫星数据预处理方法，包括重力梯度观测数据预处理中的潮汐/非潮汐效应改正、粗差探测方法、外部校准方法等(吴云龙, 2010)；分析了星敏感器噪声特性，提出了联合两个或多个星敏感器姿态数据求解最佳姿态四元数方法(郭泽华等, 2021)；提出了一种基于地面重力的卫星在轨检校方法，从地面先验重力数据的时空精度、重力梯度仪观测噪声等预处理中的外部检校环节开展了分析研究(瞿庆亮等，2021).总体而言，高精度重力梯度观测数据L1级构建的系统方法与技术，受限于国外原始数据和具体技术文献公开，国内少有学者在此研究领域开展较为系统的研究工作.随着我国推进自主重力卫星任务的发展需要，突破国外机构对重力卫星L0-L1级数据预处理的技术封锁，实现自主卫星重力梯度观测数据L1级构建能力，具有迫切技术需求(许厚泽, 2001; 宁津生, 2002).

本文面向我国发展的梯度卫星的任务需要，系统研究并初步实现了自主的卫星重力梯度观测数据L1级构建方法，主要包括加速度计电压数据转换、多星敏感器联合姿态数据的角速度重建、卫星重力梯度分量构建等数据处理技术环节.研究工作可为下一步我国推进实施民用重力梯度测量卫星任务提供自主的原始数据处理技术支撑与储备.

1 卫星重力梯度观测数据构建原理

从重力梯度测量卫星原始观测数据构建重力梯度分量，主要以重力梯度仪和星敏感器两种关键载荷得到的观测数据为基础(Rispens and Bouman, 2009).在卫星重力梯度测量中，观测量同时包括了离心加速度和角加速度.通过观测量对称性和非对称性的表达式，将角加速度分离.通过卫星搭载的星敏感器提供重力梯度仪的初始加速度矢量，对其积分可求得重力梯度仪的角速度分量，从而得到离心加速度，最终可将重力梯度从观测值中分离出来(Rummel et al., 2011).

为构建重力梯度分量，需要将其他信号由加速度计观测值中分离.根据卫星重力梯度测量原理(Stummer et al., 2012)，重力梯度仪中加速度计观测的加速度值ai为：

(1)

图1 梯度仪坐标系内梯度仪的6个加速度计位置Fig.1 Position of six accelerometers of gradiometer in gradiometer reference frame

GOCE卫星所搭载的6个加速度计都是通过刚性连接固定在梯度仪内，因此它们在卫星质心处的非保守力加速度相同，加速度计对距离矢量ri+rj≈0，ij∈{14，25，36}，可通过建立共模加速度确定.

定义共模加速度acij和差分加速度adij：

(2)

非保守力加速度d在作差分时相互抵消：

(3)

Ad=[ad,14ad,25ad,36],

(4)

(5)

利用V与Ω2的对称性可得：

AdL-1+(AdL-1)T=-2 (V-Ω2),

(6)

(7)

由式(7)可提取角加速度:

(8)

(9)

(10)

2 L1级观测数据构建方法

重力梯度测量卫星观测数据L1级构建主要包括电压数据转换、姿态数据重建、角速度重建和重力梯度分量构建等四个步骤，具体计算方法如下：

图2 单个加速度计电极结构Fig.2 Electrode structure of single accelerometer

2.1 加速度计电压数据转换

重力梯度仪中单个加速度计参考系(Accelerometer Reference Frame, ARF)的电极结构如图2所示，在加速度计的较小壁上各有两个电极(Y1,Y2,Z1,Z2)，在加速度计的较大壁上各有四个电极(X1,X2,X3,X4)，其中两个相对的电极形成一对.检测质量为扁平立方体(尺寸为4 cm×4 cm×1 cm)，在1 g环境下进行地面试验标定，利用较大壁上的四对电极完成对检测质量的悬浮.

Frommknecht等(2011)对控制电压-加速度转换的处理步骤做出了详细阐述.其关键环节是对控制电压观测数据施加静电增益(即增益因子)，再转化为加速度观测数据.需要注意的是，每个电极对增益因子均略有不同.将加速度计中8个极板的电压观测值通过与适当的静电增益重组矩阵相乘，将控制电压Vci转换为线性加速度与角加速度.

(11)

(12)

其中，VP为极化电压；ε0为真空介电常数；m为检测质块的质量；e为检测质块与极板间的距离.

2.2 多星敏感器联合重建姿态数据

根据郭泽华等(2021)提出的多星敏感器联合算法，构建噪声分布加权矩阵，可获得卫星最佳惯性姿态四元数.联合解算得到的角速度不会受到单个星敏感器的视轴测量精度较低的影响，能够有效抑制由于坐标系变换(SSRF-GRF)导致的精度较低角速度误差传播到其他分量(Siemes, 2018b).将星敏感器测得的四元数建模，计算得到最佳四元数q*为：

(13)

2.3 角速度重建

卫星角速度与姿态四元数的精度直接影响地球重力场模型反演精度.星敏感器、重力梯度仪观测数据确定角速度时，会出现两者角速度噪声分别在高频、低频增大的影响.为了精确地确定卫星在空间中的惯性角速度，其重建过程中应充分考虑到EGG与STR的误差特性.基于Stummer等(2011)提出的维纳滤波重建角速度方法，可将EGG与STR的角速度根据其精度在频域内联合.表1、表2分别为EGG和STR的噪声模型，STR和EGG三个角速度分量具有不同的噪声功率频谱密度，如图3所示.

表1 EGG角速度噪声模型Table 1 EGG angular velocity noise model

表2 STR角速度噪声模型Table 2 STR angular velocity noise model

图3 STR和EGG角速度噪声的PSD1/2Fig.3 PSD1/2 of STR and EGG angular velocity noise

由于重力梯度仪和星敏感器角速度的精度与频率有关，而功率谱密度P(f)可表示在频率f处的精度，可根据公式(14)、(15)计算维纳滤波的权重(Stummer et al., 2012)：

(14)

(15)

所有频率的权重之和等于1，由此可得到：

(16)

如图4所示，实线、虚线分别表示STR、EGG的维纳滤波权重.由于维纳滤波是在时域中联合角速度，滤波系数可通过离散傅里叶逆变换(Inverse Discrete Fourier Transform, IDFT)从权重中获得：

(17)

(18)

图4 EGG和STR维纳滤波权重Fig.4 EGG and STR Wiener filter weight

2.4 重力梯度数据构建

基于上述公式得到的加速度计观测数据和角速度观测数据，综合式(6)和(7)：

(19)

(20)

卫星星体在3个方向上所受到的非保守力可通过加速度计测量的共模加速度表示，其共模加速度存在以下关系：

E{ac,ij-ac,kl}=0;ij=14,25,36;ij≠kl

(21)

即共模加速度间存在以下线性组合关系：

ac,14,x-ac,25,x=0,

ac,14,x-ac,36,x=0,

ac,25,y-ac,36,y=0,

ac,25,y-ac,14,y=0,

ac,36,z-ac,25,z=0,

ac,36,z-ac,14,z=0.

(22)

通过式(3)得到差分加速度后，只有角速度保留在差分加速度中.为了从式(7)中得到重力梯度，应从差模加速度中分离出重力梯度.即V的主对角线重力梯度分量计算式为：

(23)

(24)

(25)

非对角线重力梯度分量计算式为：

(26)

(27)

(28)

式中，ωx、ωy、ωz为卫星的角速度.高精度的卫星角速度ωx、ωy和ωz为解算梯度V过程中的关键参数.

3 计算结果与分析

3.1 数据

考虑到GOCE卫星已有官方发布的观测数据，可作为本文构建的算法结果进行比对.本文选择GOCE卫星关键载荷的原始数据进行计算，数据长度为2013年10月8日至10月14日(共7天)，ESA(2006)提供了GOCE观测数据内容和格式的详细描述，文中所用的数据文件类型如下：EGG_NOM：加速度计控制电压数据；STR_VC2, STR_VC3：星敏感器姿态数据；AUX_EGG_DB: 重力梯度仪臂长及SSRF至GRF转换矩阵.

3.2 加速度计电压数据转换

本文针对GOCE卫星6个加速度计(A1—A6)展开分析计算，通过对电压数据转换得到加速度数据，并与ESA发布的L1b数据进行对比分析.表3为解算的加速度计A1—A6增益因子.

从表3可以看出，如果不考虑漂移，Gy、Gz的增益因子变化很小，Gx的增益因子变化则比较明显.Gx的增益因子约为(3.64～4.51)×10-4m·s-2·V-1，Gy、Gz的增益因子约为(9.67～9.89)×10-7m·s-2·V-1.以加速度计A1计算结果为例，图5为加速度计A1各个轴线性加速度与ESA发布的L1b数据中NLA_A1(Acceleration Nominal Linear, NLA)差值，图6为加速度计A1的计算值Cal_A1与ESA中NLA_A1各轴线性加速度功率谱密度(Power Spectral Density，PSD).表4给出了A1计算值与NLA_A1各个轴加速度的统计特性.

从图5可以看出，计算所得的A1线性加速度值与官方发布加速度数据，在x、y、z轴上分量ax,ay,az的差值均在非常小的范围内.图5中黑线为各轴差值的均值线，x、y、z轴对应均值分别为-8.6×10-14m·s-2、-3.2×10-11m·s-2、1.9×10-13m·s-2，各轴最大差值分别为-12.7×10-14m·s-2、-4.2×10-11m·s-2、2.5×10-13m·s-2.由于NLA_A1加速度分量ax,ay,az的量级分别为10-8m·s-2，10-6m·s-2，10-7m·s-2，其差值量级体现出很强的相似性.如图6所示为A1各轴加速度功率谱密度，表明计算结果与NLA_A1中对应数据在功率谱密度上均呈现一致性(Welch, 1967)，其中ax、az分别在低频段2～7 mHz、高频段300～500 mHz内精度稍低，其余频段各轴精度几乎相同.表4中各轴角速度的统计特性也体现了计算结果与GOCE官方发布数据在各个轴上高度相似性，其标准差均非常接近，差异极小，达到了1.1×10-14～2.1×10-12m·s-2量级，精度结果十分相近.

表3 加速度计A1—A6增益因子(单位：m·s-2·V-1)Table 3 Gain factors of accelerometers A1—A6 (unit：m·s-2·V-1)

表4 加速度计A1与NLA_A1各个轴加速度的统计Table 4 Statistics of triaxial acceleration of accelerometer A1 and NLA_A1

图5 加速度计A1三轴线性加速度差值Fig.5 Triaxial line acceleration difference of accelerometer A1

图6 A1三轴线性加速度计算结果与NLA_A1的PSD1/2对比Fig.6 PSD1/2 comparison between triaxial acceleration calculation results and NLA_A1

各个加速度计观测的加速度值通过式(2)、(3)进一步计算出共模与差分加速度，结合式(19)—(21)逆校准矩阵ICM，得到校正共模加速度(Calibrated Common Mode, CCM)和校正差分加速度(Calibrated Differential Mode, CDM).图7为重力梯度仪x、y、z轴校正后共模与差分加速度平方根功率谱密度.其中差分加速度在30～200 mHz内呈现白噪声特性，CDM-Y25与CDM-Z36在该频段范围内精度最高，达到了2×10-11、5×10-11m·s-2·Hz-1/2量级，其余各差分加速度约在10-9m·s-2·Hz-1/2量级.共模加速度呈现出随着频率升高其精度越高的特点，各轴共模加速度在测量带宽(Measurement Band Width, MBW)范围内(5～100 mHz)趋势一致且精度相似.

重力梯度仪由6个三轴加速度计组成，每个加速度计都有2个超灵敏轴(Ultra Sensitive, US)和1个低灵敏度轴(Less Sensitive, LS).由于所有6个加速度计超灵敏轴均固定安装在GRF中的x方向上(如图1)，对x轴方向共模加速度作差, 可检验计算共模加速度是否达到重力梯度仪的设计精度.由误差传播公式：

(29)

(30)

式中，σa_comb为共模加速度差值的标准差；σa_c/d_US为共模/差分加速度的标准差；σa_US为单个US轴标准差.将式(30)代入式(29)可得：

σa_comb≈σa_US.

(31)

图7 x、y、z轴上加速度计对差分(a)与共模(b)加速度PSD1/2Fig.7 x, y and z axis accelerometers couples differential mode (a) and common mode (b) acceleration PSD1/2

图8为x轴共模加速度作差后的功率谱密度，三条曲线为US轴加速度的真实噪声，在测量带宽5～100 mHz内各US轴噪声达到了10-10～10-11m·s-2·Hz-1/2范围内的良好指标，符合GOCE卫星重力梯度仪的设计要求.

图8 x轴共模加速度的差值PSD1/2Fig.8 The difference of x axis common mode acceleration PSD1/2

3.3 角速度重建

图9为GRF下多星敏感器联合前后各轴角速度平方根的功率谱密度.从图9a与图9b可以看出，与单个STR相比，多个STR联合后显示出y、z轴的角速度分量wy、wz的总体精度提升明显；在10～100 mHz处逐步下降了约1个量级，精度达到10-6rad·s-1·Hz-1/2量级；在3～30 mHz内wy的精度最高，而在30 mHz以上wx精度最高，wy与wz精度几乎相同，且wx与wy、wz均呈现出精度变化趋势相似.通过对图9a与图9b的分析可得，联合多星敏感器能有效抑制由于SSRF-GRF坐标系变换导致精度较低的角速度噪声传播到其他分量.

图9 GRF下多星敏感器联合前后各轴角速度的PSD1/2Fig.9 PSD1/2 of angular velocity of each axis before and after multiple star sensor combination under GRF

图10为基于维纳滤波的重建角速度与EGG各轴角速度平方根功率谱密度，结果表明在100 mHz频率范围内，重建角速度相较于EGG角速度的精度有明显提升.在测量带宽5～100 mHz内，角速度分量ωx、ωy、ωz平方根功率谱密度均约在32 mHz处精度改进最大，其平方根功率谱密度最大改进值范围是(5.21～6.56)×10-11rad·s-1·Hz-1/2，其中ωy的精度提升改进最大约为6.56×10-11rad·s-1·Hz-1/2.

图10 基于维纳滤波重建角速度与EGG各轴角速度的PSD1/2Fig.10 Angular velocity reconstructed based on Wiener filter and EGG angular velocity PSD1/2

3.4 重力梯度分量构建结果

式(6)计算结果中仍含有离心角速度Ω，应被减去才能得到重力梯度分量Vij；再将校准差分加速度与恢复得到的角速度代入式(23)—(28)计算得到重力梯度分量；最后基于拉普拉斯迹准则进行检验.图11为计算所得重力梯度分量与GOCE任务发布的重力梯度分量EGG_GGT (Gravity Gradients Tensor, GGT)的平方根功率谱密度分析，表5给出了重力梯度各分量计算值与GOCE任务重力梯度分量的统计特性.

从图11中可以看出，6个梯度分量的平方根功率谱密度均约在1.85×10-4Hz处达到峰值，即在该频率处的精度达到最低，Vxy约为1.7×105mE·Hz-1/2，其余分量达到了(1.5～5.0)×106mE·Hz-1/2，其原因是受到STR姿态信号噪声的影响.STR姿态噪声以每转周期(cycle per revolution, cpr)的频率(1 cpr≈1.85×10-4Hz)表现出明显的重复性，即在kcpr频率处(k为整数)6个梯度分量精度均受到的影响最大，并且在相应的频率下向较高的频率呈1/f幅度减小.图11中6个重力梯度分量中对应的平方根功率谱密度均显示出整体较为吻合，趋势较为一致.在5 mHz以上频率范围呈现出一致性，精度量级相当.在30 mHz以上频率范围，计算得到的Vxx、Vyy、Vzz、Vxz的噪声约为10 mE·Hz-1/2,Vxy、Vxz的噪声约为200 mE·Hz-1/2，明显体现了加速度计US轴和LS轴不同敏感度对重力梯度分量的影响，即Vxx、Vyy、Vzz和Vxz的测量精度较高，而Vxy和Vyz的测量精度较低.Vxx、Vyy、Vzz、Vxy对应的平方根功率谱在0.02～5 mHz低频率段内，其计算分量精度略低于EGG_GGT中对应分量(如图11a—d所示).

图11 计算所得梯度分量与L1b_GGT的PSD1/2Fig.11 PSD1/2 of the calculated gradient component and L1b_GGT

由表5可得，计算得到的Vxx与官方发布的GGT_Vxx相比较，其标准差分别为6791.54 mE、6714.88 mE，两者相差约为76.7 mE；而Vyy、Vzz、Vxy的标准差均在20.4～34.2 mE范围内.如图11e、f所示，计算得到的Vxz、Vyz与GGT对应分量在全频段内的精度呈现出极强的相似性，Vxz、Vyz与GGT对应分量的标准差分别为11.5 mE、17.5 mE.测量带宽范围内，重力梯度对角线分量在5 mHz处噪声约为1 E·Hz-1/2，0.1 Hz处约为10 mE·Hz-1/2，并在20～40 mHz范围内降低到8 mE·Hz-1/2.图12中计算得到梯度分量的迹与GGT相比，在低于10 mHz范围内受到的噪声影响更大，约低于1个量级.GOCE任务设计重力梯度张量迹在20～100 mHz内精度为11 mE·Hz-1/2，本文计算结果在该范围内约为10 mE·Hz-1/2，达到设计精度.

表5 梯度各分量计算值与L1b_GGT的统计Table 5 Calculation values of gradient components and L1b_GGT statistics

图12 计算所得梯度迹与L1b_GGT迹的PSD1/2Fig.12 PSD1/2 of calculated gradient trace and L1b_GGT trace

4 结论

本文针对重力梯度测量卫星L1级数据处理，以欧洲空间局发布的GOCE卫星数据资料为参考，初步实现了重力梯度测量卫星L1级数据全流程构建计算方法.建立了重力梯度仪观测电压数据到加速度数据转换，提出了多星敏感器联合法，有效地抑制了坐标系变换导致的精度较低角速度分量噪声传播，重建了高精度卫星姿态角速度，并最终构建了符合精度要求的重力梯度分量.

(1)通过解算电压增益因子得到的加速度观测数据与欧空局官方数据在各轴标准差达到1.1×10-14～2.1×10-12m·s-2量级，在测量带宽5～100 mHz内，共模加速度计算精度达到了梯度仪超灵敏轴10-10～10-11m·s-2·Hz-1/2的精度要求.

(2)联合多星敏感器姿态数据能有效抑制SSRF-GRF转换过程中导致的噪声传播，角速度重建精度达到10-6rad/s量级，其中分量wy、wz相比联合前约提升1个量级；基于维纳滤波重建角速度方法，将EGG与STR的角速度根据其精度在频域内进行联合，在测量带宽5～100 mHz内角速度分量wy精度最大提升了6.56×10-11rad·s-1·Hz-1/2.

(3)构建了重力梯度各分量，其中Vxx、Vyy、Vzz在5 mHz处约为1 E·Hz-1/2，0.1 Hz处约为10 mE·Hz-1/2，与欧空局发布结果的标准差相差为20.4～76.7 mE；在20～100 mHz频率范围内，重力梯度张量的迹约为10 mE·Hz-1/2，符合官方精度要求，验证了本文构建方法的有效性.

(4)本文初步系统实现了重力梯度测量卫星L1级数据构建计算方法，能够为我国未来开展重力卫星任务提供自主的数据处理技术积累和科学参考.