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

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

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

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

1.1 點源穩(wěn)定電場邊值問題

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

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

圖1 2D 結(jié)構(gòu)模型Fig.1 Model of 2D structure

1.2 電阻率法有限元逼近

1.2.1 分塊等級網(wǎng)格剖分

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

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

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

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

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

1.2.2 單元分析

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

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

其中：

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

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

類似地，可得積分項

其中：

其中：

1.2.3 總剛合成

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

只有電源點對應(yīng)位置非零(坐標(biāo)原點為第1 個網(wǎng)格節(jié)點)；U 為未知的電位向量。

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

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

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

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

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

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

1.3.1 行索引存貯格式(CSR)

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

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

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

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

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

1.3.2 CSR 格式總剛直接集成

表1 離散波數(shù)及相應(yīng)的權(quán)系數(shù)Table 1 Discrete wavenumbers and corresponding weights

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

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

整數(shù)Ni，為與節(jié)點i 關(guān)聯(lián)的節(jié)點個數(shù)；

整型數(shù)組Ci，為按升序依次存儲關(guān)聯(lián)節(jié)點編號。

其中：

圖5 關(guān)聯(lián)節(jié)點示意圖Fig.5 Schematic diagram of associated node

2 數(shù)值計算結(jié)果與分析

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

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

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

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

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

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

2.1 模型一[28]

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

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

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

從圖6 可以看出：3 種不同的網(wǎng)格剖分，在遠離電源點時(r＞3)都能夠得到非常精確的結(jié)果，相對誤差控制在0.5%以內(nèi)；徑向采取均勻剖分時(λ=1)，數(shù)值解在點源附近的誤差較大，超過4%；局部加密等級網(wǎng)格在不增加計算工作量的前提下，顯著提高了奇點附近的精度，且2 種不同的等級網(wǎng)格(λ=2, 3)得到幾乎完全相同的結(jié)果。經(jīng)計算，基于等級網(wǎng)格有限元解最大相對誤差為 0.65%，而平均相對誤差僅為0.30%。從計算時間和模擬精度2 個方面考慮，此例計算結(jié)果都與文獻[28]中采用1 600×1 600(2 563 201個節(jié)點，2 560 000 個單元)均勻矩形剖分時，外推瀑布式多網(wǎng)格法的計算結(jié)果相當(dāng)。文獻[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。數(shù)值模擬結(jié)果見圖8。從圖8 可看出：采用等級網(wǎng)格的有限元法能在200 m 范圍內(nèi)得到高精度的計算結(jié)果，平均相對誤差為0.29%；而200 m 以外相對誤差顯著增加，這是因為Fourier離散波數(shù)是基于均勻地層精確解利用200 m 內(nèi)的極距序列最優(yōu)化選取的，適用范圍為1～200 m[29]。除此之外，可得到與模型一完全類似的結(jié)論：(1) 局部加密等級網(wǎng)格可顯著提高奇點附近的模擬精度(x=1 附近相對誤差從4.60%顯著減為0.10%；(2) 網(wǎng)格參數(shù)λ=2和λ=3 對應(yīng)的視電阻率誤差曲線幾乎重合，即對此類電法正演問題網(wǎng)格參數(shù)λ 取為2 即可。

2.3 模型三[21]

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

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

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

3 結(jié)論

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

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

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

(4) 局部加密等級網(wǎng)格為結(jié)構(gòu)化網(wǎng)格，避開了通常非結(jié)構(gòu)化網(wǎng)格計算時繁瑣的網(wǎng)格剖分及總剛集成過程，有限元程序?qū)崿F(xiàn)簡便。

(5) Intel MKL 的PARDISO 為一種高效、并行直接稀疏求解器，能在普通PC 機上5 s 內(nèi)求解直流電阻率法正演有限元離散得到的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] 底青云, 王妙月. 穩(wěn)定電流場有限元法模擬研究[J]. 地球物理學(xué)報, 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]. 地球物理學(xué)報, 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] 羅延鐘, 孟永良. 關(guān)于用有限單元法對二維構(gòu)造作電阻率法模擬的幾個問題[J]. 地球物理學(xué)報, 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]. 地球物理學(xué)報, 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] 阮百堯, 徐世浙. 電導(dǎo)率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬[J]. 地球科學(xué): 中國地質(zhì)大學(xué)學(xué)報, 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] 阮百堯, 熊彬. 電導(dǎo)率連續(xù)變化的三維護電阻率測深有限元模擬[J]. 地球物理學(xué)報, 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]. 地球物理學(xué)報, 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] 熊彬, 阮百堯. 電位雙二次變化二維地電斷面電阻率測深有限元數(shù)值模擬[J]. 地球物理學(xué)報, 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]. 北京: 科學(xué)出版社,1994: 131-205.XU Shizhe. The finite element method in geophysics[M].Beijing: Science Press, 1994: 131-205.

[17] 徐世浙, 劉斌, 阮百堯. 電阻率法中求解異常電位的有限單元法[J]. 地球物理學(xué)報, 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] 湯井田, 公勁喆. 三維直流電阻率有限元-無限元耦合數(shù)值模擬[J]. 地球物理學(xué)報, 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] 陳小斌, 胡文寶. 有限元直接迭代算法及其在線源頻率域電磁響應(yīng)計算中的應(yīng)用[J]. 地球物理學(xué)報, 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] 湯井田, 王飛燕, 任政勇. 基于非結(jié)構(gòu)化網(wǎng)格的2.5D 直流電阻率自適應(yīng)有限元數(shù)值模擬[J]. 地球物理學(xué)報, 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] 任政勇, 湯井田. 基于局部加密非結(jié)構(gòu)化網(wǎng)格的三維電阻率法有限元數(shù)值模擬[J]. 地球物理學(xué)報, 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] 陳傳淼. 科學(xué)計算概論[M]. 北京: 科學(xué)出版社, 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. 直流電阻率三維正演的代數(shù)多重網(wǎng)格方法[J]. 地球物理學(xué)報, 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 正演的外推瀑布式多重網(wǎng)格法[J]. 地球物理學(xué)報, 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 逆變換離散波數(shù)的最優(yōu)化選取[J]. 中南大學(xué)學(xué)報(自然科學(xué)版), 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譜多指數(shù)反演[J]. 物理學(xué)報, 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.