王 鈞，朱國(guó)維

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

1 引 言

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

2 方法原理

高密度電法基于可能潛在的目標(biāo)地質(zhì)體與周?chē)鷰r土層的電阻率差異來(lái)進(jìn)行探測(cè)分析，通過(guò)觀察人工建立的地下穩(wěn)定電流場(chǎng)的傳導(dǎo)分布規(guī)律來(lái)獲取對(duì)應(yīng)地區(qū)的地質(zhì)情況(圖1)。在三維笛卡爾坐標(biāo)系下，假設(shè)任一點(diǎn)(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為測(cè)量電極。由于實(shí)際工作需要，選擇A、B、M、N的相對(duì)位置來(lái)設(shè)置探測(cè)裝置，現(xiàn)場(chǎng)探測(cè)得到的電阻率受多因素的影響，因此其為視電阻率[14]，例如圖2的單極-偶極裝置[15]：

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

此時(shí)，測(cè)量的視電阻率

(4)

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

2.1 有限差分法

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

(5)

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

(6)

式(6)中，Cm為對(duì)應(yīng)差分系數(shù)，將此格式代入式(2)中即可得差分控制方程。可以將全部結(jié)點(diǎn)建立的差分控制方程寫(xiě)成矩陣形式[19]：

[C]·[U]=[Q]

(7)

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

2.2 有限元法

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

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

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

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

2.3 光滑約束最小二乘法

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

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

(10)

(11)

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

3 模型與算例

以單極-偶極裝置為例，X方向電極數(shù)為39個(gè)、Y方向電極數(shù)為11個(gè)，各方向電極距為1 m，即X方向探測(cè)區(qū)間為0～38 m、Y方向探測(cè)區(qū)間為0～10 m。采用的模型電阻率值分別有低阻值20Ω·m、背景值100Ω·m，以及高阻值1 000 Ω·m。正演建模時(shí)網(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層，最大探測(cè)擬深為13.2 m，取8個(gè)不同深度面進(jìn)行三維可視化成像。

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

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

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

如圖6所示，模型六和模型七的異常體電阻率值與模型四和模型五相反，擴(kuò)散假象的認(rèn)識(shí)與上述分析一致，方形的低阻體尺寸擴(kuò)大后對(duì)淺層與深層擴(kuò)散影響變得較多，綜合模型四、五、六、七可見(jiàn)，異常體X向與Y向的擴(kuò)展比例也是成像的一個(gè)影響因素，體現(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°的正方形，目的在于研究三維高密度電法對(duì)于邊棱的成像能力(圖7)?？梢钥闯?，隨著深度的變大，邊棱逐漸圓滑化。

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

模型十和模型十一分別為同一標(biāo)準(zhǔn)上引入凸凹地形因素(圖8)。可見(jiàn)地形的起伏變化會(huì)對(duì)視電阻率分布產(chǎn)生影響，地表凸地形頂點(diǎn)區(qū)域呈現(xiàn)出較低阻特征，而在兩邊小范圍內(nèi)出現(xiàn)較高阻，地表凹地形得到相反的結(jié)果，原因在于等電位面在凸處變寬，電流線變稀疏，而等電位面在凹處變窄，電流線變密集。而且地形效應(yīng)將嚴(yán)重影響局部電阻率分布，使之呈現(xiàn)一定的“起伏狀”，體積邊界效應(yīng)也更加明顯。

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

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

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

4 實(shí)例探測(cè)

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

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

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

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

5 結(jié) 論

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