偽隨機激電法的有限單元法數值模擬研究

楊曉弘，曾凡秋

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

0 前言

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

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

1 偽隨機激電法基本原理

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

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

(1)

由于該方法場源包含多個頻率，而且能同時接收和測量這些頻率信號，所以該方法具有一次可以測量多個頻率信號，減弱了由于干擾和電流變化等的影響，提高了測量精度，特別是相對精度，可以觀測弱信號，發送電流小，儀器輕便，信息量豐富，可自動去除感應耦合等特點[1]。

2 數值模擬的基本原理與方法

2．1 數值模擬的基本原理

偽隨機激電法也是頻率域激電方法，對其進行數值模擬的出發點是Pelton定義的Cole-Cole模型[4-6]。Cole-Cole模型可以用式(2)來表示

(2)

其中ρ0為頻率為零時巖、礦石的電阻率；m為充電率；c為頻率相關系數；τ為時間常數。當激發電流的頻率較低的時候可以忽略電磁效應的情況，這時可以用Cole-Cole模型得到的復電阻率或復電位來替換大地模型的地質體，通過求解大地模型得到不同頻率下地表節點上的復電位值，就可以獲得偽隨機激電法數值模擬的參數[7-8]。

當激發電流頻率較低的時候，可以在忽略電磁效應的情況下用穩定電流場的邊值問題和變分問題對以上問題進行研究。三維構造中雙點電源電場模型的電位邊值問題為[9-10]

(3)

式中u為電位；σ為三維地下空間介質的電導率；I為電流強度；Ω為模型區域；n為邊界的外法向方向；Γs為模型的地面邊界；γA、γB為測點到電源點A、B的距離；?！逓闊o窮遠邊界。

與上述邊值問題(3)等價的變分問題為

(4)

2．2 有限單元法

2.2.1 單元剖分

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

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

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

2.2.2 單元積分

經過線性插值后，對變分方程的單元積分變為

(5)

{Ne}T{φe}]dV

(6)

其中N為形函數。

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

2.2.3 總體合成及求變分

(7)

令δF(φ)=0,

則有

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

即

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

2.2.4 解方程及計算視幅頻率

解形如{K}{φ}=0的方程就得到模型中各個節點上的復電位值，偽隨機激電法的基本參數都可以通過相應的定義計算得到。

3 模型計算

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

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

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

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

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

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

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

圖6是不同的頻率相關系數c2所對應的幅頻率曲線圖。

圖6 不同頻率相關系數對應的幅頻率曲線(m=0.8,τ=50 s)Fig．6 Curves of Fs in different correlation coefficient(m=0.8,τ=50 s)

從圖6可知，隨著頻率相關系數c2的增大，幅頻率極值也相應變大，但在極值兩翼上，幅頻率表現出無規律的變化，幅頻率的極值頻率基本上與頻率相關系數無關，故頻率相關系數與幅頻率除了在極值頻率附近外，在其他頻率段上相關性不大。

4 結論

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

參考文獻：

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

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

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

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

[5] 傅良魁. 應用地球物理教程[M]．北京：地質出版社，1991.

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

[7] 羅延鐘,孟永良. 關于用有限單元法計算二維構造點電源場的幾個問題[J].地球物理學報, 1986,29(1): 32-35.

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

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

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