林小穩(wěn), 柯式鎮(zhèn), 賀秋麗, 許巍, 舒心, 張琪

(1.中國(guó)石油大學(xué)(北京)地球物理與信息工程學(xué)院, 北京 102249; 2.中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249; 3.中國(guó)石油集團(tuán)測(cè)井有限公司技術(shù)中心, 陜西 西安 710077)

0 引 言

激發(fā)極化法(IP)在地球物理勘探領(lǐng)域已經(jīng)被廣泛應(yīng)用。復(fù)電阻率法(CR)又稱(chēng)頻譜激電法(SIP),是20世紀(jì)70年代發(fā)展起來(lái)的激電分支方法[1]。在對(duì)巖石復(fù)電阻率研究中通過(guò)等效電路建立模型模擬巖石的頻散特性,使用最為廣泛的是Cole-Cole模型。利用實(shí)驗(yàn)室?guī)r心測(cè)量數(shù)據(jù),選擇合適的算法反演出Cole-Cole模型中的電頻譜參數(shù),進(jìn)而對(duì)解釋模型進(jìn)行標(biāo)定是復(fù)電阻率測(cè)井資料解釋的核心問(wèn)題。國(guó)內(nèi)外學(xué)者對(duì)模型參數(shù)的反演作了許多研究,Jaggar等于1988年提出了基于馬奎特算法的反演[2]。2001年,J.Xiang等提出了基于多重最小二乘估計(jì)的直接反演算法,該方法不需要任何中間變量的介入,也不用對(duì)參數(shù)初始值進(jìn)行估計(jì),可反演得到一階Cole-Cole模型中的4個(gè)參數(shù)[3]。2005年曹中林引入遺傳算法對(duì)一階Cole-Cole模型以及和形式的二階模型進(jìn)行了反演,并與馬奎特算法進(jìn)行了對(duì)比[4-5]。2006年童茂松等引入模擬退火算法對(duì)一階Cole-Cole模型、基于乘積形式的二階模型及Dias模型進(jìn)行了反演,并提出基于乘積形式的二階復(fù)Cole-Cole模型比一階模型更具有優(yōu)勢(shì)[6]。在實(shí)際工作中如何選擇合適的模型以及反演算法,國(guó)內(nèi)外學(xué)者對(duì)此沒(méi)有作過(guò)深入系統(tǒng)的研究,也沒(méi)有給出過(guò)明確的結(jié)論。

本文在100 Hz～10 MHz頻率范圍內(nèi)測(cè)量了多組人造砂巖樣品的復(fù)電阻率頻譜曲線并進(jìn)行了反演。在反演復(fù)電阻率頻譜時(shí),使用局部尋優(yōu)的最小二乘法以及在全局尋優(yōu)中的2種算法——單點(diǎn)搜索的模擬退火算法和并行群體搜索的遺傳算法。本文利用這3種算法結(jié)合基于乘積形式的二階Cole-Cole模型對(duì)理論模型和實(shí)測(cè)數(shù)據(jù)分別進(jìn)行了反演,綜合對(duì)比了這3種算法的優(yōu)缺點(diǎn),驗(yàn)證了二階Cole-Cole乘積模型與實(shí)測(cè)數(shù)據(jù)的適配性,為反演算法的選擇提供了指導(dǎo)。

1 復(fù)電阻率測(cè)井

復(fù)電阻率測(cè)井方法是基于地層巖石電阻率的頻散特性而發(fā)展起來(lái)的一種新的測(cè)井方法。它充分利用了巖石的頻散信息,改變了僅考慮探測(cè)空間范圍影響而忽視巖石頻散性質(zhì)的傳統(tǒng)認(rèn)識(shí),在一定程度上克服了地層水礦化度的影響,可以更加直觀地反映含油性,提高對(duì)油層的識(shí)別能力。利用復(fù)電阻率測(cè)井?dāng)?shù)據(jù)可以按照?qǐng)D1所示的流程計(jì)算含油氣飽和度[7],在這一流程中,利用巖心實(shí)驗(yàn)數(shù)據(jù)反演出Cole-Cole模型中的頻譜參數(shù)對(duì)模型進(jìn)行標(biāo)定是最為基礎(chǔ)也是最為核心的一步,反演所得參數(shù)的準(zhǔn)確性將影響模型標(biāo)定的準(zhǔn)確性,最終影響儲(chǔ)層評(píng)價(jià)的準(zhǔn)確性。

圖1 復(fù)電阻率測(cè)井解釋流程圖

2 Cole-Cole模型

