張錢江, 戴世坤, 陳龍偉, 強建科, 李昆, 趙東東

1 中南大學有色金屬成礦預測與地質環(huán)境監(jiān)測教育部重點實驗室, 長沙 410083 2 中南大學地球科學與信息物理學院, 長沙 410083 3 桂林理工大學地球科學學院, 桂林 541004

?

多源條件下直流電阻率法有限元三維數(shù)值模擬中一種近似邊界條件

張錢江1,2,3, 戴世坤1,2*, 陳龍偉3, 強建科1,2, 李昆1,2, 趙東東1,2

1 中南大學有色金屬成礦預測與地質環(huán)境監(jiān)測教育部重點實驗室, 長沙 410083 2 中南大學地球科學與信息物理學院, 長沙 410083 3 桂林理工大學地球科學學院, 桂林 541004

在多源直流電阻率法有限元三維數(shù)值模擬中,傳統(tǒng)混合邊界條件由于剛度矩陣與源點位置相關,無法實現(xiàn)線性方程組的快速回代求解,目前常用齊次邊界條件或無限元邊界進行替代,雖然實現(xiàn)了快速回代求解,但同時也降低了數(shù)值模擬的精度.為了實現(xiàn)快速回代求解,并確保數(shù)值模擬的計算精度,本文提出了一種近似邊界條件方法.其核心思想是將與場源位置相關的邊界系數(shù)矩陣從剛度矩陣中分離出來,使得分離后的剛度矩陣與場源位置無關.并將邊界系數(shù)矩陣與邊界處一次電場向量的乘積移到線性方程組右端源項中,當場源位置發(fā)生改變時,只需要重新計算右端源項就可以實現(xiàn)快速回代求解.理論模型數(shù)值計算表明,在水平地形條件下,本文邊界條件數(shù)值精度優(yōu)于混合邊界條件;在起伏地形條件下,與齊次邊界條件相比,本文邊界條件數(shù)值結果與混合邊界條件吻合度更高.關鍵詞 直流電阻率法; 有限單元法; 邊界條件; 數(shù)值模擬

1 引言

直流電阻率法作為地球物理勘探中最古老的勘探方法之一,早已廣泛應用于礦產(chǎn)資源勘探,水文地質調查,工程勘察等領域(徐世浙, 1994; White et al., 2001; 李金銘, 2005; Rucker et al., 2010; Loke et al., 2013).隨著勘探難度的加大和方法應用的深入,以多源激發(fā)和大面積觀測為代表的三維精細勘探成為研究熱點(Loke and Barker, 1996; Nyquist and Roth, 2005; 吳小平等, 2015),從而對多源條件下數(shù)值模擬的計算精度和計算效率提出了更高的要求.截斷邊界問題和線性方程組的求解問題是影響多源條件直流電阻率法有限單元法三維數(shù)值模擬計算精度和計算效率的兩個重要因素.

針對截斷邊界問題,目前最常用的做法是采用混合邊界條件(Dey and Morrison, 1979),該算法假設前提為: (1) 點電源位于水平地形上; (2) 模型區(qū)域Ω內(nèi)部的電性不均勻性對無窮邊界?！奚系碾娢环植疾划a(chǎn)生影響,其電位是點電源電位.因此,在均勻介質背景模型中采用混合邊界條件可以獲得較高的數(shù)值精度,但是在層狀介質背景模型中由于假設條件并不嚴格,邊界區(qū)域依然存在截斷邊界的影響.此外,由于剛度矩陣與場源位置相關,當場源位置發(fā)生改變時需要重新計算剛度矩陣,從而導致多源條件數(shù)值模擬的計算效率低下.為了解決這個問題,Zhao和Yedlin(1996)采用求取異常電位的方法(Lowry et al., 1989),對二次場電位施加Dirichlet邊界條件:ψ=0(ψ為二次場電位),由于ψ按照c/r2規(guī)律衰減,該方法具有非常高的數(shù)值精度.但是,異常電位法需要非常精確的一次場電位,因而很難應用到起伏地形模型計算中.阮百堯等(2001)認為混合邊界條件中邊界距離源點位置r較遠,可以看成常數(shù),從而使得剛度矩陣在源點移動過程中保持不變.強建科和羅延鐘(2007)在上述觀點基礎上采用Choleshy對剛度矩陣進行分解,當場源位置發(fā)生改變時,通過回代求解提升了正演數(shù)值模擬的計算效率.黃俊革等(2003)采用齊次邊界條件(Neumann boundary),認為當r→∞時,混合邊界條件中cos(r,n)U/r→0,因此,只需要選取適當?shù)倪吔缇嚯x,截斷邊界的影響就可以忽略.Blome等(2009)采用無限元法處理邊界問題,該方法與Dirichlet和混合邊界條件相比縮小了計算范圍,并減少了計算量.湯井田和公勁喆(2010)在有限元-無限元法中,提出了一種全新的最優(yōu)的無限元形函數(shù),然后將其與非結構化四面體有限元相結合,使得電位在無限域內(nèi)連續(xù)并在無限遠處衰減為零,在相同計算范圍下該方法計算精度優(yōu)于齊次邊界條件,并在近似測區(qū)大小的計算范圍內(nèi)與混合邊界條件計算精度相當.

