李 曼， 林文東，2

(1.東華理工大學核工程與地球物理學院，江西 南昌330013;2.遼寧省有色地質(zhì)局103 隊，遼寧 丹東118008)

地球物理解釋工作的目的是使解釋結果逼近真實地質(zhì)模型，正則化是得到穩(wěn)定解的重要手段。正則化方法中穩(wěn)定因子的設計與正則化因子的選擇是兩項重要研究內(nèi)容。

Tikhonov 等(1977)提出了正則化方法

式中，P(α，m)是總的目標函數(shù)，α 是正則化因子，φd(m)是數(shù)據(jù)目標函數(shù)，φm(m)是模型約束目標函數(shù)，也稱穩(wěn)定因子。

穩(wěn)定因子的主要功能是對模型解的空間進行限制，以減少多解性，求得穩(wěn)定解(Zhdanov，2002)，穩(wěn)定因子的設計是正則化方法中的重要研究內(nèi)容。在地球物理學中，可以使用的穩(wěn)定因子有很多，根據(jù)實際問題的不同，可以選擇的穩(wěn)定因子也有不同。本文為實現(xiàn)陡變邊界反演，研究了最小支持穩(wěn)定因子。

最小支持穩(wěn)定因子最先是由Last 等(1983)提出并運用于重力數(shù)據(jù)的反演解釋;Mehanee 等(2002)在進行了最小支持穩(wěn)定因子提高塊狀結構分辨率的研究;Candansyar(2002)在其二維大地電磁的反演的文章中進行了不同穩(wěn)定因子的比較;在一篇關于地震應用的文獻中，Portniaguine 等(2004)對聚焦參數(shù)的取值進行了分析，提出了聚焦參數(shù)依據(jù)計算機字長的選擇方法;Zhdanov 等(2004)在3D 重力張量的正則化反演中使用了該穩(wěn)定因子，使反演結果更加接近真實值;劉小軍等(2007)在二維大地電磁數(shù)據(jù)的聚焦反演算法中運用了類似的正則化技術，取得了較好的反演效果;Vignoli 等(2012)在一維條件下對最小支持穩(wěn)定因子中的聚焦參數(shù)的選擇取進行了討論。

要想得到穩(wěn)定解，就需要使得方程的右邊兩項達到較好的平衡，因此選擇恰當?shù)恼齽t化因子同樣是十分必要的。早期地球物理學學者采用一些經(jīng)驗算法確定，隨著越來越多的學者對這一方面關注，提出大量的正則化因子選取方法。其中，以Hansen 所提出的L 曲線法最為著名。該方法是基于數(shù)據(jù)誤差水平未知的啟發(fā)式后驗選取方法，在雙對數(shù)坐標下形成單調(diào)遞減的曲線，類似于字母“L”，故稱為L 曲線法。Hansen(1992)首次提出L 曲線準則方法;此后幾年，Hansen 等(1993)對該法進行了系統(tǒng)的研究，完善了L 曲線法的理論;Hansen 等(1999)給出了較實用的L 曲線法計算技術。

在地球物理方面，Li 等(1999)首次在他的文章中將“L 曲線法”引入地球物理反演中;Farquharson 等(2000)在他們的文章中對廣義交叉原理(GCV)和“L 曲線法”在非線性反演問題中選擇正則化因子的對比進行了論述;近年來，該方法得到廣泛關注，王振杰等(2004，2005)研究了用L 曲線法確定半?yún)?shù)模型中的平滑因子和一種解算病態(tài)問題的兩步解法;趙烈加等(2007)研究了基于L 曲線準則的譜反演算法;姜祝輝等(2011)在合成孔徑雷達資料反演海面風場中也使用了正則化方法;于勝杰等(2012)研究了選權擬合法進行GPS 水汽層析解算;陶肖靜等(2012)研究了半?yún)?shù)模型中影響正則化因子的因數(shù);李倫等(2013)研究了基于正則化方法的高頻地波雷達海浪方向譜反演，這些研究中均使用了L 曲線準則確定正則化因子。

本文通過簡單的點與直線關系對“L 曲線法”進行了改進;該方法不需要Hansen 經(jīng)典算法計算曲線曲率，計算簡間、高效，同時可以有效的去除“假”正則化因子，在求解實際地球物理問題中應用效果顯著。文章使用修正的L 曲線法自動選取正則化參數(shù)，采用最小支持穩(wěn)定因子作為模型約束目標函數(shù)，同時給出了形式上的模型參數(shù)表示方法，采用正則化的共軛梯度法對目標函數(shù)極小化。將最小支持穩(wěn)定因子應用到算例分析中，對理論模型進行試算，同時對聚焦參數(shù)的取值進行了討論，對比分析后，證明了本文所采用的最小支持穩(wěn)定因子表示方法及參數(shù)選擇算法的正確高效，最小支持穩(wěn)定因子具有更好的聚焦能力、有利于實現(xiàn)陡變邊界反演。

