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

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

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

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

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

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

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

1 邊值問題

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

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

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

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

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

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

2 有限體積正演算法

2.1 控制方程離散

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

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

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

圖2 二維離散區(qū)域的內(nèi)部節(jié)點示意圖Fig.2 Internal nodes of two-dimensional discretized region

2.2 邊界條件處理

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

圖3 剖分網(wǎng)格上邊界節(jié)點示意圖Fig.3 Diagram of nodes located on the top edge of mesh

而

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

其中：

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

圖4 剖分網(wǎng)格的下邊界節(jié)點示意圖Fig.4 Diagram of nodes located on the bottom edge of mesh

且

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

其中：

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

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

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

且

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

其中：

2.3 線性方程組求解

將計算區(qū)域的所有內(nèi)部節(jié)點和邊界節(jié)點進行格心型有限體積法處理，最終得出波數(shù)域中的線性方程組為

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 算法測試

3.1.1 均勻半空間模型

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

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

3.1.2 層狀介質(zhì)模型

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

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

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

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

3.2 模型試算分析

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

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

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

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

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

4 結(jié)論

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

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

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