在線性方程組求解中,多源向量求解問題的快速計算一直是研究的難點和熱點.線性方程組的求解主要包括兩大類:迭代求解法(Saad, 1996; Saad and Van der Vorst, 2000)和直接求解法(Davis, 2006).由于迭代求解法具有消耗計算資源低和在許多情況下收斂較快的優(yōu)點,在直流電阻率法三維數(shù)值模擬中占有主導地位(Dey and Morrison, 1979; Ajiz and Jennings, 1984； Spitzer, 1995; Zhang et al.,1995; Zhou and Greenhalgh, 2001; Li and Spitzer, 2002; Wu et al., 2003; 吳小平和汪彤彤, 2003).隨著計算機計算能力的飛速發(fā)展,稀疏矩陣的直接求解法開始出現(xiàn)(Liu, 1992; Schenk et al., 2001).該方法雖然需要耗費大量的時間對稀疏矩陣進行分解,并占用巨額的內(nèi)存資源,但是只需要分解一次就能實現(xiàn)多源向量的快速回代求解.因此,在多源向量求解問題中和需要大量正演計算的反演問題中具有很大的應用潛力,并已經(jīng)被應用到電磁法三維勘探中(Oldenburg et al., 2008; Arunwan and Siripunvaraporn., 2010; Yang and Oldenburg, 2012; Schwarzbach and Haber, 2013; Grayver et al., 2013; 韓波等, 2015; 楊軍等, 2015). Blome等(2009)將基于“Pardiso”的直接求解法應用到直流電阻率法有限元三維數(shù)值模擬中,并與預條件共軛梯度法(PCG)進行了對比研究,研究結果表明,在多源向量求解問題中直接解法優(yōu)勢更大.

針對多源條件直流電阻率法有限元三維數(shù)值模擬中常用邊界條件存在的問題,本文提出了一種全新的近似邊界條件方法: (1) 從電位與電場的相互關系出發(fā)給出半無限邊界上電位滿足的邊界條件; (2) 將邊界系數(shù)矩陣從剛度矩陣中分離出來,并將其與邊界處一次電場向量的乘積移到線性方程組右端源項中; (3) 采用直接解法求解器對分離后的剛度矩陣進行三角矩陣分解; (4) 當場源位置發(fā)生改變時,通過重新計算右端源項實現(xiàn)快速回代求解.文中詳細論述了近似邊界條件的計算方法,并給出了均勻背景模型和層狀介質模型中邊界處一次電場的計算表達式,最后通過理論模型驗證了本文算法的有效性.

2 有限單元法基本原理

雙點電源置于三維無限半空間Ω的表面Γs時,電位滿足微分方程(Xu, 1994)

(1)

其中σ為介質電導率,U為電位,I為點源供電電流,δ為狄拉克函數(shù),A和B分別為點源供電點所在位置,A點供正電流,B點供負電流(后文同).

在地面Γs上，電位滿足的邊界條件為

(2)

在半無限邊界Γ∞上,目前常用的邊界條件有

(3)

其中n為外法向方向,rA和rB分別為點電源A和B距離邊界的距離,f(x,y,z)為已知函數(shù),在數(shù)值模擬中通常取f(x,y,z)=0.

三維直流電阻率法雙點電源邊值問題由式(1),(2)和(3)組成.基于變分原理,可以得到與邊值問題等價的變分問題(以混合邊界條件為例)

(4)

采用四面體單元進行空間離散,離散方案為:首先將模型進行六面體單元剖分,然后在六面體中進行四面體單元的剖分(熊彬等, 2003).六面體單元剖分成六個四面體單元的方法參考文獻(Zhou and Greenhalgh, 2001),該剖分方法具有很好的數(shù)值對稱性和數(shù)值精度.在四面體單元中采用雙線性插值函數(shù):

(5)

其中Ni為插值形函數(shù),與四面體自然坐標Li相同(徐世浙, 1994),Ui和σi分別為節(jié)點上的電位值和電導率值.

