基于Fractal模型的復(fù)電阻率法2.5D有限元數(shù)值模擬

龍秀潔,陳漢波,莫亞軍,區(qū)小毅,盧勝輝

(1.廣西壯族自治區(qū)地球物理勘察院,廣西 柳州 545005；2.桂林理工大學(xué) 地球科學(xué)學(xué)院,廣西 桂林 541006)

0 引言

復(fù)電阻率法是一種重要的地球物理勘探方法，具有電性參數(shù)多的特點(diǎn)，在金屬礦勘探領(lǐng)域起著重要的作用，較其他物探方法具有抗干擾能力強(qiáng)，多參數(shù)對(duì)比解釋可提供更豐富的異常信息的優(yōu)點(diǎn)。近年來，該方法被廣泛應(yīng)用于石油地震勘探、水文地質(zhì)調(diào)查等領(lǐng)域[1]。

國內(nèi)外學(xué)者對(duì)于復(fù)電阻率法數(shù)值模擬做了大量的研究，取得眾多成果。如國內(nèi)的羅延鐘等[2-3]、張桂清等[4]、劉菘等[5]。國外的研究者，如G.W. Hohmann等[6]、A. Weller等[7]、D.W. Oldenburg等[8]等。然而，大多對(duì)于頻譜激電法的研究均基于Core-Core模型，且沒有考慮復(fù)電導(dǎo)率呈線性變化的情況。因此，本文在前人的基礎(chǔ)上，基于異常復(fù)電位法，采用Fractal模型[9]研究復(fù)電導(dǎo)率呈線性變化的復(fù)電阻率法2.5D有限元數(shù)值模擬，并對(duì)其相應(yīng)的異常響應(yīng)特征進(jìn)行分析，以期能為復(fù)電阻率法野外實(shí)際工作提供理論依據(jù)。

1 復(fù)電阻率法邊值問題

1.1 Fractal 模型

Fractal 模型是由B.R.P. Rocha提出用來描述巖、礦石激發(fā)極化現(xiàn)象的數(shù)學(xué)模型，B.R.P. Rocha等[9]對(duì)大量巖、礦石標(biāo)本及露頭進(jìn)行測(cè)量，結(jié)果表明，巖、礦石的激發(fā)極化響應(yīng)所引起的復(fù)電阻率可用以下公式進(jìn)行表述：

(1)

其中：ρ0為巖、礦石未極化的直流電阻率；m為介質(zhì)的極化率(無量綱)；δr為礦物顆粒電阻率百分比；rh=1/(1+iωτ0)；u=iωτ(1+v)；v=(iωτf)-η；τ為雙層振蕩相關(guān)的松弛時(shí)間常數(shù)；τ0為采樣松弛時(shí)間常數(shù)；τf為分形松弛時(shí)間；η為介質(zhì)的分形幾何直接相關(guān)的參數(shù)(無量綱)。

式(1)對(duì)于定量描述巖礦石的激發(fā)極化現(xiàn)象提供了一種新的手段，相對(duì)于常用的core-core模型來說，其提供了更多的巖石物性參數(shù)，從而可以更詳細(xì)地解釋巖礦石激發(fā)極化效應(yīng)所產(chǎn)生的原因及影響因素。圖1展示了不同η的復(fù)電阻率的振幅和相位隨著頻率的變化規(guī)律。計(jì)算參數(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ù)電祖率ρ(ω)幅值隨著頻率的增大而遞減，在低頻時(shí)，且η取較小值時(shí)，復(fù)電祖率幅值趨向于ρ0。而相位φ的絕對(duì)值與η成正比，在低頻及高頻時(shí)趨向于0。

1.2 基于Fractal模型的復(fù)電阻率變分問題

點(diǎn)源條件下，復(fù)電阻率法2.5D變分問題如下所示：

(2)

(3)

式中：I為電流大小。

2 有限單元法分析

2.1 有限單元剖分

采用矩形網(wǎng)格進(jìn)一步細(xì)分為4個(gè)三角形網(wǎng)格剖分方式對(duì)研究區(qū)域進(jìn)行剖分，如圖2所示：將方程(2)對(duì)研究區(qū)域Ω和無窮邊界Γ∞的積分分為各個(gè)單元的積分之和：

(4)

在單元內(nèi)對(duì)復(fù)電位及復(fù)電導(dǎo)率進(jìn)行線性插值。線性插值的母單元與子單元如圖3所示。則單元內(nèi)的復(fù)電位及復(fù)電導(dǎo)率則可表示為

圖3 線性插值的母單元(a)及子單元(b)

(5)

σ=N1σ1+N2σ2+N3σ3

，

(6)

(7)

2.2 單元分析

方程(4)中的第1項(xiàng)

(8)

其中：

(9)

方程(4)積分第2項(xiàng)

(10)

其中

(11)

方程(4)積分第3項(xiàng)

(12)

其中:

(13)

方程(4)積分第4項(xiàng)

(14)

其中：

(15)

方程(4)積分第5項(xiàng)

(16)

