鐘藝文，劉　羽

(桂林理工大學 a.地球科學學院，b.機械與控制工程學院，桂林　541004)

?

網(wǎng)格剖分及邊界對MT Occam反演的影響研究

鐘藝文a，劉羽b*

(桂林理工大學 a.地球科學學院，b.機械與控制工程學院，桂林541004)

摘要：在大地電磁有限元模擬計算中網(wǎng)格剖分及邊界放置是否得當，影響著最終的反演結(jié)果?；贠ccam反演法從網(wǎng)格剖分及邊界方面，設計不同模型，討論不同網(wǎng)格對反演精度的影響。研究表明，對于二維介質(zhì)體，只要將網(wǎng)格邊界置于適當倍數(shù)趨膚深度的距離處，便可忽略截斷邊界的影響，無需放置在無窮遠。在同一頻段及不均勻剖分的前提下，粗細網(wǎng)格均可獲得良好的反演結(jié)果。細網(wǎng)格的反演精度高于粗網(wǎng)格，但未能很好地圈定異常體的邊界。粗網(wǎng)格剖分下，MT Occam法對低阻異常體反演效果較好，而對高阻異常體效果較差。

關鍵詞：有限單元法；網(wǎng)格剖分；邊界；MT Occam反演

0引言

有限單元法因理論完善、能模擬復雜結(jié)構(gòu)等優(yōu)點被廣泛應用于地球物理場的數(shù)值模擬中。Coggon[1]最先將有限單元法引入大地電磁領域，經(jīng)Wannamaker[2]、Rodi[3]電磁場的有限元模擬的發(fā)展與完善，有限單元法進入實用階段。在國內(nèi)陳樂壽[4]首先介紹了有限單元法在大地電磁正演計算中的應用，此后徐世浙[5]等對限單元法進行了深入研究和發(fā)展。如今有限單元法已成為二、三維大地電磁數(shù)值模擬的重要方法，但隨著計算機的高速發(fā)展，人們對反演的精度要求越來越高。因此許多學者[6-8]從算法精度方面(如插值函數(shù)、線性方程組的求解)進行改進研究，而在區(qū)域剖分方面的研究相對較少[10-13]。區(qū)域剖分是插值函數(shù)選取和解線性方程組的基礎，剖分是否合理所產(chǎn)生的誤差影響遠遠大于算法方面計算精度的高低，即使有再高精度的求解系統(tǒng)，網(wǎng)格剖分得不合理也未必能得到合理的結(jié)果。區(qū)域剖分單元通常有四邊形和三角形兩種形式，前者易于剖分具有通用性，后者可以模擬復雜地形及不規(guī)則地質(zhì)體。這里基于Occam反演法，區(qū)域剖分引用的是Wannamaker等[2]將矩形單元二次剖分成4個三角形網(wǎng)格的方案。對于網(wǎng)格剖分中邊界的選取傳統(tǒng)上都選擇放置于無窮遠處，但實際模擬區(qū)域并未能取到無窮遠，所以在網(wǎng)格的邊界處勢必存在截斷邊界的影響。對于網(wǎng)格邊界的選取，在追求獲取高精度結(jié)果的同時還得兼顧計算時間和內(nèi)存需求。作者在前人[10-14]研究基礎上，同樣用趨膚深度δ作為網(wǎng)格邊界選取的依據(jù)，依次改變網(wǎng)格的邊界，對模型的正反演計算進行比較，總結(jié)網(wǎng)格邊界置于何處方可忽略截斷邊界的影響。正演是反演的基礎，作者先研究網(wǎng)格剖分對正演的影響；在網(wǎng)格剖分疏密對正演影響不大的基礎下，繼而研究總結(jié)網(wǎng)格剖分對反演的影響，為后續(xù)解釋工作提供依據(jù)。

1Occam反演法理論

地球物理中觀測數(shù)據(jù)是有限的、離散的，同時存在觀測誤差，因此造成地球物理反演具有不適定性。在眾多的反演方法中，Occam反演是能夠克服此不足的方法之一。Constable等[15-17]提出了Occam法反演理論，Occam反演法尋求的是最光滑反演模型，因此引入了粗糙度(R)項以壓制非數(shù)據(jù)產(chǎn)生的冗余構(gòu)造。Occam反演算法在計算反演模型m的過程中，不僅要求正演響應數(shù)據(jù)與實測數(shù)據(jù)d的擬合差在設定精度內(nèi)，同時要求粗糙度盡可能小。因此，使用拉格朗日乘子μ來平衡模型光滑和數(shù)據(jù)擬合程度，最終所構(gòu)造的目標泛函如式(1)所示。

