時間域激發極化法三維有限元正演

駱紅梅

（福建省地質測繪院，福建 福州 350011）

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

1 正演方法及驗證

1.1 正演方法

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

1.2 精度驗證

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

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

1.3 高低阻異常體模型試算

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

圖2 模型示意圖

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

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

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

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

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

1.4 速度對比

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

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

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

2 結論

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