卫星重力梯度恢复地球重力场的频谱分析

蔡 林，周泽兵，祝 竺，高 防，许厚泽

1 华中科技大学电子与信息工程系，武汉 430074

2 华中科技大学物理学院，武汉 430074

3 华中科技大学地球物理研究所，武汉 430074

4 中国科学院测量与地球物理研究所，武汉 430077

卫星重力梯度恢复地球重力场的频谱分析

蔡 林1，3，周泽兵2，3＊，祝 竺2，3，高 防2，3，许厚泽3，4

1 华中科技大学电子与信息工程系，武汉 430074

2 华中科技大学物理学院，武汉 430074

3 华中科技大学地球物理研究所，武汉 430074

4 中国科学院测量与地球物理研究所，武汉 430077

根据仪器功率谱密度和重力位系数阶方差的定义，本文建立了卫星重力梯度测量噪声功率谱密度与重力场模型的误差阶方差的直接对应关系，并基于此讨论了重力梯度测量精度、卫星轨道高度以及运行时间对地球重力场恢复精度的影响.相比于传统的基于最小二乘法评估卫星载荷噪声对地球重力场恢复精度的影响而言，本文提出的方法简单、直接，有助于快速设计和确定卫星重力测量计划的有关参数.

卫星重力梯度测量，地球重力场，频谱分析

1 引 言

卫星重力测量技术公认为是当今获取全球高精度高分辨率地球重力场及其时变信息最有效的技术手段之一.于2000年7月15日和2002年3月17日分别发射升空的CHAMP卫星［1］和GRACE卫星［2］采用的是卫卫跟踪（SST）观测技术，测量了地球重力场的中长波分量.2009年3月17日成功发射的GOCE卫星［3－4］结合了高低卫－卫跟踪模式和卫星重力梯度测量（SGG）两种技术模式.相对于CHAMP和GRACE而言，GOCE能更为精确地测量地球重力场的中短波分量，对地球的位系数测量能达到200阶以上.

地球重力场的恢复精度取决于卫星的测量模式、卫星高度、轨道倾角、卫星载荷的测量精度等参数.目前在卫卫跟踪（SST）测量的数据处理和误差分析中广泛采用的是利用基于轨道摄动分析理论的时域法，包括动力学积分法［5－6］、能量守恒法［7－8］和短弧边值法［9］等.对于卫星重力梯度数据的处理和评估，一般是基于将重力梯度观测值作为卫星位置或者时间的函数，在此基础上利用相应的空域最小二乘法［10－11］或时域最小二乘法［10，12－13］来确定重力场模型并进行误差分析.以上这些方法本质上都是利用最小二乘法评估重力场恢复精度［14－15］，能够在得到重力位系数的同时给出相应重力场恢复的精度，以及可以把不同测量数据组合进行联合求解，得到优化的结果.这类方法不仅需要消耗大量计算资源，且难以直接评估某些单一参数的影响.另外，当前和未来的重力场研究都面临着解算模型的分辨率越来越高、阶数越来越大的问题.如果使用传统方法研究误差分析和恢复能力，则在阶数较大时，每增加一阶都会使计算量急剧增长并带来严重的数值不稳定性问题［16］.所以，高效、稳定的评估重力场精度的方法对高分辨率重力场模型的研究有着重要意义.

为了快速评估和确定卫星重力测量计划有关参数，以及评估不同类型载荷噪声功率谱的影响，本文根据仪器功率谱密度和重力位模型阶方差的定义，讨论了重力梯度测量噪声与地球重力场的频谱关系，建立起地球重力场恢复精度与仪器噪声水平的直接表达式，这有助于项目设计初期选取和确定相应载荷指标，以及分析部分载荷指标调整的影响.最后，本文利用该方法对卫星重力梯度测量中仪器技术指标、卫星高度和运行时间与恢复地球重力场的精度之间的关系作了模拟分析.

2 理论分析

地球引力场扰动位函数T的球谐展开如下［10］：

其中GM为地心引力常数，R为地球平均半径，l和m分别为球谐展开的阶和次，与为相对于正常重力位的扰动球谐系数，（cosθ）为归一化的勒让德函数，r、θ和λ分别为地心参考系下观测点的半径、余纬和经度，r＝R＋H，H表示卫星轨道高度.

对（1）式求导并利用相关坐标关系，可得扰动位函数T在局部指北坐标系下沿z方向的梯度分量Tzz为

