大地電磁反演中改進的自適應(yīng)正則化因子選取

向 陽，于 鵬，陳 曉，唐 睿

（同濟大學(xué) 海洋地質(zhì)國家重點實驗室，上海200092）

如何解決反演過程受初始模型影響的局限并提高解的穩(wěn)定性，一直是地球物理反演的關(guān)鍵問題．吉洪諾夫正則化方法［1］已被廣泛接受和應(yīng)用于解決這種不適定性反問題．正則化反演通過在目標函數(shù)中加入模型約束，進而降低反演的多解性．正則化因子是數(shù)據(jù)擬合泛函和模型穩(wěn)定泛函之間的折衷系數(shù)，正則化因子過大表示強調(diào)模型約束，數(shù)據(jù)總體是欠擬合的［2］；正則化因子過小表示反演主要擬合數(shù)據(jù)，反演結(jié)果容易受到數(shù)據(jù)噪音的影響產(chǎn)生虛假構(gòu)造，總體是過擬合的，因此正則化因子的選取直接影響到反演結(jié)果［3］．如何選擇合適的正則化因子，使反演既滿足觀測數(shù)據(jù)的充分擬合，同時還能保證模型的穩(wěn)定性，一直是正則化反演中研究的熱點，許多學(xué)者對此進行了深入研究．

前人研究較多的是經(jīng)驗定值方法．其中包括無偏風(fēng)險估計預(yù)測方法［4－5］，該方法能給出估計值和真值之間的無偏的誤差估計，但需要已知數(shù)據(jù)的噪音方差，而廣義交叉驗證方法［6－9］，則不需要估計數(shù)據(jù)的噪音方差便可進行計算，因此應(yīng)用更為廣泛．Hansen等［10－11］提出了L曲線法，選取L曲線拐點對應(yīng)的值為最佳正則化因子，但需要進行多次試算才能得到一個最佳選擇，且這個最佳也不是一個嚴格的取值，是一個相對的最佳［12］．

為了避免上述方法計算費時的問題，出現(xiàn)了許多正則化因子隨反演過程可自適應(yīng)調(diào)整的算法．在Occam 反演中［13－14］選取逼近數(shù)據(jù)擬合期望的最大的正則化因子為每次迭代的選擇，該方法在每次反演計算中需額外正演多次，運算量較大．吳小平等［15］提出一種給定正則化因子初值，每次迭代直接按照衰減因子進行自適應(yīng)調(diào)整的方法，但反演結(jié)果受正則化因子初值和衰減因子的選值影響很大．Zhdanov等［16－17］提出一種正則化因子呈階梯狀下降的自適應(yīng)方法，選取數(shù)據(jù)擬合泛函和模型穩(wěn)定泛函的比值為正則化因子初值，并選取經(jīng)驗性的衰減因子，在數(shù)據(jù)擬合效率下降時對當(dāng)前正則化因子進行調(diào)整，具有較高的反演精度和計算效率．陳小斌等［18］提出一種完全自適應(yīng)方案（CMD 方案），選取數(shù)據(jù)擬合泛函與模型穩(wěn)定泛函的比值為每次迭代反演的正則化因子，該方案在整個迭代反演過程中始終強調(diào)數(shù)據(jù)擬合和模型穩(wěn)定的平衡，具有較高穩(wěn)定性．

為能在保持反演穩(wěn)定的基礎(chǔ)上進一步提高計算效率，本文以共軛梯度法求解帶光滑約束的MT 反問題為框架，引入了幾種正則化方式，包括L 曲線法、Occam 自適應(yīng)、CMD 自適應(yīng)以及Zhdanov自適應(yīng)算法，通過比較分析這些方法在不同初始模型條件下的反演情況，提出一種比較合理且在初始模型偏差較大時也能保證反演穩(wěn)定的自適應(yīng)正則化算法．

1 MT共軛梯度法光滑反演

Constable等［13］提出的Occam 反演是一種帶平滑約束的非線性最小二乘反演，該方法主要尋找滿足數(shù)據(jù)擬合的模型中令粗糙度極小的模型作為迭代結(jié)果，所得的反演結(jié)果比較光滑，對初始模型依賴較小，算法穩(wěn)定性高．

反演的目標函數(shù)為［3］

