基于ILU-BICGSTAB算法的起伏地形MT二維正演模擬

郭勝男， 熊 彬， 羅天涯， 吳延強， 陳 欣

(桂林理工大學 地球科學學院， 桂林 541006)

基于ILU-BICGSTAB算法的起伏地形MT二維正演模擬

郭勝男， 熊 彬， 羅天涯， 吳延強， 陳 欣

(桂林理工大學 地球科學學院， 桂林 541006)

針對起伏地形對異常體的探測效果的影響，利用有限單元法，基于MTLAB平臺編寫了大地電磁二維正演程序，為更好地模擬實際地形，采用矩形網格中進一步三角細化的剖分方式，并利用不完全LU分解的預條件BICGSTAB算法提高計算速度與精度，最后通過建立多種地電模型分析討論了起伏地形對大地電磁場的影響，驗證了該方法的實用性與準確性。

大地電磁； 有限單元； 起伏地形； 正演

0 引言

大地電磁測深是一種利用天然電磁場進行電磁測深的方法[1]，因具有較大的穿透深度和分辨能力等優點，從而廣泛應用于天然氣勘探和地殼、上地幔研究中。而在實際勘探中，地形條件因地而異，往往都是起伏不定，若忽略了起伏地形對大地電磁測深的影響，必將影響反演解釋效果，因此，對于復雜地形大地電磁場的正演研究是很有必要的。

國外很早就開展了起伏地形對大地電磁場影響的研究工作，Wannamaker等[2]首先采用了矩形中進一步三角細化單元剖分的有限單元法，分析討論了復雜地形對視電阻率及相位的影響；Franke等[3]為了更好地模擬實際地形，使用了非結構化任意三角的剖分方法；Myung等[4]采用矢量有限單元法完成了起伏地形下的MT三維正演問題。在國內，徐世浙等[5-6]使用有限單元法及邊界單元法完成了起伏地形大地電磁二維正演模擬，并針對地形對MT的影響給出了合理的解釋；陳小斌[7]提出有限元直接迭代法，分析談論了起伏地形對大地電磁測深的影響，并對視電阻率產生畸變的原因作出新的解釋；劉小軍等[8]利用求解二次場的有限單元法完成了起伏地形二維正演數值模擬；柳建新等[9]利用自適應有限元法模擬了起伏地形對大地電磁場的影響。

筆者在前人研究的基礎上，采用矩形進一步三角細化的網格剖分方式，系統方程求解使用ILU的預條件BICGSTAB算法，實現了起伏地形二維大地電磁場的模擬。

1 理論基礎

1.1 二維大地電磁正演邊值問題

兩種極化模式下二維大地電磁場所滿足的邊值問題為[10-11]：

(1)

TM模式下，

(2)

上邊界AB直接取在地面上，并以該出的u為1單位，左右邊界AD和BC分別取距離局部不均勻體足夠遠的地方，CD為下邊界，取在離地面足夠遠的地下。

TE模式下，

(3)

上邊界AB離地面足夠遠，使異常場在此邊界上為零，下邊界及左右邊界與TM模式相同。

1.2 二維大地電磁正演變分問題

MT二維正演的變分問題為式(4)。

(4)

1.3 有限單元法

采用的區域剖分方式是在矩形單元基礎上引入兩條對角線，將其細分為4個小三角形，然后在小三角形上進行線性插值，這樣做的優點是既能適應傾斜物性界面，又能節省存儲空間，增加計算效率。

將式(4)分解為各單元積分之和。

(5)

對式(5)中的各個單元及邊界分別積分，消去中間節點，得4×4矩陣，分別擴展為全體節點矩陣，然后求和，對泛函取極值，得線性方程組式(6)。

ku=0

(6)

將上邊界條件加入，即可求得各節點的u值，隨后計算輔助場[11]，便可得到地表的視電阻率及阻抗相位。

1.4 預處理BICGSTAB算法

BICGSTAB算法的主要思想是基于雙邊Lanczos算法和基于殘差正交子空間的迭代方法，與傳統的經典迭代法相比具有計算速度快、精度高、穩定性強等優點，從而廣泛應用于大型稀疏復系數矩陣的求解。預條件矩陣采用ILU技術進行分解，其基本求解思路為[12]：

