柳建新,黃朝宇,3,汪強(qiáng)強(qiáng),佟利群,劉海飛,劉 昕,趙瑩杰

(1.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)沙 410083; 2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院,長(zhǎng)沙 410083;3.湖南省國(guó)土空間調(diào)查監(jiān)測(cè)所,長(zhǎng)沙 410129;4.陜西鐵道工程勘察有限公司,西安 710000;5.內(nèi)蒙古自治區(qū)呼倫貝爾市應(yīng)急管理局,呼倫貝爾 021000)

0 引言

目前我國(guó)礦產(chǎn)資源的保有儲(chǔ)量嚴(yán)重不足,許多老礦山因資源枯竭而面臨倒閉,迫切需要在第二找礦空間 (500 m～2 000 m) 開(kāi)展深部找礦理論與方法的研究,以解決礦山接替資源問(wèn)題[1-2]。我國(guó)大中型危機(jī)礦山深邊部具有很大的找礦潛力,充分利用礦山原有的成礦模式及豐富的開(kāi)采資料,深入研究探測(cè)深度大、分辨率高的地球物理新方法與新技術(shù),是目前開(kāi)展深部找礦的一種有效途徑。

井中物探方法作為地球物理勘探方法的一個(gè)重要分支,主要用來(lái)解決井周地質(zhì)問(wèn)題,諸如尋找井旁、井底盲礦體,確定其空間位置、形態(tài)、產(chǎn)狀,追蹤和圈定礦體范圍,研究鉆孔間礦體的連續(xù)性等,其突出優(yōu)點(diǎn)就是能夠把場(chǎng)源或測(cè)量裝置借助鉆孔放入地下深處,使其接近探測(cè)對(duì)象,因此發(fā)現(xiàn)井旁隱伏礦體的能力往往比地面物探方法要強(qiáng)[3]。目前在危機(jī)礦山深邊部找礦中,主要使用的井中物探方法有井中激電法、井中瞬變電磁法、井中磁測(cè)及井中重力等[4],其中井中激電法因其觀測(cè)參數(shù)多、觀測(cè)裝置靈活以及對(duì)礦質(zhì)資源反應(yīng)靈敏等特點(diǎn),在中深部礦產(chǎn)資源勘探中具有明顯的優(yōu)勢(shì)[5]。

井中激電法的探測(cè)效果是與其正反演解釋水平密切相關(guān)的。近年來(lái),井中激電法的正反演解釋方法得到了快速的發(fā)展。J.P.Busby等[6]模擬并分析了井中三極觀測(cè)裝置的多個(gè)極化體的電阻率和極化率響應(yīng);Klaus Spitzer等[7]采用跨孔偶極和單孔偶極直流激電觀測(cè)方法,對(duì)魁北克西北部Casa Berardi金礦進(jìn)行了井中激電探測(cè),并采用三維電阻率和極化率正演模擬解釋來(lái)描述斷層構(gòu)造;Kulessa Bernd等[8]利用鉆孔電阻率數(shù)據(jù)的三維反演對(duì)冰川下排水條件進(jìn)行時(shí)移成像;呂玉增等[9]對(duì)三維地-井、井-地IP的正反演做了系統(tǒng)研究,編制了地-井、井-地IP三維快速正反演計(jì)算程序,開(kāi)發(fā)了地-井五方位IP人機(jī)交互正演擬合反演解釋軟件;F.J.Morrow等[10]基于直流電阻率穿越和跨孔電阻率層析成像法,研究了潮汐對(duì)淺層無(wú)限制含水層鹽分界面的影響;李長(zhǎng)偉等[11]實(shí)現(xiàn)了三維井中激發(fā)極化法正反演,并對(duì)快速迭代求解技術(shù)進(jìn)行了研究;熊彬等[12]基于有限差分和異常電位算法實(shí)現(xiàn)了時(shí)間域井中激發(fā)極化法的三維數(shù)值模擬,采用一維非零元素行壓縮存儲(chǔ)模式對(duì)所形成的大型稀疏矩陣進(jìn)行存儲(chǔ),節(jié)約了計(jì)算所需內(nèi)存空間,同時(shí)引入不完全Cholesky分解穩(wěn)定化雙共軛梯度(ICBU)法求解有限差分線性方程組,提高了求解的效率;白澤等[13]采用了有限差分法和不完全Cholesky共軛梯度法(ICCG)對(duì)井地電成像的點(diǎn)源和線源在產(chǎn)生的地表電勢(shì)異常進(jìn)行了三維正演模擬,并對(duì)地層電阻率數(shù)據(jù)進(jìn)行阻尼最小二乘法反演;H.D.Wondimu等[14]在加拿大安大略省北部一個(gè)巖漿硫化物礦床的一組19個(gè)鉆孔中進(jìn)行了梯度和充電法調(diào)查,并對(duì)數(shù)據(jù)反演得到三維電阻率和極化率模型,圈定了礦化帶的輪廓;趙榮春等[15]采用基于柱坐標(biāo)下的放射狀網(wǎng)格剖分方式,通過(guò)異常電位法對(duì)井中激電觀測(cè)進(jìn)行了有限元數(shù)值模擬,對(duì)井中激電觀測(cè)井旁異常體所引起的激電異常特征進(jìn)行了研究。這些研究成果推動(dòng)了井中激電法理論與方法的發(fā)展,也為筆者開(kāi)展井中正反演方法的研究奠定了基礎(chǔ)。

