韓騎， 胡祥云*， 程正璞， 楊炳南， 蔡建超， 韋偉

1 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 4300742 貴州省地礦局， 貴州銅仁 554300

?

自適應(yīng)非結(jié)構(gòu)有限元MT二維起伏地形正反演研究

韓騎1， 胡祥云1*， 程正璞1， 楊炳南2， 蔡建超1， 韋偉1

1 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 4300742 貴州省地礦局， 貴州銅仁 554300

在山區(qū)進(jìn)行MT勘探時(shí)，用規(guī)則網(wǎng)格有限元方法模擬起伏地形會(huì)受到限制.本文采用非結(jié)構(gòu)三角網(wǎng)格可以有效地模擬任意二維地質(zhì)結(jié)構(gòu)，如起伏地形、傾斜巖層和多尺度構(gòu)造等.正演引入自適應(yīng)有限元方法，其在網(wǎng)格剖分過(guò)程中能根據(jù)單元誤差自動(dòng)細(xì)化網(wǎng)格，保證了正演結(jié)果的精度.將自適應(yīng)有限元與Occam算法結(jié)合，且引用并行處理技術(shù)提高正反演計(jì)算速度.通過(guò)對(duì)比兩個(gè)理論模型，討論了地形對(duì)MT正演響應(yīng)的影響；其次進(jìn)行了不同地電模型帶地形反演展示了本文算法的正確性和適用性；最后將該方法應(yīng)用于實(shí)測(cè)MT數(shù)據(jù)處理，證明了自適應(yīng)非結(jié)構(gòu)有限元方法是復(fù)雜地形下處理MT數(shù)據(jù)的有力工具.

大地電磁； 自適應(yīng)有限元； 地形； 反演

1 引言

大地電磁測(cè)深法（MT）是通過(guò)天然交變電磁場(chǎng)研究地球模型電性結(jié)構(gòu)的一種地球物理方法.大地電磁測(cè)深方法具有較大的穿透深度和分辨能力，它不僅是研究地殼與上地幔深部構(gòu)造的強(qiáng)有力工具（李冉等，2014；呂慶田等，2011），而且可以結(jié)合地震數(shù)據(jù)（劉彥等，2012；詹艷等，2014）和重磁數(shù)據(jù)（張交東等，2012）進(jìn)行綜合解釋分析，已成為礦產(chǎn)勘查（Hu et al.，2013；陸桂福和吳新剛，2014）、地?zé)豳Y源（Wu et al.，2012；高景宏等，2010）和油氣勘探（夏訓(xùn)銀等，2012）的有效方法.

地形對(duì)大地電磁測(cè)深的影響是地球物理工作者難以避免的問(wèn)題.特別是我國(guó)山地區(qū)域地形起伏劇烈，大地電磁數(shù)據(jù)采集受各種因素的限制，當(dāng)測(cè)點(diǎn)布置在山溝時(shí)，垂直山體方向受高電阻體對(duì)電場(chǎng)的排斥作用的影響；若山體不規(guī)則，地形影響更為復(fù)雜（張翔等，1999）.有限元方法是模擬MT正演響應(yīng)的有效手段，許多研究者利用矩形網(wǎng)格（Rodi，1976；陳樂(lè)壽，1981；童孝忠等，2009；張昆等，2008）和矩形對(duì)分三角網(wǎng)格（Wannamaker et al.，1987；Zyserman et al.，1999；劉小軍等，2007；柳建新等，2009）模擬地形正演響應(yīng)，結(jié)果表明地形對(duì)TM模式的影響遠(yuǎn)遠(yuǎn)大于TE模式（王緒本等，1999），并且視電阻率曲線(xiàn)表現(xiàn)為畸變、靜態(tài)平移或者二者的混合，沒(méi)有一定規(guī)律可循.但矩形及其對(duì)分三角網(wǎng)格對(duì)復(fù)雜地質(zhì)結(jié)構(gòu)（起伏地形和交錯(cuò)地層等）模擬能力有限，且上述規(guī)則網(wǎng)格在對(duì)目標(biāo)區(qū)域進(jìn)行加密時(shí)，其背景區(qū)域也會(huì)受到不必要細(xì)化，增加了計(jì)算時(shí)間.