Cole-Cole模型是20世紀(jì)40年代初期K.S.Cole和R.H.Cole首先提出并用來(lái)描述電解質(zhì)的介電極化效應(yīng)[8],后人在此基礎(chǔ)上不斷研究,逐漸完善了這一模型。1978年,Pelton等通過(guò)研究認(rèn)為,巖石的激電效應(yīng)與Cole-Cole模型相符合,并將其引入到地球物理勘查中,由此奠定了復(fù)電阻率法的基礎(chǔ)。Pelton認(rèn)為巖石復(fù)電阻率隨頻率的變化(即復(fù)電阻率頻譜)可表示[9]為

(1)

式中,ρ(ω)為復(fù)電阻率,Ω·m;ρ0為零頻電阻率,表征導(dǎo)電性強(qiáng)弱,Ω·m;m為極化率,表征激電效應(yīng)強(qiáng)弱;τ為時(shí)間常數(shù),表征激電效應(yīng)特性;c為頻率相關(guān)系數(shù),表征頻譜或時(shí)間譜陡緩;ω為角頻率。

一階Cole-Cole模型的復(fù)電阻率曲線實(shí)部表現(xiàn)為單斜坡,虛部為單波谷的特征,而大多數(shù)實(shí)測(cè)曲線卻表現(xiàn)為多斜坡、多波谷特征。為了更好地計(jì)算巖石復(fù)電阻率的頻散特性,Pelton等用2個(gè)Cole-Cole模型描述飽和鹽水巖石的頻散現(xiàn)象,一個(gè)表征界面極化,另一個(gè)表征鹽水偶極子極化[10]。二階Cole-Cole模型有并聯(lián)和串聯(lián)2種形式,其等效電路圖(二階)如圖2和圖3所示。本文所介紹的反演算法基于二階(N=2)并聯(lián)形式的Cole-Cole模型,且其反演得到的參數(shù)極化率mi、頻率相關(guān)系數(shù)ci要在0

(2)

式中各參數(shù)的意義與一階Cole-Cole模型中參數(shù)的意義相同。

圖2 二階串聯(lián)Cole-Cole模型等效電路

圖3 二階并聯(lián)Cole-Cole模型等效電路

3 實(shí)驗(yàn)測(cè)量系統(tǒng)

實(shí)驗(yàn)室對(duì)巖石復(fù)電阻率頻譜的測(cè)量有電極法和線圈法。本文采用如圖4所示的電極法測(cè)量[11]。

圖4 巖心復(fù)電阻率測(cè)量系統(tǒng)

實(shí)驗(yàn)測(cè)量系統(tǒng)結(jié)構(gòu)由信號(hào)發(fā)生器、功率放大器、取樣移相電路、測(cè)量電極系、巖心夾持器、鎖相放大器、數(shù)字存儲(chǔ)示波器和計(jì)算機(jī)組成。各個(gè)儀器之間采用通用并行接口總線相連,以計(jì)算機(jī)作為控制中心,控制各個(gè)儀器的工作,與控制軟件一起構(gòu)成完整的掃頻測(cè)量系統(tǒng)[12]。

4 反演算法

最小二乘法[13]讓實(shí)測(cè)數(shù)據(jù)和估計(jì)數(shù)據(jù)之間的距離平方和最小,這是最小二乘法使用的逼近原則。

模擬退火算法的思想來(lái)源于對(duì)固體退火降溫過(guò)程的模擬,即將固體加溫至充分高,再讓其充分冷卻。模擬退火算法是一種啟發(fā)式隨機(jī)搜索方法,但它在搜索策略上與傳統(tǒng)的隨機(jī)搜索方法不同,它不僅引入了適當(dāng)?shù)碾S機(jī)因素,而且還引入了物理系統(tǒng)退火過(guò)程的自然機(jī)理[14]。模擬退火算法從某一較高初溫出發(fā),伴隨溫度參數(shù)不斷下降,結(jié)合概率突跳特性在解空間中隨機(jī)尋找目標(biāo)函數(shù)的全局最優(yōu)解。

遺傳算法是一種自適應(yīng)全局優(yōu)化概率搜索算法[15]。遺傳算法模擬自然進(jìn)化中優(yōu)勝劣汰、適者生存的原理,與傳統(tǒng)的搜索策略不同,遺傳算法不是簡(jiǎn)單的單點(diǎn)搜索,而是并行群體搜索,可以同時(shí)處理群體中的多個(gè)個(gè)體[16]。

5 反演結(jié)果分析

利用最小二乘法、模擬退火算法以及遺傳算法在100 Hz～10 MHz頻率范圍內(nèi)對(duì)二階并聯(lián)形式的Cole-Cole理論模型在不同的信噪比時(shí)進(jìn)行反演,反演出公式中的7個(gè)參數(shù)和復(fù)電阻率曲線。表1給出了理論模型參數(shù)和反演初始參數(shù);圖5所示是在信噪比為30時(shí)反演所得曲線。使用這3種算法對(duì)理論模型(信噪比為30)進(jìn)行了20次反演,反演所得參數(shù)及曲線平均相對(duì)誤差如表2所示。

