王新宇，王 程，毛玉蓉，嚴(yán)良俊，周 磊，高文龍

(1.長江大學(xué) 油氣資源與勘探技術(shù)教育部重點實驗室，湖北 武漢 430100；2.非常規(guī)油氣省部共建協(xié)同創(chuàng)新中心，湖北 武漢 430100；3.中煤科工集團西安研究院有限公司，陜西 西安 710077)

直流電阻率法從地表到孔中、低維到高維、粗放到精細快速發(fā)展，在有色金屬、煤田等礦產(chǎn)資源勘探和地質(zhì)調(diào)查中應(yīng)用廣泛、效果明顯[1-2]。直流電阻率法實測資料的合理解釋離不開有效的三維正反演技術(shù)，并且直流電阻率法三維正反演一直是地球物理學(xué)界的研究熱點[3-8]。

國內(nèi)外學(xué)者對電阻率法正演做了大量研究。20 世紀(jì)70 年代末期，國外學(xué)者A.Dey 等[9]將吸收邊界條件的思想應(yīng)用于電阻率法三維數(shù)值模擬，有效減少網(wǎng)格數(shù)目和邊界效應(yīng)的影響，取得良好的模擬效果；R.C.Fox 等[10]研究了起伏地形對電阻率法和激發(fā)極化法觀測數(shù)據(jù)的影響；T.Lowry[11]、ZHAO Shengkai[12]等將電位分解成背景值和異常值的疊加，提出了異常電位法，克服場源的奇異性，有效地提高計算精度；C.Rücker 等[13]采用非結(jié)構(gòu)四面體網(wǎng)格有限元法實現(xiàn)了復(fù)雜地形下電阻率法三維模擬；強建科等[14]將三棱柱單元引入起伏地形條件下電阻率法有限元三維正演中，并研究了起伏地形對電阻率法正演響應(yīng)的影響特征；湯井田等[15]采用三棱柱網(wǎng)格，推導(dǎo)了無限元三維單元映射函數(shù)，取得較高的計算精度。在井-地電阻率法數(shù)值模擬方面，王志剛等[16]基于積分方程法離散泊松方程，實現(xiàn)了井-地電阻率法三維并行正演算法；部分學(xué)者[17-18]基于有限差分法實現(xiàn)了線源井-地電法三維正演，并取得了良好的數(shù)值模擬效果；李長偉等[19]基于三棱柱網(wǎng)格實現(xiàn)了井-地電法有限元三維正演；王智等[20]對井-地電阻率法三維正演的結(jié)果采用歸一化總水平導(dǎo)數(shù)法，一定程度增強了異常體的響應(yīng)特征。

上述研究在進行正演模擬時，均未考慮鉆孔因素(井液電阻率和鉆孔孔徑)對視電阻率數(shù)據(jù)的影響。而考慮鉆孔的地電模型，網(wǎng)格剖分勢必變得困難。為解決上述問題，采用有限元法，結(jié)合混合網(wǎng)格的空間離散模式，解決小孔徑條件下“空間離散困難”的問題。對于鉆孔所在地層，采用三棱柱網(wǎng)格對鉆孔進行快速離散；對于其他地層，采用非結(jié)構(gòu)化四面體網(wǎng)格進行離散，以減少有限元的計算量。針對不同的地電模型，采用合適的網(wǎng)格類型進行離散，并開展高精度、快速的三維正演模擬。

1 混合網(wǎng)格有限元三維正演算法

1.1 有限元邊值問題及變分問題

直流電阻率法有限元三維正演中，異常電位us滿足的邊界問題[21]為：

式中：σ為地下介質(zhì)電導(dǎo)率，S/m；σ′為異常電導(dǎo)率，σ′=σ?σ0，S/m；σ0為均勻介質(zhì)電導(dǎo)率，S/m；us為異常電位，V，由地下的電性不均勻體產(chǎn)生；u0為正常電位，V，是點電流源在均勻半空間或全空間產(chǎn)生的電位，可由解析解求得；n為 模型邊界的外法向方向；r和r分別為邊界上的點到點電流源的距離和有向距離，m；Γ0和Γ∞分別為計算區(qū)域 Ω的地面邊界和無窮遠邊界。

則異常電位滿足的變分問題[21]為：

式中：dΩ為體積分；dΓ為計算區(qū)域的無窮遠邊界面積分。

1.2 混合網(wǎng)格單元分析

將計算區(qū)域 Ω離散為互不重疊的四面體和三棱柱單元，四面體和三棱柱節(jié)點編號規(guī)則如圖1 所示：

圖1 四面體和三棱柱的網(wǎng)格節(jié)點編號Fig.1 Nodes number of tetrahedral and triangular prism element

則第e個單元內(nèi)任意一點的異常電位為：