从上式可看出，若能精确测定重力梯度张量分量Tzz，则可以恢复高精度的地球重力场，得到相应的位系数和.同时（2）式也意味着，和的不确定度与取决于Tzz的测量噪声.根据球谐函数的正交性，由（2）式得到Tzz的测量噪声总功率为［17］

如果已知梯度测量噪声功率谱密度Szz（f）（单位：E2／Hz），则第l阶m次对应频带引入的噪声功率为

其中flm为第l阶m次重力位系数对应的测量频率，由傅里叶级数和球谐系数的对应关系来确定［18－19］，Δf为flm对应的测量频率间隔，取决于测量频率的分辨率和信号频率的稳定度.重力梯度测量水平或者分辨率取决于测量噪声的功率谱密度，工程上习惯用噪声功率谱密度的开方来表示，其单位为E／Hz1／2.若忽略实际重力卫星无法覆盖极轨区域的影响，同时假定测量频带内梯度测量噪声为白噪声（例如对于GOCE计划而言，在测量频带内重力梯度仪噪声谱密度可看作是白噪声）.在此假设条件下，再考虑到l阶内独立球谐系数¯Clm和¯Slm的个数为（2l＋1），可以得出第l阶对应频带内引入的梯度测量噪声功率为

把（5）和（6）式带入（4）式，整理得到

当卫星的运行时间为T时，忽略重力信号的时变影响（即忽略信号频率的涨落影响），则每阶次位系数所对应的梯度测量的频率间隔为Δf＝1／T［13］，若梯度测量噪声在测量频带内为白噪声，对每阶次的影响都为第l阶的总贡献为由（7）式可得到：

为了分析重力位系数的恢复精度对应的大地水准面起伏和重力异常，我们使用大地水准面累积误差σN和重力异常累积误差σΔg评定其影响，表达式分别为［10］

和

3 结果分析

根据（8）、（9）和（10）式，我们对不同的梯度测量噪声功率谱密度、卫星高度和卫星运行时间的影响进行分析和讨论.

3.1 梯度测量精度影响分析

当卫星运行时间和高度一定时，重力场恢复精度主要依赖于重力梯度的测量精度.重力梯度测量精度依赖于梯度仪以及卫星平台角速度的测量精度.由（8）式可知，若其它因素不变时，重力梯度测量精度提高N倍，重力场的恢复精度会相应提高N倍.若卫星轨道高度H为250km，运行时间T为6个月，梯度测量水平分别为4、7和10mE／Hz1／2，得到相应的误差阶方差估计如图1所示.需要说明的是，GOCE计划建议立项前最先确定梯度测量精度为4mE／Hz1／2［3］，GOCE计划第三次研讨会将其调整为7mE／Hz1／2［20］，GOCE卫星在轨运行时重力梯度测量水平约为10mE／Hz1／2甚至更差［20］.另外，本文利用Kaula准则评定恢复重力场位系数的最大阶数［21］，其具体表达式如下：

图1 不同梯度测量噪声情况下的误差阶方差估计Fig.1 Errors degree variance with different precision of gravity gradient measurements

从图1中可以看到，基于Kaula准则，梯度测量水平为4、7和10mE／Hz1／2所对应的恢复重力场位系数的最大阶数分别为245、231和221阶，这表明提高重力梯度的测量精度可有效提高恢复地球重力场的精度.

相应的大地水准面累积误差估计如图2a所示，当阶数为100阶时，梯度测量水平为4、7和10mE／Hz1／2所对应的大地水准面精度分别约为0.22、0.39和0.56cm.从150阶左右开始大地水准面精度随阶数升高而明显地变差，当阶数为200阶时，对应的大地水准面精度分别为2.11、3.69和5.27cm.相应的重力异常累积误差估计如图2b所示，梯度测量水平为4、7和10mE／Hz1／2，对应前50阶的重力异常累积误差分别为0.0038、0.0067和0.0095mGal，前200阶的累积误差分别为0.59、1.04和1.48mGal.表1给出不同梯度测量噪声情况下大地水准面累积误差以及重力异常累积误差与恢复阶数的关系.

表1 不同梯度测量噪声情况下的大地水准面累积误差以及重力异常累积误差估计比较Table 1 Cumulative geoid height errors and gravity anomaly errors with different precision of gravity gradient measurements

图2 不同梯度测量噪声情况下的大地水准面累积误差（a）以及重力异常累积误差（b）估计Fig.2 Cumulative geoid height errors（a）and cumulative gravity anomaly errors（b）with different precision of gravity gradient measurements