對變分問題(4)各項單元分析后所得的結果相加,再擴展成全體節(jié)點組成的矩陣或列陣.由全部單元相加,得

(6)

(7)

其中K=K1+K2,為對稱、正定的大型稀疏矩陣.

對于線性方程組(7),我們設計了電阻率為100 Ωm的均勻半空間模型,對常用不完全Cholesky分解共軛梯度法(ICCG)(Wu et al., 2003)和Pardiso求解器的求解性能進行了測試,測試結果見表1.模型測試中,模型網(wǎng)格做均勻剖分,共軛梯度法迭代終止條件為均方根誤差ε≤10-5,測試計算機配置為CPU-InterXeon(R)E5-2687W,主頻3.10GHz,內(nèi)存256.00GB.

表1 不同維數(shù)線性方程組求解參數(shù)Table 1 Parameters of linear equations with different dimensions

從表1我們可以看出,單個場源的正演計算中,ICCG具有很大的優(yōu)勢,且隨著維數(shù)的增加優(yōu)勢越大.但當場源較多時,Pardiso只需要耗費非常短的時間就能完成線性方程組回代求解,考慮到多源條件下直流電阻率法三維反演的計算效率需求,我們對兩種解法在反演計算中求解線性方程組總耗時進行粗略估算,以共軛梯度法反演為例(Ellis and Oldenberg, 1994; Zhang et al., 1995;吳小平和徐果明, 2000):假設ICCG和Pardiso求解一次線性方程組總時間分別為tICCG和tPardiso, Pardiso求解中回代求解時間為tbackward,激發(fā)場源個數(shù)Ns為50個,反演迭代次數(shù)Niter為10次,每次迭代中共軛梯度內(nèi)循環(huán)Ncg為8次.則共軛梯度反演過程中兩種求解方法求解方程組總耗時計算表達式分別為

(8)

×Niter.

(9)

根據(jù)計算表達式(8)和(9)我們可以估算出兩種解法在不同維數(shù)反演計算中求解線性方程組總耗時以及求解效率比,計算結果見表2.從估算結果中我們可以看出,在多源條件直流電阻率法三維共軛梯度反演計算中,當計算機硬件滿足計算要求時,采用直接解法具有更高的計算效率.此外,在實際計算過程中當條件數(shù)較差時將大大增加迭代求解法的計算時間,而直接求解法不受這一因素的影響.

表2 兩種解法在共軛梯度反演中 求解方程組總耗時估算Table 2 Total estimated time cost in resolving linear equations with conjugate gradient inversion in the two methods

3 邊界條件算法

由式(7)可知,采用有限單元法得到的線性方程組可以表示為

K1U+K2U=P.

(10)

常用邊界條件中(如式(3)所示):Dirichlet邊界條件通常將?！奕〉米銐蜻h,近似認為U|?！?0,由于邊界區(qū)域取得足夠大,需要增加大量的計算節(jié)點,從而降低了計算效率.該方法也可以近似認為U|?！逓榻財噙吔缟弦淮螆鲭娢恢?這樣?！蘧筒槐厝〉煤艽?從而降低了計算工作量.徐世浙(1994)在位場延拓有限單元法數(shù)值模擬中詳細論述了這種方法的加載方法.

K1U=P.

(11)

混合邊界條件相比其他兩種邊界條件,具有更高的數(shù)值計算精度,但由于K2與源點位置相關,當源點位置發(fā)生變化時需要重新計算邊界項系數(shù)矩陣K2,因此在多源條件下無法實現(xiàn)線性方程組的快速回代求解.針對這個問題,阮百堯等(2001)將K2中與源點位置相關的r看成常數(shù),使得邊界項系數(shù)矩陣K2在源點移動過程中保持不變.這種假設在測線較長,源點位置移動較大時存在較大的數(shù)值誤差.

綜合上述三種邊界條件方法,針對多源條件直流電阻率法三維數(shù)值模擬的計算效率和計算精度問題,本文提出了一種近似邊界條件方法:該方法取具有解析解或數(shù)值濾波解的背景模型計算出邊界一次電場E0,并作為邊界上的值,然后將其與邊界單元積分的乘積向量移項到右端項源項.在半無限邊界Γ∞上,從電位與電場的相互關系出發(fā),電位滿足的邊界條件為

(12)

從而可得邊界積分項變分表達式為

(13)

對邊界積分項(13)進行單元分析,將所得的結果相加,再擴展成全體節(jié)點組成的矩陣或列陣,由全部邊界單元相加得

(14)

