激電二維數(shù)值模擬研究

陳永凌，蔣首進，謝 丹

（1．武警黃金十二支隊，成都 610059；2．武警黃金十支隊，昆明 65000）

激電二維數(shù)值模擬研究

陳永凌1，蔣首進2，謝 丹1

（1．武警黃金十二支隊，成都 610059；2．武警黃金十支隊，昆明 65000）

從點源二維地電問題出發(fā)，采用有限單元法進行了地電場進行數(shù)值模擬，采用自適應(yīng)三角剖分來實現(xiàn)起伏地表的模擬 ，針對雙邊三極裝置，實現(xiàn)了多種模型的正演研究；通過多種模型的正反演，總結(jié)異常產(chǎn)生的規(guī)律，為激電法的分析提供了有效的信息。

激發(fā)極化法；有限元法；雙邊三極

0 引言

有限元方程是由Coggon在20世紀70年代由電磁場能量最小原理導(dǎo)出的，通過計算點源二維地點問題的變分方程，得到其相應(yīng)的解。20世紀80年代，頗瑞德摩A．B研究了三維地電有限元數(shù)值模擬。二十一世紀初，Matesihat計算了二維模型使用為函數(shù)場源的三維場源的響應(yīng)。我國從1975年開始在電法勘探中引入有限元方法，取得了很好的效果。上世紀80年代鐘本善等［7］發(fā)表了基于點源二維地電的有限元正演數(shù)值模擬以及數(shù)值模擬過程中的若干問題，進行了有限單元正演模擬，推導(dǎo)了地電場模擬的有關(guān)方程；徐世浙［6］通過網(wǎng)格內(nèi)電導(dǎo)率連續(xù)變化，采用矩形網(wǎng)格剖分對地電場進行了數(shù)值模擬；阮百堯［8］提出電導(dǎo)率和電位雙線性插值，在此過程中采用線性變化的電導(dǎo)率分塊，有效地提高了有限元正演數(shù)值模擬的精度；阮百堯［8］針對結(jié)構(gòu)單元電導(dǎo)率等參數(shù)與實際不相符合的問題，相繼提出雙線性插值的模擬方法和連續(xù)線性變化數(shù)值模擬，并且優(yōu)化了系數(shù)矩陣，十分有效地提高了計算精度和速度。

作者在前人的基礎(chǔ)上，運用有限單元法進行了數(shù)值模擬，采用自適應(yīng)三角剖分來實現(xiàn)起伏地表的模擬，通過穩(wěn)定電流場的基本方程的導(dǎo)出、點源二維地電問題的討論、邊值問題的討論等。在這過程中采用雙二次插值基函數(shù)和Galerkin法將微分方程離散化，利用代數(shù)多重網(wǎng)格法求解方程組，波數(shù)域—空間域的傅立葉反變換，采用分段波數(shù)一維連分法直接計算反變換中的定積分［5］。

1 基本原理

1．1 點源二維地電問題

對于二維地電問題，我們通常假設(shè)平行于地質(zhì)體走向為y軸，其傾向方為x軸，向下方向為z軸，由于地質(zhì)體走向方向物性條件不變，則地電參數(shù)σ只隨x、z的方向發(fā)生變化，在y方向上不發(fā)生變化，即

這時候電位的基本微分式方程變?yōu)?/p>

由于在二維情況下地電參數(shù)σ只隨x、z的方向發(fā)生變化，在y方向上不發(fā)生變化，為簡化計算條件，我們使用傅氏變換

將坐標系空間中的電位U變換為二維波數(shù)域中的電位φ，式中k為空間波數(shù)。

對式（3）進行傅氏變換可以得到

式（4）就是在點源二維地電斷面時變換后的電位φ必須滿足的基本微分方程式。通過上述變換將地電三維問題轉(zhuǎn)換為二維來描述，其中k值只是一個參數(shù)，實際工作中k值為一常數(shù)，由不同的k值來求得相應(yīng)的φ值，再通過傅氏反變換來求取未知電位U。

1．2 穩(wěn)定電流場的邊值問題

流場的分布是在整個地下介質(zhì)空間分布的，而我們需要計算的只是其中的一部分空間，所以需要在區(qū)域邊界上給傅氏電位φ進行邊界條件的計算。

