基于Gmsh的起伏地形下井—地直流電法正演模擬

張宇哲，孟麟，王智

(長(zhǎng)江大學(xué) 電子信息學(xué)院，湖北 荊州 434023)

通訊作者： 王智(1985-)，男，湖北省武漢市人，博士，碩士生導(dǎo)師，主要研究方向?yàn)殡姶欧〝?shù)值模擬和機(jī)器學(xué)習(xí)。Email:1324385898@qq.com

0 引言

Gmsh是一個(gè)具有CAD內(nèi)核和后處理器的開源三維網(wǎng)格生成軟件，提供了一種帶有參數(shù)輸入和高級(jí)可視化功能的快速、輕便和用戶友好的網(wǎng)格剖分工具[1]，它由幾何建模、網(wǎng)格剖分、求解器和后處理4個(gè)部分組成。自1971年，Coggon[2]首次將有限單元法應(yīng)用到二維線電流場(chǎng)的正演數(shù)值模擬過程后，很多學(xué)者對(duì)利用有限單元法解決直流電法正演問題進(jìn)行了深入研究。如：Rijo[3]引入通用性網(wǎng)格，提高了正演的計(jì)算速度和精度；Dey等[4]對(duì)井—地電法和地表電法的各自裝置類型的特點(diǎn)進(jìn)行了對(duì)比分析；Wu X P[5]將不完全Cholesky共軛梯度方法應(yīng)用在有限單元法的正演模擬中，提升了運(yùn)算速率；Cardarelli等[6]利用艾米特插值進(jìn)行傅里葉反變換，提高了正演結(jié)果的精度；Tang J T等[7]使用非結(jié)構(gòu)化網(wǎng)格來擬合復(fù)雜地形的電阻率模型；Ren Z Y等[8-10]提出自適應(yīng)有限元算法自動(dòng)進(jìn)行網(wǎng)格加密，提高正演的準(zhǔn)確性；Pan K J等[11]將外推瀑布式多重網(wǎng)格法應(yīng)用在2.5D和3D直流電法正演中；Zhang Q J等[12]研究了加密—收縮的網(wǎng)格在有限元模擬中的應(yīng)用；Kuang X T等[13]對(duì)起伏地形條件下長(zhǎng)方體磁場(chǎng)無解析解奇點(diǎn)表達(dá)式進(jìn)行了研究；王智等[14]結(jié)合井—地觀測(cè)方法和改進(jìn)的異常電位法，確定了適用于起伏地表的自然邊界條件；武建平等[15]將有限單元法應(yīng)用在電磁法三維正演模擬中，取得了不錯(cuò)的結(jié)果；嚴(yán)波等[16]將對(duì)偶加權(quán)后驗(yàn)誤差估計(jì)方法應(yīng)用在自適應(yīng)有限元正演模擬中，提高了自適應(yīng)有限元網(wǎng)格加密的有效性。本文采用不規(guī)則網(wǎng)格剖分?jǐn)M合起伏地形，使用井—地聯(lián)合觀測(cè)方法來進(jìn)行起伏地形下的異常響應(yīng)的觀測(cè)，研究使用不同觀測(cè)裝置時(shí)山谷地形對(duì)下方異常響應(yīng)的影響。

1 2.5D有限元正演模擬

在三維無限半空間Ω中，點(diǎn)電源A的坐標(biāo)為(xA,yA,zA)，若要求解直流電法勘探的點(diǎn)源2.5D邊值問題，需要通過傅里葉變換將空間域電位u(x,y,z)轉(zhuǎn)化為波數(shù)域電位U(x,k,z)。波數(shù)域電位U滿足的偏微分方程[17]如下：

(1)

式中：σ是地下介質(zhì)電導(dǎo)率；I是供電電流強(qiáng)度；δ是Delta脈沖函數(shù)；ΓS和?！薹謩e為三維無線半空間的地表邊界和無窮遠(yuǎn)邊界；n代表邊界的外法向向量；r是?！奚先我庖稽c(diǎn)到點(diǎn)電源的向量；K1和K0分別為一階和零階修正貝塞爾函數(shù)；k為波數(shù)。該偏微分方程的邊界條件分別為：①地表邊界采用Newman邊界條件；②無窮遠(yuǎn)邊界采用混合邊界條件。

