Analysis of Boundary Effect Suppression in Downward Continuation of Airborne Gravimetric Data by Filling the Data Gap

 

LIU Xiao-gang1,2,3，XIAO Yun1,2，LI Ying-chun1,2， ZHAI Zhen-he1,2

(1. Xi’an Research Institute of Surveying and Mapping, Xi’an710054, China; 2. National 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)

 

Abstract  Downward continuation is one of the key steps in the processing of airborne gravimetric data, however, data at the boundary area will be distorted along with the changing trend of nearby data during downward continuation of airborne gravity data, and this will result in the increase of continuation error of surveying data at the boundary area. Therefore, the influence of boundary effect to the continuation results should be analyzed, and the data in the gap of boundary area should be simulated and filled, then the boundary area can be constrained and the boundary effect to the continuation results can also be suppressed. The results show that before and after the data gap is filled, the precision loss to the surveying error of iterative Tikhonov regularization method, Landweber regularization iteration method and gradient method are 30%, 23.5%, 40.4% and 28%, 21.3%, 38.7%, respectively. The precision are improved by 2, 2.2 and 1.7 percent, so the precision of downward continuation results are improved at a certain extent.

Key words  Airborne Gravimetry；Downward Continuation；Boundary Effect；Earth’s Gravitational Field Model(EGM)；Regularization

 

<!--[if !supportLists]-->1 <!--[endif]-->引言

在实际的工作中，航空重力测量常常是在起伏的航线上进行的，为了实现不同类型重力测量数据的融合以及为地球重力场模型的构建提供基础数据支撑，需要将航空重力测量数据向下延拓到地球表面[1-12]。另外，利用向下延拓可以把重力测量数据的观测面外推到场源体附近，突出局部场，有利于提高重力测量数据解释的可靠性。

由于航空重力测量数据的不均匀或地理因素的限制（边界区域处存在危险地带或国家边界的限制等区域），边界区域处航空重力测量数据点就比较稀少，缺少测量点的边界区域空白处会失去实际测量点的约束，受到其邻近区域函数变化趋势的影响，从而呈现出自由的变化状态，因此，在航空重力数据向下延拓时，缺少测量数据的边界区域处会随着其邻近区域数据的变化趋势而扭曲变形，向下延拓的高度越高，这种扭曲变形的程度越为严重，相应边界处延拓结果的误差也会越大，这就是通常所说的边界效应。

因此，在航空重力测量数据向下延拓时，如何有效抑制边界效应的影响，是提高航空重力测量数据向下延拓结果精度的一个重要措施。目前，在航空重力数据向下延拓的边界效应抑制理论与方法的研究上，有部分学者已有一定的研究成果。王兴涛在将航空重力测量数据向下延拓时，顾及了外区地面重力异常对空中重力异常的影响，从而有效地减小了边界效应的影响，使向下延拓数据的可用范围基本与空中相同，但是，随着延拓距离的增加，边界效应也明显增加[13]。蒋涛提出了基于矩谐分析的航空重力向下延拓方法，以重力扰动作为基本观测值，给出了基于矩谐分析的向下延拓模型和算法，比较结果表明，矩谐分析在延拓精度、稳定性和边界效应等方面都要优于直接法和基于广义岭估计的逆泊松积分法，能取得良好的向下延拓效果[14]。此外，很少见到关于边界效应问题处理方法的报道。

基于上述考虑，本文研究航空重力测量数据向下延拓时边界效应的抑制理论与方法，对缺少重力测量数据的边界区域进行模拟仿真，填充了数据空白，对缺少测量点的边界区域进行约束，以便抑制边界效应对向下延拓结果的影响。

<!--[if !supportLists]-->2 <!--[endif]-->向下延拓的基础模型

航空重力测量数据向下延拓的基本原理如图1所示：

图1  向下延拓示意图

表示观测面上的航空重力测量数据，表示延拓面上的航空重力测量数据，由观测面上的重力数据计算延拓面上的重力数据，这种换算称为航空重力测量数据的向下延拓。

本文在研究航空重力测量数据向下延拓的边界效应抑制方法时，以迭代Tikhonov正则化法、Landweber正则化迭代法和梯度法等三种数学模型作为基础模型。

迭代Tikhonov正则化法的公式为[15]

Landweber正则化迭代法的公式为[15]