式中：Pα（m）為參數(shù)泛函；α為正則化因子；m為待解模型；φ（m）為數(shù)據(jù)擬合泛函；S（m）為模型穩(wěn)定泛函；離散形式下，Wd，Wm分別為數(shù)據(jù)和模型的加權(quán)矩陣；A為正演算子；d為觀測數(shù)據(jù)；We為穩(wěn)定泛函矩陣，此處為Occam反演中的粗糙化差分矩陣；e為防止分母為零的小值；mapr為先驗信息．

由于共軛梯度法具有不用求解大型線性方程組并且收斂速度快等優(yōu)點，本文采用共軛梯度法代替文獻［14］中的牛頓法求解模型參數(shù)，整個共軛梯度法的計算過程如下［3］：

式中：n為當(dāng)前迭代次數(shù)；rn為數(shù)據(jù)擬合殘差；A（mn），d分別為正演數(shù)據(jù)和觀測數(shù)據(jù)；mn為模型參數(shù)；sn為模型穩(wěn)定項；Wen為穩(wěn)定泛函矩陣；αn為正則化因子為梯度方向；Fmn為偏導(dǎo)數(shù)矩陣為系數(shù)項為共軛梯度方向為步長．

本文反演計算中，迭代終止條件為數(shù)據(jù)擬合誤差小于給定的閾值，若達不到此條件，則當(dāng)?shù)螖?shù)達到給定的最大次數(shù)時退出反演迭代．

2 模型試驗

2．1 模型試驗一

本文所采用的8 層層狀模型與文獻［18］中一致，在320～0．000 55Hz之間取40個計算頻點，在正演響應(yīng)數(shù)據(jù)上加入標準差為5%的隨機噪音．反演模型為41層，縱向上從0．01km 到100km 按對數(shù)等間距剖分，反演迭代的擬合閾值為5%，即與所加入的噪音等級一致．反演中選取Bostick反演結(jié)果為先驗信息．

2．1．1 多種正則化因子選取方法在不同初始模型下的反演結(jié)果

實際資料計算中，初始模型是不容易選擇的，故作者選擇了電阻率分別為1，10，100，1000，10 000 Ω·m的5個跨度比較大的均勻背景為初始模型來檢驗各種正則化反演方法的穩(wěn)定性．為了直觀地表現(xiàn)不同正則化因子選取方式對反演結(jié)果的影響，本文列出了不同初始模型條件下各方法的反演結(jié)果（圖1），包括無正則化反演，L 曲線正則化反演，Occam 自適應(yīng)正則化反演，CMD 自適應(yīng)正則化反演，Zhdanov自適應(yīng)正則化反演．因篇幅所限，各種正則化反演的流程以及參數(shù)選擇參見文獻［18－20］．需要指出的是，由L 曲線法計算后確定的正則化因子為0．1，Zhdanov自適應(yīng)算法的衰減因子取0．9．

從圖1可以看出，無正則化（圖1a）和L 曲線法（圖1b）的反演結(jié)果受初始模型的影響較大，只有在初始模型與實際模型偏差較小時才能得到較為理想的結(jié)果，當(dāng)初始模型偏差較大時，反演結(jié)果很不穩(wěn)定．Occam 自適應(yīng)反演結(jié)果（圖1c）受初始模型影響也很大，這主要是由于該方法雖然給定了很大的正則化因子為初值，但在迭代初期選擇使擬合誤差下降最快的正則化因子，正則化因子衰減很快，當(dāng)初始模型偏離較大時，模型穩(wěn)定泛函發(fā)揮的作用不夠，故所得反演結(jié)果不理想．CMD 方案的反演結(jié)果（圖1d）受初始模型影響相對較小，但由于正則化因子在計算中變化范圍較小，最終的數(shù)據(jù)擬合精度不夠，在較淺和較深區(qū)域的反演結(jié)果受初始模型的影響較大．Zhdanov自適應(yīng)方法在不同的初始模型下所得到的反演結(jié)果（圖1e）總體來說比較理想，但不同初始模型反演結(jié)果的一致性還不強，作者認為這主要是因為正則化因子的初始值在初始模型偏離較大時還是相對偏小．自適應(yīng)正則化的方法在初始模型偏離較小時可以得到與傳統(tǒng)定值方法相當(dāng)?shù)姆囱萁Y(jié)果，且在實際操作中更加便捷，總體來說它們對初始模型的依賴性要小于定值的方法．

