一種適用于電磁法數(shù)值模擬的二維建模與網(wǎng)格剖分方法

何 畏， 吳文鸝， 顧觀文， 梁 萌， 陳 實(shí)， 馮 斌

(1.中國(guó)地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，廊坊 065000；2.國(guó)土資源部 地球物理電磁法探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，廊坊 065000)

一種適用于電磁法數(shù)值模擬的二維建模與網(wǎng)格剖分方法

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

(1.中國(guó)地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，廊坊 065000；2.國(guó)土資源部 地球物理電磁法探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，廊坊 065000)

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

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

0 引言

近年來(lái)電磁法處理解釋技術(shù)發(fā)展迅速，已廣泛應(yīng)用于礦產(chǎn)資源勘查與開發(fā)、環(huán)境監(jiān)測(cè)與保護(hù)、自然災(zāi)害預(yù)防與治理等領(lǐng)域[1-3]。由于電磁法三維反演技術(shù)的實(shí)用化還需進(jìn)一步地實(shí)踐[4-7]，以及一維的局限性，學(xué)者們處理解釋數(shù)據(jù)時(shí)應(yīng)用二維處理方法目前較廣泛[8-12]。如何提高數(shù)據(jù)處理效率、將先進(jìn)方法技術(shù)推廣應(yīng)用，一直是二維研究的重點(diǎn)?？梢暬换ソ：途W(wǎng)格剖分方法,是提高電法數(shù)值模擬效率和方法技術(shù)實(shí)用化的重要技術(shù)支持，可解決數(shù)據(jù)處理流程繁瑣或不清晰、模型數(shù)據(jù)配置靠手工文本編輯和無(wú)可視化功能、數(shù)據(jù)傳遞或交換僅依賴于文件格式等問(wèn)題[8-12]，因此有必要設(shè)計(jì)開發(fā)建模與網(wǎng)格剖分方法。

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

筆者針對(duì)不同電磁勘探方法模擬解釋對(duì)網(wǎng)格單元的需求，采用計(jì)算機(jī)圖形學(xué)、人機(jī)交互、拓?fù)潢P(guān)系學(xué)等技術(shù)，研究了交互建模和網(wǎng)格剖分特點(diǎn)及對(duì)模擬解釋的影響[12]，設(shè)計(jì)了基于測(cè)線信息的二維交互建模、四邊形、三角形剖分流程和方法。這些建模與網(wǎng)格剖分技術(shù)的實(shí)現(xiàn)，為電磁法正演模擬解釋方法的實(shí)際推廣應(yīng)用提供了獲取網(wǎng)格模型工具，且利用計(jì)算機(jī)可視化與交互技術(shù)使得建模與網(wǎng)格剖分具有流程清晰、可視化好、效率高等特點(diǎn)。

1 二維建模方法

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

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

多邊形搜索方法是二維交互建模的關(guān)鍵技術(shù)，基于圖形學(xué)中樹結(jié)構(gòu)(Tree Structure)原理，多邊形搜索轉(zhuǎn)化為搜索閉合曲線樹環(huán)，曲線樹的環(huán)搜索方法設(shè)計(jì)如圖2所示。

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

2 網(wǎng)格剖分設(shè)計(jì)

電磁法數(shù)據(jù)的二維模擬解釋常使用三角形或四邊形網(wǎng)格數(shù)據(jù)，因此需將基于勘探線信息構(gòu)建的矢量模型進(jìn)行相關(guān)網(wǎng)格剖分并完成屬性填充，觀測(cè)區(qū)與擴(kuò)展區(qū)設(shè)置、網(wǎng)格剖分方法流程如3所示。

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

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

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

2.1 三角形網(wǎng)格

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

圖3 觀測(cè)區(qū)、擴(kuò)展區(qū)設(shè)置與網(wǎng)格剖分流程Fig.3 Flow chart of the observation area, the extended area and the mesh generation(a)觀測(cè)區(qū)與擴(kuò)展區(qū)設(shè)置；(b)網(wǎng)格剖分流程

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

圖5 起伏地形下觀測(cè)區(qū)剖面三角網(wǎng)格剖分Fig.5 Triangular mesh generation principle in undulating topography

2.2 四邊形網(wǎng)格

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

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

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

(1)

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

3 實(shí)例與分析

為了驗(yàn)證本文建模與網(wǎng)格剖分方法的有效性，在Win7系統(tǒng)、Intel Xeon E52630 CPU @2.2GHz計(jì)算機(jī)上，利用Qt 5.1.0(32 bit)開發(fā)工具、C++語(yǔ)言，實(shí)現(xiàn)了建模與網(wǎng)格剖分方法。

3.1 交互建模、三角網(wǎng)格剖分與二維正演模擬

實(shí)驗(yàn)測(cè)線上有41個(gè)測(cè)點(diǎn)，點(diǎn)距為50 m，最小高程為0 m，最大高程為60 m，橫向模型起始坐標(biāo)點(diǎn)為0，測(cè)深深度為700 m；垂向網(wǎng)格第一層厚度為10 m，測(cè)深增長(zhǎng)率為1.1(實(shí)際可變)，擴(kuò)展區(qū)橫向擴(kuò)展第一網(wǎng)格寬度、數(shù)量、增長(zhǎng)率分別為50 m、10個(gè)、1.3，垂向擴(kuò)展網(wǎng)格數(shù)量為10(垂向擴(kuò)展增長(zhǎng)率和測(cè)深深度增長(zhǎng)率一致，垂向擴(kuò)展第一層網(wǎng)格高度與觀測(cè)區(qū)最后一層網(wǎng)格高度相同)，建模、剖分與MT二維正演模擬結(jié)果如圖7所示。

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

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