2)若k≤kmax且ε≥εmax，轉步驟3)，否則停止，輸出xk。

4)若k=1，則ρ1=r0，否則，計算βk-1=ρk-1·αk-1/ρk-2ωk-1,ρk=rk+1+βk-1(ρk-1-ωk-1Vk-1)。

2 正演模擬與分析

2.1 層狀均勻介質模型

將箱式電阻爐溫度升至500 ℃時，放入有涂層和無涂層的Q235鋼片樣品，放置一定時間后從爐內取出試樣，冷卻至室溫，稱重，以樣品氧化增重速率V(O)表征試樣的被氧化能力.V(O)值越小，說明試樣抗氧化性能越好.V(O)值的計算公式如式 (1).

首先利用具有解析解的三層層狀介質模型，檢驗本文程序及算法的準確性。模型的基本參數為：ρ1=100Ω·m，h1=1 000m，ρ2=10Ω·m，h2=500m,ρ3=100Ω·m。如圖1和圖2所示，本文程序所計算出的兩種極化模式下視電阻率、相位與解析解相當吻合，平均相對誤差僅為0.2%，而利用傳統的經典算法求出的平均相對誤差則高達0.9%， 這驗證了本文程序及算法的可行性。

2.2 地壘起伏模型

如圖3所示，為地壘起伏地形，上寬為400m，下寬為1 600m，高度為500m，研究區域包括空氣及地形線以下ρ1=100Ω·m的均勻半空間，空氣視電阻率取值1012Ω·m。圖4為地壘起伏模型下，不同測點和不同頻率(所取頻率分別為0.01Hz，1Hz，100Hz)之間的模擬結果，其中，圖4(a)、圖4(b)、

圖1 視電阻率解析解與有限元計算結果對比Fig.1 Contrast of apparent resistivity between numerical and analytical results by FEM

圖2 阻抗相位解析解與有限元計算結果對比Fig.2 Contrast of impedance phase between numerical and analytical results by FEM

圖3 地壘起伏地形Fig.3 Horst topography

圖4(c)、圖4(d)分別為TM和TE極化模式下視電阻率及阻抗相位圖。從圖4可以看到，無論是TM模式還是TE模式，地壘起伏地形使得視電阻率和阻抗相位都產生了畸變。總體來說，TM模式相比TE模式，受地形影響更為嚴重，其中，TM極化模式的視電阻率值因地壘地形而減小，且與頻率成反比，而對TE模式下MT響應則相反，即使其視電阻率值增大，且頻率越高，響應越嚴重，而對于阻抗相位，兩種模式下的影響差別不大，頻率越高，影響越大。

2.3 地塹起伏模型

如圖5所示，地塹起伏地形上、下寬度分別為1 500m和500m，高差為500m，研究區域為地形線以下的均勻半空間，視電阻率ρ=100Ω·m。從圖6可以看到：TM極化模式下的視電阻率及阻抗相位的畸變遠比TE極化模式下的視電阻率與阻抗相位嚴重，其中，TM極化模式下的視電阻率因地塹地形的影響而增高，TE則與之相反；隨著頻率的增高，TM極化模式視電阻率影響較小，而TE模式則較高；對于阻抗相位，則影響基本一致，高頻影響較大，低頻影響較小。

圖4 地壘地形MT響應Fig.4 MT response with horst topography(a)、(b) 分別為TM視電阻率、相位曲線；(c)、(d) 分別為TE視電阻率、相位曲線

圖5 地塹起伏地形Fig.5 Graben topograph

2.4 地壘和地塹起伏異常體模型

如圖7及圖8所示，分別為地壘和地塹起伏異常體組合模型。地形起伏參數如圖所示，研究區域為地形線以下的均勻半空間，其視電阻率為ρ=100Ω·m。在地壘和地塹正下方距地表2 000m處有一低阻體(ρ=10Ω·m)，長和寬分別為1 000m。圖9和圖10 分別給出了地壘和地塹起伏異常體組合模型下的TM模式及TE模式下的視電阻率-頻率和阻抗相位-頻率擬斷面圖。從圖9和圖10中可以看到,無論是TM、TE模式，MT對地下異常體的探測都受地型起伏的影響。其中，TM極化模式的視電阻率和阻抗相位的擬斷面圖受地形影響較大，完全反映不出低阻異常體的特征；而TE極化模式的視電阻率和阻抗相位的擬斷面圖受地形影響較小 ，基本能反應出低阻異常體特征，高頻反映出了地形的影響規律，對于低頻而言，則幾乎沒有影響，與前面討論分析得出的結論保持一致。