用有限單元法解決上述邊值問題時(shí)，需推導(dǎo)邊值問題對(duì)應(yīng)的變分問題，式(1)對(duì)應(yīng)的變分問題為：

(2)

得到對(duì)應(yīng)變分問題后，用三角形單元將求解區(qū)域離散化；然后在每個(gè)離散單元e上用多項(xiàng)式函數(shù)和對(duì)應(yīng)的形函數(shù)來近似待求量；最后將每個(gè)單元合并，求解大型稀疏線性方程組，得到各個(gè)節(jié)點(diǎn)處的波數(shù)域電位值U，進(jìn)行傅里葉反變換便可得到空間中的電位值u。本文選用的波數(shù)ki與傅里葉反變換系數(shù)gi選自文獻(xiàn)[17]中采用最優(yōu)化方法得到的5個(gè)波數(shù)和對(duì)應(yīng)的傅里葉反變化系數(shù)。

2 建模及觀測(cè)裝置

在Gmsh建立模型和網(wǎng)格剖分過程中，有2種建模方式：①采用Gmsh軟件左側(cè)的GUI進(jìn)行交互建模；②在geo格式文件中采用Gmsh自己的腳本語句進(jìn)行建模。本文選擇后一種方式，Gmsh腳本建模原理請(qǐng)參考本文附錄A。

井—地電阻率觀測(cè)裝置由井下供電電極和地表測(cè)量電極共同組成，由于井中的觀測(cè)電極距離地下的異常體更近，所以得到的觀測(cè)數(shù)據(jù)量更大且更精確，還可以分別移動(dòng)供電電極和測(cè)量電極的相對(duì)位置，來觀測(cè)和研究不同水平位置和不同深度的地下介質(zhì)情況。根據(jù)實(shí)際使用的觀測(cè)裝置的不同，井—地電阻率觀測(cè)裝置可以細(xì)分為井—地二極觀測(cè)裝置、井—地三極觀測(cè)裝置、井—地偶極觀測(cè)裝置等。在記錄和作圖時(shí)，井—地二極裝置(圖1a)以井中供電電極A的縱坐標(biāo)為記錄點(diǎn)縱坐標(biāo)，地表觀測(cè)電極M的橫坐標(biāo)為記錄點(diǎn)橫坐標(biāo)；井—地三極裝置(圖1b)以井中供電電極A的縱坐標(biāo)為記錄點(diǎn)縱坐標(biāo)，地表觀測(cè)電極MN中點(diǎn)的橫坐標(biāo)為記錄點(diǎn)橫坐標(biāo)。

圖1 觀測(cè)裝置示意[18]Fig.1 Example of observation device[18]

3 數(shù)值模擬及分析

地質(zhì)模型的建立和網(wǎng)格剖分是使用有限單元法解決地球物理正演問題的重要環(huán)節(jié)。用Gmsh軟件對(duì)地質(zhì)模型進(jìn)行不規(guī)則的網(wǎng)格剖分時(shí)進(jìn)行局部加密，能更精確地?cái)M合各種地形情況，提高正演結(jié)果的準(zhǔn)確性。2.5D有限單元正演數(shù)值模擬程序采用CSR格式來壓縮儲(chǔ)存大型稀疏矩陣；選用Eigen庫中的BiCGSTAB(biconjugate gradient stabilized method)求解器解正演過程中的線性方程組。首先，為了驗(yàn)證Gmsh在正演問題中的適用性和有限單元正演數(shù)值模擬算法的正確性，選用具有解析解的層狀模型進(jìn)行正演模擬；然后，使用井—地二極和井—地三極電阻率觀測(cè)裝置進(jìn)行觀測(cè)，對(duì)兩種山谷地形下異常體的異常響應(yīng)進(jìn)行正演模擬，討論山谷地形對(duì)下方異常體的異常響應(yīng)的影響。

