地球物理学报  2021, Vol. 64 Issue (12): 4437-4448   PDF    
卫星重力梯度观测数据L1级构建方法
吴云龙1,2,3, 郭泽华1,2, 肖云4, 马林5     
1. 地震大地测量重点实验室, 中国地震局地震研究所, 武汉 430071;
2. 防灾科技学院, 河北三河 065201;
3. 引力与固体潮国家野外科学观测研究站, 武汉 430071;
4. 西安测绘研究所, 西安 710054;
5. 航天东方红卫星有限公司, 北京 100094
摘要:高精度重力梯度观测数据L1级构建的系统方法是推进我国自主重力卫星任务重要的基础数据处理技术.本文以GOCE卫星L1级数据预处理技术和关键载荷原始数据为参考,面向我国发展的梯度测量卫星的任务需要,系统研究并初步实现了卫星重力梯度观测数据L1级构建方法,主要包括加速度计电压数据转换、多星敏感器联合姿态数据的角速度重建、卫星重力梯度分量构建等技术内容.计算结果表明,加速度计超灵敏轴精度为10-10~10-11m·s-2·Hz-1/2,达到重力梯度仪设计精度要求;多星敏感器联合解算最佳姿态角速度wywz在10~100 mHz内精度约提升1个量级,其精度约达到10-5 rad·s-1·Hz-1/2量级,能够有效抑制低精度角速度分量在坐标系转换中导致的噪声传播;基于维纳滤波方法恢复的角速度在5~100 mHz频段内的平方根功率谱密度提升了(5.21~6.56)×10-11 rad·s-1·Hz-1/2,显示了基于高精度角速度解算重力梯度分量的必要性;构建重力梯度各分量计算值与全球重力场和海洋环流探测器(GOCE)官方公布的重力梯度分量精度相当,其梯度张量的迹在20~100 mHz频段范围内约为10 mE·Hz-1/2,验证了本文构建方法的有效性.研究工作可为下一步我国推进实施民用重力梯度测量卫星任务提供自主的原始数据处理技术支撑与储备.
关键词: 卫星重力梯度      预处理      星敏感器      角速度重建      重力梯度数据构建     
L1 level construction method of satellite gravity gradiometry observations
WU YunLong1,2,3, GUO ZeHua1,2, XIAO Yun4, MA Lin5     
1. Key Laboratory of Earthquake Geodesy, Institute of Seismology, China Earthquake Administration, Wuhan 430071, China;
2. College of Disaster Prevention Science and Technology, Hebei Sanhe 065201, China;
3. National Observation and Research Station of Gravitation and Earth Tide, Wuhan 430071, China;
4. Xi'an Research Institute of Surveying and Mapping, Xi'an 710054, China;
5. DFH Satellite Co. Ltd, Beijing 100094, China
Abstract: The L1 level construction method of high-precision satellite gravity gradiometry observations is an important basic data processing technology for promoting the national gravity satellite mission. Base on the GOCE satellite pre-processing technology and raw data,the L1 level construction method of satellite gravity gradiometry observation is studied and initially implemented in this work,which aim for the needs of the national gravity gradiometry satellite mission. The main steps include the conversion of the accelerometer voltage data,the reconstruction of the angular rate based on the combination of the star sensors,and the construction of the gravity gradient of the satellite. The results indicate that the accuracy of the accelerometer ultra-sensitive axis is 10-10~10-11 m·s-2·Hz-1/2,which achieve the designed accuracy requirements of the gradiometer. The optimal determination of the angular rate wy,wz is improved by about 1 order of magnitude in the range of 10~100 mHz at about 10-5 rad·s-1·Hz-1/2,which is effective in suppressing noise propagation caused by low precision angular velocity components in coordinate system conversions. The square root power spectrum of the angular rate reconstruction based on the Wiener filtering method in the 5~100 mHz band is enhanced by (5.21~6.56)×10-11 rad·s-1·Hz-1/2,which indicates the necessity of solving the gravity gradient based on high accuracy angular rate. Lastly,the calculated values for each component of the gravity gradient are of comparable accuracy to the official GOCE gravity gradient components. The trace of the gradient tensor is about 10 mE·Hz-1/2 in the frequency range of 20~100 mHz,which validates our method. This work provides the technical support and reserve of independent raw data processing for promoting the national civil gravity gradiometry satellite mission.
Keywords: Satellite gravity gradiometry    Pre-process    Star sensor    Angular rate reconstruction    Gravity gradient construction    
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)

