2. 西安测绘研究所, 陕西 西安 710054;
3. 大地测量与地球动力学国家重点实验室, 湖北 武汉 430077;
4. 信息工程大学地理空间信息学院, 河南 郑州 450052
2. Xi'an Research Institute of Surveying and Mapping, Xi'an 710054, China;
3. State Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China;
4. Institute of Geospatial Information, Information Engineering University, Zhengzhou 450052, China
航空重力测量技术,因其可以在沙漠、沼泽、冰川、原始森林、陆海交界等一些难以开展地面重力测量的区域进行作业,快速经济地获取精度良好、分布均匀、大面积的地球重力场中高频信息,从而成为地球重力场研究的最为热门的领域之一[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]
式中,Δgh(r, θ, λ)表示空中某一点的航空重力数据;Δg0(R, θ′, λ′)表示地面上某一点的重力数据;
由于在实际情况下空中和地面重力异常均是离散值,式(1)可写成如下形式的观测方程
式中,gh为空中重力异常Δgh所组成的向量;X为地面区域的重力异常Δg0所组成的向量;A为相应的系数阵。在航空重力测量中,空中重力异常Δgh是已知观测值,地面重力异常为待估参数,这是Poisson积分方程的逆问题。
式(2)的传统离散化形式为
式中,矩阵Aij的表达式为
式中,
对系数矩阵A作如下奇异值分解
式中,U、V分别由矩阵A的左、右特征向量组成;
考虑到式(5),式(2)的最小二乘解及其谱分解形式为[19]
式中,λi为A的奇异值;ui、vi分别为其相应的左右奇异向量;rank(A)表示矩阵A的秩;而l=gh为误差方程的自由项。
式(6)即为传统的最小二乘法延拓模型。
2 改进最小二乘法延拓模型采用式(4)所示的Poisson积分离散化形式在某些情况下不能得到正确的结果,甚至误差很大。主要原因是,与实际结果相比,处于计算点向径方向的流动点对计算点的贡献较大,如果以该点的核函数值作为整个格网的核函数值,将会带来较大误差。因此,有学者给出了修改向径方向流动点加权系数的离散化形式。矩阵A的对角线元素为[27-30]
矩阵A的非对角线元素为
式中,ψ0表示选取积分半径;NC表示积分半径内的测量点个数。
该离散化公式的推导思路是先将积分半径内的其他格网对计算点的延拓贡献值求出,再将整个积分范围内的理论贡献值与之相减,即得到向径方向的流动点对计算点的贡献。
将式(6)中的矩阵A用式(7)和式(8)来替换,就可以得到改进的最小二乘法延拓模型。
3 Tikhonov正则化法延拓模型航空重力测量数据向下延拓时,观测方程是病态的,其系数矩阵A的奇异值单调地趋向零,由式(6)难以获得稳定的解。因此,需对病态方程进行正则化处理,以抑制高频测量噪声对延拓结果的影响。在众多正则化方法中,以Tikhonov正则化法的应用最为广泛,所采用的准则如下
式中,α>0为正则化参数;‖X‖表示X的范数。
根据式(9)的约束条件,可得Tikhonov正则化解为
此处,Xα的下标表示相应于正则化参数α的解。σ02为单位权方差。
Tikhonov正则化解是有偏的,其偏差为
其中,X为X的真值。
Tikhonov正则化解的平均均方误差为
将式(5)代入式(10)、(12),可得Tikhonov正则化解的谱分解式为[19]
通过实际计算表明,采用式(13)来获得Tikhonov正则化解的速度要优于式(10),这是因为谱分解式避免了许多大型矩阵的乘法和求逆运算,节约了大量计算时间。
正则化解Xα与其真值X之差的范数的数学期望为
上式表明,Tikhonov正则化解的误差由两部分组成:第一项为观测误差所引起的估值误差,随α的增大单调减小;第二项为正则化所引起的估值误差,随α的增大单调增大,当α=0时,该项误差为0。选择合适的正则化参数α可起到平衡这两项误差的作用,使得它们的和为最小,这正是通过正则化算法可获得精确、稳定解的原因。
4 改进Poisson积分迭代法延拓模型如何正确地确定正则化参数,则是整个正则化算法的关键,参数选择的优劣直接影响了最终延拓结果的精度。在模拟试验中,最佳正则化参数的确定,是在延拓面数据已知的情况下根据最小延拓误差来估算的。而对于实际的航空重力测量数据向下延拓,延拓面并不一定存在已知的测量数据。因此,在这种情况下如何确定正则化参数,存在较大的难度。为了减弱延拓模型的病态性,提高延拓结果的精度,本文提出Poisson积分迭代法这一延拓思路,进而建立相应的延拓模型。
由于R/r也是调和函数,因此有
在式(5)两端同时乘以Δg0(R, θ, λ),有
将式(1)减去式(16),通过变换可以得到
对式(17)进行改化,可以得到
由于式(18)两端都含有Δg0(R, θ, λ)项,因此,该式可以进一步改化成迭代形式。作为第一次近似,命
再用下式进行迭代计算
式(20)则为本文建立的Poisson积分迭代法向下延拓模型,该式可以进一步修改成矩阵计算的形式
式中,矩阵A可以用式(4)来计算(本文称之为新建模型Ⅰ,即Poisson积分迭代法向下延拓模型),也可以用式(7)和式(8)来计算(新建模型Ⅱ,即改进Poisson积分迭代法向下延拓模型)。当前后两次计算值之间的差异小于ε时,就可以终止迭代。
5 数值试验与结果分析 5.1 试验数据说明采用中国某地区实测航空重力测量数据,分辨率是5′×5′,飞行高度为3400 m。其等值线图如图 1所示。
相应的地面重力测量数据,其等值线图如图 2所示。
航空和地面重力测量数据的统计结果如表 1所示。
mGal | ||||
数据类型 | 最大值 | 最小值 | 平均值 | 标准差 |
航空重力测量数据 | 86.950 | -42.800 | 16.214 | 27.692 |
地面重力测量数据 | 101.210 | -45.140 | 18.056 | 29.537 |
根据数字高程模型SRTM30,采用Global Mapper软件获得了该地区5′×5′的地形数据,其等值线图如图 3所示。
地面重力测量数据对应的地形数据统计结果如表 2所示。
以地面点的平均地形高度1 476.115 m作为向下延拓的基准面,则延拓高度是1 923.885 m。
根据图 1、图 2与图 3的对比可以看出,航空和地面重力数据的变化趋势,与地形的起伏密切相关。地形起伏较为剧烈的地区,重力数据也有相应的变化。
5.2 对比试验结果首先,利用传统的最小二乘法延拓模型,将航空重力测量数据向下延拓到地面平均高程面上,其等值线图如图 4所示。将延拓结果与地面重力测量数据进行比较,精度统计结果列于表 3中。
mGal | ||||
模型类别 | 最大值 | 最小值 | 平均值 | 标准差 |
传统最小二乘法延拓模型 | 73.69 | -32.01 | 12.92 | 20.65 |
改进最小二乘法延拓模型 | 14.45 | -20.05 | 1.49 | 5.60 |
Tikhonov正则化法延拓模型 | 15.41 | -17.25 | 2.01 | 5.26 |
Poisson积分迭代法延拓模型Ⅰ | 15.10 | -18.50 | 1.82 | 5.39 |
Poisson积分迭代法延拓模型Ⅱ | 15.14 | -18.42 | 1.84 | 5.39 |
其次,利用改进的最小二乘法延拓模型,将航空重力测量数据向下延拓到地面平均高程面上,其等值线图如图 5所示。将延拓结果与地面重力测量数据进行比较,精度统计结果列于表 3中。
再次,在使用Tikhonov正则化法延拓模型之前,需要确定其正则化参数。根据地面已知重力测量数据,计算了延拓误差与正则化参数的关系,如图 6所示,并确定了最优正则化参数α=0.031。采用该正则化参数,将航空重力测量数据向下延拓,延拓结果等值线图如图 7所示,精度统计结果列入表 3中。
最后,利用本文建立的Poisson积分迭代法延拓模型将航空重力测量数据向下延拓,延拓结果的等值线图如图 8、图 9所示,精度统计结果列于表 3中。
根据表 3的统计结果,相比较于传统的最小二乘法延拓模型,改进的最小二乘法延拓模型的结果精度提高了约15 mGal,采用Tikhonov正则化法对延拓模型的病态性进行修正以后,延拓结果精度进一步提高了约0.34 mGal,本文新建立的Poisson积分迭代法延拓模型精度要优于改进的最小二乘法延拓模型,精度提高了约0.21 mGal,但稍逊于Tikhonov正则化法延拓模型,精度略低于0.13 mGal。
通过Poisson积分迭代法延拓模型Ⅰ和Ⅱ的结果可以看出,不论采用传统还是改进的Poisson积分离散化公式,延拓结果几乎一致,这说明本文建立的Poisson积分迭代法延拓模型不仅可以修正传统Poisson积分离散化误差给延拓模型带来的病态性影响,而且对延拓模型本身的病态性有一定的抑制作用(精度优于改进的最小二乘法延拓模型)。
6 结束语航空重力测量技术,因其可以快速经济有效地获取分布均匀、精度良好、大面积的重力场中高频信息,在我国陆地和陆海交界区域的重力测量中得到了很好的应用。在国内外学者研究成果的基础上,本文利用实测航空和地面重力数据,对传统最小二乘法、改进最小二乘法、Tikhonov正则化法等延拓模型进行了数值分析,最后建立了Poisson积分迭代法和改进Poisson积分迭代法延拓模型。本文的研究成果,可应用于我国航空重力标量和矢量测量数据的处理中,从而为获得更高精度、更高分辨率的地球重力场产品提供一定的技术支撑。
[1] |
孙中苗. 航空重力测量理论、方法及应用研究[D]. 郑州: 信息工程大学, 2004. SUN Zhongmiao. Theory, Methods and Applications of Airborne Gravimetry[D]. Zhengzhou: Information Engineering University, 2004. |
[2] |
孙中苗, 夏哲仁, 石磐.
航空重力测量研究进展[J]. 地球物理学进展, 2004, 19(3): 492–496.
SUN Zhongmiao, XIA Zheren, SHI Pan. Progress and Advance in Airborne Gravimetry[J]. Progress in Geophysics, 2004, 19(3): 492–496. DOI:10.3969/j.issn.1004-2903.2004.03.002 |
[3] |
赵德军. 航空矢量重力测量的理论与方法[D]. 郑州: 信息工程大学, 2005. ZHAO Dejun. Theory and Methods of Airborne Vector Gravimetry System[D]. Zhengzhou: Information Engineering University, 2005. |
[4] | GOLI M, NAJAFI-ALAMDARI M, VANÍČEK P. Numerical Behaviour of the Downward Continuation of Gravity Anomalies[J]. Studia Geophysica et Geodaetica, 2011, 55(2): 191–202. DOI:10.1007/s11200-011-0011-8 |
[5] |
刘晓刚, 李姗姗, 吴星.
卫星重力梯度数据的向下延拓[J]. 大地测量与地球动力学, 2011, 31(1): 132–137.
LIU Xiaogang, LI Shanshan, WU Xing. Downward Continuation of Satellite Gravity Gradient Data[J]. Journal of Geodesy and Geodynamics, 2011, 31(1): 132–137. |
[6] |
刘晓刚, 明锋, 常宜峰, 等.
航空重力数据向下延拓的频域正则化迭代法[J]. 测绘科学与工程, 2014, 34(2): 9–14, 19.
LIU Xiaogang, MING Feng, CHANG Yifeng, et al. Downward Continuation of Airborne Gravimetry Data by Regularization Iterative Methods in Frequency Domain[J]. Geomatics Science and Engineering, 2014, 34(2): 9–14, 19. |
[7] | HEISKANEN W A, MORITZ H. Physical Geodesy[M]. San Francisco: W. H. Freeman and Company, 1967. |
[8] |
张嗥, 陈琼, 丛明日.
航空重力测量数据向下延拓中空间协方差函数特性研究[J]. 测绘科学, 2006, 31(4): 51–53.
ZHANG Hao, CHEN Qiong, CONG Mingri. Research on Space Covariance Function's Characteristic During Downward Continuation of Airborne Gravity Measurement Data[J]. Science of Surveying and Mapping, 2006, 31(4): 51–53. DOI:10.3771/j.issn.1009-2307.2006.04.015 |
[9] |
翟振和, 孙中苗.
基于配置法的局部重力场延拓模型构建与应用分析[J]. 地球物理学报, 2009, 52(7): 1700–1706.
ZHAI Zhenhe, SUN Zhongmiao. Continuation Model Construction and Application Analysis of Local Gravity Field Based on Least Square Collocation[J]. Chinese Journal of Geophysics, 2009, 52(7): 1700–1706. DOI:10.3969/j.issn.0001-5733.2009.07.004 |
[10] |
张松堂, 陆银龙.
基于少量地面控制的航空重力测量[J]. 海洋测绘, 2009, 29(3): 7–11.
ZHANG Songtang, LU Yinlong. Airborne Gravimetry by A Few Ground Controls[J]. Hydrographic Surveying and Charting, 2009, 29(3): 7–11. DOI:10.3969/j.issn.1671-3044.2009.03.003 |
[11] | FORSBERG R. Downward Continuation of Airborne Gravity Data[EB/OL]. [2016-10-10]. http://74.125.155.132/scholar?q=cache:K4rJHdrG6LMJ,2010. |
[12] |
王兴涛, 夏哲仁, 石磐, 等.
航空重力测量数据向下延拓方法比较[J]. 地球物理学报, 2004, 47(6): 1017–1022.
WANG Xingtao, XIA Zheren, SHI Pan, et al. A Comparison of Different Downward Continuation Methods for Airborne Gravity Data[J]. Chinese Journal of Geophysics, 2004, 47(6): 1017–1022. DOI:10.3321/j.issn:0001-5733.2004.06.013 |
[13] |
周波阳, 罗志才, 许闯, 等.
航空重力数据向下延拓的FFT快速算法比较[J]. 大地测量与地球动力学, 2013, 33(1): 64–68, 73.
ZHOU Boyang, LUO Zhicai, XU Chuang, et al. Comparison Among Fast Fourier Transform Algorithms for Downward Continuation of Airborne Gravity Data[J]. Journal of Geodesy and Geodynamics, 2013, 33(1): 64–68, 73. |
[14] |
黄谟涛, 欧阳永忠, 刘敏, 等.
海域航空重力测量数据向下延拓的实用方法[J]. 武汉大学学报(信息科学版), 2014, 39(10): 1147–1152.
HUANG Motao, OUYANG Yongzhong, LIU Min, et al. Practical Methods for the Downward Continuation of Airborne Gravity Data in the Sea Area[J]. Geomatics and Information Science of Wuhan University, 2014, 39(10): 1147–1152. |
[15] |
黄谟涛, 宁津生, 欧阳永忠, 等.
联合使用位模型和地形信息的陆区航空重力向下延拓方法[J]. 测绘学报, 2015, 44(4): 355–362.
HUANG Motao, NING Jinsheng, OUYANG Yongzhong, et al. Downward Continuation of Airborne Gravimetry on Land Using Geopotential Model and Terrain Information[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(4): 355–362. DOI:10.11947/j.AGCS.2015.20130751 |
[16] |
刘敏, 黄谟涛, 欧阳永忠, 等.
顾及地形效应的重力向下延拓模型分析与检验[J]. 测绘学报, 2016, 45(5): 521–530.
LIU Min, HUANG Motao, OUYANG Yongzhong, et al. Test and Analysis of Downward Continuation Models for Airborne Gravity Data with regard to the Effect of Topographic Height[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(5): 521–530. DOI:10.11947/j.AGCS.2016.20150453 |
[17] | ILK K H, KUSCHE J, RUDOLPH S. A Contribution to Data Combination in Ill-posed Downward Continuation Problems[J]. Journal of Geodynamics, 2002, 33(1-2): 75–99. DOI:10.1016/S0264-3707(01)00056-4 |
[18] | KERN M. An Analysis of the Combination and Downward Continuation of Satellite, Airborne and Terrestrial Gravity Data[R]. Calgary: University of Calgary, 2003. |
[19] |
王兴涛, 石磐, 朱非洲.
航空重力测量数据向下延拓的正则化算法及其谱分解[J]. 测绘学报, 2004, 33(1): 33–38.
WANG Xingtao, SHI Pan, ZHU Feizhou. Regularization Methods and Spectral Decomposition for the Downward Continuation of Airborne Gravity Data[J]. Acta Geodaetica et Cartographica Sinica, 2004, 33(1): 33–38. DOI:10.3321/j.issn:1001-1595.2004.01.006 |
[20] |
成怡, 郝燕玲, 刘繁明.
航空重力测量数据向下延拓及其影响因素分析[J]. 系统仿真学报, 2008, 20(8): 2190–2194.
CHEN Yi, HAO Yanling, LIU Fanming. Analysis of Effect and Downward Continuation of Airborne Gravity Data[J]. Journal of System Simulation, 2008, 20(8): 2190–2194. |
[21] |
郝燕玲, 成怡, 孙枫, 等.
Tikhonov正则化向下延拓算法仿真实验研究[J]. 仪器仪表学报, 2008, 29(3): 605–609.
HAO Yanling, CHEN Yi, SUN Feng, et al. Simulation Research on Tikhonov Regularation Algorithm in Downward Continuation[J]. Chinese Journal of Scientific Instrument, 2008, 29(3): 605–609. DOI:10.3321/j.issn:0254-3087.2008.03.030 |
[22] |
邓凯亮, 暴景阳, 黄谟涛, 等.
航空重力数据向下延拓的Tikhonov正则化法仿真研究[J]. 武汉大学学报(信息科学版), 2010, 35(12): 1414–1417.
DENG Kailiang, BAO Jingyang, HUANG Motao, et al. Simulation of Tikhonov Regulation Algorithm in Downward Continuation of Airborne Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2010, 35(12): 1414–1417. |
[23] |
邓凯亮, 黄谟涛, 暴景阳, 等.
向下延拓航空重力数据的Tikhonov双参数正则化法[J]. 测绘学报, 2011, 40(6): 690–696.
DENG Kailiang, HUANG Motao, BAO Jingyang, et al. Tikhonov Two-parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 690–696. |
[24] |
蒋涛, 李建成, 王正涛, 等.
航空重力向下延拓病态问题的求解[J]. 测绘学报, 2011, 40(6): 684–689.
JIANG Tao, LI Jiancheng, WANG Zhengtao, et al. Solution of Ill-posed Problem in Downward Continuation of Airborne Gravity[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 684–689. |
[25] |
孙文, 吴晓平, 王庆宾, 等.
航空重力数据向下延拓的波数域迭代Tikhonov正则化方法[J]. 测绘学报, 2014, 43(6): 566–574.
SUN Wen, WU Xiaoping, WANG Qingbin, et al. Wave Number Domain Iterative Tikhonov Regularization Method for Downward Continuation of Airborne Gravity Data[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(6): 566–574. DOI:10.13485/j.cnki.11-2089.2014.0080 |
[26] |
刘晓刚, 李迎春, 肖云, 等.
重力与磁力测量数据向下延拓中最优正则化参数确定方法[J]. 测绘学报, 2014, 43(9): 881–887.
LIU Xiaogang, LI Yingchun, XIAO Yun, et al. Optimal Regularization Parameter Determination Method in Downward Continuation of Gravimetric and Geomagnetic Data[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(9): 881–887. DOI:10.13485/j.cnki.11-2089.2014.0160 |
[27] | VANÍČEK P, SUN W, ONG P, et al. Downward Continuation of Helmert's Gravity[J]. Journal of Geodesy, 1996, 71(1): 21–34. DOI:10.1007/s001900050072 |
[28] | SUN W, VANÍČEK P. On Some Problems of the Downward Continuation of the 5'×5' Mean Helmert Gravity Disturbance[J]. Journal of Geodesy, 1998, 72(7-8): 411–420. DOI:10.1007/s001900050180 |
[29] |
管斌, 孙中苗, 刘晓刚.
重力异常向上延拓中Poisson积分离散化方法比较[J]. 测绘科学与工程, 2014, 34(6): 12–17.
GUAN Bin, SUN Zhongmiao, LIU Xiaogang. Comparison of Discretization Methods of Poisson Integral in Upward Continuation of Gravity Anomaly[J]. Geomatics Science and Engineering, 2014, 34(6): 12–17. |
[30] |
刘晓刚, 孙中苗, 管斌, 等.
Poisson积分离散化改进公式在航空重力向下延拓中的应用[J]. 测绘科学与工程, 2017, 37(1): 5–9.
LIU Xiaogang, SUN Zhongmiao, GUAN Bin, et al. Application of Improved Poisson Integral Discretization Formula in Downward Continuation of Airborne Gravimetry Data[J]. Geomatics Science and Engineering, 2017, 37(1): 5–9. |