三維電場直流電阻率邊界單元法地形改正

段長生，王緒本，呂玉增，李長偉，黃修保

（1．成都理工大學，成都 610059；2．桂林理工大學，桂林 541004；3．贛中南地質礦產勘查研究院，南昌 330029）

三維電場直流電阻率邊界單元法地形改正

段長生1，3，王緒本1，呂玉增2，李長偉2，黃修保3

（1．成都理工大學，成都 610059；2．桂林理工大學，桂林 541004；3．贛中南地質礦產勘查研究院，南昌 330029）

電法觀測易受到地形的影響，給資料解釋造成困難。這里對數值計算中地面網格節點的立體角（ω）公式進行了嚴格地推導，利用邊界單元法對三維直流電場的地形影響進行了計算并編寫了程序。90°角域山脊地形算例表明，解析解與數值解的電位及視電阻率相對誤差小于5%，驗證了算法及程序的正確性。實例計算了3D山谷地形2D觀測剖面的地形效應，當旁側剖面與主剖面距離大于24 m時，地形效應最大值小于5%。

邊界單元法；地形改正；立體角；球面三角形

0 引言

20世紀50年代以來，電法勘探在礦產資源勘查和工程勘察等方面的應用越來越廣泛，但地形起伏較大的地區，不僅采集難度加大，資料畸變也嚴重，給解釋造成困難，因此地形影響引起了有關學者們的重視［1－4］。起伏地形不但可以引起假異常，還會掩蓋地下由礦體或目標物引起的真異常，如果不能正確地認識和消除，就會導致完全錯誤的解釋結果。

對起伏地形資料的解釋，人們一般采用經驗判斷、帶地形直接反演和比值法消除地形的影響。比值法是用實測視電阻率除以“地形校正因子（地形模型視電阻率／地形模型電阻率）”，從而得到改正后的平地形下的視電阻率。比值法中對地形校正因子的計算的方法很多，早期的有基于物理實驗的定性解釋法；基于數學公式的解析計算方法（比如保角變換法［5］）。但由于定性解釋的局限性以及數學公式的復雜性，這些方法的應用受到一定限制。隨著計算機技術的發展和地球物理解釋方法的進步，眾多學者研究了適合于起伏地形的2D、3D數值模擬和解釋技術，從而使地形影響問題的模擬、識別和校正成為可能［6］。主要的數值計算方法有：積分方程法、邊界元方法、有限差分法和有限元方法等。其中的邊界單元法利用格林函數將三維邊值問題轉換成二維計算，計算量減少、計算速度提高，有利于大型數值計算以及快速反演的實現。

作者從區域邊界點立體角的定義出發，依據計算球面三角形面積的原理，利用地面網格點的三維坐標給出了該地形改正公式中的ω（立體角）的計算方法［7］。由于該法僅利用已有網格節點坐標，并不局限剖分單元的形式，因而可以用于多種剖分方式下節點的立體角計算。據此作者編寫了地形改正程序，給出了90°角域山脊地形算例，計算結果表明，解析解與數值解的電位及視電阻率相對誤差小于5%，驗證了算法及程序的正確性。實例計算了3D山谷地形2D剖面的地形效應，在設計模型下，當測線與主剖面的距離小于10 m時，視電阻率受地形影響仍然較大，中心值大于1．1Ω·m，而當距離大于24 m時，中心值小于1．05Ω·m，地形效應最大值小于5%，幾乎可以忽略。

1 邊界元模擬法

1．1 邊值問題

假定地下介質是均勻的，電阻率為1Ω·m，地表A點電流強度為1 A，則電位滿足的基本方程為［8］

其中：ωA是電源點A對地下區域所張的立體角；δ是狄拉克函數。

1．2 邊界積分方程

其中：ωp是計算點p的對地下區域所張的立體角；up、u0p分別是計算點p的電位、正常電位；是與計算單元的法向向量的夾角。

1．3 邊界元法

1．3．1 單元積分

式（2）中第二項利用邊界單元法求解，首先用三角形單元將計算區域劃分為m個小區域，單元內采用線性插值則