圖1 不同初始模型下的不同正則化反演方法的反演結(jié)果Fig．1 Inversion results of different regularized methods by different initial models

2．1．2 本文改進的自適應(yīng)正則化因子選取方法的反演結(jié)果

作者認為一種行之有效的正則化因子選取方案應(yīng)該充分發(fā)揮模型穩(wěn)定泛函的約束作用，以減小反演結(jié)果對初始模型的依賴性．針對初始模型偏離較大或初始模型不易選擇的情況，結(jié)合上述模型試驗結(jié)果，作者以Zhdanov自適應(yīng)算法［3］為基本框架，提出了一種改進的自適應(yīng)正則化方案．

在Zhdanov自適應(yīng)算法中，選取數(shù)據(jù)擬合泛函和模型穩(wěn)定泛函的比值為初始值，作者認為在初始模型m1偏差較大時，這樣選取的正則化因子初值是偏小的，不能較好地發(fā)揮模型穩(wěn)定的約束效果，應(yīng)該選取更大的正則化因子α1為初值以增大模型穩(wěn)定泛函的作用，使其在初始模型偏離真值較大時也能充分穩(wěn)定，即在反演迭代的初期，充分發(fā)揮先驗信息的約束作用，而在反演迭代的后期，隨著正則化因子的減小進一步突出數(shù)據(jù)擬合的作用，得到一個較為精確的解．基于此，作者將這一比值乘以一個倍數(shù)k予以放大

Zhdanov自適應(yīng)正則化反演算法中，當(dāng)數(shù)據(jù)擬合效率下降時，即迭代前后數(shù)據(jù)擬合相對變化量Δφ＝（‖Wdrn‖2－‖Wdrn＋1‖2）／‖Wdrn‖2很小時，則按照一定的衰減因子q對當(dāng)前的正則化因子進行調(diào)整，Zhdanov指出衰減因子q的經(jīng)驗取值區(qū)間為［0．5，0．9］，但并沒有提出明確的衰減機制．基于此，作者試圖尋找一種可以兼顧反演穩(wěn)定性和計算效率的衰減機制．

作者同樣按照Δφ來衡量反演過程的迭代效率，并將迭代效率分為3種情況，低效率，高效率和過渡階段，在不同階段，選取不同的衰減方案．在滿足式4（a）時，數(shù)據(jù)擬合下降緩慢，說明當(dāng)前的正則化因子不再適應(yīng)當(dāng)前的反演迭代，應(yīng)選擇一個較小的衰減因子來提高反演精度；在滿足式4（b）時，數(shù)據(jù)擬合下降較快，應(yīng)選擇一個較大的衰減因子來保持反演穩(wěn)定性；在滿足式4（c）時，數(shù)據(jù)擬合下降一般，為了兼顧反演精度和穩(wěn)定性，作者選?。踧l，0．9］這一衰減區(qū)間．

作者設(shè)計了6種參數(shù)組合，選擇電阻率分別為1，10，100，1 000，10 000Ω·m 的均勻背景為初始模型，反演結(jié)果和正則化因子迭代曲線見圖2．結(jié)合圖1e、圖2，可以看出隨著正則化因子初值以及衰減因子的增大，反演的一致性逐漸變強，反演結(jié)果對初始模型的依賴性逐漸變?。@也驗證了作者對Zhdanov自適應(yīng)正則化反演結(jié)果（圖1e）跳躍的分析：直接選擇數(shù)據(jù)擬合泛函和模型穩(wěn)定泛函的比值為正則化因子初值在初始模型偏差較大的情況下是偏小的．在k取5，10，100時，即使衰減因子較?。▓D2b，2d，2f），模型穩(wěn)定泛函的約束作用還是可以得到較好發(fā)揮，這說明通過增大正則化因子初值使穩(wěn)定泛函在反演初期發(fā)揮較大貢獻，可以減小反演結(jié)果對初始模型的依賴性．此外，在正則化因子初值相同的情況下，選擇較大的衰減因子（圖2c，2e），可以在反演過程中保持模型穩(wěn)定泛函的約束作用，反演結(jié)果的一致性優(yōu)于選擇較小的衰減因子時的方案．

