航空重力测量数据向下延拓的改进Poisson积分迭代法

刘晓刚，孙中苗，管 斌，范昊鹏,4

1. 地理信息工程国家重点实验室，陕西 西安 710054； 2. 西安测绘研究所，陕西 西安 710054； 3. 大地测量与地球动力学国家重点实验室，湖北 武汉 430077； 4. 信息工程大学地理空间信息学院，河南 郑州 450052

航空重力测量技术，因其可以在沙漠、沼泽、冰川、原始森林、陆海交界等一些难以开展地面重力测量的区域进行作业，快速经济地获取精度良好、分布均匀、大面积的地球重力场中高频信息，从而成为地球重力场研究的最为热门的领域之一[1-2]。

航空重力测量技术主要包括标量测量技术和矢量测量技术。航空重力标量测量技术已经非常成熟，而航空重力矢量测量技术正在国家高分辨率对地观测系统重大专项的支持下开展研究。我国的航空重力矢量测量仪初样机已经研制成功，其集成试验在经历了飞机选型论证及加改装、试验区勘选及测量、系统安装集成及测试标定等准备工作后，2015年6—10月，西安测绘研究所组织国防科技大学、第一测绘导航基地等单位在内蒙古地区进行了我国首次航空重力矢量测量仪集成飞行试验，获得了第一手的航空和地面重力矢量测量数据。2017年6月，西安测绘研究所组织上述单位在山西地区进行了我国第二次航空重力矢量测量仪集成飞行试验。因此，这些原始测量数据的预处理及其专业处理的成果，是评判我国首套航空重力矢量测量仪系统的研制是否满足要求的重要标准，而向下延拓则是其专业处理中非常重要的一个环节。

航空重力测量数据向下延拓结果的好坏，直接影响到其进一步应用，如不同类型重力测量数据的融合、全球或区域地球重力场模型的构建、全球或区域(似)大地水准面的精化、水下重力匹配辅助导航中重力基准图的生成等。将空中重力测量数据向下延拓到地球表面采用的是逆Poisson积分方程，属于典型的不适定问题[3-6]。为了减弱方程的病态性对于延拓结果的影响，国内外学者提出了很多解决思路，如最小二乘配置法[7-11]、直接代表法[12]、迭代法[13]、解析法[14-16]和正则化法[17-26]等。

目前使用的大多数向下延拓方法都是在Poisson积分方程的基础上发展而来的，而局部区域空中格网重力测量数据的向下延拓，需要先将Poisson积分进行离散化处理，从而达到离散求和的目的。本文在对传统向下延拓模型进行数值分析的基础上，提出了Poisson积分迭代法这一延拓思路，并与传统最小二乘法、改进最小二乘法、Tikhonov正则化法等延拓模型进行了精度比较。

1 传统最小二乘法延拓模型

根据Poisson积分公式，将地面重力异常数据向上延拓，有[7]

(1)

由于在实际情况下空中和地面重力异常均是离散值，式(1)可写成如下形式的观测方程

gh=AX

(2)

式中，gh为空中重力异常Δgh所组成的向量；X为地面区域的重力异常Δg0所组成的向量；A为相应的系数阵。在航空重力测量中，空中重力异常Δgh是已知观测值，地面重力异常为待估参数，这是Poisson积分方程的逆问题。

式(2)的传统离散化形式为

(3)

式中，矩阵Aij的表达式为

(4)

对系数矩阵A作如下奇异值分解

A=UΛVT

(5)

式中，U、V分别由矩阵A的左、右特征向量组成；U=[u1u2…un]，V=[v1v2…vn]，U∈Rn×n,V∈Rn×n，且UTU=UUT=I，VTV=VVT=I；Λ∈Rn×n为对角阵，其对角元素是矩阵A的奇异值λ1、λ2、…、λn，且λi按递减的次序排列。

考虑到式(5)，式(2)的最小二乘解及其谱分解形式为[19]

(6)

式中，λi为A的奇异值；ui、vi分别为其相应的左右奇异向量；rank(A)表示矩阵A的秩；而l=gh为误差方程的自由项。

式(6)即为传统的最小二乘法延拓模型。

2 改进最小二乘法延拓模型