其中：

(17)

方程(4)積分第6項(xiàng)

(18)

式中：

(19)

2.3 總體合成

在單元內(nèi)，將式(9)、式(11)、式(13)、式(15)、式(17)、式(19)的積分結(jié)果進(jìn)行相加后，再擴(kuò)展成由全體節(jié)點(diǎn)組成的矩陣或列陣：

(20)

2.4 求解線性方程組

對(duì)方程(20) 進(jìn)行求變分：

(21)

(22)

(23)

求解方程(22)便可得波數(shù)域中的異常復(fù)電位值。本文采用不完全LU分解的穩(wěn)定雙共軛梯度算法來求解線性方程組，具體求解過程請(qǐng)?jiān)斠娢墨I(xiàn)[10]，在此不再贅述。求解出異常復(fù)電位后，再通過傅里葉反變換便可得到三維空間中的異常復(fù)電位。

2.5 復(fù)電阻率的計(jì)算

求解出三維空間中的異常復(fù)電位，便可按照以下公式計(jì)算復(fù)電阻率及相位值：

(24)

(25)

3 模型計(jì)算

3.1 模型1

為了檢驗(yàn)文中算法及所編制的程序的正確性，設(shè)置一個(gè)兩層水平層狀極化介質(zhì)模型，厚度為6 m，電阻率為100 Ω。Fractal模型參數(shù)如表1所示。供電頻率f=0.125 Hz，將2.5D有限元數(shù)值解與一維數(shù)字濾波算法所計(jì)算出來的解析解進(jìn)行對(duì)比，如圖4所示。由圖4可知，復(fù)電阻率的實(shí)分量和虛分量數(shù)值解與解析解基本一致，從而驗(yàn)證了本文所述算法及所編制的程序的正確性。

表1 模型1參數(shù)

圖4 2.5維有限元數(shù)值解與一維數(shù)字濾波解析解對(duì)比

3.2 模型2

模型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)范圍與異常體的埋藏位置基本吻合。在計(jì)算的頻率范圍內(nèi)，復(fù)電阻率的相位及虛分量均為負(fù)值，偶極裝置下的復(fù)電阻率振幅和相位均呈“八”字型異常特征，特別值得注意的是復(fù)電阻率的虛分量，其幅值隨著頻率的增大而減小，且遠(yuǎn)小于實(shí)分量。

圖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

3.3 模型3

為了進(jìn)一步驗(yàn)證本文算法的有效性，設(shè)計(jì)一個(gè)極化半空間中有雙異常體模型。異常體大小均為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的不同測(cè)量裝置下的復(fù)電阻率法2.5D有限元數(shù)值模擬結(jié)果。由圖9～12可看到，復(fù)電阻率的各分量都很好地反映出異常體的存在，且異常響應(yīng)范圍與異常體的埋藏位置基本吻合。在計(jì)算的頻率范圍內(nèi)，4種測(cè)量裝置下的視復(fù)電阻率的虛分量均為負(fù)值；對(duì)稱四極裝置下，復(fù)電阻率的實(shí)部和虛部呈現(xiàn)雙“柱體”型異常特征，如圖9所示；如圖10所示，偶極裝置下的復(fù)電阻率實(shí)部和虛部均呈雙“八”字型異常特征；從圖11可知，三極測(cè)量裝置，復(fù)電阻率的實(shí)部和虛部的異常特征與對(duì)稱四極類似，不過虛部的右側(cè)部分的值大于左側(cè)；如圖12可知，二極裝置下，可以看到復(fù)電阻率的異常位置與模型非常吻合，其縱向分辨率較于三極和對(duì)稱四極測(cè)量裝置高。同樣，我們可以觀察到，在4種測(cè)量裝置下，復(fù)電阻率的實(shí)部隨著頻率的增大，產(chǎn)生的變化不大；不同于實(shí)部分量，復(fù)電阻率的虛分量，其幅值隨著頻率的增大而減小，且遠(yuǎn)小于實(shí)分量。

圖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

4 結(jié)論

本文實(shí)現(xiàn)了復(fù)電阻率法2.5D有限元數(shù)值模擬，采用異常復(fù)電法，避免了源的奇異性的影響，提高了數(shù)值解的精度；引入Fractal 模型對(duì)頻率域激發(fā)極化響應(yīng)進(jìn)行研究，可以為更好解釋激發(fā)極化效應(yīng)提供一種新的手段；對(duì)復(fù)電導(dǎo)率及復(fù)電位進(jìn)行線性插值，可以模擬更為真實(shí)的野外地質(zhì)情況；分析了所設(shè)計(jì)的地電模型的復(fù)電祖率異常響應(yīng)特征，可為野外實(shí)際勘探提供理論依據(jù)。數(shù)值模擬結(jié)果驗(yàn)證了本文基于Fractal 模型的CR2.5D有限元數(shù)值模擬算法的正確性和精確性，表明引入Fractal 模型研究復(fù)電導(dǎo)率呈線性變化的CR響應(yīng)特征是可行的、有效的。