顧觀文，武曄，石硯斌

(1.防災(zāi)科技學(xué)院 地球科學(xué)學(xué)院，河北 廊坊 065201; 2.河北省地震動力學(xué)重點(diǎn)實(shí)驗(yàn)室，河北 廊坊065201)

0 引言

大地電磁測深法(magnetotelluric sounding, MT)具有施工方便、勘探效率高、成本低 (相對于地震勘探)、勘探深度大等優(yōu)點(diǎn)，目前已被廣泛應(yīng)用于資源勘查、能源勘探及深部構(gòu)造探測等方面。對于大地電磁數(shù)據(jù)的反演和解釋，由于受限于理論方法及計算設(shè)備的運(yùn)算能力，早期以一維或二維反演技術(shù)為主。20 世紀(jì)90年代以來，隨著計算設(shè)備運(yùn)算能力的提高及數(shù)值計算方法的進(jìn)步，大地電磁三維正反演技術(shù)開始逐步發(fā)展，不同的三維正演方法(積分方程、有限差分、有限元等)及其計算技術(shù)取得了巨大進(jìn)展[1-15]。目前在實(shí)際中得到應(yīng)用的三維反演技術(shù)主要是基于有限差分法三維正演的反演方法，特別是國內(nèi)實(shí)測大地電磁資料的三維解釋基本上都采用基于有限差分法的三維反演技術(shù)[16-19]。不同于有限差分法，有限單元法在模擬起伏地形以及復(fù)雜地質(zhì)體的電磁響應(yīng)方面具有明顯優(yōu)勢，特別是近些年發(fā)展迅速的矢量有限元法，由于其能有效地解決傳統(tǒng)節(jié)點(diǎn)有限元法存在的偽解問題，目前已成為復(fù)雜地形和復(fù)雜地質(zhì)體三維電磁響應(yīng)模擬的主要方法。但有限單元法也存在一些不足，運(yùn)算量大、計算時間長是導(dǎo)致基于有限元法的大地電磁三維反演技術(shù)實(shí)用化進(jìn)程相對滯后(相對于基于有限差分法的三維反演技術(shù))的主要因素。為此，開展基于有限元的MT快速三維正演算法研究，提高三維正演計算效率，對于大地電磁三維反演技術(shù)的實(shí)用性具有重要的意義。

MT三維正演模擬最終需要求解復(fù)數(shù)大型線性方程組，如何快速、準(zhǔn)確地求解此線性方程成為大地電磁三維正演實(shí)現(xiàn)中的關(guān)鍵工作，此線性方程的求解時間約占整個正演計算時間的80%[20]。目前求解大型線性方程組主要有兩種方法，一類為Krylov子空間迭代法，另一類為直接解法。迭代算法從一個初始解出發(fā)，通過逐步修正來逼近真實(shí)解，在每一步的計算過程中只需要計算矩陣向量乘積或矩陣的轉(zhuǎn)置與向量的乘積，因此其內(nèi)存需求較小。但是，迭代算法可能存在不收斂的問題：一方面，方程組階數(shù)的增加導(dǎo)致矩陣條件數(shù)增大；另一方面，求解低頻電磁場問題時，方程的定解條件變?nèi)?，也加大了矩陣的病態(tài)性，使迭代算法收斂變慢，導(dǎo)致過多的迭代求解步驟，從而增加了計算時間。采用迭代法求解線性方程組時，需要引入散度校正技術(shù)以改善在低頻范圍(尤其在0.01 Hz及以下)求解的收斂性[5,11,21-22]。直接解法對計算機(jī)的內(nèi)存容量要求較高，但直接法求解穩(wěn)定，對于給定的方程組，直接解法總能獲得其精確解，尤其是對于大地電磁的低頻正演模擬問題，無需做散度校正也可獲得高精度的數(shù)值解。值得關(guān)注的是，近些年直接求解器獲得了非常大的突破，產(chǎn)生了有效的直接法求解器，如基于LU矩陣分解法和OpenMP的PARDISO(parallel direct solver)求解器[23]，基于波前法和MPI的MUMPS(multifrontal massive parallel solver)求解器[24-25]?；谥苯咏夥ㄔ诮鉀Q大規(guī)模電磁模擬問題方面取得了明顯效果[13,26-28]。直接求解法得到了越來越多的重視，有逐步取代迭代求解器的發(fā)展趨勢[29]?；赑ARDISO直接求解器的直接求解法具有計算精度高和無需散度校正等優(yōu)勢，同時PARDISO求解器支持復(fù)數(shù)雙精度數(shù)據(jù)類型、64位整數(shù)索引、并可搭配共享內(nèi)存環(huán)境下的OpenMP并行方式，可以方便地應(yīng)用在電磁有限元分析后的大型稀疏矩陣的求解中。