(1)

(2)

(3)

其中：Nm是模型參數(shù)個數(shù)；Npt是約束項項數(shù)。

將式(1)中的問題線性化，引用公式F[mk+△]=F[mk]+JK△,△=(mk+1-mk)表示的是模型修改量，mk是第k迭代的解，Jk為雅克比矩陣，即F[mk]相對于mk的偏導數(shù)。因此可以得到模型的迭代公式(式(4))。

mk+1(μ)=[μ(?T?+TT)+(WJk)TWJK]-1

(4)

2模型計算

2.1網(wǎng)格邊界的影響

在二維結(jié)構(gòu)中，在地面使用右旋制，令x軸為走向方向向北，y軸垂直x軸水平向東，z軸垂直向下。如圖1所示，假設u為所要求的電磁場，在不同模式下，令上邊界AB處的u都為常數(shù)“1”，值得一提的是TM模式下AB邊直接取在地面，而TE模式下，為了消除地下電磁場的變化對地面電磁場的影響，AB邊必須距離地面一定的距離。左右邊界AC、BD應取在距橫向不均勻體足夠遠的地方，使得在該處電磁場沿深度的分布可以被看成是與地下一維介質(zhì)或均勻半空間介質(zhì)時的情況相類似。對于下邊界，要求異常場在CD邊上的場值為“0”，CD邊以下可看作是均勻半空間，所以CD邊滿足第三類邊界條件。綜合以上所討論u滿足的邊界條件如式(5)所示。

(5)

圖1　研究區(qū)域及坐標Fig.1　Study area and coordinates

在做邊界截取研究時，網(wǎng)格尺寸對計算結(jié)果有很大影響，所以在不改變網(wǎng)格尺寸的情況下，只對外邊界做調(diào)整?；赪annamaker等[2]在PW2D的做法，TE模式下，將研究區(qū)域總橫向距離的一半作為AB邊的距離。MT Occam反演不依賴初始模型，通常取均勻半空間，電阻率取100 Ω·m。正演網(wǎng)格數(shù)為59(橫向)×17(縱向)，Occam要求反演網(wǎng)格是正演網(wǎng)格的子集，所以正反演網(wǎng)格必須是對齊的。反演網(wǎng)格設計：橫向按等間距1 000 m剖分并取兩個間距作為網(wǎng)格單元的寬，在此基礎上首尾兩網(wǎng)格單元的寬分別由±1 km、±3 km、±9 km、±27 km、±81 km、±243 km、±729 km七個間距組成；縱向網(wǎng)格由地面往下除前三個是按等間距剖分外，以下的都是按常數(shù)比例增長剖分，縱向最后一層的厚度則是由多個縱向剖分間距組成。頻率分別為：0.937 5 Hz，0.625 Hz，0.468 8 Hz，0.312 5 Hz，0.234 4 Hz，0.156 3 Hz，0.117 2 Hz，0.078 13 Hz，0.058 59 Hz，0.039 06 Hz，0.029 3 Hz，0.019 53 Hz，0.014 65 Hz，0.009 77 Hz，0.007 33 Hz，0.004 88 Hz，0.003 66 Hz，0.002 44 Hz，0.001 83 Hz，0.001 22 Hz。

作者采用二維異常嵌入體模型做模擬計算，模型如圖 2所示。若嵌入體是低阻體，等于1 Ω·m。頻率取0.390 6 Hz，取背景電阻率為100 Ω·m，依據(jù)公式可知，趨膚深度大約是8 km。按上述的網(wǎng)格剖分左右邊界已置于足夠遠。根據(jù)吳娟[10]、湯井田等[9]研究，只要下邊界置于1 δ距離以上，就不會影

響計算精度。因縱向剖分的層數(shù)足夠多，下邊界所在的地方也可以認為是足夠遠。二維介質(zhì)沒有解析解，所以認為在此邊界條件下計算得的是相當精確的數(shù)值解，所以把計算得到的這組數(shù)據(jù)作為參考解。

圖2　二維地質(zhì)體模型Fig.2　Two-dimensional geological model

首先在以上的網(wǎng)格剖分為基礎，將正演網(wǎng)格的左右邊界分別置于120 km、40 km、24 km處。在測點O處觀測得到的視電阻率如表1所示。

