基于改進(jìn)粒子群算法的疊前彈性阻抗反演

張　進(jìn), 陳曉琦, 邢　磊, 安玲芳, 趙　璽, 安振芳

(中國海洋大學(xué)　海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，青島　266100)

?

基于改進(jìn)粒子群算法的疊前彈性阻抗反演

張進(jìn), 陳曉琦, 邢磊, 安玲芳, 趙璽, 安振芳

(中國海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，青島266100)

摘要：彈性阻抗反演是結(jié)合聲阻抗反演與AVO反演的疊前地震反演技術(shù)，能夠克服疊后波阻抗反演的缺陷，反映振幅隨偏移距變化的信息，已經(jīng)廣泛應(yīng)用于地震巖性識(shí)別和流體特征的獲取。常規(guī)的線性迭代彈性阻抗反演方法存在依賴初始模型、容易陷入局部極值等缺陷。針對(duì)這一問題，提出了一種基于改進(jìn)粒子群算法的彈性阻抗非線性反演方法，并利用該算法對(duì)勝利油田某工區(qū)地震資料進(jìn)行了彈性阻抗反演，獲得了多個(gè)彈性參數(shù)剖面，與實(shí)際鉆井結(jié)果相符，該方法為復(fù)雜油氣藏的勘探開發(fā)提供了一種有效可行的途徑。

關(guān)鍵詞：彈性阻抗反演； 粒子群算法； 非線性反演

0引言

常規(guī)疊后波阻抗反演技術(shù)基于地震波垂直入射的假設(shè)，其得到的反射振幅為共中心點(diǎn)道集疊加平均的結(jié)果，故不能反映反射振幅隨炮檢距或入射角的變化。為獲取潛在的AVO效應(yīng)，Connolly[1]將縱波反射系數(shù)隨炮檢距的變化引入地震道正、反演問題中，為孔隙流體和巖性識(shí)別提出了彈性阻抗EI(Elastic Impedance)的概念，能夠提供振幅隨著偏移距變化的信息，解決了在較大炮間距入射角情況下的縱波正、反演問題，具有足夠的精度和良好的保真性，從EI提出至今，彈性阻抗反演取得了突破性的發(fā)展；Cambois[2]對(duì)AVO反演和彈性波阻抗(EI)反演方法進(jìn)行了比較，得出彈性阻抗的反演結(jié)果更理想，在抗噪能力方面比AVO反演更有優(yōu)勢(shì)；Mallick等[3]討論了Connolly彈性阻抗公式的應(yīng)用及限制條件；Whitcombe D N[4]引入縱波速度、橫波速度、密度三個(gè)參量，對(duì)Connolly提出的彈性阻抗進(jìn)行了歸一化；Whitcombe 等[5]在前人對(duì)于流體因子和巖性預(yù)測(cè)等方面研究成果的基礎(chǔ)上，提出擴(kuò)充彈性波阻抗(EEI)，可用于巖性和流體預(yù)測(cè)的;馬勁風(fēng)[6]基于Subhashis Mallick[7]和Yanghua Wang[8]提出了廣義彈性阻抗(GEI)的概念，并給出使用彈性模量表示的Zoeppritz 方程的簡(jiǎn)化公式；倪逸[9]通過對(duì)目前幾種彈性波阻抗計(jì)算方法的分析，提出了一種基于范數(shù)動(dòng)態(tài)可調(diào)的改進(jìn)方法；印興耀等[10]提出利用Connolly彈性阻抗方程從三個(gè)角度反演結(jié)果中提取縱、橫波速度和密度參數(shù)的方法，可以對(duì)地下儲(chǔ)層的展布及含油氣性進(jìn)行預(yù)測(cè)；王保麗等[11]提出了一種新的基于佐普里茲Gary近似方程的彈性阻抗反演公式，能夠減小間接計(jì)算產(chǎn)生的累計(jì)誤差。