式中,下标i∈{1, 2, …, 6}分别代表图 1中加速度计A1—A6;ri为第i个加速度计与卫星质心的距离矢量;V为重力梯度张量;Ω2ri为卫星在惯性空间中旋转产生的离心加速度,为卫星角加速度产生的欧拉项,两者均反映了惯性加速度;d表示非保守力加速度,如空气阻力和太阳辐射.

图 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)

卫星在太空运行过程中,梯度仪观测数据会受到各种误差的影响,如比例因子误差、质心偏离、加速度计轴不垂直等,导致共模加速度影响差分加速度精度.为了消除以上误差影响,需对观测加速度进行校准.通过逆校准矩阵(Inverse Calibration Matrices, ICM)应用于共模和差分加速度进行校准,再通过式(7)将角加速度从差分加速度中分离出来:

(4)

(5)

利用VΩ2的对称性可得:

(6)

利用的偏对称性可得:

(7)

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

(8)

(9)

(10)

式中,分别为角加速度xyz轴的分量;LxLyLz分别为距离矢量Lxyz轴的分量.

2 L1级观测数据构建方法

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

2.1 加速度计电压数据转换

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

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

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

(11)

式中,下标i∈{1, 2, …, 6}分别代表图 1中加速度计A1—A6,为线性加速度,为8个极板的控制电压观测值,Gread, xGread, yz分别为xyz极板的转换常系数,gix1, …, giz2为加速度计各个电极对的增益因子,其表达式为:

(12)

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

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

GOCE卫星中搭载了3个星敏感器(STR1、STR2和STR3),其提供了卫星的高精度姿态信息.单个星敏感器姿态数据计算得到角速度含有较大噪声,且在其进行坐标系转换时会传递到其他角速度分量(Bandikova and Flury, 2014).星敏感器坐标系(Star Sensor Reference Frame,SSRF)中z轴角速度分量ωzSSRF精度低于xy轴角速度分量ωxSSRFωySSRF精度,由星敏感器坐标系转换到梯度仪坐标系(Gradiometer Reference Frame,GRF)时,ωzSSRF的误差将传播至GRF中角速度分量ωxGRFωyGRFωzGRF中,进而影响整个重力梯度分量的精度(Siemes et al., 2019b).

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

(13)

式中,q*为联合后得到的最佳姿态四元素;qmeasuredIRF→SSRFx分别为由惯性坐标系(Inertial Reference Frame, IRF)到SSRFx的实测姿态四元数;RGRF→SSRFx为GRF到SSRFx的旋转矩阵;分别为梯度仪坐标系中STRx的噪声;其中x∈{1, 2, 3}为不同星敏感器的标识.

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/2 Fig. 3 PSD1/2 of STR and EGG angular velocity noise

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

(14)

(15)

式中,H(f)EGGH(f)STR分别为维纳滤波中EGG和STR角速度的权重;P(f)EGGP(f)STR分别为梯度仪和星敏感器在频率f处的精度.

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

(16)

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

(17)

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

式中,hnEGGhnSTR分别为梯度仪和星敏感器角速度滤波系数.通过时域中的卷积获得恢复的角速度:

(18)

式中,ωn为重建角速度;ωnSTR为星敏感器角速度;ωnEGG为重力梯度仪角速度.

2.4 重力梯度数据构建

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

(19)

(20)

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

(21)

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

(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 加速度计A1—A6增益因子(单位:m·s-2·V-1) Table 3 Gain factors of accelerometers A1—A6 (unit: m·s-2·V-1)

表 3可以看出,如果不考虑漂移,GyGz的增益因子变化很小,Gx的增益因子变化则比较明显.Gx的增益因子约为(3.64~4.51)×10-4 m·s-2·V-1GyGz的增益因子约为(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三轴线性加速度差值 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
表 4 加速度计A1与NLA_A1各个轴加速度的统计 Table 4 Statistics of triaxial acceleration of accelerometer A1 and NLA_A1

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

各个加速度计观测的加速度值通过式(2)、(3)进一步计算出共模与差分加速度,结合式(19)—(21)逆校准矩阵ICM,得到校正共模加速度(Calibrated Common Mode, CCM)和校正差分加速度(Calibrated Differential Mode, CDM).图 7为重力梯度仪xyz轴校正后共模与差分加速度平方根功率谱密度.其中差分加速度在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)趋势一致且精度相似.

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

重力梯度仪由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)可得:

