最小阻尼二乘算法

㈠ 非线性回归方程为什么不过样本中心

关于回归方程，存在线性回归方程和非线性回归方程，其目的在于最大程度的拟合一组数据，数据由自变量和因变量组成。当判断线性回归方程的相关系数r的绝对值落在0.75-1之间时，可以认为自变量与因变量相关，用最小二乘法建立回归方程。而对于非线性回归方程，不妨考虑将自变量进行合理转化换元转化为线性回归方程。当然，既然是拟合，多数情况下点都不落在直线上，样本中心点除外。网络知道非线性回归方程为什么不过样本中心?
怎么判断用线性回归还是非线性回归？
有奖励写回答共9个回答
假面
如果是你希望，就带上XX的假面
聊聊关注成为第21480位粉丝
优先选择线性回归，因为线性回归容易处理。也可以选择非线性回归。非线性回归很复杂，而线性回归的方法基本上前人已经完善的差不多了。
处理可线性化处理的非线性回归的基本方法是，通过变量变换，将非线性回归化为线性回归，然后用线性回归方法处理。
假定根据理论或经验，已获得输出变量与输入变量之间的非线性表达式，但表达式的系数是未知的，要根据输入输出的n次观察结果来确定系数的值。
(1)最小阻尼二乘算法扩展阅读：
多数非线性模型中，参数必须限制在有意义的区间内。指的是在迭代过程中对参数的限制。分为线性约束和非线性约束。线性约束中将参数乘以常数 但这个常数不能为其他参数或者自身。非线性约束中至少有一个参数和其他参数相乘或者相除或者进行幂运算。
线性回归模型经常用最小二乘逼近来拟合，但他们也可能用别的方法来拟合，比如用最小化“拟合缺陷”在一些其他规范里（比如最小绝对误差回归），或者在桥回归中最小化最小二乘损失函数的惩罚。
相反，最小二乘逼近可以用来拟合那些非线性的模型。因此，尽管“最小二乘法”和“线性模型”是紧密相连的，但他们是不能划等号的。

㈡ 最小二乘最优化反演方法

1.对称四极装置电阻率测深一维线性反演的基本思路

·第1步：建立目标函数

电阻率测深数据的反演方法可归结为使下面目标函数趋于极小：

地球物理数据处理基础

这里，NS是供电极距数；ρ_ai是第i个极距的实测视电阻率；M（M_j，j＝1，2…，NM）是模型参数（地层的电阻率或厚度），NM是预测模型参数个数；ρ_ci是由预测模型正演计算所得的第i个极距的理论视电阻率；α是范数，当α＝2时即为最小二乘法。

由于地球物理反问题存在多解性。因此在求解时，为减少反演的多解性，常在目标函数式（9－48）中对模型加入先验信息，如加入模型先验信息的目标函数变为

地球物理数据处理基础

式中：M_b为背景模型；C为光滑度矩阵；λ为拉格朗日系数。即要求预测模型正演结果与实测数据最接近，而且要求预测模型与给定的背景模型最接近并且光滑。显然拉格朗日

系数λ可以调控对预测模型的要求偏重于哪边。

·第2步：线性化

对式（9－48）利用泰勒展开，将模型的解在初始模型M处展开，并忽略二次及以上高次项，得

地球物理数据处理基础

要上式趋于极小，则对于各供电极距i（i＝1，2，…，NS）要满足下面的线性方程：

地球物理数据处理基础

·第3步：计算视电阻率对模型参数的偏导数 得到偏导数矩阵。

·第4步：解线性方程组（9－50），得到预测模型M的修改量ΔM，建立新的预测模型。

·第5步：计算新模型的正演理论视电阻率，与实测视电阻率进行对比，如果精度满足要求，则新的预测模型即为反演结果；否则回到第3步重新计算偏导数矩阵，重新计算模型修改量，直到精度满足要求。

上述反演过程中有三个问题需要解决：①偏导数的求解；②模型参数的处理，模型参数中有不同量纲的层电阻率和厚度，尤其电阻率参数变化范围很大，如果直接由上述方法求解，不但会导致方程（9－50）严重病态，而且电阻率和厚度参数的修改也不会正确，从而会导致反演方法不收敛；③线性方程组（9－50）的求解，线性方程组（9－50）往往是不相容方程组，求解方法不同，收敛速度和收敛性也不同。下面将对这三个问题的处理方法逐一进行介绍。

2.偏导数的计算

对称四极装置电阻率测深的视电阻率对模型参数的偏导 可以从数值滤波正演过程中直接给出。但一般反演中，偏导数通常用差分方法来计算。这里介绍用差分法来求对称四极装置电阻率测深的视电阻率对模型参数的偏导

取ΔM_j＝0.1M_j，则

地球物理数据处理基础

3.模型参数的无量纲化处理

为了解决视电阻率数据和电阻率模型参数变化范围大的问题，不少曾用它们的对数值来进行计算，但是这并没有解决层厚度参数与电阻率参数量纲不同的问题，因此，在反演过程中，出现了视电阻率测深曲线拟合很好，但预测模型与真实模型相差很大的情况，所以必须对方程（9－50）无量纲化。

将式（9－51）的偏导数代入方程（9－50）中可得

地球物理数据处理基础

其中：i＝1，2，…，NS。

上式两端同时除以ρ_ci（M₁，M₂，…，M_j，…，M_NM）可得

地球物理数据处理基础

这样方程组（9－52）可写为矩阵形式

地球物理数据处理基础

上式方程组中，左端系数矩阵A中各系数A_ij、未知向量x中各变量x_j以及右端向量B中各系数B_i都无量纲，从而解决了参数变量量纲不同的问题。

解线性方程组（9－54），则新模型参数Μ^*｛M^*_j，j＝1，2，…，NM｝为

地球物理数据处理基础

另外为了防止模型参数修改过量，实际过程中又作如下规定：x_j>1.5时取x_j＝1.5，x_j<－0.5时取x_j＝－0.5，即每次修改量不超过原有模型参数值的一半，以保证收敛稳定。

4.线性方程组求解

线性方程组（9－54）的求解问题，由于A为NS×NM阶矩阵，且对地球物理反演问题而言A多半是接近奇异的，所以不能用常规方法求解。解这个非方阵线性方程组最常用的方法有：

（1）奇异值分解法：它的基本思想是建立在奇异值分解定理上，即任意NS×NM阶矩阵A均可分解为A＝UWV^T，这里U为NS×NS阶正交阵和V为NM×NM阶正交阵，

地球物理数据处理基础

其中：δ₁，…，δ_r为矩阵A的奇异值，r是矩阵的秩。当A非奇异时，奇异值较大，方程组（9－54）有广义逆解X＝VW^－1U^T，这里

地球物理数据处理基础

当A接近奇异时，有的奇异值就较小，此时由于W^－1中系数过大，上述解的误差就较大。为了解决这个问题，维根斯（Wiggins）提出用最接近的矩阵R来代替A，R＝UW_eV^T，而

地球物理数据处理基础

W中小的奇异值在这里便被零代替了。因此有较精确的广义逆解X＝VW^－1_eU^T。

（2）最小二乘法：在方程组（9－54）两端左乘A^T，则方程变为A^TAx＝A^TB，当A非奇异时，有唯一解X＝（A^TA）^－1A^TB；当A接近奇异时，用上式求解误差就会很大，导致反演不收敛。为了解决这个问题，马奎特提出给系数矩阵A^TA的对角元素上加上一个常数，即A^TAx＋λI＝A^TB，以改善方程的条件，这就是着名的马奎特法又称阻尼最小二乘法。

分析奇异值（维根斯）分解方法和阻尼最小二乘（马奎特）方法，可以认为奇异值分解（维根斯）方法是将导致小的奇异值的方程取消来改善奇异值，而马奎特方法则是通过在A^TA的对角元素上加常数来增大A的奇异值，两者的思想刚好相反。

㈢ 最优化选择法数学原理

2.2.1 目标函数

设观测异常以ΔZ_k表示，k为观测点序号，k=1，2，…，m，m为观测点数。

设所选用的地质体模型的理论异常以 Z 表示，Z 是模型体参量和观测点坐标的函数，即

Z=f（x_k，y_k，z_k，b₁，b₂，…，b_n）

式中：x_k，y_k，z_k为观测点的坐标；b₁，b₂，…，b_n为模型体的参量，如空间位置、产状、物性等，参量的个数为n。

模型体的初始参量用

，

，

，…，

表示。

理论曲线与实测曲线之间的符合程度，是以各测点上理论异常与实测异常之差的平方和（即偏差平方和）来衡量的，用φ表示，即

地球物理数据处理教程

