樸英哲，李桐林，劉永亮

1.金策工業(yè)綜合大學資源探測工學系，朝鮮 平壤 999093

2.吉林大學地球探測科學與技術學院，長春 130026

0 引言

大地電磁反演的非唯一性是眾所周知的［1－3］，而Occam反演是克服該缺陷的方法之一。在大地電磁（magnetotellurics，MT）數(shù)據解釋中，Constable等［2］初次將Occam反演這種說法用于一維反演。Occam反演是光滑模型反演，具有較好的穩(wěn)定收斂性和結果可靠性。另外，Occam反演能夠引入先驗信息來去掉對已知構造邊界處的光滑度或者加上對同類電性單元的電阻率差的限制?；谶@些優(yōu)點，Occam反演廣泛應用于 MT二維解釋［3－4］。

Siripunvaraporn 等［5－7］研 究 出 了 數(shù) 據 空 間Occam反演方法，并且成功地應用于MT二維、三維反演。他們在數(shù)據空間方法中，又結合了CG的優(yōu)點而研究出了 DCGOCC［8］。張羅磊等［9］還提出了光滑模型與尖銳邊界相結合的反演方法。對Occam反演普遍而詳細的論述可見于文獻［2－3，5－7］中。

Occam反演不但被應用于MT數(shù)據解釋而且被應用于多種地球物理勘探反演解釋當中［10－17］。MT二維數(shù)據Occam反演方法［3］是其中最有代表性的例子。拉格朗日乘子是介于模型光滑和數(shù)據擬合間的折衷參數(shù)，每次迭代反演為了求取適當?shù)睦窭嗜粘俗有枰M行多次正演計算，尤其在接近收斂時更是如此。為此，不少研究人員提出了直接求取拉格朗日乘子的方法。

吳小平等［18］提出了每次迭代以固定的比率減少拉格朗日乘子的方法，還指出雖然這種反演的結果非最光滑模型，但因為觀測數(shù)據是反演解釋的第一手資料，而模型光滑作為反演約束條件僅是穩(wěn)定迭代的手段，只有使理論數(shù)據與實際數(shù)據盡可能一致才能分辨所有的構造特征，尤其對精確數(shù)據的反演更是如此。但這種方法忽視了因觀測數(shù)據存在噪聲而產生多余構造的可能性，以及對精確數(shù)據來說，調整正演與觀測數(shù)據之間的擬合差期望值是更合理的。吳小平的求取方法速度快，在3DMT正演需要較長時間的情況下，反演多采取了該方法［8，19］。MT三維反演的拉格朗日乘子選取方法與吳小平等的方法有一些差別之處，就是隨著擬合差的變化而減小或增大拉格朗日乘子。

張羅磊等［9］提出的方法是根據反演目標泛函中數(shù)據誤差部分和模型粗糙度部分占用的比重選取初始值，每次迭代隨模型粗糙度的變化，在總體粗糙度沒有減少時，以粗糙度變化率來減少乘子。這種方法也快，但容易陷入局部極小值。

另外，關于與拉格朗日乘子類似的Tikhonov正規(guī)化因子的選取方法有幾種，如基于離差原理（discrepancy principle）的 方 法［20］、廣 義 交 叉 驗 證（generalized cross validation）方 法［21］、L－曲 線（L－curve）法［22］和 U－曲線（U－curve）法［23］。這些方法從每次迭代需要反復計算正演的這一點上來看與deGroot－Hedlin等［3］的 Occam2DMT 相似（以下將deGroot－Hedlin等［3］的方法稱為 Occam2DMT）。

筆者首先介紹了Occam反演和Occam2DMT的拉格朗日乘子求取方法，然后提出了改進的方法。接著通過幾種模型實驗證明了改進的方法比原方法更有效，比吳小平等［18］的方法更穩(wěn)定。

1 Occam反演中拉格朗日求取方法改進

1.1 Occam反演方法及拉格朗日乘子求取方法

Occam反演是在一定的擬合誤差標準下求使模型粗糙度最小的解。因此，反演的目標泛函以模型粗糙度、擬合差及拉格朗日乘子構成：

其中：m為模型向量；‖▽m‖2為模型粗糙度；d為觀測數(shù)據向量；F為正演算子；W為利用數(shù)據標準差進行規(guī)一化的矩陣；‖Wd－WF（m）‖為數(shù)據擬合差（以下用X表示）；X＊為擬合差的期望值；μ為拉格朗日乘子。

這里模型向量的泛函是非線性，為此，在初始模型m1附近作線性化，用迭代方法求解。第二次迭代模型近似為