其中，α为正则化参数，用于平衡不稳定性及光滑性；n为迭代次数；为的傅里叶变换；，u、v分别表示空域变量x、y对应的频域变量。

梯度法截取到二阶项的公式为[16,17]

其中，表示延拓距离，r表示观测面上某一点的地心向径，表示延拓面上对应点的地心向径，表示观测面上某一点的地心余纬和地心经度。

<!--[if !supportLists]-->3 <!--[endif]-->边界效应抑制方法

为了抑制边界效应对向下延拓精度的影响，本文拟采用三种思路：一是采用已有的高精度重力场模型，对边界区域的数据进行模拟后填充，或者采用地面上已有的重力测量数据，向上延拓到测量高度后对数据空白处进行填充；二是采用一些插值方法（如拉格朗日插值、牛顿插值、分段线性插值和三次样条插值等），通过向外插值的方式扩大测量范围，使边界点变成内部点，从而对缺少测量点的边界区域进行插值约束；三是采用半余弦函数扩边法，将航空重力测量数据在平面四个方向上进行扩边处理，从而增加测量数据范围，减弱边界效应的影响。基本思路如图2所示：

图2  解决边界效应的基本思路图

由于时间关系，本文在具体数值试验时，目前仅仅测试了数据空白填充这种方法。

<!--[if !supportLists]-->4 <!--[endif]-->数值实验与结果分析

<!--[if !supportLists]-->4.1 <!--[endif]-->航空重力测量数据仿真

采用2160阶的EGM2008地球重力场模型，对中国某地区的航空重力测量数据和地面重力数据进行仿真。

将航空重力测量飞机的飞行高度设为3000m，分辨率为5′×5′。图3和图4分别表示观测面理论航空重力测量数据和延拓面理论重力数据。

 

图3  观测面理论航空重力测量数据等值线图（mGal）

 

图4  延拓面理论重力数据等值线图（mGal）

为了验证本文向下延拓算法的有效性，在观测面理论航空重力测量数据中加入均值为0、方差为9mGal的高斯白噪声，其等值线图如图5所示。表1为几种数据的统计结果。

 

图5  含有高斯白噪声的观测面航空重力测量数据等值线图（mGal）

表1  数据统计结果（mGal）

统计值 观测面高度	最大值	最小值	平均值	标准差
3000m（不含噪）	61.597	-51.053	-2.399	21.135
3000m（含噪）	68.799	-53.366	-2.383	21.306
0m	74.977	-61.050	-2.242	23.364

 

<!--[if !supportLists]-->4.2 <!--[endif]-->数据空白填充前实验结果

（1）迭代Tikhonov正则化法

在对迭代Tikhonov正则化方法的延拓精度进行测试时，首先需要对延拓误差与正则化因子的关系进行验证。当迭代次数分别取10次、20次、30次和40次时，延拓误差与正则化因子的关系如图6所示。

 

图6  延拓误差与正则化因子关系图（mGal）

从图6可以看出，对于不同的迭代次数，迭代Tikhonov正则化法在曲线的轨点处取得了最小的延拓误差，而且迭代次数越多时，所对应的最佳正则化因子也越大。

表2表示迭代Tikhonov正则化法在取得最小延拓误差时对应的正则化因子和迭代次数。因此，表2的统计结果进一步证实了图6所得出的结论。

表2  迭代Tikhonov正则化法最小延拓误差对应的正则化因子和迭代次数（mGal）

正则化因子（α）	迭代次数（n）	最大值	最小值	平均值	均方差
2.581	10	15.113	-12.549	0.051	4.002
5.411	20	14.711	-12.297	0.061	3.927
8.251	30	14.579	-12.208	0.064	3.902

 

图7给出了迭代Tikhonov正则化法在取得最小延拓误差（即α=8.251，n=30）时对应的延拓结果的等值线图。

 

图7  迭代Tikhonov正则化法延拓结果等值线图（mGal）

从图4和图7的比较来看，迭代Tikhonov正则化方法延拓结果对延拓面上理论重力数据的逼近效果很好。

（2）Landweber正则化迭代法

在对Landweber正则化迭代法的延拓精度进行测试时，首先需要对延拓误差与正则化因子及迭代次数的关系进行验证。具体计算时，取正则化因子分别为α=0.5、1.0、1.5、2.0，延拓误差与迭代次数的关系如图8所示：

 

