基于局部加密等級網格的2.5D 直流電法有限元模擬

胡宏伶，肖曉，潘克家 ，湯井田，謝維

(1. 湖南師范大學 數學與計算機科學學院，湖南 長沙，410081;2. 中南大學 地球科學與信息物理學院，有色金屬成礦預測教育部重點實驗室，湖南 長沙，410083;3. 中南大學數學與統計學院, 湖南 長沙，410083;4. 成都理工大學 油氣藏地質及開發工程國家重點實驗室, 四川 成都，610059)

自1971 年Coggon[1]將有限單元法(finite element method，FEM)應用到直流電阻率法數值模擬以來，FEM 以其理論完備、邊界處理能力強和通用性強等優點，在直流電法正演中得到了廣泛應用[2-18]。直流電法數值模擬問題具有如下特點：地下往往存在復雜的電性不均勻結構，電性差異可達幾個數量級；研究區域具有無界性，需要進行截斷邊界處理；點電源處電位為無窮大，具有奇性。這些都將導致傳統有限元法模擬精度及效率降低。為解決上述問題，提高有限元的計算效率及點源處的模擬精度，湯井田等[19]提出一種3D 直流電阻率有限元—無限元耦合數值模擬方法，有效地降低了截斷邊界的影響。為提高源點附近的精度，陳小斌等[20]采用組合網格技巧，對源點附近進行局部加密，并利用有限元直接迭代算法求解線源頻率域大地電磁正演問題。湯井田等[21-24]采用非結構化網格的自適應有限元，通過加密奇點(點電源)處網格密度以減小有限元的計算誤差，提高了有限元模擬的速度和精度。采用局部加密的非結構化網格，能提高奇點附近有限元的模擬精度，但網格剖分及總體剛度矩陣集成繁瑣，編程較復雜，且網格剖分具有很大的不確定性。非結構化網格破壞了FEM 固有的許多超收斂結構，在提高局部精度的同時降低了FEM 的整體收斂速度[25]。另一種處理方式是采用較密的均勻剖分，并利用多網格法求解超大規模有限元方程。為解決通常幾何多重網格法粗網格難以模擬復雜電性差異的問題，Moucha 等[26]提出一種電導率的6 參數化表示法，并給出基于9 點格式有限差分的多網格法，求解2D 直流電法正演問題；魯晶津等[27]提出直流電阻率3D 正演的代數多重網格方法。潘克家等[28-29]引入外推瀑布式多網格法求解2.5D/3D 直流電阻率法數值模擬問題，大大提高了直流電法正演的計算速度。Li等[30]提出一種3D 直流電法的多源系統有限元模擬方案。Wang 等[31]給出非結構化網格下各向異性3D 直流電阻率法有限元模擬。在此，本文作者擬結合局部加密策略和結構化網格的優點，提出一種基于λ-等級網格的2.5D 直流電法有限元模擬方案。

1 電阻率法2.5D 有限元逼近

1.1 點源穩定電場邊值問題

假設地下為2D 構造，以電源點為坐標原點建立合適坐標系使z 軸平行于構造走向，結合z 方向的Fourier 變換，點源波數域中電位函數 U (x, y , k)在Ω內滿足如下2D 變系數非齊次亥姆霍茨(Helmholtz)方程邊值問題[16]：

考慮到本文將采用局部加密的等級網格，截斷邊界?！奕槿鐖D1 所示以原點為圓心、r∞為半徑的半圓，則電位u 滿足的橢圓邊值問題(1)可進一步簡化為

圖1 2D 結構模型Fig.1 Model of 2D structure

1.2 電阻率法有限元逼近

1.2.1 分塊等級網格剖分

直流電阻率法正演實際上為求解具有點源奇點的橢圓邊值問題。陳傳淼[25]認為解決奇點問題最好的辦法是局部加密網格法。

考慮半圓型求解區域Ω={(r,θ)|0＜ r＜1,0＜ θ＜π}。 首先，將Ω分為k個子扇形Ωi(i=1,2, …, k)，確保每個子扇形的圓心角αi＜π/2。然后，在每個子扇形上，以r 方向的剖分為基準，進行局部加密的三角形剖分。

設先將半徑方向作均勻剖分 rk*=kh0(k=0, 1, …,N)，步長 h0=r∞/N。為了在O 點附近局部加密，可取區間J=[0, rl*]，并在其上取m 個等級網格節點。

其中：等級網格參數λ ＞1。顯然，rm=rl*=lh0。而小區間Ij=(rj,rj+1)的長度為