表1 理論模型參數(shù)及反演初始參數(shù)

表2 理論模型反演結(jié)果平均相對(duì)誤差表

圖6 理論模型反演迭代收斂圖

反演時(shí)對(duì)于最小二乘法設(shè)置了2組初值(見(jiàn)表1)。第1組設(shè)置恰當(dāng);第2組設(shè)置不恰當(dāng)以作為第1組的對(duì)比。反演結(jié)果表明,使用第1組初始值時(shí)最小二乘法可以得到準(zhǔn)確的解;使用第2組初始值時(shí)最小二乘法所得的結(jié)果誤差很大,這恰恰反映了最小二乘法對(duì)初始值具有依賴(lài)性,容易陷入局部極值而得不到準(zhǔn)確結(jié)果的局限性。模擬退火算法進(jìn)行反演,使用第2組初始值初始值依然可以得到準(zhǔn)確的反演結(jié)果,這也驗(yàn)證了作為全局尋優(yōu)的模擬退火算法,擺脫了對(duì)初始值的依賴(lài)性,避免陷入局部極值;針對(duì)群體開(kāi)始搜索的遺傳算法,作為一種并行的全局尋優(yōu)算法,同時(shí)對(duì)模型的7個(gè)參數(shù)進(jìn)行并行搜索,也可以得到準(zhǔn)確的反演結(jié)果。對(duì)于這3種算法,最小二乘法耗時(shí)最短,遺傳算法耗時(shí)最長(zhǎng)。算法的運(yùn)行時(shí)間受算法運(yùn)行平臺(tái)配置、算法內(nèi)部?jī)?yōu)化等影響,本文比較的是模擬退火算法和遺傳算法的收斂速度。圖6所示是同時(shí)循環(huán)2 000次時(shí)模擬退火算法與遺傳算法的迭代收斂圖。由圖6可知模擬退火算法的收斂速度慢于遺傳算法,這是因?yàn)槟M退火算法是通過(guò)搜索鄰域逐漸搜索全局最優(yōu)解,對(duì)整個(gè)搜索空間的狀況了解不多,不便于使搜索過(guò)程進(jìn)入最有希望的搜索區(qū)域,從而使得模擬退火算法的運(yùn)算效率不高;遺傳算法針對(duì)群體開(kāi)始搜索,運(yùn)算效率較高,但遺傳算法對(duì)局部搜索的并不充分。在實(shí)現(xiàn)過(guò)程中,模擬退火算法更加簡(jiǎn)單,容易實(shí)現(xiàn);遺傳算法過(guò)程繁瑣,計(jì)算過(guò)程中有較多的參數(shù),實(shí)現(xiàn)過(guò)程比較復(fù)雜。

通過(guò)對(duì)理論模型反演,驗(yàn)證了模擬退火算法和遺傳算法在反演工作中的可行性。利用第3部分所述的實(shí)驗(yàn)測(cè)量系統(tǒng)測(cè)量了不同飽和度人造巖心的復(fù)電阻率數(shù)據(jù),利用所測(cè)數(shù)據(jù)中第1組的實(shí)驗(yàn)實(shí)測(cè)數(shù)據(jù)的反演結(jié)果曲線圖如圖7所示,所選巖心的孔隙度為0.27,礦化度為1 kmg/L,含水飽和度為100%。反演時(shí)與理論模型反演相同,也設(shè)置了2組初始值,對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行了20次反演,反演得出的7個(gè)參數(shù)和曲線平均相對(duì)誤差如表3所示。由實(shí)測(cè)數(shù)據(jù)的反演結(jié)果可得出與前文一致的結(jié)論,二階乘積形式的Cole-Cole模型可以很好的與實(shí)測(cè)數(shù)據(jù)相匹配。通過(guò)對(duì)這3種算法的綜合對(duì)比,可以看出全局尋優(yōu)的模擬退火算法因其不依賴(lài)初始值的設(shè)置,可以找到全局最優(yōu)解,過(guò)程簡(jiǎn)單明了,容易實(shí)現(xiàn)的特性在實(shí)際工作中更有優(yōu)勢(shì)。

圖7 實(shí)測(cè)數(shù)據(jù)反演曲線圖

算法ρ0m1c1m2c2曲線平均相對(duì)誤差/%最小二乘法(曲線1)469.40.52480.1000.58790.99976.324×e-060.94114.7最小二乘法(曲線2)970.230.77240.82920.5640.99756.295×e-060.943074.2模擬退火算法484.450.54220.1010.60460.97806.337×e-060.94013.7遺傳算法500.030.55580.1030.62030.99796.349×e-060.93944.3

6 結(jié) 論

