時(shí)間域激發(fā)極化法三維有限元正演

駱紅梅

（福建省地質(zhì)測(cè)繪院，福建 福州 350011）

時(shí)間域激發(fā)極化法三維正演國(guó)內(nèi)外有不少的研究，大多都是在直流電正演的基礎(chǔ)上引入等效電阻率或COLECOLE 模型來做等效計(jì)算，然后采用有限元技術(shù)來實(shí)現(xiàn)激發(fā)極化場(chǎng)的正演[1]。但是三維有限元的計(jì)算速度和龐大的內(nèi)存需求一直是困擾激發(fā)極化法三維反演的主要問題[2]。國(guó)內(nèi)外近幾年針對(duì)以上問題發(fā)表了不少文章，主要是采用定半帶寬存儲(chǔ)稀疏矩陣，直接采用LDLT 分解法求解方程組；采用一維非零元素壓縮存儲(chǔ)模式，然后運(yùn)用預(yù)條件共軛梯度法（PCG）求解方程組。

1 正演方法及驗(yàn)證

1.1 正演方法

本文采用三維有限元數(shù)值模擬的方法，通過詳細(xì)推導(dǎo)三維地電場(chǎng)變分問題所滿足的變分方程，在直流電法三維有限元正演算法[3]的基礎(chǔ)上加入COLE-COLE 模型，利用數(shù)字濾波求解γ 函數(shù)，可實(shí)現(xiàn)時(shí)間域激發(fā)極化法三維正演計(jì)算。在正演計(jì)算時(shí)，采用六面體矩形網(wǎng)格對(duì)變分方程進(jìn)行離散，發(fā)射源附近適當(dāng)加密；有限元形成的大型稀疏對(duì)稱方程組采用CSR 存儲(chǔ)格式進(jìn)行存儲(chǔ)，大大節(jié)約了內(nèi)存空間?？紤]到現(xiàn)在的計(jì)算機(jī)大多都是多核多線程的，為了使計(jì)算機(jī)資源得到最大限度的利用，使求解方程組的速度得到提高，所以采用并行求解技術(shù)。而MKL 庫(kù)中有許多現(xiàn)成的并行求解算法可以利用，其中并行求解器PARDISO，不僅采用CSR 存儲(chǔ)格式，而且為多線程多核心并行計(jì)算，可以大大提高時(shí)域激發(fā)極化發(fā)三維正演的計(jì)算效率和速度。綜上所述步驟就可以實(shí)現(xiàn)時(shí)域激發(fā)極化發(fā)有限元三維正演。

1.2 精度驗(yàn)證

首先計(jì)算兩正負(fù)點(diǎn)源在均勻半空間中產(chǎn)生的場(chǎng)值，并與解析解做對(duì)比，來驗(yàn)證程序的正確性。圖1 為計(jì)算點(diǎn)到源點(diǎn)距離隨相對(duì)誤差分布圖。從圖可以看出隨著計(jì)算點(diǎn)與源點(diǎn)距離的減小誤差逐漸增大，在距離源點(diǎn)4m 處，相對(duì)誤差小于1.6%。所以當(dāng)計(jì)算點(diǎn)距離源點(diǎn)大于4m 的地方，計(jì)算結(jié)果是正確，精度是可靠的。

圖1 測(cè)點(diǎn)離源的距離的誤差分布圖

1.3 高低阻異常體模型試算

采用中梯裝置，建立如圖2 的高、低阻異常模型。高阻異常體參數(shù)：圍巖電阻率10Ω·m、極化率0.01、時(shí)間常數(shù)3.0、頻率相關(guān)系數(shù)0.1；異常體電阻率分別為20Ω·m、50Ω·m、100Ω·m、500Ω·m，極化率0.3，頻率相關(guān)系數(shù)0.3，時(shí)間常數(shù)5.0。

圖2 模型示意圖

從圖3 高阻異常視極化電阻率分布圖可以看出，在視極化電阻率圖上，高阻異常體在它的正上方會(huì)產(chǎn)生一個(gè)視極化電阻率高阻異常，且隨著異常體與圍巖電阻率差異的增大，它的視極化電阻率異常幅值越大。

圖3 高阻異常視極化電阻率分布圖

低阻異常體與圍巖參數(shù)：圍巖電阻率200Ω·m、極化率0.01、時(shí)間常數(shù)3.0、頻率相關(guān)系數(shù)0.1；異常體電阻率10Ω·m、50Ω·m、100Ω·m、150Ω·m，極化率0.3、頻率相關(guān)系數(shù)0.1、時(shí)間常數(shù)分別為3.0。從圖4 低阻異常視極化電阻率分布圖可以看出，低阻高極化異常體在它的正上方產(chǎn)生低的視極化電阻率異常，且隨著異常體與圍巖的電阻率差異越大，它的視極化電阻率異常越明顯，與高阻異常有相同的規(guī)律。

高、低阻模型的時(shí)間域有限元三維模型的響應(yīng)特征說明本文采用的正演方法是合理的，有效的。

圖4 低阻異常視極化電阻率分布圖

1.4 速度對(duì)比

為了對(duì)比程序速度上的優(yōu)勢(shì)，特意編寫了用定半帶寬存儲(chǔ)，直接采用LDTD 分解法求解方程組的程序和用CSR 存儲(chǔ)，用預(yù)條件共軛梯度法（PCG）解方程的程序來作對(duì)比。在同一臺(tái)計(jì)算機(jī)（雙核四線程，內(nèi)存2G）上，對(duì)三個(gè)程序的計(jì)算速度進(jìn)行對(duì)比，見表1。

表1 不同計(jì)算方法CPU 計(jì)算時(shí)間對(duì)比表

可以看出本文所采用的CSR 存儲(chǔ)格式和Pardiso 解方程的方法比其他兩種方法的計(jì)算速度高。相比定半帶寬存儲(chǔ)模式，內(nèi)存需求少；相比PCG 法，計(jì)算速度優(yōu)勢(shì)明顯，且隨著剖分節(jié)點(diǎn)數(shù)和計(jì)算機(jī)核數(shù)的增加，計(jì)算速度優(yōu)勢(shì)越明顯。

2 結(jié)論

通過詳細(xì)推導(dǎo)三維地電場(chǎng)變分問題所滿足的變分方程，并對(duì)方程進(jìn)行有限元離散分析，得到要求解的大型稀疏對(duì)稱方程組。采用一維非零元素壓縮存儲(chǔ)的CSR 模式和PARDISO 并行求解器求解方程。設(shè)計(jì)了均勻半空間模型，用解析解與數(shù)值解相互擬合，證明了該方法是正確的、精度是可靠的。采用中梯裝置試算了高、低阻異常模型，計(jì)算了異常體的視極化電阻率，證明該方法是合理的、有效的。最后將三種計(jì)算方法的速度進(jìn)行了比較，在保證精度和較小的存儲(chǔ)需求的前提下，本文的方法顯著地提高了計(jì)算速度，為解決激發(fā)極化法三維反演問題提供了一定的基礎(chǔ)。