常用的反演方法一般為局部線性迭代方法或者將非線性問題線性化，具有依賴初始值、易陷入局部極值的缺陷。對(duì)于精度要求較高的EI反演來說，采用非線性優(yōu)化方法，其解的性質(zhì)、狀態(tài)均有良好的表現(xiàn)。近年來，一些全局優(yōu)化特性并且通用性好的搜索算法(如遺傳算法、模擬退火算法、粒子群算法、人工神經(jīng)網(wǎng)絡(luò)算法等)，已經(jīng)廣泛應(yīng)用在地震反演中。其中，粒子群算法對(duì)初始模型的依賴性小，參數(shù)設(shè)置少，收斂速度快，能夠保證地震資料反演結(jié)果的精度。這里在基本粒子群算法中加入模擬退火算子用于反演彈性阻抗，通過對(duì)理論模型的分析和實(shí)際地震資料的反演，證明該方法的合理性和有效性。

1彈性阻抗的基本原理

與AVO的理論基礎(chǔ)相同，彈性阻抗是基于平面波在兩種半無限空間彈性介質(zhì)分界面上的反射和透射所滿足的Zoeppritz方程及其不同簡(jiǎn)化形式，是聲阻抗(Acoustic impedance, AI)的推廣。聲阻抗表示入射角為零時(shí)的特例，而彈性阻抗是關(guān)于入射角的函數(shù)，可以標(biāo)定和反演非零偏移數(shù)據(jù)體，其精度已經(jīng)在廣泛的應(yīng)用中得到證明。

Connolly根據(jù)Zoeppritz方程的三項(xiàng)Aki & Richards簡(jiǎn)化公式推導(dǎo)出含有縱、橫波速度、密度和入射角參數(shù)的彈性阻抗公式，如式(1)所示。

(2)

常規(guī)彈性阻抗反演流程如圖1所示。

圖1　彈性阻抗反演流程示意圖Fig．1　The procedures of elastic impedance inversion

反演步驟為：

1)由于EI是關(guān)于縱、橫波速度，密度和入射角的函數(shù)，在反演之前需要將地震數(shù)據(jù)的偏移數(shù)據(jù)體轉(zhuǎn)化為角度數(shù)據(jù)體。

2)根據(jù)測(cè)井曲線或者巖心數(shù)據(jù)，計(jì)算出井旁道EI偽測(cè)井曲線。

3)角度子波的提取。子波提取是地震反演的關(guān)鍵之一，可直接影響反演精度，在提取之前要先對(duì)多口井的測(cè)井曲線進(jìn)行層位標(biāo)定。

4)對(duì)控制點(diǎn)上的EI曲線進(jìn)行內(nèi)插，利用角道集部分疊加資料與井旁道相應(yīng)入射角的彈性阻抗建立波阻抗剖面的低頻模型，進(jìn)行低頻成分的補(bǔ)充。

5)用非線性反演方法計(jì)算帶限的彈性阻抗剖面，將低頻模型的彈性波阻抗與帶限的彈性波阻抗剖面融合，最終可以獲得與實(shí)際情況相符合的彈性波阻抗剖面。

在得到不同角度下的彈性阻抗之后，可根據(jù)變換后的Connolly方程反演得到縱波速度Vp、橫波速度Vs、密度ρ等巖性參數(shù)：

(3)

與聲阻抗AI相比，彈性阻抗不存在聲阻抗垂直入射的前提，能夠克服由于疊加導(dǎo)致的有效信息損失的缺陷，并且EI由于含有AVO信息，其對(duì)油氣飽和度更為敏感，故可以更全面、直觀地反應(yīng)與油氣相關(guān)的信息，為油氣藏的識(shí)別提供依據(jù)，降低油氣檢測(cè)的多解性。

地震資料的彈性阻抗反演實(shí)際上是一個(gè)最優(yōu)化求解的非線性問題[12]，常規(guī)的線性算法雖收斂速度快，但易陷入局部最優(yōu)。粒子群算法是一種基于迭代的優(yōu)化工具，具有較強(qiáng)的搜索能力，容易理解，易于實(shí)現(xiàn)，近年來在模式分類、函數(shù)尋優(yōu)、系統(tǒng)控制、以及工程應(yīng)用等多個(gè)領(lǐng)域得到了成功地應(yīng)用[13]。作者嘗試在基本粒子群算法中加入模擬退火算子，并用其進(jìn)行彈性阻抗的反演，取得了良好的效果。

2粒子群算法

2.1基本粒子群算法

