基于2.5D 有限元的高密度電法不同裝置勘探效果研究

趙榮春，呂玉增，張 智，韋柳椰，潘泓序

(桂林理工大學 地球科學學院，廣西 桂林 541004)

高密度電法起源于20 世紀70 年代末期的陣列電法探測思想，是在80 年代從美國和日本發展起來的一種電阻率方法[1-3]。高密度電法兼具電剖面和電測深的特點，憑借工作效率高、反映地電信息豐富、施工簡便等優點，已經被廣泛應用于堤壩勘查、考古、工程物探、水文地質勘探等各個方面[4-6]。

高密度電法基本理論與常規電阻率法相同，然而在諸多方面又體現出了較常規電法的優越性，高密度電阻率法具有多種排列裝置，如：溫納排列、偶極-偶極排列、施倫貝謝爾排列、三極排列等排列裝置，不同裝置因其自身的排列方式不同，其在不同測量環境下表現出的探測效果往往有著很大的差異性。因此，在實際工作環境中針對不同探測對象選擇合適的排列裝置就顯得尤為重要，為探究高密度電法不同排列裝置的探測效果，前人在實踐應用及數值模擬方面都做了諸多研究。在實踐應用中，楊佳鳴[7]、李文忠[8]等對具有明顯層狀結構的斷崖剖面，分別采用溫納、偶極-偶極、微分和施倫貝謝爾等裝置進行對比測量，發現施倫貝謝爾裝置對地層分層有較好的反映；郭清石[9]、孫忠輝[10]等采用數值模擬及模型實驗的方法模擬了高密度電法對溶洞、溶槽的探測，對比發現不同裝置對低阻溶洞都具有很好的分辨能力；王剛等[11]采用不同排列裝置分別對地質剖面及輸水管道進行勘察，對不同探測裝置的分辨特點進行了分析。在電阻率法數值模擬方面，R.C.Fox 等[12]利用有限單元法，模擬了2D 構造下使用偶極-偶極裝置測量時地形對探測效果的影響，結果表明地形會對探測帶來顯著影響；K.M.Befus[13]針對電阻率、極化率在復雜條件下的成像問題進行了研究，并在實際探測中進行了成功應用；歐陽永永等[14]通過建立高阻異常體和隱伏直立低阻體模型，分析認為溫納裝置對于高阻異常體的探測優于其他裝置，而對隱伏直立低阻體的探測溫納和三極裝置的探測效果更好；王智等[15]采用GMSH 進行幾何建模，對研究區域進行非結構化四面體網格剖分，實現了三維井地電阻率法的正演模擬；鄭冰等[16]基于有限差分法對高密度電法進行了正演分析，研究發現當探測水平相鄰低阻異常體時，偶極-偶極、三極、對稱四極裝置對異常體的探測效果優于其他裝置；吳小平等[17]應用非結構化的網格剖分，對具有復雜地形條件下的三維電阻率進行了高精度的數值計算；C.Rücker[18]、T.Günther[19]等采用非結構化網格剖分實現了三維電阻率有限元正反演，提高了對復雜模型的處理能力。

在前人諸多的研究中，基于野外實際環境的研究往往會受限于地質環境本身的單一性，因此，無法對多種地質背景及異常形態進行對比研究，在數值模擬方面則缺乏對高密度電法測量進行較為系統的分析。高密度電法因其廣泛的適用性，其二維和三維的正演問題仍然是電法研究的一個重要方向，實際工作中針對不同探測對象如何選擇合適的測量裝置以獲得更好的探測效果是一個值得探索的問題。因此，筆者從點源電位在三維構造中所滿足的微分方程出發，推導了2.5D 電位所滿足的變分問題，采用Delaunay 三角化算法實現非結構化網格剖分，實現了對高密度電法的正演，在此基礎上結合野外實際應用，設計了生產實踐中常見的地質模型，通過設置不同的電性參數，采用不同的高密度電法排列裝置進行正反演計算，并對其探測的效果進行比較和分析，以期為高密度電法勘探的實際應用提供參考和優化策略。

