梁 卓， 熊 彬， 徐志峰， 蘭懷慷， 羅天涯， 李祖強(qiáng), 崔楊洋

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

0 引言

天然形成的洞穴(喀斯特溶洞、陷落洞等)與人工洞穴像地道、采空區(qū)等對地面建筑可能造成嚴(yán)重危害，所以查明地下空洞的位置和形態(tài)十分必要。電阻率法已廣泛應(yīng)用于洞體探測，并獲得了不同程度的效果。

在電阻率法研究中，由于無法得到復(fù)雜地電條件下電場的解析解，需要使用數(shù)值模擬方法求解。有限差分法、有限單元法等需要進(jìn)行網(wǎng)格剖分，并求解大型線性方程組，計(jì)算量較大，往往還需要借助非結(jié)構(gòu)化網(wǎng)格，而邊界單元法[1-6]可以較好地解決地形二、三維正演問題[7-11]，不僅適合求解有界區(qū)域的問題，而更適合無界區(qū)域問題。筆者使用邊界單元方法快速計(jì)算點(diǎn)電流源在二度體地電模型上激發(fā)的電場，獲得二度體條件下的復(fù)雜地電條件的異常特征，為復(fù)雜地電異常的解釋提供借鑒。

由于電位滿足的積分方程的核是一個(gè)以為奇點(diǎn)的函數(shù)，在使用邊界單元法離散該積分方程時(shí)，在包含點(diǎn)電流源的單元上，由于電位的梯度變化很大，運(yùn)用較小的離散單元尺度法，同時(shí)采用二次場近似方法。在相同單元尺度條件下，我們比較了用總場法和二次場近似法求解場源單元內(nèi)的電位；計(jì)算結(jié)果表明，在適當(dāng)?shù)臈l件下，兩種方法都能達(dá)到要求的精度，在相同的網(wǎng)格條件下，對于不同的地電模型各有優(yōu)勢。

1 二度體的電位控制方程

圖1 2.5D邊值問題模型示意圖Fig.1 2.5 Dimension boundary value problem

如圖1所示，二度體沿z軸方向延伸，大地是電阻率為ρ1的均勻各向同性介質(zhì)，地表起伏且在地表A點(diǎn)放置電流強(qiáng)度為I的供電電極，地表邊界Γ1以及由封閉曲線Γ2共同構(gòu)成區(qū)域Ω2，在區(qū)域Ω1中，有由封閉的曲線構(gòu)成的區(qū)域Ω1，區(qū)域Ω1的介質(zhì)的電阻率為，則電位滿足的偏微分方程的邊值問題[12]如下：

在區(qū)域Ω1中：2φ1=0

(1)

在區(qū)域Ω2中：2φ2=0

(2)

(3)

(4)

(5)

其中δs(A)為面Dirac函數(shù)。

考慮到地電結(jié)構(gòu)為二度體，點(diǎn)源產(chǎn)生的電位關(guān)于垂直x軸的過點(diǎn)電流源的位置的平面是對稱的，即φ(x,y,z)=φ(-x,y,z)，則對電位的x變量做余弦變換，即：

(6)

(7)

式中:ω1p、ω2p是點(diǎn)P對區(qū)域Ω1、Ω2所張立體角，r是點(diǎn)p到點(diǎn)電流源的距離。從方程(7)可知，當(dāng)異常體的電阻率和背景電阻率相等時(shí)，方程(7)為起伏界面時(shí)電位的積分方程。

2 積分方程的離散

積分方程(7)直接求解非常困難，為此采用邊界單元法求解該積分方程。將Γ1∪Γ2剖分成n個(gè)單元，設(shè)單元節(jié)點(diǎn)i上的電位為Ui，為了避免點(diǎn)源在節(jié)點(diǎn)上的奇異性，將點(diǎn)電流源放在某個(gè)單元的內(nèi)部。

1)算法1。設(shè)每個(gè)單元上的總電位U是線性變化的，在某個(gè)單元內(nèi)對總電位采用線性插值，對節(jié)點(diǎn)i有式(8)的關(guān)于離散總電位的方程。

(8)