粒子群算法(Particle Swarm Optimization，PSO)來源于鳥群和魚群的群體運(yùn)動(dòng)行為，最早由Kennedy博士和Eberhart博士[14-15]提出，其數(shù)學(xué)描述如下：

vid(t+1)=wvid(t)+c1r1(pid(t)-xid(t))+

c2r2(gd(t)-xid(t))

(4)

1≤i≤m1≤d≤n

(5)

其中：w為慣性權(quán)重因子；c1和c2為加速因子，均為正常數(shù)，c1為調(diào)節(jié)使粒子飛向自身最好位置方向的步長(zhǎng)，c2為調(diào)節(jié)粒子飛向全局最好位置的步長(zhǎng)；t是當(dāng)前迭代次數(shù)；r1、r2是[0,1]內(nèi)的隨機(jī)數(shù)。

(6)

可設(shè)定上限速度Vmax來確保粒子不會(huì)由于速度過大而錯(cuò)過最優(yōu)解，根據(jù)式(6)對(duì)粒子的速度進(jìn)行限制，隨機(jī)產(chǎn)生粒子的速度和初始位置，然后按照式(4)、式(5)進(jìn)行迭代，直到找到符合條件的解。粒子群算法實(shí)現(xiàn)方便、收斂速度快，參數(shù)設(shè)置少，是一種高效的搜索方法，但“早熟”收斂和易陷入局部極值是其最主要的缺點(diǎn)。

2.2加入模擬退火算子的隨機(jī)粒子群算法

在上述的基本粒子群算法中，當(dāng)w=0，粒子的進(jìn)化方程就變?yōu)槭?7)。

(7)

(8)

(9)

(10)

(11)

為使粒子j以較大概率趨于最優(yōu)解，作者提出加入模擬退火算子的隨機(jī)粒子群算法(ASPSO)：以當(dāng)前歷史最優(yōu)位置gd為初始狀態(tài)，選取初始溫度T=T0，采用經(jīng)典退火過程Tk=akT0，根據(jù)公式(12)產(chǎn)生下一個(gè)狀態(tài)：

(12)

(13)

模擬退火方法本身具有很好的全局收斂性，因而采用該方法生成微粒j并根據(jù)式(12)、式(13)進(jìn)行狀態(tài)更新，不會(huì)對(duì)隨機(jī)粒子群算法的全局收斂性產(chǎn)生負(fù)面影響。

2.3改進(jìn)粒子群算法性能測(cè)試

為了檢測(cè)粒子全優(yōu)化算法的性能，用Generalized Rastrigin函數(shù)對(duì)其進(jìn)行測(cè)試，Rastrigin函數(shù)是一個(gè)具有大量局部極優(yōu)的典型復(fù)雜多峰函數(shù)(圖2)，在運(yùn)算過程中容易使算法陷入局部極優(yōu)，卻不能得到全局最優(yōu)解，其表達(dá)式為式(14)。

-5.12≤xi≤5.12

(14)

將基本PSO算法與ASPSO算法分別運(yùn)行50次，迭代終止條件為迭代誤差達(dá)到1e-7，或迭代次數(shù)達(dá)到10 000次，測(cè)試并對(duì)比其收斂速度和精度，基本PSO算法與ASPSO算法反演Generalized Rastrigin函數(shù)(包括2維和5維)的參數(shù)設(shè)定如下：

基本PSO算法：m=30，w=0.9，c1=c2=1.8，vmax=5.12，xmax=5.12，xmin=-5.12。

ASPSO算法：m=30，w=0.5，c1=c2=1.8，vmax=5.12，xmax=5.12，xmin=-5.12，采用經(jīng)典退火方式，初始溫度T0=100 000，降溫系數(shù)a=0.95。

圖2　Generalized Rastrigin函數(shù)Fig．2　Generalized Rastrigin function

兩種算法基于Generalized Rastrigin函數(shù)的計(jì)算結(jié)果如表1所示，從表1中可以看出，ASPSO算法在收斂速度和精度上較普通PSO算法都有很大改善。對(duì)于復(fù)雜多峰Rastrigin函數(shù)，ASPSO算法運(yùn)算效率是普通PSO算法6倍～20倍，并且解的質(zhì)量有很大提高。