1 2.5D 有限元數值模擬

有一定走向的構造稱為二維構造，可以用2 個坐標變量描述二維構造的地電特征，但是點電源產生的電場是三維的，因此，點源二維構造的電場其實是三維的，可采用傅里葉變換，將三維電場變換成帶參數的2.5D 問題，進而使用有限單元法求解。

1.1 穩定電流場電位邊值問題

穩定電流場中單點源存在的情況下電導率分布均勻的3D 邊值問題可表示為[20]：

當地質體的電性分布沿某一方向不變時，三維問題可以通過傅里葉變換轉化為2.5D 問題[21]。

將方程(2)應用于方程(1)，可得：

用有限單元法求解方程(3)所示的邊值問題，需先將其變換成對應的變分問題[20]，之后進行數值求解。

1.2 有限單元法

為求解方程(3)所對應的變分問題，本文采用有限單元法(Finite Element Method，FEM)[22]在Matlab 中進行數值求解。有限單元法求解的基本思想是將研究區域離散成一組單元，其中的未知電位V由與節點相連的簡單線性插值函數近似表示，因此，必須創建相應的網格進行計算。本文采用Delaunay 三角化算法實現非結構化網格剖分，并對研究區域進行局部加密[23-24]。經過相應的積分變換，方程(3)可表示為[22]：

將方程(4)所示的單個的有限元矩陣結合在一起創建一個線性方程組可表示為：

解線性方程組(5)可得到各節點的電位。當方程(5)求解完畢后需將波數域的計算結果通過反傅里葉變換式(6)將其轉換到實數域中，以獲得測點的真實電位V。

最優化算法中，一般采用數值計算下式代替式(6)的計算，計算離散波數和加權系數的方法可以在文獻[25]中找到，這里不再重復。在本文的計算中，將使用高斯求積和拉蓋爾積分的組合方法[26]進行求解。

1.3 數值模擬精度分析

為驗證正演程序的正確性，分別采用數值模擬與球形異常體解析式進行計算[27]。

低阻球體模型參數設置為：球心距地表h=5 m，球體半徑ra=2 m，低阻球體的電阻率ρ2為50 Ω·m，圍巖的電阻率ρ1為1 000 Ω·m，采用施倫貝謝爾裝置觀測，取極距為15 m 作剖面，如圖1 所示。分析可見，數值計算結果與解析解基本吻合，解析解與數值解平均相對誤差不超過1.5%，正演計算具有較高的精度。

圖1 數值解與解析解對比Fig.1 Comparison between numerical and analytical solutions

2 模型算例與結果分析

2.1 排列裝置及測線布設介紹

高密度電法勘探主要有溫納、偶極-偶極、施倫貝謝爾、二極、三極等排列裝置，以此為基礎，又演變出了十幾種裝置類型，本文將針對實踐中應用較為廣泛的溫納α、溫納β、施倫貝謝爾、偶極-偶極裝置進行研究，在模型計算中統一設定測線長為400 m，電極距為10 m，共布設40 個電極。在反演計算方面，目前市面上有眾多反演軟件及相關的算法，然而瑞典RES2DINV 反演軟件應用最為普遍，因此，本文將以此對數值結果進行最小二乘法反演計算，在反演計算中，設置阻尼系數最大值為0.1，最小值為0.03，最大迭代次數為5 次，收斂極限為5%，不使用圓滑約束。導出反演數據用Surfer繪制反演剖面圖，正演剖面中以1/4 最大電極距近似表示深度。

本文中各排列裝置電極排列如圖2 所示。在溫納α、溫納β排列中，A、B、M、N電極間距離始終固定相等，隔離系數t在測量中逐漸增大；在施倫貝謝爾測量中，AM=BN，系數t逐漸變大，系數m隨著t的變大進行調整，m的變化范圍為1～5；在偶極-偶極裝置測量中，AB=MN，系數t逐漸變大，系數m隨著t的變大而調整，變化范圍為1～3。

圖2 不同排列裝置電極分布Fig.2 Electrode distributions in different array configurations

2.2 水平層狀介質模型

