偽隨機(jī)激電法的有限單元法數(shù)值模擬研究

楊曉弘，曾凡秋

(湖南省有色地質(zhì)勘查研究院，長(zhǎng)沙 410015)

0 前言

長(zhǎng)期以來(lái)電法勘探，包括激電法和電磁法在發(fā)現(xiàn)銅、鉛、鋅等金屬礦產(chǎn)方面發(fā)揮了巨大的作用。然而，區(qū)分礦與非礦多年來(lái)一直是困擾礦產(chǎn)地球物理工作者的難題。地球物理學(xué)家何繼善院士經(jīng)過(guò)多年的理論研究與生產(chǎn)實(shí)踐，提出了偽隨機(jī)激電法理論。偽隨機(jī)激電法的場(chǎng)源是采用2n系列偽隨機(jī)信號(hào)波形電流作為電法的激勵(lì)場(chǎng)源，接收經(jīng)過(guò)介質(zhì)響應(yīng)后的偽隨機(jī)多頻電場(chǎng)信號(hào)和磁場(chǎng)信號(hào)，經(jīng)處理、分析后提取地質(zhì)信息，而達(dá)到電法勘探的目的。該激勵(lì)場(chǎng)源可以應(yīng)用于CSAMT、MT、IP、SIP、CR等。偽隨機(jī)信號(hào)具有易于大功率發(fā)送的優(yōu)點(diǎn)，與方波相比，其帶寬是可控的。在實(shí)踐測(cè)量中可以通過(guò)控制信號(hào)帶寬來(lái)控制測(cè)量效率，從而完全改變了人工場(chǎng)源頻率域電法勘探效率低這一缺點(diǎn)[1]。

不論是哪種地球物理勘探方法，在資料的解釋方面都存在問(wèn)題的多解性，由于物理模擬比較昂貴和數(shù)字計(jì)算又不太適合物性復(fù)雜分布的地質(zhì)情況，所以數(shù)值模擬成為大多數(shù)地球物理工作者有力的工具[2-3]。作者試圖通過(guò)對(duì)偽隨機(jī)激電法進(jìn)行正演數(shù)值模擬，來(lái)研究偽隨機(jī)激電法的一些特點(diǎn)和屬性，在偽隨機(jī)激電法數(shù)值模擬方面進(jìn)行一些有益的嘗試。

1 偽隨機(jī)激電法基本原理

2n系列偽隨機(jī)、多頻信號(hào)電磁法理論是由中國(guó)工程院院士何繼善所發(fā)明和命名，激發(fā)信號(hào)隨著n的不同，在時(shí)間域上具有不同的波形。但它們的頻譜有一個(gè)共同的特點(diǎn)，就是它們的頻譜在頻率為2kω0(k=0, 1, 2, …,n-2,n-1；ω0為基波頻率)共n個(gè)頻率上大小基本相等，在這些頻率上的功率之和占了信號(hào)總平均功率的大部分，其余諧波含有的功率只占信號(hào)總平均功率的小部分，2n系列偽隨機(jī)信號(hào)的名稱也因此得來(lái)。

用p(2,n,t)表示偽隨機(jī)n頻波，在周期[0,T)內(nèi)可以表示為：

(1)

由于該方法場(chǎng)源包含多個(gè)頻率，而且能同時(shí)接收和測(cè)量這些頻率信號(hào)，所以該方法具有一次可以測(cè)量多個(gè)頻率信號(hào)，減弱了由于干擾和電流變化等的影響，提高了測(cè)量精度，特別是相對(duì)精度，可以觀測(cè)弱信號(hào)，發(fā)送電流小，儀器輕便，信息量豐富，可自動(dòng)去除感應(yīng)耦合等特點(diǎn)[1]。

2 數(shù)值模擬的基本原理與方法

2．1 數(shù)值模擬的基本原理

偽隨機(jī)激電法也是頻率域激電方法，對(duì)其進(jìn)行數(shù)值模擬的出發(fā)點(diǎn)是Pelton定義的Cole-Cole模型[4-6]。Cole-Cole模型可以用式(2)來(lái)表示

(2)