3.2 卫星轨道高度影响分析

轨道高度是重力卫星的关键参数之一，其选择直接影响到卫星重力恢复精度.当卫星运行时间和重力梯度测量精度一定时，不考虑大气密度变化等因素时，降低轨道高度有利于提高恢复重力场的精度.假设在重力梯度测量卫星的梯度测量噪声不变的情况下，梯度测量水平为7mE／Hz1／2，运行时间为6个月，卫星轨道高度H分别为200、250和300km，得到相应的误差阶方差估计如图3所示，重力场位系数的最大恢复阶次数分别为287、230和191阶，即轨道高度每降低50km，恢复重力场系数的最大阶数提高约40～50阶.因此，降低轨道高度可以大幅提高地球重力场恢复的精度和空间分辨率.

相应的大地水准面累积误差估计如图4a所示，当阶数小于100阶时，恢复的大地水准面精度在不同卫星轨道高度情况下差异不大，分别为0.31、0.39和0.58cm，这意味着降低卫星轨道高度对阶数小于100阶的大地水准面精度的恢复影响较小.随着阶数增加，降低卫星轨道高度对大地水准面精度的恢复有非常大的影响，当阶数为200阶时，对应的大地水准面精度分别为0.96、3.69和15.16cm，分别相差约4倍.相应的重力异常累积误差估计如图4b所示，从中可得重力异常累积误差估计随卫星轨道高度增加而变差.由重力异常累积误差估计同样可以得到，阶数较低时，卫星轨道高度对重力异常的作用不明显，当轨道高度为200、250和300km时，前50阶的重力异常累积误差分别为0.0050、0.0067和0.0090mGal；而阶数较高时，卫星轨道高度的作用随阶数增加越来越明显，前200阶的重力异常累积误差分别为0.25、1.04和4.35mGal，也是分别相差约4倍.表2给出了不同高度情况下的大地水准面累积误差以及重力异常累积误差估计与恢复阶数的关系.由此可看出，在卫星实际运行条件下，降低卫星运行轨道高度有利于提高地球重力场中短波段系数的恢复精度和空间分辨率.

图3 不同高度情况下的误差阶方差估计Fig.3 Errors degree variance with different altitudes

3.3 卫星运行时间影响分析

重力场的恢复精度还取决于卫星的运行时间.由（8）式可知，其他因素不变时，卫星运行时间延长N倍，静态重力场的恢复精度会提高倍.当然，考虑到重力场的时变效应，需选择合适的运行时间进行模拟分析.若卫星轨道高度为250km，梯度测量水平为7mE／Hz1／2，运行时间T为1、3和6个月得到相应的误差阶方差估计如图5所示，对应恢复重力场位系数的最大阶数分别为206、222和230阶，表明增加运行时间可提高恢复地球重力场的精度.

表2 不同高度情况下的大地水准面累积误差以及重力异常累积误差估计比较Table 2 Cumulative geoid height errors and cumulative gravity anomaly errors with different altitudes

相应的大地水准面累积误差估计和重力异常累积误差估计如图6a和6b所示，从中可以看出，当阶数为200阶时，运行时间为1、3和6个月所对应的大地水准面精度分别为9.04、5.22和3.69cm，重力异常累积误差估计分别为2.54、1.47和1.04mGal.大地水准面和重力异常恢复精度随着卫星运行时间增加而变好，表3给出了多种运行时间情况下的大地水准面累积误差以及重力异常累积误差估计，当重力梯度测量噪声和卫星轨道高度一定时，重力场恢复的空间分辨率依赖于卫星的运行时间.值得强调的是，若考虑到实际情况中重力场的时变效应，卫星运行时间进一步增加对地球重力场恢复精度和空间分辨率的改善不会太显著.

图4 不同高度情况下的大地水准面累积误差（a）以及重力异常累积误差（b）估计Fig.4 Cumulative geoid height errors（a）and cumulative gravity anomaly errors（b）with different altitudes

表3 不同运行时间下的大地水准面累积误差以及重力异常累积误差估计比较Table 3 Cumulative geoid height errors and cumulative gravity anomaly errors with different operation duration

图5 不同运行时间下的误差阶方差估计Fig.5 Errors degree variance with different mission duration

4 讨论和结论