采用井中激電法解決危機(jī)礦山的中深部找礦問(wèn)題,首先需要獲取足夠多的反映深部礦體的信息,而這種深部信息的獲取必須充分利用礦體及其外圍的發(fā)射和接收空間。筆者針對(duì)單孔井中激電法進(jìn)行了初步研究,首先設(shè)計(jì)了兩種多電極系觀測(cè)方法,以盡可能多地獲取井旁地電信息。其次系統(tǒng)研究了井中激電法的二維有限元正演模擬方法以及多種觀測(cè)方式的聯(lián)合反演方法,主要包括點(diǎn)源場(chǎng)邊值和變分問(wèn)題、有限元解法、構(gòu)建反演區(qū)域、聯(lián)合反演方程的建立及約束方式等。最后為驗(yàn)證方法的有效性,編制了Windows界面下的反演軟件,對(duì)地電模型進(jìn)行反演試算,反演效果較好,驗(yàn)證了該方法的可行性。

1 單孔井中激電陣列觀測(cè)方法

在礦區(qū)存在鉆孔的情況下,可以充分利用地面和鉆孔的觀測(cè)空間,獲取大量井旁地電信息,以提高井旁目標(biāo)體的探測(cè)效果?？紤]到探測(cè)分辨率和野外觀測(cè)效率問(wèn)題,筆者在地表三極觀測(cè)裝置的基礎(chǔ)上,設(shè)計(jì)了兩種針對(duì)單個(gè)鉆孔情況下的井中激電陣列觀測(cè)方法。

為提高探測(cè)的縱向分辨率,設(shè)計(jì)的井中激電陣列觀測(cè)方法如圖1(a)所示。觀測(cè)原理為:將供電電極B(負(fù)極)作為無(wú)窮遠(yuǎn)極,放置于離井口較遠(yuǎn)的位置,其余供電電極A1、A2、…、An(均為正極,相鄰電極以等算數(shù)間隔或等對(duì)數(shù)間隔方式排列,可根據(jù)勘探要求進(jìn)行設(shè)計(jì))及測(cè)量電極M1和M2置于井中,在井中由下至上進(jìn)行供電和測(cè)量。即A1供電M1M2測(cè)量,A2供電M1M2測(cè)量,…,An供電M1M2測(cè)量,第一個(gè)測(cè)深點(diǎn)測(cè)量結(jié)束;供電電極A1、A2、…、An及測(cè)量電極M1和M2整體向上移動(dòng)一個(gè)點(diǎn)位,重復(fù)A1供電M1M2測(cè)量,A2供電M1M2測(cè)量,…,An供電M1M2測(cè)量,第二個(gè)測(cè)深點(diǎn)測(cè)量結(jié)束;再繼續(xù)上移一個(gè)點(diǎn)位,再繼續(xù)供電和測(cè)量,如此往復(fù),直到整口井測(cè)量結(jié)束(或者僅測(cè)量某異常深度段),這樣可獲取井旁不同深度范圍內(nèi)的電性信息。

