付德俊

(廣東核力工程勘察院， 廣州 510800)

0 引言

直流電阻率法(DC)是利用地下介質(zhì)的電阻率參數(shù)來實(shí)現(xiàn)對地下地質(zhì)結(jié)構(gòu)的分辨和研究，該方法被廣泛應(yīng)用到環(huán)境與工程勘察中，其已成為目前常用的物探方法之一。

常規(guī)直流電阻率法將觀測電極布設(shè)在地表，來觀測地表電位的變化，并采用特定觀測裝置來刻畫出需要的參數(shù)(如視電阻率等)，顯示出地下地質(zhì)體在水平或垂直方向存在電性差異，從而實(shí)現(xiàn)對地質(zhì)體的定性分析。但是，對于特殊地質(zhì)情況，地表觀測到電位或視電阻率數(shù)據(jù)很難對深部地質(zhì)體進(jìn)行約束，較難地刻畫地下地質(zhì)體的邊界范圍，無法滿足中深部工程勘察需求，為了彌補(bǔ)這一缺陷，開展井地直流電阻率法反演成像研究具有重要的研究意義，其彌補(bǔ)了縱向上深部數(shù)據(jù)對地下地質(zhì)體約束。

井地直流電阻率常用的觀測裝置可以分為井-地、地-井-地、井-井以及井-地-井等，為了研究不同觀測裝置對地下地質(zhì)體的分辨能力，具有重要意義。目前，國內(nèi)、外研究學(xué)者進(jìn)行很多相關(guān)的研究工作，Zhou等[1-2]采用解析偏導(dǎo)數(shù)矩陣完成了2D/3D井間DC快速成像研究；Wikinson等[3]重點(diǎn)分析了礦井DC反演成像技術(shù)的研究工作；國內(nèi)學(xué)者同樣在這方面做了大量的研究工作；王志剛等[4-6]實(shí)現(xiàn)積分方程法的Born 近似井地DC三維快速成像研究；呂玉增等[7]實(shí)現(xiàn)了井間三維直接成像研究；岳建華等[8]開展了巷道井地DC快速成像研究；徐凱軍等[9]開展了基于共軛梯度算法的三維井地電阻率成像研究；藍(lán)澤鸞等[10]開展了基于TV約束的2.5維井地直流電阻率反演研究。

筆者采用變分原理推導(dǎo)2.5維井地DC滿足邊值問題，構(gòu)建了基于二階最大平滑因子的井地DC正則化反演目標(biāo)函數(shù)，采用共軛梯度算法實(shí)現(xiàn)了該目標(biāo)函數(shù)最優(yōu)化求解，最終完成不同觀測方式的井地2.5維直流電阻率穩(wěn)定、高效的反演研究。研究結(jié)果表明，不同觀測方式的反演對異常體的約束存在明顯的差異，靠近井所在的位置對異常體的約束較為明顯，但是若沒有地表數(shù)據(jù)的約束可能出現(xiàn)對異常體過渡約束，使異常體的范圍出現(xiàn)明顯的偏差。

1 井地2.5維直流電阻率正演

1.1 邊值問題

在波數(shù)域，井地2.5維DC滿足的邊值問題[11]如式(1)所示。

(1)

式中:區(qū)域Ω為二維研究區(qū)域；Γs、?！奘嵌S區(qū)域的邊界，Γs為區(qū)域Ω的地表邊界，?！逓閰^(qū)域Ω的地下邊界；σ為電導(dǎo)率；r為點(diǎn)電源到邊界的距離；n為無窮遠(yuǎn)邊界的外法向方向；k為波數(shù)；K0、K1分別為零階、一階第二類貝塞爾函數(shù)。

將公式(1)轉(zhuǎn)化為變分問題：

Iδ(A)U]dΩ+

(2)

對計(jì)算區(qū)離散、構(gòu)造插值函數(shù)，式(2)離散為：

(3)

其中，P為除供電點(diǎn)為0.5以外全為零的列向量。

對式(3)求變分并令其為零，得到線性方程組

KU=P

(4)

通過求解線性式(4)，得到波數(shù)域的電位U，然后通過傅氏反變換得到電位：

(5)

式中:r為點(diǎn)的位置；km是波數(shù)；gm是加權(quán)系數(shù)[11]。

表1 傅氏變換系數(shù)km和gm

1.2 2.5維井地DC觀測方式

常用的井地DC觀測裝置可以分解成一下四種：①井-地；②地-井-地；③井-地-井；④井-井測量裝置，其野外布置方式圖1所示。

圖1 2.5維井地DC觀測裝置示意圖Fig.1 Schematic diagram of2.5-D well ground DC observation device(a)井-地;(b)地-井-地; (c)井-地-井;(d)井-井

