基于格點型有限體積法的直流電測深2.5維正演模擬

歐陽黎明，趙煊煊，童孝忠，謝維，李歡

(1. 中南大學 地球科學與信息物理學院，湖南 長沙，410083；2. 湖南省國土空間調查監測所，湖南 長沙，410129；3. 山西省煤炭地質物探測繪院有限公司，山西 晉中，030600)

直流電測深法勘探是根據地殼中不同巖礦石之間的電阻率或電導率差異，通過在地表觀測天然或人工建立的電場分布來尋找不同巖體結構、調查地質構造以及解決各類工程地質問題[1]。直流電測深勘探方法在探測幾十米至幾百米深度的目標體時具有獨特的優勢，如勘探儀器操作簡單、野外勘測成本低廉及定量反演解釋效果較好等。為了利用直流電測深法查明地下構造及巖體的空間展布特征，這需要對直流電測深勘探數據進行定性解釋與定量解釋，其中包括正演數值模擬與反演成像。直流電測深法勘探的正反演解釋需要以勘探地質對象的實際電性結構變化為基礎。當前，直流電測深法的一維、二維正反演技術是實測資料處理與解釋的主要技術[2-3]。

直流電測深正演數值模擬是基于地電模型的電性分布求解場分布，目前，成熟的直流電測深正演數值模擬方法主要有有限差分法、有限單元法和有限體積法。有限差分法(finite difference method)是一種近似的數值模擬方法，其正演計算的基本原理是采用差商來代替微商或偏微商，將偏微分方程定解問題轉化為適定的線性方程組來求解，這可很好地近似表征地下介質的電性結構[4-7]。有限單元法(finite element method)是求解地球物理波場變分問題的一種數值模擬近似方法，通過將偏微分方程定解問題轉換為變分方程，然后將被積分區域離散化，再按照一定的規律和需要以及異常體的響應精度進行網格離散，最終也能將偏微分方程定解問題轉化為適定的線性方程組來求解[8]。有限單元法已廣泛應用于直流電測深正演模擬，取得了較好的數值計算效果[9-12]。

有限體積法(finite volume method)也稱控制體積法，其基本求解思路是將整個計算區域剖分為有限個連續的單元，并且在每個離散網格節點的周圍都按一定的規則形成一個互不重復的控制體，最后將控制方程對每一個控制體進行積分運算，從而將偏微分方程定解問題轉化為適定的線性方程組來求解[13]。有限體積法的求解思路簡單，通過離散控制方程和定解條件將微分方程轉化為代數方程，遵循因變量在控制體積內的守恒原理。無論是常微分方程、偏微分方程、各種類型的二階線性方程還是高階線性或非線性方程，均可利用有限體積法將方程求解轉化為線性方程組，而后利用計算機求其數值解[14]。近年來，有限體積法在地球物理勘探中逐漸得到應用，眾多學者研究了地球物理正演模擬問題[15-17]。通常，有限體積法構造控制體存在格點型(vertex-centered)和格心型(cell-centered)共2 種方式[18]。格心型方式將單一的網格單元作為控制體單元，并假設單元內電位為常數，這需要進一步近似地面-空氣界面節點的直流電測深離散電位[19]，必然會給直流電測深的電位模擬帶來計算誤差。而格點型方式以網格節點為中心形成控制體單元，并將未知電位置于網格離散節點處，從而可容易地計算地表處的電位，減小數值計算誤差。

本文從點電源電位所滿足的邊值問題出發，建立點源2.5維直流電測深的格點型有限體積正演算法，以便為后續的時間域激電和復電阻率法的正反演研究提供數值計算工具。通過試算存在解析解的均勻半空間地電模型和層狀介質地電模型，驗證格點型有限體積正演算法的準確性和穩定性，并分析直流電測深2.5維正演模擬的誤差來源。利用格點型有限體積數值算法計算典型二維地電模型的視電阻率響應，并總結直流電測深的異常響應規律和特點，以便為實測數據的定性解釋提供指導。

1 邊值問題