表1　兩側(cè)邊界變化下視電阻率觀測結(jié)果Tab.1　Apparent resistivity result when changing the both sides of mesh boundary

注：R120、R40、R24分別是左右邊界置于120 km、40 km、24 km處觀測的視電阻率；ε為誤差百分比。

從表1可知，TE模式下，當網(wǎng)格邊界從1 000 km外不斷縮減到120 km處時，所觀測到的視電阻率沒有變化。當縮減至40 km(5δ)、24 km時，相對高頻所觀測到的視電阻率的相對誤差在1 %左右。隨著頻率繼續(xù)減小，相對誤差增加到了5%左右。對TM模式作了相關研究，發(fā)現(xiàn)不管邊界取在120 km還是40 km或者24 km處，所觀測到的視電阻率相對誤差都維持在小于1 %的數(shù)量級。以上實驗表明，只要將兩側(cè)邊界置于5 δ以上的距離處，便可以當作不受截斷邊界影響。雖然在相對低頻處會存在一些小的偏差，但整體上可以反映淺部及深部的電性分布信息。

其次，在正演網(wǎng)格數(shù)59(橫向)×17(縱向)基礎上，反演網(wǎng)格數(shù)為24×9+13×3，兩側(cè)邊界分別置于120 km、40 km、24 km處，研究邊界的選取對反演結(jié)果的影響。

Occam反演所尋求的目標擬合差是“1”，反演規(guī)定的最大迭代次數(shù)是20次。由表2的數(shù)據(jù)可知，反演的最終擬合差精度都相當高，但是Occam反演結(jié)果尋求的是最光滑模型，所以還得考慮粗糙度這方面。隨著左右邊界的位置不斷拉近，反演模型的粗糙度會逐漸增加，反演不能正常收斂，以至于反演的迭代次數(shù)都達到限定次數(shù)時反演才結(jié)束，因此就會增加反演時間。綜合上述討論，截斷邊界對正反演結(jié)果都有一定的影響，左右邊界都應取在大于5δ距離處，但也無需置于無窮遠。

表2　反演取不同網(wǎng)格邊界的反演結(jié)果Tab.2　Inversion results when inversion take different boundary of grid

2.2粗細網(wǎng)格對反演結(jié)果的影響

反演最終結(jié)果是每個單元上的電阻率值，求解方程組后將節(jié)點場值作垂向偏導，采用差分公式計算偏導數(shù)，由此可知網(wǎng)格剖分粗細與偏導數(shù)計算存在某種關系。以下將從網(wǎng)格的粗細做研究，依然采用如圖2所示模型，均勻半空間中有一個矩形異常體，可取1 Ω·m或者1 000 Ω·m。網(wǎng)格剖分成兩種形式：粗網(wǎng)格，正演網(wǎng)格數(shù)都是59(橫向)×35(縱向)，反演網(wǎng)格數(shù)為303，測點數(shù)為21，頻點數(shù)是20；細網(wǎng)格，正演網(wǎng)格數(shù)是103(橫向)×35(縱向)，反演網(wǎng)格數(shù)為889，測點數(shù)為21，頻點數(shù)是20。

反演是一個不斷正演的過程，為了量化影響反演的因素，首先討論了不同網(wǎng)格剖分對正演的影響。由圖3可知，低阻模型的低頻段粗細網(wǎng)格正演的數(shù)據(jù)有些許的偏差，整體上粗網(wǎng)格與細網(wǎng)格都獲得基本一致的正演的結(jié)果。由此說明，只要網(wǎng)格橫縱向剖分的尺寸在某個范圍內(nèi)網(wǎng)格剖分的疏密對正演計算影響不大。