图8  延拓误差与迭代次数关系图（mGal）

从图8可以看出，对于不同的正则化因子，Landweber正则化迭代法在曲线的轨点处取得了最小的延拓误差。但是当正则化因子α=2.0时，延拓误差突然增大，因此可以得知，选取的正则化因子不宜过大，以小于2.0为宜。

当迭代次数分别取10次、20次、30次和40次时，延拓误差与正则化因子的关系如图9所示。

 

图9  延拓误差与正则化因子关系图（mGal）

从图9可以看出，当迭代次数越多时，在取得最小延拓误差处对应的正则化因子越小，与迭代Tikhonov正则化法取得的结论正好相反，这是由于Landweber正则化迭代法的表达式和迭代Tikhonov正则化法的表达式存在的显著不同引起的。另外，在正则化因子大于2.0时，可以看到延拓误差急剧增大，这与图8显示的结果相符，从而说明了在采用Landweber正则化迭代法进行航空重力测量数据向下延拓时，正则化因子应该小于2.0。

表3表示Landweber正则化迭代法在延拓结果取得最小延拓误差时对应的正则化因子和迭代次数。

表3  Landweber正则化迭代法最小延拓误差对应的正则化因子和迭代次数（mGal）

正则化因子（α）	迭代次数（n）	最大值	最小值	平均值	均方差
0.321	10	13.465	-11.591	0.091	3.705
0.171	20	13.809	-11.882	0.085	3.779
0.111	30	14.159	-11.727	0.071	3.806
0.091	40	13.864	-12.259	0.089	3.827

 

下图10给出了Landweber正则化迭代法在取得最小延拓误差（即α=0.321）时对应的延拓结果的等值线图。

 

图10  Landweber正则化迭代法延拓结果等值线图（mGal）

从图4和图10的比较来看，Landweber正则化迭代法延拓结果对延拓面上的理论重力数据的逼近效果很好。

（3）梯度法

图11给出了利用梯度法将观测面上含有噪声的航空重力测量数据向下延拓3000m后延拓结果的等值线图。

 

图11  梯度法延拓结果等值线图（mGal）

从图4和图11的比较来看，梯度法延拓结果对延拓面上的理论磁异常数据的逼近效果比较好。

梯度法延拓结果的精度统计情况见表4。

表4  梯度法延拓结果精度统计（mGal）

最大值	最小值	平均值	均方差
16.459	-15.804	0.141	4.211

 

<!--[if !supportLists]-->4.3 <!--[endif]-->数据空白填充后实验结果

采用2160阶EGM2008地球重力场模型，对试验区周边1°范围内的航空重力测量数据进行仿真，从而对迭代Tikhonov正则化法、Landweber正则化迭代法和梯度法等三种航空重力测量数据向下延拓基础模型的边界效应抑制效果进行分析。

（1）迭代Tikhonov正则化法

根据表2的统计结果，将数据空白填充后迭代Tikhonov正则化法取得的延拓结果列于表5中。

表5  数据空白填充后迭代Tikhonov正则化法延拓结果统计结果（mGal）

正则化因子（α）	迭代次数（n）	最大值	最小值	平均值	均方差
2.581	10	15.113	-12.549	-0.021	3.945
5.411	20	14.324	-13.810	-0.011	3.868
8.251	30	14.236	-13.716	-0.008	3.842

 

与表2的统计结果进行比较后可以看出，在对数据空白区进行填充后，得到的延拓结果精度普遍比未填充结果的精度提高约0.06mGal的水平，也就是说，在数据空白未填充之前，迭代Tikhonov正则化法对测量数据进行延拓造成的精度损失为测量误差的30%，而在数据空白填充之后，造成的精度损失为测量误差的28%，提高了大约2个百分点，有一定的提升效果。

α=8.251、n=30时对应的延拓结果的等值线图如图12所示。

 

图12  数据空白填充后迭代Tikhonov正则化法延拓结果等值线图（mGal）

（2）Landweber正则化迭代法

根据表3的统计结果，将数据空白填充后Landweber正则化迭代法取得的延拓结果列于表6中。

表6  数据空白填充后Landweber正则化迭代法延拓结果统计结果（mGal）