目的在于求得一组关于模型体参量的修改量δ₁，δ₂，…，δ_n，来修改模型体给定的初值参量，即

地球物理数据处理教程

于是求出关于模型体参量的一组新值，而由这组新参量形成的模型体的理论异常与实测异常之间的偏差平方和将取极小，即是

地球物理数据处理教程

代入式（2.2.1）中将使φ值获得极小，这时b_i即为我们的解释结果，这称为最小二乘意义下的最优化选择法。

我们称φ为目标函数，用它来衡量理论曲线与实测曲线的符合程度。最优化方法的关键在于求取使φ值获得极小参量的改正值δ_i，而f通常是b_i的非线性函数，因而该问题归结为非线性函数极小的问题。

2.2.2 求非线性函数极小的迭代过程

从上已知f为b_i的非线性函数，那么要求它与实测值之间的偏差平方和φ为极小的问题就称为非线性极小问题，或称为非线性参数的估计问题。如果是线性问题，参数估计比较简单，通常进行一次计算即可求出参数的真值，而对非线性问题，参数估计却要复杂得多，为了求解，通常将函数在参数初值邻域内展成线性（忽略高次项），即所谓的线性化，然后再求得改正量δ_i（i=1，2，…，n），由于这是一种近似方法，因而不可能使φ一次达到极小，而需要一个迭代过程，通过反复计算而逐步逼近函数φ的极小值。

图2.1 不同埋深时的重力异常

为了说明这个求极小的迭代过程，可以举一个单参量的例子，即假如我们要确定引起重力异常Δg_k的场源地质体的深度，假设场源为一个已知体积和密度的球体模型，如图2.1所示，那么φ就是球心埋深z的函数，如果球心埋深的真值为h，我们首先取初值为z^（0），这时函数

地球物理数据处理教程

式中：Δg_k为实测异常；g（z）是球心埋深为z的理论重力异常；φ随z的变化情况示于图2.2 中，要求使φ获极小的z，即要求使

地球物理数据处理教程

的根。由于z^（0）和φ（z^（0））不能一次求出φ的极小来，通常采用迭代的办法，如图2.3所示，例如用牛顿切线法迭代求根，根据下式

地球物理数据处理教程

得到一个更近似于根的值z^（1），但不等于h，因此需进一步再用上式，将z^（1）作为新的初值z^（0），可得到新的z^（1）更接近于h，如此反复下去可以使z值无限接近于h，当满足精度要求时，我们认为它近似等于h了，停止迭代，这时的z^（1）就作为h值。

图2.2 函数φ（z）随z变化示意图

图2.3 用牛顿切线法求φ′（z）=0的根示意图

2.2.3 单参量非线性函数的极小问题

单参量不仅是讨论多参量的基础，而且往往在求多参量极小时要直接用到单参量极小的方法，因此有必要作一介绍。

求单参量极小的方法很多，上面用到的牛顿切线法就是其中之一，在此我们介绍一种用得较多的函数拟合法，以及精度较高的DSC-Powell方法。

2.2.3.1 函数拟合法

2.2.3.1.1 二次函数拟合法

A.不计算导数的情况

设取三个参量值x₁、x₂、x₃，它们对应的φ 值就应为φ₁、φ₂、φ₃，过三个点（x₁，φ₁；x₂，φ₂；x₃，φ₃）作二次抛物线，应有下式

地球物理数据处理教程

联立φ₁、φ₂、φ₃的方程式，即可得出系数A、B、C来。

当A＞0时，应有极小点存在，我们设极小点为d，那么根据极小的必要条件有

地球物理数据处理教程

将A、B的表达式代入即得

地球物理数据处理教程

当x₁、x₂、x₃为等距的三点时，上式可简化为

地球物理数据处理教程

B.计算导数的情况

设已知两个点的参量值x₁和x₂对应的函数值φ₁、φ₂，并已求得x₁点的一阶导数值φ′（x₁），可用下列方法求极小点d：

地球物理数据处理教程

联立φ₁、φ₂、φ′（x₁）三个方程即可得A、B、C，代入极小点的表达式即可求得极小点。

为了简化起见，不妨设x₁为坐标原点（即x₁=0），设x₂=1，于是上面各式简化成：

φ′（x₁）=B

φ₁=C

φ₂=A+B+C

A=φ₂-φ′（x₁）-φ₁

则

地球物理数据处理教程

2.2.3.1.2 三次函数拟合法

取两个点的参量值x₁和x₂，及相应的φ₁和φ₂值，并已得到该两点的一阶导数值φ′（x₁）和φ′（x₂），我们选用一个三次多项式

φ=Ax³+Bx²+Cx+D

代入上面给出的4个条件，同样，为了简化起见，不妨设x₁为坐标原点（即x₁=0），设x₂=1，则有

φ₁=D

φ₂=A+B+C+D

φ′（x₁）=C

φ′（x₂）=3A+2B+C

联立求解，可定出4个系数A、B、C、D，按照求极小的必要条件

φ′=3Ax²+2Bx+C=0

当二阶导数

φ″=6Ax+2B＞0

时有极小存在，极小点d就为

地球物理数据处理教程

为了计算方便，令

v=φ′（x₁）

u=φ′（x₂）

S=-3（φ₁-φ₂）=3（A+B+C）

Z=s-u-v=B+C

W²=Z²-vu=B²-3AC

于是极小点d就可用下列形式表示：

地球物理数据处理教程

2.2.3.2 DSC-Powell 法

该法为比较细致的单参量探测法，精度比较高，计算工作量较大，大致可分为两部分来完成，其探测（迭代）过程如图2.4所示。

2.2.3.2.1 确定极小值所在的区间

采用的是一种直接探测法，做法可归纳如下。

第一步：给定探测方向x、初值点x₀和初始步长Δx，计算φ（x₀）和φ（x₀+Δx），若φ（x₀+Δx）≤φ（x₀），转向第二步；若φ（x₀+Δx）＞φ（x₀），则取-Δx为步长Δx，转向第二步。

第二步：计算x_k+1=x_k+Δx，计算φ（x_k+1）。

第三步：如果φ（x_k+1）≤φ（x_k），以2Δx为新步长代替Δx，且用k代表k+1，转向第二步。

如果φ（x_k+1）＞φ（x_k），则以x_m表示x_k+1，以x_m-1表示x_k，将上步的x_k作为x_m-2，并计算

地球物理数据处理教程

第四步：在4个等距点（x_m-2、x_m-1、x_m+1、x_m）中，去掉四点中离φ（x）最小点最远的那一点，即或是x_m，或是x_m-2，剩下的三点按顺序以x_α、x_b、x_c表示，其中x_b为中点，那么（x_α，x_c）区间即为极小值所在的区间。

2.2.3.2.2 用二次函数拟合法求极小点

将上面已确定的等距的 x_α、x_b、x_c三点及 φ 值，用二次函数拟合法即用公式（2.2.3）求得极小点，令为x^*点。再将x_α、x_b、x_c、x^*四点中舍去φ值最大的点，剩下的点重新令为α、b、c，则又得三点和它们相应的φ值，用公式（2.2.2）求其极小点x^*，如此反复使用公式（2.2.2），逐步缩小极小值的区间，一直到两次求得的极小点位置差小于事先给定的精度为止，x^*点即为极小点。

图2.4 DSC-Powell法示意图

2.2.4 广义最小二乘法（Gauss 法）

重磁反问题中的最优化方法，一般是指多参量的非线性最优估计问题，理论模型异常z=f（

，b₁，b₂，…，b_n）是参数b_i（i=1，2，3，…，n）的非线性函数，其中

=（x，y，z）为测点的坐标。由前已知ΔZ_k（k=1，2，…，m）表示在第k个观测点

上的实测异常，现在要寻求与观测异常相对应的理论模型的参量值b_i（i=1，2，…，n），使理论异常与实测异常的偏差平方和

地球物理数据处理教程

为极小。

设b_i的初值为

，则上述问题，即是要求修正量δ_i，使

地球物理数据处理教程

代入φ中，使φ获得极小。

高斯提出了首先将f函数线性化的近似迭代方法，即将f在

处按台劳级数展开取其线性项。

地球物理数据处理教程

式中

，

地球物理数据处理教程

当

和

给出后，

和

均可直接计算出来。将台劳展开式代入式（2.2.6）中，目标函数φ为

地球物理数据处理教程

要求

使φ取得极小，根据极小的必要条件

地球物理数据处理教程

将上式化为

地球物理数据处理教程

写成方程组形式

地球物理数据处理教程

式中：

，

（i，j=1，2，…，n）

再写成矩阵形式，有

地球物理数据处理教程

即

地球物理数据处理教程

其中

A=P^TP

地球物理数据处理教程