2 共軛梯度正則化反演

構(gòu)建正則化反演目標(biāo)函數(shù)[12]：

φ(m,α)=φd(m)+αφm(m)

(6)

其中：φ(m,α)為構(gòu)建的正則化目標(biāo)函數(shù)；α是正則化因子；φd(m)為數(shù)據(jù)約束項(xiàng)；φm(m)為模型約束項(xiàng)；m為模型參數(shù)。數(shù)據(jù)約束項(xiàng)與模型約束項(xiàng)表達(dá)式分別為：

φd(m)=Wd(A(m)-dobs)TWd(A(m)-dobs)

φm(m)= (m-mapr)T(m-mapr)

(7)

其中：A表示為正演算子；dobs表示為觀測數(shù)據(jù)；mapr為已知地下結(jié)構(gòu)的先驗(yàn)信息；Wd為數(shù)據(jù)權(quán)重矩陣。為了獲得最優(yōu)的解，必須對正則化目標(biāo)函數(shù)求最小，即

φ(m,α)=φd(m)+αφm(m)→min

(8)

筆者采用共軛梯度算法求解公式的最優(yōu)化問題，其采用的基本思想為最速下降法[13]：

mn+1=mn+Δm

(9)

初始梯度迭代向量表達(dá)式為：

(10)

首次迭代方向?yàn)樘荻认蛄颗c首次迭代共軛梯度向量的線性組合

(11)

第n+1次迭代的方法為：

(12)

采用線性搜索求解目標(biāo)函數(shù)最小來獲取最速下降迭代步長：

(13)

通過簡化，可獲取到下降步長為式(14)。

(14)

3 算例分析

3.1 算法驗(yàn)證

圖2展示了電導(dǎo)率為0.01 s/m均勻半空間的理論電位曲線與數(shù)值模擬的電位曲線對比圖，相對誤差范圍在2%以內(nèi),對比結(jié)果表明本文正演算法正確可靠。

圖2 數(shù)值解與解析解電位對比曲線圖Fig.2 Comparision of FEM numerical and one- dimensional analytical solutions

圖3 地電模型斷面圖Fig.3 Sectional view of the electric model

圖4 井-地反演斷面圖Fig.4 Sketch map of borehole-surface dc resistivity

圖5 地-井-地反演斷面圖Fig.5 Sketch map of surface-borehole- surface dc resistivity

圖6 井-地-井反演斷面圖Fig.6 Sketch map of well - to - well dc resistivity

圖7 井-井反演斷面圖Fig.7 Sketch map ofcrosshole dc resistivity

3.2 不同觀測方式反演對比分析

圖3為地電模型，均勻半空間下存在埋深以及大小相同而不同電導(dǎo)率的兩個(gè)異常體，異常體1的電導(dǎo)率為0.005 s/m ，異常體2的電導(dǎo)率為0.02 s/m，背景電導(dǎo)率為0.01 s/m，異常體尺寸為4 m×4 m，異常體的頂部埋深為2 m，設(shè)計(jì)了如上所述的四種井地直流電阻率的觀測方式，圖(4)～圖(7)的虛線表示井設(shè)計(jì)的位置。理論正演得到的合成數(shù)據(jù)(直接采用觀測電位作為數(shù)據(jù)來源)加入3%的高斯白噪聲后作為反演輸入數(shù)據(jù)，此次反演網(wǎng)格與正演網(wǎng)格一致，網(wǎng)格大小為42×20，并向外進(jìn)行了網(wǎng)格延拓以降低截?cái)噙吔鐜淼臄?shù)值誤差。采用Pole-Pole采集裝置，對于井-地測量方式總電極數(shù)為20個(gè)，地表布設(shè)12個(gè)，井中布設(shè)8個(gè)；地-井-地測量方式總電極數(shù)為24個(gè)，地表布設(shè)16個(gè)，井中布設(shè)8個(gè)；井-地-井測量方式總電極數(shù)為32個(gè)，地表布設(shè)16個(gè)，井中布設(shè)16個(gè)；井-井測量方式總電極數(shù)為16個(gè)，地表布設(shè)0個(gè)，井中布設(shè)16個(gè)。最后采用共軛梯度法實(shí)現(xiàn)了不同觀測裝置的井地電阻率2.5維二階平滑模型正則化反演研究。

