一種適用于電磁法數值模擬的二維建模與網格剖分方法

何 畏， 吳文鸝， 顧觀文， 梁 萌， 陳 實， 馮 斌

(1.中國地質科學院 地球物理地球化學勘查研究所，廊坊 065000；2.國土資源部 地球物理電磁法探測技術重點實驗室，廊坊 065000)

一種適用于電磁法數值模擬的二維建模與網格剖分方法

何 畏1,2， 吳文鸝1,2， 顧觀文1,2， 梁 萌1,2， 陳 實1,2， 馮 斌1,2

(1.中國地質科學院 地球物理地球化學勘查研究所，廊坊 065000；2.國土資源部 地球物理電磁法探測技術重點實驗室，廊坊 065000)

為滿足電磁法二維數值模擬解釋對交互建模與可視化網格剖分的需求，利用計算機圖形學、人機交互、拓撲關系學等技術，設計了基于測線、測點等信息的二維矢量建模、四邊形和三角形網格剖分方法與流程，并開發形成了軟件模塊。該方法與軟件可提供起伏地形下大地電磁、可控源音頻大地電磁法二維數值模擬的三角形、四邊形網格(下邊界隨地形變化或平地形)剖分實用化工具。將該方法應用到多種模型試驗中，取得了良好的建模與網格剖分效果。

大地電磁； 可控源音頻大地電磁； 二維建模； 網格剖分； 四邊形； 三角形

0 引言

近年來電磁法處理解釋技術發展迅速，已廣泛應用于礦產資源勘查與開發、環境監測與保護、自然災害預防與治理等領域[1-3]。由于電磁法三維反演技術的實用化還需進一步地實踐[4-7]，以及一維的局限性，學者們處理解釋數據時應用二維處理方法目前較廣泛[8-12]。如何提高數據處理效率、將先進方法技術推廣應用，一直是二維研究的重點。可視化交互建模和網格剖分方法,是提高電法數值模擬效率和方法技術實用化的重要技術支持，可解決數據處理流程繁瑣或不清晰、模型數據配置靠手工文本編輯和無可視化功能、數據傳遞或交換僅依賴于文件格式等問題[8-12]，因此有必要設計開發建模與網格剖分方法。

電磁法數據處理與解釋軟件現狀發展并不同步，國外軟件研發、實際資料處理應用較早。在電法儀器配套使用的軟件開發方面，①美國Zonge公司開發了與多功能電法儀器GDP-32配套的軟件[13]，主要針對各電法子方法進行數據預處理、一維/二維反演，如帶地形二維反演TDIP模塊、CSAMT(遠區數據)/AMT二維反演模塊，美國AGI公司開發了與其高密度IP電法儀器的軟件模塊[14]，用于IP的二維反演成像EarthImager 2D模塊，該類與儀器配套軟件通用性不強，不具有二維交互建模與網格剖分功能；②通用類軟件以EMIGMA、WinGLink為代表：加拿大PetRos EiKon公司開發的地球物理解釋平臺EMIGMA[15]，具有多種電法的一維三維正反演功能，不具備二維正反演功能；③Geosystem與WesternGeco公司開發地球物理解釋軟件WinGLink[16]，具有MT/AMT二維反演功能，不具有二維交互建模與可視化網格剖分功能。國內電法二維處理解釋軟件方面[17]，以MTSoft(成都理工大學)、電磁勘探資料處理與解釋系統(中國石油大學)、MT-Pioneer(中國地震局地質研究所)、電法數據處理系統GeoElectro(吉林大學)以及電法工作站WEM(中國地質科學院物化探所)等為代表，其中MTSoft、電磁勘探資料處理與解釋系統、MT-Pioneer是針對大地電磁的數據預處理與反演開發的軟件系統，GeoElectro具有TDIP二維反演功能，WEM具有TDIP、MT、CSAMT、TEM、SIP五種電法方法的二維正反演[18-19]，這些軟件的二維正演模塊含有二維簡單體建模與網格剖分功能。

