一維大地電磁Occam反演拉格朗日乘子的搜索

張君濤，周 軍，王緒本，夏時斌，鐘紅梅

（1.成都理工大學(xué)地球探測與信息技術(shù)教育部重點實驗室，成都 610059；2.四川省核工業(yè)地質(zhì)調(diào)查院，成都 610061）

一維大地電磁Occam反演拉格朗日乘子的搜索

張君濤1，周 軍1，王緒本1，夏時斌1，鐘紅梅2

（1.成都理工大學(xué)地球探測與信息技術(shù)教育部重點實驗室，成都 610059；2.四川省核工業(yè)地質(zhì)調(diào)查院，成都 610061）

在大地電磁反演中，Occam法因其在反演穩(wěn)定性和模型分辨率等方面的優(yōu)勢，得到廣泛應(yīng)用。但由于其每次迭代都需要不斷地搜索拉格朗日乘子，因而拉格朗日乘子的搜索效率對Occam法反演的運算速度起著至關(guān)重要的作用。為提高拉格朗日乘子的搜索效率，這里提出將拉格朗日乘子的搜索從自然數(shù)域中轉(zhuǎn)到以10為底的對數(shù)域中進行，同時在區(qū)間最小值的搜索中采用二次函數(shù)極小值搜索法。通過大量一維大地電磁反演驗證，該方法將拉格朗日乘子的搜索效率提高了20%～50%，大大提高了Occam反演的運算速度，方便了其在高維反演中的應(yīng)用。

大地電磁；Occam反演；拉格朗日乘子；對數(shù)；二次函數(shù)搜索

0 引言

在大地電磁反演中，利用目標函數(shù)梯度信息的最優(yōu)化方法占據(jù)著主要的地位［1－3，9］，包含Occam法［4－5］、非線性共軛梯度法（NLCG）［7］、擬牛頓法等反演方法。其中Occam法因其在反演穩(wěn)定性和模型分辨率等方面的優(yōu)勢，在現(xiàn)今的大地電磁反演中得到了廣泛地應(yīng)用，但由于該方法在每次反演迭代中，都要搜索一個滿足數(shù)據(jù)誤差最小或者模型粗糙度最大的拉格朗日乘子，因而需要進行多次正演，增加了反演的計算量，限制了其在高維反演的應(yīng)用。

很多學(xué)者在優(yōu)化Occam反演的計算量上做了大量的研究，Siripunvaraporn［6－7］等研究了數(shù)據(jù)空間Occam反演，并且成功的應(yīng)用到了MT二維、三維反演，極大的加強了Occam法在二維、三維反演中的應(yīng)用；吳小平［11］等在討論了Occam反演中數(shù)據(jù)擬合差隨拉格朗日乘子變化的基礎(chǔ)上，采用拉格朗日乘子在一定步長下逐次遞減的求取方法，每次迭代只需一次正演，極大地提高了計算速度，但初始拉格朗日乘子及減小比例的設(shè)定對反演穩(wěn)定性及結(jié)果影響較大。在對拉格朗日乘子搜索方式的改進上，也有很多人提出了改進方式［12－13］，但依然有很多地方值得探討。

基于上述考慮，作者對Occam法反演中提高拉格朗日乘子的搜索效率上做了進一步研究，結(jié)合拉格朗日乘子在Occam方法反演迭代中的變化特性，提出了對該參數(shù)應(yīng)當采用對數(shù)變換，并結(jié)合二次搜索方法來開展搜索計算。將上述方法應(yīng)用到了一維Occam反演的拉格朗日乘子搜索中，試圖提高拉格朗日乘子搜索的計算速度。這里首先介紹了Occam法反演的原理，然后詳細介紹了對數(shù)二次法搜索原理，并進行一維模型算例分析，說明了其快速搜索的優(yōu)勢。

1 Occam反演基本原理

一維Occam反演的目標函數(shù)包含數(shù)據(jù)目標函數(shù)和模型粗糙度函數(shù)：

其中：μ為拉格朗日乘子；大括號內(nèi)的為數(shù)據(jù)目標函數(shù)：d為實測數(shù)據(jù)；F為正演算子；W為數(shù)據(jù)誤差加權(quán)矩陣；X＊為擬合差期望值；‖?zm‖2為模型粗糙度函數(shù)。