式中：P称为雅可比（Jacobi）矩阵，是理论模型函数对参量的一阶导数矩阵。A为正定对称矩阵，实际计算时，当实测异常值已给出，模型体的初值

已选定后，A和

即可计算出，求解方程（2.2.7）即可求出

，从而可得

。

上面推导出的方程（2.2.7）是将f线性化所得，因而只有当f为真正的线性函数时，

才是真正的极小点

，即一步到达极小；当f为非线性函数时，台劳式线性化仅为近似式，近似程序视

的大小而定，当|δ_i|较大时，二次以上项忽略的误差就大，反之就小，所以对于非线性函数

不能简单地作为极小点

，一般将

作为新的初值

再重复上述做法，再解方程（2.2.7）又得到新的

，反复迭代下去，直到满足精度要求为止（例如|δ_i|小到允许误差）。

在高斯法应用中常常出现一种困难，即迭代过程不稳定，当

过大时，台劳展开的高次项太大而不能忽略时，就可能发生这样的情况，即用方程（2.2.7）求得的解，得到的参量

所对应的φ值大于

所对应的φ值，那么它将不能稳定地收敛于φ的极小值，即是出现了发散的情况，一般说来当f非线性程度越明显时，越易出现发散的情况。

因此高斯法的一种改进形式如下，即不直接把

作为校正值，而将它作为校正方向，记为

，而在该方向上用单变量求极小的方法寻找在这个方向上的极小点，即寻找一个α，使目标函数φ（

）为极小，取

作为新的初值，再继续迭代（0＜α＜1）。

把这个改进的方法称为广义最小二乘法，它使迭代过程的稳定性有所改善，即使这样当初值取得不好时，也有可能出现不收敛。

2.2.5 最速下降法

从前述已知，我们的目的是要求目标函数的极小，高斯法是利用将f函数线性化，建立一个正规方程（2.2.7）来求取修正量的，最速下降法是另一类型方法，它直接寻找φ函数的下降方向来求取修正量，所以它又称为直接法，而高斯法又称为间接法。

从目标函数φ出发来寻找其下降方向

地球物理数据处理教程

始终是大于或等于0，因此它一定有极小存在，我们首先考虑初值点

的一个邻域内，将φ在

处台劳展开取至线性项，有

地球物理数据处理教程

希望寻找使Φ下降的方向，即要找新点

，使φ（

）＜φ（

）

即要求φ（

）-φ（

）＞0，

且越大越好，那么可得

地球物理数据处理教程

地球物理数据处理教程

式中

表示φ函数对

的各分量的导数所组成的向量，即梯度向量。

要使上式取极大，有

地球物理数据处理教程

上式说明了φ值下降最快的方向

，应该是与梯度方向

相反的方向，即负梯度方向，那么修正量就应在负梯度方向上来求取。下面讨论从

出发，沿负梯度方向上求取极小点的方法，除了用前面介绍过的方法外，在此再介绍一种近似计算方法。

要求从

出发，沿-

方向的极小点，即要求λ使φ

为-

方向上极小点。根据极小必要条件，有

地球物理数据处理教程

如果φ为二次函数时，λ可以直接解出，在重磁反问题中φ为非二次函数，且函数形式较复杂，一般无法直接解出λ，而采用近似法，先将φ（

）台劳展开，取至线性项，即

地球物理数据处理教程

假设粗略认为φ的极小值为零，则极小点的λ应有

地球物理数据处理教程

这个方法计算简单，但误差较大，特别是

远离真正极小点

时，φ值较大，上式的假设不适合，当接近真极小点

附近时，可以采用。但在重磁反问题中，由于实测值Z_k中含有干扰成分，所以即使到了

附近，φ值仍不会为零，因而上述计算λ的方法不能直接采用，可将上述计算的λ作为一个区间估计值，再用其他方法计算［0，λ］之间真正的λ值。

从上所述可将最速下降法叙述如下：从初值

出发，沿着φ（

）的负梯度方向-

（

）寻找极小点

，然后又从

出发，沿着φ（

）的负梯度方向-

（

）寻找极小点

，一直迭代下去，直到找到

为止。

由于这个方法是沿着初值点的最快下降方向，在该方向上如果采用单方向求极小的方法得到该方向上的极小点，那么又称“最优”、“最速”下降法。但需要指出的是，所谓“最速”是就初值点的邻域而言，所谓“最优”是指在初值点的负梯度方向上，所以它的着眼点是就局部而言，就初值点邻域而言，而对整体往往是既非“最优”，又非“最速”，而是一条曲折的弯路，难怪有人称它为“瞎子下山法”，如图2.5所示，当φ的等值面为拉长的椭球时更是如此。但它有一个十分可贵的优点，即在迭代的每一步都保证φ值下降，所以它是稳定收敛的，在φ函数复杂时，计算工作量较大些，对于大型计算机比较适用。

图2.5 最速下降法迭代过程示意图

图2.6 修正量的方向

2.2.6 阻尼最小二乘法（Marguardt）

比较上述两种方法可知，Gauss法修正量的步长大，当φ近于二次函数，可以很快收敛，但当φ为非二次函数，初值又给得不好时，常常引起发散。而最速下降法却能保证稳定的收敛，但修正量的步长小，计算工作量大。当φ的等值面为拉长的椭球时，Gauss法的修正量

和最速下降法的修正量

之间的夹角γ可达80°～90°，如图2.6所示。

对于φ为二次函数的情况下，高斯法的修正量

方向是指向φ的极小点，而最速下降法修正量

的方向是垂直于通过

点的φ函数等值面的切平面。因而当φ为比较复杂的函数时，有可能使

出现发散而失败。

阻尼最小二乘法是在Gauss法和最速下降法之间取某种插值，它力图能以最大步长前进，同时又能紧靠负梯度方向，这样既能保证收敛又能加快速度。它的基本思想是：在迭代过程的每一步，最好尽量使用Gauss法修正量方向

，以使修正步长尽可能地增大，如当这种情况下不能收敛时，再逐步改用接近最速下降的方向

，同时缩小步长，以保证收敛，下面以

表示由阻尼最小二乘法得出的修正量。

实现上述思想只要将方程

地球物理数据处理教程

改变为

地球物理数据处理教程

就能实现了。式中

为我们所要求的修正量，即称Marguardt修正向量，I为单位矩阵，λ是用来控制修正方向和步长的任意正数，又称阻尼因子，它起到阻止发散的作用，方程（2.2.9）中

显然是λ的函数，即

地球物理数据处理教程

通过这一改变后，即原来的正规方程（2.2.7）系数矩阵的主对角线上加一正数，从而使条件数得到了改善。如果原来A是奇异的，而A+λI可成为正定的，设原来A的最大特征值和最小特征值为μ_max和μ_min，则条件数就发生了如下变化：

地球物理数据处理教程

使病态条件数改善，对于计算来说，是十分有利的。

从方程（2.2.7）可看出，右端项为

地球物理数据处理教程

而φ的负梯度向量

的第i个分量

地球物理数据处理教程

所以

，即方程（2.2.7）、（2.2.9）的右端项

的方向即为负梯度方向，值为负梯度值的一半。

在方程（2.2.9）中，当λ=0时，即是（2.2.7）方程，这时

就是

；当λ→∞时，δ₀→

，而

是负梯度方向，这时

就是最速下降方向，所以阻尼最小二乘法的修正量

，是最速下降修正量

和Gauss法修正量

之间的某种插值，λ就是这种插值的权系数。

Marguardt向量

具有以下三个特性：

（1）当λ越来越大时，

的长度越来越小，且

地球物理数据处理教程

‖

‖表示

向量的范数，也即是它的长度。

（2）当λ由零逐渐增大时，

的方向逐渐由Gauss法的方向

转向最速下降法方向

，λ越大，

方向越接近

方向。

（3）对λ＞0的任意正数，

（满足方程（2.2.9））使φ在半径为‖

‖的球面上取得极小。

图2.7Δ₀（λ）随λ的变化情况示意图

以上三个性质说明，当λ逐渐增大时，

的方向由

向

靠近，它的大小‖

‖逐渐减小，λ→∞时，‖

‖→0，如图2.7所示。因此在迭代的任何一步，我们总可以找到充分大的λ，来保证稳定的收敛，因为当φ 不下降时，就加大λ向

靠，一直到使φ下降为止，从而保证收敛。性质（3）说明在跨出同样的步长时，以

（λ）方向最好，这就保证了该法的优越性。在实际计算时，总是在保证收敛的前提下，取较小的λ，以获得较大的步长前进。

下面介绍阻尼最小二乘法的迭代步骤，即实际计算过程。

（1）给出模型体参量初值

，计算φ（