近年來(lái)許多學(xué)者（Key and Weiss，2006；Franke et al.，2007；Li and Key，2007；Li and Pek，2008）采用自適應(yīng)非結(jié)構(gòu)三角網(wǎng)格有限元方法計(jì)算電磁場(chǎng).非結(jié)構(gòu)三角網(wǎng)格能夠真實(shí)地模擬各種復(fù)雜地質(zhì)結(jié)構(gòu)，該方法利用Delaunay三角劃分原理（Shewchuk，2002）對(duì)目標(biāo)區(qū)域進(jìn)行細(xì)密網(wǎng)格剖分，而背景區(qū)域則粗糙網(wǎng)格剖分.自適應(yīng)有限元對(duì)計(jì)算結(jié)果進(jìn)行誤差評(píng)估（Ovall，2006），誤差較大的單元網(wǎng)格自動(dòng)細(xì)分直至有限元計(jì)算結(jié)果滿(mǎn)足給定的誤差范圍.自適應(yīng)非結(jié)構(gòu)有限元方法因其強(qiáng)大的適應(yīng)性，越來(lái)越受到地球物理研究者的青睞，其方法技術(shù)也日臻完善.現(xiàn)已成功實(shí)現(xiàn)了海洋二維MT正演模擬（Key and Weiss，2006；Franke et al.，2007；趙慧等，2014），海洋2.5D可控源正演模擬（Key and Ovall，2011）.在此基礎(chǔ)上，本文采用非結(jié)構(gòu)三角網(wǎng)格剖分、自適應(yīng)有限元與Occam算法結(jié)合的策略來(lái)研究陸地MT二維帶地形反演.首先討論了不同地形對(duì)MT正演的影響；其次建立3個(gè)典型地電模型進(jìn)行帶地形反演，并與平地反演結(jié)果進(jìn)行對(duì)比；最后對(duì)實(shí)測(cè)MT數(shù)據(jù)進(jìn)行反演并解釋.

2 方法原理

（1）

（2）

事先給出一次場(chǎng)對(duì)應(yīng)的電導(dǎo)率為σ0，一次場(chǎng)為均勻半空間模型其電磁場(chǎng)值由解析解算出，二次場(chǎng)則是由σ-σ0的差異引起，采用自適應(yīng)有限元解出場(chǎng)值，最后將一次場(chǎng)和二次場(chǎng)相加算出總場(chǎng)值.該方法可通過(guò)提高二次場(chǎng)解的精度來(lái)提高整體解的精度（Rodi，1976），與總場(chǎng)法求解的結(jié)果對(duì)比驗(yàn)證了該方法的正確性（趙廣茂等，2008；劉小軍等，2007）.這種方法最大優(yōu)勢(shì)在于有限元網(wǎng)格只需要精確地模擬二次場(chǎng)的變化，因此可采用包含更少單元的網(wǎng)格（Key and Ovall，2011），提高計(jì)算速度.利用Helmholtz方程得到二次場(chǎng)TE和TM偏微分方程表達(dá)式如下：

（3）

（4）

式中ω為角頻率，μ=4π×10-7H·m-1為磁導(dǎo)率常數(shù).對(duì)研究區(qū)域進(jìn)行三角剖分后生成包含N個(gè)頂點(diǎn)的N維向量空間V，式（3）和（4）對(duì)應(yīng)的變分問(wèn)題可以統(tǒng)一寫(xiě)成：

B（v,u）=F（v）,

（5）

（6）

（7）

二次場(chǎng)與一次場(chǎng)相加求得總場(chǎng)值，則MT阻抗為

（8）

MT視電阻率和相位角為

（9）

相對(duì)于有限元的解un，解的誤差G（u-un）更被研究者重視，其中u是（5）式的真實(shí)解，G為誤差評(píng)估的對(duì)偶函數(shù).自適應(yīng)網(wǎng)格劃分的實(shí)質(zhì)是計(jì)算出每個(gè)三角單元有限元解的誤差，對(duì)誤差較大的三角單元再次細(xì)分，直到誤差達(dá)到允許范圍.本文采用對(duì)偶加權(quán)殘差（DWR）法（Ovall，2006）計(jì)算單元誤差，屬于一種后驗(yàn)誤差處理技術(shù)，通過(guò)分級(jí)函數(shù)對(duì)殘差進(jìn)行加權(quán).選擇一個(gè)合適G函數(shù)十分重要，DWR給出了一個(gè)比較合適的G函數(shù)，單元誤差計(jì)算過(guò)程中多次利用有限元數(shù)值方法進(jìn)行求解，所求誤差是逼近真實(shí)解的近似值.后驗(yàn)誤差方法在工程學(xué)領(lǐng)域是一項(xiàng)非常成熟的技術(shù)，近幾十年被引入到電磁場(chǎng)模擬中，已應(yīng)用于海洋可控源電磁場(chǎng)（Li and Key，2007；Key and Ovall，2011）、海洋二維MT電磁場(chǎng) （Franke et al，2007；Key and Weiss，2006；趙慧等，2014）和二維電阻率各向異性介質(zhì)大地電磁場(chǎng)（Li and Pet，2008）模擬.

反演問(wèn)題采用Occam算法（Constable et al.，1987），該算法自發(fā)表以來(lái)由于其穩(wěn)定性在電磁法反演領(lǐng)域得到了廣泛的應(yīng)用（Ghaedrahmati et al.，2014；Siripunvaraporn and Sarakorn，2011；Sudha et al.，2014；何梅興等，2011）.Occam反演實(shí)質(zhì)為帶光滑約束的最小二乘法，其泛函U表達(dá)式為：

（10）