重新排列剖分節點為{r0, r1, …, rm=rl*, …,rm+N-l=rN*}，并對扇形Ωi進行局部加密的三角剖分。將半徑為rj的圓弧作j 等分得到j+1 個節點，并將它們連接為折線。然后，用直線依次連接相鄰兩圓弧上的節點，得到1 個三角形網格。幾個子區域上的網格就構成了所需的分塊局部加密等級網格，見圖2。為了簡便，本文實際計算中直接取l=N，即在整個求解區域采取λ-等級網格。

1.2.2 單元分析

在單元e(逆時針編號依次為i, j, m)上進行三角形線性插值：

圖2 λ=2 時的分塊局部加密等級網格Fig.2 Piecewise locally refined graded mesh when λ=2

其中：

首先將式(3)中的積分分解成各單元e 上的積分。若單元e 內電導率為常數 σe，利用式(6)可得

其中：Ue=(ui, uj, um)T，為單元e 的u 值列向量；K1e=(kst)3×3為3 階對稱矩陣；

類似地，可得積分項

其中：

其中：

1.2.3 總剛合成

設區域Ω 剖分為n 個節點，將矩陣K1e，K2e和K3e相加可得3 階單元剛度矩陣Ke，并擴展成n 階方陣，然后疊加得

只有電源點對應位置非零(坐標原點為第1 個網格節點)；U 為未知的電位向量。

令式(15)變分為0，得到線性方程組

由式(10)，(12)及(14)知總體剛度矩陣K 為高度稀疏的對稱正定矩陣，可利用PARDISO 等直接稀疏求解器快速求解得到波數域電位 U (x , y , k)。得到波數域中對應一系列不同波數ki時的電位 U (x , y , ki)后，利用Fourier 逆變換即可得到感興趣的過點電源剖面(z=0)上的總電位：

其中：ki和gi分別為離散波數及相應的權系數。采用文獻[32]中由差分進化算法[33]計算得到離散波數，見表1。利用單源電場的電位，通過組合不同的采集裝置，即可得到不同裝置對應的視電阻率曲線，以用于電法勘探解釋。

與通常采取的矩形計算區域相比，本文選取半圓型截斷邊界，簡化了原問題(1)中截斷邊界上的混合邊界條件，使得貝塞爾函數不必代入式(13)的邊界積分進行計算，而只需計算 K0(kr∞)和 K1(kr∞)各1 次，大大簡化了剛度矩陣的計算過程。另一方面，半圓形求解區域非常適合如圖2 所示的結構化分塊等級網格剖分，且能自動對奇點(點電源)局部加密，在保證精度的同時提高計算速度。若對矩形區域采取通常的結構化局部加密策略，即x 方向和y 方向分別采取由密到疏的剖分(如對數均勻或λ-等級網格)，則會出現如圖3 所示非常扁長的矩形單元，網格質量下降，嚴重影響有限元模擬精度[25]。

圖3 矩形域上局部加密網格Fig.3 Locally refined mesh on rectangular domain

1.3 壓縮存儲格式及總剛直接集成

1.3.1 行索引存貯格式(CSR)

此種存儲格式按照行的順序，依次將對稱稀疏矩陣的上三角非零元素存放在一段連續存儲空間上。假設A 為n 階對稱方陣，m 表示其上三角部分非零元的個數，數組下標從1 開始，則這種存儲結構由以下3個數組構成：

雙精度數組a，即按行存儲A 的上三角部分的非零元素，長度為m；

整型數組Ja，即依次存儲A 中元素的列號，長度為m；

整型數組Ia，即存儲A 每行第1 個非零元在數組a 中的位置，長度為n+1。其中，Ia(n+1)=m+1，即第i行非零元個數為Ia(i+1)-Ia(i)。

圖4 所示為1 個對稱CSR 格式存儲的示例。此種存儲格式可用于Intel MKL的PARDISO直接稀疏求解器，并且存儲效率非常高。對圖3 所示的矩形網格采用矩形雙線性元，有限元離散形成的總體剛度陣每行最多5 個非零元。而對圖2 所示的等級網格，采用三角形線性元，總體剛度陣每行最多4 個非零元。

1.3.2 CSR 格式總剛直接集成

表1 離散波數及相應的權系數Table 1 Discrete wavenumbers and corresponding weights

圖4 CSR 格式存儲示例Fig.4 Example of CSR storage format