根据仪器功率谱密度和重力场模型阶方差的定义，本文直接建立了梯度测量噪声功率谱对位系数恢复精度的表达式，并在此基础上讨论了卫星重力梯度测量中仪器技术指标、卫星轨道高度和卫星运行时间对重力场恢复的影响.结果分析表明，重力梯度测量精度越高、卫星轨道高度越低以及运行时间越长，恢复地球重力场精度越高.其结果与国外学者用数值仿真方法得到的结果符合很好［10，22－23］.相对于传统的研究卫星重力测量误差分析和恢复能力论证的方法，本文提出的方法具有简单和直接等优点，且不需要消耗大量计算资源.该方法特别适合项目建议初期对部分观测量或者载荷指标进行快速评估和确定，也可以对卫星重力测量系统数值模拟仿真中参数选择和优化设计提供指导.

需要指出的是，在本文的推导中采用的是极圆轨道，而实际重力测量卫星的轨道是偏心率很小的椭圆，且轨道的倾角略大于或小于90°，这会导致数据采集不能覆盖极区以及恢复的重力场模型精度与理想情况有一定偏差.此外，本文仅使用Vzz进行估计，未考虑Vxx和Vyy以及重力梯度卫星中卫卫高低跟踪模式的贡献.本文讨论中，假设了测量频带范围内梯度测量噪声功率谱密度与频率无关，即测量频带内白噪声的假设.下一步将对这些问题做进一步分析，特别是实际梯度测量中色噪声的影响分析.

图6 不同运行时间下的大地水准面累积误差（a）以及重力异常累积误差（b）估计Fig.6 Cumulative geoid height errors（a）cumulative gravity anomaly errors（b）with different operation duration

致 谢 感谢中国科学院测量与地球物理研究所钟敏研究员、武汉大学测绘学院姜卫平教授和徐新禹博士的有益讨论，感谢审稿专家提出的修改意见和建议.

（References）

［1］ Reigber C，Balmino G，Schwintzer P，et al.A high－quality global gravity field model from CHAMP GPS tracking data and accelerometry（EIGEN－1S）.Geophysical Research Letters，2002，29：1692.

［2］ Jekeli C.The determination of gravitational potential differences from satellite－to－satellite tracking.Celestial Mechanics and Dynamical Astronomy，1999，75（2）：85－101.

［3］ European Space Agency.Gravity field and steady－state ocean Circulation Exploration Reports for Mission Selection.The Four Candidate Earth Explorer Core Missions，SP－1233（1）.European Space Agency，1999.

［4］ Rummel R，Colombo O L.Gravity field determination from satellite gradiometry.Bulletin Geodesique，1985，59（3）：233－246.

［5］ Reigber C.Gravity field recovery from satellite tracking data.∥SansòF，Rummel R eds.Theory of Satellite Geodesy and Gravity Field Determination.Lecture Notes in Earth Sciences，Vol.25.New York：Springer－Verlag，1989.

［6］ 周旭华，许厚泽，吴斌等.用GRACE卫星跟踪数据反演地球重力场.地球物理学报，2006，49（3）：718－723.Zhou X H，Hsu H，Wu B，et al.Earth＇s gravity field derived from GRACE satellite tracking data.Chinese J.Geophys.（in Chinese），2006，49（3）：718－723.

［7］ Han S C，Jekeli C，Shum C K.Efficient gravity field recovery using in situ disturbing potential observables from CHAMP.Geophysical Research Letters，2002，29：1789.

［8］ 徐天河，杨元喜.基于能量守恒方法恢复CHAMP重力场模型.测绘学报，2005，34（1）：1－6.Xu T H，Yang Y X.CHAMP gravity field recovery using energy conservation.Acta Geodaetica et Cartographic Sinica（in Chinese），2005，34（1）：1－6.

［9］ Mayer－Güerr T，Ilk K H，Eicker A，et al.ITG－CHAMP01：A CHAMP gravity field model from short kinematic arcs over a one－year observation period.Journal of Geodesy，2005，78（7－8）：462－480.

［10］ Rummel R，van Gelderen M，Koop R，et al.Spherical Harmonic Analysis of Satellite Gradiometry.Publications on Geodesy，New Series 39，Netherlands Geodetic Commission，1993.

［11］ Schuh W D.Tailored Numerical Solution Strategies for the Global Determination of the Earth＇s Gravity Field.Technical Report，Graz：Technical University Graz，1996.

［12］ Schrama E J O.Gravity field error analysis：Applications of GPS receivers and gradiometers on low orbiting platforms，NASA TM－100769，GSFC Greenbelt Md.20771，1990.

