激電二維數值模擬研究

陳永凌，蔣首進，謝 丹

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

激電二維數值模擬研究

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

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

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

激發極化法；有限元法；雙邊三極

0 引言

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

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

1 基本原理

1．1 點源二維地電問題

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

這時候電位的基本微分式方程變為

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

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

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

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

1．2 穩定電流場的邊值問題

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

邊界條件同常分為三類：

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

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

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

由于狄里希萊條件會產生計算值比解析解的值小；諾依曼條件會產生計算值比解析解的值大。混合邊界條件具有能在遠處邊界面上保持電位的物理特性，且又不需要首先段設一個“正常”的電場分布的優點，這樣可以減少邊界附近放延展網格的數量，從而減少了計算的工作量。同時因為在邊界上保持電位的連續，沒有電性的突變，對地下邊界的反射作用也消除。由于混合邊界條件具有較高的模擬精度，較好的反演效果，所以一般情況下在計算過程中采用混合邊界條件。

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

從上面的推導中得到了點源二維目標函數

的極小值相對應。

其中：φ為目標函數；D為研究區域；Γ為研究區域D的邊界。

1．4 網格剖分示意圖

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

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

網格剖分的原則一般來講有以下幾點：

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

2）每個單元內介質的電導率為常數，在不同電導率介質的分界面或者求解區的邊界上，使用三角形的一個邊互相連接，以模擬邊界線的近似折線取代邊界線。

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

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

圖1 網格剖分示意圖Fig．1 The grid subdivision schemes

2 模型計算

本節主要對固定點源測深采用了有限元進行了正演模擬，就多個模型實例進行數值模擬的結果分析總結。由于考慮到數據量的問題，因此在本章模型計算后成正演曲線時，取了其中的測深供電點中的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），我們可以看出數值模擬較好地體現了地下有的異常體。

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），我們可以得到數值模擬較好地體現了地下有的異常體。

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

圖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 純地形

依據圖8中的凸地形，進行純地形響應研究。

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

圖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 純地形

依據圖13中的負地形模型，進行純地形響應研究。

從模型正演測深點曲線得到（圖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 結論

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

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

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

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

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

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

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

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

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

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

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

［8］ 阮百堯，徐世浙．電導率分塊線性變化二維地電斷面電阻率測深有限元數值模擬［J］．地球科學，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－），男，碩士，主要研究方向為礦產地球物理，E-mail：476811938＠qq．com