其中K3為邊界系數(shù)矩陣.

令(14)式的變分為零,并綜合線性方程組(10)可得新的線性方程組:

K1U+K3E=P.

(15)

我們知道,如果以截斷邊界上的近似場值作為邊界條件,那么截斷邊界造成的邊界反射就能得到相應的壓制.在直流電阻率法勘探中,總場為一次場和二次場的疊加,二次場由一次場激發(fā),而一次場占主導,因此將一次場代替總場作為截斷邊界上的邊界條件切實可行.取截斷邊界上電場值E近似等于一次電場值E0(邊界處一次電場計算方法見附錄),并將其與邊界系數(shù)矩陣K3的乘積向量移到右端源項中,得到如下線性方程組：

K1U=P-K3E0.

(16)

采用Pardiso_64位求解器對線性方程組(16)進行求解,當源點位置發(fā)生改變時,通過計算邊界處一次場電場E0得到新的源向量,從而實現(xiàn)多源向量的快速回代求解.

由直流電阻率法三維點源控制方程(1)可知,剛度矩陣K1中每行非零元素之和為零,存在如下等式：

K1c=0,

(17)

其中c為任意常數(shù)向量.

因此,線性方程組(16)為病態(tài)方程組,求解得到的電位U存在多解,

(18)

其中Ur為真實電位值.

這種多解現(xiàn)象不影響兩點之間的電位差計算結果,以及對視電阻率數(shù)據(jù)的反演結果.當需要計算精確電位值時,可以選取參考點消除U中的常數(shù)向量c.參考點為具有相對準確數(shù)值的點(如邊界處,參考值選為一次場電位值).

4 算例及分析

為了驗證本文邊界條件算法的計算精度,共選取了三個理論模型進行測試:均勻半空間模型、水平層狀介質模型和起伏地形模型.算例中如無特別說明x軸為南北方向,y軸為西東方向,z軸垂直向上,坐標原點位于模型地表中心點.算例中,線性方程組的求解采用Pardiso_64位求解器,測試計算機配置為CPU-Inter Core(TM) i7-4770,主頻3.40 GHz,內(nèi)存32.00 GB;本文邊界條件中電位參考點為偶極源中心點(該點電位大多數(shù)情況下理論解為0).

4.1 均勻半空間模型

設計電阻率為100 Ωm的均勻半空間模型,如圖1所示.其中,模型網(wǎng)格剖分Nx×Ny×Nz為101×101×51,網(wǎng)格總節(jié)點個數(shù)為520251;水平和垂直方向均采用均勻網(wǎng)格剖分,網(wǎng)格間距為10 m;偶極點電源沿y軸方向布設,供電電流IA=10 A,IB=-10 A,A和B坐標分別為(0,-20,0)和(0,20,0).

圖1 均勻半空間三維模型Fig.1 Homogeneous half space model

圖2 均勻半空間模型數(shù)值解和解析解對比紅色實線為混合邊界條件,粉紅色點線為齊次邊界條件,藍色虛線為本文邊界條件,綠色虛線為解析解.(a) z=0平面圖; (b) x=0 剖面圖; (c) z=0, x=0測線.Fig.2 Comparison between the numerical and analytic solutions of the half homogeneous space modelRed solid lines represent the Mixed boundary condition; Pink dotted lines represent the Neumann boundary condition; Blue dashed lines represent the new boundary condition; Green dashed lines represent the Analytical solution. (a) Horizontal plane at z=0; (b) The vertical section at x=0; (c) The observation line at z=0 and x=0.

4.2 水平層狀介質模型

如圖3所示為水平層狀介質兩層模型,第一層厚度為100 m,電阻率為100 Ωm;第二層電阻率為1000 Ωm.模型網(wǎng)格剖分Nx×Ny×Nz為101×101×51,網(wǎng)格總節(jié)點個數(shù)為520251;采用均勻網(wǎng)格剖分,間距為10 m;坐標原點位于地表中心點,供電電流IA=10 A,IB=-10 A,A和B坐標分別為(0,-20,0)和(0,20,0).

圖3 兩層層狀介質模型Fig.3 Two-layered model

在上述模型邊界各增加7個節(jié)點,按照1.5倍網(wǎng)格間距進行擴邊處理,如圖5所示為擴邊處理后數(shù)值解與解析解對比(去除擴邊區(qū)域),經(jīng)過簡單擴邊處理后,兩種邊界條件算法在邊界區(qū)域數(shù)值精度都得到了明顯的提升,混合邊界條件和本文邊界條件平均數(shù)值誤差分別降低到2.784%和0.59%,本文邊界條件仍然具有更好的數(shù)值結果.

