各向異性介質(zhì)空間-波數(shù)混合域直流電阻率法三維數(shù)值模擬研究

戴世坤， 凌嘉宣*， 陳輕蕊， 張瑩， 李昆

1 中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室， 長沙 410083 2 中南大學(xué)有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室， 長沙 410083 3 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長沙 410083 4 西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院， 成都 610500

0 引言

直流電阻率法作為地球物理勘探中較為成熟的方法，被廣泛應(yīng)用(Ren et al., 2018a)，但其在數(shù)值模擬研究中，由于場源奇異性、大差異電導(dǎo)率對比、復(fù)雜幾何結(jié)構(gòu)、非均勻性、各向異性和大尺度規(guī)模等問題對算法和計算機體系結(jié)構(gòu)提出了更高的要求，為了解決這些問題，多年來很多學(xué)者進行了探索(Wait，1990；Al Hagrey，1994；阮百堯和熊彬，2002；熊彬和阮百堯，2002；Li and Spitzer，2002；沈松金和郭乃川，2008；沈松金等，2009；陳桂波等，2009；王威和吳小平，2010；王威，2013；Wang et al.，2013；Yin et al.，2016；Ren et al.，2018b；殷長春等，2018；張錢江等，2016；任政勇等，2018).目前，針對地下介質(zhì)為各向異性結(jié)構(gòu)的直流電阻率法數(shù)值方法研究包括積分方程法(Eskola and Hongisto，1981，1997；Eskola，1988，1992；Eloranta，1988；Flykt et al.，1996；Li and Uren，1997；陳桂波等，2009)，有限差分法(Dey and Morrison，1979；Lowry et al.，1989；Spitzer，1995；Zhao and Yedlin，1996；Wang and Fang，2001；Hou et al.，2006；Chen et al.，2017)，有限單元法(Coggon，1971；Bibby，1978；Pridmore et al.，1981；徐世浙，1988；Zhou and Greenhalgh， 2001; Zhou et al.，2009；熊彬和阮百堯，2002；Li and Spitzer，2002，2005；王威和吳小平，2010；王威，2013；Wang et al.，2013；Yan et al.，2016；Yang et al.，2018；任政勇等，2018).相對于其他兩種方法，有限單元法因有完善的理論和可用非結(jié)構(gòu)化網(wǎng)格剖分的特性，更適合模擬地形和地下復(fù)雜結(jié)構(gòu)介質(zhì)，因此被深入研究和廣泛應(yīng)用，這為解決復(fù)雜各向異性介質(zhì)模型下三維數(shù)值模擬問題提供巨大動力(Wang et al.，2013).為了探索和分析地下各向異性介質(zhì)的特性，近年來部分學(xué)者采用基于有限單元法的直流電阻率法對地下三維各向異性介質(zhì)進行了研究.Li 和 Spitzer(2005)推導(dǎo)了各向異性下基于總場和二次場的邊值條件，并利用混合邊界條件實現(xiàn)了地下介質(zhì)為不同各向異性形態(tài)的三維數(shù)值模擬；Zhou等(2009)利用基于“Gaussian quadrature grid”方法的有限單元法實現(xiàn)各向異性介質(zhì)的響應(yīng)計算，該方法不需要與表面地形匹配的恒定單元的網(wǎng)格就取得較好的效果；王威和吳小平(2010)利用二次場方法和SSOR預(yù)條件共軛梯度算法實現(xiàn)了電阻率任意三維各向異性模型的地電響應(yīng)；Wang等(2011, 2013)采用非結(jié)構(gòu)化網(wǎng)格有限單元法實現(xiàn)了任意各向異性三維直流電阻率法的三維數(shù)值模擬；Yang等(2018)采用非結(jié)構(gòu)化網(wǎng)格自適應(yīng)有限單元法模擬任意各向異性介質(zhì)的響應(yīng)，該方法采用基于梯度恢復(fù)理論的局部自適應(yīng)網(wǎng)格細化技術(shù)，結(jié)合圓形掃描直流測量技術(shù)，獲得了更多關(guān)于地球不同深度的電阻率和電各向異性特性的信息.Ren等(2018c)在考慮地形和各向異性介質(zhì)條件下，利用面向目標(biāo)的自適應(yīng)有限單元法實現(xiàn)了三維直流電法的數(shù)值模擬.