筆者針對不同電磁勘探方法模擬解釋對網格單元的需求，采用計算機圖形學、人機交互、拓撲關系學等技術，研究了交互建模和網格剖分特點及對模擬解釋的影響[12]，設計了基于測線信息的二維交互建模、四邊形、三角形剖分流程和方法。這些建模與網格剖分技術的實現，為電磁法正演模擬解釋方法的實際推廣應用提供了獲取網格模型工具，且利用計算機可視化與交互技術使得建模與網格剖分具有流程清晰、可視化好、效率高等特點。

1 二維建模方法

二維建模是網格剖分的基礎，其正確、合理的建模方式對網格剖分尤為重要。本建模方法的設計主要依賴于計算機圖形學、人機交互和可視化技術，包括兩個步驟：①建模數據準備，其數據來源于電磁勘探觀測區內某測線信息(含測線長度、所有測點坐標、測深深度等)，形成帶地形的二維模型邊界區域，可將已知地質或地震剖面作為參考背景資料；②人機交互建模，人機交互式鼠標勾畫曲線多邊形，基于無向圖搜索環多邊形方法[20-23]，獲取多邊形模型，設置或修改多邊形屬性，形成二維矢量模型。二維建模流程如圖1所示。

圖1 二維建模與網格剖分流程Fig.1 Flow chart of the two-dimensional modeling and mesh generation

多邊形搜索方法是二維交互建模的關鍵技術，基于圖形學中樹結構(Tree Structure)原理，多邊形搜索轉化為搜索閉合曲線樹環，曲線樹的環搜索方法設計如圖2所示。

圖2是將獲取多邊形轉化為搜索樹形環的方法，其父子樹根節點表示環的起始曲線，通過曲線首尾節點的索引關系搜索孩子節點，逐級搜索，最終形成閉合樹形環。該搜索方法涉及名詞或術語：①曲線(Curves Tree)，平面內離散節點連成非閉合、無相交的線段集；②合理環(Reasonable Ring)，由曲線集形成閉合環，但不一定是目標環；③正確環(Correct Ring)，所有合理環中，除了最外圍邊界線段之外，存在一個環內不包含其它線段的環；④最小面積法(Minimum Area Method)，正確環的特點是所有合理環中面積最小的環，即最小面積法的用途是判定環是否為正確環。

2 網格剖分設計

電磁法數據的二維模擬解釋常使用三角形或四邊形網格數據，因此需將基于勘探線信息構建的矢量模型進行相關網格剖分并完成屬性填充，觀測區與擴展區設置、網格剖分方法流程如3所示。

根據不同電磁勘探方法，在二維矢量模型區域進行觀測區三角形、四邊形網格剖分，以及擴展區網格剖分，這些網格單元屬性由網格中心點與二維矢量多邊形模型的包含關系決定，其填充原理是網格中心點坐標在多邊形內時，將該多邊形的屬性(電阻率或極化率物性、顏色等屬性)賦值給該網格，否則將其屬性賦值為空氣層網格單元標志(如物性設為負數，可視化顏色為白色或不可見狀態)。

圖2 環樹形搜索方法流程圖Fig.2 Flow chart of ring tree search method

圖4(a)是四邊形ABCD(線段AB、CD方向為鉛垂方向，線段AD、BC任意方向)，該形狀的四邊形可滿足激發極化法(時間域激發極化法(TDIP)、頻率域激發極化法(SIP))正演初始模型特點。四邊形剖分參數確定四邊形AB、CD邊為網格垂向高度，測點間點距和網格橫向剖分方式(測區下邊界隨地形或平地形)確定AD、BC邊。如圖4(b)所示，將該四邊形兩對角線相交，交點為點O，形成三角形Δ1-ABO、Δ2-BCO、Δ3-CDO、Δ4-DOA(三角形號和節點均按逆時針排列)，以該種方式形成的三角形單元可用于大地電磁(MT)二維正演初始三角網格。

2.1 三角形網格

