1 引 言

在实际工作中，航空重力与磁力测量通常是在起伏的航线上进行的。然而重力与磁力资料的定量解释方法一般要求测量数据分布在一个平面上，因此需要将实测资料换算到一个平面上。

重力与磁力测量数据的向下延拓直接关系到这些数据的处理及反演方法的应用效果。利用向下延拓可以把重力与磁力测量数据的观测面外推到场源体附近，突出局部场，有利于提高重力与磁力测量数据解释的可靠性。然而，向下延拓是一个典型的不适定问题，主要表现为计算的不稳定性。随着向下延拓深度的增大，会对重力与磁力测量数据中的高频干扰信号起着显著的放大作用，导致不能分辨有效信号。

目前处理不适定问题的主要方法是采用正则化技术，国内外学者也作了大量的研究。文献^[1]采用正则化技术研究了位场数据的向下延拓，并推导出了稳定的频率域延拓算子；文献^[2]讨论了正则化方法位场向下延拓的4个频率响应公式；文献^[3]介绍了解决非线性病态性问题的改进Landweber迭代法的收敛分析；文献^[4]评述了向下延拓及其基于Tikhonov正则化算法的正则化过程；文献^[5]研究了航空和地表重力数据向下延拓的几种正则化方法；文献^[6]基于Poisson积分方程，给出了向下延拓的正则化算法；文献^[7]采用正则化方法中的Landweber迭代法来处理不适定问题；文献^{[8, 9]}结合Tikhonov正则化算法和移去-恢复技术对航空重力测量数据进行了向下延拓的仿真试验；文献^[10]根据逆Poisson积分方程，采用L曲线法确定Tikhonov正则化参数；文献^[11]研究了位场数据的正则化向下延拓法，并考虑了几种有统计学意义的正则化方案；文献^[12]研究了迭代正则化方法(包括迭代Tikhonov法、Landweber迭代法和截断奇异值分解法)在位场下延中的应用；文献^[13]比较3种正则化参数的选取方法；文献^[14]推导了积分迭代法对应的正则化滤波子函数；文献^[15]探讨了3种空间域迭代解法在位场向下延拓中的异同和各自优势；文献^[16]提出了Tikhonov双参数正则化法；文献^[17]提出将广义岭估计用于求解航空重力向下延拓病态问题，研究了求解逆Poisson积分问题的3种正则化方法。

在向下延拓的正则化方法研究中，正则化参数的确定是一项非常重要的任务，参数选择的优劣直接影响了最终延拓结果的精度。目前，许多重力与磁力测量数据向下延拓方法都使用了正则化技术。然而，最佳正则化因子的确定，在模拟试验中是在延拓面数据已知的情况下根据最小延拓误差来估算的。而对于实际的重力与磁力测量数据向下延拓，延拓面并不一定存在已知的测量数据。因此，在这种情况下如何确定最优正则化参数，是应该重点研究的问题。总体来说，在目前的重力与磁力测量数据向下延拓的正则化技术研究中，有关正则化参数确定方法的研究较少。因此，最优正则化参数的确定问题研究具有重要意义。

本文根据观测面和延拓面重力与磁力测量数据的Poisson积分平面近似关系，结合快速傅里叶变换算法，将其转换到频率域进行计算。为了克服计算的不稳定性并提高计算结果的精度，引入Landweber正则化迭代法，在此基础上采用L曲线法研究磁力数据向下延拓中最优正则化参数的确定。

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

图 1 向下延拓示意图 Fig. 1 Downward continuation sketch map

图选项

u₀(x,y)表示观测面上的重力与磁力测量数据，u_h(ξ,η)表示延拓面上的重力与磁力测量数据，计算观测面以下至场源以上空间z=h平面的重力与磁力测量数据，这种换算称为重力与磁力测量数据向下延拓。

根据重力与磁力测量数据向上延拓公式，得到观测面重力与磁力测量数据u₀(x,y)与向下延拓面重力与磁力数据u_h(x,y)之间的平面近似关系为^[19]

点击浏览公式

对上式进行傅里叶变换，经过整理，得到向下延拓的基本原理公式

点击浏览公式

从式(2)可以看出，由于向下延拓算子e^2πfh的不稳定性，会对重力与磁力测量数据中的高频噪声有着显著的放大作用，在延拓深度较大时导致延拓结果精度不高，为了解决这种不适定问题，需要引入正则化因子，因此，研究了Landweber正则化迭代法数学模型。

式(1)可以修改为

点击浏览公式

由式(2)可知，重力与磁力测量数据向下延拓即为求解u_h(x,y)，并使得目标函数

点击浏览公式

点击浏览公式

点击浏览公式

对式(6)两端进行傅里叶变换，并由数学归纳法，经过调整，即可得到重力与磁力测量数据向下延拓的Landweber正则化迭代法公式为^{[3, 14, 15]}

点击浏览公式

利用式(7)进行重力与磁力测量数据的向下延拓时，如果是采用仿真数据，则可以根据观测面和延拓面的仿真数据来求解其最优正则化参数，但在处理实测数据时，最优正则化参数的确定就需要采用其他方法，本文使用的是L曲线法。

在实测重力与磁力测量数据的噪声水平未知的情况下，可以采用L曲线法来确定正则化参数。根据式(3)，L曲线法采用对数尺度描述拟合误差||Ku_h^α－u₀||₂与正则解||u_h^α||₂随参数α变化的对比曲线，它的特征是曲线呈明显的大写字母“L”形，并根据其轨点来确定正则化参数α^*[18]。令

点击浏览公式

点击浏览公式

点击浏览公式

将拟合误差||Ku_h^α－u₀||₂与正则解||u_h^α||₂转换到频率域，有