為提高探測(cè)的橫向分辨率,設(shè)計(jì)的井中激電陣列觀測(cè)方法如圖1(b)所示。觀測(cè)原理為:將供電電極B(負(fù)極)作為無(wú)窮遠(yuǎn)極,放置于離井口較遠(yuǎn)的位置,其余供電電極A1、A2、…、An(均為正極)布設(shè)于過(guò)井口且垂直異常走向的同一測(cè)線上,而測(cè)量電極M1、M2、…、Mn由下至上置于井中,地面電極供電井中電極進(jìn)行差分測(cè)量。即A1供電M1M2、M2M3、…、Mn-1Mn同時(shí)測(cè)量;A2供電M1M2、M2M3、…、Mn-1Mn同時(shí)測(cè)量;直到An供電M1M2、M2M3…、Mn-1Mn同時(shí)測(cè)量,至此觀測(cè)結(jié)束。這種觀測(cè)方式有利于分辨異常體位于鉆孔的哪一側(cè),實(shí)際應(yīng)用中盡量使供電電極排列大于異常體的橫向分布范圍。

2 井中激電2.5維有限元正演模擬

2.1 波數(shù)域點(diǎn)源場(chǎng)的邊值問(wèn)題和變分問(wèn)題

當(dāng)點(diǎn)源A置于地下時(shí),波數(shù)域電位V的邊值問(wèn)題為[16]

(1)

其中:Ω為積分區(qū)域;Γs為地表邊界;?！逓榻?cái)噙吔?σ為介質(zhì)的電導(dǎo)率;f=Iδ(x-xA)δ(z-zA)/2為點(diǎn)源項(xiàng)。當(dāng)點(diǎn)源距離截?cái)噙吔巛^遠(yuǎn)時(shí),第三類邊界條件D為

式中:rp和rq分別為點(diǎn)源p及其相對(duì)地表的鏡像點(diǎn)q至邊界積分點(diǎn)的距離;cos(rp,n)和cos(rq,n)分別為矢徑rp和rq與積分區(qū)域Ω的外法向量n的夾角余弦;λ為波數(shù);K0和K1分別為第二類零階和一階修正貝塞爾函數(shù)。

利用泛函分析,將方程(1)轉(zhuǎn)化為等價(jià)的變分問(wèn)題

δF(V)=0

(2)

2.2 有限單元法

采用有限元法求解變分方程(2),求解過(guò)程如下:

1)單元剖分?？紤]到地形因素,將整個(gè)研究區(qū)域剖分成有限個(gè)三角形,具體如圖2所示。方程(2)對(duì)區(qū)域Ω 的積分將轉(zhuǎn)化為對(duì)各三角形單元e和邊界單元Γe的積分之和:

(3)

圖2 地電模型的網(wǎng)格剖分示意圖

2)線性插值。若三角形頂點(diǎn)的局部編號(hào)記為1、2、3,其空間坐標(biāo)記為(x1,z1)、(x2,z2)、(x2,z2),電位值記為V1、V2、V3,電導(dǎo)率值記為σ1、σ2、σ3。在三角形單元內(nèi),任意一點(diǎn)的電位和電導(dǎo)率均采用線性插值,即:

(4)