分析剛度矩陣的組裝過程得知：只有2 個節點相互關聯(出現在同一單元中)，才會對總體剛度矩陣有貢獻。在總剛K 的第i 行，非零元的列標正好為與節點i 關聯的節點編號。為此，為每個節點i 設計1 個名為“關聯節點”的數據結構，用于存儲非零元的列標。

整數Ni，為與節點i 關聯的節點個數；

整型數組Ci，為按升序依次存儲關聯節點編號。

其中：

圖5 關聯節點示意圖Fig.5 Schematic diagram of associated node

2 數值計算結果與分析

本文程序在Windows 7(Intel? CoreTMi5-760，2.8 GHz，4.0 GB 內存)下運行，軟件平臺為Intel Parallel Studio XE 2011，包含Intel Visual Fortran Compiler XE 12.0 以及Intel Math Kernel Library 10.3。

首先考慮扇形區域Ω={0＜ r ＜1,0＜ θ＜π/3}上如下拉普拉斯方程邊值問題：

其中：邊界函數g 由精確解u=r1/2sin(θ /2)確定，精確解在原點處不可導，具有奇性。采用如前所述的局部加密的等級網格，利用有限元方法求解橢圓邊值問題(21)，計算結果如表2 所示(其中：N 為徑向剖分數；誤差表示數值解與精確解之間的最大絕對誤差)。從表2 可看出，均勻網格下對此類奇異解問題，有限元只有半階收斂；網格參數λ 增加1，收斂階增加半階；固定N，隨著λ 增加，精度顯著提高；λ 取4 時已達到線性有限元最高階即2 階收斂，故沒必要進一步增大λ。另外，此例說明對此類奇異解問題，均勻剖分收斂非常緩慢，而采用局部加密的等級網格可顯著提高有限元解的精度，即使λ=2 時128 等級剖分的計算精度(4.69×10-6)已與均勻1 024 剖分時的精度(4.20×10-6)相當，而λ=4 時128 等級剖分的誤差僅為2.45×10-7，也遠比均勻1 024 剖分時的4.20×10-6??；同時也驗證了本文方法的可行性和有限元程序的正確性。

表2 等級網格有限元解收斂性分析Table 2 Convergence analysis of finite element solutions for graded meshes

需要指出的是：以上模型問題為調和方程邊值問題，沒有考慮物性的不均勻性以及附加源的影響；同時，精確解的奇性也比較弱，只在原點導數不存在，而函數仍是連續的。然而，在直流電法數值模擬中，地下存在復雜的電性不均勻結構，電性差異甚至可達幾個數量級；同時電法正演問題為無界問題，需要進行截斷邊界處理，且在點電源處電位函數本身就不連續，為無窮大。這些都將導致有限元模擬精度及效率下降，因此，不能期望其達到如表2 所示模型問題的計算精度。

基于行索引稀疏存儲模式，利用Intel MKL 的PARDISO 直接稀疏求解器求解有限元方程(17)。因為基于Cholesky分解的直接求解器計算時間主要取決于問題的規模及系數矩陣的稀疏結構，故每次求解962 001 階有限元方程的時間差不多，約為4.5 s；而進行1 次完整的正演計算(需要解5 次方程組)，總耗時約為24.5 s，由此可知正演程序其它部分(包括5 次單元剛度陣的計算及總體剛度矩陣的集成) 僅耗時2 s。這是因為局部加密的λ-等級網格為結構化網格，生成的剛度矩陣稀疏結構簡單，且每行最多7 個非零元，實際上考慮到對稱性僅需存儲4 個元素，不僅減小計算機內存，而且可大大簡化總體剛度矩陣的集成。而采取半圓形求解區域，簡化了原邊值問題(1)截斷邊界?！奚系幕旌线吔鐥l件，使得貝塞爾函數不必參與單剛的計算，提高了總剛的計算速度。

2.1 模型一[28]

2 層水平地層：第1 層電阻率為10 ?·m，厚度為10 m；第2 層電阻率為100 ?·m。采用單極供電，供電電極位于坐標原點。圖6 所示為兩極裝置得到的視電阻率曲線，其中精確解為線性濾波法計算結果[34]。圖7 所示為主剖面上的相對誤差分布。從圖7 可看出：相對誤差關于y 軸對稱，且在y=10 左右有明顯突變，反映出介質水平分界面的位置。

圖6 兩層地電模型計算結果Fig.6 Numerical results for two-layer geoelectric model

圖7 主剖面上的相對誤差Fig.7 Relative errors on main cross-section