(31)

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

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

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

图 9 GRF下多星敏感器联合前后各轴角速度的PSD1/2 Fig. 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-11 rad·s-1·Hz-1/2,其中ωy的精度提升改进最大约为6.56×10-11 rad·s-1·Hz-1/2.

图 10 基于维纳滤波重建角速度与EGG各轴角速度的PSD1/2 Fig. 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 计算所得梯度分量与L1b_GGT的PSD1/2 Fig. 11 PSD1/2 of the calculated gradient component and L1b_GGT
表 5 梯度各分量计算值与L1b_GGT的统计 Table 5 Calculation values of gradient components and L1b_GGT statistics

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

表 5可得,计算得到的Vxx与官方发布的GGT_Vxx相比较,其标准差分别为6791.54 mE、6714.88 mE,两者相差约为76.7 mE;而VyyVzzVxy的标准差均在20.4~34.2 mE范围内.如图 11ef所示,计算得到的VxzVyz与GGT对应分量在全频段内的精度呈现出极强的相似性,VxzVyz与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,达到设计精度.

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

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

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

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

(3) 构建了重力梯度各分量,其中VxxVyyVzz在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级数据构建计算方法,能够为我国未来开展重力卫星任务提供自主的数据处理技术积累和科学参考.

References
Bandikova T, Flury J. 2014. Improvement of the GRACE star camera data based on the revision of the combination method. Advances in Space Research, 54(9): 1818-1827. DOI:10.1016/j.asr.2014.07.004
Baur O, Bock H, Höck E, et al. 2014. Comparison of GOCE-GPS gravity fields derived by different approaches. Journal of Geodesy, 88(10): 959-973. DOI:10.1007/s00190-014-0736-6
Becker S, Brockmann J M, Schuh W D. 2014. Mean dynamic topography estimates purely based on GOCE gravity field models and altimetry. Geophysical Research Letters, 41(6): 2063-2069. DOI:10.1002/2014GL059510
Bouman J, Fuchs M, Ivins E, et al. 2014. Antarctic outlet glacier mass change resolved at basin scale from satellite gravity gradiometry. Geophysical Research Letters, 41(16): 5919-5926. DOI:10.1002/2014GL060637
Cesare C, Sechi G, Catastini G. 2008. Gradiometer ground processing algorithms documentation. //Technical Note to ESA, GO-TN-AI-0067, Issue 7, Alenia Aerospazio.
Drewes H, Kuglitsch F, Adám J, et al. 2016. The geodesist's handbook 2016. Journal of Geodesy, 90(10): 907-1205. DOI:10.1007/s00190-016-0948-z
ESA. 1999. Gravity field and steady-state ocean circulation explorer mission. //Reports for Mission Selection. The Four Candidate Earth Explorer Core Missions, SP-1233(1). Noordwijk, The Netherlands: European Space Agency.
ESA. 2006. GOCE L1b Products user Handbook, SERCO/DATAMAT Consortium, ESA Tech. Note GOCE-GSEG-EOPG-TN-06-0137. Noordwijk: European Space Agency.
Floberghagen R, Fehringer M, Lamarre D, et al. 2011. Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission. Journal of Geodesy, 85(11): 749-758. DOI:10.1007/s00190-011-0498-3
Frommknecht B, Lamarre D, Meloni M, et al. 2011. GOCE level 1b data processing. Journal of Geodesy, 85(11): 759-775. DOI:10.1007/s00190-011-0497-4
Fuchs M J, Bouman J, Broerse T, et al. 2014. Observing coseismic gravity change from the Japan Tohoku-Oki 2011 earthquake with GOCE gravity gradiometry. Journal of Geophysical Research: Solid Earth, 118(10): 5712-5721. DOI:10.1002/jgrb.50381
Gruber T, Gerlach C, Haagmans R. 2012. Intercontinental height datum connection with GOCE and GPS-levelling data. Journal of Geodetic Science, 2(4): 270-280. DOI:10.2478/v10156-012-0001-y
Guo Z H, Wu Y L, Xiao Y, et al. 2021. Reconstruction method of satellite gravity gradient measurement angular velocity by combining star tracker quaternion. Geomatics and Information Science of Wuhan University (in Chinese), 46(9): 1336-1344. DOI:10.13203/j.whugis20200595
Hirt C. 2014. GOCE's view below the ice of Antarctica: Satellite gravimetry confirms improvements in Bedmap2 bedrock knowledge. Geophysical Research Letters, 41(14): 5021-5028. DOI:10.1002/2014GL060636
Knudsen P, Bingham R, Andersen O, et al. 2011. A global mean dynamic topography and ocean circulation estimation using a preliminary GOCE gravity model. Journal of Geodesy, 85(11): 861-879. DOI:10.1007/s00190-011-0485-8
Lu B, F rste C, Barthelmes F, et al. 2020. Using real polar ground gravimetry data to solve the GOCE polar gap problem in satellite-only gravity field recovery. Journal of Geodesy, 94(3): 34. DOI:10.1007/s00190-020-01361-z
Ning J S. 2002. The satellite gravity surveying technology and research of earth's gravity field. Journal of Geodesy and Geodynamics (in Chinese), 22(1): 1-5. DOI:10.3969/j.issn.1671-5942.2002.01.001
Pail R, Bruinsma S, Migliaccio F, et al. 2011. First GOCE gravity field models derived by three different approaches. Journal of Geodesy, 85(11): 819-843. DOI:10.1007/s00190-011-0467-x
Qu Q L, Chang X T, Zhu G B, et al. 2021. Calibration of satellite gravity gradients based on ground gravity data. Chinese Journal of Geophysics (in Chinese), 64(8): 2590-2598. DOI:10.6038/cjg2021O0463
Ran J J, Zhong M, Xu H Z, et al. 2015. Analysis of the gravity field recovery accuracy from the low-low satellite-to-satellite tracking mission. Chinese Journal of Geophysics, 58(10): 3487-3495. DOI:10.6038/cjg20151005
Rispens S, Bouman J. 2009. Calibrating the GOCE accelerations with star sensor data and a global gravity field model. Journal of Geodesy, 83(8): 737-749. DOI:10.1007/s00190-008-0290-1
Rummel R, Yi W Y, Stummer C. 2011. GOCE gravitational gradiometry. Journal of Geodesy, 85(11): 777-790. DOI:10.1007/s00190-011-0500-0
Siemes C, Haagmans R, Kern M, et al. 2012. Monitoring GOCE gradiometer calibration parameters using accelerometer and star sensor data: Methodology and first results. Journal of Geodesy, 86(8): 629-645. DOI:10.1007/s00190-012-0545-8
Siemes C. 2018a. Improving GOCE cross-track gravity gradients. Journal of Geodesy, 92(1): 33-45. DOI:10.1007/s00190-017-1042-x
Siemes C. 2018b. GOCE Level 1b gravity gradient processing algorithms. Technical Report ESA-EOPSM-GOCE-TN-3397, Paris, France: European Space Agency.
Siemes C, Rexer M, Schlicht A, et al. 2019a. GOCE gradiometer data calibration. Journal of Geodesy, 93(9): 1603-1630. DOI:10.1007/s00190-019-01271-9
Siemes C, Rexer M, Haagmans R. 2019b. GOCE star tracker attitude quaternion calibration and combination. Advances in Space Research, 63(3): 1133-1146. DOI:10.1016/j.asr.2018.10.030
Stummer C, Fecher T, Pail R. 2011. Alternative method for angular rate determination within the GOCE gradiometer processing. Journal of Geodesy, 85(9): 585-596. DOI:10.1007/s00190-011-0461-3
Stummer C, Siemes C, Pail R, et al. 2012. Upgrade of the GOCE Level 1b gradiometer processor. Advances in Space Research, 49(4): 739-752. DOI:10.1016/j.asr.2011.11.027
Sun W K, Okubo S. 2004. Coseismic deformations detectable by satellite gravity missions: A case study of Alaska (1964, 2002) and Hokkaido (2003) earthquakes in the spectral domain. Journal of Geophysical Research: Solid Earth, 109(B4): B04405. DOI:10.1029/2003JB002554
Sun W K. 2008. Progress and current situation of research on theory and observation of gravity change caused by seismicity and volcanism. Journal of Geodesy and Geodynamics (in Chinese), 28(4): 44-53, 71. DOI:10.3969/j.issn.1671-5942.2008.04.008
Tapley B D, Bettadpur S, Ries J C, et al. 2004. GRACE measurements of mass variability in the Earth system. Science, 305(5683): 503-505. DOI:10.1126/science.1099192
Visser P. 2017. Using the GOCE star trackers for validating the calibration of its accelerometers. Journal of Geodesy, 92(8): 833-846. DOI:10.1007/s00190-017-1097-8
Visser P N A M, Van Der Wal W, Schrama E J O, et al. 2014. Assessment of observing time-variable gravity from GOCE GPS and accelerometer observations. Journal of Geodesy, 88(11): 1029-1046. DOI:10.1007/s00190-014-0741-9
Visser P N A M, Van Den IJssel J A A. 2016. Calibration and validation of individual GOCE accelerometers by precise orbit determination. Journal of Geodesy, 90(1): 1-13. DOI:10.1007/s00190-015-0850-0
Wan X Y, Yu J H. 2013. Derivation of the radial gradient of the gravity based on non-full tensor satellite gravity gradients. Journal of Geodynamics, 66: 59-64. DOI:10.1016/j.jog.2013.02.005
Welch P. 1967. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2): 70-73. DOI:10.1109/TAU.1967.1161901
Wu Y L. 2010. Study on pre-processing of GOCE satellite gravity gradiometry data[Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.
Xu H Z. 2001. Satellite gravity missions-new hotpoint in geodesy. Science of Surveying and Mapping (in Chinese), 26(3): 1-3. DOI:10.3771/j.issn.1009-2307.2001.03.001
Xu X Y, Jiang W P, Zhang X M, et al. 2018. Ability of recovering the global gravity field of a new satellite gravimetry system. Chinese Journal of Geophysics (in Chinese), 61(6): 2227-2236. DOI:10.6038/cjg2018L0270
Zhong M, Duan J B, Xu H Z, et al. 2009. Trend of China land water storage redistribution at medi- and large-spatial scales in recent five years by satellite gravity observations. Chinese Science Bulletin, 54(5): 816-821. DOI:10.1007/s11434-008-0556-2
郭泽华, 吴云龙, 肖云, 等. 2021. 联合星象仪四元数的卫星重力梯度测量角速度重建方法. 武汉大学学报·信息科学版, 46(9): 1336-1344. DOI:10.13203/j.whugis20200595
宁津生. 2002. 卫星重力探测技术与地球重力场研究. 大地测量与地球动力学, 22(1): 1-5. DOI:10.3969/j.issn.1671-5942.2002.01.001
瞿庆亮, 常晓涛, 朱广彬, 等. 2021. 基于地面重力的卫星重力梯度检校方法. 地球物理学报, 64(8): 2590-2598. DOI:10.6038/cjg2021O0463
冉将军, 钟敏, 许厚泽, 等. 2015. 模拟分析低低跟踪模式重力卫星反演地球重力场的精度. 地球物理学报, 58(10): 3487-3495. DOI:10.6038/cjg20151005
孙文科. 2008. 地震火山活动产生重力变化的理论与观测研究的进展及现状. 大地测量与地球动力学, 28(4): 44-53, 71. DOI:10.3969/j.issn.1671-5942.2008.04.008
吴云龙. 2010. GOCE卫星重力梯度测量数据的预处理研究[博士论文]. 武汉: 武汉大学.
许厚泽. 2001. 卫星重力研究: 21世纪大地测量研究的新热点. 测绘科学, 26(3): 1-3. DOI:10.3771/j.issn.1009-2307.2001.03.001
徐新禹, 姜卫平, 张晓敏, 等. 2018. 一种新型重力测量卫星系统确定全球重力场的性能分析. 地球物理学报, 61(6): 2227-2236. DOI:10.6038/cjg2018L0270
钟敏, 段建宾, 许厚泽, 等. 2009. 利用卫星重力观测研究近5年中国陆地水量中长空间尺度的变化趋势. 科学通报, 54(9): 1290-1294. DOI:10.1007/s11434-008-0556-2