陳 恒， 嚴(yán)良俊，代小威，李鵬飛

(1.長江大學(xué) 湖北非常規(guī)油氣協(xié)同創(chuàng)新中心，湖北 武漢 430100；2.中國石油塔里木油田分公司 勘探開發(fā)研究院，新疆 庫爾勒 841000)

?

模擬退火法與阻尼最小二乘法聯(lián)合反演巖石激電譜參數(shù)

陳 恒1， 嚴(yán)良俊1，代小威1，李鵬飛2

(1.長江大學(xué) 湖北非常規(guī)油氣協(xié)同創(chuàng)新中心，湖北 武漢 430100；2.中國石油塔里木油田分公司 勘探開發(fā)研究院，新疆 庫爾勒 841000)

Cole—Cole模型是一種定量描述激電響應(yīng)譜特性的方法，即可以用4個(gè)參數(shù)來定量描述巖石激電響應(yīng)特性,這4個(gè)參數(shù)是巖礦石激電響應(yīng)特征的基本參數(shù)。由于不同巖礦石具有不同的激電譜特征，因而存在著利用這4個(gè)Cole—Cole參數(shù)來識(shí)別它們的可能性。阻尼最小二乘法是應(yīng)用地球物理學(xué)中常用的一種最優(yōu)化反演技術(shù)，廣泛地用在譜激電反演之中，具有很快的收斂速度。然而阻尼最小二乘法對初始值很依賴，非線性模擬退火反演方法能夠很好地消除這種影響。因此聯(lián)合這兩種方法反演巖石激電譜參數(shù)，經(jīng)理論和實(shí)測數(shù)據(jù)計(jì)算表明，反演能較好地?cái)M合實(shí)測曲線，取得了一定的效果。

模擬退火法；阻尼最小二乘法；激電譜；反演； Cole—Cole模型

1 引 言

巖礦石的復(fù)電阻率頻譜，一般都滿足Cole—Cole模型，可以用四個(gè)Cole—Cole參數(shù)描述。一般，巖礦石四個(gè)Cole—Cole參數(shù)是通過對巖礦石的復(fù)電阻率頻譜測量數(shù)據(jù)進(jìn)行最優(yōu)化反演來確定的。當(dāng)其頻譜能很好擬合實(shí)測的巖礦石的復(fù)電阻率頻譜時(shí)，擬合所用的Cole—Cole參數(shù)就是巖礦石的Cole—Cole參數(shù)。頻譜激電資料的譜激電反演，不論是常規(guī)的還是求真參數(shù)的，一般都采用阻尼最小二乘法最優(yōu)化技術(shù)[1]。

阻尼最小二乘法又稱馬奎特迭代法，是應(yīng)用地球物理學(xué)中常用的一種最優(yōu)化反演技術(shù)。它通過自適應(yīng)調(diào)整阻尼因子來達(dá)到收斂特性，具有更高的迭代收斂速度，但在阻尼最小二乘法的計(jì)算過程中，初值是一個(gè)很重要的因素，若選擇的初值接近真解時(shí)，收斂速度很快且能得到收斂，但如果初值遠(yuǎn)離真解時(shí)，優(yōu)化結(jié)果往往過早地陷入局部最優(yōu)解，從而得到的結(jié)果完全背離真解。要解決該問題，方法一是要通過大量的原始信息對其初值有一個(gè)大概的估計(jì)，但這在現(xiàn)實(shí)工作中是很難做到的，另外一種方法就是選擇一種合適的方法與其相結(jié)合，消除阻尼最小二乘法對初始值的依賴，并且還應(yīng)具有很快的收斂速度[2,4]。