反演的過程就是尋找使目標函數(shù)得到最小值的模型。為此Occam法將正演算子F線性化，在初始模型m0處按一階泰勒級數(shù)展開：

其中：J（m0）為正演算子F對模型參數(shù)m0的導(dǎo)數(shù)，稱為雅可比矩陣。將式（2）帶入式（1）中，由此得到的目標函數(shù)已經(jīng)簡化為一個關(guān)于模型參數(shù)m的二次函數(shù)：

其中：d0＝d－F（m0）＋Jm0，J＝J（m0），再求式（3）關(guān)于m的一階導(dǎo)數(shù)，得到梯度函數(shù)g（m）：

此時的梯度函數(shù)是一個關(guān)于模型參數(shù)的一次函數(shù)，令梯度g（m）＝0，求得的模型參數(shù)即為目標函數(shù)取得極小值的模型：

由于正演算子是非線性的，目標函數(shù)也不是二次函數(shù)，故反演過程中需要不斷迭代，反復(fù)修正初始模型m0，重復(fù)式（5）過程，最終求得目標函數(shù)的最小值。

在每次迭代中，對于不同的拉格朗日乘子u，利用式（5）可得到不同的模型參數(shù)及其相應(yīng)的絕對擬合誤差。相比其他反演方法直接輸入一個經(jīng)驗的拉格朗日乘子，Occam法反演提出在每次迭代中都要通過搜索，得到一個滿足數(shù)據(jù)擬合誤差最小的μ值（如果數(shù)據(jù)擬合誤差低于目標擬合誤差，則搜索滿足目標數(shù)據(jù)擬合差的最小模型粗糙度的μ值）。用此μ值求得模型參數(shù)m再作為下一次迭代的初始模型，進行第二次迭代。這樣反復(fù)迭代，最后找到一個滿足目標擬合差且粗糙度最小的模型，或者是達到最大迭代次數(shù)后得到的擬合差最小的模型。正是由于Occam法反演對拉格朗日乘子的搜索，使其具有很高的反演穩(wěn)定性和模型分辨率，但同時也增加了反演的運算量。

2 拉格朗日乘子的搜索方法

考慮到拉格朗日乘子的搜索效率對Occam法反演運算時間的巨大影響，作者對拉格朗日乘子的搜索方式進行改進。

所有患者采用隨機數(shù)字表法分為舌下含服組和口服組,每組各40例?？诜M患者口服酒石酸美托洛爾片(阿斯利康制藥有限公司,產(chǎn)品批號1407184)25 mg(1片);舌下含服組患者舌下含服酒石酸美托洛爾片25 mg(1片)。

對于μ值的搜索，一般采用先搜索到一個谷值區(qū)間，然后在此區(qū)間上搜索數(shù)據(jù)擬合差的最小值，搜索到最小值后進行下一次迭代；而在搜索過程中，一旦有數(shù)據(jù)擬合誤差小于目標擬合誤差就直接進入滿足擬合差的最小模型粗糙度的μ值搜索。這里也采用以上搜索策略，但做了兩點改進。

2.1 拉格朗日乘子對數(shù)化搜索

理論上μ值可以為大于“0”的任意數(shù)，但如果直接對μ進行搜索，在區(qū)間搜索過程中，μ值的修正步長是隨著μ值的變化而變化的，且初始μ值如果給的不合適，區(qū)間搜索的次數(shù)會相對較大。考慮到將原目標函數(shù)中的拉格朗日乘子μ變?yōu)?0λ后，并不會改變其原有目標函數(shù)的極值特性，只要滿足μ＝10λ其對應(yīng)的模型參數(shù)及其對應(yīng)的目標函數(shù)與數(shù)據(jù)擬合差都是一致的，因而作者大膽設(shè)想將μ值做對數(shù)變化，從而將目標函數(shù)中對μ值的搜索轉(zhuǎn)為對10λ中的λ做搜索。將μ＝10λ代入到目標函數(shù)式（1）中得式（6）。

式（6）經(jīng)過式（3）、式（4）、式（5）運算過程可得每次迭代中模型參數(shù)的求取公式：