直流電測深法勘探通常采用點電源供電，而點電源產生的電位具有三維分布的特征，即u=u(x，y，z)，其滿足如下三維變系數泊松方程邊值問題[20]：

其中：u為點電源產生的電位；δ(x)為狄拉克函數；ρ為點電源處的電荷密度；ΓS為地表-空氣界面，即邊值問題的上邊界；Γ∞為截斷邊界，即無窮遠邊界；r為點電源S到截斷邊界Γ∞的距離；σ為電導率(電阻率的倒數)，其單位為S/m；cos(r，n)為矢徑r與邊界外法線方向n的夾角余弦。

對于二維地電模型，若地下介質的電性參數沿走向y不發生變化，即σ=σ(x，z)，為了消除走向坐標y，需要對式(1)中的偏微分方程進行傅里葉變換。由于電位u(x，y，z)是實函數，并且是變量y的 偶 函 數，即u(x，y，z) =u(x，-y，z)，所 以，對u(x，y，z)進行余弦傅里葉變換，且積分區間選擇為0至+∞，則余弦傅里葉變換公式為

其中：U(λ，x，z)為空間域電位u(x，y，z)經余弦傅里葉變換后的波數域電位；λ為波數或傅里葉變換變量。

對式(1)中的偏微分方程和邊界條件同時進行傅里葉余弦變換，可得波數域電位U(λ，x，z)滿足如下帶參數(波數λ)的二維變系數亥姆霍茲方程邊值問題[21]：

式中：Q為穩恒電流密度；K1和K0分別為第二類的1階、0階修正貝塞爾函數。

2 有限體積正演算法

2.1 控制方程離散

采用矩形網格單元將直流電測深2.5維正演的二維計算區域進行離散，如圖1所示。令

圖1 二維地電模型離散化Fig.1 Discretization for two-dimensional geo-electric model

二維離散區域的內部節點示意圖如圖2 所示。取二維地電模型離散化后網格內部的任意一個節點，按照有限體積法的基本思想，在控制容積中對式(3)中的偏微分方程兩端進行積分：

圖2 二維離散區域的內部節點示意圖Fig.2 Internal nodes of two-dimensional discretized region

2.2 邊界條件處理

上邊界節點滿足第二類邊界條件方程，其控制單元形成的邊界為線段a—b—IV—V—VI—VII(見圖3)。因此，在該線段上的積分為

圖3 剖分網格上邊界節點示意圖Fig.3 Diagram of nodes located on the top edge of mesh

而

將式(10)和式(11)代入式(5)，可得上邊界節點(不含2個角點)滿足的代數方程為

其中：

下邊界節點滿足第三類邊界條件，其控制單元形成的邊界為線段I—II—III—b—a—VIII，如圖4所示，則在該線段上的積分為

圖4 剖分網格的下邊界節點示意圖Fig.4 Diagram of nodes located on the bottom edge of mesh

且

將式(13)和式(14)代入式(5)，可得下邊界節點(不含兩個角點)滿足的代數方程為

其中：

左右邊界節點均滿足第三類邊界條件方程，因此，對左右邊界節點滿足的方程，采用上面相同的方式進行處理。

對于離散區域的角點處理，以上邊界右端角點為例，給出其處理辦法。上邊界右角點控制單元形成的閉曲線為a—d—VI—VII(見圖5)。因此，在該閉曲線上的積分為

圖5 剖分網格的上邊界右端角點示意圖Fig.5 Diagram of right corner node on the bottom of mesh

且

將式(16)和式(17)代入式(5)，可得上邊界右角點滿足的代數方程為

其中：

2.3 線性方程組求解

將計算區域的所有內部節點和邊界節點進行格心型有限體積法處理，最終得出波數域中的線性方程組為

通過求解該適定線性方程組(19)，即得波數域中所有離散網格節點的電位U。對于多個供電電源，需要對方程(19)進行多次求解，然后對波數域電位進行傅里葉逆變換，便可得到點源激勵下二維地電斷面的空間域電位，并根據相應的觀測裝置計算視電阻率。