根據(jù)大地電磁場滿足的麥克斯韋方程組推導(dǎo)了大地電磁三維正演的邊值問題，利用矢量有限元法基于六面體網(wǎng)格單元對計算區(qū)域進(jìn)行離散，采用無需散度校正的直接求解方法(PARDISO)求解矢量有限元法對應(yīng)的大型線性方程組。在對三維地電模型進(jìn)行數(shù)值模擬的過程中，通過數(shù)值解與解析解對比以及本文模擬結(jié)果與國際公認(rèn)檢驗(yàn)?zāi)Ｐ偷挠嬎憬Y(jié)果對比，驗(yàn)證三維矢量有限元方法及程序的正確性。通過無需散度校正直接求解法與帶散度校正的迭代求解法對比，證明本文快速三維正演算法的有效性。

1 大地電磁三維正演算法

1.1 大地電磁三維正演的邊值問題

正演方法基于六面體網(wǎng)格剖分的矢量有限元法[11-12]，在大地電磁研究的頻率范圍內(nèi)，忽略位移電流的作用。取時諧場為e-iωt，麥克斯韋方程組的微分形式表示如下：

(1a)

(1b)

(1c)

(1d)

式中：E為電場強(qiáng)度矢量，H為磁場強(qiáng)度矢量，σ是介質(zhì)的電導(dǎo)率，μ是介質(zhì)的磁導(dǎo)率，ε是介質(zhì)的介電常數(shù)，q為自由電荷密度。對式(1a)兩邊取旋度，再將式(1b)代入到式(1a)，則將電磁場滿足的一階微分方程變?yōu)槎A微分方程，可得電場E滿足下面方程：

(2)

求解偏微分方程組(2)得到電場分量Ex、Ey、Ez后，根據(jù)式(1a)求得磁場分量Hx、Hy、Hz。

由于空氣中的電場和磁場受地形和不均勻體影響而不均勻，因此在數(shù)值模擬中必須包括空氣層和地下介質(zhì)。如圖1所示，模擬區(qū)域設(shè)為Ω，由空氣層和地下介質(zhì)兩部分組成，Ω=Ω1∪Ω2。在模擬區(qū)域中，空氣層的電導(dǎo)率一般在10-6～10-10S·m-1之間，這時空氣和地下介質(zhì)的接觸面變成內(nèi)部邊界。

圖1 帶地形三維MT數(shù)值模擬區(qū)域剖面示意[11]Fig.1 Section diagram of numerical modeling domain for 3D MT with topography[11]

將式(2)寫成

(3)

取第一類邊界條件：

E(x,y,z)|?Ω=g(x,y,z)|?Ω,

(4)

式(4)中g(shù)是邊界上的矢量電場，可以用一維或者二維MT計算值[4,30]。這樣，式(3)和式(4)構(gòu)成了大地電磁三維正演的邊值問題。

1.2 矢量有限元分析

用有限元法求解上述區(qū)域(圖1)的電磁場問題需要對研究區(qū)域離散化，即對研究區(qū)域進(jìn)行六面體網(wǎng)格剖分。如圖2a所示，沿x、y和z軸方向分別剖分成Nx、Ny和Nz段，網(wǎng)格間距分別為Δx(i)(i=1,…,Nx)、Δy(j)(j=1,…,Ny)和Δz(k)(k=1,…,Nz)。可以推導(dǎo)，每個六面體單元的內(nèi)部電場分量用六面體的12條棱邊的場值分量(如圖2b所示)通過插值求取，公式為

(5)

a—區(qū)域剖分示意; b—電場分量位置示意a—domain subdivision; b—location of electric field components圖2 矢量有限元法的區(qū)域剖分示意Fig.2 Domain subdivision of the vector finite element method

式(5)寫成矢量形式如下：

(6)

(7)

(8)

(9)

式中：坐標(biāo)轉(zhuǎn)換函數(shù)為ξ=(x-xc)/a,η=(y-yc)/b,ζ=(z-zc)/c；(xc,yc,zc)是六面體的中心坐標(biāo)，2a、2b、2c分別是六面體x、y、z方向的長度；(ξi,ηi,ζi)的取值與棱邊編號有關(guān)。

由式(3)定義矢量余函數(shù)：

