地球物理学进展  2014, Vol. 29 Issue (1): 46-50   PDF    
引用本文
刘晓刚, 闫志闯, 孙文, 周睿. 确定地球重力场模型的最小二乘配置法与调和分析法的精度评析[J]. 地球物理学进展, 2014, 29(1): 46-50, doi: 10.6038/pg20140106  
LIU Xiao-gang. Precision evaluation and analysis of least squares collocation method and spherical harmonic analysis method in the determination of the Earth’s gravity field model. Progress in Geophysics, 2014, 29(1): 46-50, doi: 10.6038/pg20140106   
确定地球重力场模型的最小二乘配置法与调和分析法的精度评析
刘晓刚1,2,3, 闫志闯4, 孙文4, 周睿4    
1. 西安测绘研究所, 西安 710054;
2. 地理信息工程国家重点实验室, 西安 710054;
3. 武汉大学地球空间环境与大地测量教育部重点实验室, 武汉 430079;
4. 信息工程大学地理空间信息学院, 郑州 450052
收稿日期: 2013-8-18 修回日期: 2013-12-12 .
基金项目：国家自然科学基金资助项目（41304022， 41174026， 41104047）、国家973资助项目（61322201， 2013CB733303）和地球空间环境和大地测量教育部重点实验室开放基金资助项目（11-01-03）联合资助.
作者简介：刘晓刚，男，1983年生，汉族，博士，助理研究员，甘肃天水人，主要从事卫星重力、磁力测量研究.（E-mail：liuxiaogang_1949@163.com）
摘要：本文给出了利用卫星重力梯度张量Tzz分量恢复地球重力场模型的最小二乘配置法与调和分析法的基本原理公式， 通过数值实验对两种方法确定地球重力场模型的精度进行了比较， 最后总结评述了两种方法在反演地球重力场模型时的优缺点.
关键词： 地球重力场模型     最小二乘配置法     调和分析法     卫星重力梯度     GOCE    
Precision evaluation and analysis of least squares collocation method and spherical harmonic analysis method in the determination of the Earth’s gravity field model
LIU Xiao-gang1,2,3,YAN Zhi-chuang4,SUN Wen4,ZHOU Rui    
1. Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China;
2. State Key Laboratory of Geo-Information Engineering, Xi’an 710054, China;
3. Key laboratory of Geo-space Environment and Geodesy of Ministry of Education, Wuhan University, Wuhan 430079, China;
4. Institute of Geospatial Information, Information Engineering University, Zhengzhou 450052, China
Abstract: Fundamental formulae of least squares collocation(LSC) method and spherical harmonic analysis method to recover the Earth’s gravity field model(EGM) from Tzz component of satellite gravity gradients(SGG) are deduced, the precisions of these two methods in the determination of the EGM are compared by numerical experiments, and the advantage and disadvantage of these two methods in the recovery of the EGM are summarized and reviewed at last.
Key words: Earth’s gravity field model(EGM)     least squares collocation(LSC)     spherical harmonic analysis     satellite gravity gradient(SGG)     GOCE    

0 引 言

自从卫星重力梯度(SGG：Satellite Gravity Gradient)技术提出以来，国内外众多学者在利用SGG数据恢复地球重力场模型的调和分析法与最小二乘配置法(LSC：Least Squares Collocation)研究方面，已经开展了大量卓有成效的工作.