4.3 起伏地形模型

隆起四方臺地形模型參數(shù)為:地形高50 m,長和寬各160 m,背景電阻率為100 Ωm,在地形正下方置入一個100 m×100 m×25 m的低阻體塊,距離頂界面120 m,電阻率為10 Ωm.如圖6所示為四方臺地形三維示意圖,地形位于地表正中央,如圖7所示為四方臺地形模型在x=0處的電阻率斷面圖.模型網(wǎng)格剖分Nx×Ny×Nz為115×115×58,網(wǎng)格總節(jié)點個數(shù)為767050;水平方向采用均勻網(wǎng)格剖分,網(wǎng)格間距為10 m;深度方向初始網(wǎng)格間距為5 m,隨著層數(shù)的增加網(wǎng)格間距逐漸遞增.模型邊界各選取7個網(wǎng)格節(jié)點按照相鄰網(wǎng)格間距的1.5倍進行擴邊處理,供電電流IA=10 A,IB=-10 A.

為了研究偶極點電源位于不同位置時截斷邊界對數(shù)值模擬的影響,本文共選取了3對偶極點電源(均位于地表):

第1對偶極源A和B均位于地形左邊水平地表上,其坐標分別為(0,-300,0)和(0,-240,0),數(shù)值結果見圖8a1—8a3;

第2對偶極源中A位于地形左邊水平地表上,B位于斜坡地形上,其坐標分別為(0,-170,0)和(0,-110,25),數(shù)值結果見圖8b1—8b3;

第3對偶極源A和B均位于地形頂端,其坐標分別為(0,-30,50)和(0,30,50),數(shù)值結果見圖8c1—8c3.

如圖8所示為偶極點電源位于不同位置時三種

圖4 兩層介質模型數(shù)值解和解析解對比紅色實線為混合邊界條件,藍色實線為本文邊界條件,綠色虛線為解析解.(a) z=0平面圖; (b) x=0剖面圖; (c) z=0, x=0測線.Fig.4 Comparison between the numerical and analytic solutions of the half homogeneous space modelRed solid lines represent the Mixed boundary condition; Blue dashed lines represent the new boundary condition; Green dashed lines represent the Analytical solution.(a) Horizontal plane at z=0; (b) The vertical section at x=0; (c) The observation line at z=0 and x=0.

圖5 擴邊處理后數(shù)值解與解析解對比Fig.5 Comparison between the numerical and analytic solutions for the two-layer model after boundary extensionSolid lines represent the Mixed boundary condition; Blue dashed lines represent the new boundary condition; Green dashed lines represent the Analytical solution.

圖6 四方臺地形三維示意圖Fig.6 A 3D view of the square frustum hill model

圖7 四方臺地形模型在x=0處的斷面圖Fig.7 The vertical section of the square frustum hill model at x=0

圖8 三種邊界條件數(shù)值模擬結果對比圖Fig.8 Comparison of numerical simulation results from three boundary conditions

邊界條件(混合邊界條件、齊次邊界條件和本文邊界條件)數(shù)值模擬結果在對數(shù)區(qū)間上的對比圖,圖8a1,8b1和8c1為地表電位投影平面圖,圖8a2,8b2和8c2為x=0處電位斷面圖,圖8a3,8b3和8c3為地表過源線電位測線圖.圖中紅線為混合邊界條件數(shù)值模擬結果,綠線為齊次邊界條件數(shù)值模擬結果,藍線為本文邊界條件數(shù)值模擬結果.在起伏地形條件下,雖然沒有解析解作為參考,但是與齊次邊界條件相比,混合邊界條件數(shù)值結果可靠性更高.對比圖中三種邊界條件數(shù)值結果可以發(fā)現(xiàn),在邊界區(qū)域中,本文邊界條件數(shù)值結果與混合邊界條件非常接近,而齊次邊界條件差異很大.

5 結論

本文在多源條件下直流電阻率法有限單元三維數(shù)值模擬中,針對目前常用邊界條件存在的問題,提出了一種新的近似邊界條件方法,該方法聯(lián)合直接解法實現(xiàn)了快速回代求解,大大提升了數(shù)值模擬的計算效率.在理論模型數(shù)值算例中,通過對比同等條件下混合邊界條件,齊次邊界條件和本文邊界條件在邊界區(qū)域的數(shù)值精度可以發(fā)現(xiàn):在水平地形條件下,由于能夠獲得較精確的邊界一次電場,本文邊界條件數(shù)值精度優(yōu)于混合邊界條件;在起伏地形條件下,相比齊次邊界條件,本文邊界條件與混合邊界條件數(shù)值結果更接近.鑒于該方法的特點,可以將其推廣應用于多源條件直流電阻率法2.5維或其他人工源電磁法數(shù)值模擬計算中.