）；给出实测场值ΔZ_k（k=1，2，…，m）；给出阻尼因子的初值λ^（0）及改变λ的比例系数v。

（2）开始迭代，λ=λ^（0）/v

（3）计算A，（A+λI）及右端项

在初值点

的值，得方程（2.2.9），（A+λI）

的系数矩阵及右端项。

（4）求解方程（2.2.9）得

。

（5）计算

及φ（

）。

（6）比较φ（

）和φ（

）。

若φ（

）＜φ（

），则该次迭代成功。判断

是否满足精度要求，若满足停止迭代，这时的

即为极小点

；若不满足精度要求，则将

作为新

，φ（

）作为新φ（

），减小λ作为新的λ^（0），转向第（2）步，继续迭代下去。

若φ（

）＞φ（

），则该次迭代失败，增大阻尼因子λ，将λ·v作为新的λ，转向第（4）步，即重新求解（A+λI）

方程，重新得到新的

。

该方法中阻尼因子λ的选择十分重要，上述选法是一种简单可行的方法，还有很多不同的选择方法，可参阅有关的书籍。

㈣ 预测模型建立

(一)参数拟合原理

在得到单井涌水量与所测量的地球物理测井各种参数之间的关系方程之后,可以发现里面还有很多待定的常数,这些常数在各种不同的地方是不一样的,为了能够确定这些系数,就需要获得这个地区的单井涌水量和对应的测井参数,然后拟合得到对应于这个地区的待定参数,这个被称为参数拟合。本程序所采用的拟合方法是改进型阻尼最小二乘法进行多参数数据拟合^[14]。下面介绍一下拟合方法的原理。

设按上述任一模型计算得到的第i个孔的单位涌水量为q_i。抽水实测单位涌水量为q_j,由前述诸个模型可见,q_j是个非线性多元变量函数,因而采用下述两种函数作为目标函数。并用最优化方法求取选定模型的待定系数是适当的。

(A+λ²K)ΔP=B

(1)目标函数取各井单位涌水量相对误差的平方和

含水层含水量预测综合物探技术

式中:λ为阻尼系数。

(2)目标函数取各井单位涌水量绝对误差的平方和

含水层含水量预测综合物探技术

选用哪种目标函数,应根据预测区各井单井涌水量的差异大小以及预测要求而定。若涌水量差异较大,而对涌水量较小者的预测精度要求较高,则宜选择相对误差的平方和作为目标函数,此时,小水量钻孔的预测精度虽然提高了,但大水量钻孔的预测精度相对降低了。若涌水量变化较小,且对涌水量较小者并不要求与大水量钻孔有相同高的预测精度,则适宜采用绝对误差的平方和作为目标函数。拟合流程见图5-4。

(二)模型构建

使用最小二乘准则,待求的模型系数a、b、c、d、e、f、g、R的值,应使得目标函数取极小值。显然,这是个非线性多元变量函数求最小二乘极小的问题,可采用最优化方法中比较有效的马奎特法(或称阻尼最小二乘法)求解,通常经过几次迭代就可求得各个模型的待定系数。

马奎特法是最优化中求最小二乘极小解比较有效的算法,它比梯度法、共轭梯度法收敛快,又比高斯牛顿法稳定,因而早已在很多其他反演解释中得到广泛应用。

经典马奎特算法中,由模型系数组成的矢量及其修正量的各元素相互间差别很大时,阻尼系数必将取得较大,这将增加迭代次数,降低运算速度,同时他还要求模型系数初值应靠近极小点,否则不易收敛,也就是说稳定性不理想。因此,我们采用加权阻尼因子的方法,即将经典马奎特方程中的单位矩阵K修改为与模型系数的大小有关的对角阵K,效果是模型系数大,阻尼小;模型系数小,阻尼大。从而使各模型系数以同等速度向极小点收敛,提高了算法的运算速度与稳定性,这就是改进的阻尼最小二乘法,其方程为

含水层含水量预测综合物探技术

图5-4 多参数拟合流程图

含水层含水量预测综合物探技术

利用上述拟合方法所求取的预测模型的待定参量a、b、c、d、e、、fg、R代入(5-61)式,便得到利用地球物理测井电阻率参量预测含水层含水量模型。

㈤ 二维断面反演的阻尼最小二乘法

二维视电阻率断面反演是利用网格剖分，在不同的网格单元内设置不同的电阻率，通过反演模型的求解过程，不断修改各单元的电阻率，最后由其分布可确定地电断面的几何特性与物理特性，进而指明所测地电断面的地质构造特征。

图2⁃1⁃74 二维断面反演剖分示意图

设在被研究的地电断面Ω域上，待求的地电模型参数为M个（即M个待求的地电阻率值ρ₁，ρ₂，…，ρ_M）。通常，我们是将Ω域按一定法则剖分，例如可按图2⁃1⁃74 进行剖分，每个网格单元对应一个电阻率。这些单元中的模型参量，可用向量ρ=（ρ₁，ρ₂，…，ρ_M）表示。实测视电阻率拟断面上，取N个视电阻率采样值（即

，i=1，2，…，N）。二维拟断面反演是不断修改地电模型的电阻率参数，使理论模型拟断面向实测拟断面逼近。在理论计算值

（i=1，2，…N）向实测视电阻率

逼近过程中，通过不断改变电阻率值ρ，使理论计算的ρ^理论与实测视电阻率ρ^实测之间的误差尽可能小（一般＜5%），以此作为衡量实测视电阻率ρ^实测和理论计算ρ^实测间拟合程度。通常采用对数型拟合方差F作为拟合视电阻率的目标函数（注意，以下的对数运算均系对电阻率的数值｛ρ_i｝_Ω·m进行），但简化记为

地电场与电法勘探

上式中

是根据初始模型参数（电阻率初值）正演计算的结果，它是地电断面参数ρ和电极距的函数，即

=ρ^实测（ρ，d_i）（d_i为与电极距有关的量）。目标函数F反映了实测拟断面数据与理论拟断面数据间的拟合程度，目标函数是模型参数的函数。二维视电阻率拟断面反演的目的就是要找到一组模型参数ρ=（ρ₁，ρ₂，…，ρ_M）使目标函数取得极小值，即

地电场与电法勘探

由于理论计算的ρ^理论（ρ，d_i）是模型函数ρ的非线性函数，故式（2⁃1⁃132）被称为非线性最小二乘问题，求取模型参数ρ拟合过程相当于数学上求多元函数极小值问题。对于非线性函数F直接求出ρ是很困难的，为此需要对非线性函数进行线性化近似处理。对假定的地电断面，给出一组模型参数初值ρ⁰=（

，

，…，

），将lnρ^理论（ρ，d_i）在初值ρ⁰附近做泰勒级数展开，将二阶及二阶以上的偏导数项略去，展开式的结果如下：

地电场与电法勘探

地电场与电法勘探

地电场与电法勘探

则式（2⁃1⁃133）可写为

地电场与电法勘探

将式（2⁃1⁃136）代入式（2⁃1⁃131），可得到目标函数F的近似表达式：

地电场与电法勘探

将式（2⁃1⁃137）右端记为�F，则非线性最小二乘问题（2⁃1⁃132）即可转换为线性最小二乘问题：

地电场与电法勘探

根据极值存在的必要条件，使

达到最小的Δρ_j（j=1，2，…，M），应满足下列方程组：

地电场与电法勘探

（k=1，2，…，M）

整理后得：

地电场与电法勘探

将上式写成矩阵形式：

A^TAΔρ=A^TAΔG（2⁃1⁃141）

式（2⁃1⁃141）即为目标函数（式2⁃1⁃138）的法方程。其中N×M矩阵A称为雅可比矩阵，其元素由式（2⁃1⁃134）来确定；ΔG为列向量，其元素为

地电场与电法勘探

求解法方程式（2⁃1⁃141），可得出模型参数的修正量Δρ，取ρ¹=ρ⁰+Δρ作为新的模型近似值，若

F（ρ¹）＜ F（ρ⁰）

且

F（ρ¹）＜ ε（ε为给定精度）（2⁃1⁃143）

则ρ¹作为二维地电断面的反演解释结果。若达不到精度，则以ρ¹取代ρ⁰重复以上过程，直至求出符合精度要求的模型参数为止。

以上求解过程的特点是将非线性最小二乘问题［式（2⁃1⁃132）］转化为求解一系列最小二乘问题［式（2⁃1⁃138）］，虽然每一步求得的Δρ只是

=min的极小元，还不能使F（ρ¹+Δρ）达到极小，但只要模型参数初值 ρ⁰选取得当，这种逐步线性化的过程是收敛的。

