航空重力数据向下延拓的迭代Tikhonov正则化法

冯淑萍，高延民

（1.西安测绘总站，陕西 西安 710054）

航空重力数据向下延拓的迭代Tikhonov正则化法

冯淑萍1，高延民1

（1.西安测绘总站，陕西 西安 710054）

根据观测面和延拓面测量数据的Poisson积分平面近似关系，结合快速傅立叶变换算法，将向下延拓转换到频率域进行计算，并采用迭代Tikhonov正则化方法，克服计算的不稳定性，提高计算结果的精度，实现了航空重力测量数据的向下延拓。最后采用模拟航空重力测量数据验证了该算法的有效性，取得了较好的延拓结果。

航空重力；向下延拓；正则化参数；迭代Tikhonov正则化法；快速傅立叶变换

在实际的工作中，航空重力测量常常是在起伏的航线上进行的，然而重力资料的定量解释方法要求测量数据分布在一个平面上，因此，需要将实测资料向下延拓到一个平面上。而向下延拓是一个典型的不适定问题[1-2]，主要表现为计算的不稳定性。随着向下延拓深度的增大，会对重力测量数据中的高频干扰信号起显著的放大作用，从而不能分辨有效信号。为了解决这个问题，本文采用迭代Tikhonov正则化方法来实现航空重力测量数据的向下延拓。

1 向下延拓的基础模型

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

图1 向下延拓示意图

图中，Δgh（ξ，η）表示观测面上的航空重力测量数据；Δg0（x，y）表示延拓面上的航空重力测量数据。计算观测面以下至场源以上z=0这个平面的航空重力测量数据称为航空重力测量数据的向下延拓。

根据航空重力测量数据向上延拓公式，可得观测面航空重力测量数据与延拓面重力数据之间的近似关系公式[3-6]：

式中，h是向下延拓深度；r是延拓面上点（x，y，0）与观测面上点（ξ，η，h）之间的距离。式（1）是第一类Fredholm积分方程，具有实对称核，且可以表示为二维卷积：

其中：

将式（2）转换到谱域里计算，其谱表达式为：

因为：

其中，f=（u2+v2）1/2；u、v分别表示空域变量x、y对应的频域变量。将式（5）代入式（4），可以得到：

通过变换可以得到：

对式（7）两端进行Fourier逆变换，得到延拓面重力数据：

从式（8）可以看出，向下延拓的基础数学模型原理比较简单，但是由于其向下延拓算子 的不稳定性对航空重力测量数据中的高频噪声有着显著的放大作用，在延拓深度较大时会导致延拓结果精度不高，要解决这种不适定问题，可引入正则化因子。因此，本文研究了迭代Tikhonov正则化法数学模型。

2 迭代Tikhonov正则化方法

式（2）可以修改为：

式中，K表示第一类Fredholm积分算子。对于求解不适定问题采用Tikhonov正则化是一个常用的方法，它指的是求解一个极小化的正则化泛函，即

式中，α为正则化参数，用于平衡不稳定性和光滑性。上式等价于如下的Euler方程：

K*为算子K的伴随算子，由于航空重力测量数据向下延拓的Fredholm算子K为对称的线性紧算子，有：

因此，式（11）可以表示为：

即

对式（14）两边同时作傅立叶变换，得到：

经过调整，可以得到：

考虑到式（5）、（6），式（16）可以变形为：

在式（17）两端进行Fourier逆变换，即可得到航空重力测量数据向下延拓的Tikhonov正则化法公式为：

由于Tikhonov正则化的饱和效应，使得正则化解与准确解的误差估计不能达到阶数最优，迭代的Tikhonov正则化对此进行了改进，其迭代形式如下：

将上式变形为：

根据数学归纳法，式（22）又可以写为：

根据式（20），式（23）可以改化为：

再根据式（5）、（6），经过化简，可以得到：

在式（24）两端进行Fourier逆变换，可得到航空重力测量数据向下延拓的迭代Tikhonov正则化法公式为[7～13]：