(10)

把矢量基函數(shù)作為權(quán)函數(shù)，采用迦遼金方法[31-32]使整個域內(nèi)的積分矢量余函數(shù)為最小，即：

(11)

(12)

利用矢量恒等式：

(13)

則式(12)右邊第一項(xiàng)分為兩項(xiàng)：

(14)

根據(jù)高斯散度定理，可知：

(15)

對于第e個單元，式(11)可化為

(16)

KeEe=Se,

(17)

式中：Se表示場源項(xiàng)；Ee表示棱邊上的電場；Ke為單元剛度矩陣，是一個12×12階的復(fù)數(shù)矩陣，可按下式解析計算[32]得出：

(18)

將每個單元電場滿足的線性方程進(jìn)行組合，可以得到整個計算域上電場滿足的線性方程組：

K·E=s,

(19)

式中：K是系統(tǒng)剛度矩陣；E是整個計算域的網(wǎng)格單元棱邊上的電場值向量；s是源向量，由計算域的上、下、左、右的邊界場值與邊界上的單元剛度矩陣計算得到。

1.3 線性方程組求解

方程(19)為復(fù)數(shù)大型線性方程組，PARDISO是針對大規(guī)模稀疏線性方程組開發(fā)的高效并行直接求解器，采用PARDISO直接求解器求解方程(19)。該求解器采用BLAS軟件包，實(shí)現(xiàn)算法中線性代數(shù)操作的并行計算。大量數(shù)值測試表明，PARDISO是目前最快的線性稀疏矩陣求解方法之一[33]。Intel公司獲得授權(quán)后，Intel?Math Kernel Library(Intel MKL)提供了PARDISO的優(yōu)化版本，計算效率高，穩(wěn)定性好。

1.3.1 PARDISO求解過程

PARDISO求解器基于LU分解，融合了向左和向右算法的優(yōu)點(diǎn)，采用超節(jié)點(diǎn)消去樹進(jìn)行填充元優(yōu)化，降低算法的時空消耗。對于復(fù)線性方程組有

Ax=b，

(20)

式中：A矩陣是一個稀疏矩陣，矩陣中的元索值多數(shù)為零，只有與當(dāng)前節(jié)點(diǎn)直接相連的節(jié)點(diǎn)處為非零元素。利用PARDISO求解方程(20)的具體步驟如下：

③ 方程求解與迭代：利用LU分解結(jié)果，求解方程。如對結(jié)果精度有進(jìn)一步要求，使用迭代法進(jìn)一步提高精度。

④ 迭代結(jié)束，釋放計算過程所占內(nèi)存。

1.3.2 壓縮稀疏矩陣存儲

PARDISO求解器采用目前廣泛流行的壓縮稀疏行格式(compressed sparse row, CSR)對稀疏矩陣進(jìn)行壓縮存儲，大大降低了矩陣元素的訪問時間和空間存儲成本。該方法以行為單位存儲每個非零數(shù)據(jù)。對于一個對稱、稀疏矩陣A，PARDISO對矩陣的存儲包括三個數(shù)組：

① 雙精度數(shù)組a—矩陣A上三角部分的非零元素。A的非零元素通過下面的ja與ia映射到a數(shù)組中。

② 整型數(shù)組ja—數(shù)組a中每個元素在矩陣A中的列號。

③ 整型數(shù)組ia—矩陣A中上三角部分每行第一個非零元素在數(shù)組a中的位置，其最后一個元素為上三角矩陣的非零元素總數(shù)加一。

以式(21)為例說明CSR存儲格式：

(21)

大地電磁三維正演問題(式(3)、式(4))經(jīng)有限元離散后，得到的剛度矩陣(式(19) 中的K)為稀疏、對稱矩陣，可采用對稱CSR格式存儲。

根據(jù)大地電磁場線性方程組特有的稀疏性，開發(fā)有針對性的數(shù)據(jù)結(jié)構(gòu)接口，將PARDISO快速求解器應(yīng)用到方程組(19)的求解，得到計算域上網(wǎng)格單元棱邊上的電場值，然后根據(jù)麥克斯韋方程組(式(1a))微分求取磁場。

1.4 視電阻率及阻抗相位的計算

根據(jù)Newman等[34]、譚捍東等[35]的研究，假設(shè)兩種線性無關(guān)的場源激發(fā)的表面電場和磁場分別為Ex1、Ey1、Hx1、Hy1、Ex2、Ey2、Hx2、Hy2，由此可以計算三維MT的張量阻抗：

