基于矢量有限元的大地電磁快速三維正演研究

顧觀文，武曄，石硯斌

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

0 引言

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

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

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

1 大地電磁三維正演算法

1.1 大地電磁三維正演的邊值問(wèn)題

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

(1a)

(1b)

(1c)

(1d)

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

(2)

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

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

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

將式(2)寫(xiě)成

(3)

取第一類邊界條件：

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

(4)

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

1.2 矢量有限元分析

用有限元法求解上述區(qū)域(圖1)的電磁場(chǎng)問(wèn)題需要對(duì)研究區(qū)域離散化，即對(duì)研究區(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)，每個(gè)六面體單元的內(nèi)部電場(chǎng)分量用六面體的12條棱邊的場(chǎng)值分量(如圖2b所示)通過(guò)插值求取，公式為

(5)

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

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

(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方向的長(zhǎng)度；(ξi,ηi,ζi)的取值與棱邊編號(hào)有關(guān)。

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

(10)

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

(11)

(12)

利用矢量恒等式：

(13)

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

(14)

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

(15)

對(duì)于第e個(gè)單元，式(11)可化為

(16)

KeEe=Se,

(17)

式中：Se表示場(chǎng)源項(xiàng)；Ee表示棱邊上的電場(chǎng)；Ke為單元?jiǎng)偠染仃?，是一個(gè)12×12階的復(fù)數(shù)矩陣，可按下式解析計(jì)算[32]得出：

(18)

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

K·E=s,

(19)

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

1.3 線性方程組求解

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

1.3.1 PARDISO求解過(guò)程

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

Ax=b，

(20)

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

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

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

1.3.2 壓縮稀疏矩陣存儲(chǔ)

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

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

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

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

以式(21)為例說(shuō)明CSR存儲(chǔ)格式：

(21)

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

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

1.4 視電阻率及阻抗相位的計(jì)算

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

(22)

張量阻抗的每一個(gè)分量可表示為

(23a)

(23b)

(23c)

(23d)

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

(24)

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

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

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

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

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

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

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

觀測(cè)頻率范圍為10 000～0.000 1 Hz，共取21個(gè)頻點(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è)為一個(gè)較大的常數(shù)，一般在106～1010Ω·m之間，本文涉及的三維模型其空氣層電阻率均取為1010Ω·m，表1、圖3為數(shù)值模擬結(jié)果。從表1可以看出，視電阻率誤差最高為0.743 7%(頻點(diǎn)為50 Hz)，其余頻點(diǎn)計(jì)算誤差均未超過(guò)0.6%，尤其在0.1 Hz以后誤差均在0.1%以下。相位計(jì)算誤差最高為0.233 56%(頻點(diǎn)為100 Hz)，其余頻點(diǎn)計(jì)算誤差均在0.2%以下。

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

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

圖3 均勻半空間模型三維正演視電阻率(a)和相位(b)與解析解對(duì)比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%(觀測(cè)頻率50 Hz)，其他頻點(diǎn)誤差均在2%以下，表明模擬精度符合要求。

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

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

圖4 H型層狀模型三維正演視電阻率(a)和相位(b)數(shù)值解與解析解對(duì)比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%(觀測(cè)頻率50 Hz)，其他頻點(diǎn)誤差均在3.5%以下，表明模擬精度符合要求。

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

2.1.2 COMMEMI 3D-2模型

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

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

圖5 K型層狀模型三維正演視電阻率(a)和相位(b)數(shù)值解與解析解對(duì)比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層為空氣層)個(gè)網(wǎng)格，采用矢量有限元法對(duì)該模型進(jìn)行三維正演模擬，并將矢量有限元模擬結(jié)果與Wannamaker等的積分方程法(integral equation，IE)三維模擬結(jié)果[1]進(jìn)行比較。圖7為觀測(cè)頻率f=0.001 Hz的模擬結(jié)果對(duì)比，圖中實(shí)心黑色圓圈是本文的無(wú)需散度校正基于PARDISO直接求解的矢量有限元法(VFE-PARDISO)計(jì)算結(jié)果，方框是積分方程法(IE)計(jì)算結(jié)果，可以看出二種方法的計(jì)算結(jié)果基本一致。

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

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

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

對(duì)上述模型(圖9)進(jìn)行三維正演模擬，計(jì)算頻率f=2 Hz時(shí)TE極化模式和TM極化模式的視電阻率和相位曲線。圖10中黑色虛線和實(shí)線分別是二維有限元(2D FE)法的TE、TM極化模式的計(jì)算結(jié)果，圓圈和方框是本文的矢量有限元法(VFE-PARDISO)的xy模式、yx模式的計(jì)算結(jié)果。從圖中可以看出，本研究算法(VFE-PARDISO)與二維有限元的視電阻率、相位的計(jì)算結(jié)果一致，從而說(shuō)明了本文研究的三維正演算法對(duì)起伏地形地電模型的計(jì)算結(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 本文矢量有限元正演算法的計(jì)算結(jié)果與IE方法的計(jì)算結(jié)果對(duì)比Fig.7 Comparison between the calculation results of vector finite element forward algorithm and IE method

a—三維地形示意； b—地形網(wǎng)格剖分和MT測(cè)點(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)計(jì)算的二維地形影響與二維有限元結(jié)果對(duì)比Fig.10 Comparision between modeling results of 3DVFEM(PARDISO) and 2DFEM for 2D ridge

3 無(wú)需散度校正的直接解法與帶散度校正的迭代解法計(jì)算對(duì)比

為了對(duì)比無(wú)需散度校正的PARDISO直接解法(VFE-PARDISO without divergence correction)和帶散度校正的BICG迭代解法(VFE-BICG with divergence correction)的計(jì)算精度和計(jì)算時(shí)間，分別采用這兩種求解方法對(duì)如圖9所示的二維山峰地形模型進(jìn)行三維正演模擬，并對(duì)比模擬結(jié)果。兩種求解方法的三維正演模擬均在曙光W560-G20工作站上完成，計(jì)算及程序編譯環(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)。整個(gè)模型區(qū)域沿x、y和z方向剖分為43×43×29(其中z方向山峰頂面以上7層為空氣層)個(gè)網(wǎng)格單元。

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

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

4 結(jié)論

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