法方程组的系数矩阵A^TA，一般病态十分严重，甚至奇异。为保证反演过程收敛，增强法方程线的数值稳定性，可采用改进的阻尼最小二乘法（又名马奎特法）

（A^TA+λS）Δρ=A^TΔG（2⁃1⁃144）

上式中λ为阻尼因子；

S为对角矩阵，其形式如下：

地电场与电法勘探

总结以上阻尼最小二乘法的反演思路，我们可将其归纳为以下四个主要步骤：

（1）给出初值；

（2）计算理论拟断面（用2.5维数值模拟方法作正演计算）；

（3）解法方程：（A^TA+λS）Δρ=A^TΔG；

（4）让ρ¹=ρ⁰+Δρ⁰作为新模型参数，重新迭代反演。

可以看出，阻尼最小二乘法的主要计算工作量是解法方程，而要解法方程组，关键在于求出系数矩阵A（雅可比矩阵），A中元素用式（2⁃1⁃134）来计算。可见雅可比矩阵的计算是反演成像计算中很重要的一步。

计算雅可比矩阵的方法有多种，这里我们不再进行讨论。

㈥ 电测深曲线的数字解释法

近年来，用电子计算机对水平层电测深曲线进行数字解释发展较快，已经提出不少方法。其中用最优化法拟合电阻率转换函数或直接拟合ρ_s的解释方法用得较广，本节以修改的阻尼最小二乘法拟合T函数为例，介绍其原理（流程图见图2⁃1⁃37）。

图2⁃1⁃37 拟合T函数流程图

（一）用数字滤波法求阻率转换函数

将计算MN→0时ρ_s的基本公式（2⁃1⁃49 a）写为

地电场与电法勘探

应用傅里叶⁃贝塞尔积分公式：

地电场与电法勘探

中的第二式，即反转定理，可将（2⁃1⁃83）变为

地电场与电法勘探

对上式作变量代换，令

r=e^x，m=e^-y（2⁃1⁃86）

并令

T（y）=T₁（e^-y），ρ′_s（x）=ρ_s（e^x）（2⁃1⁃87）

便得到

地电场与电法勘探

式中

J′₁（y-x）=J₁（e^x-y）（2⁃1⁃89）

从数字滤波理论可知，（2⁃1⁃88）的褶积运算为空间域数字滤波的表达式，如将J′₁（y-x）看成为数字滤波器的脉冲响应，则当输入信号为ρ′_s（x）时，该滤波器的输出信号便是T（y），于是将根据实测视电阻率曲线计算电阻率转换函数的运算变成为数字滤波过程。

在野外测量ρ_s（r）曲线时，并非测连续曲线，而是在r等于某些数的极距上测量，这种测量称为空间采样测量，它所测到的是一系列的离散量。

在用电子计算机进行数字滤波时，ρ′_s（x）也是离散取值的，如将Δx称为取样间隔，而将iΔx表示第i个取样点的自变量值，则电极距对ρ′_s（x）离散采样取值时，所得采样系列为：ρ′_s（0）、ρ′_s（1Δx）、ρ′_s（2Δx），……ρ′_s（iΔx）……

野外实测曲线采样时，后一个极距大约是前者的1.5 倍，相邻极距间的间距是不等的，并且一般较计算机解释所用的步长要大，故需有内插程序，以便求出等步长的ρ_s值。此外，无论采用哪一种换算方法，均需曲线首尾支达到渐近线。通常实测曲线的左支与渐近线接近，而右支往往达不到渐近值，这便需根据实测曲线作外推延长。

在实际计算时，应将（2⁃1⁃88）式改写成离散取值的形式，根据对各种地电断面的ρ′_s（x）进行频谱分析可知，ρ′_s（x）的谱ρ_s（f）是有限分布的，即它具有某一截止频率f_z。根据采样定理，当采样间隔Δx≤

时，则ρ′_s（x）可用其离散取值ρ′_s（iΔx）表示：

地电场与电法勘探

代入（2⁃1⁃88）式便得：

地电场与电法勘探

令y=kΔx，s=x-iΔx，上式可写成：

地电场与电法勘探

式中

地电场与电法勘探

a〔（k-i）Δx〕称为滤波系数，在实际计算时，其个数是有限的，用l和L分别表示滤波系数的起止编号，并略去Δx不写，则（2⁃1⁃92）式可写为

地电场与电法勘探

滤波系数与电极装置类型有关，与地层参数无关，一旦将滤波系数求出，便可适合各种地电断面。格霍什（D.P.Ghosh，1971）首次发表了他用特殊方法求得的两组滤波系数，取样间隔为Δx=

ln10，其中12个一组的l=-3，L=8，顺次为：

0.0060、-0.0783、0.3999、0.3492、0.1675、0.0858、0.0358、0.0198、0.0067、0.0051、0.0007、0.0018，其后，不同作者发表了用不同方法求得的滤波系数，多者达百余个。

（二）用最优化法求层参数

用（2⁃1⁃94）式可求得实际的电阻率转换函数，以下用T_S（k）表示。用最优化法求层参数的大致过程是：先根据实际情况，给定一组层参数（叫做初值），算出T函数的理论值T_L（k），将它与T_S（k）比较，计算二者的误差，并根据它修改层参数；再算理论值，再作比较，再修改层参数，直到计算的T_L（k）与T_S（k）之差在规定范围内为止，便将这时理论值所对应的层参数作为解释结果。现稍为详细地讨论如下。

用P表示层数数ρ₁、ρ₂……ρ_n、h₁、h₂……h_n-1，用一组P便可算出T_L（k，P），理论值与实际值的相对误差可用δ表示：

地电场与电法勘探

式中m为极距的采样数，n为层参数个数。δ称为最优化的目标函数，最优化的过程便为使δ达到极小值的过程。

T_L（k，P）在P⁰点（初始值）写成泰勒级数，略去二次及二次以上的项，得到：

地电场与电法勘探

从而可将δ写为

地电场与电法勘探

增量ΔP称为步长，在给定层参数初始值的条件下，可将δ看作是ΔP（ΔP₁、ΔP₂、……ΔP_j……、ΔP_n）的函数。（2⁃1⁃96）式表示将层参数修改ΔP，使它从P⁰变为P′=P⁰+ΔP后的拟合差，因此，应确定δ取极小值的条件以选择ΔP。而多元函数取极小值的充分条件是全部偏导数为零，即

地电场与电法勘探

地电场与电法勘探

l=1，2，……n

此式可改写为

地电场与电法勘探

上式为n个未知函ΔP₁，ΔP₂，……，ΔP_n的n个线性联立方程，解此联立方程组，便可求得使δ到达极小点的最佳步长ΔP。但由于前面在作泰勒级数展开时，只取了一次项，即将T当作线性函数，然而T不是线性函数。故这样求得的ΔP，并不能使δ达到极小点，而只是向极小点靠近一步，达到了P′点（P′=P⁰+ΔP）。为了达到极小点，可再从P′点开始再用近似线性的办法形成新的线性方程组，又求得另一个步长ΔP^（2），对参数作第二次修正，如此反复逐次迭代逼近，最后使T_L（k，P）与T_S（k）拟合得最好，并将这时的T_L（k，P）所对应的层参数作为解释结果。上述过程便为最小二乘法的最优化过程。

在最小二乘法中，为了加快逼近过程，人们希望步长尽可能大些，然而步长太大，将T函数近似看作线性函数的误差也大了，有时甚至会出现无法逼近极小点而使计算失败的情况。为了解决这个矛盾，有人提出了修改的阻尼最小二乘法（或称修改的马奎特法），它的目标函数是：

地电场与电法勘探

其中R为阻尼系数，即在计算过程中，不但考虑到拟合差，而且还对步长加以约束，使之减少

倍，这样一来，使最优化过程的稳定性增加，而对速度影响也不大。

上述用数字滤波法求T函数并用最优化法拟合T函数的过程可用前面图2⁃1⁃37所示的计算流程图表示。

（三）用数字滤波求ρ_s及与实测ρ_s的直接拟合

直接将（2⁃1⁃83）式按（2⁃1⁃86）和（2⁃1⁃87）式的关系作变量代换，可将电阻率也变换为褶积形式：

式中

地电场与电法勘探

当离散取值时，上式也可以化为

地电场与电法勘探

式中

地电场与电法勘探

b为滤波系数。

（2⁃1⁃99）式为在已知T（k）的条件下，用数字滤波求ρ′_s（i）的表达式，一方面可用（2⁃1⁃99）式作正演计算，在给定层参数后，求出T（y），然后进行数字滤波便可求出ρ_s理论曲线。