對每個(gè)節(jié)點(diǎn)均可得到上述方程，由全部節(jié)點(diǎn)得到含有n個(gè)方程n個(gè)未知數(shù)的方程組為式(9)。

(Θ-F)U1=B1

(9)

2)算法2。由于在點(diǎn)電源附近電位的梯度比較大，在一個(gè)單元內(nèi)的電位不能使用線性變化近似，為了考慮電位在點(diǎn)電源單元內(nèi)的非線性特征，將電位看成是點(diǎn)電流源在均勻空間產(chǎn)生的一次場電位和在一個(gè)單元內(nèi)為常數(shù)的二次電位C之和。

U=φ0+C

(10)

在某個(gè)節(jié)點(diǎn)i上，由式(10)得到：Ci=Ui-φ0i，則在該單元上的電位可近似表示為式(11)。

U=Ui+φ0-φ0i

(11)

將式(11)代入式(7)，得到節(jié)點(diǎn)上電位的方程為式(12)。

(12)

對每個(gè)節(jié)點(diǎn)都可以得到方程(12)，并將各離散量用列矢量表示，則得到關(guān)于電位的線性方程為式(13)。

(Θ-H)U1=B1+C1

(13)

由式(6)可知空間域的電位為式(14)。

(14)

在過點(diǎn)電源的主剖面上,x=0，得到式(15)。

(15)

上述廣義積分可由如下的6個(gè)積分系數(shù)辛普森積分公式算得：

(16)

式中km、Cm值列于表1中。

表1 積分系數(shù)

圖2 算法結(jié)果對比情況Fig.2 Comparisons of algorithm results(a)結(jié)果與精確解對比;(b)誤差曲線

3 均勻半空間結(jié)果比較

對場源單元內(nèi)電位采用總場法和二次場近似法求解的兩種情況(簡稱算法1和算法2)。假設(shè)地表水平，單元大小為1 m，大地電阻率為1 Ω·m，地表各節(jié)點(diǎn)的電位如圖2所示。從圖2中可以看出，對場源單元內(nèi)電位作總場法和二次場近似法的不同結(jié)果與精確解作對比，三者吻合較好，均驗(yàn)證了算法的正確性與可靠性。在電源點(diǎn)附近r=5 m的范圍以內(nèi)，采用總場法的算法1中相對誤差在0.005%～0.292%，而作二次場近似法的算法2中相對誤差在0.001%～0.288%，算法2精確性明顯高于算法1；r>5 m時(shí)，算法2的誤差顯著增大，并且明顯大于算法1。經(jīng)上述11個(gè)點(diǎn) 的統(tǒng)計(jì)，算法1中邊界單元法的計(jì)算誤差在0.005%～0.292%，算法2中誤差在0.001%～0.338%。這表明，對單元內(nèi)的電位作二次場近似法可以明顯提高電源點(diǎn)附近計(jì)算結(jié)果的精度。

圖3 模型示意圖 1Fig.3 Model diagram 1

圖4 視電阻率曲線及相對誤差曲線Fig.4 Apparent resistivity curves and relative error curves(a)ρ=1 Ω·m; (b)ρ=10 Ω·m;(c)ρ=1 000 Ω·m; (d)ρ=1 000 000 Ω·m

4 存在異常體時(shí)的計(jì)算結(jié)果對比

將源點(diǎn)A固定在原點(diǎn)，采用單極-偶極裝置(pole-Dipole)，大地電阻率為100 Ω·m，電流強(qiáng)度大小為1.0 A，在水平距離52.5 m下方分別賦存不同電阻率值、不同埋深和不同半徑大小的長圓柱目標(biāo)體，測點(diǎn)由左向右移動(dòng)，并計(jì)算兩種算法所得的視電阻率數(shù)據(jù)。

4.1 模型一

在水平距離52.5 m下方賦存不同電阻率值的長圓柱目標(biāo)體，半徑r=1.5 m，中心埋深h=20 m。電阻率值分別為1 Ω·m、10 Ω·m、1 000 Ω·m、1 000 000 Ω·m(按空洞模型假設(shè))。假設(shè)地表水平，單元大小為1 m，取電位電極MN=5 m，測點(diǎn)沿水平軸由左向右移動(dòng)。