附錄：截斷邊界一次場電場E0的計算

(1) 均勻背景模型

水平均勻半空間模型中,邊界M點電位表達式為

(A1)

(A2)

式中rAM和rBM分別為點電源A和B到邊界M點的距離;dxAM和dxBM為點電源A和B到邊界M點所在yz界面的垂直距離,dyAM和dyBM為點電源A和B到邊界M點所在xz界面的垂直距離,dzAM和dzBM為點電源A和B到邊界M點所在yx底界面的垂直距離.

當模型具有起伏地形時,邊界一次場電場E0仍然按照式(2)進行計算.

(2) 層狀介質模型

邊界處M點電位表達式由數(shù)值濾波解法求得,遞推表達式參考文獻(Li and Spitzer, 2002;李金銘, 2005),濾波系數(shù)參考文獻(Anderson, 1989; Guptasarma and Singh, 1997).根據(jù)電場與電位的相互關系可以求得邊界處M點電場分量,基于篇幅的限制,這里給出兩層層狀介質中邊界一次場電場E0的計算表達式:

第一層介質中

(A3)

第二層介質中

(A4)

當模型具有起伏地形時,邊界一次場電場值E0仍然根據(jù)式(A3)和(A4)計算.

致謝 感謝中國石油大學(北京)趙建國教授對本文給予的幫助,并特別對兩位匿名評審專家提出的寶貴修改意見表示感謝,最后感謝編輯們的辛苦工作.

Ajiz M A, Jennings A. 1984. A robust incomplete Choleski-conjugate gradient algorithm.Int.J.Numer.Meth.Eng., 20(5): 949-966.

Anderson W L. 1989. A hybrid fast hankel transform algorithm for electromagnetic modeling.Geophysics, 54(2): 263-266.

Arunwan T R, Siripunvaraporn W. 2010. An efficient modified hierarchical domain decomposition for two-dimensional magnetotelluric forward modeling.Geophys.J.Int., 183(2): 634-644.Blome M, Maurer H R, Schmidt K. 2009. Advances in three-dimensional geoelectric forward solver techniques.Geophys.J.Int., 176(3): 740-752.

Davis T A. 2006. Direct Methods for Sparse Linear Systems. Philadelphia: SIAM.

Dey A, Morrison H F. 1979. Resistivity modeling for arbitrarily shaped three-dimensional structures.Geophysics, 44(4): 753-780.

Ellis R G, Oldenberg D W. 1994. The pole-pole 3-D DC-resistivity inverse problem: a conjugate gradient approach.Geophys.J.Int., 119(1): 187-194.

Grayver A V, Streich R, Ritter O. 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver.Geophys.J.Int., 193(3): 1432-1446.

Guptasarma D, Singh B. 1997. New digital linear filters for Hankel J0and J1transforms.GeophysicalProspecting, 45(5): 745-762.

Han B, Hu X Y, Huang Y F, et al. 2015. 3-D frequency-domain CSEM modeling using a parallel direct solver.ChineseJ.Geophys. (in Chinese), 58(8): 2812-2826, doi: 10.6038/cjg20150816.

Huang J G, Ruan B Y, Bao G S. 2003. Finite element method for IP Modeling on 3-D Geoelectric Section.EarthScience—JournalofChinaUniversityofGeosciences(in Chinese), 28(3): 323-326.

Li J M. 2005. Geoelectric Field and Electrical Exploration (in Chinese). Beijing: Geological Publishing House, 136-212.

Li Y G, Spitzer K. 2002. Three-dimensional DC resistivity forward modelling using finite elements in comparison with finite-difference solutions.Geophys.J.Int., 151(3): 924-934.

Liu J W H. 1992. The Multifrontal Method for Sparse Matrix Solution: Theory and Practice.SIAMRev., 34(1): 82-109.

Loke M H, Barker R D. 1996. Practical techniques for 3D resistivity surveys and data inversion.GeophysicalProspecting, 44(3): 499-523.Loke M H, Chambers J E, Rucker D F, et al. 2013. Recent developments in the direct-current geoelectrical imaging method.JournalofAppliedGeophysics, 95: 135-156. Lowry T, Allen M B, Shive P N. 1989. Singularity removal: a refinement of resistivity modeling techniques.Geophysics, 54(6): 766-774.