從圖6 可以看出：3 種不同的網格剖分，在遠離電源點時(r＞3)都能夠得到非常精確的結果，相對誤差控制在0.5%以內；徑向采取均勻剖分時(λ=1)，數值解在點源附近的誤差較大，超過4%；局部加密等級網格在不增加計算工作量的前提下，顯著提高了奇點附近的精度，且2 種不同的等級網格(λ=2, 3)得到幾乎完全相同的結果。經計算，基于等級網格有限元解最大相對誤差為 0.65%，而平均相對誤差僅為0.30%。從計算時間和模擬精度2 個方面考慮，此例計算結果都與文獻[28]中采用1 600×1 600(2 563 201個節點，2 560 000 個單元)均勻矩形剖分時，外推瀑布式多網格法的計算結果相當。文獻[28]中計算時間為28 s，最大相對誤差0.87%，平均相對誤差0.22%。

2.2 模型二[5]

K 型3 層層狀地層：第1 層電阻率為50 ?·m，厚度為5 m；第2 層電阻率為100 ?·m，厚度為10 m；第3 層電阻率為20 ?·m。數值模擬結果見圖8。從圖8 可看出：采用等級網格的有限元法能在200 m 范圍內得到高精度的計算結果，平均相對誤差為0.29%；而200 m 以外相對誤差顯著增加，這是因為Fourier離散波數是基于均勻地層精確解利用200 m 內的極距序列最優化選取的，適用范圍為1～200 m[29]。除此之外，可得到與模型一完全類似的結論：(1) 局部加密等級網格可顯著提高奇點附近的模擬精度(x=1 附近相對誤差從4.60%顯著減為0.10%；(2) 網格參數λ=2和λ=3 對應的視電阻率誤差曲線幾乎重合，即對此類電法正演問題網格參數λ 取為2 即可。

2.3 模型三[21]

圖8 3 層地電模型計算結果Fig.8 Numerical results for three-layer geoelectric model

圖9 垂直電性界面模型計算結果Fig.9 Numerical results of vertical electrical interface

2D 垂直典型界面，左邊介質電阻率為1 ?·m，右邊介質電阻率為100 ?·m；單位點電源位于左邊介質上，距界面5 m?；诰植考用艿燃壘W格的有限元法得到的視電阻率及相對誤差如圖9 所示，其中圖9(b)中虛線表示分界面的位置。從圖9 可看出：采用均勻網格剖分，奇點及分界面附近的相對誤差比較大；而等級網格相對誤差在1～200 m 范圍內都比較小，控制在0.70%以內，與文獻[21]中基于非結構化網格的自適應有限元計算結果精度相當，文獻[21]中給出的誤差限為0.69%。在介質分界面處(x=5 m)，這2 種網格剖分下視電阻率的相對誤差均有體現，具有明顯的突變。值得注意的是：利用等級網格，x=1 處的相對誤差從3.94%驟降到0.12%，這進一步說明了等級網格的確能顯著提高奇點附近的精度。

3 結論

(1) 在無窮遠處截取圓形邊界，大大簡化了2.5D直流電法模擬電位滿足的橢圓邊值問題及單元剛度矩陣的計算。

(2) 局部加密等級網格能在不改變問題規模的前提下，顯著提高2.5D 直流電法正演的有限元模擬精度。

(3) 對簡單的模型問題，隨著等級網格參數λ的增大(小于4)，奇異解的計算精度顯著增大，但對直流電法正演影響不大，取為2 即可。

(4) 局部加密等級網格為結構化網格，避開了通常非結構化網格計算時繁瑣的網格剖分及總剛集成過程，有限元程序實現簡便。

(5) Intel MKL 的PARDISO 為一種高效、并行直接稀疏求解器，能在普通PC 機上5 s 內求解直流電阻率法正演有限元離散得到的100 萬階稀疏線性方程組，可用于各種地球物理正演問題的計算。

[1] Coggon J H. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1): 132-145.

[2] Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped two dimensional structures[J]. Geophysical Prospecting, 1979,27(1): 106-136.

[3] 底青云, 王妙月. 穩定電流場有限元法模擬研究[J]. 地球物理學報, 1998, 41(2): 252-260.DI Qingyun, WANG Miaoyue. The real-like 2D FEM modeling research on the field characteristics of direct electric current field[J]. Chinese Journal of Geophysics, 1998, 41(2): 252-260.

[4] LI Changwei, XIONG Bin, QIANG Jianke, et al. Multiple linear system techniques for 3D finite element method modeling of direct current resistivity[J]. Journal of Central South University,2012, 19(2): 424-432.