式中：對于四面體單元t=4，對于三棱柱單元t=6；(x,y,z)為第e個單元內(nèi)任意一點的異常電位；為第e個單元上第i個節(jié)點的異常電位。

四面體的線性插值基函數(shù)(x,y,z)可表示為：

式中：i為四面體節(jié)點索引；Ve為四面體單元體積；ai、bi、ci、di為四面體線性插值基函數(shù)系數(shù)(參照文獻[21]計算)。

正三棱柱的線性插值基函數(shù)[14]可表示為：

式中：i為三棱柱節(jié)點索引；Δ為三棱柱三角形面積；h為三棱柱單元高度(h=2s)；z0為三棱柱中心z坐標(biāo)；ai、bi、ci(i=k,m,n,k′,m′,n′)為三棱柱線性插值基函數(shù)系數(shù)。

采用有限元法對式(2)進行離散后，對每個單元進行積分得到單元系數(shù)矩陣，最終形成總體系數(shù)矩陣：

式中：Ne為混合網(wǎng)格單元總數(shù)；Ke=σ(K1e+K2e)，Ke為與 σ相關(guān)的單元系數(shù)矩陣；Pe=σ′(K1e+K2e)，Pe為與σ′相關(guān)的單元系數(shù)矩陣；K1e為式(2)中與體積分相關(guān)的單元系數(shù)矩陣的擴展矩陣；K2e為式(2)中與邊界積分相關(guān)的單元系數(shù)矩陣的擴展矩陣；K為總電位的系數(shù)矩陣；P為異常電位的系數(shù)矩陣。

對于每個四面體單元，擴展矩陣K1e和K2e表示為：

式中：Δ1為四面體或三棱柱的三角形面積。

對于每個正三棱柱單元，擴展矩陣K1e表示為：

當(dāng)三棱柱側(cè)面為無窮遠邊界時，擴展矩陣K2e為：

式中：Δ2為三棱柱外邊界四邊形的面積。當(dāng)三棱柱側(cè)面(四邊形)位于無窮遠邊界時，K2e積分采用式(13)計算；當(dāng)三棱柱底面和四面體表面(三角形)位于無窮遠邊界時，K2e積分采用式(8)計算。至此，四面體與三棱柱單元積分和邊界積分完成，由于采用標(biāo)量有限元，四面體與三棱柱共節(jié)點處的積分點重合；如圖1 所示，當(dāng)四面體單元2、3、4 節(jié)點與三棱柱單元k′、m′、n′重合時，積分點重合的單元自動滿足標(biāo)量場的一致性。令式(6)的變分為零，得到Ax=b形式的線性方程組：

1.3 線性方程組求解

本文采用改進的對稱逐步超松弛預(yù)處理共軛梯度迭代算法(SSOR-PCG)[22]求解式(14)。引入對稱逐步超松弛迭代(SSOR 預(yù)條件)的分裂矩陣作為預(yù)處理矩陣：

式中：H為預(yù)處理矩陣；D為A的對角陣；L為A的嚴(yán)格下三角矩陣；0<ω<2為松弛因子(算法細節(jié)參照文獻[22])。

1.4 視電阻率計算

在求解方程后得到各節(jié)點的異常電位值，從而得到各節(jié)點的總電位值：

式中：u為每個節(jié)點的總電位。

本文主要研究井-地、地-井工作方式的視電阻率響應(yīng)特征，圖2 為井-地、地-井工作方式的示意圖。二極裝置與三極裝置的視電阻率計算方式分別如下：

圖2 電阻率法井-地、地-井工作方式Fig.2 Surface-to-borehole and borehole-to-surface working mode of resistivity method

式中：ρs為視電阻率，Ω·m；I為電流大小，A；KAM和KAMN分別為二極裝置與三極裝置的裝置系數(shù)；r和r′分別為點電源A和其關(guān)于地表鏡像點到測點M或N的距離。

2 混合網(wǎng)格有限元電阻率法數(shù)值模擬分析

2.1 混合網(wǎng)格算法精度驗證