整個研究過程使用的是MT Occam 3.0版本程序，Occam 反演需要輸入四個文件，即DATA、MESH、MODEL、startup，DATA文件是經(jīng)正演而產(chǎn)生的，Occam反演的最終結(jié)果是反演網(wǎng)格數(shù)個電阻率，保存在迭代文件ITERXX中，這些電阻率是地下介質(zhì)的綜合反映，即某種平均下的電阻率，所以會有一定偏差。使用畫圖程序?qū)⒆罱K迭代文件中電阻率進行繪圖便可得到圖4、圖5中的圖形。由圖4可知，不管是粗網(wǎng)格還是細網(wǎng)格在埋深5 km處都出現(xiàn)低阻異常，異常體中心電阻率值與實驗值相差甚微，異常體中心位置都能與初始模型基本吻合，但細網(wǎng)格的模擬精度更優(yōu)于粗網(wǎng)格。而左右邊界方面，粗細網(wǎng)格異常體都未能很好的重現(xiàn)初始模型的邊界。細網(wǎng)格上邊界出現(xiàn)向上擴撒現(xiàn)象，而粗網(wǎng)格沒有。粗細網(wǎng)格的下邊界都出現(xiàn)異常體向下擴散延伸現(xiàn)象，細網(wǎng)格更是明顯。圖5是高阻模型下的反演結(jié)果，粗細網(wǎng)格反演的整個模擬區(qū)域都未出現(xiàn)明顯的高阻異?，F(xiàn)象，很難判斷異常體的中心位置。出現(xiàn)此情況主要是由于在正演時產(chǎn)生的數(shù)據(jù)偏差較大加上繪圖的色彩較容易混淆。結(jié)合反演數(shù)據(jù)分析可知，在原設計異常體位置處反演得出的電阻率略微的高于周邊的電阻率，因此從數(shù)據(jù)方面是可以圈定異常體。

圖3　不同網(wǎng)格下測點O的視電阻率曲線對比圖Fig.3　Comparison of apparent resistivity curve at station O of different grid(a)低阻模型；(b)高阻模型

圖4　低阻模型下不同網(wǎng)格的MT Occam反演Fig.4　The MT Occam inversion of different grid in low resistance model(a)粗網(wǎng)格；(b)細網(wǎng)格

圖5　高阻模型下不同網(wǎng)格的MT Occam反演Fig.5　The MT Occam inversion of different grid in high resistance model(a)粗網(wǎng)格；(b)細網(wǎng)格

3結(jié)論

針對Occam框架下的MT反演問題，通過將正反演網(wǎng)格的邊界置于不同位置，用不同反演網(wǎng)格作對比實驗得出的結(jié)果可知：

1)用趨膚深度作為網(wǎng)格左右邊界劃分的度量，正反演網(wǎng)格的左右邊界置于5倍趨膚深度以上但無需置于無窮遠，便可忽略截斷邊界的影響。

2)在同一頻段內(nèi)，細網(wǎng)格的模擬精度略高于粗網(wǎng)格，而細網(wǎng)格在圈定異常體范圍方面稍有欠缺。粗網(wǎng)格剖分下Occam反演法對低阻異常相對靈敏，反演結(jié)果精度相當高，而反演高阻異常體得的結(jié)果較差，與真實相差甚遠。

3)基于淺部剖分都相對比較細，而深部則是相對較粗的不均勻剖分粗細網(wǎng)格兩種形式都可獲得良好的反演效果。

參考文獻：

[1]COGGON J H.Electromagnetic and electrical modeling by the finite element method[J].Geophysics,1971,36(1):132-155.

[2]WANNAMAKER P.E，STODT J.A.A stable finite element solution for two-dimensional magnetotelluric modelling[J].Geophysical Journal of the Royal Astronomical Society，1987：277-296.

[3]RODI W L.A technique for improving the accuracy of finite element solution for magnctotelluric data [J].Geophysics，1976,44(2)：483-506.

[4]陳樂壽.有限元法在大地電磁場正演計算中的應用于改進[J].石油物探，1981,20(3):84-103.

CHEN L S.Application and improvement of finite-element method in forward calculation of geo-electromagnetic field[J].geophysical prospecting for petrole,1981,20(3):84-103.(In Chinese)

[5]徐世浙.地球物理中的有限單元法[M].北京:科學出版社，1994.

XU S Z.The finite element method in geophysics [M].Beijing:Science Press,1994.(In Chinese)

[6]陳小斌,胡文寶．有限元直接迭代算法在MT二維正演計算中的應用[J]．石油地球物理勘探，2000,35(4):487-496.

CHEN X B,HU W B.Application of finite-element direct iteration algorithm to mt 2-d forward computation[J].oil geophysical prospecting，2000，35(4):487-496.(In Chinese)

[7]馬為，陳小斌.大地電磁測深二維正演中輔助場的新算法[J].地震地質(zhì)，2008,30(2):525-533.

MA W,CHEN X B.A new algorithm for the calculation of auxiliary field in mt 2-d forward modeling [J].Seismology and Geology,2008,30(2):525-533.(In Chinese)

