最小阻尼二乘演算法

㈠ 非線性回歸方程為什麼不過樣本中心

關於回歸方程，存在線性回歸方程和非線性回歸方程，其目的在於最大程度的擬合一組數據，數據由自變數和因變數組成。當判斷線性回歸方程的相關系數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）。