可以看出式（7）只是將式（5）中的μ變?yōu)?0λ，從而將對μ的搜索轉(zhuǎn)為對λ的搜索，這將更加快速地搜索到谷值區(qū)間，并減小初始λ值對搜索效率的影響。

從圖1可以看出，圖1（a）絕對擬合誤差曲線整體上變化平緩，更近似一個二次曲線；如果我們將區(qū)間搜索步長在以10為底的對數(shù)上設(shè)為“1”，一般在左右?guī)讉€步長之內(nèi)就能搜索到谷值區(qū)間，從而擺脫了初始值的設(shè)定對搜索次數(shù)的影響。圖1（b）絕對擬合誤差曲線在最小值左右變化相當劇烈，而在μ值較大時變化相當平緩。這給區(qū)間搜索步長的確定帶來很大的不便，且初始μ值對區(qū)間搜索的次數(shù)有很大的影響。正是基于上述因素，這里更加明確將拉格朗日乘子的搜索從自然數(shù)域中轉(zhuǎn)到了以10為底的對數(shù)域進行，將拉格朗日乘子變化為μ＝10λ，從而將對μ值的搜索改為對λ的搜索，其實質(zhì)就是在以10為底的對數(shù)上進行μ值搜索。為便于進行谷值區(qū)間最小值搜索，將搜索到的谷值區(qū)間三點按λ值從大到小保存為（λ0，rms0）、（λ1，rms1）、（λ2，rms2）。

圖1 模型一首次迭代絕對擬合誤差曲線Fig.1 The absolute error curve of the first iteration of the model 1 （a）拉格朗日乘子以10為底；（b）拉格朗日乘子用自然數(shù)

2.2 二次函數(shù)搜索區(qū)間極小λ值

做了上述變化后能快速地搜索到一個谷值區(qū)間，但此區(qū)間的范圍相對自然數(shù)搜索的區(qū)間較大，如果區(qū)間最小值搜索繼續(xù)使用二分法進行，則搜索次數(shù)相對自然數(shù)搜索的谷值區(qū)間二分法，增加了很多次。同樣以圖2（a）為例，首次迭代對數(shù)搜索區(qū)間次數(shù)為4次，但用二分法進行區(qū)間最小值搜索則用了8次，共計12次；自然數(shù)區(qū)間搜索用了9次，但最小值搜索只用了5次，共14次；二者相比，對數(shù)搜索基本沒有優(yōu)勢。因此，提出利用二次函數(shù)極小值來搜索區(qū)間最小值。

構(gòu)建二次函數(shù)：

本文對于退出搜索需要滿足以下要求：

其中ε≤1。當滿足式（9）要求后，則取擬合差小的λ值及其相對應(yīng)的模型進入下一次迭代。

3 算例分析

為了驗證本方法的搜索效率，建立了四種地電模型（圖2）。特別說明模型一為陳小斌［10］在自適應(yīng)正則化反演中的八層驗證模型，這里利用此模型將本方法與傳統(tǒng)拉格朗日乘子二分法的搜索過程進行了詳細地比較。四個模型的反演目標擬合差均為0.001。為便于敘述，將作者所描述的搜索方法稱為“對數(shù)二次法搜索”，將傳統(tǒng)拉格朗日乘子二分法的搜索稱為“自然二分法搜索”。

表1為八層模型在首次迭代中的μ值變化過程，對數(shù)二分法只經(jīng)過4次區(qū)間搜索加3次最小值搜索就完成了迭代，而自然二分法經(jīng)歷了9次區(qū)間搜索加5次最小值搜索才完成迭代。而從搜索的結(jié)果來看，兩種方法搜索到的μ值基本上沒有差異，甚至本例中對數(shù)二次法搜索得到的絕對擬合誤差還更小一點?？梢妼?shù)二分法在搜索上的優(yōu)勢，這種優(yōu)勢在首次迭代中體現(xiàn)得最為明顯。

圖2 地電模型及Occam反演結(jié)果Fig.2 The geoelectric model and Occam inversion results （a）模型1；（b）模型2；（c）模型3；（d）模型4

表1 模型1首次迭代μ值搜索過程Tab.1 The searching process of Lagrange multiplierin first iteration of model 1

