阻尼最小二乘法_阻尼最小二乘法基本原理

个人意见：样本的标准差是表明离散程度,那么在标准差中是和什么进行比较的呢?答案是样本平均数.那么就说明所有的样本元素在过样本平均数的线性直线(即回归线)上侧或下侧波动.样本中心点即为样本平均数,故回归直线一定过样本中心点.『有条件详见理科数学书』

阻尼最小二乘法_阻尼最小二乘法基本原理

可线性化问题

处理可线性化处理的非线性回归的基本方法是，通过变量变换，将非线性回归化为线性回归，然后用线性回归方法处理。假定根据理论或经验，已获得输出变量与输入变量之间的非线性表达式，但表达式的系数是未知的，要根据输入输出的n次观察结果来确定系数的值。按最小二乘法原理来求出系数值，所得到的模型为非线性回归模型(nonlinear regression model)

非线性回归简介

如果回归模型的因变量是自变量的一次以上函数形式，回归规律在图形上表现为形态各异的各种曲线，称为非线性回归。[2] 这类模型称为非线性回归模型。在许多实际问题中，回归函数往往是较复杂的非线性函数。非线性函数的求解一般可分为将非线性变换成线性和不能变换成线性两大类。

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

·第1步：建立目标函数

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

地球物理数据处理基础

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

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

地球物理数据处理基础

式中：Mb为背景模型；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.偏导数的计算

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

取ΔMj＝0.1Mj，则

地球物理数据处理基础

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

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

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

地球物理数据处理基础

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

上式两端同时除以ρci（M1，M2，…，Mj，…，MNM）可得

地球物理数据处理基础

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

地球物理数据处理基础

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

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

地球物理数据处理基础

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

4.线性方程组求解

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

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

地球物理数据处理基础

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

地球物理数据处理基础

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

地球物理数据处理基础

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

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

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

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

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

3.1.1 基本原理

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

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

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

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

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

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

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

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

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

3.1.2 等效源的求取

图3.1 等效源示意图

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

（ΔZmk）i=mi｛［2（zi-zk）2-（xk-xi）2］sinA-3（zi-zk）（xk-xi）cosA｝/［（xk-xi）2+（zk-zi）2］5/2（3.1.1）

（ΔZmk）i表示第i个等效源产生在第k个观测点上的磁场垂直分量；（xi，zi）为等效源的坐标；（xk，zk）为观测点的坐标（i=1，2，…，n，k=1，2，…，m）；共有 n个等效源，m个观测点，通常m≥n；如果我们记mi的系数部分为cki，即

cki=｛［2（zi-zk）2-（xk-xi）2］sinA-3（zi-zk）（xk-xi）cosA｝/［（xk-xi）2+（zk-zi）2］5/2

则式（3.1.1）简记为

（ΔZmk）i=ckimi （3.1.2）

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

cki=（zi-zk）/［（xk-xi）2+（zk-zi）2］3/2

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

地球物理数据处理教程

按等效源的定义，等效源产生在观测点上的场值应与观测值相同，实测值为 ，即已知测点有m个，可列出m个方程式如下

地球物理数据处理教程

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

地球物理数据处理教程

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

地球物理数据处理教程

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

地球物理数据处理教程

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

地球物理数据处理教程

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

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

地球物理数据处理教程

式中：fk=ZOK-Zmk

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

地球物理数据处理教程

做法类同于第二章所述，设初值为 ， 为增量，求解式（3.1.9）得 ，用 + 作为新的初值，继续求解新的 ，迭代下去，直到满足精度要求为止。从式（3.1.2）可知，Zm为 的线性函数，因而C矩阵（即Zmk对 的一阶偏导数矩阵）为常数矩阵，所以CTC在整个迭代过程中是不变的，它的各元素仅与等效源的位置及磁化方向和测点的位置有关，而与变量 无关，这是一个突出优点，可以减少计算工作量。

根据试验发现，用最速下降法求取的等效源 各分量之间相差较小，且有一定的规律，而阻尼最小二乘法求取的等效源 各分量之间相差较大，因而用最速下降法求出的等效源作位场转换的效果通常较好，但计算比较费时间，而用阻尼最小二乘法求出的等效源作位场转换的效果，特别是向下解析延拓的效果较差，但计算较快，如图3.2所示，图中示出了对于山脊地形下有一水平圆柱体的理论异常，用等效源法进行向下解析延拓，从图中可见延拓到第二层的效果，即用两种不同的最优化方法求出等效源 ，用该 来进行向下延拓的效果，可见最速下降法在下延方面的效果优于阻尼最小二乘法，但时间却往往要超过十多倍。图中虚线表示阻尼最小二乘法的结果，点划线表示最速下降法结果，实线表示理论曲线。

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

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