3 数值实验与结果分析

3.1 航空重力测量数据仿真

采用2 160阶的EGM2008地球重力场模型，对中国某地区的航空重力测量数据和地面重力数据进行仿真。将航空重力测量飞机的飞行高度设为3 000 m，分辨率设为5'×5'。图2和图3分别表示观测面理论航空重力测量数据和延拓面理论重力数据。为了验证本课题向下延拓算法的有效性，在观测面理论航空重力测量数据中加入均值为0、标准差为3 mGal的高斯白噪声，其等值线图如图4所示。

3.2 迭代Tikhonov正则化法实验结果

在对迭代Tikhonov正则化方法的延拓精度进行测试时，首先需要对延拓误差与正则化因子的关系进行验证。表1表示迭代Tikhonov正则化法在取得最小延拓误差时对应的正则化因子和迭代次数。因此，表1的统计结果进一步证实了所得出的结论。

表1 最小延拓误差对应的正则化因子和迭代次数/mGal

图2 观测面理论航空重力测量数据等值线/mGal

图3 延拓面理论重力数据等值线/mGal

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

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

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

确定合适的正则化参数是迭代Tikhonov正则化法获取最优解的关键，如果α远大于1，将导致逼近问题对原问题过于光滑，使计算结果与原问题的解相差甚远；相反，如果α远小于1，则逼近问题未能很好地改良算子的谱，正则化迭代法的效果不明显。

图5 迭代Tikhonov正则化法延拓结果等值线/mGal

[1] 栾文贵.场位解析延拓的稳定化算法[J].地球物理学报,1983,26(3):263-274

[2] 梁锦文.位场向下延拓的正则化方法[J].地球物理学报,1989,32(5):600-608

[3] 王兴涛,夏哲仁,石磐,等.航空重力测量数据向下延拓方法比较[J].地球物理学报,2004,47(6):1 017-1 022

[4] 罗志才,宁津生,晁定波.卫星重力梯度向下延拓的谱方法[J].测绘学报,1997,26(2):168-175

[5] 孙中苗.航空重力测量理论方法及应用研究[D].郑州:信息工程大学,2004

[6] 刘晓刚,李珊珊,吴星.卫星重力梯度数据的向下延拓[J].大地测量与地球动力学,2011,31(1):132-137

[7] TIKHONOV A N, ARSENIN V Y. Solution of Certain Integral Equations of the First Kind[J]. J Assoc Comput Mach, 1962 (9):84-97

[8] 成怡,郝燕玲,刘繁明.航空重力测量数据向下延拓及其影响因素分析[J].系统仿真学报,2008,20(8):2 190-2 194

[9] 郝燕玲,成怡,孙枫,等.Tikhonov正则化向下延拓算法仿真实验研究[J].仪器仪表学报,2008,29(3):605-609

[10] 邓凯亮,暴景阳,黄谟涛,等.航空重力数据向下延拓的Tikhonov正则化法仿真研究[J].武汉大学学报(信息科学版),2010,35(12):1 414-1 417

[11] 曾小牛,李夕海,韩绍卿,等.位场向下延拓三种迭代方法之比较[J].地球物理学进展,2011,26(3): 908-915

[12] 邓凯亮,黄谟涛,暴景阳,等.向下延拓航空重力数据的Tikhonov双参数正则化法[J].测绘学报,2011,40(6):690-696

[13] 蒋涛,李建成,王正涛,等.航空重力向下延拓病态问题的求解[J].测绘学报,2011,40(6):112-122

P223

B

1672-4623（2016）10-0053-03

10.3969/j.issn.1672-4623.2016.10.015

冯淑萍，工程师，主要从事航空重力研究。

2015-07-13。

项目来源：国家自然科学基金资助项目（41304022）；国家重点基础研究发展计划资助项目（61322201，2013CB733303）；高分专项青年创新基金资助项目（GFZX04060103-5-12）。