［13］ Sneeuw N.A semi－analytical approach to gravity field analysis from satellite observations.Munich：Technical University of Munich，2000.

［14］ Tscherning C C.Computation of spherical harmonic coefficients and their error estimates using least－squares collocation.Journal of Geodesy，2001，75（1）：12－18.

［15］ Baur O，Austen G，Kusche J.Efficient GOCE satellite gravity field recovery based on least－squares using QR decomposition.Journal of Geodesy，2008，82（4－5）：207－221.

［16］ Gruber C，Novák P，Sebera J.FFT－based high－performance spherical harmonic transformation.Studia Geophysica et Geodaetica，2011，55（3）：489－500.

［17］ Ilk K H，Flury J，Rummel R，et al.Mass transport and mass distribution in the earth system.Proposal for a German Priority Research Program，GeoForschungs Zentrum Potsdam.2004.

［18］ Kazhdan M，Funkhouser T，Rusinkiewicz S.Rotation invariant spherical harmonic representation of 3Dshape descriptors.∥Proceedings of the 2003Eurographics／ACM SIGGRAPH Symposium on Geometry Processing.Aachen，Germany：Eurographics Association Press，2003：156－16.

［19］ Novotni M，Klein R.3DZernike descriptors for content based shape retrieval.∥Proceedings of the Eighth ACM Symposium on Solid Modeling and Applications.Washington，USA：ACM Press，2003：216－225.

［20］ Drinkwater M R，Haagmans R，Muzi D，et al.The GOCE Gravity Mission：ESA′s First Core Earth Explorer.Proceedings of 3rd International GOCE User Workshop，6－8 November，2006，Frascati，Italy，ESA SP－627，ISBN 92－9092－938－3，2007：1－8.

［21］ Kaula W M.Theory of Satellite Geodesy.Massachusetts：Blaisdell Publishing Company Press，1966：34－35.

［22］ Koop R.Global gravity field modelling using satellite gravity gradiometry.Publications on Geodesy，New Series 38，Netherlands Geodetic Commission，1993.

［23］ Schrama E J O.Gravity field error analysis：applications of global positioning system receivers and gradiometers on low orbiting platforms.Journal of Geophysical Research，1991，96（B12）：20041－20051.

Spectral analysis for recovering the Earth′s gravity potential by satellite gravity gradients

CAI Lin1，3，ZHOU Ze－Bing2，3＊，ZHU Zhu2，3，GAO Fang2，3，HSU Houtse3，4
1 Department of Electronics and Information Engineering，Huazhong University of Science and Technology，Wuhan 430074，China
2 School of Physics，Huazhong University of Science and Technology，Wuhan 430074，China
3 Institute of Geophysics，Huazhong University of Science and Technology，Wuhan 430074，China
4 Institute of Geodesy and Geophysics，Chinese Academy of Sciences，Wuhan 430077，China

Direct relationship between the power spectrum density of the satellite gradiometry and the coefficients of the Earth′s gravity potential is established based on the definitions of the instrument′s power spectrum density and the Earth′s gravity field potential，and then the effect of gravity gradient measurement accuracy，the altitude of the satellite，and the operation duration on recovering of Earth＇s gravity field is analyzed by using this method.This method is simple and direct compared with the others based on the least square approach，and is very useful to design and verify the parameters for the Earth gravity recovery missions.

Satellite gravity gradiometry，The Earth＇s gravity potential，Spectral analysis

10.6038／j.issn.0001－5733.2012.05.014

P228

2011－10－18，2011－11－27收修定稿

国家自然科学基金项目（11075058）资助.

蔡林，男，1984年生，在读博士研究生，主要从事卫星重力测量方法研究.E－mail：mumu0202＠gmail.com

＊通讯作者周泽兵，男，1973年生，教授，主要从事卫星重力测量和引力实验研究.E－mail：zhouzb＠mail.hust.edu.cn

蔡林，周泽兵，祝竺等.卫星重力梯度恢复地球重力场的频谱分析.地球物理学报，2012，55（5）：1565－1571，

10.6038／j.issn.0001－5733.2012.05.014.

Cai L，Zhou Z B，Zhu Z，et al.Spectral analysis for recovering the Earth′s gravity potential by satellite gravity gradients.Chinese J.Geophys.（in Chinese），2012，55（5）：1565－1571，doi：10.6038／j.issn.0001－5733.2012.05.014.

（本文编辑 胡素芳）