基于數(shù)字濾波算法的時(shí)域激電場(chǎng)Cole-Cole 模型參數(shù)最優(yōu)化計(jì)算

劉瑞澤, 孟慶鑫, 喬圣揚(yáng)

河北地質(zhì)大學(xué), 河北 石家莊 050031

0 引言

時(shí)間域譜激發(fā)極化法是通過(guò)測(cè)量人工場(chǎng)源引起的大地視時(shí)變電阻率特征, 尋找電性異常體, 并根據(jù)巖、 礦石時(shí)譜參數(shù), 評(píng)價(jià)電性異常體的一種電法勘探方法[1]。 時(shí)譜激電法能獲得反映地下不同時(shí)間地質(zhì)信息的電性譜[2], 擁有一定的應(yīng)用前景[1]。 近年來(lái), 國(guó)內(nèi)外研究者將譜激電法應(yīng)用于金屬礦物分析[3]、 不同流體飽和度巖層[4]、 巖樣滲漏率計(jì)算[5]、 工程地質(zhì)勘探中垃圾填埋場(chǎng)污染范圍[6]、 監(jiān)測(cè)地球化學(xué)場(chǎng)[7]等方面; 同時(shí), 隨著科技進(jìn)步、 儀器設(shè)備的改進(jìn), 時(shí)譜激電法在數(shù)據(jù)采集[8]、 測(cè)量高信噪比[9]等技術(shù)方面得到改善, 以及有關(guān)時(shí)譜激電新的計(jì)算方法的引入[10]。 綜上所述, 科研人員不僅拓寬時(shí)譜激電法的應(yīng)用領(lǐng)域, 而且完善了理論方法和技術(shù)手段, 這些都使時(shí)譜激電法得以發(fā)展。 上述成果仍需要推廣至其他領(lǐng)域, 也需要將其他領(lǐng)域的一些好的理論成果和技術(shù)引入到時(shí)譜激電法中; 而且, 以往主要研究的是視時(shí)變電阻率, 極化率的大小[11]等參數(shù), 在處理一些復(fù)雜介質(zhì)問(wèn)題時(shí), 這些參數(shù)反映出的電性差異不明顯, 甚至不能區(qū)分電性異常體, 而且解釋方法單一[12]。 隨著研究的不斷深入, 研究者發(fā)現(xiàn)譜激電特征參數(shù)能夠區(qū)分復(fù)雜介質(zhì)的電性差異, 這種應(yīng)用多參數(shù)聯(lián)合解釋電性差異的方法日趨普遍[13]。 為了更好的對(duì)野外觀測(cè)數(shù)據(jù)進(jìn)行解釋, 譜特征參數(shù)的最優(yōu)化計(jì)算是有必要的。 故本文將在數(shù)字濾波算法的基礎(chǔ)上, 完成時(shí)譜激電法譜特征參數(shù)的最優(yōu)化計(jì)算; 通過(guò)理論算例, 驗(yàn)證反算出的譜特征參數(shù)的精度; 進(jìn)行不同充電時(shí)間的譜特征參數(shù)最優(yōu)化計(jì)算, 分析不同充電時(shí)間對(duì)譜特征參數(shù)精度的影響。

1 時(shí)譜激電正演

自激發(fā)極化法廣泛應(yīng)用于實(shí)際生產(chǎn)后, 研究者在激電譜特性區(qū)分礦化類(lèi)型[14]、 利用在復(fù)平面中不同的復(fù)電阻率頻譜類(lèi)型定性評(píng)價(jià)異常源性質(zhì)[15]、 利用激電衰減時(shí)和衰減度等特性找地下水[16]等方面取得重要成果; 此外, 也有研究者在定量描述激電譜特性方面做過(guò)許多研究, 其中柯?tīng)?柯?tīng)柲Ｐ偷奶岢鯷17][18]提供了對(duì)激電譜定量描述的手段。 本文基于柯?tīng)?柯?tīng)柲Ｐ秃蛿?shù)字濾波算法來(lái)研究譜特征參數(shù)。

1.1 柯?tīng)?柯?tīng)柲Ｐ?/h3>

當(dāng)通過(guò)巖、 礦石的供電電流密度不大時(shí), 巖、 礦石的激電效應(yīng)具有線(xiàn)性規(guī)律, 并且在時(shí)間范圍內(nèi)可以用不變的參數(shù)來(lái)描述, 是一個(gè)線(xiàn)性不變系統(tǒng)。 Madden和Cantwell、 Pelton 等[17]借 用 了 一 種 與 Cole 和Cole[19]提出的張弛模型類(lèi)似的傳遞函數(shù), 來(lái)描述復(fù)電阻率頻譜的性質(zhì), 將等效電路模型稱(chēng)之為柯?tīng)?柯?tīng)柲Ｐ? 函數(shù)稱(chēng)之柯?tīng)?柯?tīng)栕杩贡磉_(dá)式。