圖2 不同初始模型不同自適應(yīng)正則化因子選取方案的反演結(jié)果Fig．2 Inversion results with different adaptive regularization parameters by different initial models

此外，作者統(tǒng)計了k＝10，ql分別取0．4，0．5，0．6，0．7，0．8，0．9時的反演迭代次數(shù)、數(shù)據(jù)擬合誤差和模型逼近誤差，見圖3，同樣選擇電阻率分別為1，10，100，1 000，10 000Ω·m 的均勻背景為初始模型進行反演，數(shù)據(jù)擬合和模型逼近均用相對均方根誤差來統(tǒng)計．從圖3a可以看出，隨著正則因子初值和衰減因子的增大，反演達到閾值所需的迭代次數(shù)隨之增多．分析圖3b，3c，可以看出，ql取0．4時數(shù)據(jù)擬合和模型逼近情況跳躍最大，說明穩(wěn)定性有待提高；ql取0．9時數(shù)據(jù)擬合和模型逼近跳躍較小，說明穩(wěn)定性較高，但迭代次數(shù)較多，迭代效率不高；綜合考慮以上因素，作者認為ql取0．5～0．8時本模型的反演可以兼顧模型穩(wěn)定性和反演精度．

圖3 不同初始模型不同自適應(yīng)正則化因子選取方案的迭代次數(shù)、數(shù)據(jù)擬合和模型逼近統(tǒng)計Fig．3 Statistical graph of iterative number，data misfit and model approximation with different adaptive regularization parameters with different initial models

2．2 模型試驗二

為了驗證本文改進的自適應(yīng)算法在復(fù)雜模型條件下的適應(yīng)性，作者設(shè)計了如圖4a所示的二維模型．采用k＝10，ql＝0．6，qh＝1．0參數(shù)組合對該模型進行了電場極化模式反演，初始模型分別為100 和10 000Ω·m 兩種方案，反演結(jié)果見圖4b，4c，此外，還進行了Zhdanov自適應(yīng)正則化反演（圖4d，4e）和正則化因子取定值0．1時的正則化反演（圖4f，4g）予以對比．

圖4 不同初始模型下的不同正則化方法的反演結(jié)果Fig．4 Inversion results of different regularized methods by different initial models

本文改進的自適應(yīng)算法在兩種初始模型條件下的模型逼近誤差分別為0．208和0．206，可見初始模型對反演結(jié)果影響很小，所得反演結(jié)果一致性很高；Zhdanov自適應(yīng)算法在兩種初始模型條件下的模型逼近分別為0．217和0．397，所得反演結(jié)果一致性一般，在初始模型為10 000Ω·m 時，個別測點的反演已經(jīng)出現(xiàn)不穩(wěn)定的情況，跳越較大；正則化因子取定值0．1時正則化反演在兩種初始模型條件下的模型逼近分別為0．211 和1．285，所得反演結(jié)果差異很大，初始模型為10 000Ω·m 時，已不能反演地下真實的電阻率分布，偏差很大，這也說明了經(jīng)驗定值法的局限性．可見本文改進的自適應(yīng)算法在復(fù)雜模型條件下依舊可以保持反演的穩(wěn)定性，對初始模型依賴小，而且可以保持較高的反演精度．

3 結(jié)論

正則化因子的選擇決定著正則化反演的結(jié)果，如何選擇合適的正則化因子一直是反演研究中的重點和熱點．本文通過比較多種正則化因子的選取方法，分析了各種方法在給定不同初始模型條件下的表現(xiàn)，提出了一種改進的自適應(yīng)正則化方法，得到了以下結(jié)論：

（1）本文提出的改進的自適應(yīng)正則化方法，通過選擇較大正則化因子初值和自動控制衰減因子，相比其他幾種自適應(yīng)正則化算法對初始模型的依賴小，可在不同初始模型的條件下得到精度及穩(wěn)定性都很高的解，具有明顯的優(yōu)越性．

（2）本文提出的自適應(yīng)正則化方案是基于共軛梯度法求解的光滑反演，但可推廣至其他地球物理反演方法以及高維反演．建議在不強調(diào)計算時間代價的情況下，應(yīng)選擇較大的正則化因子初值和衰減因子；在必須考慮計算效率的情況下，建議選擇較大的正則化因子初值和較小衰減因子以兼顧反演的穩(wěn)定性和計算效率．