(1) 基于乘積形式的二階Cole-Cole可以很好地匹配實(shí)測(cè)數(shù)據(jù),可以應(yīng)用于實(shí)際反演。

(2) 在初值設(shè)置合理時(shí),最小二乘法可以快速準(zhǔn)確地得到反演結(jié)果;但當(dāng)初值設(shè)置不合理時(shí)最小二乘法容易陷入局部極值,無(wú)法得到正確的反演結(jié)果,因此,使用最小二乘法時(shí)需要對(duì)初始值進(jìn)行估計(jì),容易帶來(lái)誤差。

(3) 模擬退火算法和遺傳算法這2種全局尋優(yōu)算法,不依賴(lài)初始值的設(shè)置,避免了陷入局部極值,可以準(zhǔn)確得到全局最優(yōu)解,反演效果較好。

(4) 模擬退火算法由單點(diǎn)開(kāi)始搜索領(lǐng)域逐漸搜索全局最優(yōu)解,運(yùn)算效率不高,但過(guò)程簡(jiǎn)單明了,容易實(shí)現(xiàn);遺傳算法針對(duì)群體開(kāi)始并行搜索,運(yùn)算效率較高,但實(shí)現(xiàn)過(guò)程參數(shù)較多,比較復(fù)雜。

(5) 實(shí)際工作中使用模擬退火算法較好,擺脫了對(duì)初始值的依賴(lài),容易實(shí)現(xiàn),可以提高工作效率。

參考文獻(xiàn):

[1] Lesmes D P, Frye K M. Influence of Pore Fluid Chemistry on the Complex Conductivity and Induced Polarization Responses of Berea Sandstone [J]. Journal of Geophysical Research, 2001, 106: 4079-4090.

[2] Jaggar S R, Fell P A. Forward and Inverse Cole-Cole Modeling in the Analysis of Frequency Domain Electrical Impedance Data [J]. Exploration Geophysics, 1988, 19: 463-470.

[3] Xiang J, et al. Direct Inversion of the Apparent Complex-resistivity Spectrum [J]. Geophysics, 2000, 66(5): 1399-1404.

[4] Cao Z L, et al. Inversion Study of Spectral Induced Polarization Based on Improved Genetic Algorithm [C]∥Progress in Electromagnetic Research Symposium, Hangzhou, 2005.

[5] 曹中林, 昌彥君, 等. 基于演化算法的復(fù)電阻率頻譜參數(shù)反演 [J]. 工程地球物理學(xué)報(bào), 2005, 2(1): 33-38.

[6] 童茂松, 丁柱. 巖石復(fù)電阻率頻譜模型參數(shù)的反演 [J]. 測(cè)井技術(shù), 2006, 30(4): 303-305.

[7] 柯式鎮(zhèn), 馮啟寧, 何億成, 等. 電極法復(fù)電阻率測(cè)井研究 [J]. 石油學(xué)報(bào), 2006, 27(2): 89-82.

[8] Cole K S, Cole R H. Dispersion and Absorption in Dielectrics [J]. Chem Phys, 1941, 1: 341-351.

[9] Pelton W H, Jones N B, et al. Mineral Discrimination and Removal of Inductive Coupling with Muti-frequency IP [J]. Geophysics, 1978, 43(3): 588-609.

[10] Pelton W H, Smith B D, Sill W R. Interpretation of Complex Resistivity and Dielectric Data, part II [J]. Geophysics, 1984, 29(4): 11-45.

[11] 張雷潔, 柯式鎮(zhèn), 等. 40 Hz～110 MHz巖石復(fù)電阻率頻譜實(shí)驗(yàn)測(cè)量 [C]∥中國(guó)地球物理2013——第二十專(zhuān)題論文集, 中國(guó)地球物理學(xué)會(huì)第二十九屆年會(huì), 昆明, 2013.

[12] 柯式鎮(zhèn), 馮啟寧, 孫艷茹. 復(fù)電阻率頻散模型及其參數(shù)的獲取方法 [J]. 測(cè)井技術(shù), 1999, 23(6): 416-418.

[13] Plackett R L. The Discovery of the Method of Least Squares [J]. Biometrika, 1972, 59: 239-251.

[14] Kirkpatrick S, Gelatt C D, Vecchi M P. Optimization by Simulated Annealing [J]. Science, 1983, 220(4598): 671-680.

[15] Stoffa P L, Sen M K. Non-linear Multiparameter Optimization Using Genetic Algorithm: Inversion of Plane Wave Seismograms [J]. Geophysics, 1991, 56: 1794-1810.

[16] Hu B, Tang G, Ma J W, et al. Parametric Inversion of Viscoelastic Media from VSP Data Using a Genetic Algorithm [J]. Applied Geophysics, 2007, 4(3): 194-200.