研究者提出的柯?tīng)?柯?tīng)栕杩贡磉_(dá)式為[17]:

式中:Z(0)、m、τ、c分別為頻率為零時(shí)的阻抗、極化率、 時(shí)間常數(shù)和頻率相關(guān)系數(shù), 并將其稱(chēng)為柯?tīng)?柯?tīng)枀?shù)。Z(0) 的單位為Ω,τ的單位為s,m和c無(wú)量綱。

根據(jù)上式不難得到在頻率域中的復(fù)電阻率的柯?tīng)?柯?tīng)柋磉_(dá)式[1]:

在這一復(fù)電阻率的柯?tīng)?柯?tīng)柋磉_(dá)式中,ρ0為頻率為零時(shí)的電阻率, 其余參數(shù)與阻抗表達(dá)式中的相同,也將這些參數(shù)稱(chēng)之為柯?tīng)?柯?tīng)栴l域譜特征參數(shù)。ρ0的單位為Ω,τ的單位為s,m和c無(wú)量綱。

1.2 基于數(shù)字濾波算法柯?tīng)?柯?tīng)柲Ｐ偷念l-時(shí)轉(zhuǎn)換

研究者常將柯?tīng)?柯?tīng)柲Ｐ蛻?yīng)用于頻域激發(fā)極化法中, 并用柯?tīng)?柯?tīng)栕V特征參數(shù)描述激電譜特性。研究者的研究證實(shí)了譜特征參數(shù)可以很好的定量描述多數(shù)巖、 礦石的激電譜特性, 并指出可以根據(jù)巖、 礦石的譜特征參數(shù)尋找激電異常體。 本文采用Guptasarma[20]研發(fā)的一種計(jì)算瞬變響應(yīng)的數(shù)字濾波算法, 對(duì)柯?tīng)?柯?tīng)柲Ｐ瓦M(jìn)行“頻-時(shí)” 轉(zhuǎn)換, 以此求取時(shí)間域中的電阻率柯?tīng)?柯?tīng)柲Ｐ汀?/p>

為了計(jì)算復(fù)電阻率的柯?tīng)?柯?tīng)柋磉_(dá)式ρ(ω) 在t時(shí)刻的響應(yīng)ρ(t) , 首先確定頻域復(fù)電阻率的虛分量、實(shí)分量。 研究者推導(dǎo)出的含有虛分量和實(shí)分量的復(fù)電阻率表達(dá)式為[1,17]:

本數(shù)字濾波算法只需要頻域復(fù)電阻率的實(shí)分量,即為:

再根據(jù)數(shù)字濾波的離散式[20]:

Guptasarma 提出的數(shù)字濾波算法要用到21 個(gè)濾波系數(shù), 然后通過(guò)累加的方式來(lái)求取電阻率或電位差。

轉(zhuǎn)換成時(shí)間域電阻率或電位差的公式后, 要通過(guò)建模試算, 觀察不同譜特征參數(shù)的時(shí)間域譜激電法電阻率或電位差特性曲線(xiàn)變化特征[21]。

根據(jù)前人的研究文獻(xiàn)資料[11,21], 設(shè)置多組不同的譜特征參數(shù), 依據(jù)不同物質(zhì)的充電飽和特性和常用觀測(cè)時(shí)道, 以時(shí)間作為自變量, 電阻率作為因變量, 來(lái)觀測(cè)不同譜特征參數(shù)引起曲線(xiàn)的變化特征。 通過(guò)編程的方式, 來(lái)繪制不同譜特征參數(shù)對(duì)應(yīng)的電阻率和電位差隨時(shí)間變化的曲線(xiàn)圖。

圖1 為時(shí)間常數(shù)和極化率取定值, 時(shí)間相關(guān)系數(shù)在變化的時(shí)變電阻率曲線(xiàn)圖。 在時(shí)間軸為對(duì)數(shù)比例尺、 時(shí)變電阻率為算數(shù)比例尺的坐標(biāo)軸中, 時(shí)變電阻率與時(shí)間在時(shí)間常數(shù)取值附近有中心對(duì)稱(chēng)的趨勢(shì)。 并且在遠(yuǎn)離中心對(duì)稱(chēng)點(diǎn)處, 隨著頻率相關(guān)系數(shù)數(shù)值的增加, 時(shí)變電阻率隨時(shí)間變化的曲線(xiàn)斜率逐漸減小, 而在中心對(duì)稱(chēng)點(diǎn)附近, 則隨著頻率相關(guān)系數(shù)數(shù)值的增加, 時(shí)變電阻率隨時(shí)間變化曲線(xiàn)的斜率增大。