正则化因子（α）	迭代次数（n）	最大值	最小值	平均值	均方差
0.321	10	13.484	-12.913	0.018	3.638
0.171	20	13.674	-13.167	0.011	3.714
0.111	30	14.013	-13.408	0.013	3.744
0.091	40	13.622	-13.219	0.013	3.761

 

与表3的统计结果进行比较后可以看出，在对数据空白区进行填充后，得到的延拓结果精度普遍比未填充结果的精度提高约0.07mGal的水平，也就是说，在数据空白未填充之前，Landweber正则化迭代法对测量数据进行延拓造成的精度损失为测量误差的23.5%，而在数据空白填充之后，造成的精度损失为测量误差的21.3%，提高了大约2.2个百分点，有一定的提升效果。

α=0.321、n=10时对应的延拓结果的等值线图如图13所示。

 

图13  数据空白填充后Landweber正则化迭代法延拓结果等值线图（mGal）

（3）梯度法

将数据空白填充后的航空重力测量数据向下延拓，其延拓结果列于表7中。

表7  数据空白填充后的梯度法延拓结果精度统计（mGal）

最大值	最小值	平均值	均方差
16.290	-15.847	0.136	4.162

与表4的统计结果进行比较后可以看出，在对数据空白区进行填充后，得到的延拓结果精度普遍比未填充结果的精度提高约0.05mGal的水平，也就是说，在数据空白未填充之前，梯度法对测量数据进行延拓造成的精度损失为测量误差的40.4%，而在数据空白填充之后，造成的精度损失为测量误差的38.7%，提高了大约1.7个百分点，有一定的提升效果。

延拓结果对应的等值线图如图14所示。

 

图14  数据空白填充后梯度法延拓结果等值线图（mGal）

<!--[if !supportLists]-->5 <!--[endif]-->结束语

根据迭代Tikhonov正则化法、Landweber正则化迭代法、梯度法等三种数学模型在数据空白填充后的向下延拓结果来看，采用数据空白填充这种方法来抑制边界效应，虽然在一定程度上提高了向下延拓结果的精度，但提升效果并不十分理想。因此，下一步还需要测试其它两种边界效应抑制方法，从而减弱边界效应的影响，进一步提高航空重力测量数据向下延拓的精度。

 

参考文献

<!--[if !supportLists]-->[1] <!--[endif]-->Novák P, Kern M, Schwarz K P. Numerical studies on the harmonic downward continuation of band-limited airborne gravity[J]. Studia geoph. et geod., 2001, 45: 317-345.

<!--[if !supportLists]-->[2] <!--[endif]-->Novák P, Heck B. Downward continuation and geoid determination based on band-limited airborne gravity data[J]. Journal of Geodesy, 2002, 76: 269-278.

<!--[if !supportLists]-->[3] <!--[endif]-->Wang X T, Xia Z R, Shi P, et al. A comparison of different downward continuation methods for airborne gravity data[J]. Chinese J. Geophys., 2004, 47(6): 1017-1022. (王兴涛，夏哲仁，石磐，等．航空重力测量数据向下延拓方法比较[J]．地球物理学报，2004，47(6)：1017-1022.)

<!--[if !supportLists]-->[4] <!--[endif]-->Sun Z M. Theory, methods and applications of airborne gravimetry[D]. Zhengzhou: Information Engineering University, 2004. (孙中苗. 航空重力测量理论方法及应用研究[D]. 郑州: 信息工程大学, 2004.)

<!--[if !supportLists]-->[5] <!--[endif]-->Chen Y, Hao Y L, Liu F M. Analysis of effect and downward continuation of airborne gravity data[J]. Journal of System Simulation, 2008，20(8)：2190-2194．(成怡，郝燕玲，刘繁明．航空重力测量数据向下延拓及其影响因素分析[J]．系统仿真学报，2008，20(8)：2190-2194.)

<!--[if !supportLists]-->[6] <!--[endif]-->Hao Y L, Chen Y, Sun F, et al. Simulation research on Tikhonov regularization algorithm in downward continuation[J]. Chinese Joumal of Scientific Instrumen, 2008, 29(3): 605-609. (郝燕玲，成怡，孙枫，等．Tikhonov正则化向下延拓算法仿真实验研究[J]．仪器仪表学报，2008，29(3)：605-609.)