[8]柳建新，郭榮文，童孝忠，等.不完全LU分解預處理的BICGSTAB算法在大地電磁二維正演模擬中的應用[J].中南大學學報，2009,40(2):484-491.

LIU J X，GUO Z W,TONG X Z,et al.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling [J].Journal of central university，2009,40(2):484-491.(In Chinese)

[9]湯井田,薛帥.MT有限元模擬中截斷邊界的影響[J].吉林大學學報，2013,43(1):267-274.

TANG J T,XUE S.Influence of Truncated Boundary in FEM Numerical Simulation of MT [J].Journal of Jilin university,2013,43(1):267-274.(In Chinese)

[10]吳娟，席振銖，王鶴.網(wǎng)格尺寸及邊界對大地電磁有限元正演精度的影響[J].物探化探計算技術(shù)，2012,34(1):27-32.

WU J,XI Z Z，WANG H.Effect of grid size and boundary on MT forward modeling using finite element method [J].Computing Technology of Geophysical and Geochemical Exploration,2012,34(1):27-32.(In Chinese)

[11]朱崇利，董云,王延平，等.網(wǎng)格剖分對大地電磁測深反演精度的影響[J].物探與化探，2014,38(4):737-741.

ZHU C L,DONG Y，WANG Y P，et al.The effect of grid subdivision on the accuracy of MT inversion [J].Geophysical and Geochemical Exploration,2014,38(4):737-741.(In Chinese)

[12]朱崇利.網(wǎng)格剖分對反演的影響[J].地球物理學進展，2014,29(2):889-893.

ZHU C L.The influence of grid subdivision on inversion [J].Progress in Geophysics,2014,29(2):889-893.(In Chinese)

[13]歐東新.計算機精度和網(wǎng)格大小對大地電磁有限單元法正演的影響[J].桂林工學院學報，2007,27(3):329-332.

OU D X.Effect of computer precision and grid length on MT simulating using FEM [J].Journal of Guilin Institute of Technology，2007,27(3):329-332.(In Chinese)

[14]SHARMA P S,KAIKKONEN P.An automated finite element mesh generation and element coding in 2-D electromagnetic inversion[J].Geophysics，1998,34(3):93-114.

[15]CONSTABLE S C,PARKER R L,CONSTABLE C G.Occam’s inversion:a practical algorithm for generating smooth models from electromagetic sounding Data [J] Geophysics ,1987,52(3):289-300.

[16]DEGROOT-HEDLIN C,CONSTABLE S C.Occam’s inversion to generate smooth two dimensional models from magnetotelluric data[J].Geophysics ,1990,55(12):1613-1624.

[17]DEGROOT-HEDLIN C,STEVEN CONSTABLE .Occam’s inversion and the north American central plains electrical anomaly[J].Geomag.eoelectr,1993,45:985-99.

The study of effect of mesh generation and boundary on MT Occam inversion

ZHONG Yi-wena,LIU Yub*

(Guilin University of Technology a.College of Earth Sciences,b.College of Mechanical and Control Engineering,Guilin541004,China)

Abstract:The grid subdivision is reasonable or not and the boundary select appropriately of MT modeling using finite element method have influence on the inversion results.Based on Occam inversion method,using different models from grid subdivision and boundary settings,discussed their effect to the inversion results.The result show that mesh boundary are not necessary to place at infinity only in the appropriate multiple skin depth distance can ignore the influence of the truncation boundary on two-dimensional media.At the premise of same frequency band and non-uniform subdivision,good inversion results can be obtained from both fine mesh and coarse mesh.The inversion accuracy of fine grid better than that of the coarse grid,but not good at delineating the boundary of the abnormal body.The MT Occam results of coarse grid for the abnormal body of low resistance is better,however for abnormal body of high resistance,the results is worse.

Key words:finite element method;mesh generation;boundary;MT Occam inversion

收稿日期：2015-04-22改回日期：2015-06-09

基金項目:國家自然科學基金項目(41264005)

作者簡介：鐘藝文(1990-)，男，碩士，研究方向為電磁數(shù)值模擬，E-mail：zhongyw321 @163.com。*通信作者：劉羽(1961-)，男，教授，博士，主要研究方向為數(shù)據(jù)處理及并行計算，E-mail：lewis_5709 @163.com。

文章編號：1001-1749(2016)02-0145-06

中圖分類號：P 631.3

文獻標志碼：A

DOI：10.3969/j.issn.1001-1749.2016.02.01