Nyquist J E, Roth M J S. 2005. Improved 3D pole-dipole resistivity surveys using radial measurement pairs.GeophysicalResearchLetters, 32: L21416, doi: 10.1029/2005GL024153.

Oldenburg D W, Haber E, Shekhtman R. 2008. Forward modelling and inversion of multi-source TEM data. ∥ 78thAnn. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 559-563.

Qiang J K, Luo Y Z. 2007. The resistivity FEM numerical modeling on 3-D undulating topography.ChineseJ.Geophys. (in Chinese), 50(5): 1606-1613. Ruan B Y, Xiong B, Xu S Z. 2001. Finite element method for modeling resistivity sounding on 3-D geoelectric section.EarthScience—JournalofChinaUniversityofGeosciences(in Chinese), 26(1): 73-77. Rucker D F, Loke M H, Levitt M T, et al. 2010. Electrical-resistivity characterization of an industrial site using long electrodes.Geophysics, 75(4): WA95-WA104.

Saad Y. 1996. Iterative Methods for Sparse Linear Systems. Boston: PWS Pub. Co.

Saad Y, Vorst H A D. 2000. Iterative solution of linear systems in the 20thcentury.J.Comput.Appl.Math., 123(1-2): 1-33. Schenk O, G?rtner K, Fichtner W, et al. 2001. PARDISO: a high-performance serial and parallel sparse linear solver in semiconductor device simulation.FutureGenerat.Comput.Syst., 18(1): 69-78.

Schwarzbach C, Haber E. 2013. Finite element based inversion for time-harmonic electromagnetic problems.Geophys.J.Int., 193(2): 615-634. Spitzer K. 1995. A 3-D finite-difference algorithm for dc resistivity modelling using conjugate gradient methods.Geophys.J.Int., 123(3): 903-914.

Tang J T, Gong J Z. 2010. 3D DC resistivity forward modeling by finite-infinite element coupling method.ChineseJ.Geophys. (in Chinese), 53(3): 717-728, doi: 10.3969/j.issn.0001-5733.2010.03.027.

White R M S, Collins S, Denne R, et al. 2001. A new survey design for 3D IP inversion modelling at Copper hill.ExplorationGeophysics, 32(4): 152-155.

Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method.ChineseJ.Geophys. (in Chinese), 43(3): 420-427.

Wu X P, Xiao Y F, Qi C, et al. 2003. Computations of secondary potential for 3D DC resistivity modelling using an incomplete Choleski conjugate-gradient method.GeophysicalProspecting, 51(6): 567-577.

Wu X P, Wang T T. 2003. A 3-D finite-element resistivity forward modeling using conjugate gradient algorithm.ChineseJ.Geophys. (in Chinese), 46(3): 428-432. Wu X P, Liu Y, Wang W. 2015. 3D resistivity inversion incorporating topography based on unstructured meshes.ChineseJ.Geophys. (in Chinese), 58(8): 2706-2717, doi: 10.6038/cjg20150808.

Xiong B, Ruan B Y, Luo Y Z. 2003. 3-D numerical simulation study of DC resistivity anomaly under complicated terrain.GeologyandProspecting(in Chinese), 39(4): 60-64.

Xu S Z. 1994. The Finite Element Method in Geophysics (in Chinese). Beijing: Science Press, 159-172.

Yang D K, Oldenburg D W. 2012. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit.Geophysics, 77(2): B23-B34.

Yang J, Liu Y, Wu X P. 2015. 3D simulation of marine CSEM using vector finite element method on unstructured grids.ChineseJ.Geophys. (in Chinese), 58(8): 2827-2838, doi: 10.6038/cjg20150817.

Zhang J, Mackie R L, Madden T R. 1995. 3-D resistivity forward modeling and inversion using conjugate gradients.Geophysics, 60(5): 1313-1325.

Zhao S K, Yedlin M J. 1996. Some refinements on the finite-difference method for 3-D DC resistivity modeling.Geophysics, 61(5): 1301-1307.

Zhou B, Greenhalgh S A. 2001. Finite element three-dimensional direct current resistivity modelling: accuracy and efficiency considerations.Geophys.J.Int., 145(3): 679-688.

附中文參考文獻

韓波, 胡祥云, 黃一凡等. 2015. 基于并行化直接解法的頻率域可控源電磁三維正演. 地球物理學報, 58(8): 2812-2826, doi: 10.6038/cjg20150816.

黃俊革, 阮百堯, 鮑光淑. 2003. 三維地電斷面激發(fā)極化法有限元數(shù)值模擬. 地球科學—中國地質大學學報, 28(3): 323-326.

李金銘. 2005. 地電場與電法勘探. 北京: 地質出版社, 136-212.