表2、表3分別為模型一Occam反演用對數(shù)二次法搜索與自然二分法搜索的詳細迭代過程。結(jié)果可以看出，對數(shù)二次法區(qū)間搜索往往只需要3次～4次，而3次搜索已經(jīng)是區(qū)間搜索的極限；二次最小值搜索基本上也只需要1次～3次。整體對比，每次迭代對數(shù)二次法搜索的次數(shù)都明顯少于自然二分法搜索。

表4為四種模型分別用兩方搜索法進行Occam反演的μ值變化次數(shù)對比。無一例外，對數(shù)二次法搜索次數(shù)均少于自然二分法搜索，減小比例在20% ～50%之間。

從圖2反演結(jié)果上看，本方法在提高μ值搜索效率的同時，對Occam反演的結(jié)果沒有任何影響。

表2 模型1每次迭代對數(shù)二次法搜索次數(shù)Tab.2 The search process of each iteration of model 1by logarithmic quadratic method

表3 模型1每次迭代自然二分法搜索次數(shù)Tab.3 The search process of each iteration of model 1by natural dichotomymethod

表4 四種模型μ值變化總次數(shù)對比Tab.4 The comparison of the number of Lagrange multiplier's search

4 結(jié)論

作者提出了一種Occam反演中拉格朗日乘子的搜索方法，此方法主要做了以下兩方面的改進。

1）利用關(guān)于拉格朗日乘子的絕對擬合誤差曲線的形態(tài)特征，用公式μ＝10λ對拉格朗日乘子做個變化，從而將對μ值的搜索變?yōu)閷Ζ酥档乃阉?，即將拉格朗日乘子的搜索從自然?shù)域中轉(zhuǎn)化到了以10為底的對數(shù)域中進行。在運用過程中，發(fā)現(xiàn)對λ進行區(qū)間搜索的確很快，但也變相地加大了區(qū)間最小值的搜索范圍，在實際算例中，整體的搜索效率并沒有提高多少，故而進行了第二次改進。

2）利用區(qū)間搜索到的已知三點，構(gòu)建關(guān)于λ值的擬合差二次函數(shù)，并求得此二次函數(shù)取得極小值時的λ值；用此λ值代入反演中求得擬合差，再得到一個新的谷值區(qū)間。這樣反復(fù)進行二次函數(shù)構(gòu)建，尋找區(qū)間最小值。在實際算例中，我們發(fā)現(xiàn)往往只需要1次～3次二次函數(shù)構(gòu)建就能搜索到區(qū)間的最小值，這大大地提高了區(qū)間最小值的搜索效率。

將以上兩種改進運用到拉格朗日乘子的搜索中，用模型算例驗證發(fā)現(xiàn)，此方法將搜索效率提高了20%～50%，大大減少了Occam法反演的運算時間。

本搜索方法僅僅在一維大地電磁反演中進行了驗證，有待在二維、三維中運用。

［1］ 石應(yīng)駿，劉國棟，吳廣耀，等.大地電磁測深法教程［M］.北京：地展出版社，1985.

SHI Y J，LIU G D，WU G Y，Tutorial of MT sounding［M］.Beijing：Geological Publishing House，1985.（In Chinese）

［2］ 陳樂壽，劉任，王天生.大地電磁測深資料處理與解釋［M］.北京：石油工業(yè)出版社，1989.

CHEN L S，LIU R，WANG T S.Data processing and interpretation of MT sounding［M］.Beijing：Petroleum Industry Press，1989.（In Chinese）

［3］ 王家映.地球物理反演理論［M］.北京：高等教育出版社，2002.

WANG J Y.Inverse theory in geophysics［M］.Beijing：Higher Education Press，2002.（In Chinese）

［4］ CONSTABLE S C，PARKER R L，CONSTABLE C G.Occam's inversion：A practical algorithm for generating smooth models from electromagnetic sounding data［J］.Geophysics，1987，52（3）：289－300.

［5］ DEGROOT－HEDLIN C D，CONSTABLE S C.Occam's Inversion to Generate Smooth Two Dimensional Models form MagnetotelluricData［J］.Geophysics，1990，55（12）：1613－1624.

