王 鈞，朱國維

(中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院，北京 100083)

1 引 言

高密度電法是當前應(yīng)用較廣的淺層地球物理探測技術(shù)，三維高密度成像數(shù)據(jù)主要是通過布設(shè)的多條二維測線數(shù)據(jù)平移合并或真三維觀測系統(tǒng)采集完成[1]。在實際的地質(zhì)體探測中，關(guān)注較多的是可能存在的異常體空間定位與分布范圍，而不僅僅是一條測線反演出來的某個剖面或多剖面的簡單疊加[2]。當前三維高密度電法的探究使用已經(jīng)較為廣泛，施龍青等將三維高密度電法應(yīng)用于底板與頂板的含水探測中，克服了二維情況下只能獲取巷道底板和頂板電阻率而不能得到地質(zhì)體內(nèi)部電阻率分布特征的缺點[3,4]。黃俊革等使用E-SCAN電極布設(shè)方式對常見的地質(zhì)體進行模擬，得到了低阻異常響應(yīng)較好的結(jié)果，建議先用此方式測量區(qū)域的面積，以此確定可能存在的異常體范圍，再反演出縱向的位置特征[5]。黃真萍等將三維高密度電法用于各種地電模型，提出了反演體積效應(yīng)的假象問題，認為這種方法對于異常體有較好的分布探測能力，但要綜合考慮存在的干擾因素[6]。高陽等在巖溶塌陷區(qū)使用三維高密度電法探測技術(shù)，獲得了巖溶地質(zhì)體立體分布的尺寸范圍，表明此方法對于常規(guī)電法存在的體積與旁側(cè)效應(yīng)有一定的減弱能力[7]。朱瑞等在擬建水庫壩區(qū)開展三維高密度電法地質(zhì)情況探測，并結(jié)合現(xiàn)場鉆孔資料驗證方法的可靠性，提出了為避免多解性問題，應(yīng)采取綜合探測分析方法[8]。張先林等通過時移三維高密度電法對黃土層灌溉水入滲方式進行研究，并通過Voxler軟件將反演結(jié)果進行成像處理，得到了試驗區(qū)的入滲過程與對應(yīng)規(guī)律[9]。本文將三維高密度電法正演軟件Res3dmod的模擬結(jié)果導(dǎo)入反演軟件Res3dinv中進行計算，然后結(jié)合Matlab程序?qū)崿F(xiàn)了可視化反演成像，并分析了成像效果，得到了各因素變化帶來的響應(yīng)特征及規(guī)律，最后根據(jù)實例的Voxler軟件成像結(jié)果討論了所對應(yīng)的實際地質(zhì)情況。

2 方法原理

高密度電法基于可能潛在的目標地質(zhì)體與周圍巖土層的電阻率差異來進行探測分析，通過觀察人工建立的地下穩(wěn)定電流場的傳導(dǎo)分布規(guī)律來獲取對應(yīng)地區(qū)的地質(zhì)情況(圖1)。在三維笛卡爾坐標系下，假設(shè)任一點(x0,y0,z0)處電流為+I，則有[10]：

(1)

=?·( [σ]? [U] )

(2)

圖1 電阻率法電極連接Fig.1 The electrode connection of resistivity method

(3)

ρ為電阻率(Ω·m);K為裝置系數(shù);A、B為供電電極;M、N為測量電極。由于實際工作需要，選擇A、B、M、N的相對位置來設(shè)置探測裝置，現(xiàn)場探測得到的電阻率受多因素的影響，因此其為視電阻率[14]，例如圖2的單極-偶極裝置[15]：

圖2 單極-偶極裝置示意圖Fig.2 The schematic diagram of pole-dipole device

此時，測量的視電阻率

(4)

其具有較高的橫向分辨率，采用單極供電方式，結(jié)果與其他裝置間具有轉(zhuǎn)換性[16]。

2.1 有限差分法

在實際計算中，可以將電場的總電位U總表征為正常場U正與異常場U異之和[17，18]，正常場可根據(jù)解析方程式求解：

(5)

異常場通過有限差分法迭代求解，有限差分主要依據(jù)泰勒級數(shù)展開等數(shù)學(xué)近似手段，以點差商代替點微商，其可以采用向前、向后，以及中心差分格式。以中心差分交錯網(wǎng)格為例，2M階精度為：

(6)

式(6)中，Cm為對應(yīng)差分系數(shù)，將此格式代入式(2)中即可得差分控制方程?？梢詫⑷拷Y(jié)點建立的差分控制方程寫成矩陣形式[19]：

[C]·[U]=[Q]

(7)