強建科, 羅延鐘. 2007. 三維地形直流電阻率有限元法模擬. 地球物理學報, 50(5): 1606-1613.

阮百堯, 熊彬, 徐世浙. 2001. 三維地電斷面電阻率測深有限元數(shù)值模擬. 地球科學—中國地質大學學報, 26(1): 73-77.

湯井田, 公勁喆. 2010. 三維直流電阻率有限元-無限元耦合數(shù)值模擬. 地球物理學報, 53(3): 717-728, doi: 10.3969/j.issn.0001-5733.2010.03.027.

吳小平, 徐果明. 2000. 利用共軛梯度法的電阻率三維反演研究. 地球物理學報, 43(3): 420-427.

吳小平, 汪彤彤. 2003. 利用共軛梯度算法的電阻率三維有限元正演. 地球物理學報, 46(3): 428-432.

吳小平, 劉洋, 王威. 2015. 基于非結構網(wǎng)格的電阻率三維帶地形反演. 地球物理學報, 58(8): 2706-2717, doi: 10.6038/cjg20150808.

熊彬, 阮百堯, 羅延鐘. 2003. 復雜地形條件下直流電阻率異常三維數(shù)值模擬研究. 地質與勘探, 39(4): 60-64.

徐世浙. 1994. 地球物理中的有限單元法. 北京: 科學出版社, 159-172.

楊軍, 劉穎, 吳小平. 2015. 海洋可控源電磁三維非結構矢量有限元數(shù)值模擬. 地球物理學報, 58(8): 2827-2838, doi: 10.6038/cjg20150817.

(本文編輯 胡素芳)

An approximate boundary condition for FEM-based 3-D numerical simulation with multi-source direct current resistivity method

ZHANG Qian-Jiang1,2,3, DAI Shi-Kun1,2*, CHEN Long-Wei3, QIANG Jian-Ke1,2, LI Kun1,2, ZHAO Dong-Dong1,2

1KeyLaboratoryofMetallogenicPredictionofNonferrousMetalsandGeologicalEnvironmentMonitoring,MinistryofEducation,CentralSouthUniversity,Changsha410083,China2SchoolofGeosciencesandInfo-physics,CentralSounthUniversity，Changsha410083,China3CollegeofEarthSciences,GuilinUniversityofTechnology,Guilin541004,China

In finite element-based three-dimensional numerical simulation of direct current (DC) resistivity method, applying traditional mixed boundary condition cannot accomplish solving linear equations through recursion due to computation of stiffness coupling with source positions. Currently Neumann boundary condition or infinite element boundary condition is usually used instead. Although rapidly resolving, these two methods reduce the precision of numerical simulation. To realize recursive resolving rapidly and ensure simulation precision, an approximate boundary condition is proposed to implement both fast recursive resolving and precise simulation. The key ideas underlying the method are to separate boundary coefficient matrix coupled with source positions from stiffness matrix so as to make the resultant stiffness matrix independent of source positions. Production of boundary coefficient matrix and primary electric field vectors on the boundary is then transferred to the right hand of linear equations. In doing so only right-hand source items need to be computed when source positions are changed. Synthetic tests show that the numerical simulation precision applying the newly proposed boundary condition is superior to the one using mixed boundary condition in the case of horizontal topography. Also, in the case of rugged topography, the simulation results, compared with the application of neumann boundary, are much closer to those with mixed boundary condition.

Direct current resistivity method; Finite element method; Boundary condition; Numerical simulation

10.6038/cjg20160927.

國家重大科研儀器設備研制專項資助項目(41227803)和國家自然科學基金項目(41574127,41174104,41674075)聯(lián)合資助.

張錢江,男,1985年生,博士研究生,主要從事地球物理數(shù)值模擬與反演成像研究.E-mail:qjz2013@csu.edu.cn

10.6038/cjg20160927

P631

2015-12-08,2016-04-14收修定稿

張錢江, 戴世坤, 陳龍偉等. 2016. 多源條件下直流電阻率法有限元三維數(shù)值模擬中一種近似邊界條件. 地球物理學報，59(9):3448-3458,

Zhang Q J, Dai S K, Chen L W, et al. 2016. An approximate boundary condition for FEM-based 3-D numerical simulation with multi-source direct current resistivity method.ChineseJ.Geophys. (in Chinese),59(9):3448-3458,doi:10.6038/cjg20160927.

*通訊作者 戴世坤,男,1964年生,教授,博士生導師,研究方向為電法勘探三維高效反演成像算法及其實用化軟件開發(fā).

E-mail:dskgmes@csu,edu,cn