其中：l是三角元頂點；Nl是形函數；xl、yl、zl是頂點坐標，

1．4 立體角ω的計算方法

根據式（3）求解方程組時還需要確定各網格節點的立體角。對于不光滑邊界上的孤立點p，立體角ωp定義為：以p點為球心，作一半徑為ε（無窮小）的球體，若區域內球體表面積s，則ωp＝s／ε2［1］。

在水平面（圖1中紅線）上，x、y將整個區域按照矩形剖分后再剖分為三角形，x方向剖分nx＿n次。其中一內部節點p的剖分情況如圖1所示，第n、n－nx＿n＋1、n＋1點的坐標分別為（xn，yn，zn）、（xn－nx＿n＋1，yn－nx＿n＋1，zn－nx＿n＋1）、（xn＋1，yn＋1，zn＋1）。對第n點做一半徑為ε（ε趨于無窮小）的球體，此時該點的立體角由地面以下區域內球體的表面積表示。按圖1的方法將地面以下表面積的計算分為八個部分（現以圖1中第二部分的計算為例），則在水平面上，n－nx＿n＋1、n＋1點到n點及n－nx＿n＋1點到n＋1點的距離分別為a1＝

圖1 內部節點網格示意圖Fig．1 Schematic diagram for inner node mesh

待求表面積為圖2藍線區域部分，可由θ1、θ2、α、ε來表示。

圖2 球面三角形面積計算示意圖Fig．2 Schematic diagram for Spherical triangle area calculation

同樣，可以求出其他七個部分的表面積，疊加就得到整個的表面積。從而易得p點的立體角為

2 地形校正

利用式（3）可求得網格節點電位

其中：K為裝置系數；ΔU為電位差；I為供電電流；ρT五數值計算視電阻率。

當I＝1A，大地電阻率ρ＝1Ω·m，計算得到的ρT在數值上等于校正因子αT，若實測視電阻率為ρs，則經比值法地改后的視電阻率

3 算例

3．1 算例1

計算模型為驗證模型，即有解析解的90°角域山脊地形，如圖3所示，模型的電阻率為1Ω·m。由于y坐標具有對稱性，1 A正負點源供電點的x、z坐標分別為：Axz＝（－40 m、－40 m），Bxz＝（40 m、－40 m）。在供電點中間范圍x＝－20．5 m到x＝20．5 m測量。數值解與解析解電位對比如圖4所示，數值解與解析解曲線形態一致，數值相近，表1將數值解與解析解電位和視電阻率值進行了統計，21對電位和視電阻率值數值和解析計算結果中，相對誤差均較小，其中電位最大相對誤差為3．87%，平均相對誤差為1．73%；電阻率最大相對誤差為2．34%，平均相對誤差為0．64%。

圖3 90°角域山脊地形Fig．3 90°ridge topography model

圖4 數值解與解析解電位對比圖Fig．4 Comparison between analytical and numerical solutions

3．2 算例2

模型為地形效應特征分析模型，電阻率為1 Ω·m的山谷地形，中心坐標為（0 m，0 m，－4．631 m）。x、y坐標范圍均為－40 m～40 m。沿y方向除主剖面外布置了5條線，距離主剖面依次為5 m、10 m、24 m、30 m、45 m，邊界單元法計算的視電阻率異常曲線圖示于圖5（a），圖5（b）是相應的地形示意。由圖5可見，總體來講，距離主剖面越遠地形效應越小，當測線與主剖面的距離小于10 m時，視電阻率受地形影響仍然較大，中心值大于1．1Ω· m，而當距離大于24 m時，中心值小于1．05Ω·m，地形效應最大值小于5%，幾乎可以忽略。

3．3 算例3