其中ρ0為頻率為零時(shí)巖、礦石的電阻率；m為充電率；c為頻率相關(guān)系數(shù)；τ為時(shí)間常數(shù)。當(dāng)激發(fā)電流的頻率較低的時(shí)候可以忽略電磁效應(yīng)的情況，這時(shí)可以用Cole-Cole模型得到的復(fù)電阻率或復(fù)電位來(lái)替換大地模型的地質(zhì)體，通過(guò)求解大地模型得到不同頻率下地表節(jié)點(diǎn)上的復(fù)電位值，就可以獲得偽隨機(jī)激電法數(shù)值模擬的參數(shù)[7-8]。

當(dāng)激發(fā)電流頻率較低的時(shí)候，可以在忽略電磁效應(yīng)的情況下用穩(wěn)定電流場(chǎng)的邊值問(wèn)題和變分問(wèn)題對(duì)以上問(wèn)題進(jìn)行研究。三維構(gòu)造中雙點(diǎn)電源電場(chǎng)模型的電位邊值問(wèn)題為[9-10]

(3)

式中u為電位；σ為三維地下空間介質(zhì)的電導(dǎo)率；I為電流強(qiáng)度；Ω為模型區(qū)域；n為邊界的外法向方向；Γs為模型的地面邊界；γA、γB為測(cè)點(diǎn)到電源點(diǎn)A、B的距離；?！逓闊o(wú)窮遠(yuǎn)邊界。

與上述邊值問(wèn)題(3)等價(jià)的變分問(wèn)題為

(4)

2．2 有限單元法

2.2.1 單元剖分

為方便程序的編制，用正六面體單元對(duì)大地模型進(jìn)行剖分(圖1，圖2)。

圖1 模型剖分示意圖Fig．1 Sketch of model division

圖2 剖分單元Fig．2 Element of the model

2.2.2 單元積分

經(jīng)過(guò)線性插值后，對(duì)變分方程的單元積分變?yōu)?/p>

(5)

{Ne}T{φe}]dV

(6)

其中N為形函數(shù)。

式中ξi、ηi、ζi分別是等參單元中點(diǎn)i(i=1,…，8)的坐標(biāo)。

2.2.3 總體合成及求變分

(7)

令δF(φ)=0,

則有

[K′]{φ}-{φ}=0,

即

{K}{φ}=0，其中K=[K′]-。

2.2.4 解方程及計(jì)算視幅頻率

解形如{K}{φ}=0的方程就得到模型中各個(gè)節(jié)點(diǎn)上的復(fù)電位值，偽隨機(jī)激電法的基本參數(shù)都可以通過(guò)相應(yīng)的定義計(jì)算得到。

3 模型計(jì)算

假設(shè)有一個(gè)大小為500 m×100 m×100 m三維均勻半空間大地模型，模型的剖面圖如圖3所示。

圖3 模型斷面示意圖Fig．3 Sketch of the earth model

采用偽隨機(jī)信號(hào)激發(fā)時(shí)，首先研究極化異常體的時(shí)間常數(shù)τ對(duì)幅頻率的影響情況。當(dāng)測(cè)量電極在極化異常體的上方地面中心位置固定的時(shí)候，選用的激勵(lì)電流頻率范為 0.000 1 Hz～10 000 Hz，極化異常體的m2=0.8、c2=0.25；大地介質(zhì)中的m1=0.04、c1=0.25、τ1=1 s，采用不同的τ2值時(shí)就可得到對(duì)應(yīng)的幅頻率曲線圖(圖4)。

圖4 不同時(shí)間常數(shù)對(duì)應(yīng)的幅頻率曲線(m=0.8,c=0.25)Fig．4 Curves of Fs in different time constant(m=0.8,c=0.25)

從圖4可以看出，幅頻率在某個(gè)頻率上取得其極值，而隨著頻率的無(wú)限增大或無(wú)限減少，幅頻率趨向于“0”值，即異常體來(lái)不及激化或激化已經(jīng)趨于穩(wěn)定。隨著時(shí)間常數(shù)τ值的增大，幅頻率極值頻率向低頻方向移動(dòng)。幅頻率的極值頻率雖然隨著時(shí)間常數(shù)的不同而不同，但幅頻率的極值卻相同，故時(shí)間常數(shù)不會(huì)影響幅頻率的極值大小。時(shí)間常數(shù)也不影響幅頻率曲線的形態(tài)，而只影響其位置。通過(guò)對(duì)比野外實(shí)測(cè)幅頻率曲線，其極值頻率能在一定程度上反映勘探區(qū)內(nèi)的綜合時(shí)間常數(shù)的大小，而大部分的金屬硫化物的時(shí)間常數(shù)都很大，無(wú)激電效應(yīng)或激電效應(yīng)很小的圍巖的時(shí)間常數(shù)一般較小，故通過(guò)比較幅頻率的極值頻率大小能對(duì)異常源作一定的定性解釋。

