基于最小支持的2.5維直流電阻率正則化反演研究

李 曼， 林文東，2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中eps 是計(jì)算機(jī)浮點(diǎn)相對精度，在其文章中，β 按如下公式取值:

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

2 L 曲線法則及修正

2.1 Hansen 的L 曲線法

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

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

圖1 L 曲線Fig.1 L-curve

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

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

對于一般的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。從而可以得出:

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

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

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

2.2 修正的L 曲線法

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

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

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

針對這一情況，對L 曲線準(zhǔn)則進(jìn)行了改進(jìn)，原理基于點(diǎn)與直線的關(guān)系以及點(diǎn)到直線的距離，具體算法如下:

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

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

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

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

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

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

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

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

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

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

3 算例分析

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

姜祝輝，黃思訓(xùn)，何然，等. 2011. 合成孔徑雷達(dá)資料反演海面風(fēng)場的正則化方法研究［J］.物理學(xué)報(bào)，60(6):1-8.

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

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

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

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

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

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

趙烈加，劉衛(wèi)，陳小華，等. 2007.基于L 曲線準(zhǔn)則的T2 譜反演算法研究［J］.石油天然氣學(xué)報(bào)，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.