(22)

張量阻抗的每一個分量可表示為

(23a)

(23b)

(23c)

(23d)

式中：下標(biāo)1、2表示極化模式，下標(biāo)x、y表示x、y分量，E表示電場，H表示磁場，Z表示阻抗；式23b和23c定義的響應(yīng)分別稱為XY和YX模式響應(yīng)，按照式24可求出三維介質(zhì)的視電阻率和相位：

(24)

式中:i=X,Y;j=X,Y。

2 三維正演結(jié)果的正確性驗(yàn)證

為了驗(yàn)證本文三維正演算法的正確性，分別對水平地形地電模型和起伏地形地電模型進(jìn)行正確性驗(yàn)證。對于水平地形條件下的驗(yàn)證，采用均勻半空間、典型層狀模型和國際電磁學(xué)術(shù)討論會COMMEMI研究小組Zhdanov等設(shè)計的COMMEMI 3D-2模型[36]，對于起伏地形條件下的驗(yàn)證，采用Wannamaker等設(shè)計的目前國際上普遍認(rèn)可的二維山峰模型[37]。

2.1 水平地形條件下的模型驗(yàn)證

2.1.1 均勻半空間模型和層狀模型

1) 模型網(wǎng)格剖分及觀測頻率

以下的均勻半空間模型、H型和K型層狀模型均采用相同網(wǎng)格剖分，并選取相同的觀測頻率。沿x、y和z方向剖分為41×41×35(其中z方向地面以上10層為空氣層)個網(wǎng)格。

觀測頻率范圍為10 000～0.000 1 Hz，共取21個頻點(diǎn)，分別為(單位：Hz)：10 000，8 000，5 000，2 000，1 000，500，200，100，50，10，5，2，1，0.5，0.1，0.05，0.01，0.005，0.001，0.000 5，0.000 1。

2) 均勻半空間模型

大地的電阻率設(shè)定為100 Ω·m，空氣通常認(rèn)為是絕緣體，但為了求解方程2，需要模型的電阻率為有限值，空氣層的電阻率設(shè)為一個較大的常數(shù)，一般在106～1010Ω·m之間，本文涉及的三維模型其空氣層電阻率均取為1010Ω·m，表1、圖3為數(shù)值模擬結(jié)果。從表1可以看出，視電阻率誤差最高為0.743 7%(頻點(diǎn)為50 Hz)，其余頻點(diǎn)計算誤差均未超過0.6%，尤其在0.1 Hz以后誤差均在0.1%以下。相位計算誤差最高為0.233 56%(頻點(diǎn)為100 Hz)，其余頻點(diǎn)計算誤差均在0.2%以下。

從圖3給出的視電阻率、相位的對比曲線可以看出，在10 000～0.000 1 Hz頻段范圍，三維矢量有限元正演(vector finite element with PARDISO solver(VFE-PARDISO)曲線與解析解曲線重合。

表1 均勻半空間模型矢量有限元解與解析解對比

圖3 均勻半空間模型三維正演視電阻率(a)和相位(b)與解析解對比Fig.3 Comparison of 3D forward apparent resistivity (a) and phase (b) of homogeneous half space model with analytical solution

3) H型層狀模型

H型模型參數(shù)設(shè)置如表2所示。表3、圖4為H型層狀模型三維數(shù)值模擬結(jié)果，可以看出，矢量有限元數(shù)值解與解析解的最大誤差為1.8%(觀測頻率50 Hz)，其他頻點(diǎn)誤差均在2%以下，表明模擬精度符合要求。

表2 H型層狀模型參數(shù)

表3 H型層狀模型矢量有限元解與解析解對比

圖4 H型層狀模型三維正演視電阻率(a)和相位(b)數(shù)值解與解析解對比Fig.4 Comparison of apparent resistivity (a) and phase (b) of 3D forward modeling of H-type layered model with analytical solutions

(4) K型層狀模型

K型模型參數(shù)設(shè)置如表4所示。表5、圖5為K型層狀模型三維數(shù)值模擬結(jié)果，可以看出，矢量有限元數(shù)值解與解析解的最大誤差為5.5%(觀測頻率50 Hz)，其他頻點(diǎn)誤差均在3.5%以下，表明模擬精度符合要求。

表4 K型層狀模型參數(shù)

2.1.2 COMMEMI 3D-2模型