圖4～圖7分別展示了井-地、地-井-地、井-地-井以及井-井四種觀測方式的反演結(jié)果。此次反演最大迭代次數(shù)為20次，圖8、圖9分別展示數(shù)據(jù)誤差與迭代次數(shù)的衰減圖以及正則化因子α與迭代次數(shù)的變化曲線。從反演結(jié)果中可以看出，反演結(jié)果都能夠較好體現(xiàn)異常所在的位置，反演結(jié)果較為明顯。但是對于不同觀測方式其反演結(jié)果得到細(xì)節(jié)不一樣，對異常體邊界的識(shí)別程度不一致。從圖4的反演結(jié)果可以看出，靠井所在位置的異常體，其反演的結(jié)果呈現(xiàn)的范圍更加集中，而遠(yuǎn)離井的異常體反映的區(qū)域相對較大。圖5展示了地-井-地的觀測裝置的反演結(jié)果，當(dāng)井位于兩個(gè)異常體的中間時(shí)，井所采集的數(shù)據(jù)對兩個(gè)異常體都存在約束，兩個(gè)異常體都能較為明顯分辨。對比圖4～圖5可發(fā)現(xiàn)，地-井-地的觀測方式相對井-地方式提高提高了異常體的分辨率，其對井旁左右的異常體都有較好地約束，而井-地觀測方式對靠近井的異常體約束權(quán)重相對大。

圖6井-地-井觀測方式的反演結(jié)果表明，異常體的邊界相對前兩種的反演結(jié)果更為集中，反演的結(jié)果與實(shí)際模型的吻合度高，反演得到的電阻率值與真實(shí)電阻率值相差不大。圖7的井-井觀測方式的反演結(jié)果表明，其反演得到的異常體范圍與真實(shí)結(jié)果相差較大，由于缺乏地表數(shù)據(jù)的約束導(dǎo)致異常反演趨向井所在的位置，井附近的細(xì)節(jié)較為明顯，總體上異常體相應(yīng)位置得到體現(xiàn)。從圖6～圖7的對比發(fā)現(xiàn)，地表觀測數(shù)據(jù)能夠較為全面進(jìn)行約束，避免了過渡依賴井中數(shù)據(jù)，避免假異常的存在。

綜合分析圖4～圖7的結(jié)果可以表明，井中數(shù)據(jù)能夠較好彌補(bǔ)地表觀測數(shù)據(jù)對縱向約束減弱的缺陷，同時(shí)能夠提高反演非唯一性。從反演的結(jié)果表明，井-地-井的觀測方式對異常體的邊界圈定相對其他三種觀測方式較好，但是存在一些小的細(xì)節(jié)異常，地-井-地觀測方式能夠較好對井兩旁異常體都有較好的約束效果，而遠(yuǎn)離井的異常體其約束能力逐漸降低，這一現(xiàn)象在井-地測量方式得到體現(xiàn)。

圖8、圖9分別展示了數(shù)據(jù)誤差隨迭代次數(shù)的增加而逐漸減小，正則化因子α隨迭代次數(shù)的增加而逐漸減小的這一規(guī)律。當(dāng)模型達(dá)到一定時(shí)，即隨著迭代次數(shù)的增加人為的減少模型約束權(quán)重，使反演結(jié)果更加傾向?qū)?shù)據(jù)的擬合來降低反演的數(shù)據(jù)誤差，最終達(dá)到穩(wěn)定的反演結(jié)果。

4 結(jié)論

不同觀測方式的井地電阻率2.5維共軛梯度反演結(jié)果表明，我們開發(fā)的算法正確可靠，反演算法穩(wěn)定高效。根據(jù)以上的研究可以得出以下幾點(diǎn)認(rèn)識(shí)：

1)井中數(shù)據(jù)能夠較好彌補(bǔ)地表觀測數(shù)據(jù)對縱向約束減弱的缺陷，能夠有效地降低反演非唯一性，同時(shí)井-地-井觀測方式所采用電極數(shù)量大于其他觀測方式，其獲得的地下異常體信息較為豐富，其反演的異常范圍更為集中，其對異常體的圈定要優(yōu)于其他觀測方式。

圖8 數(shù)據(jù)誤差與迭代次數(shù)曲線圖Fig.8 Curve of between iterationserror and iterations

圖9 正則化因子α與迭代次數(shù)變化曲線圖Fig.9 The graph of the change of regularization factor and iteration number

2)由于地表測量電極的缺失，導(dǎo)致井-井觀測方式在橫向上對異常體的約束減少，數(shù)據(jù)量的減少很大程度上增加反演的非唯一性，導(dǎo)致反演精度差，使得井的左右兩側(cè)出現(xiàn)假異常，但也能大致地反應(yīng)異常的范圍，能夠給實(shí)際情況提供參考。

3)基于平滑約束的二階最大平滑穩(wěn)定因子的共軛梯度法迭代反演方法，能夠較為穩(wěn)定獲取反演結(jié)果，同時(shí)模型約束項(xiàng)增加提高反演的穩(wěn)定性，其為后期開發(fā)更為穩(wěn)定模型約束項(xiàng)提供借鑒。