另一方面，与前面拟合T曲线相似，用拟合ρ_s曲线来解反问题。即在给定层参数初值后，计算ρ_s理论曲线，并以之与实测的ρ_s拟合，逐次迭代后二者达到最佳拟合，从而求出层参数。

由于视电阻率与地层参数间的函数关系及其偏导数比较复杂，拟合ρ_s曲线较拟合T曲线要慢得多，但拟合T曲线的结果精度要低些，因此，有人提出先拟合T曲线，求得较佳参数作为初值，再进一步拟合ρ_s曲线，这样可兼顾速度与精度。

㈦ 马夸特反演法(阻尼最小二乘法)

对于线性系统:

d=Gm (3.40)

G为M×N矩阵，G的秩为r，当min(M，N)＞r，称为混定问题^[1]。解混定问题的常用方法是马夸特反演法，又称为“阻尼最小二乘法”。

建立如下目标函数:

ψ=(d-Gm)_T(d-Gm)+ε·m^Tm (3.41)

ε是一个大于零的数，称为阻尼因子，一般通过实验确定。对式(3.41)求对m的极小，得

(G^TG+ε·I)m=G^Td (3.42)

其中I为单位矩阵，由式(3.42)可得模型的解:

m=(G^TG+ε·I)^-1G^Td (3.43)

如果ε为零，上式就是常规的最小二乘法的解，如果ε不为零，相当于在G^TG的对角线增加一个大于零的数，这极大地改善了常规最小二乘法方程组的系数矩阵的条件数(条件数越小越好)^[1]:

地球物理反演教程

其中:令G^TG为A，λ为A^TA的特征值。马夸特反演法相当于加大原来的系数矩阵对角线的元素，在方程病态时能减小系数矩阵的条件数，使方程组稳定求解，因此在地球物理反演中有广泛的应用。

如果仅仅从解方程的角度也可以应用马夸特反演法，对任何如下线性方程组:

Am=B (3.45)

如果A的条件数差，都可以变为如下形式以获得稳定的解:

(A^TA+ε·I)m=A^TB (3.46)

㈧ 等效源法

在场论中，场和场源具有唯一的对应关系，但在实际中，观测到的场只是整个场的一部分，再加上观测误差和随机干扰，使其对观测值的解释不可避免地出现多解性，即是说有多种可能的场源分布与观测场在一定的误差范围内对应。与观测场对应的场源，如果不是真正的场源，我们称它为等效源。

在一般的资料解释中，多解性带来解释的困难，这里利用多解性来对位场进行转换，以突出有用信息，增加解释手段，对复杂异常的解释尤其重要。

3.1.1 基本原理

与观测场对应的多个场源中，选择一组最简单的场源，例如按一定位置分布或不固定位置的点荷、线荷，对磁异常还有磁偶极子、偶极线等，用最优化方法确定它们的质量或磁量，使它们产生在观测面上的场值与实测场值相吻合，利用这组等效源就能很方便地作各种位场变换，即是说观测值的各种变换值可以等价地认为是等效源产生的各种变换值，例如曲化平，向上、向下延拓，异常不同分量之间的换算，求变换磁化方向的磁异常（包括化到磁极），以及垂向一次、二次导数等等，均可用等效源的正演计算来求得。

因此，可以归纳出等效源法的特点如下：

（1）把各种繁杂的位场转换，变成一个简单的正演计算，计算过程简单，便于统一处理。

（2）不丢掉边部测点，条件好的情况下，可适当外推。

（3）对地形起伏较大的观测面，作位场转换的效果仍较好。

（4）由于等效源产生的场仅在观测面内与真实场源产生的场在一定误差范围内吻合，因而用等效源进行位场的转换，特别是向下延拓时就必然只限于一定范围之内，而不是整个空间，这就是该法的局限。

从上可知，该法应该分成以下两步进行：

（1）选择等效源模型，并用最优化方法求取等效源的质量或磁荷量。

（2）用等效源来计算各种位场转换值。

3.1.2 等效源的求取

图3.1 等效源示意图

以二维剖面垂直磁异常转换为例，如图3.1所示，选择等效源模型为磁偶极子，磁化方向的倾角为A，那么第i个等效源产生在第k个观测点上的垂直磁场值为

（ΔZ_mk）_i=m_i｛［2（z_i-z_k）²-（x_k-x_i）²］sinA-3（z_i-z_k）（x_k-x_i）cosA｝/［（x_k-x_i）²+（z_k-z_i）²］^5/2（3.1.1）

（ΔZ_mk）_i表示第i个等效源产生在第k个观测点上的磁场垂直分量；（x_i，z_i）为等效源的坐标；（x_k，z_k）为观测点的坐标（i=1，2，…，n，k=1，2，…，m）；共有 n个等效源，m个观测点，通常m≥n；如果我们记m_i的系数部分为c_ki，即

c_ki=｛［2（z_i-z_k）²-（x_k-x_i）²］sinA-3（z_i-z_k）（x_k-x_i）cosA｝/［（x_k-x_i）²+（z_k-z_i）²］^5/2

则式（3.1.1）简记为

（ΔZ_mk）_i=c_kim_i （3.1.2）

当取不同的等效源模型时，如取点磁荷，则式（3.1.2）的形式不变，仅是c_ki的表达式不同，点磁荷对c_ki的表达式变为

c_ki=（z_i-z_k）/［（x_k-x_i）²+（z_k-z_i）²］^3/2

按照场的叠加原理，第k个观测点上的磁场垂直分量应为所有等效源产生在该点上的场值之累加求和，即为

地球物理数据处理教程

按等效源的定义，等效源产生在观测点上的场值应与观测值相同，实测值为

，即已知测点有m个，可列出m个方程式如下

地球物理数据处理教程

上述方程组未知数为m₁，m₂，…，m_n共n个，而方程式有m个，因m≥n，所以为一超定方程组，一般求其最小二乘意义下的近似解，即目标函数 φ 为极小的解 m₁，m₂，…，m_n。

地球物理数据处理教程

从而归结为求解多元函数的极值问题，可用第二章所述的最优化方法求解。需要注意的是，当等效源位置固定时，Z_mk为自变量m₁，m₂，…，m_n的线性函数，可用直接解线性方程组的方法求解，但是我们发现，此线性方程组的系数矩阵C^TC的条件不好，易于奇异，所以我们采用以下最优化方法解此方程：

地球物理数据处理教程

以往我们采用最速下降法和阻尼最小二乘法求取等效源。最速下降法的做法是，给定一个初值

，迭代过程由下式给出

地球物理数据处理教程

式中

为

点处最速下降方向，即负梯度方向。

地球物理数据处理教程

α^（l）为第l次迭代的最佳步长。

由于Z_mk是

的线性函数，所以φ为标准的二次型函数，可以解析的求出最佳步长的表达式

地球物理数据处理教程

式中：f_k=Z_OK-Z_mk

用阻尼最小二乘法求取等效源，只需将式（3.1.6）改为下列形式，求解下列方程即可：

地球物理数据处理教程

做法类同于第二章所述，设初值为

，

为增量，求解式（3.1.9）得

，用

+

作为新的初值，继续求解新的

，迭代下去，直到满足精度要求为止。从式（3.1.2）可知，Z_m为

的线性函数，因而C矩阵（即Z_mk对

的一阶偏导数矩阵）为常数矩阵，所以C^TC在整个迭代过程中是不变的，它的各元素仅与等效源的位置及磁化方向和测点的位置有关，而与变量

无关，这是一个突出优点，可以减少计算工作量。

根据试验发现，用最速下降法求取的等效源

各分量之间相差较小，且有一定的规律，而阻尼最小二乘法求取的等效源

各分量之间相差较大，因而用最速下降法求出的等效源作位场转换的效果通常较好，但计算比较费时间，而用阻尼最小二乘法求出的等效源作位场转换的效果，特别是向下解析延拓的效果较差，但计算较快，如图3.2所示，图中示出了对于山脊地形下有一水平圆柱体的理论异常，用等效源法进行向下解析延拓，从图中可见延拓到第二层的效果，即用两种不同的最优化方法求出等效源

，用该

来进行向下延拓的效果，可见最速下降法在下延方面的效果优于阻尼最小二乘法，但时间却往往要超过十多倍。图中虚线表示阻尼最小二乘法的结果，点划线表示最速下降法结果，实线表示理论曲线。

根据计算发现，在地形起伏较大、向下延拓较深时，即使最速下降法也将引起较大的误差。其等效源位置的选择，较大地影响延拓效果，当等效源位置太浅时，易引起下延结果的振荡，而位置太深时，易引起求取等效源的系数矩阵奇异，下面简单给予说明。

为讨论方便，我们以点荷作等效源模型来讨论，这时式（3.1.2）中的c_ki应为