自適應(yīng)有限元方法進(jìn)行電磁場(chǎng)數(shù)值模擬的主要優(yōu)勢(shì)在于可根據(jù)每次有限元計(jì)算的結(jié)果自動(dòng)進(jìn)行網(wǎng)格精煉，確保了計(jì)算結(jié)果的精度，避免了以往正反演過(guò)程中網(wǎng)格繪制的繁瑣工作，圖1給出了自適應(yīng)有限元與Occam算法結(jié)合的計(jì)算流程.首先給出反演初始模型作為反演細(xì)化的初始網(wǎng)格和電阻率模型參數(shù)，一般來(lái)說(shuō)，目標(biāo)區(qū)域網(wǎng)格劃分細(xì)致，選擇小三角單元；背景區(qū)域則采用大三角單元.隨后進(jìn)行內(nèi)存空間分配，即采用多少線(xiàn)程參與計(jì)算，每個(gè)線(xiàn)程負(fù)責(zé)多少測(cè)點(diǎn)的自適應(yīng)有限元運(yùn)算.本文設(shè)置單元誤差閾值0.02，采用有限元方法計(jì)算每個(gè)單元的場(chǎng)值，由DWR方法估計(jì)單元局部誤差，將局部誤差大于誤差閾值的單元進(jìn)行細(xì)化，得到新的三角網(wǎng)格.重復(fù)以上過(guò)程直至所有單元誤差達(dá)到誤差閾值或者已達(dá)到最大網(wǎng)格細(xì)化次數(shù)，則網(wǎng)格精煉停止，輸出正演結(jié)果.設(shè)置反演收斂精度1.5，計(jì)算數(shù)據(jù)均方根誤差RMS，若未達(dá)到反演收斂精度，則更新模型重新開(kāi)始自適應(yīng)有限元計(jì)算過(guò)程直至達(dá)到反演收斂精度或者達(dá)到最大迭代次數(shù)，輸出反演結(jié)果.

圖1 自適應(yīng)有限元算法反演流程

3 模型計(jì)算

3.1 地形對(duì)MT響應(yīng)影響研究

為了研究地形對(duì)MT勘探的影響，在100 Ωm的均勻半空間中設(shè)置一10 Ωm的矩形低阻體，在圖2a中設(shè)置包含山谷、山峰起伏地形，最大高差300 m，圖2b采用水平地形.從左往右依次設(shè)置11個(gè)測(cè)點(diǎn)，間距300 m，頻率范圍為0.1～10000 Hz，共22個(gè)頻點(diǎn).

圖2 地形影響研究模型

采用自適應(yīng)有限元方法計(jì)算各測(cè)點(diǎn)的視電阻率和相位響應(yīng)值，選取代表典型地形的三個(gè)測(cè)點(diǎn)MT04、MT05和MT06正演結(jié)果進(jìn)行比較.如圖3所示，紅色曲線(xiàn)為T(mén)M模式，藍(lán)色曲線(xiàn)為T(mén)E模式.圖3a中TE模式視電阻率和相位曲線(xiàn)相對(duì)于平地模型（圖3b）波動(dòng)不大；但是對(duì)于TM模式，在山峰處（MT04）視電阻率曲線(xiàn)低頻陡然下降，山坡處（MT05）視電阻率曲線(xiàn)整體向上平移，且在高頻處上揚(yáng)，山谷處（MT06）視電阻率曲線(xiàn)高頻處上揚(yáng).不同地形對(duì)正演響應(yīng)產(chǎn)生的干擾不同，視電阻率曲線(xiàn)會(huì)產(chǎn)生畸變、平移或二者的結(jié)合，沒(méi)有十分確切的規(guī)律可循，整體來(lái)說(shuō)地形對(duì)視電阻率的影響大于相位角，對(duì)TM模式的影響遠(yuǎn)遠(yuǎn)大于TE模式，推測(cè)其原因是兩種極化模式中，磁場(chǎng)一般變化平緩，TE模式電場(chǎng)平行于構(gòu)造走向，TM模式電場(chǎng)垂直走向，垂直于走向的電場(chǎng)會(huì)產(chǎn)生感應(yīng)電場(chǎng)，該感應(yīng)電場(chǎng)在空間上的分布受地形起伏的影響嚴(yán)重，而在TE模式中只有與走向平行的電場(chǎng).因此，實(shí)現(xiàn)帶地形反演消除因地形產(chǎn)生的數(shù)據(jù)異常十分重要.

3.2 典型地電模型正反演研究

在山區(qū)進(jìn)行MT工作主要以尋找金屬礦為目的，金屬礦一般呈現(xiàn)低阻性質(zhì).如圖4所示，在100 Ωm均勻半空間埋藏兩個(gè)形態(tài)不一，埋深不同的10 Ωm的低阻體，1號(hào)低阻體中心埋深400 m，呈塊狀；2號(hào)低阻體中心埋深750 m，呈薄板狀.測(cè)線(xiàn)垂直走向方向分布，全長(zhǎng)2000 m，其中包含40個(gè)測(cè)點(diǎn)，頻率范圍0.1～10000 Hz，共60個(gè)頻點(diǎn)，地形采用一組實(shí)測(cè)高程數(shù)據(jù)，最大高差達(dá)250 m.

圖3 MT04，MT05，MT06測(cè)點(diǎn)視電阻率和相位響應(yīng)

圖4 金屬礦模型（其中1、2為低阻體編號(hào)）