<!--[if !supportLists]-->[7] <!--[endif]-->Zhai Z H, Li N. Downward continuation of airborne gravity data based on land data[J]. Science of Surveying and Mapping, 2009, 34(Suppl.): 29-30. (翟振和, 李楠. 基于地面重力基准数据的航空重力数据向下延拓[J]. 测绘科学, 2009, 34(Suppl.): 29-30.)

<!--[if !supportLists]-->[8] <!--[endif]-->Deng K L, Bao J Y, Huang M T, 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. (邓凯亮，暴景阳，黄谟涛，等．航空重力数据向下延拓的Tikhonov正则化法仿真研究[J]．武汉大学学报（信息科学版），2010，35(12)：1414-1417.)

<!--[if !supportLists]-->[9] <!--[endif]-->Deng K L, Huang M T, Bao J Y, 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. (邓凯亮，黄谟涛，暴景阳，等．向下延拓航空重力数据的Tikhonov双参数正则化法[J]．测绘学报，2011，40(6)：690-696.)

<!--[if !supportLists]-->[10] <!--[endif]-->Jiang T, Li J C, Wang Z T, et al. Solution of ill-posed problem in downward continuation of airborne gravity[J]. Acta Geodaetica et Cartographica Sinica, 2011，40(6)：112-122. (蒋涛，李建成，王正涛，等．航空重力向下延拓病态问题的求解[J]．测绘学报，2011，40(6)：112-122.)

<!--[if !supportLists]-->[11] <!--[endif]-->Zhou B Y, Luo Z C, Xu C, et al. Comparison among fast fourier transform algorithms for downward continuation of airborne gravity data[J]. J. Geodesy and Geodynamics, 2013, 33(1): 64-68. (周波阳, 罗志才, 许闯, 等. 航空重力数据向下延拓的FFT快速算法比较[J]．大地测量与地球动力学，2013，33(1)：64-68.)

<!--[if !supportLists]-->[12] <!--[endif]-->Sun W, Wu X P, Wang Q B, 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. (孙文, 吴晓平, 王庆宾, 等．航空重力数据向下延拓的波数域迭代Tikhonov正则化方法[J]．测绘学报, 2014, 43(6): 566-574.)

<!--[if !supportLists]-->[13] <!--[endif]-->Wang X T, Shi P. Regularization methods and spectral decomposition for the downward continuation of airborne gravity data[J]. Acta Geodaetica et Cartographica Sinica, 2004 , 33(1): 33-38. (王兴涛，石磐．航空重力测量数据向下延拓的正则化算法及其谱分解[J]．测绘学报，2004，33(1)：33-38．)

<!--[if !supportLists]-->[14] <!--[endif]-->Jiang T, Dang Y M, Zhang C Y, et al. Downward continuation of airborne gravity data based on rectangular harmonic analysis[J]. Acta Geodaetica et Cartographica Sinica, 2013 , 42(4): 475-480. (蒋涛, 党亚民, 章传银, 等. 基于矩谐分析的航空重力向下延拓[J]. 测绘学报, 2013 , 42(4): 475-480.)

<!--[if !supportLists]-->[15] <!--[endif]-->Zeng X N, Li X H, Han S Q, et al. A comparison of three iteration methods for downward continuation of potential fields[J]. Progress in Geophys., 2011, 26(3): 908-915. (曾小牛, 李夕海, 韩绍卿等. 位场向下延拓三种迭代方法之比较[J]. 地球物理学进展, 2011, 26(3): 908-915.)

<!--[if !supportLists]-->[16] <!--[endif]-->Liu X G. Theory and methods of the Earth’s gravity field model recovery from GOCE data[D]. Zhengzhou: Information Engineering University, 2011.（刘晓刚. GOCE卫星测量恢复地球重力场模型的理论与方法[D]．郑州：信息工程大学，2011.）

<!--[if !supportLists]-->[17] <!--[endif]-->Zhong B, Liu H L, Luo Z C, et al. Reduction and gridded processing of satellite gravity gradient data[J]. J. Geodesy and Geodynamics, 2011, 31(3): 79-84. (钟波, 刘华亮, 罗志才, 等. 卫星重力梯度数据的归算与格网化处理[J]. 大地测量与地球动力学, 2011, 31(3): 79-84.)

上一篇：全球及局部海洋扰动重力反演的快速解析方法

下一篇：水准数据粗差探测的组合分析方法