圖5為起伏地形下觀測區二維剖面三角網格剖分原理示意圖，在平面坐標系XOY內，9個實心黑點為理論觀測點，測線上地形起伏綠色邊界由9個觀測點連接形成封閉曲線表達；該測線上9個測點有7個不同高程值，在這些不同高程處做與OX的平行張，并與地形邊界線求交，其交點如圖5中空心點所示；從左到右，由觀測點、交點形成不等間距的四邊形網格，這些四邊形單元按圖4(b)四邊形三角化法形成4個三角形，這些三角形單元經建模屬性填充后即可作為MT二維正演的初始模型。向左、向下、向右擴展區的三角形屬性填充同觀測區內方法一致，其相應多邊形區域屬性由相關邊界的觀測區模型屬性確定，地形邊界線以上空氣層擴展區的三角形屬性設為空氣層標志(物性值為負數)。

圖3 觀測區、擴展區設置與網格剖分流程Fig.3 Flow chart of the observation area, the extended area and the mesh generation(a)觀測區與擴展區設置；(b)網格剖分流程

圖4 三角形、四邊形單元Fig.4 Triangular and quadrilateral element(a)四邊形；(b)四邊形三角化

圖5 起伏地形下觀測區剖面三角網格剖分Fig.5 Triangular mesh generation principle in undulating topography

2.2 四邊形網格

四邊形網格剖分涉及CSAMT法中模型下邊界隨地形(和地形邊界變化一致)、平地形(為電磁法二維正演模擬需求預留該功能)兩種情形，其網格橫向間距為測點距離(圖6)。

圖6(a)為下邊界隨地形四邊形網格剖分原理圖，實現下邊界隨地形的方法是,每一層垂向網格單元的橫向邊和由地形邊界形成第一層四邊形單元上頂邊的變化一致(橫向坐標、斜率、距離)，四邊形橫向邊長為兩測點間距，垂向長度可根據用戶剖分需要設置等長度或按一定比例增長方式。圖6(b)為建模區域下邊界為平地形網格剖分原理圖，除第一層外的每一層垂向上四邊形網格的右下角節點坐標(左下角類似)按公式(1)計算。

圖6 起伏地形下觀測區剖面四邊形網格剖分原理Fig.6 Quadrilateral mesh generation principle in undulating topography(a)下邊界隨地形；(b)下邊界平地形

(1)

式中：rbp.z為右下角節點網格垂向坐標值；t_v[j]、t_v[j+1]為第j、j+1個高程值；z_h[i]為第i層網格垂向高度；H為測深深度；min_t_v為最小高程。整個剖面的四邊形網格坐標計算方式可以按式(1)方法循環計算，計算原則是凹凸地形隨著網格垂向增加逐漸變平緩，第一層四邊形的左下、右上坐標為第二層四邊形的左上、右上坐標，依次計算，即可獲得下邊界平地形的四邊形網格。

3 實例與分析

為了驗證本文建模與網格剖分方法的有效性，在Win7系統、Intel Xeon E52630 CPU @2.2GHz計算機上，利用Qt 5.1.0(32 bit)開發工具、C++語言，實現了建模與網格剖分方法。

3.1 交互建模、三角網格剖分與二維正演模擬

實驗測線上有41個測點，點距為50 m，最小高程為0 m，最大高程為60 m，橫向模型起始坐標點為0，測深深度為700 m；垂向網格第一層厚度為10 m，測深增長率為1.1(實際可變)，擴展區橫向擴展第一網格寬度、數量、增長率分別為50 m、10個、1.3，垂向擴展網格數量為10(垂向擴展增長率和測深深度增長率一致，垂向擴展第一層網格高度與觀測區最后一層網格高度相同)，建模、剖分與MT二維正演模擬結果如圖7所示。

圖7 二維交互建模、三角剖分與正演模擬Fig.7 Two-dimensional interactive modeling and triangulation and forward simulation(a)鼠標交互勾畫地電模型;(b)二維模型構建;(c)測區三角剖分;(d)測區、擴展區三角剖分;(e)TE模式二維正演視電阻率擬斷面;(f)TE模式二維正演相位擬斷面;(g)TM模式二維正演視電阻率擬斷面;(h)TM模式二維正演相位擬斷面