以圖4模型進(jìn)行網(wǎng)格剖分，初始網(wǎng)格劃分為2035個(gè)三角單元（圖5a），對(duì)單元誤差取對(duì)數(shù)值顯示，其分布如圖5b所示.采用四線(xiàn)程對(duì)初始網(wǎng)格細(xì)化，以其中一個(gè)線(xiàn)程為例，該線(xiàn)程負(fù)責(zé)最后10個(gè)測(cè)點(diǎn)的自適應(yīng)正演運(yùn)算，經(jīng)過(guò)8次細(xì)化，生成18854個(gè)三角單元（圖5c），單元誤差達(dá)到誤差閾值以?xún)?nèi)（圖5d）則停止計(jì)算.網(wǎng)格主要在電性分界處即低阻體邊界和測(cè)點(diǎn)附近進(jìn)行細(xì)化，圖5e顯示測(cè)點(diǎn)附近網(wǎng)格細(xì)化信息.其余三個(gè)線(xiàn)程同時(shí)進(jìn)行自適應(yīng)細(xì)化過(guò)程，共耗時(shí)123 s.

利用上述模型的正演數(shù)據(jù)添加4%的隨機(jī)噪聲進(jìn)行Occam反演，對(duì)于目標(biāo)區(qū)域采用小三角單元細(xì)致剖分，周?chē)尘皡^(qū)域則粗糙剖分.反演初始模型如圖6，其中包含5958個(gè)電阻率為100 Ωm三角單元.

帶地形反演結(jié)果（圖7a）顯示1號(hào)和2號(hào)低阻體中心位置，可圈定出1號(hào)低阻體的輪廓，2號(hào)低阻體邊界較模糊.圖7b為相同數(shù)據(jù)平地反演結(jié)果，相對(duì)于帶地形反演，低阻體出現(xiàn)中心位置偏移及范圍擴(kuò)大的情況，且在近地表處存在虛假異常也稱(chēng)為靜態(tài)效應(yīng).地形是引起靜態(tài)效應(yīng)的原因之一，常規(guī)處理中，反演前會(huì)進(jìn)行靜態(tài)效應(yīng)校正，用以消除淺層假異?，F(xiàn)象，常用方法為整體平移數(shù)據(jù)曲線(xiàn)，這需要豐富的經(jīng)驗(yàn)且平移結(jié)果不一定準(zhǔn)確，帶地形反演不能完全消除但可以有效地抑制靜態(tài)效應(yīng).

為了驗(yàn)證該算法對(duì)淺部薄層低阻礦床的反演效果，我們建立圖8a所示的地電模型，模型分為兩層，上層電阻率為100 Ωm，層厚約870 m；下層為1000 Ωm，上層存在10 Ωm的低阻體，中心埋深450 m.測(cè)點(diǎn)、地形、頻點(diǎn)設(shè)置與圖4模型一致.反演結(jié)果（圖8b）較好地顯示薄層低阻體的中心位置，可識(shí)別下層高阻帶，且地層分界明顯.

在我國(guó)地形起伏的山區(qū)中，構(gòu)造發(fā)育復(fù)雜，常見(jiàn)褶皺、斷層分布.建立相應(yīng)的綜合褶皺模型（圖9a）進(jìn)行反演，其中以向斜為主，連接2個(gè)小型背斜，10和100 Ωm地層互層分布.測(cè)點(diǎn)、地形、頻點(diǎn)設(shè)置與圖4模型一致.反演結(jié)果（圖9b）顯示該褶皺的大致輪廓，地層分層清晰，可看出中部向斜構(gòu)造.

表1給出了上述三個(gè)模型反演耗時(shí)，RMS誤差及粗糙度.由于本文模型運(yùn)算采用四線(xiàn)程，模型每更新一次，便重新開(kāi)始自適應(yīng)細(xì)化，導(dǎo)致反演時(shí)間大大增加.在實(shí)際應(yīng)用中，處理大量的數(shù)據(jù)需要更多的時(shí)間，如何優(yōu)化算法是亟待解決的問(wèn)題.

表1 反演耗時(shí)、誤差及粗糙度列表

4 實(shí)測(cè)資料解釋

工區(qū)位于中國(guó)西南山地礦區(qū)，區(qū)內(nèi)地形起伏劇烈.該區(qū)地表水系發(fā)育，礦區(qū)出露地層比較單一，主要為中三疊統(tǒng)許滿(mǎn)組及第四系，但構(gòu)造發(fā)育復(fù)雜，主要為褶皺、斷裂和節(jié)理.數(shù)據(jù)采集使用加拿大phoenix公司V8多功能電法儀，通過(guò)不極化電極測(cè)量Ex、Ey兩道正交的電場(chǎng)水平分量數(shù)據(jù)，磁傳感器測(cè)量Hx、Hy、Hz三道正交的磁場(chǎng)分量數(shù)據(jù)，從而實(shí)現(xiàn)五個(gè)電磁場(chǎng)分量的電磁測(cè)深數(shù)據(jù)采集.采集頻率范圍設(shè)置為10400～0.35 Hz，在低頻處截屏到1 Hz，探測(cè)深度約地下2 km以?xún)?nèi)電性結(jié)構(gòu).數(shù)據(jù)噪聲干擾不強(qiáng)，在高頻有一天然死頻帶，采用D+方法進(jìn)行處理.受地形影響，TM和TE曲線(xiàn)首端較多分離，擬斷面圖“掛面條”現(xiàn)象嚴(yán)重.