非線性模擬退火反演方法很好地克服了初始值設(shè)定這一缺點(diǎn)。每修改一次模型，將會(huì)隨機(jī)產(chǎn)生一個(gè)新的狀態(tài)模型，如果新模型能夠被系統(tǒng)接受，則進(jìn)行‘降溫’處理，產(chǎn)生下一個(gè)新模型。這種接受新模型的方式使其成為一種全局最優(yōu)算法[2]。該方法能跳出局部最優(yōu)的“陷阱”，避免了在反演迭代過程中陷入目標(biāo)函數(shù)的局部最小值，從而得到全局最優(yōu)解，并得到理論證明和實(shí)際應(yīng)用驗(yàn)證[3]。模擬退火方法雖然可以達(dá)到全局最小，但計(jì)算量太大，尤其是解的個(gè)數(shù)太多和解空間范圍太大時(shí)，該方法很難收斂[5]。

利用阻尼最小二乘法和模擬退火方法各自的優(yōu)缺點(diǎn)，國內(nèi)外利用這兩種方法進(jìn)行反演的例子很多。劉建軍等利用阻尼最小二乘法與模擬退火法結(jié)合實(shí)現(xiàn)非線性模型參數(shù)估計(jì)，并將此算法運(yùn)用于植物光合作用光響應(yīng)新模型對水稻的光響應(yīng)曲線的擬合，得到的擬合參數(shù)與DPS擬合值極為接近[3]。王偉男等討論了對于觀測數(shù)據(jù)的非線性擬合，目前廣為采用的是基于蒙特卡羅法的模擬退火方法和阻尼最小二乘法[6]。王芳等基于模擬退火算法中的全局收斂思想，建立了一種解非線性反演問題的新策略[7]。

2 基本原理與方法

2.1 復(fù)電阻率的Cole—Cole表達(dá)式

大多數(shù)巖，礦石的激電譜可以用4個(gè)Cole—Cole參數(shù)來很好的定量描述，一般都滿足Cole—Cole模型，可以根據(jù)巖礦石的這4個(gè)參數(shù)來評價(jià)激電異常。其復(fù)電阻率的Cole—Cole表達(dá)式為：

(1)

式中，ρ0是零頻率時(shí)的電阻率，單位是Ω·m；m是極化率，無量綱；c為復(fù)電阻率隨頻率變化程度的參數(shù)，稱為頻率相關(guān)系數(shù)，無量綱；τ為時(shí)間常數(shù)[1]，單位為s。

2.2 阻尼最小二乘法基本原理

阻尼最小二乘法是在最小二乘法和最速下降法的基礎(chǔ)上發(fā)展起來的，它結(jié)合了最小二乘法與最速下降法的優(yōu)點(diǎn)，并在兩種方法之間取某種折衷，力圖以最大的步長，同時(shí)又靠近最速下降方向，以保證穩(wěn)定收斂，并加快收斂速度。這種方法又稱馬奎特法[8,9]。

假設(shè)要反演的一組實(shí)測地球物理數(shù)據(jù)向量用D表示，設(shè)其有n個(gè)數(shù)據(jù)；設(shè)正演模型的已知模型參數(shù)向量為X(頻點(diǎn)的頻率)，它也有n個(gè)數(shù)；設(shè)待求參數(shù)向量為P(Cole—Cole參數(shù))，有m個(gè)數(shù)[1]。

正演模型是P和X的函數(shù)，設(shè)為G(P,X)，顯然應(yīng)當(dāng)有D=G(P,X)。

對于正演模型G(P,X)，按照臺(tái)勞級(jí)數(shù)展開并略去所有二階以上的項(xiàng)可以得到：

G(P,X) =G(P0,X)+

(2)

上式可以寫成矩陣形式：ΔG=AΔP

(3)

其中，ΔG是數(shù)據(jù)差列向量；A是雅克比矩陣；ΔP是參數(shù)差列向量；最小二乘法的擬合優(yōu)度用x2值衡量，根據(jù)上式即可得到目標(biāo)函數(shù)：

(4)

上式是Δpj的多元函數(shù)，取得極小值的充分條件是ε(P0)對所有的Δpj的偏導(dǎo)數(shù)為零。即：

=0

(5)

進(jìn)一步化簡可以得到:

(ATA)ΔP=ATΔG

(6)

由式(6)可以得到：