圖1 時(shí)間常數(shù)和極化率取定值, 不同的時(shí)間相關(guān)系數(shù)對(duì)應(yīng)的時(shí)變電阻率曲線(xiàn)圖Fig.1 The time constant and polarizability are fixed values,and the time-varying resistivity curves corresponding to different time correlation coefficients

圖2 為極化率和時(shí)間相關(guān)系數(shù)一定, 時(shí)間常數(shù)在變化的時(shí)變電阻率曲線(xiàn)圖。 在時(shí)間軸為對(duì)數(shù)比例尺、時(shí)變電阻率為算數(shù)比例尺的坐標(biāo)軸中, 不同時(shí)間常數(shù)所對(duì)應(yīng)的時(shí)變電阻率與時(shí)間的關(guān)系曲線(xiàn)的變化趨勢(shì)是一致的, 不同時(shí)間常數(shù)所對(duì)應(yīng)的時(shí)變電阻率隨時(shí)間的變化曲線(xiàn)可以通過(guò)左右平移得到。

圖2 極化率和時(shí)間相關(guān)系數(shù)取定值, 不同的時(shí)間常數(shù)對(duì)應(yīng)的時(shí)變電阻率曲線(xiàn)圖Fig.2 The polarizability and time correlation coefficient are taken a fixed value, and the time-varying resistivity curve corresponding to different time constants

圖3 為時(shí)間常數(shù)和時(shí)間相關(guān)系數(shù)一定, 極化率在變化的時(shí)變電阻率曲線(xiàn)圖。 在時(shí)間軸為對(duì)數(shù)比例尺、時(shí)變電阻率為算數(shù)比例尺的坐標(biāo)軸中, 時(shí)變電阻率時(shí)譜有如下特征: 結(jié)果恒為正值, 其數(shù)值隨時(shí)間增大單調(diào)地增加; 在長(zhǎng)時(shí)間供電時(shí), 時(shí)變電阻率趨于ρ(∞), 在短時(shí)間供電時(shí), 時(shí)變電阻率趨于ρ(∞)(1-m) 值。

圖3 時(shí)間常數(shù)和極化率取定值, 不同的極化率對(duì)應(yīng)的時(shí)變電阻率曲線(xiàn)圖Fig.3 The time constant and the polarizability are fixed values, and the time-varying resistivity curves corresponding to different polarizabilities

圖4 為不同譜特征參數(shù)對(duì)應(yīng)的充放電曲線(xiàn)圖。 在橫坐標(biāo)為充電時(shí)間, 縱坐標(biāo)為電位差的充放電曲線(xiàn)中, 不同的譜特征參數(shù)對(duì)應(yīng)著不同的充放電曲線(xiàn)。 不同的譜特征參數(shù)對(duì)應(yīng)著不同的極化異?，F(xiàn)象, 可以通過(guò)譜特征參數(shù)來(lái)區(qū)分極化異常。

圖4 不同譜特征參數(shù)對(duì)應(yīng)的充放電曲線(xiàn)圖Fig.4 Charge and discharge curves corresponding to different spectral characteristic parameters

2 時(shí)間域激電場(chǎng)視譜特征參數(shù)最優(yōu)化計(jì)算

根據(jù)上文中提到的通過(guò)數(shù)字濾波算法, 將頻域復(fù)變電阻率表達(dá)式進(jìn)行“頻-時(shí)” 轉(zhuǎn)換后, 得到時(shí)間域電阻率表達(dá)式, 利用編程軟件, 完成時(shí)譜特征參數(shù)最優(yōu)化計(jì)算。 其中, 根據(jù)之前的研究[10], 選擇常用的理論算例(主要為時(shí)間常數(shù)和相關(guān)系數(shù)的選擇), 來(lái)進(jìn)行譜特征參數(shù)的最優(yōu)化計(jì)算, 并通過(guò)理論算例來(lái)驗(yàn)證該最優(yōu)化算法計(jì)算的精度。

2.1 目標(biāo)函數(shù)的構(gòu)建

本文的目標(biāo)函數(shù)采用最小二乘原理來(lái)進(jìn)行構(gòu)建。首先譜特征參數(shù)最優(yōu)化計(jì)算目標(biāo)函數(shù)中的元素(采用了比值的形式, 目的是只對(duì)τ,c進(jìn)行最優(yōu)化計(jì)算, 減少參數(shù), 降低多解性) 為:

式中:mona,moffa分別表示充電和放電時(shí)的視充電率, 而視充電率可通過(guò)一次場(chǎng)、 二次場(chǎng)電位差的比值獲得。

譜特征參數(shù)最優(yōu)化計(jì)算的目標(biāo)函數(shù)為:

2.2 視譜特征參數(shù)最優(yōu)化計(jì)算的算法步驟

步1: 給出τ1,c1, 0 ≤ε≤1;

k =1。

步2: 如果 ▽?duì)?τk)2+ ▽?duì)?ck)2≤ε, 則停;