選擇工區(qū)中A3測(cè)線(xiàn)中部進(jìn)行反演.測(cè)線(xiàn)由西南向東北方向布設(shè)，最大高差242 m，地形起伏明顯.工區(qū)布置如圖10所示（中間測(cè)線(xiàn)為A3），圖左下角紅色點(diǎn)表示反演起始測(cè)點(diǎn)，右上角紅色點(diǎn)為終點(diǎn).A3測(cè)線(xiàn)所在位置出露地層為中三疊統(tǒng)許滿(mǎn)組第三段（T2xm3），其包含的主要巖礦石電性見(jiàn)表2.圖10顯示A3測(cè)線(xiàn)中部穿越一大型向斜構(gòu)造和一背斜構(gòu)造，測(cè)線(xiàn)末端發(fā)育一系列小型褶皺，與兩大褶皺主軸近似平行，測(cè)線(xiàn)中部和末端分別穿越小型河流.

圖6 反演初始模型

圖7 金屬礦模型反演結(jié)果（其中1、2是低阻體編號(hào)）

圖8 薄層低阻體模型

圖9 綜合褶皺模型

圖10 工區(qū)地形地質(zhì)及工程布置圖

圖11 A3測(cè)線(xiàn)反演結(jié)果

圖12 錳礦實(shí)測(cè)MT數(shù)據(jù)反演結(jié)果

圖13 電阻率測(cè)深曲線(xiàn)

A3測(cè)線(xiàn)中部共40個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)間距約40 m，全長(zhǎng)1600 m，頻率范圍1～10400 Hz，共53個(gè)頻點(diǎn)，反演共耗時(shí)31230 s，RMS為3.74.結(jié)果如圖11所示，大致可劃分為三層電性層，黑色虛線(xiàn)表示地層線(xiàn)，淺地表處為低阻層，隨深度增加電阻率逐漸上升，變化范圍10～1000 Ωm，與表2所列T2xm3地層電阻率變化情況一致；反演結(jié)果顯示背斜和向斜構(gòu)造，背斜上部被剝蝕，整體電性表現(xiàn)為高阻，在其內(nèi)部包含低阻體.工區(qū)地質(zhì)資料顯示該處斷層、裂隙、水系發(fā)育豐富，推測(cè)該低阻體成因可能是地表水沿裂隙滲漏使該處呈現(xiàn)低阻性質(zhì).

表2 測(cè)區(qū)主要巖礦石電阻率統(tǒng)計(jì)表

另一組數(shù)據(jù)為尋找錳礦實(shí)測(cè)電磁數(shù)據(jù)，錳礦成礦地質(zhì)體主要由炭質(zhì)、粉砂質(zhì)粘土巖夾條帶狀、塊狀菱錳礦組成，其綜合電性特征表現(xiàn)為低電阻率特征，其上覆地層為含礫砂巖、細(xì)晶白云巖等高電阻率巖性.

測(cè)線(xiàn)共25個(gè)測(cè)點(diǎn)，全長(zhǎng)1200 m，0.35～10400 Hz共60個(gè)頻點(diǎn).其反演結(jié)果如圖12所示，地層分層明顯，淺部地層為高阻，厚度約200～300 m，隨深度增加電阻率逐漸減小，在200 m深度處電阻率有增大趨勢(shì)，之后又逐漸減小.該測(cè)線(xiàn)處分布有三口鉆孔，如圖中黑色圓柱所示，分別位于測(cè)線(xiàn)前段（ZK01）、中段（ZK02）和后段（ZK03）.ZK01見(jiàn)礦深度214 m；ZK02見(jiàn)礦深度261 m；ZK03見(jiàn)礦深度501 m.反演結(jié)果中鉆孔處視電阻率隨深度變化曲線(xiàn)如圖13，括號(hào)內(nèi)的數(shù)字代表鉆孔所在水平位置，箭頭指向見(jiàn)礦深度.三條測(cè)深曲線(xiàn)變化趨勢(shì)一致，電阻率隨深度迅速減小后又緩慢回升，見(jiàn)礦深度并不在視電阻率最小處，而是在視電阻率緩慢回升帶，結(jié)合反演結(jié)果和鉆孔資料用紅色虛線(xiàn)在圖12中圈定出成礦構(gòu)造帶大致范圍.

5 結(jié)論

本文利用非結(jié)構(gòu)自適應(yīng)有限元算法研究地形對(duì)MT正演數(shù)據(jù)的影響，該方法能夠真實(shí)地模擬起伏地形、不規(guī)則異常體及復(fù)雜形態(tài)的地質(zhì)構(gòu)造.研究結(jié)果表明地形對(duì)TE模式影響很小，對(duì)TM模式影響遠(yuǎn)大于TE模式，與平地視電阻率正演曲線(xiàn)相比，TM模式在山峰處低頻部分下降，斜坡處曲線(xiàn)整體向上平移，在山谷處高頻部分上揚(yáng).