［6］ SIRIPUNVARAPORN W，EGBERTG.An efficient data－subspaceinversion method for 2－D magnetotelluricdata［J］.Geophysics，2000，65（3）：791－803.

［7］ SIRIPUNVARAPORN W，UYESHIMA M，EGBERT G.Three dimensional inversion for network－magnetotelluricdata［J］.Earth Planets Space，2004，56（9）：893－902.

［8］ RODI W L，MACKIE R L.Nonlinear conjugate gradients algorithm for 2－D magnetotelluricinversion［J］.Geophysics，2001，66（3）：174－187.

［9］ 韓波，胡祥云，何展祥，等.大地電磁反演方法的數(shù)學(xué)分類［J］.石油地球物理勘探，2012，47（1）：177－187.

HAN B，HU X Y，HE ZH X，et al.Mathematical classification of magnetotelluric inversion methods［J］.Oil Geophysical Prospecting，2012，47（1）：177－187.（In Chinese）

［10］陳小斌，趙國澤，湯吉，等.大地電磁自適應(yīng)正則化反演算法［J］.地球物理學(xué)報，2005，48（4）：937－946.

CHEN X B，ZHAO G Z，TANG J，et al.An adaptive regularized inversion algorithm for magnetotelluric data ［J］.Chinese J Geophys，2005，48（4）：937－946.（In Chinese）

［11］吳小平，徐果明.大地電磁數(shù)據(jù)的OCCAM反演改進［J］.地球物理學(xué)報，1998，41（4）：547－554.

WU X P，XU G M.Improvement of OCCAM＇s inversion for MT date［J］.Chinese J Geophys，1998，41 （4）：547－554.（In Chinese）

［12］樸英哲，李桐林，劉永亮.在大地電磁二維Occam反演中求取拉格朗日乘子方法改進［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，2014，44（2）：660－667.

PAK Y C，LI T L，LIU Y L.Improvement of choosing lagrange multiplier on MT 2DOccam inversion ［J］.Journal of Jilin University：Earth Science Edition，2014，44（2）：660－667.（In Chinese）

［13］劉俊峰，鄧居智，陳輝，等.一種用于Occam反演中搜索拉格朗日乘子的方法［J］.工程地球物理學(xué)報，2013，10（3）：344－350.

LIU J F，DENG J Z，CHEN H，et al.A method used for searching lagrange multiplier in Occam inversion ［J］.Chinese journal of engineering geophysics，2013，10（3）：344－350.（In Chinese）

Lagrange multiplier's search in one－dimensional magnetotelluric Occam's inversion

ZHANG Jun－tao1，ZHOU Jun1，WANG Xu－ben1，XIA Shi－bin1，ZHONG Hong－mei2
（1.Key Lab of Earth Exploration ＆information Techniques of Ministry of Education，Chengdu university of technology，Chengdu 610059，China；2.Sichuan Nuclear Industry Geological Survey Institute，Chengdu 610059，China）

In magnetotelluric inversion，because of its stability in inversion and the advantage of model resolution and so on，Occam method is widely used.But in each iteration，Occam＇s inversion needs to constantly search Lagrange multiplier.The searching efficiency of Lagrange multiplier to Occam＇s inversion speed plays a vital role.To improve the search efficiency of Lagrange multiplier，a new searching method was presented in this paper.First Lagrange multiplier of transfer from the natural number field to search for the 10logs base domain.Second in the search of interval minimum search method using quadratic function.Through a number of one－dimensional magnetotelluric inversions，Lagrange multiplier search efficiency is increased by 20%－50%，and the operation speed of Occam＇s inversion is greatly improved.

magnetotelluric；Occam inversion；Lagrange multiplier；logarithmic；quadratic function

P 631.3

：A

10.3969／j.issn.1001－1749.2015.06.03

1001－1749（2015）06－0687－06

2014－11－18改回日期：2014－12－14

地質(zhì)調(diào)查項目（12120113095100）；重大儀器專項（2011YQ05006007）；四川省科技支撐計劃項目（2013FZ0084）；四川省教育廳項目（13ZB0062）

張君濤（1989－），男，碩士，主要研究方向為電磁法數(shù)據(jù)處理與正反演解釋，E－mail：Zjuntao＠foxmail.com。