邊界條件同常分為三類：

其中：第一類邊界條件又叫“狄里希萊條件”；第二類邊界條件又叫“諾依曼條件”；第三類邊界條件又叫“混合邊界條件”。

上面第三類邊界條件為修正后的第三類邊界條件，式中λ為邊界修正參數(shù)，η由式（5）確定。

式中：lλ為λ點的供電電流強度；rλ為測點到場源距離為λ點的點矢徑為邊界處法線向量。

由于狄里希萊條件會產(chǎn)生計算值比解析解的值?。恢Z依曼條件會產(chǎn)生計算值比解析解的值大?；旌线吔鐥l件具有能在遠處邊界面上保持電位的物理特性，且又不需要首先段設(shè)一個“正?！钡碾妶龇植嫉膬?yōu)點，這樣可以減少邊界附近放延展網(wǎng)格的數(shù)量，從而減少了計算的工作量。同時因為在邊界上保持電位的連續(xù)，沒有電性的突變，對地下邊界的反射作用也消除。由于混合邊界條件具有較高的模擬精度，較好的反演效果，所以一般情況下在計算過程中采用混合邊界條件。

1．3 地電二維有限元法基本原理

從上面的推導(dǎo)中得到了點源二維目標函數(shù)

的極小值相對應(yīng)。

其中：φ為目標函數(shù)；D為研究區(qū)域；Γ為研究區(qū)域D的邊界。

1．4 網(wǎng)格剖分示意圖

在使用有限元法求解研究區(qū)域的時候，通常需要將研究區(qū)域剖分成許多連續(xù)的小單元。網(wǎng)格剖分種類繁多，通常有以下幾種剖分類型：①非結(jié)構(gòu)化三角網(wǎng)格剖分；②矩形三角剖分；③四邊形剖分；④三角形剖分。四邊形剖分相對其他剖分較簡單，且收斂速度較快，但四邊形剖分對地形起伏及復(fù)雜模型模擬較為粗超，不能較好地模擬其真實的形態(tài)，不適應(yīng)于起伏地形條件及復(fù)雜條件的地電模型，且在不同方向上電性結(jié)構(gòu)變化也不相同；非結(jié)構(gòu)化三角網(wǎng)格剖分具有三角單元的特征，對各種約束條件滿足，但由于網(wǎng)格數(shù)量非常大，導(dǎo)致計算速度較慢；三角網(wǎng)格剖分具有滿足地形起伏條件，可以很好地模擬異常體形態(tài)，滿足地電模型的約束條件和三角單元的特性，且網(wǎng)格節(jié)點數(shù)量適合，計算速度較快，因而成為一種常用的網(wǎng)格剖分方法。

結(jié)合以上結(jié)論，采用了如圖1所示的方式進行網(wǎng)格剖分，解決了網(wǎng)格大小與精度以及計算量之間存在的一系列矛盾，也較好地適應(yīng)了地形的起伏。

網(wǎng)格剖分的原則一般來講有以下幾點：

1）各個單元的節(jié)點只能與相鄰的單元節(jié)點相重合，但是不能成為相鄰的單元的邊上或者邊內(nèi)的點。

2）每個單元內(nèi)介質(zhì)的電導(dǎo)率為常數(shù)，在不同電導(dǎo)率介質(zhì)的分界面或者求解區(qū)的邊界上，使用三角形的一個邊互相連接，以模擬邊界線的近似折線取代邊界線。

3）三角形的三條邊或者三個頂點的差別不能太大。

4）在某些區(qū)域物性比變化劇烈，點位參數(shù)的二階或高階導(dǎo)數(shù)數(shù)值較大的區(qū)域，單元剖分應(yīng)剖分細密一些，而在電場變化平緩的地方，網(wǎng)格可以分稀疏些［5－6］。

圖1 網(wǎng)格剖分示意圖Fig．1 The grid subdivision schemes

2 模型計算

本節(jié)主要對固定點源測深采用了有限元進行了正演模擬，就多個模型實例進行數(shù)值模擬的結(jié)果分析總結(jié)。由于考慮到數(shù)據(jù)量的問題，因此在本章模型計算后成正演曲線時，取了其中的測深供電點中的7個點成圖，其位置分別為：A 1＝－50 m、A 2＝－30 m、A 3＝－10 m、A 4＝0 m、A 5＝10 m、A 6＝30 m、A 7＝50 m。