將自適應(yīng)有限元與Occam算法結(jié)合進(jìn)行帶地形反演運(yùn)算，通過(guò)理論模型試算和分析，結(jié)果證明了該算法能清晰識(shí)別出異常體的中心位置且分層性良好，與平地反演結(jié)果相比，不會(huì)在淺部產(chǎn)生虛假異常，能有效抑制靜態(tài)效應(yīng).應(yīng)用于電磁實(shí)測(cè)數(shù)據(jù)反演結(jié)果顯示了地層分層信息、背斜構(gòu)造、見(jiàn)礦層位，探測(cè)的地質(zhì)結(jié)構(gòu)和圈定的成礦構(gòu)造帶均與已有地質(zhì)資料相符合，表明該算法具有較強(qiáng)的實(shí)用性，可用于實(shí)測(cè)資料解釋分析.

在此基礎(chǔ)上，將繼續(xù)研究帶地形反演條件下如何正確進(jìn)行靜態(tài)校正，優(yōu)化算法縮短反演計(jì)算時(shí)間，并根據(jù)其他已知信息進(jìn)行模型約束反演和聯(lián)合反演，最終將非結(jié)構(gòu)網(wǎng)格引入起伏地形下大地電磁三維正反演研究中.

致謝 感謝中國(guó)地質(zhì)科學(xué)院呂慶田研究員和中南大學(xué)湯井田教授對(duì)文章的悉心指導(dǎo)，對(duì)手稿提出的意見(jiàn)豐富了文章的內(nèi)容.

Chen L S. 1981. Application and improvement of finite element method in forward calculation of geo-electromagnetic field.GeophysicalProspectingforPetrole（in Chinese）, （3）: 84-103.

Coggon J H. 1971. Electromagnetic and electrical modeling by the finite element method.Geophysics, 36（1）: 132-155.

Constable S C, Parker R L, Constable C G. 1987. Occam′s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data.Geophysics, 52（3）: 289-300.

Franke A, B?rner R U, Spitzer K. 2007. Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography.GeophysicalJournalInternational, 171（1）: 71-86.

Gao J H, Tong T G, Qiang J K, et al. 2010. A magnetotelluric study of geothermal resources in Kaifeng depression, He′nan Province.GeophysicalandGeochemicalExploration（in Chinese）, 34（4）: 440-443.

Ghaedrahmati R, Moradzadeh A, Fathianpour N, et al. 2014. Investigating 2-D MT inversion codes using real field data.ArabianJournalofGeosciences, 7（6）: 2315-2328.

He M X, Hu X Y, Ye Y X, et al. 2011. 2.5D controlled source audio-frequency magnetotellurics Occam inversion.ProgressinGeophysics（in Chinese）, 26（6）: 2163-2170, doi: 10.3969/j.issn.1004-2903.2011.06.033.

Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China.Geophysics, 78（3）: B111-B119.

Key K, Weiss C. 2006. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example.Geophysics, 71（6）: G291-G299.Key K, Ovall J. 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling.GeophysicalJournalInternational, 186（1）: 137-154.

Li R, Tang J, Dong Z Y, et al. 2014. Deep electrical conductivity structure of the southern area in Yunnan Province.ChineseJournalofGeophysics（in Chinese）, 57（4）: 1111-1122, doi: 10.6038/cjg20140409.

Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modeling: Part 1—An adaptive finite-element algorithm.Geophysics, 72（2）: WA51-WA62.Li Y G, Pek J. 2008. Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media.GeophysicalJournalInternational, 175（3）: 942-954.

Liu J X, He H, Liu H F, et al. 2009. 2D inversion of VES data under rolling topography and its application.ComputingTechniquesforGeophysicalandGeochemicalExploration（in Chinese）, 31（4）: 293-296.

Liu X J, Wang J L, Yu P. 2007. Secondary field-based two-dimensional magnetotelluric numerical modeling by finite element method.JournalofTongjiUniversity（NaturalScience） （in Chinese）, 35（8）: 1113-1117.

Liu Y, Lü Q T, Meng G X, et al. 2012. Joint electromagnetic and seismic inversion survey: status and prospect.ProgressinGeophysics（in Chinese）, 27（6）: 2444-2451, doi: 10.6038/j.issn.1004-2903.2012.06.019.

Lu G F, Wu X G. 2014. Application of comprehensive electrical prospecting method in the exploration of the concealed metallic deposit.MineralExploration（in Chinese）, 5（4）: 617-622.

Lü Q T, Shi D N, Tang J T, et al. 2011. Probing on deep structure of middle and lower reaches of the Yangtze metallogenic belt and typical ore concentration area: A review of annual progress of SinoProbe-03.ActaGeoscienticaSinica（in Chinese）, 32（3）: 257-268.

Ovall J S. 2006. Asymptotically exact functional error estimators based on superconvergent gradient recovery.NumerischeMathematik, 102（3）: 543-558.