1 最小支持穩(wěn)定因子及聚焦參數(shù)的選取

穩(wěn)定因子的主要功能是對模型解的空間進行限制，以減少多解性，求得穩(wěn)定解。地球物理學中常用的穩(wěn)定因子主要有最小范數(shù)穩(wěn)定因子、最大平滑穩(wěn)定因子等，但這兩者都是基于模型光滑，當模型出現(xiàn)突變時，這兩種穩(wěn)定因子計算結果仍是漸變的，不利于地質(zhì)資料解釋。最小支持穩(wěn)定因子:

式中，β 為聚焦參數(shù)。

Zhdanov(2002)給出了穩(wěn)定因子的偽二次統(tǒng)一表示:

其中函數(shù)we是關于m 的一個非線性函數(shù)，在這種情況下，由式(3)所確定的函數(shù)s(m)不是二次的，這就是稱其為“偽二次”函數(shù)的原因。

用穩(wěn)定因子的“偽二次”式(3)，可以把對應的Tiknonov 正則化方程(1)變?yōu)?

最小支持穩(wěn)定因子在正則化問題中的一個簡化的“假的二次”形式，可以使用如下表示方法以達到在地球物理問題中得到近似正則化解的目的。

對于最小支持穩(wěn)定因子SMS(m)來說，假定若mapr= 0，則

根據(jù)最小支持穩(wěn)定因子的性質(zhì)，在反演中模型會盡可能的集中，從而可以區(qū)分物性差異較大的界面。對于最小支持穩(wěn)定因子來說，聚焦參數(shù)β 的取值也對反演結果具有很大的影響。Portniaguine 等(2004)，Vignoli 等(2012)的研究表明聚焦參數(shù)的取值與計算機的精度有關，下面介紹Vignoli 的選擇方法。當m = mapr時，為避免奇點β 的取值要滿足如下:

其中eps 是計算機浮點相對精度，在其文章中，β 按如下公式取值:

其中max(m - mapr)是m - mapr向量中的最大元素。

2 L 曲線法則及修正

2.1 Hansen 的L 曲線法

正則化因子作為數(shù)據(jù)目標函數(shù)和模型約束目標函數(shù)的權重系數(shù)，若正則化因子過大則產(chǎn)生的結果主要擬合先驗模型，若正則化因子過小則產(chǎn)生的結果主要擬合觀測數(shù)據(jù)。若在這兩種情況之間取值，就可以得到折中解，該解在某種程度上不僅可以滿足精度要求，也可以保證解的穩(wěn)定性。

由Hansen 等提出的L 曲線技術是一種基于數(shù)據(jù)誤差水平未知的啟發(fā)式后驗選取正則化因子的方法，該方法以log-log 尺度來描述φd(m)與φm(m)，其正則化因子α 的取值為雙對數(shù)曲線上曲率的最大值，一般是其“隅角”所對應的位置，如圖1中α0所示位置。

圖1 L 曲線Fig.1 L-curve

在Hansen 的經(jīng)典論文中給出了通過曲線的曲率半徑計算“隅角”的方法。地球物理的反演問題是已知觀測數(shù)據(jù)求與其對應模型的過程。設d 為觀測數(shù)據(jù)向量，m 是模型參數(shù)向量，A 為把地球物理模型映射到理想數(shù)據(jù)的函數(shù)，正問題可表示為，

式中A 為正演算符，Tikhonov 正則化的基本思想是將解的范數(shù)作為先驗信息考慮，則問題的求解(假定正問題是線性的)將轉化為求解以下問題的最小值:

對于一般的m × n 矩陣A，一定存在正交陣U= (u1，u2，Λ，un)∈Rm×m和V = (v1，v2，…，Vn)∈Rn×n，可以使得

其中σi表示奇異值，且σ1≥σ2≥…≥σp≥0 是遞減的順序排列，∑= diag(σ1，σ2，…，σp)。ui為左奇異向量，相互正交，即有UTU = Im，vi為右奇異向量，相互正交，即有VTV = In。從而可以得出:

為了運算簡便，令正則化參數(shù)α = λ2。從而方程組A(m)= d 的Tikhonov 正則化所對應的解可以寫為:

對于Tikhonov 正則化方程組來說，應用奇異值分解可以得到:

求κ(η)的最大值便可以得到合適的正則化因子。

2.2 修正的L 曲線法