前人關(guān)于直流電阻率法三維有限單元法的數(shù)值模擬研究中，大都直接在空間域中將三維模型進行剖分，然后求解滿足邊值問題的大型方程組，這樣求解方程組占用內(nèi)存大，計算時間較長.本文采用一種空間波數(shù)混合域方法對直流電阻率法進行三維數(shù)值模擬，具體為對異常電位滿足的微偏微分方程沿水平方向進行二維傅里葉變換，令水平方向轉(zhuǎn)換成波數(shù)域，保留垂直方向為空間域，由此可根據(jù)地下介質(zhì)電流密度變化快慢靈活剖分網(wǎng)格(電流密度變化劇烈區(qū)域采用稠密剖分方式，電流密度變化緩慢區(qū)域采用稀疏剖分方式).這樣將空間域異常電位滿足的三維偏微分方程轉(zhuǎn)化成不同波數(shù)滿足的一維常微分方程，把一個大規(guī)模三維數(shù)值模擬問題分解為多個一維數(shù)值模擬問題，并在空間波數(shù)混合域中利用一維有限單元法求解方程組，大大減小內(nèi)存使用量和計算時間.利用壓縮算子進行迭代，最后獲得高精度近似解.該方法通過利用二維傅里葉的快速性和追趕法(Boisvert, 1991)求解常微分方程的高效性，實現(xiàn)了直流電阻率法三維數(shù)值模擬，這為研究和分析地下介質(zhì)各向異性特性提供一種新途徑.

1 方法原理

1.1 控制方程

將點電源置于各向異性介質(zhì)中，可得電位滿足的微分方程(徐世浙，1994)

(1)

(2)

(3)

其中

(4)

矩陣D是與三個歐拉角有關(guān)的變換矩陣，具體表達式為

(5)

(6)

根據(jù)疊加原理，總電位U可以分解為

U=Ua+Ub,(7)

將式(6)和式(7)代入式(1)，可得

-2Iδ(r-r0),(8)

正常電位Ub滿足

(9)

將式(9)代入式(8)并整理，得到異常電位Ua滿足的微分方程為

(10)

散射密度ja滿足

(11)

式中E為總電場，式(10)可表示為

(12)

(12)式即為各向異性介質(zhì)空間域異常電位Ua滿足的控制方程.將式(12)展開，可得

(13)

(14)

(1)令總電場E(0)的初始值等于背景電場；

(3)對波數(shù)域異常電場和電位進行二維傅里葉變換得到空間域異常電場和異常電位，然后計算總電場和總電位；

(4)引入壓縮算子對總電場進行修改，得到新的總電場E(1)；

(5)判斷前后兩次總電場的相對殘差是否滿足期望的數(shù)值精度.若滿足精度要求，輸出結(jié)果，否則取E(0)=E(1)，重復(fù)步驟(2)—(5)，直至滿足收斂條件.

1.2 邊界條件

為了確定控制方程式(14)的定解問題，還需給出z向上合適的邊界條件.設(shè)邊界區(qū)域如圖1所示，在笛卡兒坐標(biāo)系下，令z軸垂直向下為正方向，計算區(qū)域取水平地面為上邊界zmin，取地下離異常體一定遠處為下邊界zmax.

圖1 邊界條件示意圖Fig.1 Schematic diagram of boundary conditions

對于上邊界zmin，空氣與地下介質(zhì)的分界面滿足電位法向為零，由此，上邊界zmin處電位滿足

(15)

下邊界zmax離異常體一定距離，ja在下邊界處為零，由式(14)可得

(16)

該方程的通解為

(17)

因二次電位在下邊界只存在下行波e-kx，則有

(18)

對上式求導(dǎo)，可得

(19)

聯(lián)立式(15)和式(19)獲得控制方程式(14)的邊界條件：

(20)

1.3 邊值問題

綜上所述，各向異性介質(zhì)的空間-波數(shù)混合域異常電位滿足的邊值問題為

(21)

式(21)將各向異性介質(zhì)空間域異常電位滿足的三維偏微分方程轉(zhuǎn)化成不同波數(shù)滿足的一維常微分方程，由此大大減少了計算量以及對存儲的需求.因不同波數(shù)之間的常微分方程相互獨立，故分別使用有限單元法求解，并利用追趕法求解定帶寬線性方程組.

1.4 常微分方程求解

該方法保留垂直方向在空間域，這樣可以根據(jù)地下異常區(qū)域電流密度變化快慢靈活剖分，可準(zhǔn)確模擬地下異常區(qū)域.對于異常電位在空間-波數(shù)混合域中滿足的常微分方程式(21)，本文采用基于二次插值的一維有限單元法進行求解.基于變分原理(徐世浙,1994)，邊值問題式(21)等價變分問題為