2．1 二層模型

模型為分為兩層（圖2），第一層的電阻率為50 Ω·m，極化率為η＝5．0%，厚度為15 m；第二層的電阻率為100Ω·m，極化率為η＝1．0%。

圖2 二層層狀模型示意圖Fig．2 Diagram of two layered model

從雙層模型正演測深點的曲線（圖3及圖4），我們可以看出數(shù)值模擬較好地體現(xiàn)了地下有的異常體。

2．2 單正方柱體模型

模型為10 m×10 m低阻、高極化正方柱體（圖5），中心位置為（0 m，－10 m），電阻率為50Ω·m，極化率為η＝5．0%；背景電阻率為100Ω·m，極化率為η＝1．0%。

圖3 二層層狀模型固定點源測深正演模擬各點源測深視電阻率曲線Fig．3 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in two layered model

圖4 二層層狀模型固定點源測深正演模擬各點源測深視極化率曲線Fig．4 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in two layered model

圖5 單正方柱體模型示意圖Fig．5 Diagram of single square cylinder model

從單正方柱體模型正演測深點的曲線（圖6及圖7），我們可以得到數(shù)值模擬較好地體現(xiàn)了地下有的異常體。

圖6 單正方柱體模型固定點源測深正演模擬各點源測深視電阻率曲線Fig．6 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in single square cylinder model

圖7 單正方柱體模型固定點源測深正演模擬各點源測深視極化率曲線Fig．7 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in single square cylinder model

2．3 正地形模型

2．3．1 含異常體

模型為一個正地形模型（圖8），起伏高度為4 m，背景電阻率為100Ω·m，極化率為1%，正方柱體電阻率為50Ω·m，極化率為5%，大小為10 m× 10 m中心位置為（0 m，－10 m）。

從模型正演測深點的曲線（圖9及圖10）看出，電阻率明顯受到了地形的影響，表現(xiàn)為在地形起伏段的兩端出現(xiàn)了較大的電阻率值，因為①當供電點靠近山脊時，山脊部分，電流密度相對小，測量得到電位值低，根據(jù)視電阻率公式可得，此時視電阻率偏低，呈低阻；②當供電點接近山峰時（或者說供電點進入山脊部分），此時電流密度在山脊部分相對較大，測量得到的電位值偏高，視電阻率值變高，呈高阻。但是正地形模型正演測點的極化率曲線未受到影響，這是因為極化率未受到地形影響。

圖8 正地形模型示意圖Fig．8 Diagram of positive landform model

圖9 正地形模型固定點源測深正演模擬各點源測深視電阻率曲線Fig．9 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

圖10 正地形模型固定點源測深正演模擬各點源測深視極化率曲線Fig．10 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

2．3．2 純地形

依據(jù)圖8中的凸地形，進行純地形響應(yīng)研究。

模型正演測深點曲線得出（圖11及圖12），電阻率受到了地形的影響，山脊會造成一個低阻的假異常，然而極化率卻不會受到影響。

圖11 正地形模型固定點源測深正演模擬各點源測深視電阻率曲線Fig．11 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

圖12 正地形模型固定點源測深正演模擬各點源測深視極化率曲線Fig．12 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

2．4 負地形模型

2．4．1 含異常體

模型為一個負地形模型（圖13），背景電阻率為100Ω·m，極化率為1%，正方柱體電阻率為50Ω ·m，極化率為5%，大小為10 m×10 m中心位置為（0 m，－10 m）。

圖13 負地形模型示意圖Fig．13 Diagram of negative landform model

從模型正演測深點曲線得到（圖14及圖15），電阻率也受到了地形的影響，表現(xiàn)為地形起伏段的兩端出現(xiàn)了較低的電阻率值，因為①當供電點靠近山谷時，山谷部分，電流密度相對大，測量得到電位值高，根據(jù)視電阻率公式可得，此時視電阻率偏高，呈高阻；②供電點進入山谷部分，此時電流密度相對山谷部分較小，測量得到的電位值偏低，視電阻率值變低，呈低阻。但是負地形模型正演測點的極化率曲線未受到影響。