3.1 算法驗(yàn)證

水平層狀模型如圖2所示，采用二極裝置進(jìn)行觀測(cè)，得到的電位解析解和數(shù)值解的對(duì)比如圖3所示。從結(jié)果可以看出數(shù)值解和解析解的誤差主要集中在靠近供電電極的位置，其余位置的誤差較小，證明了Gmsh在正演問題中的適用性和有限單元正演數(shù)值模擬算法的正確性。

圖2 水平層狀模型Fig.2 Two-layers geo-electric cross section

圖3 電位數(shù)值解和解析解對(duì)比Fig.3 The curve for analytical solution and calculation solution of potential

3.2 計(jì)算實(shí)例

所設(shè)3個(gè)模型如圖4所示，分別為平坦地形、底部平緩的山谷地形以及底部尖銳的山谷地形下的異常體模型，井中供電電極A的變換范圍為0～20 m，間距1 m；地表觀測(cè)電極的變換范圍為0～50 m，間距1 m;山谷地形底部深度為5 m;異常體頂部埋深為8.5 m，長(zhǎng)4 m，高2 m。為了比較各種地形和裝置的觀測(cè)視電阻率異常，圍巖電阻率設(shè)定為100 Ω·m，低阻體電阻率為5 Ω·m，高阻體電阻率為5 000 Ω·m。在作圖時(shí)，超出測(cè)量點(diǎn)范圍的數(shù)據(jù)由Surfer軟件根據(jù)測(cè)量范圍內(nèi)的數(shù)據(jù)推導(dǎo)畫出。

圖4 模型示意Fig.4 Sketch map of model

采用井—地二極裝置對(duì)模型一進(jìn)行正演模擬所得視電阻率斷面如圖5所示。在圖5a中，低阻異常體的視電阻率斷面中只存在獨(dú)立的橢圓形低阻異常，是由于低阻異常體對(duì)電流的吸引作用而形成；垂直方向上低阻異常中心深度z=8.5 m，與異常體的頂部埋深一致，最低幅值為76 Ω·m，水平方向上低阻異常中心的范圍x=23～27 m，與異常體的寬度相同。圖5b顯示,高阻異常體的視電阻率斷面垂直方向上存在對(duì)稱的高阻異常和低阻異常，當(dāng)源點(diǎn)在異常體上方時(shí)，由于高阻體對(duì)電流的排斥作用，異常體上方地表的電流密度大于正常電流密度，因此呈現(xiàn)出大于背景值的高阻異常，當(dāng)源點(diǎn)在異常體下方時(shí)，由于高阻體的屏蔽電流作用，異常體上方地表的電流密度小于正常電流密度，因此呈現(xiàn)小于背景值的低阻異常；上部高阻異常的最高幅值為120 Ω·m，下部低阻異常的最低幅值為84 Ω·m，異常分界面為z=9.5 m，和異常體中心深度相同;在水平方向上，高阻、低阻異常的分布范圍都為x=23～27 m，與異常體的寬度相同。兩種情況下均可通過視電阻率斷面圖中的異常響應(yīng)準(zhǔn)確推斷異常體位置。

圖5 模型一的視電阻率斷面(井—地二極裝置)Fig.5 Apparent resistivity cross-section view of model 1(borehole-ground pole-pole device)

圖6為模型二對(duì)應(yīng)的純山谷地形的視電阻率斷面。采用井—地二極裝置觀測(cè)時(shí)，在水平方向上視電阻率等值線最為密集的地方在山谷地形和平坦地形的分界處(圖6a)，山谷地形最低處下方存在一個(gè)相對(duì)左右兩側(cè)的高阻異常，異常幅值為95 Ω·m；在垂直方向上，相對(duì)左右兩側(cè)的高阻異常的頂部和山谷地形的底部的深度平齊。采用井—地三極裝置，在水平方向從左到右依次為低阻—高阻—低阻趨勢(shì)(圖6b)，兩側(cè)的低阻脈沖異常位于山谷地形和平坦地形的分界處，最低異常幅值為46 Ω·m；高阻異常位于山谷地形最低處下方， 最高異常幅值為106 Ω·m；在垂直方向上，兩側(cè)低阻脈沖異常最低幅值等值線頂部高度和山谷地形的底部的深度相等。