ΔP=(ATA)-1ATΔG

(7)

加入阻尼因子k,其中I是單位矩陣，即得

ΔP=(ATA+kI)-1ATΔG

(8)

ΔP即是使目標(biāo)函數(shù)ε往極小方向減小的最佳修正步長，對于絕大多數(shù)地球物理問題來說，G(P,X)不是P的線性函數(shù)，因此要經(jīng)過若干次類似迭代，才能使目標(biāo)函數(shù)滿足某一精度要求，這里采用一邊迭代一邊調(diào)整阻尼因子[10]。

2.3 模擬退火法基本原理

近年來，模擬退火法已在地球物理資料的反演中被廣泛地采用，并取得了一些令人矚目的成果，它已成為目前非線性反演的重要方法之一。它的基本思想是，把物理系統(tǒng)的能量作為最優(yōu)化問題的目標(biāo)函數(shù)，其能量達(dá)到最小，系統(tǒng)處于最優(yōu)狀態(tài)。整個(gè)最優(yōu)化過程中，溫度是一個(gè)至關(guān)重要的參數(shù)[8]。據(jù)統(tǒng)計(jì)熱力學(xué)，物體中每個(gè)分子的狀態(tài)服從吉布斯概率分布，即為：

(9)

式中：E(ri)為第i個(gè)分子的能量函數(shù)；ri為第i個(gè)分子所處的狀態(tài)；k為玻爾茲曼常數(shù)；T為溫度；而ρ(ri)為第i個(gè)分子的概率密度。其中溫度T只是一個(gè)控制參數(shù)，玻爾茲曼常數(shù)k起一個(gè)尺度因子的作用，實(shí)際應(yīng)用取為1。顯然，T較小時(shí)，對ΔE起著一種放大作用，使低溫時(shí)不易接受模型修改。這相當(dāng)于物體在低溫時(shí)分子被束縛在平衡位置附近[8,9]。

2.4 模擬退火法與阻尼最小二乘法聯(lián)合反演流程

圖1 模擬退火法與阻尼最小二乘法聯(lián)合反演流程Fig.1 Flow chart of the joint inversion of simulated annealing method and damped least square method

3 反演結(jié)果分析

3.1 理論模型計(jì)算

選擇Cole—Cole模型參數(shù)ρ0=25 Ω·m，m=0.6，τ=1 s，c=0.5，頻率范圍為0.01Hz～1 kHz，用所給的參數(shù)帶入Cole—Cole模型進(jìn)行理論計(jì)算[11，12]，然后把計(jì)算出來的數(shù)據(jù)作為初始理論數(shù)據(jù)，用模擬退火法和阻尼最小二乘法聯(lián)合反演進(jìn)行復(fù)電阻率的振幅反演。本文先隨機(jī)給一個(gè)初始值，用模擬退火法進(jìn)行反演，然后把模擬退火法反演出來的參數(shù)當(dāng)作阻尼最小二乘法初始值進(jìn)行反演，反演出來的參數(shù)見表1，聯(lián)合反演結(jié)果見圖2。為了使本方法更具有說服力，筆者在計(jì)算出來的理論數(shù)據(jù)中加入±5%以內(nèi)的隨機(jī)相對噪聲，把加入噪聲的數(shù)據(jù)再來進(jìn)行聯(lián)合反演，反演出來的參數(shù)見表1，聯(lián)合反演結(jié)果見圖3。

表1 理論模型聯(lián)合反演參數(shù)

圖2 理論模型聯(lián)合反演結(jié)果Fig.2 Joint inversion result of theoretical model

圖3 加入噪聲理論模型聯(lián)合反演結(jié)果Fig.3 Joint inversion result of theoretical model after adding noise

3.2 阻尼最小二乘法單獨(dú)反演結(jié)果分析

