在大地電磁二維Occam反演中求取拉格朗日乘子方法改進(jìn)

樸英哲，李桐林，劉永亮

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

2.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

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

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

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

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

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

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

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

1 Occam反演中拉格朗日求取方法改進(jìn)

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

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

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

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

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

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

圖1 μ值和擬合差的關(guān)系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個(gè)步驟組成。

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

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

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

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

1.2 改進(jìn)的拉格朗日乘子求取方法

一般來(lái)說(shuō)，Occam反演由2個(gè)階段組成：第一階段是使擬合差減小到擬合差期望值，第二階段是在保持?jǐn)M合差為期望值的同時(shí)，搜索粗糙度最小的模型。Occam2DMT在這2個(gè)階段中，都利用上述的由3個(gè)步驟組成的方法求取μ＊。

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

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

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

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

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

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

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

經(jīng)計(jì)算，本文方法與Occam2DMT求取的μ＊的誤差很小（小于10－5），而且能夠排除多余的正演計(jì)算。本方法的算法如圖3。

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

2 模型實(shí)驗(yàn)與結(jié)果分析

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

首先，原方法和本文方法比起來(lái)，每次反演迭代的μ＊值變化完全一致，當(dāng)然反演結(jié)果也一致，只是在正演次數(shù)上有差別。表1為8個(gè)模型的原方法與調(diào)整μ0前和調(diào)整μ0后改進(jìn)方法的正演次數(shù)比較。

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

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

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

圖4 地電模型Fig.4 Synthetic models

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

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

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

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

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

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

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

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

圖7 模型6反演結(jié)果對(duì)比Fig.7 Comparison of model 6inversion results

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

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

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

圖8 5線的反演結(jié)果Fig.8 Inversion result of line 5

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

4 結(jié)論

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

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

值得提出的是，當(dāng)觀測(cè)數(shù)據(jù)含有較強(qiáng)干擾噪聲或地下電阻率較復(fù)雜時(shí)，無(wú)法求得在預(yù)定誤差范圍內(nèi)的解，此時(shí)，本文方法不能有效地減少正演次數(shù)來(lái)減少計(jì)算時(shí)間。

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

（References）：

［1］劉國(guó)棟，陳樂(lè)壽.大地電磁測(cè)深研究［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］張羅磊，于鵬，王家林，等.光滑模型與尖銳邊界結(jié)合的MT二維反演方法［J］.地球物理學(xué)報(bào)，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機(jī)群的大地電磁Occam 反演并行計(jì)算研究［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］翁愛(ài)華.Occam反演及其在瞬變電磁測(cè)深中的應(yīng)用［J］.地質(zhì)與勘探，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ù)據(jù)的Occam反演改進(jìn)［J］.地球物理學(xué)報(bào)，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有限元模擬中截?cái)噙吔绲挠绊懀跩］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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］趙廣茂，李桐林，王大勇，等.基于二次場(chǎng)二維起伏地形MT有限元數(shù)值模擬［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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.