建立如圖3 所示水平層狀介質模型，用以探究高密度電法對覆蓋層的測量效果及各裝置的垂直分辨率，設計覆蓋層厚度為20 m，電阻率為300 Ω·m，下方基巖電阻率為2 000 Ω·m。

圖3 水平層狀介質模型網格剖分Fig.3 Gridding of a horizontal layered medium model

圖4 為水平層狀介質模型視電阻率正反演剖面。綜合分析可見，無論是正演剖面還是反演剖面，4 種裝置對地層都有較好的分層效果，地層界限與理論模型20 m 的深度基本對應，覆蓋層電阻率大小接近設計值300 Ω·m。對比正反演的成果，可見視電阻率正演剖面巖層界限清晰且沒有其他干擾異常的出現，反觀反演剖面則可見一些干擾信息，如溫納α法中出現了明顯的高阻圈閉異常(圖4b)，溫納β中出現了高阻基巖的斷裂(圖4d)，這些干擾信息都有可能導致在處理實際問題時解釋出現誤判，因此，在地形平坦且地層分層結構明顯時，不經反演的實測數據對于解釋具有重要的參考價值。

圖4 水平層狀介質模型視電阻率正演與反演剖面Fig.4 Forward and inversion profiles of apparent resistivity derived from the horizontal layered medium model

2.3 單個低阻、高阻異常體模型

高密度電法常被用來尋找地下低阻充水溶洞或者高阻地下空腔(如：考古、采空區)，建立如圖5 所示正方形異常體模型，用以模擬低阻、高阻地質體的情況。設計異常體位于測線中央正下方，異常體中心位于水平方向200 m 處，異常體頂部埋深20 m，邊長為20 m，圍巖電阻率為800 Ω·m，異常體為低阻時電阻率為50 Ω·m，為高阻時電阻率為5 000 Ω·m 。

圖5 正方形異常體模型網格剖分Fig.5 Gridding of a square anomaly model

圖6 為低阻模型視電阻率正反演剖面。由反演剖面可見，4 種裝置在異常體中心位置30 m 處都有非常明顯的低阻異常現象，都可以很好地對應異常體的空間位置，相比而言，溫納β、偶極-偶極、施倫貝謝爾裝置其反演效果具有更好的收斂性，低阻異常與實際模型大小及位置都有很好的對應，而溫納α正演剖面(圖6a)則在中心異常的兩側出現了低阻假異常；從正演的成果上分析，溫納β、偶極-偶極(圖6c、圖6e)的低阻異常中心基本與實際模型位置對應，施倫貝謝爾的異常(圖6g)只在水平方向上與異常位置對應很好，而溫納α的視電阻率正演剖面(圖6a)則出現了類似于雙低阻異常的現象，與實際模型存在很大的差距，因此，在生產實踐中對于類似異常形態的解釋要加以區分。

圖6 低阻模型視電阻率正演與反演剖面Fig.6 Forward and inversion profiles of apparent resistivity derived from the low-resistivity anomaly model

圖7 為高阻模型視電阻率正反演剖面圖。由反演剖面可見4 種裝置都具有很好的探測效果，高阻異常中心基本與實際模型對應，能很好地反映出高阻異常所處的空間位置；而從正演的情況分析，4 種裝置的視電阻率正演剖面異常形態都很有各自的特點，相對而言溫納α的視電阻率異常形態與模型更匹配，探測效果更好。

圖7 高阻模型視電阻率正演與反演剖面Fig.7 Forward and inversion profiles of apparent resistivity derived from the high-resistivity anomaly model

2.4 水平雙低阻異常體模型

為探究高密度電法不同排列裝置在水平方向上的分辨率，設置雙低阻模型(圖8)，模擬當研究區域存在2 個位置相隔較近的低阻異常體時的地質結構，設置低阻體為正方形，邊長都為20 m，頂部埋深為20 m，異常體中心在水平方向上的位置分別為170、230 m，兩異常體之間距離為40 m，異常體電阻率為50 Ω·m，圍巖電阻率為800 Ω·m。