cki=（zi-zk）/［（xk-xi）2+（zk-zi）2］3/2

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

地球物理数据处理教程 图3.2 两种最优化方法效果比较

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

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

可见，所有的cki都将趋于一个很小的数，从而使系数矩阵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。由此可见，奇异值分解法大大优于马奎特方法，而且计算时间二者几乎相同。

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

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

地球物理数据处理教程

式中：（xk，yk，zk）为观测点坐标；（xi，yi，zi）为等效源的坐标（k=1，2，…，m；i=1，2，…，n）。

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

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

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

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

不管是解析延拓还是曲化平实质上都是求观测面以上或以下某一位置上的场值，当利用观测值求出等效源 以后，要求空间任意一点l（xL，zL）的场值，可以利用式（3.1.3）简单地求得

地球物理数据处理教程

式中：Zml表示等效源在l点产生的场值，即作为观测值在l点的延拓值，式中的cLi是将cki的表达式中（xk，zk）换成（xL，zL）即可。

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

一般是指将所测到的磁场垂直分量换算成水平分量，只要利用观测到的垂直分量求得等效源 以后，同样利用式（3.1.3）可求得观测面上或空间其他位置上的水平分量，只需将式（3.1.3）中的cki的表达式换成相应的水平分量的表达式，若要求的不是观测面上的水平分量，而是空间其他位置上的水平分量，只需将cki中的（xk，zk）换成所要求的那些点的坐标（xL，zL）即可得到

地球物理数据处理教程

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

地球物理数据处理教程

当等效源取为其他模型时， 换为相应的形式。

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

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

地球物理数据处理教程

式中： =micos（A-A0）

cLi即为原来cki的表达式，将其中（xk，zk）换为（xL，zL）；A角换为A0角。

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

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

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

（cLi）z′=｛［6（zi-zL）3-9（zi-zL）·（xL-xi）2］sinA+［3（xL-xi）3-12（xL-xi）（zi-zL）2］cosA｝/［（xL-xi）2+（zi-zL）2］7/2（3.1.15）

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

（cLi）z″=｛［9（zL-zi）4-72（zi-zL）·（xL-xi）2+24（zi-zL）4］sinA+［24（xL-xi）3（zi-zL）+21·（xL-xi）2（zi-zL）-60（xL-xi）·（zi-zL）3］cosA｝/［（xL-xi）2+（zi-zL）2］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是用等效源法进行某矿区磁测资料反演的结果，可见强值区域反映了矿体（强磁性体）的部位。

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

对于线性系统:

d=Gm (3.40)

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

建立如下目标函数:

ψ=(d-Gm)T(d-Gm)+ε·mTm (3.41)

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

(GTG+ε·I)m=GTd (3.42)

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

m=(GTG+ε·I)-1GTd (3.43)

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

地球物理反演教程

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

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

Am=B (3.45)

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

(ATA+ε·I)m=ATB (3.46)

（一）视纵向电导解释法

根据均匀大地表面的理论研究，瞬变电磁测深的等效深度为

地电场与电法勘探

在不均匀介质中该式为

地电场与电法勘探

比例系数ατ依赖于地电断面参数和收-发距，用实验方法确定。

与频率测深方法相同，引入等效二层断面的概念。这时等效厚度和等效电阻率确定为

因此

地电场与电法勘探

对于远区场，356≤ατ≤503，而对于近区场ατ≤189。用与频率测深相同的方法可得视电阻率断面图或视纵向电导断面图。这时纵坐标可用ctgα=h/S。可见，真纵向电导值随α角的增加而减少（见图2-4-66）。在层状介质中，Sτ（hτ）曲线具有由不同电阻率分界面引起的一系列转折曲线。拐点纵坐标对应于这些层的上盘位置。利用倾角α或渐近线可确定每一层的纵向电导值。良导层的α值大，不太良导者α值小。对比相邻点相同α角或曲线形态，最后可追索同一地质界线。该断面实质上是时间断面，然而，若至少有一个钻井控制时，则可类似于地震勘探法，对比S（h）和Sτ（hτ）曲线，将时间断面转换为深度断面或地电断面。

表示。

视纵向电导曲线直观地给出断面中赋存良导地层情况。由于随着时间推移瞬变电磁场向地层深处传播，所以当出现新的导电层时，在纵向电导曲线上便出现数值增加。我们在下面所述的“视纵向电导”指的是，在断电后的某一时刻，电磁场的有效分布范围内地电断面能被具有一定深度的均匀导电平面等效代替时，后者的纵向电导。由此可见，这样的等效导电平面和其纵向电导以及等效深度均随时间而改变。