表1　基本PSO算法和ASPSO算法基于Generalized Rastrigin函數(shù)的結(jié)果對(duì)比

3模型驗(yàn)證

在地震勘探中，無噪時(shí)的地震道記錄S(t)可以表示為反射系數(shù)r(t)與地震子波w(t)的褶積：

(15)

由方程(15)遞推反演，可以得到彈性阻抗的遞推公式：

(16)

彈性阻抗反演可以看作函數(shù)優(yōu)化求極值問題，即要求取每個(gè)角道集的反射系數(shù)序列r(t)，使其與角度子波w(t)褶積后與實(shí)際地震記錄之間的誤差函數(shù)F(r)最小，則優(yōu)化目標(biāo)函數(shù)為式(17)。

(17)

采用改進(jìn)粒子群算法進(jìn)行彈性阻抗反演的步驟如下：

1)抽角道集并提取各個(gè)角道集的子波。這里模型中抽取5°、25°、45°三個(gè)角道集作為觀測(cè)值。

2)根據(jù)混合粒子群算法生成各角道集的反射系數(shù)序列r(t)，混合粒子群的每一維就對(duì)應(yīng)地下每一地層的反射系數(shù)，任意粒子在反射系數(shù)空間進(jìn)行搜索，每搜索到的一個(gè)新位置就相當(dāng)于得到了一個(gè)反射系數(shù)序列的候選解，根據(jù)式(15)生成模型的地震道記錄S(t)。

3)根據(jù)式(17)計(jì)算誤差函數(shù)F(r)。

5)利用優(yōu)化的粒子群算法進(jìn)行迭代，優(yōu)化最優(yōu)解，直到滿足截止條件。

作者將含有100個(gè)時(shí)間樣點(diǎn)的測(cè)井曲線作為數(shù)據(jù)模型，采樣間隔1 ms?；旌狭Ｗ尤核惴ǖ膮?shù)設(shè)置為：種群大小m=30；慣性因子w=0.1；學(xué)習(xí)因子c1=c2=1.8,xmin=-1，xmax=1，vmax=1，誤差要求小于0.000 001；采用經(jīng)典退火方式，初始溫度T0=100 000，降溫系數(shù)a=0.95。

圖3　角道集反演的反射系數(shù)和合成地震記錄與真實(shí)的反射系數(shù)和實(shí)際地震記錄對(duì)比Fig．3　The inversion coefficient of 5° gathers and synthetic seismic compare with the actual coefficient and records(a)反演的反射系數(shù)與實(shí)際反射系數(shù)對(duì)比;(b)反演的地震記錄與實(shí)際地震記錄對(duì)比

僅以5°角道集反演結(jié)果(圖3)說明基于混合粒子群算法的彈性阻抗反演方法的有效性。從圖3看出，反演得到的反射系數(shù)和地震記錄與實(shí)際值都吻合得很好。根據(jù)此方法，分別反演得到5°、25°和45°角道集的反射系數(shù)序列，然后遞推得出這三個(gè)角道集的彈性阻抗曲線，如圖4所示。

得到三個(gè)角度的彈性阻抗曲線后，計(jì)算各層的縱波阻抗、橫波阻抗、密度、縱波速度、橫波速度等彈性參數(shù)(虛線)，以上參數(shù)曲線與原測(cè)井曲線對(duì)比如圖5所示。由圖5可以看出，反演得到的各彈性參數(shù)曲線與測(cè)井?dāng)?shù)據(jù)曲線吻合得很好。

圖4　反演得到的5°、25°、45°三個(gè)角度的EI曲線Fig．4　The EI inversion curves of 5°、25°、45°

虛線表示根據(jù)反演的EI計(jì)算出的各彈性參數(shù)，實(shí)線表示原測(cè)井曲線圖5　根據(jù)反演的EI計(jì)算出的各彈性參數(shù)與原測(cè)井曲線對(duì)比Fig．5　Comparison of the elastic parameters curves and the original logging curves

4應(yīng)用實(shí)例及效果分析