式(7)中，[C]是差分方程的系數(shù)矩陣，與地層電阻率相關(guān)；[U]是節(jié)點的電位矩陣；[Q]為給定地質(zhì)體分布及邊界后的常向量。計算式(7)即可得到電位的空間點值。下文的模型一至九采用此方法。

2.2 有限元法

有限元是運用變分的方法將微分方程邊值求解換成泛函極值求解，進而將求解區(qū)域劃分為各個子單元，建立插值函數(shù)與數(shù)值積分形式進行單元的集成迭代求解(圖3)。其能夠處理不規(guī)則模型，如下文的起伏地形模型十和十一。

圖3 網(wǎng)格分布示意圖Fig.3 The schematic diagram of grid distribution

此方法針對的變分問題可以表述為式(8)[20]：

式中，Ω表示區(qū)域;Γ為邊界向量;n指向邊界外法線;r為電源點到邊界點的距離。

2.3 光滑約束最小二乘法

在三維電阻率數(shù)據(jù)的反演過程中，常采用光滑約束的最小二乘優(yōu)化方法，能夠有效降低模型與實測數(shù)據(jù)之間的差異。三維光滑約束最小二乘法基于以下方程[21]:

(JTJ+λF)Δqk=JTg-λFqk

(10)

(11)

式(10)和式(11)中，Cx、Cy是水平光滑濾波系數(shù)向量，Cz為垂直光滑濾波系數(shù)向量；J是雅克比偏微分向量，JT為J的轉(zhuǎn)置向量；λ為阻尼因子；q為模型擾動向量；g為偏差向量。

3 模型與算例

以單極-偶極裝置為例，X方向電極數(shù)為39個、Y方向電極數(shù)為11個，各方向電極距為1 m，即X方向探測區(qū)間為0～38 m、Y方向探測區(qū)間為0～10 m。采用的模型電阻率值分別有低阻值20Ω·m、背景值100Ω·m，以及高阻值1 000 Ω·m。正演建模時網(wǎng)格劃分為9層：第1層(0～0.47 m)；第2層(0.47～1.09 m)；第3層(1.09～1.84 m)；第4層(1.84～2.66 m)；第5層(2.66～3.47 m)；第6層(3.47～4.68 m)；第7層(4.68～6.08 m)；第8層(6.08～7.49 m)；第9層(7.49～9.01 m)。帶地形反演的結(jié)果模型有10層，其余模型反演結(jié)果均為13層，最大探測擬深為13.2 m，取8個不同深度面進行三維可視化成像。

如圖4所示，模型一異常體設(shè)定在3、4、5層；模型二異常體設(shè)定在4、5、6層；模型三異常體設(shè)定在5、6、7層，平面尺寸均為3 m×4 m?？梢钥吹侥Ｐ鸵浑娮杪手行闹瞪疃任恢脧?.97～3.03 m，模型二電阻率中心值深度位置從1.56～3.94 m，模型三電阻率中心值深度位置從2.25～6.18 m，反演中心位置出現(xiàn)一定偏差，但延伸范圍基本能說明異常體的存在影響區(qū)間。其中處于較低深度的異常體反演形狀邊界清晰，而處于較深位置時異常體反演形狀邊界發(fā)生變形，由于高低阻間的響應(yīng)相互影響，出現(xiàn)了一定的電阻率過渡區(qū)域，且隨著深度的增加，體積效應(yīng)越發(fā)明顯，說明了高密度電法淺層的適用性，而在深層需要注意體積發(fā)散的假象。

如圖5所示，模型四與模型五異常體均設(shè)定在3、4、5層，其中模型四低阻體平面尺寸為2 m×2 m，高阻體平面尺寸為2 m×4 m；模型五低阻體平面尺寸為4 m×4 m，高阻體平面尺寸為4 m×8 m?？梢钥吹?，對于淺層同一深度不同尺寸異常體響應(yīng)的效果一致，區(qū)別在于擴散假象與高阻體尺寸大小呈正相關(guān)，而與低阻體尺寸大小呈負相關(guān)。方形的高阻體對淺層與深層的擴散影響較多。

圖5 模型四、五尺寸大小反演結(jié)果Fig.5 The size inversion results of models 4 and 5

如圖6所示，模型六和模型七的異常體電阻率值與模型四和模型五相反，擴散假象的認識與上述分析一致，方形的低阻體尺寸擴大后對淺層與深層擴散影響變得較多，綜合模型四、五、六、七可見，異常體X向與Y向的擴展比例也是成像的一個影響因素，體現(xiàn)出形體的延伸響應(yīng)特征，尺寸變大后等距異常體產(chǎn)生了更大的影響，反映了異常體不同邊距比的抗干擾能力不同。