由式（2）可知，m2為μ的函數(shù)，數(shù)據擬合差也是μ的函數(shù)。從m1得到m2時，隨μ值從0到無窮大變化，模型m2沿模型空間中的一定軌道移動。模型空間中的每個模型都對應相應的擬合差，所以可以想到模型空間中的擬合差等值線（圖1）。

Occam2DMT首先用Brent方法［24］找到一個μ2，使得數(shù)據擬合差極小，即

圖1 μ值和擬合差的關系Fig.1 Relationship betweenμand root mean square misfit

Brent方法［24］需要事先確定包含極小值的區(qū)間（a，b），為此，在初始μ0附近找到a、b及μ1（μ1∈（a，b）），令

若X（μ1）或X（μ2）小于誤差限，利用 van Wijngaarden－Dekker－Brent方法［24］搜索m2軌道與誤差限等值線交叉的最大的μ值，令μ＊為此值，即

若X（μ2）＞X＊，令μ＊＝μ2。

此μ＊值是最佳拉格朗日乘子。按上述原理，Occam2DMT求取μ＊的方法由以下3個步驟組成。

①確定極小值區(qū)間：找到滿足式（4）的a、b及μ1。若X（μ1）≥X＊，轉移到步驟②；否則轉移到步驟③。

②X極小化：在區(qū)間（a，b）中，用Brent方法搜索滿足式（3）的極小點μ2。若X（μ2）＜X＊，轉移到步驟③；否則令μ＊＝μ2，終止。

③交叉點搜索：用WDB方法搜索滿足式（5）的μ＊。

經以上3個步驟確定μ＊之后，將此值代入式（2）計算m2。以上3個步驟需要進行反復正演，正演次數(shù)直接關系到反演的計算量。

1.2 改進的拉格朗日乘子求取方法

一般來說，Occam反演由2個階段組成：第一階段是使擬合差減小到擬合差期望值，第二階段是在保持擬合差為期望值的同時，搜索粗糙度最小的模型。Occam2DMT在這2個階段中，都利用上述的由3個步驟組成的方法求取μ＊。

在第一階段的每次迭代中，除了最后迭代，擬合差都未達到期望值，所以在模型空間中m2軌道未交叉X＊等值線，未經過步驟③，搜索μ2達到目的。但在第一階段的最后迭代和第二階段的每次迭代中，都經過步驟③，因此目的是搜索μ＊的。那么此時能夠使求取μ＊的方法優(yōu)化。

實際上，步驟①的目的只是求函數(shù)X（μ）的極小值點；若在步驟①的反復計算中，有一點的X函數(shù)值小于誤差限時（圖2a），盡管未確定極小值區(qū)間，也可以直接轉移到步驟③搜索大于μ1的交叉點μ＊。如此，可以排除不必要的計算，且結果μ＊值不受任何影響。

圖2 求取拉格朗日乘子Fig.2 Choosing Lagrange multiplier

在第二階段中，雖然模型m1的相應擬合差小于誤差限，但也有步驟①a、b及μ1的相應擬合差都大于誤差限（圖2b），此時得進行步驟②。不過步驟②不必找到極小點，只需要有一個點函數(shù)值小于X＊。因此，在步驟②的反復計算中一旦有一點的X函數(shù)值小于誤差限，就可以終止X函數(shù)極小化，開始搜索交叉點。如此未影響μ＊值，也排除多余的計算。

很明顯，上述改進更符合于Occam思想，并且其求解空間與原方法一致。

下一個問題是如何設定每次迭代的初始μ0。也許μ0越接近最佳值μ＊，計算量越少。第一迭代的初 始 值由使用者預定。第一迭代以后，Occam2DMT令μ0為前次迭代的μ2。在反演的第一階段中μ2與μ＊一致，為此這種選擇是適當?shù)?。但在反演第二階段中μ2與μ＊稍微差別，所以這種選擇可能不恰當。

根據模型實例，在模型光滑階段中μ＊值的變化較小。因此令初始μ0為前次迭代的μ＊也許更合理。筆者考慮到試算模型的μ＊值變化特征，將模型光滑階段中第i次迭代的初始值設定為，其 中為前兩次迭代的μ＊值。若i＜3，令

經計算，本文方法與Occam2DMT求取的μ＊的誤差很?。ㄐ∮?0－5），而且能夠排除多余的正演計算。本方法的算法如圖3。

圖3 本文方法第i次迭代的算法Fig.3 Flow chart of the ith iteration of this inversion