3.2 四邊形網(wǎng)格剖分與二維正演模擬

本實(shí)驗(yàn)測(cè)點(diǎn)數(shù)、點(diǎn)距和測(cè)深、建模屬性等數(shù)據(jù)與前面相同，垂向網(wǎng)格第一層厚度為50 m，測(cè)深增長(zhǎng)率為1.2，擴(kuò)展區(qū)橫向擴(kuò)展第一網(wǎng)格寬度、數(shù)量、增長(zhǎng)率分別為50 m、20個(gè)、1.2，垂向擴(kuò)展網(wǎng)格數(shù)量為20個(gè)。觀測(cè)區(qū)內(nèi)網(wǎng)格模型下邊界以隨地形方式剖分，觀測(cè)區(qū)四邊形網(wǎng)格剖分(下邊界隨地形)、CSAMT二維正演模擬結(jié)果如圖8所示。

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

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

3.3 四邊形網(wǎng)格剖分

本實(shí)驗(yàn)測(cè)點(diǎn)數(shù)、點(diǎn)距和測(cè)深、建模屬性等數(shù)據(jù)同上，垂向網(wǎng)格第一層厚度為20 m，測(cè)深增長(zhǎng)率為1.0，擴(kuò)展區(qū)橫向擴(kuò)展第一網(wǎng)格寬度、數(shù)量、增長(zhǎng)率分別為60 m、5個(gè)、1.2，垂向擴(kuò)展網(wǎng)格數(shù)量為5個(gè)。在觀測(cè)區(qū)內(nèi)網(wǎng)格模型下邊界以平地形方式剖分，觀測(cè)區(qū)四邊形網(wǎng)格剖分結(jié)果如圖9所示。

圖9 起伏地形觀測(cè)區(qū)內(nèi)下邊界平地形、四邊形剖分Fig.9 Quadrilateral subdivision when lower boundary flat terrain

圖9是按公式(1)計(jì)算每層的頂點(diǎn)坐標(biāo)，下邊界逐漸和水平方向一致，形成的觀測(cè)區(qū)內(nèi)下邊界平地形的網(wǎng)格剖分和屬性填充示意圖。該類網(wǎng)格剖分方法為后續(xù)電磁法正演模擬預(yù)留網(wǎng)格剖分技術(shù)，可為需求該類網(wǎng)格剖分的正演模擬提供借鑒。

4 結(jié)論

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

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

3)設(shè)計(jì)了解釋剖面下邊界隨地形的四邊形網(wǎng)格剖分方法，開發(fā)了其相應(yīng)軟件模塊，并開展了CSAMT正演數(shù)值模擬試驗(yàn)，結(jié)果表明建模和網(wǎng)格剖分方法正確，滿足了CSAMT二維數(shù)值模擬需求；設(shè)計(jì)了下邊界平地形的四邊形網(wǎng)格剖分方法，為電磁法二維正演模擬提供了相關(guān)網(wǎng)格剖分借鑒。

[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]．北京：科學(xué)出版社，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].地球物理學(xué)報(bào)，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] 林庚浩, 馬天寶, 寧建國(guó). 三維前處理有限差分網(wǎng)格剖分技術(shù)的研究[J].計(jì)算力學(xué)學(xué)報(bào)，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] 曹丹平，周建科，印興耀. 三角網(wǎng)格有限元法波動(dòng)模擬的數(shù)值頻散及穩(wěn)定性研究[J].地球物理學(xué)報(bào)，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] 張志勇，劉慶成. 基于收縮二叉樹結(jié)構(gòu)網(wǎng)格剖分的大地電磁二維有限單元法正演研究[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] 劉云，王緒本. 大地電磁二維自適應(yīng)地形有限元正演模擬[J].地震地質(zhì)，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] 朱崇利. 網(wǎng)格剖分對(duì)反演影響[J].地球物理學(xué)進(jìn)展，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] 美國(guó)zonge公司GDP-32配套軟件信息.http://zonge.com/services/explorationists/exploration-geophysical-methods/

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

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

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

[17] 王林飛，雄盛青，何輝，等. 非地震地球物理軟件發(fā)展現(xiàn)狀與趨勢(shì)[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] 吳文鸝. 電法勘探工作站軟件系統(tǒng)簡(jiǎn)介[J]. 地質(zhì)與勘探，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] 顧觀文，武曄. 電法軟件設(shè)計(jì)中交互與可視化技術(shù)應(yīng)用[J]. 地質(zhì)與勘探，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] 楊元生. 無(wú)向圖與有向圖的全部生成樹的計(jì)算機(jī)算法[J]. 計(jì)算機(jī)學(xué)報(bào)，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] 王耘，胡樹根，孫偉寧，等. 基于節(jié)點(diǎn)預(yù)處理的環(huán)搜索方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，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] 何畏，吳文鸝，陳實(shí)，等. 基于無(wú)向圖的二維地質(zhì)建模設(shè)計(jì)與應(yīng)用研究[J]. 物探化探計(jì)算技術(shù)，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

國(guó)家重大科學(xué)儀器設(shè)備開發(fā)專項(xiàng)(2011YQ050060)；國(guó)家“863”高科技研究發(fā)展計(jì)劃重點(diǎn)課題(SS2014AA063110)；中國(guó)地質(zhì)大調(diào)查項(xiàng)目(12120113100600,DD20179374)

何畏(1981-)，男，碩士，高級(jí)工程師，主要研究方向?yàn)橛?jì)算機(jī)建模在地學(xué)中的應(yīng)用與物探數(shù)據(jù)處理方法技術(shù)研究，E-mail：hewei@igge.cn。

1001-1749(2017)06-0719-08

TP 317.4

A

10.3969/j.issn.1001-1749.2017.06.02