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

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

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

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

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

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

1 2.5D 有限元數(shù)值模擬

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

1.1 穩(wěn)定電流場電位邊值問題

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

當?shù)刭|(zhì)體的電性分布沿某一方向不變時，三維問題可以通過傅里葉變換轉(zhuǎn)化為2.5D 問題[21]。

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

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

1.2 有限單元法

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

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

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

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

1.3 數(shù)值模擬精度分析

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

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

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

2 模型算例與結(jié)果分析

2.1 排列裝置及測線布設(shè)介紹

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

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

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

2.2 水平層狀介質(zhì)模型

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

圖3 水平層狀介質(zhì)模型網(wǎng)格剖分Fig.3 Gridding of a horizontal layered medium model

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

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

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

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

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

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

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

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

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

2.4 水平雙低阻異常體模型

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

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

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

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

2.5 斷層破碎帶模型

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

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

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

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

3 結(jié)果與討論

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

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

4 結(jié)論

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

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

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

符號注釋：

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