本次反演通過實(shí)驗(yàn)實(shí)測某地區(qū)巖石1號(hào)點(diǎn)、2號(hào)點(diǎn)和3號(hào)點(diǎn)的復(fù)電阻率數(shù)據(jù)，研究復(fù)電阻率隨頻率的變化關(guān)系。實(shí)測頻率范圍是0.01～10 kHz，考慮到電磁耦合效應(yīng)，我們只研究0.01～1 kHz的實(shí)驗(yàn)數(shù)據(jù)[13,14]。本次反演用單Cole—Cole模型進(jìn)行反演,并只對幅值進(jìn)行反演。先編程用阻尼最小二乘法進(jìn)行反演Cole—Cole模型參數(shù)。在對1號(hào)點(diǎn)反演過程中，取初值ρ0=200 Ω·m，m=0.2，τ=100 s,c=0.1無法收斂，然后通過不斷變化取初值，最后當(dāng)取初值ρ0=100 Ω·m，m=0.6,τ=0.5 s,c=0.25時(shí)反演得到了巖石頻譜參數(shù)，但也不是巖石真頻譜參數(shù)。同樣2號(hào)點(diǎn)和3號(hào)點(diǎn)也是通過不斷變化取初值才最終反演出Cole—Cole模型參數(shù)。1號(hào)點(diǎn)、2號(hào)點(diǎn)和3號(hào)點(diǎn)的反演參數(shù)見表2，反演結(jié)果見圖4～圖6。這說明阻尼最小二乘法反演對初值嚴(yán)重依賴，如果初值選取不當(dāng)，反演可能無法收斂或者不能收斂到全局最小[12]。

表2 各個(gè)點(diǎn)實(shí)測曲線反演參數(shù)

圖4 1號(hào)點(diǎn)的實(shí)測曲線聯(lián)合反演結(jié)果 Fig.4 Joint inversion result of the measured curve of the 1st point

圖5 2號(hào)點(diǎn)的實(shí)測曲線聯(lián)合反演結(jié)果Fig.5 Joint inversion result of the measured curve of the 2nd point

圖6 3號(hào)點(diǎn)的實(shí)測曲線聯(lián)合反演結(jié)果Fig.6 Joint inversion result of the measured curve of the 3rd point

3.3 阻尼最小二乘法和模擬退火法聯(lián)合反演分析

由于阻尼最小二乘法對初始值嚴(yán)重依賴，因此本文通過阻尼最小二乘法和模擬退火法聯(lián)合反演的思想來消除這種影響。先編程用非線性模擬退火法反演得到Cole—Cole模型參數(shù)，并作為阻尼最小二乘法反演的初始參數(shù)，這樣阻尼最小二乘法就會(huì)快速反演得到巖石頻譜參數(shù)。1號(hào)點(diǎn)、2號(hào)點(diǎn)和3號(hào)點(diǎn)的聯(lián)合反演參數(shù)見表2。

圖4～圖6分別給出了1號(hào)點(diǎn)、2號(hào)點(diǎn)和3號(hào)點(diǎn)實(shí)測曲線，阻尼最小二乘法單獨(dú)反演結(jié)果，模擬退火法單獨(dú)反演結(jié)果，以及阻尼最小二乘法和模擬退火法聯(lián)合反演結(jié)果。將實(shí)測曲線與反演得到的電圖阻率進(jìn)行對比，結(jié)果如圖4～圖6所示。

4 結(jié) 論

本文充分利用模擬退火法不依賴初值以及能收斂到全局最小的優(yōu)點(diǎn)，編程實(shí)現(xiàn)了阻尼最小二乘法和模擬退火法聯(lián)合反演的程序，并反演得到了巖石頻譜參數(shù)，取得了一定的效果。經(jīng)過對理論和實(shí)測的數(shù)據(jù)計(jì)算表明，這種聯(lián)合反演能較好地?cái)M合實(shí)測曲線，降低了阻尼最小二乘法對初始值的依賴程度，并具有一定的抗噪能力。

[1]劉崧.譜激電法[M].武漢：中國地質(zhì)大學(xué)出版社，1997.