在圖1中沿z向進行單元剖分，每個單元的異常電位和電流密度采用節(jié)點二次插值函數(shù)表示.通過離散，對(22)式逐項單元分析、總體合成，得到擴展矩陣或列陣：

(23)

(24)

(25)

式中，Ka是五對角矩陣，Pa是源項.

對于該五對角線性方程組，采用追趕法求解，通過空間波數(shù)混合域電場與電位的關(guān)系式(28)易得空間-波數(shù)混合域異常電場，然后使用標(biāo)準(zhǔn)FFT法(Tontini et al.，2009)對空間-波數(shù)混合域的電場進行反傅里葉變換，得到空間域異常電場.再將背景電場加上異常電場得到總電場，并采用壓縮算子迭代計算獲得空間域中修改的總電場.

1.5 壓縮算子

只有當(dāng)異常體電導(dǎo)率與背景介質(zhì)電導(dǎo)率之間差異較小時，求解式(25)才能得到較為準(zhǔn)確的解(Berdichevsky and Zhdanov, 1984)，通常地下介質(zhì)電導(dǎo)率差異較大(≥101)，此時求解式(25)得到的數(shù)值解與正確解差異很大，為此，Zhdanov和Fang(1996)和Gao(2005)在求基于積分方程的解時提出了一種迭代算子，鑒于三維直流電法積分解和微分解的等效性，將該迭代算子引入基于微分解的三維直流電法數(shù)值模擬當(dāng)中，該迭代算子的表達式為

(26)

式中，α、β均為張量，E(n-1)表示第n-1次求解式得到的總電場，E(n)′表示第n次迭代計算得到的總電場，E(n)表示采用壓縮算子修改后得到的第n次迭代的總電場.算法的步驟為：

(3)將背景電場Eb加上異常電場Ea(1)′得到總電場E(1)′；

(4)利用壓縮算子式(26)得到修正后的總電場E(1)；

(5)判斷修改前后總電場的相對殘差abs(E(1)-E(0))/abs(E(0))是否滿足期望的數(shù)值精度ε0，即abs(E(1)-E(0))/abs(E(0))<ε0.若滿足精度要求，輸出結(jié)果，否則取E(0)=E(1)，重復(fù)步驟(2)—(5)，直至滿足收斂條件.

1.6 異常電場的計算

因使用壓縮算子迭代時需要求出電場，由式(11)可知各向異性介質(zhì)空間-波數(shù)混合域的異常電場分量與異常電位的關(guān)系式為

(27)

1.7 零波數(shù)處理

在空間波數(shù)域中求解得電場后，采用標(biāo)準(zhǔn)FFT法進行反傅里葉變換時，需對波數(shù)k(kx,ky)為零的情況進行處理，否則會影響計算精度，造成較大誤差.當(dāng)波數(shù)kx,ky均為零時，異常電場的表達式為

(28)

2 算例及分析

算例均采用Fortran語言編寫的各向異性介質(zhì)空間-波數(shù)混合域三維直流電阻率數(shù)值模擬程序測試，測試平臺配置為Inter(R) Core(TM) i3-4150 CPU，@3.50 GHz(4核)，16 GB RAM.

2.1 算法正確性驗證

為了驗證本文算法的準(zhǔn)確性和計算精度，對各向同性半空間存在各向異性棱柱異常體模型1(如圖2所示)進行數(shù)值計算，計算結(jié)果與自適應(yīng)有限單元法(Ren et al.，2018c)對比.該模型的計算區(qū)域為200 m×200 m×81 m，背景電阻率為100 Ωm，異常體相關(guān)參數(shù)如下：頂部埋深為10 m，大小為12 m×12 m×10 m，三個電阻率主分量為ρ1/ρ2/ρ3=10/5/2 Ωm，對應(yīng)歐拉角為α/β/γ=30°/45°/60°.異常體中心位置在地面上的投影位于原點坐標(biāo).采用均勻網(wǎng)格進行剖分，網(wǎng)格大小為1 m.采用二極裝置測量，供電電極位于原點坐標(biāo)，電流大小為1 A.在使用本文算法計算時，使用標(biāo)準(zhǔn)FFT法進行傅里葉變換，為了減小截斷效應(yīng)帶來的誤差影響，分別沿x和y軸方向向外擴邊不同距離(0 m，20 m，40 m，80 m)進行計算.圖3表示本文算法和自適應(yīng)有限單元法在地面上分別沿x軸方向測量得到的視電阻率曲線圖及相對誤差圖，從圖中可看出本文算法與自適應(yīng)有限單元法的相對誤差小于1%，驗證了本文算法的準(zhǔn)確性.同時也說明了本文算法的精度隨著擴邊距離的增加而提高.