圖8 水平雙低阻異常體模型網格剖分Fig.8 Gridding of a horizontal double low-resistivity anomaly model

圖9 為水平雙低阻模型視電阻率正反演剖面。由正演剖面可見，溫納α、溫納β、偶極-偶極(圖9a、圖9c、圖9e)這3 種裝置都顯示為以水平方向200 m 為中心的單個低阻異常形態，無法反映出預設的雙低阻異常體的存在；而施倫貝謝爾的正演剖面(圖9g)則出現了明顯的低阻雙異常現象，且低阻異常中心位置與預設模型對應很好。分析反演剖面可見，4 種裝置都出現了明顯的低阻雙異常形態，尤其是偶極-偶極與施倫貝謝爾裝置(圖9f、圖9h)，低阻異常在空間位置上與預設模型具有很好的對應關系，而溫納α與溫納β(圖9b、圖9d)的低阻異常范圍收斂效果則相對差一些，但仍能反映出低阻異常體所處的空間位置。

圖9 水平雙低阻模型視電阻率正演與反演剖面Fig.9 Forward and inversion profiles of apparent resistivity derived from the horizontal double low-resistivity anomaly model

2.5 斷層破碎帶模型

斷層破碎帶為巖石在構造運動下發生斷裂形成，在漫長的地質作用下，其上部通常被第四系沉積物所覆蓋，而破碎帶中常被巖石碎屑、土壤等所填充，常為地下水滲流通道，因此，在地質找水中破碎帶是重要的目標體。由于破碎帶處地質結構相對不穩定，也是道路施工及隧道掘進過程中所要重點查明的對象，針對破碎帶探測的研究具有非常重要的現實意義。根據破碎帶的特征，設計如圖10 所示斷層破碎帶模型，設定覆蓋層厚度為20 m，電阻率為500 Ω·m，下伏基巖電阻率為2 000 Ω·m，破碎帶電阻率為50 Ω·m，傾斜角為60°，寬度為10 m，其頂部位置位于測線水平方向230～240 m 處。

圖10 斷層破碎帶模型網格剖分Fig.10 Gridding of a fault fracture zone model

圖11 為斷層破碎帶模型視電阻率正反演剖面。分析各裝置正演剖面可見，溫納α與施倫貝謝爾裝置所示剖面(圖11a、圖11g)在破碎帶處整體呈高阻，不能反映出低阻異常的存在，而溫納β與偶極-偶極裝置的剖面圖(圖11c、圖11e)則在破碎帶處呈現出明顯的低阻異常圈閉；在反演剖面中，溫納β與偶極-偶極裝置的剖面(圖11d、圖11f)在破碎帶位置也存在明顯的低阻異常圈閉，而溫納α的反演剖面(圖11b)在破碎帶處呈現出斷開的低阻異常，且低阻的形態與破碎帶的傾斜方向大致相當，施倫貝謝爾裝置的低阻異常形態與溫納α的大致相當，相比之下施倫貝謝爾裝置的低阻異常相對較弱。

圖11 破碎帶模型視電阻率正演與反演剖面Fig.11 Forward and inversion profiles of apparent resistivity derived from the fault fracture zone model

3 結果與討論

高密度電法在實際應用中，地下介質的復雜性和多樣性使得電法數據解釋變得復雜，不同地質體的電性參數差異、地下水體的存在、巖性變化以及構造擾動等因素都會影響數據的解釋。因此，構建不同的地質模型，采用多種采集方式進行分析研究，對于提高高密度電法的分析解釋精度具有重要的作用。筆者采用非結構化網格剖分進行有限元數值模擬，解析解與數值解平均相對誤差不超過1.5%，正演計算具有很高的精度，可以有效地模擬各種類型的地質模型。在此基礎上結合實際應用，設計了生產實踐中常見的5 種地質體模型，并分別采用了溫納α、溫納β、施倫貝謝爾、偶極-偶極4 種排列裝置進行了電阻率的正反演對比研究。

