駱紅梅

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

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

1 正演方法及驗證

1.1 正演方法

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

1.2 精度驗證

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

圖1 測點離源的距離的誤差分布圖

1.3 高低阻異常體模型試算

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

圖2 模型示意圖

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

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

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

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

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

1.4 速度對比

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

表1 不同計算方法CPU 計算時間對比表

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

2 結(jié)論

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