以1 Ω·m(圖4(a))、10 Ω·m(圖4(b))、1 000 Ω·m(圖4(c))、1 000 000 Ω·m(圖4(d))不同電阻率值的地電模型為例，計(jì)算其視電阻率的數(shù)值解，并分析兩種算法和解析解的相對誤差，結(jié)果如圖4

圖5 模型示意圖2Fig.5 Model diagram 2

所示。從圖4可以看出，高阻異常體的視電阻率異常是一個(gè)不對稱的凸形曲線，低阻異常體的視電阻率異常是一個(gè)不對稱的凹形曲線。算法2與算法1相比，異常峰值點(diǎn)向右偏移。低阻異常明顯大于高阻異常，對于低阻來說，電阻率值越小，兩種算法的相對誤差越小，表明了電阻率法對低阻異常更加敏感。對于高阻來說，隨著高阻電阻率值的增大，兩種算法結(jié)果趨于一致。

4.2 模型二

在模型一的基礎(chǔ)上，中心埋深h分別為20 m、30 m、40 m，異常體電阻率為1 000 000 Ω·m，半徑r=1.5 m。

以20 m(圖6(a))、30 m(圖6(b))、40 m(圖6(c))不同埋深地電模型為例，計(jì)算不同深度下異常體視電阻率的數(shù)值解，并分析兩種算法與解析解對比的相對誤差，結(jié)果見圖6。從圖6可以看出，算法2相對誤差曲線均基本保持不變，埋深較淺時(shí)，算法1視電阻率曲線和相對誤差曲線相對算法2向左偏移，隨著埋深不斷加大，偏移現(xiàn)象逐漸消失，二者視電阻率曲線和相對誤差曲線接近重合。說明二次場近似法在隨著異常體埋藏深度變化時(shí)，計(jì)算結(jié)果比較穩(wěn)定。

圖6 視電阻率曲線及相對誤差曲線Fig.6 Appsarent resistivity curves and relative error curves(a)中心埋深h=20 m; (b)中心埋深h=30 m;(c)中心埋深h=40 m

4.3 模型三

在模型一的基礎(chǔ)上，中心埋深h=20 m，異常體半徑大小分別為1 m、2 m和3 m，異常體電阻率為1 000 000 Ω·m。

圖7 模型示意圖3Fig.7 Model diagram 3

圖8 視電阻率曲線及相對誤差曲線Fig.8 Appsarent resistivity curves and relative error curves(a)r0=1 m；(b)r0=2 m;(c)r0=3 m

以1 m(圖8(a))、2 m(圖8(b))、3 m(圖8(c))不同半徑大小的異常體地電模型為例，計(jì)算不同大小的異常體視電阻率的數(shù)值解，并分析兩種算法下計(jì)算結(jié)果與解析解對比的相對誤差，結(jié)果如圖8所示。

從圖8可以看出，兩種算法相對誤差均低于1.5%，與解析解擬合程度較好。隨著異常體半徑不斷加大，算法1和算法2的正演結(jié)果有所差異。當(dāng)r0=1 m時(shí)，算法1和算法2的計(jì)算效果接近一致，說明當(dāng)異常體較小時(shí)，對電源點(diǎn)采用總場法和二次場近似法沒有太大區(qū)別。隨著異常體半徑不斷加大，總場法(算法1)相對誤差曲線逐漸向左偏移，且精度有所提高。因此對單元內(nèi)的電位作總場法更適合計(jì)算半徑稍大的地電模型。

5 結(jié)論

我們對單元內(nèi)電位采用總場法和二次場近似法進(jìn)行對比研究，表明這兩種方法都能準(zhǔn)確計(jì)算視電阻率值，異常曲線形態(tài)相似，異常幅值不同，均能滿足精度要求且各有優(yōu)勢。在不同地電模型的正演計(jì)算中存在一定差異。當(dāng)異常體大小及埋深一定時(shí)，對于不同電阻率兩種算法具有較高精度；隨著埋深增大，二次場近似法計(jì)算結(jié)果比較穩(wěn)定；隨著半徑增大，總場法精度更高。實(shí)驗(yàn)結(jié)果為探測地下洞體提供了借鑒,對于復(fù)雜地電模型還值得進(jìn)一步深入的研究。