采用有限體積法解直流電測深正演問題，最終將形成大型、稀疏、對稱的線性方程組，常采用直接法或迭代法求解。若供電點較多，則采用直接法計算，計算效率較高，但占用內存較大；若供電點較少，則采用迭代法計算，計算效率較高，并且占用內存較少[20-22]。鑒于此，根據地電模型和計算配置情況選擇合理的求解方法，本文采用不完全LU 分解(即上三角與下三角分解)預處理的穩定雙共軛梯度算法(ILU-BICGSTAB 迭代法)[23]，實現直流電測深2.5維正演模擬的有限體積線性方程組求解。

2.4 傅里葉逆變換及視電阻率計算

當求出若干波數域中點源二維直流電測深的電位U(λ，x，z)時，采用傅里葉逆變換就可以求出三維點源空間的電位u(x，y，z)。在剖面(y=0)上，傅里葉逆變換公式為

利用數值積分，將式(20)寫成[24]

其中：N為單位長度內波周的個數，為了獲得高精度的數值解，在實際計算過程中取N=10；，為剖面上的點至點電源處的距離；λi為利用最優化原理得到的離散波數[25]；gi為相應的傅里葉逆變換系數。離散波數λi及相應的傅里葉逆變換系數gi見表1。

表1 離散波數λi及相應的傅里葉逆變換系數giTable 1 Discrete wavenumbers λi and corresponding coefficients of inverse Fourier transform gi

直流電測深2.5 維正演模擬時需計算視電阻率，可以根據電位的模擬結果換算出不同裝置的視電阻率，其計算公式為[26]

其中：K為測量裝置系數；Δu為空間域電位差；I為點源供電電流強度，一般可取為1 A。

圖6 所示為格點型有限體積法正演模擬2.5 維直流電測深響應的基本流程圖，通過正演最終獲得三維源、二維地電模型的視電阻率。

圖6 格點型有限體積正演算法流程Fig.6 Flow-process diagram of vertex-centered finitevolume forward algorithm

3 正演算法測試與模型試算

采用數值計算的MATLAB 軟件，編寫直流電測深2.5維正演模擬的格點型有限體積法程序。本文的數值測試平臺為CPU i5-12400，RAM 16G；正演程序編制系統為MATLAB R2021b，該系統可方便調用0階和1階貝塞爾函數。

3.1 算法測試

3.1.1 均勻半空間模型

為了驗證有限體積正演算法和正演程序的正確性，選取1個均勻半空間地電模型，其電導率為0.01 S/m。采用格點型有限體積法模擬2 個異性點電源產生的電位，取計算區域的長度為200 m、寬度為100 m，且將求解區域離散為100×50 的網格。正負電流源布置于(-20 m，0 m)和(20 m，0 m)處，供電電流為I=1 A，將2.5 維數值計算結果與均勻半空間電位解析解[27]進行對比，如圖7所示。從正演模擬結果可知：利用有限體積法得到的電位數值解與解析解擬合得較好，從而驗證了正演算法和MATLAB 正演程序的正確性；另外，在供電點源附近的模擬結果存在誤差，這可以通過加密點源附近網格來提高模擬精度，也可利用邊界/奇異校正處理方法來提高數值計算精度[25]。

圖7 均勻半空間模型的點電源電位數值解與解析解的對比Fig.7 Comparison of analytical and numerical solutions of the point source potentials in the homogeneous halfspace model

3.1.2 層狀介質模型

為了進一步驗證算法的準確性，構建1個水平層狀地層模型，其模型參數為σ1=0.005 S/m，σ2=0.05 S/m和h1=4m，如圖8所示。分別采用直流電測深一維正演算法[28]和本文建立的2.5維有限體積正演算法計算該2層水平地電模型的視電阻率。選用的測量裝置為溫納測深裝置，供電電極距設計為1.5～75.0 m。

圖8 2層水平層狀模型示意圖Fig.8 Schematic diagram of two-layered horizontal model