圖2 模型1示意圖Fig.2 Schematic diagram of model 1

2.2 迭代算法收斂性驗證

為了分析異常體大小，異常體埋深以及異常體電阻率與算法收斂性的關(guān)系，利用圖2模型進行驗證，計算區(qū)域仍然為200 m×200 m×81 m.驗證異常體大小與算法收斂性關(guān)系時，只改變異常體的大小，其他參數(shù)與圖2模型的一致，對比結(jié)果如表1所示.驗證異常體埋深與算法收斂性的關(guān)系時，只改變異常體的埋深，其他參數(shù)與圖2模型的一致，對比結(jié)果如表2所示.

表1 異常體大小與算法迭代次數(shù)的關(guān)系Table 1 Relationship between the iteration times of algorithm and the size of abnormal body

表2 異常體埋深與算法迭代次數(shù)的關(guān)系Table 2 Relationship between the iteration times of algorithm and the depth of abnormal body

從表1和表2可看出，在滿足精度要求的情況下，改變異常體的大小或者埋深，算法收斂時迭代次數(shù)不變，即異常體大小或埋深不影響算法的收斂性.

圖3 本文算法和自適應(yīng)有限單元法沿x軸方向測量得到的視電阻率圖(a)以及相對誤差圖(b)(Expanding 0 m表示不擴邊)Fig.3 Apparent resistivity (a) along the x-axis measured by the proposed algorithm and the adaptive finite element method and their relative errors (b) (‘Expanding 0m’ means calculation without expanding edge).

驗證異常體電阻率值與算法收斂性的關(guān)系時，由于異常體為各向異性，無法像各向同性條件下直接給出異常體電阻率為常數(shù)，即各向異性異常體的電阻率為張量，為了便于分析，令異常體三個方向的歐拉角等于0，只改變異常體沿著x軸方向的電阻率主分量的值，保持另外兩個方向的電阻率主分量不變，背景電阻率為100 Ωm，其他參數(shù)與圖2模型的一致，結(jié)果如表3所示.從中可看出，異常電阻率與背景電阻率的差異決定了算法的收斂速度，異常電阻率與背景電阻率之間對比度越大，收斂速度越慢.另外，從表3可看出高阻異常比低阻異常的收斂速度更快.改變y軸方向電阻率主分量或z軸方向電阻率主分量也得到類似的結(jié)論，故不再贅述.

表3 異常體電阻率差異與算法迭代次數(shù)的關(guān)系Table 3 Relationship between the iteration times of algorithm and difference of abnormal body resistivity

2.3 算法效率分析

目前基于各向異性介質(zhì)的直流電阻率法三維數(shù)值模擬算法的文獻中，關(guān)于算法計算效率的闡述較少.為了研究本文算法在不同剖分網(wǎng)格下計算時間的變化規(guī)律，設(shè)計均勻半空間中存在兩個異常體的模型來計算和驗證.模型2如圖4所示，計算區(qū)域為100 m×100 m×61 m，均勻半空間的電導(dǎo)率ρ0=100 Ωm，兩個異常體大小為12 m×12 m×10 m，頂面至地面距離均為10 m，兩個異常體之間距離為20 m，其中異常體A的電阻率為各向異性，電阻率主分量和歐拉角分別為ρ1/ρ2/ρ3=10/1000/1 Ωm，α/β/γ=0°/0°/0°，異常體B的電阻率為各向同性，ρB=10 Ωm.為了分析算法在不同剖分節(jié)點總數(shù)下的計算效率，令節(jié)點總數(shù)分別為101×101×81，201×201×81，301×301×81，401×401×81，501×501×81進行計算，不同節(jié)點總數(shù)下計算時間和占用內(nèi)存如圖5所示.從圖中可看出隨著剖分節(jié)點總數(shù)的增加，計算所需內(nèi)存和時間呈近線性增長，其中節(jié)點總數(shù)為101×101×81(即826281)時計算所需內(nèi)存為0.44 GB，計算時間僅為9.36 s，說明本文算法有較高的計算效率；節(jié)點總數(shù)為501×501×81(即20331081，節(jié)點數(shù)超過107)時計算所需內(nèi)存為10.37 GB，計算時間為250.20 s，說明本文算法在微型計算機下即可計算剖分節(jié)點總數(shù)超過千萬級的各向異性介質(zhì)模型，并能較快地算出結(jié)果.