圖6 地塹地形MT響應Fig.6 MT response with graben topography(a)、(b) 分別為TM視電阻率、相位曲線；(c)、(d) 分別為TE視電阻率、相位曲線

圖7 地壘起伏異常體組合模型Fig.7 Horst abnormal body composition model

圖8 地壘起伏異常體組合模型Fig.8 Horst abnormal body composition model

圖9 地壘異常體組合地形MT響應Fig.9 MT response with terrain abnormal topography(a)、(b) 分別為TM視電阻率、相位擬斷面圖；(c)、(d) 分別為TE視電阻率、相位擬斷面圖

圖10 地塹異常體組合地形MT響應Fig.10 MT response with graben abnormal topography(a)、(b) 分別為TM視電阻率、相位擬斷面圖；(c)、(d) 分別為TE視電阻率、相位擬斷面圖

2.5 地型起伏異常體組合模型

如圖11所示，為地壘、地塹起伏異常體組合模型。地形起伏參數如圖所示，研究區域為地形線以下的均勻半空間，視電阻率ρ=100Ω·m。在地形線中間下方距地表750m處有一低阻體(ρ=10Ω·m)，長和寬分別為1 000m。從圖12中可以看到，當異常體不位于地塹、地壘正下方時，無論是TM、TE模式，MT對地下異常探測均能一定程度的反映出來，其中，TE極化模式下的低阻異常體特征反映較TM效果更好，故在山區做MT時，應盡量選擇TE模式。

圖11 地壘、地塹起伏異常體組合模型Fig.11 Horst and graben abnormal body composition model

3 結論

我們采用矩形單元進一步細分四個小三角形的剖分方法，與不完全LU預處理BICGSTAB算法結合起來研究地形對MT正演數據的影響，結果表明：

圖12 地壘、地塹異常體組合地形MT響應 Fig.12 MT response with terrain and graben abnormal topography(a)TM模式下地壘；(b) TE模式下地壘； (c)地塹異常體組合視電阻率-頻率；(d) 阻抗相位-頻率擬斷面圖

1)該方法能很好地模擬各種起伏地形，驗證了該方法的實用性及準確性，為后期的反演提供借鑒。

2)TM極化模式下的視電阻率及阻抗相位受地形的影響遠遠大于TE極化模式，且頻率越低，影響越大，而TE則與之相反。

3)當異常體位于起伏地形正下方時，TE模式下的視電阻率及阻抗相位相比TM模式受地形影響較小，且在高頻時，地形響應的規律與上述結論基本一致，能明顯反映出異常體的基本特征；而當異常體不位于地形正下方時，則均能反應地下異常體特征。

4)在開展大地電磁測深工作時，應盡量選擇平坦地區，若無法避免起伏山區，則應選擇TE模式下測出的數據作為資料解釋參考。

[1] 程志平. 電法勘探教程[M]. 北京：冶金工業出版社，2012.CHENGZP.Electricalprospectingtutorial[M].Beijing:MetallurgicalIndustrypress, 2012. (InChinese)

[2]WANNAMAKERPE,STODTJA.Astablefiniteelementsolutionfortwo-dimensionalmagnetotelluricmodeling[J].Geophysics, 1987, 88(1): 277-296.

[3]FRANKEA,RALPHB. 2Dfiniteelementmodelingofplane-wavediffusivetime-harmonicelectromagneticfieldsusingadaptiveunstructuredgrids[J].ProceedingoftheWorkshop, 2004(S2): 1-6.

[4]MYUNGJN,HEEJK. 3Dmagnetotelluricmodelingincludingsurfacetopography[J].GeophysicalProspecting, 2007,55(2): 277-287.

