中深部礦產(chǎn)資源勘探井中陣列激電方法研究

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

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

0 引言

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

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

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

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

1 單孔井中激電陣列觀測方法

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

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

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

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

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

當(dāng)點(diǎn)源A置于地下時,波數(shù)域電位V的邊值問題為[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)時,第三類邊界條件D為

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

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

δF(V)=0

(2)

2.2 有限單元法

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

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

(3)

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

2)線性插值。若三角形頂點(diǎn)的局部編號記為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)單元積分。對剖分區(qū)域任意一個三角形單元進(jìn)行積分,有

(5)

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

i≥j

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

i≥j

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

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

(6)

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

KV=S

(7)

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

(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為測點(diǎn)M和N處的電位;Kc為裝置系數(shù),在全空間下Kc=4π/(1/rAM-1/rA*N+1/rA*M-1/rA*N),rA*M和rA*N分別為供電電極A到測量電極M和N的距離,rA*M和rA*N分別為供電電極A關(guān)于地表的鏡像點(diǎn)源A*到測量電極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)的。為提高反演過程的穩(wěn)定性和減少反演的多解性,必須在模型空間引入某種約束[18],這里對模型同時施加總體光滑約束Φ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為對角矩陣(對于模型參數(shù)已知的網(wǎng)格節(jié)點(diǎn),則在相應(yīng)的對角線上給定較大的值(如100),否則設(shè)置為1,若僅施加背景約束,則D為單位矩陣)。

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

(14)

(15)

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

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

(16)

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

(17)

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

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

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

ηa=Aη

(18)

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

(19)

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

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

(20)

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

4 反演算例分析

4.1 模型算例

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

圖3 模型反演算例

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

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

為進(jìn)一步檢驗(yàn)文中井中激電法的探測效果,在長500 cm×寬500 cm×深200 cm的水槽中,開展了井中激電法物理模擬實(shí)驗(yàn),水槽實(shí)驗(yàn)照片如圖4(a)所示。在水槽中注滿清水并添加5 kg食鹽,選取多層鋁板作為低阻高極化模型,電極布置與鋁板的相對位置如圖4(b)所示,按文中第2節(jié)描述的地表供電井中接收的陣列觀測方法采集數(shù)據(jù)。圖4(c)和圖4(d)分別為采集的視電阻率和視極化率擬斷面圖,從圖4中可以看出,在接近鋁板的接收電極處出現(xiàn)低阻高極化異常帶,但對鋁板的傾向的反映不明顯。采用井中約束反演算法對水槽模擬實(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)證了文中井地激電法探測井旁異常體的有效性。

5 結(jié)論

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

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

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

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