下面讨论由野外实测值如何确定视纵向电导和等效深度的方法。以S—MN装置（根据互换原理也可用于AB—s装置中）为例说明这一方法的实质。由（2-4-66）式可确定到深度h对应的瞬变时间t的纵向电导S。在实测的Eφ（t）曲线上取微分：

地电场与电法勘探

式中K=3PM/πr3，F（m—）=m—/（1+4m—2）5/2。然后再取一个新的辅助函数：

地电场与电法勘探

事先计算函数值：

地电场与电法勘探

形成数表或列线图。借助于列线图或数表（见前页表2-4-4），根据实测曲线的φ实（m—）值确定m—值，再由（2-4-66）式确定：

地电场与电法勘探

确定薄板参数的上述方法可适用于任一收-发距和时间。但须知，我们引入随时间增大而匀速“下沉”的等效导电平面，并用它代替均匀大地，从而求得均匀大地表面的瞬变场。如前所述，对于任意水平层状地电断面，仍然可以用此法求得瞬变场的近似解，但此时的等效平面随时间t增大而“下沉”的速度不再是常数，而是分别与各层介质的电导率有关。这样我们便可以在层状介质中引入视值，即

地电场与电法勘探

图2-4-65给出了ρτ（τ1/h1）和Sτ（τ1/h1）曲线的对比图。由图可看出，Sτ曲线达到渐近线的时间远早于ρτ曲线。这给野外工作带来很大方便。尤其在电源功率不够大的条件下开展工作时，这一优点更为突出。利用（2-4-76）式的关系可将（2-4-77）式变换为Sτ（hτ）关系曲线，且由一系列测点上的Sτ（hτ）曲线绘出其断面图。

图2-4-65 ρτ与Sτ曲线对比图

根据理论分析和实际经验可对Sτ（hτ）曲线和断面图作如下认识。Sτ（hτ）图接近由电测井资料建立起来的真S（h）曲线。如果围岩介质在电性上是均匀的，则Sτ（hτ）曲线为相对于h轴以固定角度α倾斜的直线，故（二）改进的阻尼最小二乘可行方向解释法

图2-4-66 Sτ（hτ）断面图

1.原理

改进的阻尼最小二乘法（简称LMF法）是我们所熟知的，为了尽量减少多解性影响，可引入可行方向算法。

瞬变电磁测深曲线拟合的最优化问题可以归结如下。

给定初值x（0），修正量为Δx，修正后的参量为x（1），则x（1）=x（0）+Δx。参量修正后，应使目标函数φ（x）为最小值，并设定优化应满足的约束条件，即

地电场与电法勘探

其中，φ（x）为目标函数，x=（x1……xm）T，为 xi的下界，为 xi的上界，i=1，……，m。根据阻尼最小二乘法的原理，可得法方程：

地电场与电法勘探

式中，A为雅可比矩阵，I为单位阵，α为阻尼因子。由方程组（2-4-80）解出的xi的值未必满足（2-4-79）式，为使优化参数满足约束条件，现引入可行方向法如下：

①对于可行域内的x（i0），利用LMF法求出Δxi。

②令

地电场与电法勘探

其中ε为给定小正数。

③令

地电场与电法勘探

这样求得的xi必在可行域内，即满足（2-4-79）式。

将利用可行方向法求得的x（1）i值作为下次迭代新的初值，如此迭代运算，直到求出xi，……，xm的最佳值。

2.反演方法

（1）目标函数和法方程的建立。为了在不同的采样时间上，使瞬变电磁测深曲线有相近似的拟合程度，取晚期视电阻率的相对误差为目标函数。

地电场与电法勘探

其中，ρτ（i）为第i个采样道的ρτ实测值；ρτ（i，P）为用模型参数P计算的第i个采样道上的ρτ理论值；n为总的采样道数；P为理论模型参数P=（ρ1……ρm，h1，……hm-1 ）T；m为地层层数；ρi和hi分别为第i层的电阻率和厚度。经台劳展开，取=0，得

地电场与电法勘探

其中A为雅可比矩阵

地电场与电法勘探

将（2-4-84）式的对角线项加上阻尼因子α，即形成LMF的法方程

地电场与电法勘探

约束条件为

地电场与电法勘探

按前述方法，对（2-4-85）、（2-4-86）式求解，所得的极小可行解必在可行域内。

（2）雅可比矩阵元素的计算。对于瞬变电磁测深问题，晚期视电阻率的理论函数ρτ（i，P）的一阶偏导数虽然存在，但难以用公式直接给出，但可采用简单差商公式代替一阶偏导数