其中:V=(V1,V2,V3)T,σ=(σ1,σ2,σ3)T,N=(N1,N2,N3)T,Ni=(aix+biz+ci)/2Δ,(i=1,2,3)為形函數(shù),且a1=z2-z3,a2=z3-z1,a3=z1-z2,b1=x3-x2,b2=x1-x3,b3=x2-x1,c1=x2z3-x3z2,c2=x3z1-x1z3,c3=x1z2-x2z1,Δ=(a1b2-a2b1)/2為三角形單元的面積。

3)單元積分。對(duì)剖分區(qū)域任意一個(gè)三角形單元進(jìn)行積分,有

(5)

其中:K1e為三角形單元積分的剛度矩陣,K1e=[k1eij]=[k1eij],i、j=1、2、3,并且有

i≥j

式中:α=(σ1+σ2+σ3)/12Δ,β=λ2Δ/60。K2e為無(wú)窮遠(yuǎn)邊界線單元積分的剛度矩陣,K2e=[k2eij]=[K2eji],i、j=1、2、3,假如邊界單元落在三角形12邊上,即i、j=3時(shí),k2eij=0,則邊界積分的下三角陣的非零元素為

i≥j

式中:γ=D·l12/12,l12為節(jié)點(diǎn)1和2之間的距離。對(duì)于場(chǎng)源項(xiàng)Se=[si],i=1,2,3,如果點(diǎn)源與三角形某個(gè)節(jié)點(diǎn)i重合,則si=0.5,否則si=0。

4)總體合成。對(duì)所有三角單元積分結(jié)果求和,即得到泛函F(V)的數(shù)值表達(dá)式

(6)

其中:K是由全部三角單元和邊界單元的(K1e+K2e)相加組成的M×M階對(duì)稱系數(shù)矩陣,各項(xiàng)元素與模型電阻率分布和網(wǎng)格剖分在關(guān);V是由所有M個(gè)三角網(wǎng)格節(jié)點(diǎn)上的波數(shù)域電位組成的列矢量;S則是由Se相加組成的與場(chǎng)源有關(guān)的列矢量。對(duì)方程(6)中泛函F(V)求變分,并令其為零,進(jìn)而可得線性方程組

KV=S

(7)

利用一維變帶寬壓縮存儲(chǔ)的喬里斯基分解法求解方程(7),即可得到所有網(wǎng)格節(jié)點(diǎn)的波數(shù)域電位V,再通過(guò)傅氏逆變換

(8)

便可將波數(shù)域電位V轉(zhuǎn)換成主剖面y=0上的空間域電位U(x,0,z)。

2.3 計(jì)算視電阻率和視極化率

利用式(8)的近似傅氏逆變換公式

(9)

可將前面得到的波數(shù)域網(wǎng)格節(jié)點(diǎn)電位V轉(zhuǎn)換成空間域電位U。

根據(jù)計(jì)算公式

(10)

可計(jì)算出三極裝置的視電阻率ρc。其中UM和UN為測(cè)點(diǎn)M和N處的電位;Kc為裝置系數(shù),在全空間下Kc=4π/(1/rAM-1/rA*N+1/rA*M-1/rA*N),rA*M和rA*N分別為供電電極A到測(cè)量電極M和N的距離,rA*M和rA*N分別為供電電極A關(guān)于地表的鏡像點(diǎn)源A*到測(cè)量電極M和N的距離。

(11)

3 井中激電約束反演

3.1 電阻率聯(lián)合約束反演

電阻率的線性反演方程通?？杀硎緸閇17]:

AΔm=Δd

(12)

其中:A為偏導(dǎo)數(shù)矩陣;Δm為模型參數(shù)的改正向量;Δd為數(shù)據(jù)殘差矢量。

實(shí)際中方程(12)通常是欠定且病態(tài)的。為提高反演過(guò)程的穩(wěn)定性和減少反演的多解性,必須在模型空間引入某種約束[18],這里對(duì)模型同時(shí)施加總體光滑約束Φms和屬性約束Φmz,可表示為

(13)