(Colombo，1981)将快速傅里叶变换(FFT：Fast Fourier Transform)算法引入到全球重力异常的调和分析中，并建立了全球重力场模型化理论；(吴晓平，1991)利用调和分析方法讨论了重力场元在地面和空间的谱分布特征和向下延拓问题，分析了各类数据求定重力场的最高分辨率和精度；(罗志才，1996，1997)研究了利用SGG数据确定地球重力场模型的调和分析方法，导出了由重力梯度张量的球函数展开系数确定扰动位球谐系数的实用解算模型；(张传定，2000)完善了LSC理论，提出了最小二乘复配置理论，得到了各类单定边值问题的解析解、调和分析解、最小二乘解和最小二乘复配置解，建立了超定边值问题的最小方差解、最小二乘解和最小二乘复配置解；(Tscherning，2001)研究了利用LSC法估计全球重力场位系数及其误差估计的方法，分析了不同数据类型(大地水准面高、重力异常、引力位二阶梯度)以及不同轨道类型和数据长度对求解位系数的影响，并指出使用“移去-恢复”方法可以得到更好的结果；(Migliaccio et al，2004a，2004b)提出了在空域法中采用Wiener轨道滤波器对有色噪声导致的沿轨时间序列相关的梯度观测数据进行滤波预处理，给出了一种基于调和分析解算位系数的迭代解法；(李迎春，2004)研究了利用SGG张量求解扰动位系数的调和分析方法，导出了实用解算数学模型，给出了迭代最小二乘方法的数学模型和解算途径，利用引力位的协方差函数导出了位系数与引力位二阶导数的协方差函数模型，并推导出LSC法的实用解算公式；(张传定等，2004)利用球面到“轮胎”面的映射关系，将高精度调和分析步骤分三步走模式，克服了传统调和分析中平滑因子不精确等引起的误差；(张嗥等，2005)通过空间扰动位协方差函数特性，得出SGG数据与扰动位系数的相关协方差函数，利用LSC法，推导出由SGG数据直接解算引力位系数的函数表达式；(Kotsakis，2007)研究了在尽量不改变精度的前提下，通过消除LSC法内在的平滑效应来增强其恢复地球重力场能力的方法；(徐新禹，2008；钟波，2010)基于调和分析方法推导出卫星轨道面扰动位T和梯度Tzz分量的谱组合所对应谱权的具体形式；(吴星, 2009; 吴星等，2009a，2009b)研究了广义轮胎调和分析方法、点质量调和分析方法和LSC调和分析方法，并通过数值实验验证了这些方法的有效性和实用性；(吴星等，2010)推导了扰动引力梯度数据与球谐系数之间的协方差阵和自协方差阵，并有效利用FFT技术将其降阶，得到了引力梯度径向分量LSC调和分析的完整计算公式；(刘晓刚, 2011; 刘晓刚等，2012a，2012b)建立了利用扰动重力梯度张量Tzz分量和Txx+Tyy、Tzz-Txx-Tyy组合分量确定地球重力场的调和分析法模型，进一步推导了扰动重力梯度张量对角线三分量的自协方差和互协方差函数的级数展开式，推导了单分量、组合分量与重力位系数之间协方差函数的实用计算公式，给出了利用单分量和组合分量解算地球重力场模型的LSC法基本原理公式.

从上述文献可以看出，虽然国内外众多学者对LSC法与调和分析法恢复地球重力场模型的理论开展了大量的研究，但对于两种方法之间优劣性的比较，目前尚未见诸于文献报道，因此，本文将以SGG张量Tzz分量为例，给出其恢复地球重力场模型的LSC法与调和分析法的基本原理公式，并采用模拟数据实验结果对两种方法的优缺点做出评价，从而为读者在选择不同的地球重力场模型确定方法时提供一定的参考.

1 最小二乘配置法的基本原理

利用SGG张量Tzz分量解算地球重力场模型的LSC法的基本原理公式为(罗志才，1996；张传定，2000；李迎春，2004；吴星等，2009b，2010；刘晓刚等，2012b)

其中， J表示位系数的列向量，C J,Tzz表示位系数与 T zz分量的协方差矩阵，其计算公式为

其中，fM表示地球引力常量；R是地球平均半径；r′=R+h是卫星轨道上任意一点的地心向径，h表示卫星高度；θ、λ分别表示地心余纬和地心经度；(C *nm,S nm)表示完全正常化地球扰动引力位系数；n和m分别表示球谐系数的阶和次；P nm(cosθ)表示完全正常化缔合勒让德函数.CJ,Tzz表示Tzz分量之间的协方差矩阵，其计算公式为

T zz表示全球规则平均扰动重力梯度张量Tzz分量形成的列向量.

2 调和分析法的基本原理

利用SGG张量Tzz分量解算地球重力场模型的调和分析法基本原理公式为(罗志才等，1997；张传定，2000；李迎春，2004；徐新禹，2008；吴星等，2009b；钟波，2010；刘晓刚, 2011; 刘晓刚等， 2012a，2012b)

3 两种方法的精度比较

在具体数值实验中，利用EGM2008重力场模型生成了GOCE卫星轨道上的SGG张量Tzz分量数据(刘晓刚等，2010a)，通过延拓(刘晓刚等，2011)和格网化(刘晓刚等，2010b)方法将重力梯度数据归算为250 km平均轨道高度处45′×45′的格网平均值，并对两极数据空白区进行填充，得到全球范围内的重力梯度格网平均值数据.