dk=(-▽?duì)?τk),-▽?duì)?ck))。

步3: 如果k=1, 則利用不精確搜索方式求a1,否則,*(dk),k: =k+1, 轉(zhuǎn)步2。

進(jìn)行視譜特征參數(shù)最優(yōu)解計(jì)算時(shí), 首先設(shè)置初始時(shí)間常數(shù)τ1和相關(guān)系數(shù)值c1, 以及擬合精度, 將初始值和實(shí)測(cè)數(shù)據(jù)代入目標(biāo)函數(shù), 然后根據(jù)上述步驟對(duì)目標(biāo)函數(shù)迭代求解, 最終滿(mǎn)足擬合精度時(shí)的時(shí)間常數(shù)和相關(guān)系數(shù)就是本文要求的視譜特征參數(shù)。 根據(jù)上述算法步驟, 進(jìn)行計(jì)算機(jī)編程。

步4: (τk+1,ck+1) = (τk,ck) +ak

2.3 理論算例驗(yàn)證

根據(jù)之前研究者使用柯?tīng)?柯?tīng)柲Ｐ瓦M(jìn)行激發(fā)極化法仿真研究時(shí)[18][10], 得出了25 組CC 參數(shù)模型(c范圍從0.1 到0.6,τ范圍從0.001 到100 s), 每組譜特征參數(shù)模型對(duì)應(yīng)于充電過(guò)程中4 個(gè)時(shí)間點(diǎn)(0.92、0.5、 0.3 和0.2 s) 的測(cè)量數(shù)據(jù)。 依據(jù)理論算例, 基于2.1、 2.2 構(gòu)造的目標(biāo)函數(shù)和算法步驟通過(guò)編程軟件, 最優(yōu)化計(jì)算出的譜特征參數(shù)結(jié)果如表1。

表1 理論算例模型譜特征參數(shù)與最優(yōu)化算法計(jì)算出的譜特征參數(shù)結(jié)果的對(duì)比Table 1 Comparison of the spectral characteristic parameters of the theoretical calculation example model and the results of the spectral characteristic parameters calculated by the optimization algorithm

續(xù)表1

表1 的左側(cè)為具體譜特征參數(shù)(τ、c) 試算模型數(shù)值, 表1 的右側(cè)是25 組理論算例模型的譜特征參數(shù)最優(yōu)化計(jì)算結(jié)果。 從視譜特征參數(shù)最優(yōu)化算法計(jì)算出的譜特征參數(shù)τ、c的結(jié)果來(lái)看, 反算的譜特征參數(shù)于所設(shè)的理論算例高度一致。 通過(guò)計(jì)算反算譜特征參數(shù)與理論的譜特征參數(shù)的相對(duì)誤差發(fā)現(xiàn), 它們的相對(duì)誤差都小于1‰。 這說(shuō)明譜特征參數(shù)的最優(yōu)化計(jì)算結(jié)果與理論值有很高的擬合精度, 該算法對(duì)理論算例來(lái)說(shuō)是可行的。

3 結(jié)論

(1) 基于數(shù)字濾波算法對(duì)柯?tīng)?柯?tīng)柲Ｐ瓦M(jìn)行頻時(shí)變換得到了時(shí)域表達(dá)式, 根據(jù)該表達(dá)式繪制了不同時(shí)刻和譜特征參數(shù)條件的時(shí)譜曲線(xiàn), 分析可知, 不同譜特征參數(shù)對(duì)應(yīng)的時(shí)譜曲線(xiàn)有明顯的不同, 可以根據(jù)譜特征參數(shù)所反應(yīng)的時(shí)變特性來(lái)區(qū)分極化異常體。

(2) 本文構(gòu)建了比值形式的目標(biāo)函數(shù), 該目標(biāo)函數(shù)需要反算的視譜特征參數(shù)少, 在當(dāng)前主流儀器觀測(cè)時(shí)道較少的情況下, 具有降低反算結(jié)果的多解性。

(3) 通過(guò)理論算例試算, 證明本文編寫(xiě)的視譜特征參數(shù)最優(yōu)化算法是可行的, 而且反算的結(jié)果精確度較高。

(4) 本文只驗(yàn)證了試?yán)碚撍憷淖V特征參數(shù), 之后的研究, 要將該方法應(yīng)用于簡(jiǎn)單地質(zhì)模型進(jìn)行驗(yàn)證,最終目的是將該方法應(yīng)用于野外實(shí)測(cè)數(shù)據(jù)的處理, 以計(jì)算的譜特征參數(shù)來(lái)尋找異常區(qū)域。 同時(shí)也發(fā)現(xiàn), 數(shù)字濾波算法轉(zhuǎn)換后的柯?tīng)?柯?tīng)柲Ｐ陀?jì)算量比較大, 接下來(lái)的研究將降低此算法的計(jì)算量大的問(wèn)題。