2 模型實驗與結果分析

為了比較本文方法和原方法以及吳小平等［18］的方法，構建了如圖4所示的8種地電模型。模型1和模型2是電阻率100Ω·m均勻半空間中存在矩形異常體。模型1異常體的電阻率為1 000Ω·m，頂面埋深為7km；模型2異常體的電阻率為10、1 000Ω·m，頂面埋深為10km。模型3是電阻率為100Ω·m的圍巖中存在電阻率10Ω·m傳導性侵入 巖（deGroot－Hedlin 等［3］）。 模 型 4 和 5 與deGroot－Hedlin等［4］反演的模型相同。模型6是與Siripunvaraporn等［7］相似的鄰近的不同電阻率塊體模型（圍巖電阻率為100Ω·m，異常體電阻率分別為10、1 000Ω·m）。對于所有模型利用趨膚深度來進行網格剖分和劃分網格邊界［25］。大地電磁場受地形影響［26］，筆者將模型7和模型8設定為起伏地形模型。對于起伏地形模型，為了保證正演精度，斜坡附近采用更周密網格。在起伏地形的坡角處，不同研究人員的輔助場計算方法不同［27］；為此，在坡角處未安置測點。測量數(shù)據為TE、TM 2種極化模式下的視電阻率和阻抗相位（模型1：6個測點、6個頻點（0.486～0.002Hz的對數(shù)間隔）；模型2—模型6：11個測點、16個頻點（1～0.001Hz的對數(shù)間隔），模型7、模型8：11個測點、16個頻點（100～0.1Hz的對數(shù)間隔）），并在數(shù)據中加入了2%的隨機噪聲。反演初始模型為1Ω·m的均勻半空間，μ0值為5，X＊為1.1。

首先，原方法和本文方法比起來，每次反演迭代的μ＊值變化完全一致，當然反演結果也一致，只是在正演次數(shù)上有差別。表1為8個模型的原方法與調整μ0前和調整μ0后改進方法的正演次數(shù)比較。

表1 正演次數(shù)比較Table 1 Comparison of forward modeling number

可以看到所有模型正演次數(shù)均減少了20%～50%。調整μ0后正演次數(shù)小于調整前（除模型1），不過其效果并不大。

其次，比較本文方法（調整μ0后）和吳小平等的方法。吳小平等［18］方法的關鍵是如何選取減小拉格朗日乘子的比率λ及其下限μmin，通常取決于經驗。本文確定Occam2DMT反演迭代的μ＊值是最佳選擇，使比率和下限符合μ＊值的變化。

圖4 地電模型Fig.4 Synthetic models

模型1—4中μ＊的變化和模型5—8中μ＊的變化稍有差別（圖5）。因此，把8種模型分成2組。對第一組（模型1—4，圖5a），令

對于第一組模型，吳小平等［18］的方法比本文方法收斂快，反演結果雖非最光滑模型，但其粗糙度和最光滑模型（Occam2DMT反演結果）的粗糙度差別很小。這種效果是由于適當選取了λ和μmin。

但是，用同樣的參數(shù)反演第二組模型，結果分散。因此，對第二組（模型5—8，圖5b），令

2個方法的迭代次數(shù)不相等，所以用反演所需總的時間來對比2種方法（圖6）。

圖5 μ＊隨迭代次數(shù)的變化及λ選擇Fig.5 Variation ofμ＊along iteration number and choosingλ

圖6 第二組模型數(shù)據擬合差變化比較Fig.6 Comparison of root mean square misfit variation for the second model group

對于所有的模型，本方法穩(wěn)定收斂。相比之下，吳小平等的方法收斂速度慢。本方法大約100s時已收斂到誤差限內，此后計算是找最光滑模型，但吳小平等［18］的方法大約250s時才收斂（圖6a），甚至未收斂（圖6b）。對于模型7和模型8，吳小平等方法雖收斂到誤差限，但其收斂很不穩(wěn)定且收斂時間比本方法更長（圖6c，d）。

對于模型6的反演結果，本文方法和吳小平的方法都找到設計異常體，不過在吳小平的方法反演剖面中地表附近（水平距離13～15km）出現(xiàn)了電阻率較低的異常體（圖7）。這是因為模型6的反演沒收斂，同時也說明了最光滑模型和非最光滑模型的區(qū)別。

圖7 模型6反演結果對比Fig.7 Comparison of model 6inversion results

