龍秀潔,陳漢波,莫亞軍,區(qū)小毅,盧勝輝
(1.廣西壯族自治區(qū)地球物理勘察院,廣西 柳州 545005;2.桂林理工大學(xué) 地球科學(xué)學(xué)院,廣西 桂林 541006)
復(fù)電阻率法是一種重要的地球物理勘探方法,具有電性參數(shù)多的特點,在金屬礦勘探領(lǐng)域起著重要的作用,較其他物探方法具有抗干擾能力強,多參數(shù)對比解釋可提供更豐富的異常信息的優(yōu)點。近年來,該方法被廣泛應(yīng)用于石油地震勘探、水文地質(zhì)調(diào)查等領(lǐng)域[1]。
國內(nèi)外學(xué)者對于復(fù)電阻率法數(shù)值模擬做了大量的研究,取得眾多成果。如國內(nèi)的羅延鐘等[2-3]、張桂清等[4]、劉菘等[5]。國外的研究者,如G.W. Hohmann等[6]、A. Weller等[7]、D.W. Oldenburg等[8]等。然而,大多對于頻譜激電法的研究均基于Core-Core模型,且沒有考慮復(fù)電導(dǎo)率呈線性變化的情況。因此,本文在前人的基礎(chǔ)上,基于異常復(fù)電位法,采用Fractal模型[9]研究復(fù)電導(dǎo)率呈線性變化的復(fù)電阻率法2.5D有限元數(shù)值模擬,并對其相應(yīng)的異常響應(yīng)特征進行分析,以期能為復(fù)電阻率法野外實際工作提供理論依據(jù)。
Fractal 模型是由B.R.P. Rocha提出用來描述巖、礦石激發(fā)極化現(xiàn)象的數(shù)學(xué)模型,B.R.P. Rocha等[9]對大量巖、礦石標本及露頭進行測量,結(jié)果表明,巖、礦石的激發(fā)極化響應(yīng)所引起的復(fù)電阻率可用以下公式進行表述:
(1)
其中:ρ0為巖、礦石未極化的直流電阻率;m為介質(zhì)的極化率(無量綱);δr為礦物顆粒電阻率百分比;rh=1/(1+iωτ0);u=iωτ(1+v);v=(iωτf)-η;τ為雙層振蕩相關(guān)的松弛時間常數(shù);τ0為采樣松弛時間常數(shù);τf為分形松弛時間;η為介質(zhì)的分形幾何直接相關(guān)的參數(shù)(無量綱)。
式(1)對于定量描述巖礦石的激發(fā)極化現(xiàn)象提供了一種新的手段,相對于常用的core-core模型來說,其提供了更多的巖石物性參數(shù),從而可以更詳細地解釋巖礦石激發(fā)極化效應(yīng)所產(chǎn)生的原因及影響因素。圖1展示了不同η的復(fù)電阻率的振幅和相位隨著頻率的變化規(guī)律。計算參數(shù)為:ρ0=100 Ω·m、m=0.5、δr=1.0、τ=10-6s、τ0=10-12s、τf=100、η分別取0.05、0.20、0.35、0.5、0.65、0.80、0.95。由圖1可以看得到復(fù)電祖率ρ(ω)幅值隨著頻率的增大而遞減,在低頻時,且η取較小值時,復(fù)電祖率幅值趨向于ρ0。而相位φ的絕對值與η成正比,在低頻及高頻時趨向于0。
點源條件下,復(fù)電阻率法2.5D變分問題如下所示:
(2)
(3)
式中:I為電流大小。
采用矩形網(wǎng)格進一步細分為4個三角形網(wǎng)格剖分方式對研究區(qū)域進行剖分,如圖2所示:將方程(2)對研究區(qū)域Ω和無窮邊界?!薜姆e分分為各個單元的積分之和:
(4)
在單元內(nèi)對復(fù)電位及復(fù)電導(dǎo)率進行線性插值。線性插值的母單元與子單元如圖3所示。則單元內(nèi)的復(fù)電位及復(fù)電導(dǎo)率則可表示為
圖3 線性插值的母單元(a)及子單元(b)
(5)
σ=N1σ1+N2σ2+N3σ3
,
(6)
(7)
方程(4)中的第1項
(8)
其中:
(9)
方程(4)積分第2項
(10)
其中
(11)
方程(4)積分第3項
(12)
其中:
(13)
方程(4)積分第4項
(14)
其中:
(15)
方程(4)積分第5項
(16)
其中:
(17)
方程(4)積分第6項
(18)
式中:
(19)
在單元內(nèi),將式(9)、式(11)、式(13)、式(15)、式(17)、式(19)的積分結(jié)果進行相加后,再擴展成由全體節(jié)點組成的矩陣或列陣:
(20)
對方程(20) 進行求變分:
(21)
(22)
(23)
求解方程(22)便可得波數(shù)域中的異常復(fù)電位值。本文采用不完全LU分解的穩(wěn)定雙共軛梯度算法來求解線性方程組,具體求解過程請詳見文獻[10],在此不再贅述。求解出異常復(fù)電位后,再通過傅里葉反變換便可得到三維空間中的異常復(fù)電位。
求解出三維空間中的異常復(fù)電位,便可按照以下公式計算復(fù)電阻率及相位值:
(24)
(25)
為了檢驗文中算法及所編制的程序的正確性,設(shè)置一個兩層水平層狀極化介質(zhì)模型,厚度為6 m,電阻率為100 Ω。Fractal模型參數(shù)如表1所示。供電頻率f=0.125 Hz,將2.5D有限元數(shù)值解與一維數(shù)字濾波算法所計算出來的解析解進行對比,如圖4所示。由圖4可知,復(fù)電阻率的實分量和虛分量數(shù)值解與解析解基本一致,從而驗證了本文所述算法及所編制的程序的正確性。
表1 模型1參數(shù)
圖4 2.5維有限元數(shù)值解與一維數(shù)字濾波解析解對比
模型2 為極化半空間中有一異常體。具體幾何參數(shù)如圖5所示,F(xiàn)ractal模型參數(shù)如表2所示。圖6、圖7為頻率為0.125、0.25、0.5、1.0 Hz的復(fù)電阻率法2.5D有限元數(shù)值模擬結(jié)果。由圖6、圖7可看到,復(fù)電阻率的各分量都很好地反映出異常體的存在,且異常響應(yīng)范圍與異常體的埋藏位置基本吻合。在計算的頻率范圍內(nèi),復(fù)電阻率的相位及虛分量均為負值,偶極裝置下的復(fù)電阻率振幅和相位均呈“八”字型異常特征,特別值得注意的是復(fù)電阻率的虛分量,其幅值隨著頻率的增大而減小,且遠小于實分量。
圖5 模型2斷面示意
表2 模型2的Fractal 模型參數(shù)
a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz
a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz
為了進一步驗證本文算法的有效性,設(shè)計一個極化半空間中有雙異常體模型。異常體大小均為12 m×6 m,埋深維4 m,兩者相距12 m。異常體電導(dǎo)率為0.01 s/m,背景介質(zhì)的電導(dǎo)率為0.001 s/m。具體幾何參數(shù)如圖8所示,F(xiàn)ractal模型參數(shù)如表3所示。圖9~12為頻率為0.125、0.25、0.5、1.0 Hz的不同測量裝置下的復(fù)電阻率法2.5D有限元數(shù)值模擬結(jié)果。由圖9~12可看到,復(fù)電阻率的各分量都很好地反映出異常體的存在,且異常響應(yīng)范圍與異常體的埋藏位置基本吻合。在計算的頻率范圍內(nèi),4種測量裝置下的視復(fù)電阻率的虛分量均為負值;對稱四極裝置下,復(fù)電阻率的實部和虛部呈現(xiàn)雙“柱體”型異常特征,如圖9所示;如圖10所示,偶極裝置下的復(fù)電阻率實部和虛部均呈雙“八”字型異常特征;從圖11可知,三極測量裝置,復(fù)電阻率的實部和虛部的異常特征與對稱四極類似,不過虛部的右側(cè)部分的值大于左側(cè);如圖12可知,二極裝置下,可以看到復(fù)電阻率的異常位置與模型非常吻合,其縱向分辨率較于三極和對稱四極測量裝置高。同樣,我們可以觀察到,在4種測量裝置下,復(fù)電阻率的實部隨著頻率的增大,產(chǎn)生的變化不大;不同于實部分量,復(fù)電阻率的虛分量,其幅值隨著頻率的增大而減小,且遠小于實分量。
圖8 模型3斷面示意
表3 模型3的Fractal 模型參數(shù)
a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz
a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz
a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz
a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz
本文實現(xiàn)了復(fù)電阻率法2.5D有限元數(shù)值模擬,采用異常復(fù)電法,避免了源的奇異性的影響,提高了數(shù)值解的精度;引入Fractal 模型對頻率域激發(fā)極化響應(yīng)進行研究,可以為更好解釋激發(fā)極化效應(yīng)提供一種新的手段;對復(fù)電導(dǎo)率及復(fù)電位進行線性插值,可以模擬更為真實的野外地質(zhì)情況;分析了所設(shè)計的地電模型的復(fù)電祖率異常響應(yīng)特征,可為野外實際勘探提供理論依據(jù)。數(shù)值模擬結(jié)果驗證了本文基于Fractal 模型的CR2.5D有限元數(shù)值模擬算法的正確性和精確性,表明引入Fractal 模型研究復(fù)電導(dǎo)率呈線性變化的CR響應(yīng)特征是可行的、有效的。