圖6 模型六、七高低阻變化反演結(jié)果Fig.6 The inversion results of high and low resistance changes in models 6 and 7

模型八和模型九正演模型的Z向切面分別是梯形與偏轉(zhuǎn)45°的正方形，目的在于研究三維高密度電法對于邊棱的成像能力(圖7)?？梢钥闯?，隨著深度的變大，邊棱逐漸圓滑化。

圖7 模型八、九邊界形狀反演結(jié)果Fig.7 The boundary shape inversion results of models 8 and 9

模型十和模型十一分別為同一標準上引入凸凹地形因素(圖8)?？梢姷匦蔚钠鸱兓瘯σ曤娮杪史植籍a(chǎn)生影響，地表凸地形頂點區(qū)域呈現(xiàn)出較低阻特征，而在兩邊小范圍內(nèi)出現(xiàn)較高阻，地表凹地形得到相反的結(jié)果，原因在于等電位面在凸處變寬，電流線變稀疏，而等電位面在凹處變窄，電流線變密集。而且地形效應(yīng)將嚴重影響局部電阻率分布，使之呈現(xiàn)一定的“起伏狀”，體積邊界效應(yīng)也更加明顯。

圖8 模型十、十一地形凸凹反演結(jié)果Fig.8 The terrain convex and concave inversion results of models 10 and 11

此外，模型一反演最小單元體積為0.45個單位，最大為625.694個單位，平均靈敏度為9.603 1，反演結(jié)果模型有13層，誤差為0.23 %；將模型一X向極距增大一倍，Y向極距不變，同一標準反演結(jié)果擬深也增大一倍，最小單元體積為0.9個單位，最大為2 188.685個單位，平均靈敏度為6.956 4，誤差為0.7 %，反演模型變?yōu)槭邔印＿@說明為了提高勘探深度，可以適當增大測線長度，但若不增加電極，分辨率將會下降。在模型一正演的結(jié)果上加入5 %噪聲后，反演的誤差為3.93 %，這說明了誤差并不是越小越好，因為正反演以均質(zhì)為前提模擬異常，這并不符合實際復(fù)雜地質(zhì)條件，而且在進行迭代計算時，應(yīng)當依據(jù)的是后兩次的變化幅度，這符合計算的穩(wěn)定變化特征。

應(yīng)用赤道-偶偶極、單極-偶極、偶極-偶極、溫納-施倫貝爾，以及溫納裝置對同一模型進行正反演，得到的數(shù)據(jù)量從大到小為上述裝置從左至右。探測深度從大到小分別為單極-偶極、赤道-偶偶極、偶極-偶極、溫納-施倫貝爾和溫納裝置。各裝置反映的異常體特征相近，但響應(yīng)的深淺水平位置及邊界形狀的分辨率有一定差異。

4 實例探測

以大柳塔煤礦52 304工作面附近地裂縫探測為例，本次試驗使用E60D高密度電法儀和溫納裝置，僅采用了49個電極，電極距為10 m。三維高密度野外工作布線采用圖9模式。

圖9 三維測線布設(shè)Fig.9 The layout of 3D survey line

由圖10可見，探測區(qū)地裂縫在地表淺部顯現(xiàn)出高阻特性，有一定的布展和延伸趨勢，與現(xiàn)場地面狀況相符合；此外存在的低阻區(qū)域可能與地層賦水情況相關(guān)。同時可看到電法固有的體積邊界效應(yīng)，呈現(xiàn)“發(fā)散”狀態(tài)對于異常體位置及范圍的圈定是不利的，容易造成“多解性”。

圖10 地裂縫反演結(jié)果Fig.10 The inversion results of ground fissure

5 結(jié) 論

通過三維高密度電法正反演數(shù)值模擬，得到了對應(yīng)異常體的響應(yīng)規(guī)律，注意到高密度電法存在的體積邊界效應(yīng)，受到諸如測線電極距、地質(zhì)體分布、裝置形式、網(wǎng)格剖分，邊界條件以及噪聲阻尼等參數(shù)的疊加影響。因此探測反演結(jié)果以平面來看為梯形范圍等形式截取可能是比較符合實際的響應(yīng)范圍。在進行高密度電法探測時，需因地制宜選擇合適裝置，若測線布設(shè)長度及電極密度達到要求，三維高密度電法成像后能夠以一定分辨率立體展現(xiàn)地質(zhì)異常體的延展范圍與特征，因此該方法在工程與環(huán)境地球物理領(lǐng)域是一種有效的探測方法。此外，為避免電法自身存在的缺陷問題以及復(fù)雜地質(zhì)環(huán)境的交互影響，多方法綜合探測或聯(lián)合反演是解決三維高密度電法多解不收斂假象的途徑。