c_ki=（z_i-z_k）/［（x_k-x_i）²+（z_k-z_i）²］^3/2

当等效源理论z_i取得过大时，取极限情况有

地球物理数据处理教程

图3.2 两种最优化方法效果比较

图3.3 山脊地形磁异常下延效果图

1—理论值；2—奇异值分解法结果；3—马奎特法结果

可见，所有的c_ki都将趋于一个很小的数，从而使系数矩阵C很容易奇异，求出的解极不可靠。

为此，目前采用了公认的解病态方程组效果较好的奇异值分解法（附录B给出了该方法的基本原理）直接求解式（3.1.4），取得了很好的效果（图3.3所示）。于是，可将等效源的位置放置较深，以改善下延的效果。对于比较复杂的异常，如场值中包函有较多的高频成分的情况，还可采用双层排列等效源，一层分布较浅，另一层分布较深，如图3.3中的“·”所示。

山脊地形，如图3.3所示，坡角为26 °34′的山脊，最低层下面5个长度单位处，埋有一个截面为倾斜板的水平柱体，磁化倾角为45 °。排列了两层等效源，如图中“·”所示。同时用奇异值分解法和马奎特方法计算等效源参量求出最低层的延拓值，示于图上方。图中，实线为理论正演值，点划线为奇异值分解法的结果。尽管该模型的物体埋深较浅，效果仍较好，仅在两个山坡的中部地区，有一点小的波动。在下延较深的极值部分相差也甚小，最大相对误差δ_max为7.8%，均方差σ为2.52nT，整个曲线基本上与理论值符合。而马奎特方法的效果却较差，如图中虚线所示，极值减小很多，最大相对误差δ_max为23.6%，两个山坡中部地区有较大的波动，均方差σ为11.05nT。

图3.4 大台阶地形磁异常下延效果图

1—理论值；2—马奎特法结果；3—奇异值分解法结果

又如大台阶地形，这种地形通常下延效果较差。如图3.4所示，地形为一个63°26′的大台阶，高差为60个长度单位，在台阶的底层下10个长度单位处，埋有一个截面为倾斜板的水平柱体，磁化倾角为45 °。在该物体的顶部和近底部排列两层等效源，如图3.4 中“·”所示。利用奇异值分解法直接求解方程组（3.1.4），得到等效源参量，再计算向下延拓值，将其与理论正演值对比，示于图的上方。图中实线为台阶最低层的正演理论值，点划线为奇异值分解法计算的下延值，二者几乎完全重合，即使是右边下延深度相当大的部位，也相差很小，最大相对误差Δ_max为 9.1%，均方差 σ 为0.34nT，可见下延效果较好。

为了对比，用马奎特方法求解方程（3.1.9），得到该模型等效源参量计算的下延值，如图中虚线所示。显然在右边部分下延较深处，出现了较大的振荡现象，与理论正演值相差很大，最大相对误差Δ_max为133.3%，均方差 σ 为5.69nT。由此可见，奇异值分解法大大优于马奎特方法，而且计算时间二者几乎相同。

以上两例都表明，用奇异值分解法求等效源参量的分布，是比较有效的。这是因为等效源方程组是一个易于病态的方程组，用奇异值分解法求解病态方程能改善其结果的不稳定现象，从而改进了下延效果。试验结果还表明，适当加深等效源的位置，效果更好；采用双层等效源排列对高频成分引起的影响也有改善。

对于三维问题，等效源的分布在平面内，对不固定位置的等效源分布在某一空间范围内，基本原理和做法与二维类似，只是c_ki的表达式中加入了y方向的有关部分，例如取点荷为等效源模型时，有

地球物理数据处理教程

式中：（x_k，y_k，z_k）为观测点坐标；（x_i，y_i，z_i）为等效源的坐标（k=1，2，…，m；i=1，2，…，n）。

对三维问题主要困难是等效源较多，计算比较费时间，因而它的关键在于如何安排等效源的个数，以保证对整个位场的特征能较完整地反映出来。

3.1.3 利用等效源进行各种位场转换

如果已求得了等效源，将它当作等效的场源，那么对于观测值的各种位场转换值，在一定空间范围和一定误差范围内可以等价地认为与等效源所产生的相应场值相同，因而利用等效源作一系列的正演计算即可得到各种转换值，下面举几个例子来说明。

3.1.3.1 向上、 向下解析延拓及曲化平

不管是解析延拓还是曲化平实质上都是求观测面以上或以下某一位置上的场值，当利用观测值求出等效源

以后，要求空间任意一点l（x_L，z_L）的场值，可以利用式（3.1.3）简单地求得

地球物理数据处理教程

式中：Z_ml表示等效源在l点产生的场值，即作为观测值在l点的延拓值，式中的c_Li是将c_ki的表达式中（x_k，z_k）换成（x_L，z_L）即可。

3.1.3.2 磁异常不同分量之间的换算

一般是指将所测到的磁场垂直分量换算成水平分量，只要利用观测到的垂直分量求得等效源

以后，同样利用式（3.1.3）可求得观测面上或空间其他位置上的水平分量，只需将式（3.1.3）中的c_ki的表达式换成相应的水平分量的表达式，若要求的不是观测面上的水平分量，而是空间其他位置上的水平分量，只需将c_ki中的（x_k，z_k）换成所要求的那些点的坐标（x_L，z_L）即可得到

地球物理数据处理教程

式中：

即为水平分量表达式，例如当等效源模型取为偶极子时，X分量H_x中的

为

地球物理数据处理教程

当等效源取为其他模型时，

换为相应的形式。

3.1.3.3 求转换磁化方向后的场值

这时要求选用等效源模型体必须包括有磁化方向倾角，可选为偶极子、偶极线等，不能选用点荷、线荷。可计算转换磁化方向后原测点的场值，也可计算任意点上的转换磁化方向后的场值，因此我们设所要计算场值的点坐标为l（x_L，z_L），设A为磁偶极子磁化方向的倾角（A可取为实测地区的实际磁化方向的估计值），设转换磁化方向的磁倾角为A₀，那么转换磁化方向为A₀后的l点上的场值应为：

地球物理数据处理教程

式中：

=m_icos（A-A₀）

c_Li即为原来c_ki的表达式，将其中（x_k，z_k）换为（x_L，z_L）；A角换为A₀角。

若需要化磁极，取A₀=90°即可，若需要化到磁性层的顺层方向，可取A₀为磁性层的倾角。

3.1.3.4 求观测值的垂向一次、 二次导数

这时仍用式（3.1.3），只需将系数c_ki的表达式换成等效源模型公式对z求一次、二次垂向导数公式即可，例如等效源为偶极子时，其一次垂向导数的c_Li为

（c_Li）_z′=｛［6（z_i-z_L）³-9（z_i-z_L）·（x_L-x_i）²］sinA+［3（x_L-x_i）³-12（x_L-x_i）（z_i-z_L）²］cosA｝/［（x_L-x_i）²+（z_i-z_L）²］^7/2（3.1.15）

其二次垂向导数的系数c_Li为

（c_Li）_z″=｛［9（z_L-z_i）⁴-72（z_i-z_L）·（x_L-x_i）²+24（z_i-z_L）⁴］sinA+［24（x_L-x_i）³（z_i-z_L）+21·（x_L-x_i）²（z_i-z_L）-60（x_L-x_i）·（z_i-z_L）³］cosA｝/［（x_L-x_i）²+（z_i-z_L）²］^9/2（3.1.16）

除了以上列出的转换以外，还可以作一些其他转换。从上可知，所有的转换都变得十分简单，都可归纳为用式（3.1.3）作正演计算。

3.1.4 等效源法的其他应用

从上可知，等效源法不仅可用于位场转换，也可用于正演、反演等。这里介绍几个理论模型的应用实例。

3.1.4.1 复杂情况下的化极效果

在小比例尺大范围的化磁极资料处理中，会遇到如下的情况，即在一个大范围内，不同位置的地质体的磁化方向不同。这种情况对于常用的化磁极的方法（包括频率域的化极方法）来说是难以实现的。考虑到等效源法各自等效源之间具有一定的独立性，模拟变磁化方向比较容易，于是用等效源理论模型进行了变磁化方向的化极计算，并取得了满意的结果。

图3.5示出了变磁化方向化极的效果。图3.5（a）中为三个不同磁化方向的板状体（磁倾角分别为30 °、45 °、60 °）所产生的磁异常。这里排列了两层等效源，如图中“·”所示，让它们的磁倾角在三个板状体附近接近于真实倾角，通过拟合Δz曲线所得的等效源再进行化极计算得到如图3.5（b）所示的虚线。图中实线为这三个板状体垂直磁化理论曲线。两者对比可见，除极值部分有较小差异外，曲线基本吻合，最大相对误差为4.2%，均方差为3.32nT。

