基于混合網格有限元的直流電阻率法三維正演研究

2022-06-07 10:34:44王新宇毛玉蓉嚴良俊高文龍

煤田地質與勘探 2022年5期

王新宇，王程，毛玉蓉，嚴良俊，周磊，高文龍

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

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

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

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

1 混合網格有限元三維正演算法

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

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

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

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

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

1.2 混合網格單元分析

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

圖1 四面體和三棱柱的網格節點編號Fig.1 Nodes number of tetrahedral and triangular prism element

則第e個單元內任意一點的異常電位為：

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

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

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

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

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

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

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

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

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

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

當三棱柱側面為無窮遠邊界時，擴展矩陣K2e為：

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

1.3 線性方程組求解

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

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

1.4 視電阻率計算

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

式中：u為每個節點的總電位。

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

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

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

2 混合網格有限元電阻率法數值模擬分析

2.1 混合網格算法精度驗證

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

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

圖4 層狀模型的網格剖分Fig.4 Grid discretization of layered model

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

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

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

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

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

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

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

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

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

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

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

圖9 不同鉆孔參數下視電阻率響應曲線Fig.9 Apparent resistivity response curves of different borehole parameters

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

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