其中:m為反演模型參數(shù)向量;mz為已知屬性模型參數(shù)向量(若網(wǎng)格節(jié)點(diǎn)的模型參數(shù)已知,則在相應(yīng)的節(jié)點(diǎn)上給定近似真值,否則用背景模型參數(shù)代替),C為光滑度矩陣,可采用模型網(wǎng)格節(jié)點(diǎn)間的距離構(gòu)建[19-20],D為對(duì)角矩陣(對(duì)于模型參數(shù)已知的網(wǎng)格節(jié)點(diǎn),則在相應(yīng)的對(duì)角線上給定較大的值(如100),否則設(shè)置為1,若僅施加背景約束,則D為單位矩陣)。

當(dāng)在同一斷面采用 種觀測(cè)裝置進(jìn)行數(shù)據(jù)采集時(shí),結(jié)合式(12)和式(13),在最小二乘意義下構(gòu)造目標(biāo)函數(shù)Φ為

(14)

(15)

采用共軛梯度法解方程(15)[21],得到模型參的修正量Δm,將其代入式(16),

m(k)=m(k-1)+Δm

(16)

便得到新的預(yù)測(cè)模型參數(shù)向量m(k)。經(jīng)過(guò)多次迭代,直至實(shí)測(cè)數(shù)據(jù)和模擬數(shù)據(jù)之間的平均均方誤差

(17)

滿足要求或迭代次數(shù)滿足終止條件為止,電阻率反演過(guò)程結(jié)束。

3.2 極化率聯(lián)合約束反演

在完成電阻率反演的基礎(chǔ)上,再進(jìn)行極化率反演。當(dāng)?shù)叵陆橘|(zhì)存在激發(fā)極化且相對(duì)電阻率很小時(shí),視極化率ηa和極化率η之間可近似為線性關(guān)系[22]

ηa=Aη

(18)

其中:A為偏導(dǎo)數(shù)矩陣(aij=?lnρci/?lnρj),在電阻率反演中已經(jīng)得到。當(dāng)在同一斷面采用n種觀測(cè)裝置進(jìn)行數(shù)據(jù)采集時(shí),在最小二乘意義下構(gòu)造極化率反演的目標(biāo)函數(shù)Φ,并在目標(biāo)函數(shù)中引入光滑和已知屬性約束,有

(19)

其中極化率反演的數(shù)據(jù)加權(quán)系數(shù)fk根據(jù)最后一次電阻率反演迭代的數(shù)據(jù)擬合誤差給定,ηz為井旁已知極化率參數(shù)向量,其余參量含義與上述相同。

將式(19)兩端對(duì)η求導(dǎo)并令其等于零,即可得到極化率的線性反演方程

(20)

采用共軛梯度法求解方程(20),便得到地下介質(zhì)的極化率。

4 反演算例分析

4.1 模型算例

假定在地下隱伏巖體中賦存三條低阻高極化脈狀礦體,巖層及礦脈均傾向于小點(diǎn)號(hào)方向。地表剖面長(zhǎng)為800 m,井深為800 m,地表布設(shè)42根電極,井中布設(shè)160根電極。電阻率和極化率模型如圖3(a)～圖3(b)所示?；谠摰仉娔Ｐ蛯?duì)上述觀測(cè)裝置的激電數(shù)據(jù)分別進(jìn)行二維有限元數(shù)值模擬,井中觀測(cè)方式的視電阻率和視極化率的模擬結(jié)果分別如圖3(c)～圖3(d)所示,從圖3中可看出,在礦脈位置異常呈低阻高極化特征,但位置偏下及難于推斷礦脈位于井的哪一側(cè)。地井觀測(cè)方式的視電阻率和視極化率的模擬結(jié)果分別如圖3(e)～圖3(f)所示,相對(duì)井中觀測(cè)方式而言,這種方式則更好地反映了斷面的電性的分布特征,易于推斷異常大致形態(tài)和埋深,但仍存在向下偏移的現(xiàn)象。

圖3 模型反演算例