一維正演和2.5維正演計算的視電阻率曲線對比見圖9。從圖9 可知視電阻率數值解與解析解較吻合，這說明本文建立的2.5維正演算法模擬的視電阻率精度滿足實用要求。在點源附近和供電極距大于70 m時的誤差較大，但其相對誤差均在2%以內。直流電測深2.5維有限體積正演的誤差主要來自3個方面：1) 供電點源附近奇異性問題帶來的誤差；2) 截斷邊界條件和有限體積離散帶來的誤差；3) 離散傅里葉逆變換數值積分帶來的誤差。

圖9 2層地電模型數值解與解析解對比Fig.9 Comparison of analytical and numerical solutions for two-layered geo-electrical model

3.2 模型試算分析

在直流電測深2.5維正演模擬過程中，無窮遠截斷邊界條件的近似處理會產生數值計算誤差。當勘測點距離截斷邊界越近時，誤差必然會越大，只有足夠遠的截斷邊界方能滿足直流電測深勘探的實際要求精度。但足夠遠的邊界處理方式又勢必增加網格單元的數量，增加求解有限體積線性方程組所需的時間。為了實現無窮遠邊界的仿真模擬，將直流電測深計算區域剖分為目標區域和網格外延區域，進而使第三類齊次截斷邊界不受異常體的影響。目標區域為地質構造體的賦存部分，也是直流電測深勘探數據的地表采集區域，其剖分方式為均勻網格；而網格外延區域用于實現無窮遠邊界仿真，其剖分方式為非均勻網格，網格剖分的步長按大于1的倍數等比遞增，在保證正演計算精度的情況下，減少網格剖分單元數，進而提高正演計算效率。

為了分析典型模型的直流電測深的視電阻率響應，構建1個二維簡單異常體地電模型。在電導率為0.01 S/m的均勻半空間中，存在1個電導率為0.1 S/m 的矩形異常體(長4 m、寬2 m)，其頂部距離地表2 m，具體的模型參數如圖10所示。

圖10 二維地電模型示意圖Fig.10 Schematic diagram of a two-dimensional geoelectrical model

研究區域x和z方向的剖分范圍分別為[-35，35] m 和[-30，0] m，采用1 m 的等間隔單元進行剖分，合計剖分71×31=2 201個節點。外延區域的最小剖分單元邊長為2 m，且外延節點數為10 個。選用的測量裝置為溫納測深裝置，電極距設置為1 m。圖11(a)所示為電阻率ρs有限體積數值解等值線圖。根據封閉的等值線圖，可以定性判別出低阻異常體(或高導異常體)的電阻率分布特征，并能準確判別出異常體的空間展布特征。圖11(b)所示為Res2dmod 軟件[29]計算所得電阻率等值線圖，這與有限體積法的計算結果一致，且兩者的電阻率絕對誤差最大值僅為3.4 Ω·m(見圖11(c))。

圖11 低電阻率異常體的正演模擬結果Fig.11 Forward modeling results of low-resistivity anomalous body

4 結論

1) 從點電源電位滿足的變系數泊松方程邊值問題出發，建立了直流電測深2.5維正演模擬的格點型有限體積算法。格點型有限體積法利用散度定理可直接將變系數亥姆霍茲方程邊值問題轉換成適定線性方程組來求解，原理簡單，易于編程實現。

2) 通過模擬均勻半空間模型中2個異性點電源產生的電位和層狀介質模型的直流電測深視電阻率響應，其結果驗證了格點型有限體積正演算法和MATLAB 正演程序的準確性，這為下一步開展直流電測深2.5維反演成像提供了數值計算工具。

3) 采用格點型有限體積法對典型二維地電模型進行了試算，將試算結果與有限差分數值計算結果進行對比分析，進一步驗證了正演算法的有效性。本文建立的格點型有限體積正演算法不僅完善了目前2.5 維直流電阻率測深數值模擬方法，提高了數據處理和解釋精度，而且該方法具有通用性，對類似地球物理正演問題同樣具有實用價值。