對不同的幾組λ、μmin反復進行試算，對第二組模型吳小平等方法的結果也沒有改善。其原因是每次迭代都要減小拉格朗日乘子。其實從圖5可以看出，適當?shù)睦窭嗜粘俗硬⒎鞘且宦蓽p小的。吳小平等指出該方法對其參數(shù)取值無太嚴格要求，但他的例子是一維的，且二維反演的非惟一性比一維反演嚴重。這也是吳小平等方法效果不好的原因。

3 野外數(shù)據反演結果

用本文方法處理了本溪—集安地區(qū)大地電磁第5線數(shù)據（圖8）。使用的儀器是加拿大鳳凰公司的V5－2000系列電磁儀。測線的方向大致從北西到南東。在測線上共有29個測點，測點間平均距離為大約5km。測量數(shù)據中采取了13個頻點（360.0～0.7Hz的對數(shù)間隔）。根據采集物性樣本的電阻率測量結果，該地區(qū)的地下電阻率范圍為數(shù)百到數(shù)萬Ω·m?？紤]到地形起伏與趨膚深度，把剖面網格單元大小設定為寬250～750m，高100～2 000m。反演初始模型為1 000Ω·m的均勻介質，μ0值為5，X＊為3.5。

圖8 5線的反演結果Fig.8 Inversion result of line 5

本文方法和Occam2DMT方法的反演總迭代都為8次，結果模型的粗糙度都為29.97。不過正演次數(shù)分別為52和78次，計算時間分別為26h20 min和35h30min，反演結果一致，如圖8所示。

4 結論

雖然Occam2DMT具有收斂穩(wěn)定性和結果可靠性的優(yōu)點，但由于其利用了Press等的算法，在靠近解時正演次數(shù)增大，所以其計算時間較長。為了縮短反演時間，減少正演次數(shù)很重要。以固定的比率減小或增大拉格朗日乘子的方法，雖因在每次迭代中只需一次進行正演而很快，但在實際應用中，若使用者不進行人為調整反演參數(shù)，容易造成收斂失敗或虛偽構造。筆者在求取拉格朗日乘子的一維搜索中排除多余的正演計算，使得Occam反演速度加快。本文通過對幾種模型計算和野外數(shù)據反演，得出如下結論。

1）本文方法總能獲得與Occam2DMT方法一致的解，在反演的擬合差下降到滿足期待值階段和光滑模型階段，計算效率有明顯提高。根據模型實驗，可以減少正演次數(shù)20%～50%。

值得提出的是，當觀測數(shù)據含有較強干擾噪聲或地下電阻率較復雜時，無法求得在預定誤差范圍內的解，此時，本文方法不能有效地減少正演次數(shù)來減少計算時間。

2）本文方法的反演結果是粗糙度最低的模型。結果不但可靠，而且其收斂性也很穩(wěn)定。

（References）：

［1］劉國棟，陳樂壽.大地電磁測深研究［M］.北京：地震出版社，1984.Liu Guodong，Chen Leshou.The Study of Magneto－telluric Sounding［M］.Beijing：Seismological Publishing House，1984.

［2］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.

［3］deGroot－Hedlin C D，Constable S C.Occam’s Inversion to Generate Smooth Two Dimensional Models from Magnetotelluric Data［J］.Geophysics，1990，55（12）：1613－1624.

［4］deGroot－Hedlin C D，Constable S C.Inversion of Magnetotelluric Data for 2DStructure with Sharp Resistivity Contrasts［J］.Geophysics，2004，69（1）：78－86.

［5］Siripunvaraporn W，Egbert G.An Efficient Data－Subspace Inversion Method for 2DMagnetotelluric Data［J］.Geophysics，2000，65（3）：791－803.

［6］Siripunvaraporn W，Uyeshima M，Egbert G.Three－Dimensional Inversion for Network－Magnetotelluric Data［J］.Earth Planets Space，2004，56：893－902.

［7］Siripunvaraporn W，Egbert G，Lenbury Y，et al.Three－Dimensional Magnetotelluric Inversion：Data－Space Method［J］.Phys Earth Planet Inter，2005，150：3－14.

［8］Siripunvaraporn W，Sarakorn W.An Efficient Data Space Conjugate Gradient Occam’s Method for Three－Dimensional Magnetotelluric Inversion［J］.Geophys J Int，2011，186，567－579，doi：10.1111/j.1365－246X.2011.05079.x.

［9］張羅磊，于鵬，王家林，等.光滑模型與尖銳邊界結合的MT二維反演方法［J］.地球物理學報，2009，52（6）：1625－1632.Zhang Luolei，Yu Peng，Wang Jialin，et al.Smoothest Model and Sharp Boundary Based Two－Dimen－sional Magnetotelluric Inversion［J］.Chinese Journal of Geophysics，2009，52（6）：1625－1632.