為提高兩種觀測(cè)方式的探測(cè)效果,采用文中聯(lián)合反演方法將兩種觀測(cè)方式的激電數(shù)據(jù)聯(lián)合起來(lái)進(jìn)行反演,當(dāng)對(duì)井旁地質(zhì)結(jié)構(gòu)信息不了解時(shí),則在反演中只施加光滑和背景先驗(yàn)約束信息,電阻率和極化率的反演結(jié)果分別如圖3(g)～圖3(h)所示,地層傾向及異常形態(tài)和埋深得到了較好地反映,但也不同程度地出現(xiàn)了冗余信息。當(dāng)對(duì)地層結(jié)構(gòu)及電性信息有所了解時(shí),可將其作為約束引入到反演中,反演時(shí)對(duì)上覆三層電性層給定準(zhǔn)確的電阻率和極化率值,式(15)中約束矩陣對(duì)應(yīng)的對(duì)角線元素設(shè)置約束值為100,而下覆兩層電性層給定均勻電性值2 500,約束值為1,電阻率和極化率的反演結(jié)果分別如圖3(i)～圖3(j)所示,與僅施加先驗(yàn)約束信息的反演結(jié)果相比,地層結(jié)構(gòu)、電性及異常形態(tài)均得到了明顯的改善。在反演中,可根據(jù)對(duì)地層結(jié)構(gòu)和電性信息的了解程度,通過(guò)人機(jī)交互改變約束值和阻尼的大小,改善反演效果。

4.2 水槽模擬實(shí)驗(yàn)數(shù)據(jù)反演

為進(jìn)一步檢驗(yàn)文中井中激電法的探測(cè)效果,在長(zhǎng)500 cm×寬500 cm×深200 cm的水槽中,開(kāi)展了井中激電法物理模擬實(shí)驗(yàn),水槽實(shí)驗(yàn)照片如圖4(a)所示。在水槽中注滿清水并添加5 kg食鹽,選取多層鋁板作為低阻高極化模型,電極布置與鋁板的相對(duì)位置如圖4(b)所示,按文中第2節(jié)描述的地表供電井中接收的陣列觀測(cè)方法采集數(shù)據(jù)。圖4(c)和圖4(d)分別為采集的視電阻率和視極化率擬斷面圖,從圖4中可以看出,在接近鋁板的接收電極處出現(xiàn)低阻高極化異常帶,但對(duì)鋁板的傾向的反映不明顯。采用井中約束反演算法對(duì)水槽模擬實(shí)驗(yàn)數(shù)據(jù)進(jìn)行反演計(jì)算,得到電阻率與極化率的二維斷面分別如圖4(c)和圖4(d),從中可以看出,在X方向0 cm～30 cm,埋深20 cm～60 cm范圍內(nèi)呈低電阻率和高極化率特征,異常形態(tài)與實(shí)際模型吻合度較高,驗(yàn)證了文中井地激電法探測(cè)井旁異常體的有效性。

5 結(jié)論

筆者針對(duì)井旁隱伏礦體勘查問(wèn)題,對(duì)單孔井地陣列激電探測(cè)方法進(jìn)行了研究,得到如下結(jié)論:

1)設(shè)計(jì)了井中和地-井多電極系陣列觀測(cè)方法,該方法具有觀測(cè)效率高、獲取的信息量豐富、縱向和橫向分辨率高等優(yōu)點(diǎn)。

2)研究了井中激電法的二維有限元正演模擬方法以及多種觀測(cè)方式的聯(lián)合反演方法,并對(duì)兩種電阻率數(shù)據(jù)進(jìn)行聯(lián)合反演解釋,反演效果良好,地層傾向、異常形態(tài)以及目標(biāo)體埋深均得到較好地反映,可以嘗試將觀測(cè)方法和反演解釋方法在野外生產(chǎn)實(shí)踐中應(yīng)用。

3)本項(xiàng)研究可進(jìn)一步拓展到單孔井中、地井多電極系三維陣列觀測(cè)及反演解釋中。