圖14 負地形模型固定點源測深正演模擬各點源測深視電阻率曲線Fig．14 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in negative landform model

圖15 負地形模型固定點源測深正演模擬各點源測深視極化率曲線Fig．15 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in negative landform model

2．4．2 純地形

依據(jù)圖13中的負地形模型，進行純地形響應(yīng)研究。

從模型正演測深點曲線得到（圖16及圖17），電阻率受到了地形的影響，山谷會造成造成一個高阻的假異常，然而極化率缺不會受到影響。

圖16 正地形模型固定點源測深正演模擬各點源測深視電阻率曲線Fig．16 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

圖17 正地形模型固定點源測深正演模擬各點源測深視極化率曲線Fig．17 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

4 結(jié)論

本文針對固定點源測深裝置，設(shè)計了不同模型進行數(shù)值模擬，通過正演數(shù)值模擬，得到新的認識：

1）從雙層模型、單正方柱體模型正演測深點的曲線，我們可以看出數(shù)值模擬較好地體現(xiàn)了地下有的異常體。

2）從正地形的模型（含異常體）正演測深點的曲線得到，電阻率明顯受到了地形的影響，表現(xiàn)為在地形起伏段的兩端出現(xiàn)了較大的電阻率值；負地形的模型（含異常體）正演測深點曲線，電阻率也受到了地形的影響，表現(xiàn)為地形起伏段的兩端出現(xiàn)了較低的電阻率值。但正、負地形模型正演測點的極化率曲線未受到影響，這是因為極化率未受到地形影響。

3）從純地形的模型正演測深點的曲線得到，電阻率也明顯受到了地形的影響，表現(xiàn)為正地形會產(chǎn)生一個假的低阻異常，負地形會出現(xiàn)假的高阻異常；極化率未受到影響。

［1］ SMITH J T，BOOKER J R．Rapid Inversion of Twoand Three－Dimensional Magnetotelluric data［J］．Geophys Res．，1991，96：3905－3922．

［2］ 毛先進．邊界積分方程和以其為基礎(chǔ)的電阻率層析成像［D］．長沙：中南工業(yè)大學(xué)，1997．

［3］ 熊彬，阮百堯．電位雙二次變化二維地電斷面電阻率測深有限元數(shù)值模擬IJ］．地球物理學(xué)報，2002，45：285－294．

［4］ 羅延鐘，孟水良．關(guān)于用有限單元法計算二維構(gòu)造點電源場的幾個問題［J］，地球物理學(xué)報，1986，29（6）：37－45．

［5］ 強建科，羅延鐘，熊彬．固定點源測深激電異常研究［J］．地球物理學(xué)進展．2005，20（4）：1176－1183．

［6］ 徐世浙．地球物理中的有限單元法［M］．北京：科學(xué)出版社，1994．

［7］ 周熙襄，鐘本善．點源二維電法正演的有限元單元法［J］．物探化探計算技術(shù)，1983，5（3）：19－40．

［8］ 阮百堯，徐世浙．電導(dǎo)率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬［J］．地球科學(xué)，1998，23（3）：303－307．

Research about two－dimensional IP numerical simulation

CHEN Yong-lin1，JIANG Shou-jing2，XIE Dan1
（1．The 12th corps of Gold troops of CAPF，Chengdu 610059，China；2．The 10th corps of Gold troops of CAPF，Kunming 65000，China）

In this paper，starting from the question of point source and dimensional geoelectric field，we use finite element method to simulate geoelectric field，triangle subdivision algorithm to rolling surface，and various models to complete forward simulation according to the characteristic of bilateral three-pole device．By means of forward simulation and Inversion of various models，we have summarized some features about abnormity to offer some useful information for analysis of Induced polarization．

induced polarization method；finite element method；bilateral three-pole device

P 631．3＋24

A

10．3969／j．issn．1001-1749．2014．05．05

1001-1749（2014）05-0541-06

2014-02-25 改回日期：2014-07-23

陳永凌（1988－），男，碩士，主要研究方向為礦產(chǎn)地球物理，E-mail：476811938＠qq．com