網格剖分及邊界對MT Occam反演的影響研究

鐘藝文，劉　羽

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

?

網格剖分及邊界對MT Occam反演的影響研究

鐘藝文a，劉羽b*

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

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

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

0引言

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

1Occam反演法理論

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

(1)

(2)

(3)

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

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

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

(4)

2模型計算

2.1網格邊界的影響

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

(5)

圖1　研究區域及坐標Fig.1　Study area and coordinates

在做邊界截取研究時，網格尺寸對計算結果有很大影響，所以在不改變網格尺寸的情況下，只對外邊界做調整?；赪annamaker等[2]在PW2D的做法，TE模式下，將研究區域總橫向距離的一半作為AB邊的距離。MT Occam反演不依賴初始模型，通常取均勻半空間，電阻率取100 Ω·m。正演網格數為59(橫向)×17(縱向)，Occam要求反演網格是正演網格的子集，所以正反演網格必須是對齊的。反演網格設計：橫向按等間距1 000 m剖分并取兩個間距作為網格單元的寬，在此基礎上首尾兩網格單元的寬分別由±1 km、±3 km、±9 km、±27 km、±81 km、±243 km、±729 km七個間距組成；縱向網格由地面往下除前三個是按等間距剖分外，以下的都是按常數比例增長剖分，縱向最后一層的厚度則是由多個縱向剖分間距組成。頻率分別為：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，依據公式可知，趨膚深度大約是8 km。按上述的網格剖分左右邊界已置于足夠遠。根據吳娟[10]、湯井田等[9]研究，只要下邊界置于1 δ距離以上，就不會影

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

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

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

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

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

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

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

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

表2　反演取不同網格邊界的反演結果Tab.2　Inversion results when inversion take different boundary of grid

2.2粗細網格對反演結果的影響

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

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

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

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

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

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

3結論

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

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

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

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

參考文獻：

[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].地震地質，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]吳娟，席振銖，王鶴.網格尺寸及邊界對大地電磁有限元正演精度的影響[J].物探化探計算技術，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]朱崇利，董云,王延平，等.網格剖分對大地電磁測深反演精度的影響[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]朱崇利.網格剖分對反演的影響[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]歐東新.計算機精度和網格大小對大地電磁有限單元法正演的影響[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-)，男，碩士，研究方向為電磁數值模擬，E-mail：zhongyw321 @163.com。*通信作者：劉羽(1961-)，男，教授，博士，主要研究方向為數據處理及并行計算，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