采用一維層狀地層模型(圖3) 驗證混合網(wǎng)格算法的正確性。計算區(qū)域大小為1 000 m×1 000 m×1 000 m，第一層地層的厚度為6 m，電阻率為50 Ω·m；第二層地層的厚度為10 m，電阻率為100Ω·m；基底電阻率為50 Ω·m。圖4a 為混合網(wǎng)格剖分示意圖，頂部兩層地層均采用三棱柱網(wǎng)格剖分，基底地層采用四面體網(wǎng)格剖分，最終生成170 360 個三棱柱，294 820 個四面體以及57 948 個邊界三角形和3 200 個邊界四邊形。圖4b 為非結(jié)構(gòu)四面體網(wǎng)格剖分示意圖，最終生成473 222 個四面體，10 348 個邊界三角形。發(fā)射源置于地表(0 m,0 m,0 m)處，發(fā)射電流為1 A，在x方向1～100 m 的范圍內(nèi)布置100 個等間隔測點。采用SSORPCG 算法求解線性方程組，混合網(wǎng)格求解耗時11.02 s，四面體網(wǎng)格求解耗時10.75 s；由解析解和數(shù)值解的對比結(jié)果(圖5)可知：混合網(wǎng)格數(shù)值解的最大誤差為0.15 %，四面體網(wǎng)格數(shù)值解的最大誤差為1.56%(誤差計算方法如下式)，2 種方法均滿足計算精度。但混合網(wǎng)格的計算精度要略高于四面體網(wǎng)格，這是由于在較薄的地層離散網(wǎng)格時會產(chǎn)生質(zhì)量較差的四面體網(wǎng)格，從而影響計算精度，這也證實本文混合網(wǎng)格算法的可行性、有效性。

圖3 一維層狀模型Fig.3 1D layered model

圖4 層狀模型的網(wǎng)格剖分Fig.4 Grid discretization of layered model

圖5 不同方法的結(jié)果對比Fig.5 Comparison of forward modeling results of different methods

式中：ρa為混合網(wǎng)格或四面體網(wǎng)格三維數(shù)值解，Ω·m；ρb為一維解析解，Ω·m；E為數(shù)值解誤差，%。

2.2 鉆孔模型井-地觀測方式數(shù)值模擬分析

為研究鉆孔因素對井-地電阻率法視電阻率響應(yīng)的影響特征，設(shè)計鉆孔模型(圖6 為鉆孔模型的網(wǎng)格剖分示意圖)。鉆孔孔徑為10 cm，鉆孔深度為50 m，計算區(qū)域大小為1 000 m×1 000 m×1 000 m，鉆孔內(nèi)存在電阻率為3 Ω·m的井液，圍巖電阻率為100 Ω·m。在鉆孔周圍采用三棱柱網(wǎng)格進行離散，井底下方采用四面體網(wǎng)格進行離散，混合網(wǎng)格最終生成265 145 個節(jié)點，497 200 個三棱柱，76 587 個四面體以及55 496 個三角形和6 160 個四邊形。發(fā)射源分別位于井口(0 m,0 m,0 m)，井中(0 m,0 m,5 m)、(0 m,0 m,10 m)、(0 m,0 m,15 m)處，發(fā)射電流均為1 A，采用井-地觀測方式計算地表測點的視電阻率響應(yīng)。數(shù)值模擬結(jié)果如圖7所示，將視電阻率的色標(biāo)限制在圍巖電阻率的上下5%范圍內(nèi)，超出此范圍，即可認(rèn)為鉆孔因素對視電阻率數(shù)據(jù)的影響不可忽略。由圖7a 可以看出，當(dāng)發(fā)射源位于地表時，鉆孔附近視電阻率受鉆孔因素影響最大，但隨著測點收發(fā)距的增加，視電阻率響應(yīng)逐漸趨于圍巖電阻率。如圖7b 所示，當(dāng)發(fā)射源位于地下5 m 時，地表視電阻率受鉆孔因素的影響減小，但仍然不能忽略。當(dāng)發(fā)射源位于地下10 m(圖7c)，15 m(圖7d)時，地表視電阻率響應(yīng)和圍巖電阻率非常接近，可以忽略鉆孔因素對視電阻率數(shù)據(jù)的影響。因此，在實際井-地電阻率法資料處理解釋中，根據(jù)發(fā)射源位置，合理考慮鉆孔因素的影響是必要的。

圖6 鉆孔模型的局部網(wǎng)格剖分Fig.6 Local grid of borehole models

圖7 發(fā)射源位于不同深度的視電阻率響應(yīng)Fig.7 Apparent resistivity response of point current source at different depth

2.3 鉆孔模型地-井觀測方式數(shù)值模擬分析

為研究鉆孔因素(井液電阻率和鉆孔孔徑)對地-井觀測方式視電阻率響應(yīng)曲線的影響，設(shè)計鉆孔地電模型(圖8)。鉆孔深度為100 m，圍巖采用層狀地層模型，該模型的地層厚度及電阻率參數(shù)見表1，異常體的大小為30 m×20 m×20 m，電阻率為10 Ω·m。發(fā)射源位于地表井口(0 m,0 m,0 m)處，發(fā)射電流為1 A，接收電極位于井中2～100 m，電極間距2 m，以三極裝置方式進行接收。分別設(shè)計如下2 個模型：①鉆孔孔徑為50 mm 時，計算不同井液電阻率(3、10、50 Ω·m)情況下的視電阻率響應(yīng)曲線；② 井液電阻率為3Ω·m時，計算不同鉆孔孔徑(0、50、100、150 mm)情況下的視電阻率響應(yīng)曲線。