首先，向45′×45′的全球SGG张量Tzz分量格网平均值数据中分别加入标准差为0.0 mE、0.1 mE、1.0 mE和10.0 mE的零均值白噪声，利用LSC法来解算重力场模型，其对应的重力场模型阶误差RMS、累计大地水准面和重力异常误差如下图 1所示，下表 1表示加入不同误差的Tzz数据恢复240阶重力场模型对应的精度统计结果.

图 1 LSC法恢复的重力场模型阶误差RMS、累计大地水准面和重力异常误差 (a)重力场模型阶误差RMS；(b)累计大地水准面误差；(c)累计重力异常误差. Fig. 1 The degree error RMS of the EGM, the cumulative geoid error and the cumulative gravity anomaly error recovered by LSC method (a)Degree error RMS of the EGM; (b)Cumulative geoid error; (c)Cumulative gravity anomaly error.

表 1 LSC法恢复240阶重力场模型的精度统计结果 Table 1 The precision statistical results of the EGM in 240 degree recovered by LSC method

其次，为了比较调和分析法与LSC法解算结果的精度，利用含有噪声的45′×45′全球SGG张量Tzz分量格网平均值数据，基于调和分析法恢复地球重力场模型，其对应的重力场模型阶误差RMS、累计大地水准面和重力异常误差如下图 2所示，下表 2表示加入不同误差的Tzz数据恢复240阶重力场模型对应的精度统计结果.

图 2 调和分析法恢复的重力场模型阶误差RMS、累计大地水准面和重力异常误差 (a)重力场模型阶误差RMS；(b)累计大地水准面误差；(c)累计重力异常误差. Fig. 2 The degree error RMS of the EGM, the cumulative geoid error and the cumulative gravity anomaly error recovered by spherical harmonic analysis method (a)Degree error RMS of the EGM; (b)Cumulative geoid error; (c)Cumulative gravity anomaly error.

表 2 调和分析法恢复240阶重力场模型的精度统计结果 Table 2 2 The precision statistical results of the EGM in 240 degree recovered by spherical harmonic analysis method

从图 1和表 1可以看出，含有10.0 mE误差的数据对LSC法解算结果的精度影响较大，而含有0.1 mE和1.0 mE误差的数据，其解算结果在2～180阶次内与不含误差的数据结果比较一致，这说明LSC法在解算含有误差的数据时具有一定的优势.

从图 2和表 2可以看出，当观测数据含有不同的误差时对解算结果的精度影响较大，随着加入误差的不断增大，其恢复的重力场模型有效阶次不断降低，精度也不断下降.在Tzz数据中分别加入0.1 mE、1 mE和10 mE的零均值白噪声时，其恢复的重力场模型有效阶数分别约为225、180和120阶.

通过将图 1、表 1与图 2、表 2的结果比较来看，调和分析法的解算结果精度要优于LSC法，这主要是由下列因素引起的：

1)当协方差矩阵存在病态性问题时，其求逆也带来了一定的计算误差，即使采用了正则化方法也不能完全避免.

2)配置方法是重力场数据处理的重要手段，但它毕竟是一种统计近似，只能给出信号的估值，而不能给出其真值.

3)本文的全球协方差模型是采用直至2160阶的EGM2008模型计算所得，与实际的情况也存在一定的差异，因此，先验协方差模型的不准确性，也导致计算结果精度的下降.

4 结束语

LSC法与调和分析法的共同点是：

在具体实现时需要首先将卫星轨道高度处的SGG数据归算到同一轨道高度上，接着，需要将归算后的观测数据进行格网化处理，以便得到规则的格网点值或平均值，因此，两种方法的缺陷是不可避免地存在观测数据的归算和格网化误差；如果采用FFT技术，两种方法都可以实现地球重力场模型的快速求解；此外，由于GOCE卫星观测在两极有球冠半径约6.7°的数据空白区，因此，两极数据的缺失会导致所恢复的地球重力场模型精度受到影响，因此，两种方法需要采用其它重力测量手段获得的或参考重力场模型的计算数据来对两极空白区进行填充，从而尽可能降低极区空白问题的影响.

LSC法与调和分析法的不同点是：