[5] 李長偉, 阮百堯, 呂玉增, 等. 點源場井-地電位測量三維有限元模擬[J]. 地球物理學報, 2010, 53(3): 729-736.LI Changwei, RUAN Baiyao, Lü Yuzeng, et al. Threedimensional FEM modeling of point source borehole-ground electrical potential measurements[J]. Chinese Journal of Geophysics, 2010, 53(3): 729-736.

[6] Li Y G, Spitzer K. Three-dimensional DC resistivity forward modelling using finite elements in comparison with finite-difference solutions[J]. Geophysical Journal International,2002, 151(31): 924-934.

[7] 羅延鐘, 孟永良. 關于用有限單元法對二維構造作電阻率法模擬的幾個問題[J]. 地球物理學報, 1986, 29(6): 613-621.LUO Yanzhong, MENG Yongliang. Some problems on resistivity modeling for two-dimensional structures by the finite element method[J]. Chinese Journal of Geophysics, 1986, 29(6):613-621.

[8] Mundry E. Geoelectrical model calculations for two-dimensional resistivity distributions[J]. Geophysical Prospecting, 1984, 32(1):124-131.

[9] 強建科, 羅延鐘. 三維地形直流電阻率有限元法模擬[J]. 地球物理學報, 2007, 50(5): 1606-1613.QIANG Jianke, LUO Yanzhong. The resistivity FEM numerical modeling on 32D undulating topography[J]. Chinese Journal of Geophysics, 2007, 50(5): 1606-1613.

[10] 阮百堯, 徐世浙. 電導率分塊線性變化二維地電斷面電阻率測深有限元數值模擬[J]. 地球科學: 中國地質大學學報, 1998,23(3): 303-307.RUAN Baiyao, XU Shizhe. FEM for modeling resistivity sounding on 2D geoelectric model with line variation of conductivity within each block[J]. Earth Science: Journal of China University of Geoscience, 1998, 23(3): 303-307.

[11] 阮百堯, 熊彬. 電導率連續變化的三維護電阻率測深有限元模擬[J]. 地球物理學報, 2002, 45(1): 131-138.RUAN Baiyao, XIONG Bin. A finite element modeling of three-dimensional resistivity sounding with continuous conductivity[J]. Chinese Journal of Geophysics, 2002, 45(1):131-138.

[12] TANG Jingtian, WANG Feiyan, XIAO Xiao, et al. 2.5D DC resistivity modeling considering flexibility and accuracy[J].Journal of Earth Science, 2011, 22(1): 124-130.

[13] WU Xiaoping. A 3D finite-element algorithm for DC resistivity modeling using shifted incomplete cholesky conjugate gradient method[J]. Geophysical Journal International, 2003, 154(3):947-956.

[14] 吳小平, 汪彤彤. 利用共軛梯度算法的電阻率三維有限元正演[J]. 地球物理學報, 2003, 46(3): 428-432.WU Xiaoping, WANG Tongtong. A 3D finite-element resistivity forward modeling using conjugate gradient algorithm[J].Chinese Journal of Geophysics, 2003, 46(3): 428-432.

[15] 熊彬, 阮百堯. 電位雙二次變化二維地電斷面電阻率測深有限元數值模擬[J]. 地球物理學報, 2002, 45(2): 285-295.XIONG Bin, RUAN Baiyao. A numerical simulation of 2D geoelectric section with biquadratic change of potential for resistivity sounding by the finite element method[J]. Chinese Journal of Geophysics, 2002, 45(2): 285-295.

[16] 徐世浙. 地球物理中的有限單元法[M]. 北京: 科學出版社,1994: 131-205.XU Shizhe. The finite element method in geophysics[M].Beijing: Science Press, 1994: 131-205.

[17] 徐世浙, 劉斌, 阮百堯. 電阻率法中求解異常電位的有限單元法[J]. 地球物理學報, 1994, 37(S2): 511-515.XU Shizhe, LIU Bin, RUAN Baiyao. The finite element method for solving anomalous potential for resistivity surveys[J].Chinese Journal of Geophysics, 1994, 37(S2): 511-515.

[18] Zhou B, Greenhalgh S A. Finite element three-dimensional direct current resistivity modeling: Accuracy and efficiency considerations[J]. Geophysical Journal International, 2001, 145:679-688.

[19] 湯井田, 公勁喆. 三維直流電阻率有限元-無限元耦合數值模擬[J]. 地球物理學報, 2010, 53(3): 717-728.TANG Jingtian, GONG Jinzhe. 3D DC resistivity forward modeling by finite-infinite element coupling method[J]. Chinese Journal of Geophysics, 2010, 53(3): 717-728.