地电场与电法勘探

其中ρτ（i，P）是给定一组P1，…，P2m-1求出的晚期视电阻率理论值，ρτ（i，P′）为给定一组Pj，…Pj+δPj，…，P2m-1求出的晚期视电阻率理论值，δPj为适当选取的差分步长。根据数值计算的实验结果，取δPj=（0.01～0.001）Pj比较合适。

（3）反演算例。设计模型为水平地层KH型四层断面。表2-4-5给出了该模型的真实参数、约束条件和反演结果，其中第二层和第三层具有S等值性，第二层和第三层的总的纵向电导为

S2+S3=80/200+10/200=0.45，而反演结果S2+S3=73.9/183.6+12.3/266.5=0.45。反演结果说明了方法的有效性。

表2-4-5 KH型四层断面的反演数据表 除上述反演解释方法外，还有一种根据“烟圈”理论提出的解释方法见文献（方文藻等，1993）。

关于二维、三维地电体的瞬变电磁测深反演方法，虽已取得一些初步成果，但目前仍处于试验研究阶段。

最优化选择法数学原理

2.2.1 目标函数

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

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

Z=f（xk，yk，zk，b1，b2，…，bn）

式中：xk，yk，zk为观测点的坐标；b1，b2，…，bn为模型体的参量，如空间位置、产状、物性等，参量的个数为n。

模型体的初始参量用 ， ， ，…， 表示。

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

地球物理数据处理教程

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

地球物理数据处理教程

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

地球物理数据处理教程

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

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

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

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

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

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

地球物理数据处理教程

式中：Δgk为实测异常；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.不计算导数的情况

设取三个参量值x1、x2、x3，它们对应的φ 值就应为φ1、φ2、φ3，过三个点（x1，φ1；x2，φ2；x3，φ3）作二次抛物线，应有下式

地球物理数据处理教程

联立φ1、φ2、φ3的方程式，即可得出系数A、B、C来。

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

地球物理数据处理教程

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

地球物理数据处理教程

当x1、x2、x3为等距的三点时，上式可简化为

地球物理数据处理教程

B.计算导数的情况

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

地球物理数据处理教程

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

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

φ′（x1）=B

φ1=C

φ2=A+B+C

A=φ2-φ′（x1）-φ1

则地球物理数据处理教程

2.2.3.1.2 三次函数拟合法

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

φ=Ax3+Bx2+Cx+D

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

φ1=D

φ2=A+B+C+D

φ′（x1）=C

φ′（x2）=3A+2B+C

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

φ′=3Ax2+2Bx+C=0

当二阶导数

φ″=6Ax+2B＞0

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

地球物理数据处理教程

为了计算方便，令

v=φ′（x1）

u=φ′（x2）

S=-3（φ1-φ2）=3（A+B+C）

Z=s-u-v=B+C

W2=Z2-vu=B2-3AC

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

地球物理数据处理教程

2.2.3.2 DSC-Powell 法

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

2.2.3.2.1 确定极小值所在的区间

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

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

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

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

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

地球物理数据处理教程

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

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

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

图2.4 DSC-Powell法示意图

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

重磁反问题中的最优化方法，一般是指多参量的非线性最优估计问题，理论模型异常z=f（ ，b1，b2，…，bn）是参数bi（i=1，2，3，…，n）的非线性函数，其中 =（x，y，z）为测点的坐标。由前已知ΔZk（k=1，2，…，m）表示在第k个观测点 上的实测异常，现在要寻求与观测异常相对应的理论模型的参量值bi（i=1，2，…，n），使理论异常与实测异常的偏差平方和

地球物理数据处理教程

为极小。

设bi的初值为 ，则上述问题，即是要求修正量δi，使

地球物理数据处理教程

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

高斯提出了首先将f函数线性化的近似迭代方法，即将f在 处按台劳级数展开取其线性项。

地球物理数据处理教程

式中 ，

地球物理数据处理教程

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

地球物理数据处理教程

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

地球物理数据处理教程

将上式化为

地球物理数据处理教程

写成方程组形式

地球物理数据处理教程

式中： ， （i，j=1，2，…，n）

再写成矩阵形式，有

地球物理数据处理教程

即地球物理数据处理教程

其中

A=PTP

地球物理数据处理教程

式中：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，

且越大越好，那么可得

地球物理数据处理教程 地球物理数据处理教程

式中 表示φ函数对 的各分量的导数所组成的向量，即梯度向量。

要使上式取极大，有

地球物理数据处理教程