圖6 模型二的視電阻率斷面Fig.6 Apparent resistivity cross-section view of model 2

模型二正演結(jié)果的視電阻率斷面如圖7所示。采用井—地二極觀測(cè)裝置，低阻異常體的視電阻率斷面圖(圖7a)上存在一個(gè)獨(dú)立的橢圓形低阻異常，和平坦地形下的低阻異常相比，該異常的范圍明顯縮小，其中心深度位于z=8.5 m，與異常體的頂部埋深一致，最低幅值為54 Ω·m。山谷地形底部存在相對(duì)左右兩側(cè)的高阻異常，異常中心深度和山谷地形底部深度一致，幅值為83 Ω·m；低阻異常中心的范圍為x=23～27 m，與異常體的寬度相同。高阻異常體的視電阻率斷面圖(圖7b)中高阻、低阻并存，和平坦地形下的異常響應(yīng)相比，上部分高阻異常的范圍有所縮小;高阻異常和低阻異常不完全對(duì)稱,上部高阻異常的最高幅值為128 Ω·m，下部低阻異常的最低幅值為66 Ω·m,異常分界線為z=9.5 m，和異常體中心深度相同;在水平方向上，高阻、低阻異常的分布范圍都在x=23～27 m中，與異常體的寬度相同。

山谷地形下異常體模型斷面圖中，低阻異常體所產(chǎn)生的低阻異常區(qū)域基本和低阻異常體模型重合，高阻異常區(qū)域中心位于山谷底部最深處；高阻異常體產(chǎn)生的高阻異常區(qū)域和低阻異常區(qū)域基本對(duì)稱，分界線和異常體中心埋深一致，據(jù)此可以判斷山谷地形下方為低阻異常體還是高阻異常體。

圖7c、d為采用井—地三極觀測(cè)裝置所得。圖7c中可見，低阻異常體的視電阻率斷面存在山谷地形下方的高阻異常和獨(dú)立的橢圓形低阻異常，高阻異常中心位于山谷地形的最低處，最高幅值為105 Ω·m；低阻異常中心深度z=9.5 m，與異常體中心的埋深平齊，高、低阻異常分界線為z=8.5 m，與異常體的頂部埋深一致；在水平方向上，低阻異常中心的范圍與異常體的寬度相同，位于山谷地形和平坦地形的分界處的低阻脈沖異常依然明顯存在。圖7d顯示高阻異常體的視電阻率斷面中同樣存在高阻和低阻異常，高阻異常中心深度z與異常體的頂部埋深一致，其電阻率最高幅值為180 Ω·m，低阻異常最低為10 Ω·m，高、低阻異常的分界線為z=10.5 m，與異常體的底部深度相同；在水平方向上，兩側(cè)的低阻脈沖異常位于山谷地形和平坦地形的分界處，高、低阻異常中心的范圍均為x=23～27 m，與異常體的寬度相同。

圖7 模型二的正演結(jié)果Fig.7 Apparent resistivity cross-section view of the forward result of model 2

圖8為模型三對(duì)應(yīng)的純山谷地形的視電阻率斷面。采用井—地二極裝置觀測(cè)時(shí)，視電阻率斷面存在方向向下的低阻脈沖異常(圖8a),從左到右依次為高阻—低阻—高阻趨勢(shì)，等值線最密集的位置在山谷地形和平坦地形的分界處;低阻異常位于山谷地形下方，最低幅值為76 Ω·m，向下的低阻脈沖異常的底部和山谷地形的底部的深度平齊。采用井—地三極裝置觀測(cè)時(shí)，視電阻率斷面與井—地二極裝置觀測(cè)得到的相反，存在方向向上的脈沖異常(圖8b)，異常范圍x=10～40 m，與山谷地形范圍一致;低阻異常位于山谷地形正下方，其頂部和山谷地形的底部的深度平齊，低阻異常最低幅值為25 Ω·m。