圖7中灰色線條代表三角形邊，其中圖7(a)是基于測線相關信息和剖分參數，利用筆者設計的測點不同高程之間插值地形邊界節點形成的初始網格模型和人機交互式勾畫模型效果的疊加圖；圖7(b)是本文設計的曲線樹法搜索環形成的剖面模型，該剖面包含藍色與橘色覆蓋層、黃色基巖和紅色塊區域異常，電阻率分別為500 Ω·m、1 000 Ω·m、200 Ω·m、10 Ω·m；圖7(c)是觀測區內基于四邊形交叉對角線形成的三角網格剖分和屬性填充的效果圖；圖7(d)是觀測區、擴展區內基于四邊形交叉對角線形成的三角網格剖分和屬性填充的效果圖；正演結果符合理論模型的響應規律。從圖7中可以看出，建模與三角網格剖分方法原理設計合理、技術流程清晰、可視化效果明顯，實現了交互建模與可視化網格剖分的目標，為大地電磁二維數值模擬方法的實用化提供了基礎性工具。

3.2 四邊形網格剖分與二維正演模擬

本實驗測點數、點距和測深、建模屬性等數據與前面相同，垂向網格第一層厚度為50 m，測深增長率為1.2，擴展區橫向擴展第一網格寬度、數量、增長率分別為50 m、20個、1.2，垂向擴展網格數量為20個。觀測區內網格模型下邊界以隨地形方式剖分，觀測區四邊形網格剖分(下邊界隨地形)、CSAMT二維正演模擬結果如圖8所示。

圖8(a)是每層網格下邊界起伏和地形起伏變化一致；圖8(b)、圖8(c)分別是圖8(a)的YX模式二維正演視電阻率擬斷面、YX模式二維正演相位擬斷面，正演結果符合理論模型的響應規律。從圖8可以看出，下邊界隨地形網格剖分方式都在二維矢量模型上實現了正確的剖分和屬性填充，可按CSAMT二維正演數據格式要求輸出網格單元屬性進行二維正演數值模擬計算。

圖8 起伏地形測區內四邊形剖分與正演模擬Fig.8 Quadrilateral subdivision in undulating topography and forward simulation(a)下邊界隨地形、四邊形剖分;(b)YX模式二維正演視電阻率擬斷面;(c)YX模式二維正演相位擬斷面

3.3 四邊形網格剖分

本實驗測點數、點距和測深、建模屬性等數據同上，垂向網格第一層厚度為20 m，測深增長率為1.0，擴展區橫向擴展第一網格寬度、數量、增長率分別為60 m、5個、1.2，垂向擴展網格數量為5個。在觀測區內網格模型下邊界以平地形方式剖分，觀測區四邊形網格剖分結果如圖9所示。

圖9 起伏地形觀測區內下邊界平地形、四邊形剖分Fig.9 Quadrilateral subdivision when lower boundary flat terrain

圖9是按公式(1)計算每層的頂點坐標，下邊界逐漸和水平方向一致，形成的觀測區內下邊界平地形的網格剖分和屬性填充示意圖。該類網格剖分方法為后續電磁法正演模擬預留網格剖分技術，可為需求該類網格剖分的正演模擬提供借鑒。

4 結論

1)針對目前電磁法勘探在多種領域的應用需求，我們結合勘探線剖面特點設計了二維交互建模與網格剖分方法與流程，為電磁法二維剖面交互模擬解釋提供了流程清晰、可視化好、效率高等特點的建模和網格剖分手段，具有較好的實際應用價值，同時為電磁法三維數值模擬解釋的建模與網格剖分提供了方法技術參考。

2)設計了基于可視化建模的三角形網格剖分方法，在起伏地形進行邊界求交、插值，加密了網格剖分密度，為更好地提高模擬效果奠定了基礎，并開發了二維建模與三角剖分軟件模塊，滿足了大地電磁二維數值模擬需求。