圖4 模型2示意圖Fig.4 Schematic diagram of model 2

圖5 本文算法在不同剖分節(jié)點下占用內(nèi)存和計算時間Fig.5 The memory and time of the algorithm in different nodes

圖6 模型3示意圖Fig.6 Schematic diagram of model 3

2.4 模擬各向異性介質(zhì)下的響應(yīng)特征

通常，在野外為了觀測到地下結(jié)構(gòu)的各向異性特性，選擇在不同的位置發(fā)射場源，并在同一個位置觀測.為了模擬野外觀測系統(tǒng)，采用半空間存在一電阻率各向異性棱柱體模型進行計算和分析，模型3如圖6所示，背景為均勻各向同性半空間，電阻率ρ0=100 Ωm.異常體大小為100 m×100 m×100 m，異常體中心點在地面上的投影為坐標(biāo)原點，異常體頂面至地面距離為100 m，異常體電阻率主分量和歐拉角分別為ρ1/ρ2/ρ3=1/10/1000 Ωm，α/β/γ=30°/45°/60°.采用二極裝置測量，分別在點(200,0,0)和點(0,200,0)位置供電，供電電流大小為10 A，在原點(0,0,0)測量.通過測量可知，供電電極在(200,0,0)處時，原點(0,0,0)觀測到的視電阻率值為98.11 Ωm；供電電極在(0,200,0)處時，原點(0,0,0)觀測到的視電阻率值為95.10 Ωm.由此說明了該觀測系統(tǒng)能反應(yīng)出地下結(jié)構(gòu)為各向異性的特性.

圖7 模型6示意圖圖中帶箭頭實線表示收發(fā)距R=50 m，分別在相同收發(fā)距測量視電阻率.Fig.7 Schematic diagram of model 6The solid line with arrow indicates the receiving and transmitting distance R=50 m, the apparent resistivity is measured at the same receiving and transmitting distance around the source point.

為了進一步模擬和分析地下介質(zhì)各向異性特性，仍然以各向同性半空間中存在一各向異性棱柱體模型(如圖7所示)進行計算.計算區(qū)域為2000 m×2000 m×400 m，背景電阻率為ρ0=100 Ωm，異常體大小為100 m×100 m×100 m，異常體頂面至地面距離為100 m，異常體中心在地面上的投影為坐標(biāo)原點，在該處放置供電電極，向地下供電電流大小為10 A，保持供電電極的位置不變，在距離供電電極為50 m的環(huán)形圓位置測量.下面對異常體電阻率取不同值的情況進行計算和分析：

(1)異常體為單軸各向異性時，令它的電阻率主分量分別為ρ1/ρ2/ρ3=1000/10/10 Ωm，ρ1/ρ2/ρ3=10/1000/10 Ωm，ρ1/ρ2/ρ3=10/10/1000 Ωm，三種情況進行計算，此時這三種情況下歐拉角均等于零.

(2)異常體為三軸各向異性時，令它的電阻率主分量為ρ1/ρ2/ρ3=10/1000/1 Ωm，其旋轉(zhuǎn)的歐拉角分別為(a)α=0°/45°/90°/135°,β/γ=0°/0°; (b)β=0°/45°/90°/135°,α/γ=0°/0°; (c)γ=0°/30°/60°/90°,α/β=0°/0°，對這三種情況進行計算.

圖8—9表示環(huán)形測量收發(fā)距R=50 m得到的視電阻率極向圖.極化圖中各點到中心點(供電電極)的徑向長度表示視電阻率的大小，矢徑方向表示測點的方位.為了更直觀看出旋轉(zhuǎn)不同歐拉角后得到的視電阻率值之間的差異，中心點根據(jù)實際情況取不同的起始值.

圖8 異常體為各向同性或為單軸各向異性時，在收發(fā)距R=50 m環(huán)形測量得到的視電阻率極化圖其中ρ1/ρ2/ρ3=10/10/10 Ωm表示異常體電導(dǎo)率為各向同性.Fig.8 Apparent resistivity is obtained by ring measurement at R=50 m in the case of isotropic anomalous or uniaxial anisotropic anomalous body. ρ1/ρ2/ρ3=10/10/10 Ωm means that the conductivity of the abnormal body is isotropic.