模型三正演結(jié)果的視電阻率斷面如圖9所示。圖9a、b為采用井—地二極觀測(cè)裝置所得。低阻異常體的視電阻率斷面(圖9a)中存在一個(gè)獨(dú)立的橢圓形低阻異常,異常中心深度與異常體的頂部埋深一致，位置與異常體中心水平位置相同，最低幅值為48 Ω·m。高阻異常體的視電阻率斷面(圖9b)中高阻、低阻異常并存，高阻異常的最高幅值為108 Ω·m，低阻異常的最低幅值為60 Ω·m，異常分界線和異常體中心深度相同；高阻異常的中心由于山谷地形的影響被阻斷，低阻異常的中心位于x=25 m處，與異常體的水平中心位置平齊。

圖8 模型三的視電阻率斷面Fig.8 Apparent resistivity cross-section view of model 3

圖9 模型三的正演結(jié)果Fig.9 Apparent resistivity cross-section view of the forward result of model 3

圖9c、d 為采用井—地三極觀測(cè)裝置所得。在低阻異常體的視電阻率斷面(圖9c)中，垂直方向上有存在一個(gè)橢圓形低阻異常和一個(gè)脈沖異常。在垂直方向上，橢圓形低阻異常中心位于z=9.5 m，與異常體中心埋深相同；純山谷地形顯示的低阻脈沖異常頂部應(yīng)在z=5 m處，由于存在低阻異常體的影響，低阻脈沖異常的頂部位置下降到z=12 m處；在水平方向上，低阻異常和低阻脈沖異常中心均為x=25 m，最低幅值為15 Ω·m，低阻脈沖異常的分布范圍仍與山谷地形范圍一致。高阻異常體的視電阻率斷面(圖9d)中高阻、低阻異常并存，在垂直方向上，上面的高阻異常中心深度z=8.5 m，與異常體的頂部埋深一致，其最高幅值為115 Ω·m，下面是低阻異常最低幅值為10 Ω·m，高、低阻異常的分界線z=9.5 m，與異常體的中心深度相同；在水平方向上，高阻異常在x=25 m被阻斷，下面的低阻脈沖異常形態(tài)相對(duì)于純山谷地形的的異常響應(yīng)情況來說有所變化，但影響范圍仍然在x=10～40 m，與山谷地形范圍一致。

4 結(jié)論

采用Gmsh軟件進(jìn)行建模和網(wǎng)格剖分，使用井—地二極和井—地三極電阻率裝置進(jìn)行觀測(cè)，用有限單元法對(duì)2種山谷地形下異常體的異常響應(yīng)進(jìn)行正演模擬后，對(duì)正演模擬結(jié)果進(jìn)行分析，得到了以下結(jié)論：

1)底部尖銳的純山谷地形比底部平緩的純山谷地形產(chǎn)生的異常響應(yīng)更加明顯和突出，對(duì)下方異常體的異常響應(yīng)的影響也更大。尖銳山谷地形產(chǎn)生的影響主要集中在尖銳處水平2 m內(nèi)，且異常幅度大，平緩山谷地形的影響范圍較大而異常幅度小。

2)井—地三級(jí)裝置觀測(cè)結(jié)果比井—地二級(jí)裝置觀測(cè)結(jié)果的異常響應(yīng)幅值更大，異常形態(tài)更復(fù)雜。就水平分辨率來說，井—地三極裝置的分辨率比井—地二極更強(qiáng)；2種裝置的垂向分辨率相差不大，但是異常形態(tài)各有特點(diǎn)，在實(shí)際工程中應(yīng)注意區(qū)分。

3)采用不規(guī)則網(wǎng)格剖分?jǐn)M合起伏地形和使用井—地聯(lián)合觀測(cè)方法來進(jìn)行起伏地形下的地質(zhì)情況勘探能得到較好的結(jié)果，同時(shí)也證明有限元軟件Gmsh在地球物理有限單元法正演建模和網(wǎng)格剖分方面有良好的應(yīng)用價(jià)值。