［10］劉羽，王家映，孟永良.基于PC機群的大地電磁Occam 反演并行計算研究［J］.石油物探，2006，45（3）：311－315.Liu Yu，Wang Jiaying，Meng Yongliang.PC Cluster Based Magnetotelluric 2－D Occam’s Inversion Parallel Calculation［J］.GPP，2006，45（3）：311－315.

［11］Parker R L.Geophysical Inverse Theory［M］.New Jersey：Princeton University Press，1994.

［12］deGroot－Hedlin C D.Inversion for Regional 2－D Resistivity Structure in the Presence of Galvanic Scatterers［J］.Geophys J Int，1995，122：877－888.

［13］LaBrecque D J，Ward S H.Two－Dimensional Cross－Borehole Resistivity Model Fitting［J］.Geotechnical and Environmental Geophysics，1990，1：51－57.

［14］Songkhun Boonchaisuk，Chatchai Vachiratienchai，Weerachai Siripunvaraporn.Two－Dimensional Direct Current（DC）Resistivity Inversion：Data Space Occam’s Approach［J］.Physics of the Earth and Planetary Interiors，2008，168：204－211.

［15］翁愛華.Occam反演及其在瞬變電磁測深中的應用［J］.地質與勘探，2007，42（5）：74－76.Weng Aihua.Occam’s Inversion and Its Application to Transient Electromagnetic Method［J］.Geology and Prospecting，2007，42（5）：74－76.

［16］Huang Z X，Su W，Peng Y J，et al.Rayleigh Wave Tomography of China and Adjacent Regions［J］.J Geophys Res，2003，108（B2）：ARTN 2073.

［17］Greenhalgh S A，Bing Z，Green A.Solutions，Algorithms and Inter－Relations for Local Minimization Search Geophysical Inversion［J］.J Geophys Eng，2006，3：101－113.

［18］吳小平，徐果明.大地電磁數(shù)據的Occam反演改進［J］.地球物理學報，1998，41（4）：547－554.Wu Xiaoping，Xu Guoming.Improvement of Occam’s Inversion for MT Data［J］.Chinese Journal of Geophysics，1998，41（4）：547－554.

［19］Newman G A，Alumbaugh D L.Three Dimensional Magnetotelluric Inversion Using Non－Linear Conjugate Gradients［J］.Geophys J Int，2000，140：410－424.

［20］Pereverzev S.Morozov’s Discrepancy Principle for Tikhonov Regularization of Severely Ill－Posed Problem in Finite－Dimensional Subspaces［J］.Numerical Functional Analysis and Optimization，2000，21（7）：901－916.

［21］Haber E，Oldeburg D.A GCF Based Method for Nonlinear Ill－Posed Problems［J］.Computational Geosciences，2000，4（1）：41－63.

［22］Hansen P C，Leary D P.The Use of the L－Curve in the Regularization of Discrete Ill－Posed Problems［J］.SIAM Journal on Scientific Computing，1993，14（6）：1487－1503.

［23］Stando D K，Rudnicki M.Regularization Parameter Selection in Discrete Ill－Posed Problems：The Use of the U－Curve［J］.International Journal of Applied Mathematics and Computer Science，2007，17（2）：157－164.

［24］Press H W，Teukolsky A S，Vetterling T W，et al.Numerical Recipes in Fortran 77［M］.New York：Cambridge University Press，1997.

［25］湯井田，薛帥.MT有限元模擬中截斷邊界的影響［J］.吉林大學學報：地球科學版，2013，43（1）：267－274.Tang Jingtian，Xue Shuai.Influence of Truncated Boundary in FEM Numerical Simulation of MT［J］.Journal of Jilin University：Earth Science Edition，2013，43（1）：267－274.

［26］趙廣茂，李桐林，王大勇，等.基于二次場二維起伏地形MT有限元數(shù)值模擬［J］.吉林大學學報：地球科學版，2008，38（6）：1055－1059.Zhao Guangmao，Li Tonglin，Wang Dayong，et al.Secondary Field－Based Two－Dimensional Topographic Numerical Simulation in Magnetotellurics by Finite Element Method［J］.Journal of Jilin University：Earth Science Edition，2008，38（6）：1055－1059.

［27］Li Shenghui，Booker J R，Aprea C.Inversion of Magnetotelluric Data in the Presence of Strong Bathymetry/Topography［J］.Geophysical Prospecting，2008，56：259－268.