圖8 鉆孔模型混合網(wǎng)格離散Fig.8 Mixed grid discretization of borehole model

表1 鉆孔模型的地層參數(shù)Table 1 Formation parameters of borehole model

圖9a 展示了地-井觀測方式在不同井液電阻率情況下的視電阻率響應(yīng)曲線。由圖9a 可知，淺部測點的視電阻率數(shù)據(jù)受井液電阻率影響較大，深部測點的視電阻率數(shù)據(jù)基本不受井液電阻率影響；而且隨著井液電阻率的增大，視電阻率曲線所受影響逐漸減小。實際鉆孔中的井液電阻率均為低阻，通過分析井液電阻率為3 Ω·m時的視電阻率響應(yīng)可知：當(dāng)測點收發(fā)距為10 m 左右時，井液的影響才可以忽略不計。因此，在實際地-井電阻率法勘探中，鉆孔因素應(yīng)加以重視。圖9a 中的視電阻率曲線很難反映地層的電性特性和界面信息，這是由于地下介質(zhì)存在明顯的電阻率差異，高阻地層對電流具有排斥作用、低阻異常和低阻地層對電流具有吸引作用，地層與異常體之間的電流場相互耦合，使得視電阻率響應(yīng)變得尤為復(fù)雜，很難對地下介質(zhì)的電阻率信息進行有效判別。因此，在復(fù)雜地層環(huán)境下，地-井觀測方式的視電阻率曲線變化較為復(fù)雜，給直流電阻率法資料的處理解釋工作帶來巨大困難。

圖9 不同鉆孔參數(shù)下視電阻率響應(yīng)曲線Fig.9 Apparent resistivity response curves of different borehole parameters

現(xiàn)將所有地層的電阻率設(shè)為100Ω·m，異常體的電阻率為10Ω·m，在忽略鉆孔的情況下，研究地下低阻異常體對視電阻率曲線的影響特征。圖9b 展示了地-井觀測方式的視電阻率響應(yīng)曲線，可以看出，由于低阻異常體對電流的吸引作用，低阻異常體與圍巖之間的電流場相互耦合，導(dǎo)致視電阻率響應(yīng)在一定區(qū)域內(nèi)出現(xiàn)大于圍巖電阻率的假異常。

圖9c 展示了地-井觀測方式在不同鉆孔孔徑情況下的視電阻率響應(yīng)曲線。相比于井液電阻率，視電阻率響應(yīng)更易受鉆孔孔徑大小影響。且隨著鉆孔孔徑的增大，視電阻率曲線在淺部畸變更為嚴(yán)重。而圖9a 和圖9c 在200 m 以下視電阻率曲線基本不受鉆孔因素影響，在實際電阻率資料解釋中收發(fā)距較小測點的視電阻率響應(yīng)需要加以重視，收發(fā)距較大的測點在資料信噪比較高的情況下可以直接進行處理解釋。而對于電阻率法測井，發(fā)射源與測點均在井中，可以直接借助本文算法確定最佳收發(fā)距，以消除鉆孔參數(shù)的影響，從理論上確保測井?dāng)?shù)據(jù)的有效性。

3 結(jié) 論

a.提出基于四面體和三棱柱混合網(wǎng)格有限元的直流電阻率法三維正演算法。通過一維層狀模型的數(shù)值驗證，混合網(wǎng)格算法可以在保證網(wǎng)格質(zhì)量的同時，提高計算精度。

b.為有效處理鉆孔網(wǎng)格剖分，采用混合網(wǎng)格算法對井-地、地-井方式下的地電模型進行正演模擬。數(shù)值結(jié)果表明，井液電阻率或鉆孔孔徑的變化會對收發(fā)距較小測點處的視電阻率響應(yīng)產(chǎn)生嚴(yán)重的影響，在實際井-地、地-井電阻率法勘探過程中，鉆孔參數(shù)應(yīng)加以重視。

c.對于地-井觀測方式，由于測點位于地下，復(fù)雜的地層信息和鉆孔參數(shù)會導(dǎo)致地-井觀測方式的視電阻率響應(yīng)出現(xiàn)假異常。

d.在實際勘探中，可利用本文算法，結(jié)合鉆孔信息進行正演模擬分析，選擇合適的收發(fā)距，以保證實際直流電阻率資料的合理解釋。為更好地解決以上因素帶來的問題，下一步將展開基于混合網(wǎng)格下的鉆孔參數(shù)約束的直流電法三維反演工作，為直流電阻率法精細化解釋提供理論基礎(chǔ)。