图3.5（a）不同磁化方向的板状体所产生的磁异常

图3.5（b）用等效源法化极计算曲线与垂直磁化理论曲线对比

1—图（a）中三个板状体垂直磁化理论曲线；2—用等效源法化极计算所得的Δz曲线

又选用了一个极端的情况，即当磁化倾角为零的化极模型。一般说来这种情况采用常用的化极方法效果最差，而用等效源法也取得了较为满意的结果。

图3.6（a）为一个零磁化倾角的板状体的正演Δz曲线。等效源分布为两层排列，一层分布在物体顶部深度内，另一层分布在近底部的深度上，如图中“·”所示。通过拟合Δz曲线所得的等效源再进行化磁极计算所得的结果如图3.6（b）中虚线。图中实线为理论正演垂直磁化曲线。两者几乎没有明显的差异，其均方差为0.06nT。

图3.6（a）零磁化倾角时板状模型正演曲线

图3.6（b）用等效源法化极计算与垂直磁化时曲线对比

1—垂直磁化时（a）图模型理论正演曲线；2—用等效源法化极计算所得的曲线

从上述结果可见，用等效源法进行化极计算的效果较好，它不受磁倾角大小和变磁化倾角的影响，这对于实际区域物探资料的化极处理，将是很方便的。

3.1.4.2 等效源反演效果

如前所述，既然等效源与真正场源在观测面内的场值相同，那么是否可以利用等效源来大概地圈定出真正场源的范围呢？

图3.7为用等效源的参量分布特征大致圈定场源范围的示意图。真正场源为一个二维倾斜板。在所研究的剖面范围内，设计一些均匀分布的等效源，拟合观测面内场值Δz，从而求得等效源的参量m，如图中的数值所示。可以看出，数值较大的参量集中在倾斜板的内部，在板的上顶面具有明显的数值突变现象，板的外部数值极小，而进入板内数值急剧增大，在板的左右边界虽不如上顶面那样明显地突变，但也是向着板的中部，数值逐渐变大。据此，大致可以圈出一个范围来，这个范围由下向上越来越清楚。下界面不如上界面明显，这是由于下界面的变化对场值的影响不如上顶面敏感。

图3.7Δz等效源反演参量分布与场源对比图

图3.8Δg等效源反演参量分布与场源对比图

图3.9 某矿区剖面用等效源反演效果图

1—实测值；2—等效源拟合曲线；3—地形；4—钻孔见矿位置

除了计算磁异常Δz以外，还可以计算重力异常Δg，图3.8示出了一个Δg的例子，真实场源为一个正方形水平柱体。在研究的范围内，设计了均匀分布的一些等效源，利用拟合观测面内Δg值的办法，求得等效源的参量

，如图中数值所示。可以看出，上顶面同样具有非常明显的突变，左右边界也具有较清楚的突变，仅仅是下底面不太明显，但仍然可以看出一些界面的迹象。在正方形内部具有参量的极大值，利用参量突变，基本上可以划出物体的边界。

对于一些复杂的不规则形体或复杂的地质条件，也可以考虑用等效源法圈定其边界。

图3.9是用等效源法进行某矿区磁测资料反演的结果，可见强值区域反映了矿体（强磁性体）的部位。

㈨ 自动反演解释方法

大地电磁测深曲线反演解释的任务是定量地求出实测视电阻率曲线所对应的地电断面参数，目前常用的是曲线自动拟合反演解释法。解释步骤和以前提及的相同，即首先给出一个初始模型参数，然后用其计算视电阻率理论曲线并和实测曲线进行对比，如果二者差别较大，则修改初始模型的参数，重新计算相应的理论曲线再作对比，直至实测曲线和理论曲线拟合最好，即二者之差满足给定的误差要求，这时理论曲线所对应的地电断面参数即为实测曲线的解释结果。

（一）反演解释的最优化问题

视电阻率曲线的反演解释，可归结为求满足方程：

地电场与电法勘探

的地电断面参数：

地电场与电法勘探

这里ρ_Ti表示实测视电阻率曲线在第i个周期点上的离散采样值，N表示采样数或周期点数，ρ_Li（λ₁，λ₂，…，λ_m）表示某一理论曲线在相应周期点上的视电阻率值，它是由给定的地电断面层参数ρ₁，ρ₂，…，ρ_n，h₁，h₂，…，h_n-1通过正演理论计算求得的，这里将地电断面层参数统一用向量：

λ=［λ₁，λ₂，…，λ_m］^T（3⁃2⁃57）

来表示，m表示地电断面层参数的总数，对于n层断面，m=2n-1。λ^*是理论曲线和实测曲线拟合最好时，理论曲线所对应的地电断面参数，即待求实测曲线的反演解释结果。

显然，Ψ是地电断面层参数λ的函数，称为评价函数或目标函数。使理论曲线和实测曲线相拟合，就是使二者在最小二乘意义下误差最小。曲线的反演解释问题，归结为求目标函数最小点所对应的地电断面参数。数学上称求目标函数极小点问题为最优化问题，所谓“最优”，是指在最小二乘意义下误差最小。

根据多元函数极值理论，使Ψ函数取得极小的必要条件是：

地电场与电法勘探

写成微商矩阵形式：

地电场与电法勘探

等式右端表示零向量。∇Ψ（λ）表示函数Ψ的梯度向量。梯度向量等于零的点，称为函数的稳定点，它可以是极大点、极小点或鞍点。

在稳定点对目标函数作台劳级数展开，并取至二次导数项，有

地电场与电法勘探

这里λ^*表示稳定点的坐标，且Δλ_j=λ_j-

，Δλ_k=λ_k-

。定义Ψ（λ）的二阶导数矩阵〔海色矩阵〕Q为

地电场与电法勘探

代入式（3⁃2⁃60），得：

地电场与电法勘探

这里引出向量

Δλ^T=［Δλ₁，Δλ₂，…，Δλ_m］.

如果式（3⁃2⁃62）中

Ψ（λ）-Ψ（λ^*）＞ 0，

则稳定点λ^*同时是局部极小点，这时式（3⁃2⁃62）中有：

Δλ^TQ（λ^*）Δλ ＞ 0（3⁃2⁃63）

因此，稳定点为局部极小点的充分条件应是该点处海色矩阵为正定的。

在m维空间中极小点附近目标函数的等值面：

Ψ（λ）=C，

它是m维空间的超椭球面。当m=2时是以极小点为中心的椭圆簇。因为二维空间中式（3⁃2⁃61）中的海色矩阵变为

地电场与电法勘探

矩阵Q正定的条件是：

地电场与电法勘探

由此构成的等值线方程：

地电场与电法勘探

是以（

，

）为中心的椭圆，即在极小点附近Ψ函数的“等高线”为同心椭圆簇。

（二）曲线自动拟合反演解释算法

从前面的讨论可以看出，视电阻率曲线的反演解释首先需要求出目标函数的稳定点，即求解方程组（3⁃2⁃59），这个方程组是非线性方程组，直接求解比较困难。这种类型的最优化问题称为非线性最小二乘问题。

目前，求解非线性最小二乘问题的方法，通常是给出一个初始模型：

地电场与电法勘探

然后对模型参数进行校正，使其逐步逼近极小点坐标并以此作为问题的近似解。为此必须求出模型的校正向量并对模型参数进行校正：

地电场与电法勘探

式中α称为步长因子（通常0＜α≤1），它决定沿校正方向的实际校正值。要使校正过程中目标函数是逐次下降的，即要求：

Ψ（λ^（k+1））＜ Ψ（λ^（k））

直至目标函数小于给定的质量控制指标δ（δ＞0）：

Ψ（λ^*）＜ δ

并以该模型参数作为问题的近似解。

显然，求解非线性最小二乘问题的关键在于求取校正向量。

目前比较常用求校正量（或修改量）的方法有：梯度法（亦称最速降法）、高斯⁃牛顿法（亦称最小二乘法）、马夸特法（亦称阻尼最小二乘法）和广义逆矩阵法等。它们的具体算法可参考有关文献（石应骏等，1985；陈乐寿，王光锷，1990），这里就不再一一介绍了。

另外，20世纪80年代初发展起来的连续介质一维反演方法，也是用计算机进行自动反演的。所谓连续介质系指地下岩层的电阻率是逐渐变化的。如果对一个模型若将电性层划分得足够薄，以致每层的厚度均可视为趋于零，那么每一深度皆对应一个电阻率值，这便构成了连续介质。目前，大地电磁测探的二维反演已日趋完善（杨长福，徐世淅，2005），三维反演也取得了重要进展（谭捍东等，2003）。