[20] 陳小斌, 胡文寶. 有限元直接迭代算法及其在線源頻率域電磁響應計算中的應用[J]. 地球物理學報, 2002, 45(1): 119-130.CHEN Xiaobin, HU Wenbao. Direct iterative finite element(DIFE) algorithm and its application to electromagnetic response modeling of line current source[J]. Chinese Journal of Geophysics, 2002, 45(1): 119-130.

[21] 湯井田, 王飛燕, 任政勇. 基于非結構化網格的2.5D 直流電阻率自適應有限元數值模擬[J]. 地球物理學報, 2010, 53(3):708-716.TANG Jingtian, WANG Feiyan, REN Zhengyong. 2.5D DC resistivity modeling by adaptive finite-element method with unstructured triangulation[J]. Chinese Journal of Geophysics,2010, 53(3): 708-716.

[22] 任政勇, 湯井田. 基于局部加密非結構化網格的三維電阻率法有限元數值模擬[J]. 地球物理學報, 2009, 52(10):2627-2634.REN Zhengyong, TANG Jingtian. Finite element modeling of 3D DC resistivity using locally refined unstructured meshes[J].Chinese Journal of Geophysics, 2009, 52(10): 2627-2634.

[23] REN Zhengyong, TANG Jingtian. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method[J]. Geophysics, 2010, 75(1): 7-17.

[24] TANG Jingtian, REN Zhengyong, WANG Feiyan, et al. 3D direct current resistivity forward modeling by adaptive multigrid finite element method[J]. Journal of Central South University of Technology, 2010, 17(3): 587-592.

[25] 陳傳淼. 科學計算概論[M]. 北京: 科學出版社, 2007:120-262.CHEN Chuanmiao. Introduction to scientific computing[M].Beijing: Science Press, 2007: 120–262.

[26] Moucha R, Bailey R C. An accurate and robust multigrid algorithm for 2D forward resistivity modelling[J]. Geophysical Prospecting, 2004, 52(3): 197-212.

[27] 魯晶津, 吳小平, Spitzer K. 直流電阻率三維正演的代數多重網格方法[J]. 地球物理學報, 2010, 53(3): 700-707.LU Jingjin, WU Xiaoping, Spitzer K. Algebraic multigrid method for 3D DC resistivity modeling[J]. Chinese Journal of Geophysics, 2010, 53(3): 700-707.

[28] 潘克家, 湯井田, 胡宏伶, 等. 直流電阻率法2.5D 正演的外推瀑布式多重網格法[J]. 地球物理學報, 2012, 55(8):2769-2778.PAN Kejia, TANG Jingtian, HU Hongling, et al. Extrapolation cascadic multigrid method for 2.5D direct current resistivity modeling[J]. Chinese Journal of Geophysics, 2012, 55(8):2769-2778.

[29] PAN Kejia, TANG Jingtian. 2.5D and 3D DC resistivity modelling using an extrapolation cascadic multigrid method[J].Geophysical Journal International, 2014, 197(3): 1459-1470.

[30] LI Changwei, XIONG Bin, QIANG Jianke, et al. Multiple linear system techniques for 3D finite element method modeling of direct current resistivity[J]. Journal of Central South University,2012, 19(2): 424-432.

[31] Wang W, Wu X, Spitzer K. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids[J]. Geophysical Journal International, 2013, 193(2):734-746.

[32] 潘克家, 湯井田. 2.5D 直流電法正演中Fourier 逆變換離散波數的最優化選取[J]. 中南大學學報(自然科學版), 2013, 44(7):2819-2826.PAN Kejia, TANG Jingtian. Optimized selection of discrete wavenumbers for inverse Fourier transform in 2.5D DC resistivity modeling[J]. Journal of Central South University(Science and Technology), 2013, 44(7): 2819-2826.

[33] 潘克家, 陳華, 譚永基. 基于差分進化算法的核磁共振T2譜多指數反演[J]. 物理學報, 2008, 57(9): 5956-5961.PAN Kejia, CHEN Hua, TAN Yongji. Multi-exponential inversion of T2spectrum in NMR based on differential evolution algorithm[J]. Acta Physica Sinica, 2008, 57(9): 5956-5961.

[34] Guptasarma D, Singh B. New digital linear filters for Hankel J0and J1transforms[J]. Geophysical Prospecting, 1997, 45(5):745-762.