圖8表示情況(1)下異常體為單軸各向異性，改變不同方向上電阻率主分量得到的極化圖，從圖中可看出，異常體為各向同性和各向異性時的視電阻率值不同，所以在實際應(yīng)用中不可忽略地下介質(zhì)的各向異性特性.另外，不同方位的電阻率單軸各向異性在同一測點計算得到的視電阻率值也不同.特別地，在此情況下z軸方向各向異性在地面上測量得到的視電阻率值大于沿著x軸或y軸方向各向異性測量得到的視電阻率值.

圖9表示情況(2)下異常體為三軸各向異性時，異常體分別沿x軸，y軸及z軸方向旋轉(zhuǎn)不同歐拉角得到的視電阻率極向圖，其中環(huán)形測量收發(fā)距R=50 m.其中圖9a表示沿x軸方向歐拉角α分別取0°/45°/90°/135°，而其他兩個歐拉角為0°得到的視電阻率極化圖；圖9b表示沿y軸方向歐拉角β分別取0°/45°/90°/135°，而其他兩個歐拉角為0°得到的視電阻率極化圖；圖9c表示沿z軸方向歐拉角γ分別取0°/45°/90°/135°，而其他兩個歐拉角為0°得到的視電阻率極化圖.從圖9可看出，異常體沿著同一個方向旋轉(zhuǎn)不同的歐拉角，得到的視電阻率不同，沿x方向歐拉角分別為45°和135°得到的視電阻率曲線關(guān)于x軸對稱；沿著y方向歐拉角分別為45°和135°得到的視電阻率曲線關(guān)于y軸對稱；而沿z方向歐拉角分別為0°和90°得到的視電阻率曲線關(guān)于45°和135°的軸線對稱；沿x軸旋轉(zhuǎn)得到的視電阻率值大于沿著y軸和z軸旋轉(zhuǎn)得到的視電阻率數(shù)值；另外，無論異常體沿著哪個軸方向旋轉(zhuǎn)，相同收發(fā)距的環(huán)形上測量得到的視電阻率曲線均為橢圓形.

圖9 異常體為三軸各向異性時，在收發(fā)距R=50 m環(huán)形測量得到的視電阻率極化圖(a) α=0°/45°/90°/135°, β=0°, γ=0°; (b) α=0°,β=0°/45°/90°/135°,γ=0°; (c) α=0°, β=0°,γ=0°/45°/90°/135°.Fig.9 Apparent resistivity is obtained by ring measurement at R=50 m when the resistivity of the anomalous body is triaxial anisotropic

3 結(jié)論

本文將各向異性介質(zhì)分成各向同性背景介質(zhì)和各向異性異常介質(zhì)，并提出了一種適用于計算各向異性介質(zhì)的空間波數(shù)混合域直流電法三維數(shù)值模擬方法.該方法的策略是沿水平方向進行二維傅里葉變換，將空間域異常電位滿足的三維偏微分方程轉(zhuǎn)化成不同波數(shù)滿足的一維常微分方程，這樣將一個大規(guī)模三維數(shù)值模擬問題分解為多個一維數(shù)值模擬問題，具有高度并行性.該方法充分利用不同波數(shù)問題之間高度并行性，采用一維有限單元法求解不同波數(shù)滿足的常微分方程，并引入壓縮算子進行迭代，實現(xiàn)了直流電法三維數(shù)值模擬.

通過與自適應(yīng)有限單元法對比驗證精度，說明將各向異性介質(zhì)分成各向同性背景介質(zhì)和各向異性異常介質(zhì)是可行的；分析半空間下存在兩個異常體模型的計算效率，說明了本文算法在微型計算機下能較快算出剖分節(jié)點總數(shù)超過千萬的模型的結(jié)果；異常體為各向異性情況下，按不同方向旋轉(zhuǎn)歐拉角的響應(yīng)特征結(jié)果表明，異常體為各向異性比異常體為各向同性的響應(yīng)特征更復(fù)雜；各向異性異常體沿著不同的方向旋轉(zhuǎn)得到的視電阻率不同，即使在同一方向上旋轉(zhuǎn)不同角度，得到的視電阻率差異也很大.特別說明的是，本文算法具有高度并行性以及可用于帶地形的復(fù)雜模型的計算，后續(xù)將研究它基于OpenMP或GPU平臺的并行算法以及帶地形的復(fù)雜模型的計算.

附錄A

(A1)

(A2)

(A3)

(A4)

(A5)

xxe、yye、zze對應(yīng)為單元內(nèi)三個方向的電流密度.由于z=zmax表示最后一個節(jié)點，因此，式(22)中的最后一項可以擴展為

(A6)

式中，

(A7)