[2]羅延鐘,張桂青.頻率域激電法原理[M].北京:地質(zhì)出版社,1988.

[3]劉建軍，陳明峰.阻尼最小二乘法與模擬退火法結(jié)合實(shí)現(xiàn)非線性模型參數(shù)估計(jì)[J].井岡山大學(xué)學(xué)報(bào),2010,31(6):10-14.

[4]武菊，任鵬.基于LM和SA的混合優(yōu)化算法[J].內(nèi)江師范學(xué)院學(xué)報(bào)，2010，25 (8):32-34.

[5]嚴(yán)良俊，胡文寶.大地電磁測深資料的二次函數(shù)逼近非線性反演[J].地球物理學(xué)報(bào)，2004,47(5):935-940.

[6]王偉男，童茂松，陳國華.泥質(zhì)砂巖的物理性質(zhì)及其測井應(yīng)用[M].北京：石油工業(yè)出版社,2004 .

[7]王芳，陳生昌.解非線性反演問題的新策略[J].石油物探,1999,38(3):76-85.

[8]王家映.地球物理反演理論[M].北京：高等教育出版社，2002.

[9]姚姚.地球物理反演基本理論與應(yīng)用方法[M].武漢：中國地質(zhì)大學(xué)出版社,2002.

[10]薛嘉慶.最優(yōu)化原理與方法[M].北京：冶金工業(yè)出版社，1992.

[11]李鵬飛，嚴(yán)良俊，謝興兵,等.基于最小二乘的巖石頻譜參數(shù)反演初始值選取的討論[J].地震工程學(xué)報(bào),2014，36(3):569-574.

[12]向葵，胡文寶，嚴(yán)良俊,等.川黔地區(qū)頁巖復(fù)電阻率的頻散特性[J].石油地球物理勘探，2014，49(5):1 013-1 019.

[13]周新鵬，鄒長春，喬祎,等.巖礦石電性參數(shù)測定實(shí)驗(yàn)研究[J].工程地球物理學(xué)報(bào)，2015,12(2):205-213.

[14]李鵬飛，嚴(yán)良俊，余剛.南方頁巖巖芯的復(fù)電阻率測試與分析[J].工程地球物理學(xué)報(bào)，2014,11(3):383-386.

The Joint Inversion of Simulated Annealing Method and Damped Least Square Method for IP Spectrum Parameters of Rock

Chen Heng1,Yan Liangjun1,Dai Xiaowei1,Li Pengfei2

(1.HubeiCooperativeInnovationCenterofUnconventionalOilandGas,YangtzeuUniversity,WuhanHubei430100,China；2.ResearchInstituteofExplorationandDevelopment,TarimOilfieldCompany,PetroChina,KorlaXinjiang841000,China)

The Cole—Cole model is a quantitative method for describing the characteristics of the excitation response spectrum, that is, the spectral characteristics of the IP response can be quantitatively described with 4 parameters, and these 4 parameters are the basic parameters of the characteristics of rock IP response. Because different rocks have different IP spectral characteristics, it is possible to identify them by using the four Cole—Cole parameters. Damped least square method is one of the commonly used optimization inversion techniques in the application of geophysics, which is widely used in the spectral inversion and has a fast convergence rate. However, the damped least squares method is very dependent on the initial value, and the nonlinear simulated annealing inversion method can eliminate the effect of this method. So the two methods were combined to inverse the parameters of the rock. The calculation of the theory and the measured data showed that the joint inversion can better fit the measured curves, and good results had been obtained.

simulated annealing method; damped least squares method; IP spectrum; inversion; Cole—Cole model

1672—7940(2016)02—0170—05

10.3969/j.issn.1672-7940.2016.02.006

國家自然科學(xué)基金項(xiàng)目(編號(hào)：U1562109,41274082)；國家“973”項(xiàng)目(編號(hào)：2013CB228605)

陳 恒(1990-)，男，碩士研究生，主要研究方向?yàn)殡姶欧碧?。E-mail：1033256984@qq.com

P631.3

A

2015-11-10