模型為地改效果模型，凹地形下方存在一個低阻異常球體，球體半徑4 m，埋深12 m，異常球體電阻率為0．05Ω·m，區域內電阻率為1Ω·m，如圖6（a）所示，圖中只畫出了主剖面的地形，實際地形是繞對稱軸旋轉面成的凹陷。在主剖面x坐標為－39 m，39 m處分別供正負點電源。首先用有限單元法計算模型主剖面中間梯度法的視電阻率異常，如圖6（b）；然后用邊界單元法計算由地形引起的視電阻率異常，如圖6（c）；最后用比值法進行地形影響改正，并將改正結果與水平大地有限元計算的結果作比較，如圖6（d）所示。顯然圖6（b）中的正異常是由于地形引起的，直接解釋將會得到錯誤的結果。圖6（d）表明二者相差較小，經計算的平均誤差僅為2．721%。

表1 數值解與解析解電位和視電阻率值統計表Tab．1 Statistics values of the potential and apparent resistivity of numerical solution and the analytical solution

圖5 地形及邊界單元法視電阻率異常曲線圖Fig．5 Topography and abnormal curve of apparent resistivity boundary element method

圖6 模型及數值計算對比圖Fig．6 Model and the comparison chart of numerical calculation

網格大小為61×27，在Intel（R）Core（TM）2 Duo CPU，2．66 GHz，3．00 GB內存，Windows XP下運行時間不到5 s，數值計算效率較高。

4 結語

作者從點源場的邊值問題出發，通過格林公式建立邊界積分方程。利用區域剖分網格節點的坐標，采用求解球面三角形面積的方法，推導出了求解邊界單元法地形改正公式中立體角的計算方法，將此方法運用到電阻率三維邊界單元法數值模擬中，實現了對積分方程的求解。算例說明該方法誤差較小，且由于邊界元法有類似于降維的效果，可快速處理較大的網格數據，因此可作為數值模擬理論計算或者實際生產中消除地形效應的一種方法。

［1］ 湯洪志，劉慶成，龔育齡，等．邊界單元法在高密度電阻率法二維地形改正中的應用效果［J］．物探與化探，2001，25（6）：457－479．

［2］ 陳伯舫，侯作中，范國華，等．有限差分法計算三維地形影響的電磁感應［J］．地震學報，1998，20（5）：541－544．

［3］ 黃蘭珍，田憲漠．電阻率法地形改正及其在工程地質勘查總的應用［J］．物探化探計算技術，1997，19（3）：238－241．

［4］ 李正容，何繼善，溫佩琳．電阻率和CSAMT法地形改正［D］．長沙：中南大學，1991．

［5］ 徐世浙，王柏鈞．保角變換在電法勘探中的應用［M］．北京：地質出版社，1977．

［6］ 劉樹才，何昭友，劉志新，等．適合地形起伏的二維有限單元數值模擬技術［J］．物化探計算技術，2005，27（2）：131－134．

［7］ 《數學手冊》編寫組．數學手冊［M］．北京：人民教育出版社，1979．

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

Topographic correction of 3D electric field in DC resistivity by boundary element method

DUAN Chang-sheng1，2，WANG Xu-ben1，2，LV Yu-zeng3，LI Chang-wei3，HUANG Xiu-bao2
（1．Chengdu University of Technology，Chengdu 610059，China；2．Guilin University of Technology，Guilin 541004，China；3．Institute of Geological and Mineral Exploration in Central and South Jiangxi Province，Nanchang 330029，China）

It is difficult to interpret DC electrical prospecting data which were easily Influenced by the terrain．In this paper，the solid angle of the ground grid nodes for the numerical calculation of formula were rigorously deduced，3D topography effect of the DC electric field were calculated by the boundary element method using fortran Language．The results of ridge terrain of 900 angular domain show the relative error of potential and apparent resistivity are less than 5 percent and verify the correctness of the algorithm and program．3D topography effect were analyzed in valley terrain by 2D sections．In the design mode，when the lateral distance from the cross section to the main cross－section greater than 24m，the maximum value of the topographic effect less than 5%．

boundary element method；topographic correction；solid angle；spherical triangle

P 631．3＋22

A

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

1001-1749（2014）05-0529-06

2014-05-12 改回日期：2014-07-25

國家高技術研究發展計劃（863計劃）項目（2009AA06Z108）；廣西自然科學基金項目（2013GXNSFAA019277）

段長生（1974－），男，博士，主要從事電磁法正演模擬和反演成像，E-mail：49119428＠qq．com。