Rodi W L. 1976. A technique for improving the accuracy of finite element solutions for magnetotelluric data.GeophysicalJournalInternational, 44（2）: 483-506.

Shewchuk J R. 2002. Delaunay refinement algorithms for triangular mesh generation.ComputationalGeometry, 22（1-3）: 21-74.

Siripunvaraporn W, Sarakorn W. 2011. An efficient data space conjugate gradient Occam′s method for three-dimensional magnetotelluric inversion.GeophysicalJournalInternational, 186（2）: 567-579.

Sudha, Tezkan B, Siemon B. 2014. Appraisal of a new 1D weighted joint inversion of ground based and helicopter-borne electromagnetic data.GeophysicalProspecting, 62（3）: 597-614.

Tong X Z, Liu J X, Guo R W. 2009. Solution strategies for complex 2D/3D magnetotelluric forward modeling based on the finite element method.ComputerizedTomographyTheoryandApplications（in Chinese）, 18（1）: 47-54.

Wang X B, Li Y N, Gao Y C. 1999. Two dimensional topographic responses in magneto-telluric sounding and its correction methods.ComputingTechniquesforGeophysicalandGeochemicalExploration（in Chinese）, 21（4）: 327-332.

Wannamaker P E, Stodt J A, Rijo L. 1987. A stable finite element solution for two-dimensional magnetotelluric modelling.GeophysicalJournalInternational, 88（1）: 277-296.

Wu G J, Hu X Y, Huo G P, et al. 2012. Geophysical exploration for geothermal resources: An application of MT and CSAMT in Jiangxia, Wuhan, China.JournalofEarthScience, 23（5）: 757-767.

Xia X Y, Wang J X, Liu J C, et al. 2012. The application of magnetotelluric sounding to lithofacies division in oil and gas exploration of northern nanpangjiang area.GeophysicalandGeochemicalExploration（in Chinese）, 36（2）: 174-179.

Zhan Y, Zhao G Z, Wang L F, et al. 2014. Deep electric structure beneath the intersection area of West Qinling orogenic zone with North-South Seismic tectonic zone in China.ChineseJ.Geophys. （in Chinese）, 57（8）: 2594-2607, doi: 10.6038/cjg20140819.

Zhang J D, Yang X Y, Liu C Z, et al. 2012. The fine deep structure of the northern margin of the Dabie Orogenic Belt from gravity-magnetic-electrical-seismic combination survey.ChineseJ.Geophys. （in Chinese）, 55（7）: 2292-2306.

Zhang K, Wei W B, Ye G F. 2008. The actualization of 2D finite-element magnetotelluric forwarding on Matlab.SeismologicalandGeomagneticObservationandResearch（in Chinese）, 29（5）: 83-88.

Zhang X, Hu W B, Yan L J, et al. 1999. Effects and correction of topography in magnetotelluric sounding.JournalofJianghanPetroleumInstitute（in Chinese）, 21（1）: 37-41.

Zhao G M, Li T L, Wang D Y, et al. 2008. Secondary field-based two-dimensional topographic numerical simulation in magnetotellurics by finite element method.JournalofJilinUniversity（EarthScienceEdition） （in Chinese）, 38（6）: 1055-1059.

Zhao H, Liu Y, Li Y G. 2014. Adaptive finite element forward modeling for two-dimensional marine magnetotelluric fields.OilGeophysicalProspecting（in Chinese）, 49（3）: 578-585.

Zyserman F I, Guarracino L, Santos J E. 1999. A hybridized mixed finite element domain decomposed method for two dimensional magnetotelluric modelling.Earth,PlanetsandSpace, 51（4）: 297-306.

附中文參考文獻(xiàn)

陳樂(lè)壽. 1981. 有限元法在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用及改進(jìn). 石油物探, （3）: 84-103.

高景宏, 佟鐵鋼, 強(qiáng)建科等. 2010. 開(kāi)封凹陷區(qū)地?zé)豳Y源大地電磁測(cè)深研究. 物探與化探, 34（4）: 440-443.

何梅興, 胡祥云, 葉益信等. 2011. 2.5維可控源音頻大地電磁法Occam反演理論及應(yīng)用. 地球物理學(xué)進(jìn)展, 26（6）: 2163-2170, doi: 10.3969/j.issn.1004-2903.2011.06.033.

李冉, 湯吉, 董澤義等. 2014. 云南南部地區(qū)深部電性結(jié)構(gòu)特征研究. 地球物理學(xué)報(bào), 57（4）: 1111-1122, doi: 10.6038/cjg20140409. 柳建新, 何歡, 劉海飛等. 2009. 起伏地形垂直電測(cè)深二維反演及應(yīng)用. 物探化探計(jì)算技術(shù), 31（4）: 293-296.

劉小軍, 王家林, 于鵬. 2007. 基于二次場(chǎng)的二維大地電磁有限元法數(shù)值模擬. 同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版）, 35（8）: 1113-1117.

劉彥, 呂慶田, 孟貴祥等. 2012. 大地電磁與地震聯(lián)合反演研究現(xiàn)狀與展望. 地球物理學(xué)進(jìn)展, 27（6）: 2444-2451, doi: 10.6038/j.issn.1004-2903.2012.06.019.