通過對水平層狀介質模型的分析，驗證了高密度電法在探測覆蓋層時的有效性，各裝置均展現出了良好的分層效果，然而反演剖面中顯示了一些干擾信息，如高阻圈閉異常和高阻基巖的斷裂，這些干擾信息可能導致解釋誤判，需要引起注意；綜合單個低阻和高阻模型的情況對比分析可見，溫納β、偶極-偶極、施倫貝謝爾裝置無論是正演剖面還是反演剖面，其電阻率的異常形態基本都是一致的，只有電阻率數值上的差別，而溫納α對于單個低阻體的探測其效果較差，而在探測單個高阻體時又具有很好的探測效果；在對水平雙低阻異常體模型的分析中，發現不同排列裝置在水平方向上的分辨率不同，施倫貝謝爾裝置在探測雙低阻異常體時表現出較好的效果；高密度電法不同裝置在對斷層破碎帶的探測中，不同排列裝置其正反演剖面上展現出了不同的探測效果，尤其是溫納α與施倫貝謝爾的正反演結果其形態上差異較大，因此，對于高密度電法數據的解釋必須結合正反演剖面進行綜合解釋才能得出更準確的結論。

4 結論

a.高密度電法具有多種電極排列裝置，如溫納α、溫納β、施倫貝謝爾、偶極-偶極等。在針對不同的地質目標和測區環境時，常表現出明顯不同的探測效果和分辨能力。因此，在實際勘探工作中，應基于預先調查和人工經驗，針對不同的探測對象和目標層位，選擇多種合適的排列裝置進行測量。對采集的多個排列裝置數據進行綜合分析對比，可獲得更為全面準確的解釋結果。

b.同一排列裝置在正演和反演計算所得的視電阻率異常形態上存在一定差異。正演剖面能更直觀展現目標體的空間位置，但也可能出現一些與實際情況不符的異常，如溫納α對單個低阻立方體及低阻破碎帶的正演效果較差。而反演剖面雖能較好揭示目標體異常特征，但也可能產生一些干擾異常信息，如溫納α出現高阻圈閉異常等。因此，在數據解釋時，需要結合正演和反演2 種結果綜合分析，避免片面性導致的誤判。特定目標體的解釋應優先參考能更好反映其異常特征的計算結果。

c.在探測未知區域單個異常體時，綜合考慮探測精度和工作效率，溫納β、施倫貝謝爾排列具有更佳效果。當預計區域存在多個異常體時，應優先采用具有較好水平分辨能力的施倫貝謝爾和偶極-偶極排列。對于低阻破碎帶等類型異常，溫納α和偶極-偶極排列將獲得更好的探測結果。而在存在明顯分層界限的地層中，上述4 種排列均可體現出較好的分層效果。

符號注釋：

a為電極間距；dA為點電源A到球心D的距離；F為電源項和邊界條件；gi為加權系數；H為一個L×L階的稀疏矩陣；H1i和H2i為僅依賴于節點坐標和單元形狀的有限元矩陣；i為離散波數的個數；I為電源；k為波數；ki為離散波數；K0(kr)和K1(kr)分別為零階和一階貝塞爾函數；L為節點總數；ls和l∞分別為變換后的地面及無窮遠邊界；n為勒讓德函數階數；n′為邊界的外法向向量；N為離散波數的上限；Ne為網格單元的個數；r為測點至點源的距離；rM為球心D到測量極M的距離；RAM為供電電極A到測量電極M的距離；t為隔離系數；為A點電源對M極施加的電位；V為電位；V? 為經傅里葉變換后波數域電位；V?′為包含每個節點的電位；δ(x0)δ(y0)δ(z0)為以點源為中心的Dirac 函數；(x0,y0,z0) 為點源的位置；σ(x,y,z)為電導率分布函數；Ωs為傅氏變換的剖面；Ωt為研究區域；Γs為地面邊界；Γ∞為無窮遠邊界；ρ1為圍巖的電阻率；ρ2為低阻球體的電阻率；ps為視電阻率；μ為ρ2/ρ1的值；Pn(cos∠ADM)為n階勒讓德函數；?為矢量微分算子。