對于Hansen 的L 曲線準則來說，若知道L 曲線的參數(shù)表達式，即函數(shù)η(α)和ρ(α)的精確表達式，采用曲率公式求曲率的最大值可以得到合適的正則化因子。但在實際應用中，很難得到精確的η(a)和ρ(α)表達式，此時可以采用三次樣條插值的方法來近似的確定L 曲線的“隅角”，從而選取適當正則化因子，其主要思想為:將節(jié)點個數(shù)逐步增加，形成的參數(shù)樣條可以更好的逼近L 曲線，選取的最優(yōu)正則化因子為最終樣條上曲率最大的點。

曲率公式(15)中含有兩個一階導數(shù)和兩個二階導數(shù)，實際應用中常常會導致計算結果不穩(wěn)定，Hansen 對其進行了簡化，得到了只含有一個一階導數(shù)的公式，同時可以證明L 曲線是單調(diào)遞減的曲線，公式在理論上可以很好的確定正則化因子，但在實際應用中，L 曲線可能出現(xiàn)兩個“拐點”的情況，如圖2 所示。

圖2 實際應用中可能出現(xiàn)的L 曲線形狀Fig.2 The shape of the L curve may occur in practical application

針對這一情況，對L 曲線準則進行了改進，原理基于點與直線的關系以及點到直線的距離，具體算法如下:

(1)給定初始的正則化因子αi(至少四個，水平部分兩個，豎直部分兩個)，及其對應的(pi，qi)

(2)關于βi= logαi對pi和qi分別做三次樣條插值，分別將插值函數(shù)記為p(β)和q(β);

(3)找到L 曲線即曲線(p(β)，q(β))的兩個端點，確定一條直線，寫出直線表達式Ax + By + C= 0;

(4)在L 曲線上任選一點，(x0，y0)，將其帶入直線表達式，若Ax0+By0+C ＜0，則保留，若Ax0+By0+ C ≥0，則舍棄，直到所有L 曲線上的點都選完;

(5)從選取的點中任取一點，(x1，y1)，帶入點到直線公式計算出每一點到直線的距離，計算距離最大值的點，設其為(p(β0)，q(β0))，記α0= 10β0;

(6)求解α0所對應的正則化問題，記正則化解為mα0，由正則化解可以得到新的點對(p0，q0);

(7)將新的點對(p0，q0)和α0加入到步驟(1)的點對和參數(shù)中;

(8)將步驟(7)中得到的點對和參數(shù)重復步驟(2)到步驟(7)的過程，直至收斂。

相對于一般的L 曲線準則來說，修正算法增加了步驟(3)，(4)和(6)，去掉了“假的”正則化因子，可以使選擇正則化因子更符合實際地球物理問題，為后續(xù)的研究工作打下堅實基礎。

圖3 組合模型示意圖Fig 3 Diagram of multiple model

3 算例分析

在ρ0=100 Ω·m 背景電阻率條件下，存在兩個截面為正方形(邊長4 m)的水平棱柱，其中高阻體，低阻體ρ2=50 Ω·m，相距2 m，異常體頂部埋深2 m，如圖3 所示。地表采用30個電極，電極距1 m，采用溫納-裝置進行高密度電阻率測量。

采用正演計算結果(加入0.5%的白噪聲)作為反演數(shù)據(jù)集，運用共軛梯度法、修正的L 曲線法算法，分別最小范數(shù)穩(wěn)定因子、最大平滑穩(wěn)定因子及最小支持穩(wěn)定因子對進行反演試算。反演結果(圖4，5，6，7)。

圖4 穩(wěn)定因子為最小范數(shù)穩(wěn)定因子時的反演結果Fig.4 The inversion results of the minimum norm stabilizing factor

圖5 穩(wěn)定因子為最大平滑穩(wěn)定因子(一階)時的反演結果Fig.5 The inversion results of the maximum smoothness stabilizing factor (first-order derivative)

從圖4 ～圖7 中可以看出，反演結果與實際的地電模型吻合程度很好。從反演斷面圖可以看出，在使用修正的L 曲線法自動選取正則化因子的情況下，采用不同的穩(wěn)定因子，得到的反演結果都能夠較好的反應真實的地下地質(zhì)結構。在同樣地誤差精度和迭代次數(shù)情況下，(1)最小范數(shù)穩(wěn)定因子和最大平滑穩(wěn)定因子(一階、二階)的反演斷面圖得出，反演結果能與異常體吻合，但邊界為平滑過度。(2)相比于最小范數(shù)穩(wěn)定因子和最大平滑穩(wěn)定因子，最小支持穩(wěn)定因子的反演結果能夠更好的反應陡變異常體邊界，異常體深度對應優(yōu)于最小范數(shù)與最大平滑穩(wěn)定因子。

圖6 穩(wěn)定因子為最大平滑穩(wěn)定因子(二階)時的反演結果Fig.6 The inversion results of the maximum smoothness stabilizing factor (second-order derivative)