1)LSC法优点是建立的法方程都是呈块状对角形式，因此方程的求逆等计算将非常简便、快捷；该方法也可以综合各类观测数据确定地球重力场模型，在数据处理时能顾及观测量的误差，并且能够给出待估量的误差协方差信息.该方法的缺点是要求有可靠的先验协方差模型，在恢复高阶次地球重力场模型时，需要解算超大型协方差矩阵的逆矩阵，由于大型矩阵求逆存在的病态性问题，从而影响了数值解的稳定性.

2)调和分析方法的优点是将每个位系数表示为扰动引力梯度张量各分量的全球积分形式，因此位系数可以单个求解，并且其数学模型严密，计算量小，也比较容易实现.该方法的缺点是存在积分离散化误差，也不能提供用于精度评估的参数误差协方差信息.

参考文献
[1] Colombo O L. 1981. Numerical Methods for Harmonic Analysis on the Sphere[R].   Columbus: Department of Geodetic Science and Surveying, The Ohio State University.
[2] Kotsakis C. 2007. Least-squares collocation with covariance-matching constrains[J].   Journal of Geodesy, 81(10): 661-677.
[3] Li Y C. 2004. Theories and method of the recovery of the Earth’s gravity field using satellite gravity gradients[D].   Zhengzhou: The PLA Information Engineering University.
[4] Liu X G, Wu X P, Wang B J, et al. 2010a. Study on gridding methods of satellite gravity gradient data[J].   Journal of Geodesy and Geodynamics (in Chinese), 30(6): 60-65.
[5] Liu X G, Wu X P, Zhao D M, et al. 2010b. Non-singular expression of the disturbing gravity gradients[J].   Geodaetica et Cartographica Sinica (in Chinese), 39(5): 450-457.
[6] Liu X G. 2011. Theory and methods of the Earth’s gravity field model recovery from GOCE data(in Chinese) [D].   Zhengzhou: The PLA Information Engineering University.
[7] Liu X G, Li S S, Wu X. 2011. Downward continuation of satellite gravity gradient data[J].   Journal of Geodesy and Geodynamics (in Chinese), 31(1): 132-137.
[8] Liu X G, Zhang L P, Wang J, et al. 2012a. Establishment of harmonic analysis models of single component and composite components of disturbed gravity gradient tensor[J].   Journal of Geodesy and Geodynamics (in Chinese), 32(5): 103-108, 112.
[9] Liu X G, Wu X P, Wang K. 2012b. Construction of the least squares collocation models for single component and composite components of disturbed gravity gradients[J]. Chinese J. Geophys.  , 55(5): 1572-1580.
[10] Luo Z C. 1996. Theories and method of the determination of the earth’s gravity field using satellite gravity gradients(in Chinese)[D]. Wuhan: Wuhan Technical University of Surveying and Mapping.
[11] Luo Z C, Chao D B, Ning J S. 1997. Spherical harmonic analysis on gravity gradient tensors[J].   Journal of Wuhan Technical University of Surveying and Mapping (in Chinese), 22(4): 346-349.
[12] Migliaccio F, Reguzzoni M, Sansò F. 2004a. Space-wise approach to satellite gravity field determination in the presence of coloured noise[J].   Journal of Geodesy (in Chinese), 78: 304-313.
[13] Migliaccio F, Reguzzoni M, Sansò F, et al. 2004b. GOCE: dealing with large attitude variations in the conceptual structure of the space-wise approach[C]. // Proc.   Second International GOCE User Workshop “GOCE, The Geoid and Oceanography”, ESA-ESRIN, Frascati, Italy.
[14] Tscherning C C. 2001. Computation of spherical harmonic coefficients and their error estimates using least-squares collocation[J].   Journal of Geodesy (in Chinese), 75(1): 12-18.
[15] Wu X. 2009. Research of methods of spherical harmonic analysis of the Earth’s gravity field(in Chinese)[D]. Zhengzhou: Information Engineering University.
[16] Wu X, Zhang C D, Zhao D M. 2009a. Generalized torus harmonic analysis of satellite gravity gradients component[J].   Geodaetica et Cartographica Sinica (in Chinese), 38(2): 101-107.
[17] Wu X, Zhang C D, Zhao D M. 2009b. Point-mass harmonic analysis method based on spherical boundary value problem[J]. Chinese J. Geophys.   (in Chinese), 52(12): 2993-3000.
[18] Wu X, Zhang C D, Liu X G. 2010. Least squares collocation harmonic analysis of the radial satellite gravity gradients[J].   Geodaetica et Cartographica Sinica (in Chinese), 39(5): 471-477, 483.
[19] Wu X P. 1991. Some analyses of the determination of the Earth’s gravity field from the satellite gravity gradient and terrestrial data[J].   Journal of the PLA Institute of Surveying and Mapping (in Chinese), (1): 19-32.
[20] Xu X Y. 2008. Study of determining the Earth’s gravity field from satellite gravity gradient and satellite-to-satellite tracking data (in Chinese)[D]. Wuhan: Wuhan University.
[21] Zhang C D. 2000. Satellite gravimetry: foundation, modeling methods, and data processing algorithms (in Chinese)[D]. Zhengzhou: Information Engineering University.
[22] Zhang C D, Xu H Z, Wu X. 2004. On the torus approach for gobal spherical harmonic computation[C].   Hubei: Hubei Science and Technology Publishing House, 302-314.
[23] Zhang H, Chen Q. 2005. Least squares collocation in calculating potential coefficient with satellite gravity gradiometry data[J].   Geomatics and Spatial Information Technology (in Chinese), 28(6): 34-35.
[24] Zhong B. 2010. Study on the determination of the Earth’s gravity field from satellite gravimetry mission GOCE[D].   Wuhan: Wuhan University.
[25] 李迎春. 2004. 利用卫星重力梯度测量数据恢复地球重力场的理论与方法(in Chinese)[D].   郑州: 解放军信息工程大学.
[26] 刘晓刚, 吴晓平, 王宝军,等. 2010a. 卫星重力梯度数据的网格化方法研究[J].   大地测量与地球动力学, 30(6): 60-65.
[27] 刘晓刚, 吴晓平, 赵东明,等. 2010b. 扰动重力梯度的非奇异表示[J].   测绘学报, 39(5): 450-457.
[28] 刘晓刚. 2011. GOCE卫星测量恢复地球重力场模型的理论与方法[D].   郑州: 解放军信息工程大学.
[29] 刘晓刚, 李珊珊, 吴星. 2011. 卫星重力梯度数据的向下延拓[J].   大地测量与地球动力学, 31(1): 132-137.
[30] 刘晓刚, 张丽萍, 王俊,等. 2012a. 扰动重力梯度张量单分量和组合分量调和分析法模型的建立[J].   大地测量与地球动力学, 32(5): 103-108, 112.
[31] 刘晓刚, 吴晓平, 王凯. 2012b. 扰动重力梯度张量单分量和组合分量最小二乘配置法模型的建立[J].   地球物理学报, 55(5): 1572-1580.
[32] 罗志才. 1996. 利用卫星重力梯度数据确定地球重力场的理论与方法[D]. 武汉: 武汉测绘科技大学.
[33] 罗志才, 晁定波, 宁津生. 1997. 重力梯度张量的球谐分析[J].   武汉测绘科技大学学报, 22(4): 346-349.
[34] 吴星. 2009. 卫星重力梯度数据处理理论与方法[D].   郑州: 信息工程大学.
[35] 吴星, 张传定, 赵东明. 2009a. 卫星重力梯度分量的广义轮胎调和分析[J].   测绘学报, 38(2): 101-107.
[36] 吴星, 张传定, 赵东明. 2009b. 基于球面边值问题的点质量调和分析方法[J].   地球物理学报, 52(12): 2993-3000.
[37] 吴星, 张传定, 刘晓刚. 2010. 卫星重力径向梯度数据的最小二乘配置调和分析[J].   测绘学报, 39(5): 471-477, 483.
[38] 吴晓平. 1991. 利用卫星重力梯度与地面数据确定地球重力场的分析[J].   解放军测绘学院学报, (1): 19-32.
[39] 徐新禹. 2008. 利用卫星重力梯度及卫星跟踪卫星数据确定地球重力场的研究[D]. 武汉: 武汉大学测绘学院.
[40] 张传定. 2000. 卫星重力测量-基础、模型化方法与数据处理算法[D]. 郑州: 信息工程大学.
[41] 张传定, 许厚泽, 吴星. 2004. 地球重力场调和分析中的“轮胎”问题[C].   湖北: 湖北科学技术出版社, 302-314.
[42] 张皞, 陈琼. 2005. 卫星重力梯度数据解算位系数的最小二乘配置法[J].   测绘与空间地理信息, 28(6): 34-35.
[43] 钟波. 2010. 基于GOCE卫星重力测量技术确定地球重力场的研究[D].   武汉: 武汉大学.