上式说明了φ值下降最快的方向 ，应该是与梯度方向 相反的方向，即负梯度方向，那么修正量就应在负梯度方向上来求取。下面讨论从 出发，沿负梯度方向上求取极小点的方法，除了用前面介绍过的方法外，在此再介绍一种近似计算方法。

要求从 出发，沿- 方向的极小点，即要求λ使φ 为- 方向上极小点。根据极小必要条件，有

地球物理数据处理教程

如果φ为二次函数时，λ可以直接解出，在重磁反问题中φ为非二次函数，且函数形式较复杂，一般无法直接解出λ，而采用近似法，先将φ（ ）台劳展开，取至线性项，即

地球物理数据处理教程

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

地球物理数据处理教程

这个方法计算简单，但误差较大，特别是 远离真正极小点 时，φ值较大，上式的假设不适合，当接近真极小点 附近时，可以采用。但在重磁反问题中，由于实测值Zk中含有干扰成分，所以即使到了 附近，φ值仍不会为零，因而上述计算λ的方法不能直接采用，可将上述计算的λ作为一个区间估计值，再用其他方法计算［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）方程，这时 就是 ；当λ→∞时，δ0→ ，而 是负梯度方向，这时 就是最速下降方向，所以阻尼最小二乘法的修正量 ，是最速下降修正量 和Gauss法修正量 之间的某种插值，λ就是这种插值的权系数。

Marguardt向量 具有以下三个特性：

（1）当λ越来越大时， 的长度越来越小，且

地球物理数据处理教程

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

（2）当λ由零逐渐增大时， 的方向逐渐由Gauss法的方向 转向最速下降法方向 ，λ越大， 方向越接近 方向。

（3）对λ＞0的任意正数， （满足方程（2.2.9））使φ在半径为‖ ‖的球面上取得极小。

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

以上三个性质说明，当λ逐渐增大时， 的方向由 向 靠近，它的大小‖ ‖逐渐减小，λ→∞时，‖ ‖→0，如图2.7所示。因此在迭代的任何一步，我们总可以找到充分大的λ，来保证稳定的收敛，因为当φ 不下降时，就加大λ向 靠，一直到使φ下降为止，从而保证收敛。性质（3）说明在跨出同样的步长时，以 （λ）方向最好，这就保证了该法的优越性。在实际计算时，总是在保证收敛的前提下，取较小的λ，以获得较大的步长前进。

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

（1）给出模型体参量初值 ，计算φ（ ）；给出实测场值ΔZk（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） 方程，重新得到新的 。

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

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

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

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

地电场与电法勘探

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

地电场与电法勘探

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

地电场与电法勘探

对上式作变量代换，令

地电场与电法勘探

并令

地电场与电法勘探

便得到

地电场与电法勘探

式中

地电场与电法勘探

从数字滤波理论可知，（2-1-88）的褶积运算为空间域数字滤波的表达式，如将J′1（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）是有限分布的，即它具有某一截止频率fz。根据采样定理，当采样间隔 Δx≤fz时，则 ρ′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）式可求得实际的电阻率转换函数，以下用TS（k）表示。用最优化法求层参数的大致过程是：先根据实际情况，给定一组层参数（叫做初值），算出T函数的理论值TL（k），将它与TS（k）比较，计算二者的误差，并根据它修改层参数；再算理论值，再作比较，再修改层参数，直到计算的TL（k）与TS（k）之差在规定范围内为止，便将这时理论值所对应的层参数作为解释结果。现稍为详细地讨论如下。

用P表示层数数ρ1、ρ2……ρn、h1、h2……hn-1，用一组P便可算出TL（k，P），理论值与实际值的相对误差可用δ表示：

地电场与电法勘探

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

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

地电场与电法勘探

从而可将δ写为

地电场与电法勘探

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

地电场与电法勘探 地电场与电法勘探

此式可改写为

地电场与电法勘探

上式为n个未知函ΔP1，ΔP2，……，ΔPn的n个线性联立方程，解此联立方程组，便可求得使δ到达极小点的最佳步长ΔP。但由于前面在作泰勒级数展开时，只取了一次项，即将T当作线性函数，然而T不是线性函数。故这样求得的ΔP，并不能使δ达到极小点，而只是向极小点靠近一步，达到了P′点（P′=P0+ΔP）。为了达到极小点，可再从P′点开始再用近似线性的办法形成新的线性方程组，又求得另一个步长ΔP（2），对参数作第二次修正，如此反复逐次迭代逼近，最后使TL（k，P）与TS（k）拟合得最好，并将这时的TL（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曲线，这样可兼顾速度与精度。