采用式(4)所示的Poisson积分离散化形式在某些情况下不能得到正确的结果，甚至误差很大。主要原因是，与实际结果相比，处于计算点向径方向的流动点对计算点的贡献较大，如果以该点的核函数值作为整个格网的核函数值，将会带来较大误差。因此，有学者给出了修改向径方向流动点加权系数的离散化形式。矩阵A的对角线元素为[27-30]

(7)

矩阵A的非对角线元素为

(8)

式中，ψ0表示选取积分半径；NC表示积分半径内的测量点个数。

该离散化公式的推导思路是先将积分半径内的其他格网对计算点的延拓贡献值求出，再将整个积分范围内的理论贡献值与之相减，即得到向径方向的流动点对计算点的贡献。

将式(6)中的矩阵A用式(7)和式(8)来替换，就可以得到改进的最小二乘法延拓模型。

3 Tikhonov正则化法延拓模型

航空重力测量数据向下延拓时，观测方程是病态的，其系数矩阵A的奇异值单调地趋向零，由式(6)难以获得稳定的解。因此，需对病态方程进行正则化处理，以抑制高频测量噪声对延拓结果的影响。在众多正则化方法中，以Tikhonov正则化法的应用最为广泛，所采用的准则如下

(9)

根据式(9)的约束条件，可得Tikhonov正则化解为

(10)

Tikhonov正则化解是有偏的，其偏差为

(11)

Tikhonov正则化解的平均均方误差为

MSE(Xα)=D(Xα)+bias(Xα)bias(Xα)T=

(12)

将式(5)代入式(10)、(12)，可得Tikhonov正则化解的谱分解式为[19]

(13)

通过实际计算表明，采用式(13)来获得Tikhonov正则化解的速度要优于式(10)，这是因为谱分解式避免了许多大型矩阵的乘法和求逆运算，节约了大量计算时间。

(14)

上式表明，Tikhonov正则化解的误差由两部分组成：第一项为观测误差所引起的估值误差，随α的增大单调减小；第二项为正则化所引起的估值误差，随α的增大单调增大，当α=0时，该项误差为0。选择合适的正则化参数α可起到平衡这两项误差的作用，使得它们的和为最小，这正是通过正则化算法可获得精确、稳定解的原因。

4 改进Poisson积分迭代法延拓模型

如何正确地确定正则化参数，则是整个正则化算法的关键，参数选择的优劣直接影响了最终延拓结果的精度。在模拟试验中，最佳正则化参数的确定，是在延拓面数据已知的情况下根据最小延拓误差来估算的。而对于实际的航空重力测量数据向下延拓，延拓面并不一定存在已知的测量数据。因此，在这种情况下如何确定正则化参数，存在较大的难度。为了减弱延拓模型的病态性，提高延拓结果的精度，本文提出Poisson积分迭代法这一延拓思路，进而建立相应的延拓模型。

由于R/r也是调和函数，因此有

(15)

在式(5)两端同时乘以Δg0(R,θ,λ)，有

(16)

将式(1)减去式(16)，通过变换可以得到

(17)

对式(17)进行改化，可以得到

(18)

由于式(18)两端都含有Δg0(R,θ,λ)项，因此，该式可以进一步改化成迭代形式。作为第一次近似，命

(19)

再用下式进行迭代计算

sinθ′dθ′dλ′

(20)

式(20)则为本文建立的Poisson积分迭代法向下延拓模型，该式可以进一步修改成矩阵计算的形式

(21)

式中，矩阵A可以用式(4)来计算(本文称之为新建模型Ⅰ，即Poisson积分迭代法向下延拓模型)，也可以用式(7)和式(8)来计算(新建模型Ⅱ，即改进Poisson积分迭代法向下延拓模型)。当前后两次计算值之间的差异小于ε时，就可以终止迭代。

5 数值试验与结果分析

5.1 试验数据说明

采用中国某地区实测航空重力测量数据，分辨率是5′×5′，飞行高度为3400 m。其等值线图如图1所示。

相应的地面重力测量数据，其等值线图如图2 所示。

航空和地面重力测量数据的统计结果如表1所示。

表1 数据统计结果

图1 航空重力测量数据等值线图(单位：mGal)Fig.1 Isoline chart of airborne gravimetry data (Unit：mGal)

图2 地面重力测量数据等值线图(单位：mGal)Fig.2 Isoline chart of ground gravimetry data (Unit：mGal)