圖5是m對(duì)幅頻率值的影響情況。當(dāng)充電率從0.1變化到0.9時(shí)就能得到不同m2值所對(duì)應(yīng)的幅頻率Fs的曲線圖。從該曲線圖中可以看出，幅頻率隨著充電率m2的增大，其極值也變大，極值頻率緩慢地向高頻移動(dòng)。這與實(shí)際情況相符，充電率與極化率、幅頻率在發(fā)現(xiàn)異常上是等效的，且其值也是正相關(guān)。充電率的大小影響幅頻率曲線的形態(tài)，而對(duì)于其位置影響很小。

圖5 不同充電率m對(duì)應(yīng)的幅頻率曲線(τ=50 s,c=0.25)Fig．5 Curves of Fs in different charge rate m(τ=50 s,c=0.25)

圖6是不同的頻率相關(guān)系數(shù)c2所對(duì)應(yīng)的幅頻率曲線圖。

圖6 不同頻率相關(guān)系數(shù)對(duì)應(yīng)的幅頻率曲線(m=0.8,τ=50 s)Fig．6 Curves of Fs in different correlation coefficient(m=0.8,τ=50 s)

從圖6可知，隨著頻率相關(guān)系數(shù)c2的增大，幅頻率極值也相應(yīng)變大，但在極值兩翼上，幅頻率表現(xiàn)出無(wú)規(guī)律的變化，幅頻率的極值頻率基本上與頻率相關(guān)系數(shù)無(wú)關(guān)，故頻率相關(guān)系數(shù)與幅頻率除了在極值頻率附近外，在其他頻率段上相關(guān)性不大。

4 結(jié)論

在忽略電磁效應(yīng)的情況下，通過(guò)用Cole-Cole模型的相應(yīng)替換大地模型中的參數(shù)，與用有限單元法模擬偽隨機(jī)激電法中的模型相應(yīng)。通過(guò)分析偽隨機(jī)激電法正演模擬的結(jié)果可知，模擬結(jié)果正確可靠，表明運(yùn)用有限單元法模擬偽隨機(jī)激電法中的參數(shù)是正確和適用的，為偽隨機(jī)激電法的數(shù)值模擬進(jìn)行了一些有益的嘗試。

參考文獻(xiàn)：

[1] 何繼善.雙頻激電法[M].北京：高等教育出版社，2006.

[2] 徐世浙.地球物理中的有限單元法[M].北京：科學(xué)出版社，1994.

[3] 黃俊革.三維電阻率/極化率有限元正演模擬與反演成像[D].長(zhǎng)沙：中南大學(xué)，2003.

[4] 傅良魁. 電法勘探教程[M].北京：地質(zhì)出版社，1983.

[5] 傅良魁. 應(yīng)用地球物理教程[M]．北京：地質(zhì)出版社，1991.

[6] 李金銘. 電法勘探方法發(fā)展概況[J]. 物探與化探, 1996,20(4):93-95.

[7] 羅延鐘,孟永良. 關(guān)于用有限單元法計(jì)算二維構(gòu)造點(diǎn)電源場(chǎng)的幾個(gè)問(wèn)題[J].地球物理學(xué)報(bào), 1986,29(1): 32-35.

[8] 周熙襄，鐘本善. 電法勘探數(shù)值模擬技術(shù)[M]. 成都： 四川科學(xué)技術(shù)出版社, 1986.

[9] COGGON J H.Electromagnetic and electrical modeling by the finite-element method [J].Geophysics, 1971,36(1):132-155.

[10] 李大潛. 有限元素法在電法測(cè)井中的應(yīng)用[M]. 北京：石油工業(yè)出版社，1980.