圖7 穩(wěn)定因子為最小支持穩(wěn)定因子時的反演結果Fig.7 The inversion results of the minimum support stabilizing factor

4 結論

地球物理反問題存在多解性和不穩(wěn)定性，正則化是得到穩(wěn)定解的重要手段之一。正則化方法中穩(wěn)定因子的設計與正則化因子的選擇是兩項重要研究內(nèi)容。本文開展了L 曲線法自動選取正則化因子的修正算法及最小支持穩(wěn)定因子的研究工作，應用非線性共軛梯度法對2.5 維直流電阻率反問題進行了試算，結果表明:

(1)針對實際地球物理問題，運用了點與直線關系和點與直線距離對L 曲線法進行了修改，提出了一種修正的L 曲線法，該方法計算簡單、高效，有效了避免選擇過程落入“假”正則化因子區(qū)間。

(2)最小支持穩(wěn)定因子相對于最小范數(shù)穩(wěn)定因子和最大平滑穩(wěn)定因子具有更好的識別陡變異常體邊界的能力。

姜祝輝，黃思訓，何然，等. 2011. 合成孔徑雷達資料反演海面風場的正則化方法研究［J］.物理學報，60(6):1-8.

李倫，吳雄斌，龍超，等. 2013. 基于正則化方法的高頻地波雷達海浪方向譜反演［J］.地球物理學報，56(1):219-229.

劉小軍，王家林，陳冰，等.2007.二維大地電磁數(shù)據(jù)的聚焦反演算法探討［J］.石油地球物理勘探，42(3):338-342.

陶肖靜，朱建軍，田玉淼. 2012.半?yún)?shù)模型中影響正則化因子的因素分析［J］. 武漢大學學報·信息科學版，37(3):298-301.

王振杰，歐吉坤，曲國慶.2004. 用L 曲線法確定半?yún)?shù)模型中的平滑因子［J］. 武漢大學學報·信息科學版，29(7):651-653.

王振杰，歐吉坤，柳林濤. 2005.一種解算病態(tài)問題的方法—兩步解法［J］.武漢大學學報·信息科學版，30(9):821-824.

于勝杰，柳林濤. 2012. 利用選權擬合法進行GPS 水汽層析解算［J］.武漢大學學報·信息科學版，37(2):183-186.

趙烈加，劉衛(wèi)，陳小華，等. 2007.基于L 曲線準則的T2 譜反演算法研究［J］.石油天然氣學報，29(6):82-86.

Candansayar M E. 2002. Two-dimensional inversion of magnetotelluric data with CG and SVD:A comparisonal study of different stabilizers［J］. SEG Technical Program Expanded Abstracts，21(1):669-672.

Farquharson C G，Oldenburg D W. 2000. Automatic estimation of the trade-off parameter in nonlinear inverse problems using the GCV and L-curve criteria［J］. SEG Expanded Abstracts，19(1):265-268.

Hansen P C. 1992. Analysis of discrete ill-posed problems by means of the L-curve［J］.Siam Review，34(4):561-580.

Hansen P C，O’Leary D P. 1993. The use of the L-curve in the regularization of discrete ill-posed problems［J］. SIAM J. Sci. Comput，14:1487-1503.

Hansen P C. 1998. Rank-Deficient and Discrete Ill-Posed Problems［M］. Philadelphia:SIAM.

Hansen P C. 1999. The L-curve and its use in the numerical treatment of inverse problems［M］. Denmark:Department of Mathematical Modeling.

Last B J，Kubik K. 1983. Compact gravity inversion［J］. Geophysics，48(6):713-721.

Li Y G，Oldenburg D W. 1999. 3-D inversion of DC resistivity data using an L-curve criterion［J］.SEG Expanded Abstracts，18(1):1-4.

Mehanee S，Zhdanov M S. 2002. Two-dimensional magnetotelluric inversion of blocky geoelectrical structures［J］. Journal of geophysical research，107(B4):2065.

Portniaguine O，Castagna J P. 2004. Inverse spectral decomposition［J］.SEG Expanded Abstracts，23(1):1786-1789.

Tikhonov A N，Arsenin V Y. 1977. Solution of Ill-posed Problems［M］.Washing.D.C:V.H.Winston and Sons.

Vignoli G，Deiana R，Cassiani G. 2012. Focused inversion of vertical radar profile (VRP)traveltime data［J］. Geophysics，77(1):9-18.

Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems［M］. Netherlands:Elsevier.

Zhdanov M S，Ellis R，Mukherjee S. 2004. Three-dimensional regularized focusing inversion of gravity gradient tensor component data［J］.Geophysics，69(4):925-937.