選取勝利油田某工區(qū)二維數(shù)據(jù)進(jìn)行彈性阻抗反演，對(duì)這里反演算法加以檢測(cè)。該區(qū)原始道集數(shù)據(jù)信噪比低，層間反射雜亂，不利于AVO屬性分析。為了使最終的疊前振幅數(shù)據(jù)能準(zhǔn)確反映地下界面的反射強(qiáng)度，在反演之前需要對(duì)地震資料做以下振幅保持處理：①球面補(bǔ)償；②地表一致性反褶積；③隨機(jī)噪聲衰減；④反Q濾波；⑤疊前時(shí)間偏移。

過da21井道集經(jīng)預(yù)處理后(圖6)，地震資料的品質(zhì)得到較大改善，信噪比提高，隨機(jī)噪聲得到極大壓制，目標(biāo)層段(圖6中黃色矩形所示)AVO現(xiàn)象明顯，有利于進(jìn)行疊地震屬性的反演。

圖6　過da21井預(yù)處理后道集Fig．6　Traces record after preprocessing over da21 well

通過對(duì)該區(qū)井資料的分析，發(fā)現(xiàn)該區(qū)沙二段為砂泥巖互層沉積，砂體單層厚度小，儲(chǔ)層薄，縱向疊置，橫向?qū)Ρ刃暂^差，具有橫向變化快、識(shí)別難度大的特點(diǎn)。圖7為該區(qū)沙二段砂泥巖速度統(tǒng)計(jì)結(jié)果圖，泥巖速度為2 000 m/s ～2 600 m/s，砂巖速度為2 600 m/s ～3 600 m/s，砂巖含油之后縱波速度有所下降。為此，AVO相應(yīng)表現(xiàn)為振幅隨偏移距的增大而減小，屬第一類AVO。

圖7　某工區(qū)da21井區(qū)沙二段砂、泥巖速度統(tǒng)計(jì)Fig．7　The velocity statistics of sandstone and 　mudstone of Es2 in da21 well

圖8　彈性阻抗反演得到的過da21井λρ剖面Fig．8　Theλρsection with a well location of da21

圖9　彈性阻抗反演得到的過da21井μρ剖面Fig．9　Theμρsection with a well location of da21

圖10　彈性阻抗反演得到的過da21井VP/VS剖面Fig．10　The VP/VS section with a well location of da21

5結(jié)論與認(rèn)識(shí)

彈性阻抗反演的提出，使波阻抗反演從疊后發(fā)展到疊前，可獲得多種疊后資料無法得到的信息。目前，彈性阻抗反演通常采用廣義線性反演算法，易陷入局部極優(yōu)、強(qiáng)烈依賴于初始模型。本文將粒子群算法應(yīng)用于彈性阻抗反演中，通過對(duì)每一代進(jìn)化產(chǎn)生的粒子進(jìn)行隨機(jī)攝動(dòng)，并加入模擬退火算子，改善了收斂速度和精度，用測(cè)井?dāng)?shù)據(jù)加以驗(yàn)證并應(yīng)用于勝利油田某工區(qū)實(shí)際資料的彈性阻抗反演中。提取了多種疊前地震屬性剖面。該方法提取的屬性剖面能夠反映儲(chǔ)層特征，與測(cè)井結(jié)果相符，說明基于改進(jìn)粒子群算法的彈性阻抗反演方法切實(shí)有效，可以服務(wù)于巖性識(shí)別和油氣預(yù)測(cè)。

參考文獻(xiàn)：

[1] CONNOLLY P. Elastic Impedance[J].The Leading Edge,1999,18(4):438-452.

[2] CAMBOIS G.AVO inversion and elastic impedance [J].SEG Technical Program Expanded Abstracts, 2000,70:142-145.

[3]MALLICK S.AVO and elastic impedance [J].The Leading Edge,2001,20 (10):1094-1104.

[4]WHITCOMBE D N. Elastic impedance normalization [J].Geophysics,2003,67(1):60-62.

[5]WHITCOMBE D N.Extended elastic impedance for fluid and lithology prediction [J].Geophysics,2002,67(1):63-67.

[6]馬勁風(fēng).地震勘探中廣義彈性阻抗的正反演[J].地球物理學(xué)報(bào),2003,46(3):118-124

MA J F.Forward modeling and inversion method of generalized elastic impedance in seismic exploration [J].Chinese Journal of Geophysics,2003,46(3):118-124.(In Chinese)