根据数字高程模型SRTM30，采用Global Mapper软件获得了该地区5′×5′的地形数据，其等值线图如图3所示。

图3 地形数据等值线图(单位：m)Fig.3 Isoline chart of topography data(Unit：m)

地面重力测量数据对应的地形数据统计结果如表2所示。

表2 试验区地形数据结果统计

以地面点的平均地形高度1 476.115 m作为向下延拓的基准面，则延拓高度是1 923.885 m。

根据图1、图2与图3的对比可以看出，航空和地面重力数据的变化趋势，与地形的起伏密切相关。地形起伏较为剧烈的地区，重力数据也有相应的变化。

5.2 对比试验结果

首先，利用传统的最小二乘法延拓模型，将航空重力测量数据向下延拓到地面平均高程面上，其等值线图如图4所示。将延拓结果与地面重力测量数据进行比较，精度统计结果列于表3中。

图4 传统最小二乘法延拓模型的延拓结果等值线图(单位：mGal)Fig.4 Isoline chart of the continuation results of traditional least squares model(Unit：mGal)

其次，利用改进的最小二乘法延拓模型，将航空重力测量数据向下延拓到地面平均高程面上，其等值线图如图5所示。将延拓结果与地面重力测量数据进行比较，精度统计结果列于表3中。

再次，在使用Tikhonov正则化法延拓模型之前，需要确定其正则化参数。根据地面已知重力测量数据，计算了延拓误差与正则化参数的关系，如图6所示，并确定了最优正则化参数α=0.031。采用该正则化参数，将航空重力测量数据向下延拓，延拓结果等值线图如图7所示，精度统计结果列入表3中。

图5 改进最小二乘法延拓模型的延拓结果等值线图(单位：mGal)Fig.5 Isoline chart of the continuation results of improved least squares model(Unit：mGal)

表3 延拓结果统计

图6 延拓误差与正则化参数关系图Fig.6 Relationship of continuation errorsand regularization factors

最后，利用本文建立的Poisson积分迭代法延拓模型将航空重力测量数据向下延拓，延拓结果的等值线图如图8、图9所示，精度统计结果列于表3中。

图7 Tikhonov正则化法延拓模型的延拓结果等值线图(单位：mGal)Fig.7 Isoline chart of the continuation results of Tikhonov regularization model(Unit：mGal)

图8 Poisson积分迭代法延拓模型I的延拓结果等值线图(单位：mGal)Fig.8 Isoline chart of the continuation results of Poisson integral iteration model I(Unit：mGal)

图9 Poisson积分迭代法延拓模型Ⅱ的延拓结果等值线图(单位：mGal)Fig.9 Isoline chart of the continuation results of Poisson integral iteration model Ⅱ(Unit：mGal)

根据表3的统计结果，相比较于传统的最小二乘法延拓模型，改进的最小二乘法延拓模型的结果精度提高了约15 mGal，采用Tikhonov正则化法对延拓模型的病态性进行修正以后，延拓结果精度进一步提高了约0.34 mGal，本文新建立的Poisson积分迭代法延拓模型精度要优于改进的最小二乘法延拓模型，精度提高了约0.21 mGal，但稍逊于Tikhonov正则化法延拓模型，精度略低于0.13 mGal。

通过Poisson积分迭代法延拓模型Ⅰ和Ⅱ的结果可以看出，不论采用传统还是改进的Poisson积分离散化公式，延拓结果几乎一致，这说明本文建立的Poisson积分迭代法延拓模型不仅可以修正传统Poisson积分离散化误差给延拓模型带来的病态性影响，而且对延拓模型本身的病态性有一定的抑制作用(精度优于改进的最小二乘法延拓模型)。

6 结束语

航空重力测量技术，因其可以快速经济有效地获取分布均匀、精度良好、大面积的重力场中高频信息，在我国陆地和陆海交界区域的重力测量中得到了很好的应用。在国内外学者研究成果的基础上，本文利用实测航空和地面重力数据，对传统最小二乘法、改进最小二乘法、Tikhonov正则化法等延拓模型进行了数值分析，最后建立了Poisson积分迭代法和改进Poisson积分迭代法延拓模型。本文的研究成果，可应用于我国航空重力标量和矢量测量数据的处理中，从而为获得更高精度、更高分辨率的地球重力场产品提供一定的技术支撑。