陸桂福, 吳新剛. 2014. 綜合電法勘查在隱伏金屬礦勘查中的應(yīng)用效果. 礦產(chǎn)勘查, 5（4）: 617-622.

呂慶田, 史大年, 湯井田等. 2011. 長(zhǎng)江中下游成礦帶及典型礦集區(qū)深部結(jié)構(gòu)探測(cè)——SinoProbe-03年度進(jìn)展綜述. 地球?qū)W報(bào), 32（3）: 257-268.

童孝忠, 柳建新, 郭榮文. 2009. 復(fù)雜二維/三維大地電磁的有限單元法正演模擬策略. CT理論與應(yīng)用研究, 18（1）: 47-54.

王緒本, 李永年, 高永才. 1999. 大地電磁測(cè)深二維地形影響及其校正方法研究. 物探化探計(jì)算技術(shù), 21（4）: 327-332.

夏訓(xùn)銀, 王建新, 劉俊昌等. 2012. 大地電磁測(cè)深在南盤(pán)江北部地區(qū)油氣勘探中的應(yīng)用. 物探與化探, 36（2）: 174-179.

詹艷, 趙國(guó)澤, 王立鳳等. 2014. 西秦嶺與南北地震構(gòu)造帶交匯區(qū)深部電性結(jié)構(gòu)特征. 地球物理學(xué)報(bào), 57（8）: 2594-2607, doi: 10.6038/cjg20140819.

張交東, 楊曉勇, 劉成齋等. 2012. 大別山北緣深部結(jié)構(gòu)的高精度重磁電震解析. 地球物理學(xué)報(bào), 55（7）: 2292-2306.

張昆, 魏文博, 葉高峰. 2008. 二維有限元大地電磁正演模擬在Matlab上的實(shí)現(xiàn). 地震地磁觀(guān)測(cè)與研究, 29（5）: 83-88.

張翔, 胡文寶, 嚴(yán)良俊等. 1999. 大地電磁測(cè)深中的地形影響與校正. 江漢石油學(xué)院學(xué)報(bào), 21（1）: 37-41.

趙廣茂, 李桐林, 王大勇等. 2008. 基于二次場(chǎng)二維起伏地形MT有限元數(shù)值模擬. 吉林大學(xué)學(xué)報(bào)（地球科學(xué)版）, 38（6）: 1055-1059. 趙慧, 劉穎, 李予國(guó). 2014. 自適應(yīng)有限元海洋大地電磁場(chǎng)二維正演模擬. 石油地球物理勘探, 49（3）: 578-585.

（本文編輯 何燕）

A study of two-dimensional MT inversion with steep topography using the adaptive unstructured finite element method

HAN Qi1, HU Xiang-Yun1*, CHENG Zheng-Pu1, YANG Bing-Nan2, CAI Jian-Chao1, WEI Wei1

1HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China2GuizhouGeologyandMineralsBureau，GuizhouTongren554300,China

The structured grid has its limit when simulating complicated topography, especially applying the MT method in mountainous areas where the effect of terrain cannot be ignored. Unstructured triangular grids permit efficient discretization of complex modelling domains such as those containing topography, dipping strata and multiple scale structures. To ensure the accurate result of forward modeling, we refine grids according element errors during gridding the model. We combine the adaptive finite element and the Occam method, and develop a parallel processing scheme that can efficiently improve the computational speed of forward modeling and inversion. By comparing two synthetic models, first we discuss the terrain effect by building two different models. Then forward modeling and inversion are conducted on a series of models with steep topography. The application of this approach to processing of real data from Southwest China shows that the adaptive unstructured finite element method is a powerful tool to analyze MT data of complex features such as steep surface topography.

Magnetotelluric; Adaptive finite element; Topography; Inversion

國(guó)家科技專(zhuān)項(xiàng)“深部礦產(chǎn)資源立體探測(cè)技術(shù)與試驗(yàn)”（SinoProbe-03）、國(guó)家自然科學(xué)基金項(xiàng)目（41274077，41474055）、中國(guó)地質(zhì)調(diào)查局項(xiàng)目（12120113101800）聯(lián)合資助.

韓騎，女，1990年生，碩士研究生，主要從事大地電磁二維起伏地形正反演研究.E-mail:hanqi426@gmail.com

*通訊作者 胡祥云，男，1966年生，教授，主要從事電法勘探理論及應(yīng)用研究.E-mail:xyhu@cug.edu.cn

10.6038/cjg20151228.

10.6038/cjg20151228

P631

2015-05-28，2015-11-12收修定稿

韓騎， 胡祥云， 程正璞等. 2015. 自適應(yīng)非結(jié)構(gòu)有限元MT二維起伏地形正反演研究.地球物理學(xué)報(bào)，58（12）:4675-4684,

Han Q, Hu X Y, Cheng Z P, et al. 2015. A study of two-dimensional MT inversion with steep topography using the adaptive unstructured finite element method.ChineseJ.Geophys. （in Chinese）,58（12）:4675-4684,doi:10.6038/cjg20151228.