3)設計了解釋剖面下邊界隨地形的四邊形網格剖分方法，開發了其相應軟件模塊，并開展了CSAMT正演數值模擬試驗，結果表明建模和網格剖分方法正確，滿足了CSAMT二維數值模擬需求；設計了下邊界平地形的四邊形網格剖分方法，為電磁法二維正演模擬提供了相關網格剖分借鑒。

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

[2] WANNAMAKER P E,STODT J A，RIJO L. Two-dimensional topgraphic responses in magnetotellurics modeled using finite elements[J]. Geophysics，1986,51(11)：2131-2144.

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

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

[4] WANG WEI,WU XIAOPING,KLAUS SPITZER.Three-dimensional DC anisotropic resistivity modeling using finite elements on unstructured grids[J].Geophysical Journal International,2013,193:734-746.

[5] MENG CAO,TAN HANDONG,WANG KUNPENG.3D LBFGS inversion of controlled source extremely low frequency electromagnetic data[J].Applied Geophysics,2016,13(4):689-700.

[6] 李鑫，白登海，閆永利，等. 考慮電流型畸變的大地電磁三維反演[J].地球物理學報，2016,59(6)：2302-2315.

LI X,BAI D H,YAN Y L, et al. Three-dimensional inversion of magnetotelluric resistivity model with galvanic distortion[J]. Chinese Journal of Geophysics，2016, 59(6)：2302-2315. (In Chinese)

[7] 朱成，李桐林，楊海斌，等. 帶地形頻率域可控源電磁法三維反演研究[J].石油地球物理勘探，2016,51(5)：1031-1039.

ZHU C, LI T L, YANG H B, et al. 3D controlled source electromagnetic inversion with topogragpy in the frequency domain[J]. Oil Geophysical Prospectiong，2016, 51(5)：1031-1039. (In Chinese)

[8] 林庚浩, 馬天寶, 寧建國. 三維前處理有限差分網格剖分技術的研究[J].計算力學學報，2011，28(S)：199-203.

LIN G H，MA T B，NING J G. Study on technique of finite-difference grid subdivision in three dimensional pro-processing system[J]. Chinese Journal of Computational Mechanics，2011,28(S)：199-203.(In Chinese)

[9] 曹丹平，周建科，印興耀. 三角網格有限元法波動模擬的數值頻散及穩定性研究[J].地球物理學報，2015，58(5)：1717-1730.

CAO D P，ZHOU J K，YIN X Y. The study for numerical dispersion and stability of wave motion with triangle-based finite element algorithm[J]. Chinese Journal of Geophysics，2015,58(5)：1717-1730. (In Chinese)

[10] 張志勇，劉慶成. 基于收縮二叉樹結構網格剖分的大地電磁二維有限單元法正演研究[J].石油地球物勘探，2013，48(3)：482-487.

ZHANG Z Y，LIU Q C. 2D MT numerical simulation using FEM based on bi-tree grid[J]. Oil Geophysical Prospectiong，2013,48(3)：482-487.(In Chinese)

[11] 劉云，王緒本. 大地電磁二維自適應地形有限元正演模擬[J].地震地質，2010，32(3)：382-391.

LIU Y,WANG X B. Fem using adaptive topography in 2D MT forward modeling[J].Seismology and Geology，2010，32(3)：382-391.(In Chinese)

[12] 朱崇利. 網格剖分對反演影響[J].地球物理學進展，2014，29(2)：889-894.

ZHU C L. The influence of the grid subdivision on the inversion[J].Progress in Geophysics，2014，29(2)：889-894.(In Chinese)

[13] 美國zonge公司GDP-32配套軟件信息.http://zonge.com/services/explorationists/exploration-geophysical-methods/

[14] 美國AGI公司高密度IP電法儀器配套軟件模塊信息.https://www.agiusa.com/products

[15] 加拿大PetRos EiKon公司地球物理解釋平臺EMIGMA 信息.http://www.petroseikon.com/resources/index.php

