中深部礦產資源勘探井中陣列激電方法研究

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

(1.中南大學 有色資源與地質災害探查湖南省重點實驗室,長沙 410083; 2.中南大學 地球科學與信息物理學院,長沙 410083;3.湖南省國土空間調查監測所,長沙 410129;4.陜西鐵道工程勘察有限公司,西安 710000;5.內蒙古自治區呼倫貝爾市應急管理局,呼倫貝爾 021000)

0 引言

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

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

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

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

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

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

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

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

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

2.1 波數域點源場的邊值問題和變分問題

當點源A置于地下時,波數域電位V的邊值問題為[16]

(1)

其中:Ω為積分區域;Γs為地表邊界;?！逓榻財噙吔?σ為介質的電導率;f=Iδ(x-xA)δ(z-zA)/2為點源項。當點源距離截斷邊界較遠時,第三類邊界條件D為

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

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

δF(V)=0

(2)

2.2 有限單元法

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

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

(3)

圖2 地電模型的網格剖分示意圖

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

(4)

其中:V=(V1,V2,V3)T,σ=(σ1,σ2,σ3)T,N=(N1,N2,N3)T,Ni=(aix+biz+ci)/2Δ,(i=1,2,3)為形函數,且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)單元積分。對剖分區域任意一個三角形單元進行積分,有

(5)

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

i≥j

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

i≥j

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

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

(6)

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

KV=S

(7)

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

(8)

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

2.3 計算視電阻率和視極化率

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

(9)

可將前面得到的波數域網格節點電位V轉換成空間域電位U。

根據計算公式

(10)

可計算出三極裝置的視電阻率ρc。其中UM和UN為測點M和N處的電位;Kc為裝置系數,在全空間下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關于地表的鏡像點源A*到測量電極M和N的距離。

(11)

3 井中激電約束反演

3.1 電阻率聯合約束反演

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

AΔm=Δd

(12)

其中:A為偏導數矩陣;Δm為模型參數的改正向量;Δd為數據殘差矢量。

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

(13)

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

當在同一斷面采用 種觀測裝置進行數據采集時,結合式(12)和式(13),在最小二乘意義下構造目標函數Φ為

(14)

(15)

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

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

(16)

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

(17)

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

3.2 極化率聯合約束反演

在完成電阻率反演的基礎上,再進行極化率反演。當地下介質存在激發極化且相對電阻率很小時,視極化率ηa和極化率η之間可近似為線性關系[22]

ηa=Aη

(18)

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

(19)

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

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

(20)

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

4 反演算例分析

4.1 模型算例

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

圖3 模型反演算例

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

4.2 水槽模擬實驗數據反演

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

5 結論

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

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

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

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