[7]SUBHASHIS MALLICK.A simple approximation to the P-wave reflection coefficient and its implication in the inversion of amplitude variation with offset data[J].Geophysics,1993,58(4):544-552.

[8]YANGHUA WANG.Approximations to the Zoeppritz equations and their use in AVO analysis[J].Geophysics,1999,64(6):1920-1927.

[9]倪逸.彈性阻抗計(jì)算的一種新方法[J].石油地球物理勘探,2003,38(2):147-155.

NI Y.A new method for calculation of elastic wave impedance[J].OGP,2003,38(2):147-155.(In Chinese)

[10]印興耀,袁世洪,張繁昌.從彈性波阻抗中提取巖石物性參數(shù)[C].北京：CPS/SEG 2004國際地球物理會(huì)議,2004.

YIN X Y,YUAN S H,ZHANG F C.Extraction of rock physical properties from elastic wave impedance.[C],Beijing:CPS/SEG 2004 International Geophysical Conference，2004.(In Chinese)

[11]王保麗,印興耀,張繁昌.基于Gray近似的彈性波阻抗方程及反演[J].石油地球物理勘探,2007,42(4):435-439.

WANG B l,YIN X Y,ZHANG F C.Gray approximation-based elastic wave impedance equation and inversion[J].OGP,2007,42(4):435-439.(In Chinese)

[12]楊文采.地球物理反演的理論與方法[M].北京:地質(zhì)出版社,1997.

YANG W C.Theories and methods of geophysics inversion [M].Beijing:Geological Publishing House,1997.(In Chinese)

[13]楊維,李歧強(qiáng).粒子群優(yōu)化算法綜述[J].中國工程科學(xué),2004,6(5):87-94.

YANG W,LI Q Q.Review of particle swarm optimization algorithm[J].China Engineering science,2004,6(5):87-94.(In Chinese)

[14]EBERHART R,KENNEDY J.A new optimizer using particle swarm theory[A].IEEE Proceedings of the sixth international symposium on micro machine and human science.Piscataway:IEEE Service Center,1995.39-43.

[15]KENNEDY J,EBERHART.Particle swarm optimization[A].IEEE Proceedings of international conference on neural networks[C].Piscataway:IEEE Service Center,1995.1942-1948.

[16]曾建潮,介倩,崔志華.微粒群算法[M].北京:科學(xué)出版社,2004.

ZENG J C,JIE Q,CUI Z H.Particle swarm optimization[M].Beijing:Science Press,2004.(In Chinese)

收稿日期：2015-03-30改回日期：2015-06-11

基金項(xiàng)目:國家自然科學(xué)基金(41004046,41230318)

作者簡(jiǎn)介：張進(jìn)(1978-)，男，博士，副教授，主要從事地震資料疊前反演，E-mail:zj515@ouc.edu.cn。

文章編號(hào)：1001-1749(2016)03-0353-08

中圖分類號(hào)：P 631.4

文獻(xiàn)標(biāo)志碼：A

DOI：10.3969/j.issn.1001-1749.2016.03.10

Application of particle swarm optimization in prestack elastic impedance inversion

ZHANG Jin, CHEN Xiao-qi, XING Lei, AN Ling-fang, ZHAO Xi, AN Zhen-fang

(Key Lab of Submarine Geoscience and Prospecting Techniques,Ministry of Edecation,Ocean University of China,Qingdao266100,China)

Abstract:Elastic impedance inversion is a combination technology of acoustic impedance inversion and AVO inversion, which can overcome the defects of post-stack impedance inversion and reflect the amplitude variation with offset. Elastic impedance inversion has been widely used to identify lithology and fluid. While conventional linear elastic impedance inversion methods strongly depend on the initial model and easily fall into local minimum, a nonlinear elastic impedance inversion based on improved particle swarm optimization is presented, which has been verified by model test. Finally, we applied this new inversion algorithm to the prestack seismic data in Shengli oilfield and achieved several elastic parameter profiles, which are consistent with the drilling results. The method has been proved to be a feasible and effective way for complex reservoir exploration.

Key words:elastic impedance inversion; particle swarm optimization; nonlinear inversion