COMMEMI 3D-2 模型是國際電磁學(xué)術(shù)討論會COMMEMI研究小組設(shè)計的一個三維地電模型[37]，該模型包含了多個電性差異巨大的分界面，比較復(fù)雜，已成為國內(nèi)外學(xué)者普遍認(rèn)可用于MT三維正演算法驗(yàn)證的三維模型[1,4,9,11-12,38]。圖6為COMMEMI 3D-2模型示意，背景是一個三層K型地電斷面，第一層電阻率10 Ω·m，厚度10 km；第二層電阻率100 Ω·m，厚度20 km；第三層電阻率為0.1 Ω·m。在剖面的第一層中，鑲嵌有電阻率分別為10 Ω·m的低阻棱柱體和100 Ω·m的高阻棱柱體，對稱地分布在x軸的兩側(cè)，棱柱體的規(guī)模為20 km×40 km×10 km，長軸為y方向，短軸為x方向。

表5 K型層狀模型矢量有限元解與解析解對比

圖5 K型層狀模型三維正演視電阻率(a)和相位(b)數(shù)值解與解析解對比Fig.5 Comparison of 3D forward apparent resistivity (a) and phase (b) of K-type layered model with analytical solution

圖6 COMMEMI3D-2 模型示意Fig.6 Schematic diagram of COMMEMI3D-2

將該模型沿x、y和z方向剖分為28×21×19(其中z方向地面以上8層為空氣層)個網(wǎng)格，采用矢量有限元法對該模型進(jìn)行三維正演模擬，并將矢量有限元模擬結(jié)果與Wannamaker等的積分方程法(integral equation，IE)三維模擬結(jié)果[1]進(jìn)行比較。圖7為觀測頻率f=0.001 Hz的模擬結(jié)果對比，圖中實(shí)心黑色圓圈是本文的無需散度校正基于PARDISO直接求解的矢量有限元法(VFE-PARDISO)計算結(jié)果，方框是積分方程法(IE)計算結(jié)果，可以看出二種方法的計算結(jié)果基本一致。

2.2 起伏地形條件下模型驗(yàn)證

對帶地形模型采用矩形六面體剖分，矩形網(wǎng)格x，y和z方向的間距可以是不均勻的，比如圖8所示的三維山峰地形剖分。若在地形起伏界面(空氣—地面)上加大網(wǎng)格剖分密度,使其剖分足夠精細(xì)時,大地電磁場數(shù)值模擬響應(yīng)將與理論場值相近,甚至趨于一致[39]。

為了驗(yàn)證正演程序?qū)T地形影響數(shù)值模擬效果，采用Wannamaker等設(shè)計的二維山峰模型[37]。模型背景電阻率為100 Ω·m，山峰地形如圖9所示。設(shè)二維山峰地形的走向?yàn)閤方向，山峰傾向?yàn)閥方向，z方向垂直向下。為了盡量避免對y方向的三維影響，將山峰地形沿x方向延伸3.43 km，整個模型區(qū)域(34 300 m×34 300 m×91 993 m)沿x、y和z方向剖分為43×43×29(其中z方向山峰頂面以上7層為空氣層)個網(wǎng)格單元。山峰地形采用9個縱向網(wǎng)格單元的劃分，網(wǎng)格間距為50 m。

對上述模型(圖9)進(jìn)行三維正演模擬，計算頻率f=2 Hz時TE極化模式和TM極化模式的視電阻率和相位曲線。圖10中黑色虛線和實(shí)線分別是二維有限元(2D FE)法的TE、TM極化模式的計算結(jié)果，圓圈和方框是本文的矢量有限元法(VFE-PARDISO)的xy模式、yx模式的計算結(jié)果。從圖中可以看出，本研究算法(VFE-PARDISO)與二維有限元的視電阻率、相位的計算結(jié)果一致，從而說明了本文研究的三維正演算法對起伏地形地電模型的計算結(jié)果準(zhǔn)確可靠。

a—Zxy模式正演視電阻率；b—Zyx模式正演視電阻率；c—Zxy模式正演阻抗相位； d—Zyx模式正演阻抗相位a—Zxy mode forward apparent resistivity；b—Zyx mode forward apparent resistivity；c—Zxy mode forward impedance phase； d—Zyxmode forward impedance phase圖7 本文矢量有限元正演算法的計算結(jié)果與IE方法的計算結(jié)果對比Fig.7 Comparison between the calculation results of vector finite element forward algorithm and IE method

a—三維地形示意； b—地形網(wǎng)格剖分和MT測點(diǎn)分布示意a—sketch of 3D topography； b—topography meshing and distribution of MT measurement sites圖8 三維地形及其網(wǎng)格剖分示意Fig.8 3D topography and grid