[16] Geosystem與WesternGeco公司開發的地球物理解釋軟件WinGLink信息.http://softadvice.informer.com/Geosystem_Winglink.html

[17] 王林飛，雄盛青，何輝，等. 非地震地球物理軟件發展現狀與趨勢[J].物探與化探，2011，35(6)：837-844.

WANG L F，QIONG S Q，HE H，et al. Current status and future trends of non-seismic geophysical software[J]. Geophysical and Geochemical Exploration，2011,35(6)：837 -844.(In Chinese)

[18] 吳文鸝. 電法勘探工作站軟件系統簡介[J]. 地質與勘探，2003,39(z1):147-151.

WU W L. Introduction to electrical prospecting workstation software system[J]. Geology and Prospecting，2003，39(z1):147 -151.(In Chinese)

[19] 顧觀文，武曄. 電法軟件設計中交互與可視化技術應用[J]. 地質與勘探，2004,40(z1):243-246.

GU G W, WU Y. Applications of interactive visualization technology to developing software electrical methods[J]. Geology and Prospecting，2004，40(z1):243-246.(In Chinese)

[20] 楊元生. 無向圖與有向圖的全部生成樹的計算機算法[J]. 計算機學報，1983(2)：152-154.

YANG Y S . A computerized algorithm for finding all spanning trees of directed and undirected graphs [J]. Chinese Journal of Computers，1983(2):152-154. (In Chinese)

[21] 王耘，胡樹根，孫偉寧，等. 基于節點預處理的環搜索方法[J]. 計算機輔助設計與圖形學學報，2001，13(9):774-778.

WANG Y，HU S G，SUN W N，et al. Loop detection in 2D drawing based on node pretreatment [J]. Journal of Computer-Aided Design & Computer Graphics，2001，13(9):774-778. (In Chinese)

[22] GABOW H.N.,MYERS E. W.Filding all spanning tree of directed and undirected graphs[J]. Siam Journal on Computing，1978(7):3.

[23] 何畏，吳文鸝，陳實，等. 基于無向圖的二維地質建模設計與應用研究[J]. 物探化探計算技術，2015,37(1):123-129.

HE W，WU W L，CHEN S，et al. Design and application of two dimensional geology modeling based on undirected graphs[J].Computing Techniques for Geophysical and Geochemical Exploration，2015，37(1)：123 -129.(In Chinese)

TheMethodof2DModelingandGridSubdivisionforElectromagneticNumericalSimulation

HE Wei1,2， WU Wenli1,2， GU Guanwen1,2， LIANG Meng1,2， CHEN Shi1,2， FENG Bin1,2

(1.Institute of Geophysical and Geochemical Exploration, Chinese Academy of Geological Sciences，Langfang 065000，China；2.Laboratory of Geophysical EM Probing Technologies，MLR，Langfang 065000， China)

Interactive modeling and visualization mesh generation are very important parts in the electromagnetic two-dimensional numerical simulation. Using some techniques such as computer graphics, human-computer interaction, topological relation, etc., we designed method and process about the two-dimensional vector modeling, quadrilateral mesh generation, triangular mesh generation. In method, it contains measuring lines, measuring points information, etc.. Through the study, the corresponding software module was developed. In the he method and software, some triangulation mesh generation and quadrilateral mesh generation utility tools can be provided to solve some related problems in MT and CSAMT. A large number of model tests have been carried out, also good modeling and mesh generation effects were achieved.

MT; CSAMT; two-dimensional modeling; mesh generation; quadrilateral; triangle

2017-05-17 改回日期： 2017-08-02

國家重大科學儀器設備開發專項(2011YQ050060)；國家“863”高科技研究發展計劃重點課題(SS2014AA063110)；中國地質大調查項目(12120113100600,DD20179374)

何畏(1981-)，男，碩士，高級工程師，主要研究方向為計算機建模在地學中的應用與物探數據處理方法技術研究，E-mail：hewei@igge.cn。

1001-1749(2017)06-0719-08

TP 317.4

A

10.3969/j.issn.1001-1749.2017.06.02