点击浏览公式

点击浏览公式

则对于Landweber正则化迭代法，根据式(7)，有

点击浏览公式

点击浏览公式

利用一阶和二阶差分可以得到ρ和θ的一阶和二阶导数。根据给定的α值，利用FFT算法，可快速计算得到相应的c(α)序列，绘成曲线后其极大值所对应的正则化参数即为确定的正则化参数α^*。

本文在数值试验时，以磁力测量数据为例来对最优正则化参数的确定方法进行研究。

利用球体磁场模型计算磁异常的公式为^[20]

点击浏览公式

为了增加磁力试验数据的复杂性，本文共设计了4个不同位置、不同埋深的球体磁场模型，并将其产生的磁异常进行叠加。球体磁场模型采用的基本参数如表 1所示。

球体磁场模型采用的球心坐标如表 2所示。

磁异常数据的格网分辨率为50m×50m，点数为221×221。根据球体磁场公式分别计算z=0m和z=300m两个不同观测面上的理论磁异常数据。图 2、图 3分别为观测面z=0m和z=300m的理论磁异常数据等值线图。为了验证本文确定最优正则化参数方法的有效性，在z=0m观测面理论磁异常数据中加入均值为0、方差为3nT的高斯白噪声，其等值线图如图 4所示。

图 2 观测面z=0m的理论磁异常数据等值线(单位：nT) Fig. 2 Isoline figure of theoretical magnetic anomaly data on the plane of z=0m(unit：nT)

图选项

图 3 观测面z=300m的理论磁异常数据等值线图(单位：nT) Fig. 3 Isoline figure of theoretical magnetic anomaly data on the plane of z=300m(unit：nT)

图选项

图 4 z=0m观测面加入3 nT高斯白噪声后的理论磁异常数据等值线图(单位：nT) Fig. 4 Isoline of theoretical magnetic anomaly data added with 3 nT white noise on the plane of z=0m (unit:nT)

图选项

其统计特性如表 3所示。

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

图 5 延拓误差与迭代次数的关系图 Fig. 5 Relationship figure of continuation error and iteration times

图选项

图 6 延拓误差与正则化因子的关系图 Fig. 6 Relationship figure of continuation error and regularization factors

图选项

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

从图 6可以看出，迭代次数越多时，在取得最小延拓误差处对应的正则化因子越小。另外，在正则化因子大于2.0时，可以看到延拓误差急剧增大，这与图 5显示的结果相符，说明在采用Landweber正则化迭代法进行磁场数据向下延拓时，正则化因子应该小于2.0。

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

图 7给出了Landweber正则化迭代法在取得最小延拓误差(即α=1.171)时对应的延拓结果等值线图。

图 7 延拓结果等值线图(单位：nT) Fig. 7 Isoline figure of continuation results(unit:nT)

图选项

从图 3、图 7比较看，Landweber正则化迭代法延拓结果对延拓面上理论磁异常数据的逼近效果很好。

图 8和图 9分别为迭代次数取30时Landweber正则化迭代法对应的曲率函数图和L曲线图，图 8中曲率最大值点即对应着最优正则化参数，也就是图 9中L曲线中的拐点处，有趣的是图 9中的曲线形状为倒L型。从图 8可知，α^*=0.732，将α^*代入式(7)，可以得到L曲线法选取的正则化因子对应的延拓结果，其等值线图如图 10所示。

图 8 曲率函数图 Fig. 8 Curvature function figure

图选项

图 9 L曲线图 Fig. 9 L-curve figure

图选项

图 10 L曲线法选取的正则化因子对应延拓结果的等值线图(单位：nT) Fig. 10 Isoline figure of the continuation results corresponding to regularization factor chosen by L-curve method(unit：nT)

图选项

由图 10和图 3的比较可以看出，Landweber正则化迭代法中，根据L曲线法确定的最优正则化因子，其对应的延拓结果对延拓面上的理论磁异常数据的逼近效果较好。

表 5为延拓面数据未知时利用L曲线法选取的正则化因子对应的延拓结果精度统计情况。通过与表 4比较可以看出，本文利用L曲线法来确定正则化因子，对于Landweber正则化迭代法来说精度较好。

从图 9可以看出，L曲线主要由一个“水平”的部分和一个接近“垂直”的部分组成。水平部分对应的解表示正则化参数太大，因此这部分解主要由正则化误差影响较大；垂直部分对应的解表示正则化参数太小，因此这部分解主要由观测数据误差影响较大。所以，水平部分和垂直部分分别表示过于平滑解和不太平滑解，只有曲线的拐点处则是对这两种解最好的平衡。

利用理论磁场反演实例表明：①当正则化参数小于2时，反演效果很好，当正则化参数大于2时，反演效果迅速变差；②当正则化参数处于[0.7，2]区间时，反演效率最好，因为只需10步反演就能取得精确的结果；③当正则化参数小于2时，反演效果均不错，而L曲线法确定的正则化参数(0.732)亦小于2，这也说明了L曲线法在确定正则化参数时的有效性。

本文采用L曲线法研究了Landweber正则化迭代法中最优正则化参数的确定问题，并通过模型磁测数据验证了L曲线法在确定正则化参数时的有效性。重力测量数据与磁力测量数据同属于位场数据，因此，L曲线法在重力测量数据的正则化向下延拓中也同样适用，只是其确定的正则化参数会有所不同。

对于正则化参数的确定，当观测数据噪声水平已知时，可以通过偏差原理、广义偏差原理、Arcangeli准则等方法确定；当数据噪声水平未知时，可以通过拟最优准则、广义交叉检验准则、L曲线准则等方法来确定。本文仅研究了L曲线法，还需要对其他方法的适用性进行试验。