圖9 二維山峰地形示意Fig.9 Sketch of 2D ridge

a—正演視電阻率； b—正演阻抗相位a—forward apparent resistivity； b—forward impedance phase圖10 三維矢量有限元算法(PARDISO)計算的二維地形影響與二維有限元結(jié)果對比Fig.10 Comparision between modeling results of 3DVFEM(PARDISO) and 2DFEM for 2D ridge

3 無需散度校正的直接解法與帶散度校正的迭代解法計算對比

為了對比無需散度校正的PARDISO直接解法(VFE-PARDISO without divergence correction)和帶散度校正的BICG迭代解法(VFE-BICG with divergence correction)的計算精度和計算時間，分別采用這兩種求解方法對如圖9所示的二維山峰地形模型進(jìn)行三維正演模擬，并對比模擬結(jié)果。兩種求解方法的三維正演模擬均在曙光W560-G20工作站上完成，計算及程序編譯環(huán)境如下：CPU為Intel E5-2643 (3.4G)，內(nèi)存為64GB，操作系統(tǒng)為Windows 7(64位)；編譯環(huán)境為Microsoft Visual Studio 2012(已集成Intel Parallel Studio XE 2013)。整個模型區(qū)域沿x、y和z方向剖分為43×43×29(其中z方向山峰頂面以上7層為空氣層)個網(wǎng)格單元。

兩種算法的計算結(jié)果(觀測頻率為2 Hz)對比如圖11所示，圖11中黑色虛線和實(shí)線分別是二維有限元計算的TE和TM極化模式的計算結(jié)果(視電阻率和相位)，圓圈和方框是本文的無需散度校正直接求解的矢量有限元法(VFE-PARDISO)計算的xy模式(視電阻率、相位)和yx模式(視電阻率、相位)曲線圖，菱形和三角形是帶散度校正迭代求解的矢量有限元法(VFE-BICG)計算的XY模式(視電阻率、相位)和YX模式(視電阻率、相位)曲線圖。從圖11中可以看出，本文的算法(VFE-PARDISO)和迭代求解算法(VFE-BICG)計算的結(jié)果均與二維有限元的計算結(jié)果一致，但在平地與山峰拐點(diǎn)處，PARDISO直接求解的TM結(jié)果比BICG迭代求解的結(jié)果更接近于二維有限元計算結(jié)果。在計算時間方面，無需散度校正直接求解的矢量有限元法(VFE-PARDISO)正演一次上述模型僅耗時24 s，而帶散度校正迭代求解的矢量有限元法(VFE-BICG)耗時407 s；本文的直接解法(VFE-PARDISO)與迭代解法(VFE-BICG)的計算速度比達(dá)17倍?？梢姳疚牡臒o需散度校正直接求解法與帶散度校正的迭代求解法相比，在計算精度和計算時間方面均有優(yōu)勢，特別是在計算時間方面表現(xiàn)出明顯優(yōu)勢，適合應(yīng)用于中等計算規(guī)模的三維反演算法中。

a—正演視電阻率； b—正演阻抗相位a—forward apparent resistivity； b—forward impedance phase圖11 無需散度校正的直接解法與帶散度校正的迭代解法計算結(jié)果對比曲線Fig.11 Comparison of calculation results of VFE-PARDISO without divergence correction and VFE-BICG with divergence correction

4 結(jié)論

采用基于并行直接稀疏求解器PARDISO且無需散度校正的正演方案，并利用C++語言編制正演計算程序，實(shí)現(xiàn)MT三維快速正演。為了檢驗(yàn)本文的快速正演算法及計算程序的正確性，通過數(shù)值解與解析解對比以及本文模擬結(jié)果與國際公認(rèn)檢驗(yàn)?zāi)Ｐ蚚36-37]的計算結(jié)果對比，結(jié)果表明該算法及程序在水平地形和起伏地形條件下均滿足三維正演計算的精度要求。在中等規(guī)模計算條件下，通過本文的無需散度校正直接求解法與帶散度校正的迭代求解法對比，本文的無需散度校正直接求解法在計算精度和計算時間方面均有優(yōu)勢，特別是在計算時間方面表現(xiàn)出明顯優(yōu)勢，直接解法與迭代解法的計算速度比達(dá)17倍。本文實(shí)現(xiàn)的快速三維正演算法對于促進(jìn)MT三維反演技術(shù)的實(shí)用性具有現(xiàn)實(shí)意義。