[5] 徐世浙，趙生凱. 地形對大地電磁勘探的影響[J]. 西北地震學報，1985，7(4)：69-78.XUSZ,ZHAOSK.Thetopographiceffectsonmagnetotelluricresponse[J].NorthwesternSeismologicalJournal, 1985, 7(4): 69-78. (InChinese)

[6] 徐世浙，李予國，劉斌. 大地電磁型波二維地形改正的方法與效果[J]. 地球物理學報，1997，40(6)：842-846.XUSZ,LIYG,LIUB.Themethodandefficiencyof2-Dterraincorrectionfor-polarizationofMT[J].ChineseJournalOfGeophysics, 1997, 40(6): 842-846. (InChinese)

[7] 陳小斌.MT二維正演計算中地形影響的研究[J]. 石油物探，2000,39(3)：112-120.CHENXB.OntheresearchoftheinfluenceofterraintoMT2Dforwardcomputation[J].GeophysicalProspectingForPetroleum, 2000, 39(3): 112-120. (InChinese)

[8] 劉小軍，王家林，于鵬. 基于二次場的二維大地電磁有限元法數值模擬[J]. 同濟大學學報，2007，35(8)：1113-1117.LIUXJ,WANGJL,YUP.Secondaryfield-basedtwo-dimensionalmagnetotelluricnumericalmodelingbyfiniteelementmethod[J].JournalofTongjiUniversity, 2007, 35(8): 1113-1117. (InChinese)

[9] 柳建新，孫麗影，童孝忠. 基于自適應有限元的起伏地形MT二維正演模擬[J]. 地球物理學進展，2012，27(5)：2016-2023.LIUJX,SUNLY,TONGXZ.Undulatingterrain2DMTforwardmodelingwithadaptivefiniteelement[J].ProgressInGeophysics, 2012, 27(5): 2016-2023. (InChinese)

[10]陳樂壽，王光鍔. 大地電磁測深法[M]. 北京：地質出版社，1990.CHENLS,WANGGE.Magnetotelluricsoundingmethod[M].Beijing:SciencePress, 1990. (InChinese)

[11]徐世浙. 地球物理中的有限單元法[M]. 北京：地質出版社，1994.XUSZ.Finiteelementmethodforgeophysics[M].Beijing:SciencePress, 1994. (InChinese)

[12]柳建新，蔣騰飛，童孝忠. 不完全LU分解預處理的BICGSTAB算法在大地電磁二維正演模擬中的應用[J]. 中南大學學報，2009，40(2)：485-491.LIUJX,JIANGTF,TONGXZ.ApplicationofBICGSTABalgorithmwithincompleteLUdecompositionpreconditioningtotwo-dimensionalmagnetotelluricforwardmodeling[J].JournalOfCentralSouthUniversity, 2009, 40(2): 485-491. (InChinese)

2D MT forward modeling with topography using the ILU-BICGSTAB algorithm

GUO Shengnan; XIONG Bin; LUO Tianya; WU Yanqiang; CHEN Xin

(Guilin University of Technology, College of Earth Science, Guilin 541006, China)

In allusion to the issue that complex topography effects on magnetotelluric sounding, using finite element method to compute 2D magnetotelluric responses on MATLAB platform. We used rectangular grid further refines the triangulation mesh mode in order to better simulate the actual topography. The BICGSTAB algorithm with incomplete LU decomposition for preconditioning was used to improve calculating speed and precision. This approach was verified through the calculations of various geoelectric model. Finally, we analysis topographic responses in magnetotelluric.

magnetotelluric； finite element method； topography； forward

2016-05-26 改回日期：2017-01-04

國家自然科學基金(41164004,41674075);廣西自然科學基金(2016GXNSFA038004,2013GXNSFAA019277);廣西壯族自治區桂林市“漓江學者”(2013)專項;桂林理工大學研究生創新項目(BS201601)

郭勝男(1990-)，男，碩士，主要從事電磁法數值模擬研究，E-mail：1179251693@qq.com。

熊彬(1974-)，男，博士，教授，主要從事電磁場理論及反演成像方面的教學與科研工作， E-mail：xiongbin@msn.com。

1001-1749(2017)03-0306-07

P 631.4

A

10.3969/j.issn.1001-1749.2017.03.02