［1］ Tikhonov A N．Regularization of incorrectly posed problem［M］．Dokl：Soviet Math，1963．

［2］ Oldenbury D，Li Yaoguo．Inversion for applied geophysics：a tutorial［J］．Near－surface Geophysics，2005，5：89．

［3］ Zhdanov M S．Geophysical inverse theory and regularization problems［M］．［S．l．］：Elsevier Science，2002．

［4］ Wahba G．Spline models for observational data ［M］．Philadelphia：SIAM，1990．

［5］ Lin Y Z，Wohlberg B．Application of the UPRE method to optimal parameter selection for large scale regularization problems［J］．Southwest symposium on Image Analysis and Interpretation，2008：89．

［6］ Wahba G．Practical approximate solutions to linear operator equations when the data are noisy ［J］．SIAM Journal on Numerical Analysis，1977，14（4）：651．

［7］ Golub G H，Heath M T，Wahba G．Generalized cross－validation as a method for choosing a good ridge parameter ［J］．Technometrics，1979，21：215．

［8］ Haber E，Oldenburg D．A GCV based method for nonlinear illposed problems［J］．Computational Geosciences，2000，4（1）：41．

［9］ Chung Julianne，Nagy G J．A weighted－GCV method for Lanczos－Hybrid regularization ［J］．Electronic Transactions on Numerical Analysis，2008，28：149．

［10］ Hansen P C，O’Leary D P．The use of the L－curve in the regularization of discrete ill－posed problems［J］．SIAM Journal on Numerical Analysis，1993，14（6）：1487．

［11］ Hansen P C．Regularization tools：a matlab package for analysis and solution of discrete ill－posed problems ［J］．Numerical Algorithm，1994，14（6）：1．

［12］ 吳小平，徐果明．利用共軛梯度方法的電阻率三維反演研究［J］．地球物理學(xué)報，2000，43（3）：420．WU Xiaoping，XU Guoming．Study on 3－D resistivity inversion using conjugate gradient method ［J］．Chinese Journal of Geophysics，2000，43（3）：420．

［13］ 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：289．

［14］ deGroot－Hedlin C，Constable S．Occam’s inversion to generate smooth two－dimensional models from magnetotelluric data［J］．Geophysics，1990，55（12）：1613．

［15］ 吳小平，徐果明．大地電磁數(shù)據(jù)的Occam 反演改進［J］．地球物理學(xué)報，1998，41（4）：547．WU Xiaoping，XU Guoming．Improvement of Occam’s inversion for MT data ［J］．Chinese Journal of Geophysics，1998，41（4）：547．

［16］ Kumar A，Wan L，Zhdanov M S．Regularization analysis of three－dimensional magnetotelluric inversion ［C］／／2007 SEG Annual Meeting． San Antonio： Society of Exploration Geophysicists，2007：2007－0482．

［17］ Gribenko A，Zhdanov M S．Rigorous 3D inversion of marine CSEM data based on the integral equation method ［J］．Geophysics，2007，72（2）：WA73．

［18］ 陳小斌，趙國澤，湯吉，等．大地電磁自適應(yīng)正則化反演算法［J］．地球物理學(xué)報，2005，48（4）：937．CHEN Xiaobin，ZHAO Guoze，TANG Ji，et al．An adaptive regularized inversion algorithm for magnetotelluric data ［J］．Chinese Journal of Geophysics，2005，48（4）：937．

［19］ 陳曉，于鵬，張羅磊，等．大地電磁與地震正則化同步聯(lián)合反演［J］．地震地質(zhì)，2010，32（3）：402．CHEN Xiao，YU Peng，ZHANG Luolei，et al．Regularized synchronous joint inversion of MT and seismic data ［J］．Seismology and Geology，2010，32（3）：402．

［20］ 陳曉，于鵬，張羅磊，等．地震與大地電磁測深數(shù)據(jù)的自適應(yīng)正則化同步聯(lián)合反演［J］．地球物理學(xué)報，2011，54（10）：2673．CHEN Xiao，YU Peng，ZHANG Luolei，et al．Adaptive regularized synchronous joint inversion of MT and seismic data［J］．Chinese Journal of Geophysics，2011，54（10）：2673．