《测绘学 》
构建与学术的桥梁 拉近与权威的距离
梁伟1, 徐新禹1,2
1. 武汉大学测绘学院, 湖北 武汉 430079;”>
2. 武汉大学地球空间环境与大地测量教育部重点实验室, 湖北 武汉 430079;”>
3. 国家测绘地理信息局卫星测绘应用中心, 北京 100048
收稿日期：2017-05-22；修回日期：2017-10-28
基金项目：国家自然科学基金（41774020；41210006；41404020）；国家高分专项高分遥感测绘应用示范系统项目（AH1701-4）
The Determination of an Ultra-high Gravity Field Model SGG-UGM-1 by Combining EGM2008 Gravity Anomaly and GOCE Observation Data
LIANG Wei1, XU Xinyu1,2, LI Jiancheng1,2, ZHU Guangbin3
Abstract: The theory and methods of the determination of an ultra-high gravity field model by combination of satellite observation data and gravity anomaly data are studied.And an ultra-high gravity field model named SGG-UGM-1 is computed using EGM2008 derived gravity anomaly and GOCE observation data.The block-diagonal least squares (BDLS) method for quickly estimating an ultra-high gravity field model is researched and the corresponding software module is validated by numerical experiments.OpenMP technique is introduced into the BDLS program, which improves computing efficiency dramatically.An ultra-high gravity model SGG-UGM-1 complete to degree and order 2159 is derived using the proposed calculation strategies.The fully occupied normal equation system up to degree and order 220 formed by GOCE satellite data and the block-diagonal normal equation system up to degree and order 2159 formed by EGM2008 gravity anomaly data are used for the combination.Comparison among the models SGG-UGM-1 and EGM2008, EIGEN-6C2, EIGEN-6C4, GOSG-EGM in frequency domain has been done, which shows that SGG-UGM-1 is close to the reference models and that the coefficients lower than degree 220 of SGG-UGM-1 are more accurate than that of EGM2008.The models are also validated by GPS-leveling data in China and America and airborne gravity data in Maowusu surveying area.The results show that in China the accuracy level of SGG-UGM-1 derived geoid is between EIGEN-6C2 and EIGEN-6C4, and better than GOSG-EGM and EGM2008, while in America they are almost the same.In Maowusu area the accuracy level of SGG-UGM-1 derived gravity disturbance is almost at the same accuracy level with EGM2008 and EIGEN-6C4 and better than GOSG-EGM and EIGEN-6C2.
Key words: SGG-UGM-1 ultra-high gravity field model block diagonal least squares method OpenMP parallel computing technique
构建超高阶地球重力场模型一直是物理大地测量学领域的热点研究问题。采用重力梯度测量技术的GOCE重力卫星于2009年成功发射，仅利用其观测数据恢复的全球重力场模型阶次达到了280阶[1]，模型中长波分量的精度更是有了2个数量级的提高，尤其在地面重力数据空白区，相比由GRACE和地面观测数据联合反演的EIGEN-5C等模型，在100 km空间分辨率上，GOCE任务反演的模型精度有了明显提升[2]。这一进展使得联合GOCE卫星观测数据和地面观测数据获得全球更高精度的超高阶重力场模型成为可能。反演的超高阶重力场模型可用于区域高精度大地水准面的确定[3]、全球高程基准的统一[4]、地球内部结构反演[5]等科学研究及应用。
基于球面重力异常数据反演超高阶全球重力场模型常用的方法主要有数值积分(numerical quadrature，NQ)方法、最小二乘配置(least squares collocation，LSC)方法和最小二乘(least squares，LS)方法[6]。NQ方法是基于面球谐函数的正交性，采用积分方法独立求解模型的每个系数，其计算简单，但无法直接给出模型系数的验后方差。相比NQ方法，LS方法不仅可评估模型系数的精度，而且可以联合解算多类观测数据。同时，在LS方法中，如果位系数按照特定的顺序排列，法方程矩阵会具有块对角特点，因此可以分块求解位系数，这种方法称为块对角最小二乘(block diagonal least squares，BDLS)方法[6]。该方法可避免超高维矩阵的构建及求逆，计算量大大减少，已成为目前求解超高阶重力场模型的主要方法，国际上精度得到认可的EGM96[7]、EGM2008[8]、EIGEN-6C4[9]等模型在构建时均使用了该方法。文献[10]分析了BDLS方法在重力场模型确定中的应用，并通过模拟试验比较了BDLS方法与NQ方法求解360阶重力场模型的精度。
考虑到卫星重力观测数据和重力异常数据在反演全球重力场模型时存在频谱互补性，前者对重力场的中长波信号更为敏感，后者对短波及甚短波信号更为敏感，因此联合卫星重力数据和重力异常数据可以反演得到高精度高分辨率的重力场模型，这也是当前求解超高阶重力场模型的主要方法。EGM2008模型构建时联合了GRACE卫星观测法方程和重力异常数据[8]，EIGEN-6C4模型不仅使用了这些数据，还使用了GOCE观测数据，相比EGM2008模型在一些区域的精度有明显提高[9]。文献[11]联合EGM2008模型和GOCE-TIM-R5模型解算了2160阶次重力场模型GECO；文献[12]采用全矩阵求逆的最小二乘方法解算了720阶次的GOCO05c模型，模型精度在15′的空间分辨率上(720阶次)与EIGEN-6C4相当；文献[13]基于卫星、陆地、海洋等观测数据，使用球谐分析方法确定了2159阶次的UGM08模型；文献[14]采用谱权组合方法联合EIGEN-6S卫星重力场模型和格 重力异常构建了2160阶次的重力场模型。
1 基本原理1.1 由重力异常数据快速求解模型系数的BDLS方法
格 平均重力异常观测值的观测方程为[6]
(1)
式中，Δσi=Δλ(cos θi-cosθi+1)；
；(i,j)表示第i行j列个格 ；(r,θ,λ)为球坐标[15]；GM为地球引力常数；R为参考球半径；rij为r在格 上的平均值；Pnm(cosθ)为规格化的缔合勒让德函数；IPnmi为Pnm(cosθ)sinθ的积分值；(Cnm,Cnm)为扰动位球谐系数；(n,m)为勒让德函数及模型系数的阶次；Δgij为格 重力异常平均值。
由式(1)，可以建立如下误差方程
式中，L是观测值Δgij；v是观测值的改正数；A是基于式(1)建立的设计矩阵；
在LS平差准则下，可得式(3)，其中P是观测值Δgij的权阵，式(3)即为使用LS方法求解重力场模型的基本公式
(3)
使用LS方法构建法方程时，当重力异常数据及其精度的空间分布满足下面4个条件时，法方程为稀疏矩阵，如按照一定顺序排列待求解参数，法方程矩阵N为块对角形式[16]：
(1) 格 数据分布于旋转曲面上(例如球面上)；
(2) 格 数据覆盖整个曲面且经度方向的分辨率一致；
(3) 格 数据的精度与经度无关，即数据的权重不依赖于经度的变化；
(4) 格 数据的精度关于赤道对称，即纬度θ与-θ处的数据精度相同。
例如当系数先按次排列，然后同次的系数中Cnm系数在前，Snm系数在后，同次的Cnm和Snm均按阶排列，则法方程将是严格的块对角结构。例如：6阶重力场模型系数的法方程矩阵结构如图 1所示。
图 1 法方程矩阵的结构示意图Fig. 1 The schematic diagram of the normal matrix灰色部分代表非零元素，其余为零元素
图选项
如果法方程矩阵为块对角形式，求解参数时可以不构建和求解整个法方程，而是根据法方程具体的分块结构，单独构建法方程矩阵N中的每个块对角部分及与其对应的U，形成子法方程，对其单独求解，所有子法方程解的组合与原法方程的解完全一致。这种求解参数的方法称为块对角最小二乘(BDLS)方法。求解超高阶重力场模型时，参数个数较多，法方程庞大，现有的计算条件很难满足严格LS方法的计算需求，而BDLS方法可以避免大型矩阵的构建及求逆，计算量远小于LS方法。
1.2 多源重力观测数据的最小二乘联合求解
可以分别根据观测数据L1和L2的观测模型建立它们的误差方程如式(4)所示
(4)
式中，v1、v2为观测量的改正数；A1、A2为设计矩阵、
在最小二乘准则下可以建立数据L1和L2的法方程如式(5)所示，求解法方程可以得到待求参数
(5)
式中，P1、P2分别为数据L1和L2的权阵。
同时也可联合数据L1和L2，建立它们联合观测的误差方程如下
式中，vC=[v1v2]T；AC=[A1A2]T；lC=[l1l2]T。
在最小二乘准则下可以建立它们的联合观测法方程，如式(7)所示
(7)
式中，PC为数据的权阵。
当数据L1和L2不相关时，式(7)中权阵PC=
，其中w1和w2分别为数据L1和L2的加权因子，此时有[18-19]
(8)
由式(8)可知，当数据L1和L2不相关时，使用最小二乘联合平差方法联合数据L1和L2构建的法方程和直接将各类观测数据构建的法方程按照加权因子相加形成的法方程相同[18-19]。即当观测数据不相关时，直接把法方程加权相加就可实现多源数据的最小二乘联合处理，这种方法简单有效。
在最小二乘联合平差中的一种特殊情况是数据L1和L2求解的参数x1和x2可能不完全一样，它们构建的法方程的参数和维数不同，这样不能直接将法方程相加。此时可以分别由数据L1和L2构建参数x1和x2的参数集合xC的误差方程。此时建立的参数xC的误差方程与原来参数x1和x2的误差方程的区别在于设计矩阵中增加了0元素，这些0元素是与观测量无关的参数的系数。进而建立参数xC的法方程，此时将法方程相加可以得到它们对参数xC的联合观测法方程。
2 基于OpenMP技术的块对角最小二乘软件模块设计及验证
并行计算是指同时使用多种计算资源完成计算任务，是缩短计算时间的有效手段。常用的并行计算有OpenMP多线程技术和MPI多进程技术。OpenMP可为程序中的并行计算区域开辟多个线程，每个线程独立运行，多个线程协同完成并行计算区域的计算任务[20-21]。
使用BDLS方法求解超高阶重力场模型时，需要依次构建法方程矩阵N中的块对角部分和U中与其对应的部分，并对其单独求解。法方程的每个块对角部分的构建和求解可单独进行，在程序设计时一般用DO循环完成。OpenMP技术可以为相互独立的DO循环开辟多个线程，将多个循环分配到多个线程上运行，充分利用高性能计算平台的计算资源，以提高计算效率。
图 2 模型的系数阶误差RMSFig. 2 Degree error RMS of the model coefficients
图选项
3 SGG-UGM-1超高阶重力场模型的求解
3.1 超高阶重力场模型解算的数据处理流程
图 3 SGG-UGM-1模型的解算流程Fig. 3 Calculation process of SGG-UGM-1 model
图选项
(1) 使用EGM2008模型计算参考球面(球半径为6 378 136.3 m)上5′×5′格 分辨率的重力异常数据ΔgEGM；由该数据使用BDLS方法构建并求解2159阶块对角法方程，得到2~2159阶位系数，这部分系数命名为BDLS，使用BDLS方法构建法方程时认为格 数据的权相等。
(2) 根据步骤(1)求得系数中2~100阶以及251~2159阶系数计算重力异常，然后在数据ΔgEGM中减去这部分重力异常得到残余重力异常数据ΔgR。ΔgR数据中仅有101~250阶系数的贡献，再使用ΔgR数据构建101~250阶系数的块对角法方程。
(3) 使用GOCE卫星的卫星重力梯度(satellite gravity gradiometry，SGG)数据和卫星跟踪卫星(satellite-to-satellite tracking，SST)数据构建220阶卫星观测法方程，并求解得到220阶的卫星重力场模型，其数据处理方法的描述见下节。
(4) 步骤(2)建立了101~250阶系数的法方程，步骤(3)建立了2~220阶系数的法方程，使用1.2节描述的联合平差方法，联合这两个法方程得到2~250阶系数的联合观测法方程。3.3节中有关于法方程联合过程更为详细的描述。
(5) 求解步骤(4)建立的2~250阶系数联合观测法方程，可以得到2~250阶系数，这部分系数与步骤(1)中求解的系数中251~2159阶部分一同组成2~2159阶重力场模型系数。
3.2 GOCE卫星观测法方程的构建
基于上面的观测数据，建立GOCE卫星观测法方程的简要描述如下：首先，对原始数据进行粗差探测与剔除、时变重力场信号的分离、数据内插分段等预处理；直接利用沿轨的引力梯度对角线数据(Vxx，Vyy，Vzz)建立观测方程及相应的法方程，再使用方差分量估计方法估计引力梯度张量对角线分量各自的单位权中误差，3个分量的单位权中误差分别为2.5、3.3和4.3 mE。考虑到引力梯度观测值的有色噪声特点，在最小二乘过程中引入通带为0.005~0.041 Hz的ARMA(auto regressive moving-average)滤波器。采用点域加速度方法由PKI(precise kinematic orbit)轨道观测值建立相应的SST法方程，求解SST法方程得到观测值残差，并基于残差估计了单位权方差。基于SGG各分量的单位权方差以及SST的单位权方差，得到了它们之间的加权因子，最后基于加权因子将SGG数据各分量的法方程以及SST法方程相加得到GOCE卫星观测法方程，并求解得到相应的卫星重力场模型，求解时采用Tikhonov正则化方法解决极空白问题引起的法方程不稳定性问题。文献[18]中有更为详细的数据处理方法和过程的描述。
3.3 GOCE卫星观测法方程和格 重力异常法方程的联合解算
图 4 SGG-UGM-1的法方程联合解算示意图Fig. 4 The schematic diagram of normal equation combination in SGG-UGM-1
图选项
联合观测法方程分为两部分，其中251~2159阶系数法方程与重力异常法方程中对应部分相同，而2~250阶系数法方程是使用1.2节中的最小二乘联合平差方法联合101~250阶重力异常数据块对角法方程和GOCE全阶次法方程(220阶)得到的。构建重力异常数据101~250阶块对角法方程时需要在原始重力异常数据中移除2~100阶及251~2159阶系数的贡献。
模型的220~250阶系数是联合GOCE数据和重力异常数据联合求解的，这样可以保证联合模型的这部分系数及其验后方差的频谱有更好的连续性。重力异常数据求解低阶系数的精度较差，重力异常低阶次系数法方程与卫星法方程联合可能会“污染”卫星重力的低阶系数，所以为减弱其影响，2~250阶联合法方程中只采用了重力异常的101~250阶系数块对角法方程。
4 SGG-UGM-1超高阶重力场模型的精度分析4.1 SGG-UGM-1模型的系数误差
图 5 重力场模型的系数阶误差RMSFig. 5 Degree error RMS of the gravity field models
图选项
由图 5可知，在2~70阶的频率范围内，EIGEN-6C2的误差最小，EGM2008次之，而SGG-UGM-1和GOSG-EGM较大。这是由于EIGEN-6C2与EIGEN-6C4的2~70阶系数的计算方法和使用数据类似，都是联合LAGEOS、GRACE和GOCE等数据求解的[9, 22]，EGM2008的2~70阶次系数主要是来自于GRACE观测数据的贡献；而SGG-UGM-1和GOSG-EGM的2~70阶次系数主要由GOCE数据求解，因此其与EIGEN-6C4模型的差异大于EGM2008、EIGEN-6C2。
在70~220阶，SGG-UGM-1和EIGEN-6C2与EIGEN-6C4更为接近，主要原因是两个模型在此频段内均有GOCE和地面重力数据的贡献，而GOSG-EGM仅有GOCE卫星的贡献，从170阶开始它的误差迅速增大。同时，在170~220阶部分，因为SGG-UGM-1是联合GOCE和地面重力数据解算的，因此比GOSG-EGM模型的系数误差平滑，在220阶没有出现突变现象。
图 6 重力场模型的大地水准面阶误差Fig. 6 Degree geoid error of the gravity field models
图选项
4.2 中国和美国GPS/水准数据的检核结果
表 1 中国区域的GPS/水准数据检核结果Tab. 1 Validation of the gravity field models using GPS-leveling data in China
| m | ||||
| 模型 | 最大值 | 最小值 | 平均值 | 标准差 |
| EGM2008 | 1.729 | -1.535 | 0.239 | 0.240 |
| GOSG-EGM | 0.704 | -0.793 | 0.240 | 0.169 |
| SGG-UGM-1 | 0.744 | -0.618 | 0.246 | 0.162 |
| EIGEN-6C2 | 0.808 | -0.592 | 0.243 | 0.167 |
| EIGEN-6C4 | 0.729 | -0.698 | 0.243 | 0.157 |
表选项
表 2 美国区域的GPS/水准数据检核结果Tab. 2 Validation of the gravity field models using GPS-leveling data in America
| m | ||||
| 模型 | 最大值 | 最小值 | 平均值 | 标准差 |
| EGM2008 | 0.360 | -1.396 | -0.511 | 0.284 |
| GOSG-EGM | 0.392 | -1.470 | -0.511 | 0.291 |
| SGG-UGM-1 | 0.317 | -1.407 | -0.511 | 0.280 |
| EIGEN-6C2 | 0.398 | -1.415 | -0.512 | 0.283 |
| EIGEN-6C4 | 0.397 | -1.392 | -0.512 | 0.282 |
表选项
如表 1所示，在中国区域SGG-UGM-1大地水准面的精度远优于EGM2008，略优于EIGEN-6C2和GOSG-EGM模型，略逊于EIGEN-6C4模型。如表 2所示，在美国区域SGG-UGM-1模型的精度最高，但其他几个模型的精度基本相当。SGG-UGM-1模型在美国和中国区域的精度均优于EGM2008模型，这主要是因为SGG-UGM-1模型使用了GOCE数据，GOCE数据提高了模型低阶次系数的精度，这与上文的模型系数分析的结果一致。尤其与EGM2008模型在地面重力数据空白或稀疏区对比，SGG-UGM-1模型在这些区域(如中国大陆)的精度提升更为明显[25]。SGG-UGM-1模型的大地水准面在中国区域和美国区域的精度检核结果均优于GOSG-EGM，充分体现了最小二乘联合平差方法的优势。
4.3 毛乌素测区航空重力数据的检核结果
表 3 航空重力数据的检核结果Tab. 3 Validation of the gravity field models using airborne gravity data
| mGal | ||||
| 模型 | 最大值 | 最小值 | 平均值 | 标准差 |
| EGM2008 | 4.694 | -25.565 | -10.206 | 5.556 |
| GOSG-EGM | 3.089 | -28.550 | -9.883 | 6.506 |
| SGG-UGM-1 | 5.046 | -26.715 | -10.095 | 5.792 |
| EIGEN-6C2 | 7.779 | -23.414 | -7.562 | 5.851 |
| EIGEN-6C4 | 4.179 | -25.227 | -10.118 | 5.554 |
表选项
如表 3所示，在毛乌素测区，SGG-UGM-1模型计算重力扰动的精度略逊于EGM2008和EIGEN-6C4模型，优于GOSG-EGM和EIGEN-6C2模型，总体来说几个模型的精度相当。相比GPS水准检核的结果，虽然SGG-UGM-1模型的低阶次系数相较EGM2008模型有了较大提升，但由于其高阶次系数主要是由EGM2008模型重力异常求解的，且重力扰动的误差主要来自高阶次系数，所以SGG-UGM-1航空重力扰动的精度相较EGM2008没有提升。
5 总结
(1) BDLS方法可以快速精确地确定超高阶重力场模型，采用OpenMP技术可以大幅度提高BDLS软件模块的计算效率。
(2) 采用最小二乘联合平差方法不仅可以保证求解模型在联合频段内信号和误差频谱的连续性，同时求解的模型SGG-UGM-1精度优于直接拼接获得的GOSG-EGM模型。
【引文格式】梁伟, 徐新禹, 李建成, 等. 联合EGM2008模型重力异常和GOCE观测数据构建超高阶地球重力场模型SGG-UGM-1[J]. 测绘学 ，2018，47(4)：425-434. DOI:
10.11947/j.AGCS.2018.20170269
往期精彩回顾
“2017年度中国遥感领域十大事件”揭晓
数学能力是中兴与华为的唯一区别
发表中文论文就很low？再谈“中文期刊”与“文化自信”
李德仁院士：老师教我做人做学问
权威 | 专业 | 学术 | 前沿
微信投稿邮箱 | song_qi_fan@163.com
进群请备注：姓名+单位+稿件编号
声明：本站部分文章内容及图片转载于互联 、内容不代表本站观点，如有内容涉及侵权，请您立即联系本站处理，非常感谢！