混合粒子群算法在地震波阻抗反演中的應用

劉玉敏, 高松巖

(東北石油大學 電氣信息工程學院, 黑龍江 大慶 163318)

0 引 言

隨著油田儲量的日漸枯竭以及復雜的油氣藏很難被識別, 因而從地震資料中了解地下構(gòu)造的基本情況, 對地震資料進行綜合處理解釋就顯得尤為重要。地震反演技術(shù)就是按照一定的已知的地質(zhì)規(guī)律, 將地表觀測到的地震數(shù)據(jù)進行反推, 從而推知地球內(nèi)部的結(jié)構(gòu)形態(tài)和物質(zhì)成分[1]。波阻抗是重要的巖石物理參數(shù), 能較好地提供地層速度的信息, 可直接與鉆井資料對比進行儲層巖性解釋和物性分析, 因而波阻抗反演技術(shù)在油藏描述、 儲層參數(shù)估算、 油氣橫向預測等研究工作中發(fā)揮著越來越重要的作用[2]。

地震波阻抗反演是典型的非線性化問題, 傳統(tǒng)的地震波阻抗反演常用線性化算法, 如最快下降法, 共軛梯度法, 牛頓法等[3]。線性優(yōu)化方法存在的弊端就是其目標函數(shù)大部分需要可導或可微, 這就使一些函數(shù)難以被優(yōu)化、 收斂速度慢、 容易陷入局部最優(yōu)解等問題[4,5]。因此目前許多學者都在嘗試將智能非線性優(yōu)化算法應用到地震波阻抗反演中。如文獻[6]將蟻群算法應用到地震波阻抗反演中； 文獻[7]提出了層狀模型約束下的粒子群算法, 并將其應用到波阻抗反演中； 文獻[8]將免疫算法與遺傳算法相結(jié)合應用于波阻抗反演； 文獻[9]將遺傳算法和模擬退火相結(jié)合, 應用到波阻抗反演中, 取得了良好的效果； 文獻[10]將粒子群算法用于火山巖的波阻抗反演中, 并將反演結(jié)果與神經(jīng)網(wǎng)絡方法相比較, 證實了粒子群方法的有效性。

粒子群優(yōu)化算法是一種基于群體的隨機優(yōu)化技術(shù), 是典型的非線性算法[11]。該算法具有實現(xiàn)方便、 收斂速度快、 參數(shù)設置少等優(yōu)點。目前該算法已經(jīng)廣泛應用于各個領域, 并且許多學者在不斷從優(yōu)化算法結(jié)構(gòu)、 參數(shù)等方面對傳統(tǒng)粒子群進行改進。文獻[12]將粒子群算法成功用于組合優(yōu)化問題； 文獻[13]針對粒子群優(yōu)化算法在搜索的初期收斂較快, 而在后期容易陷入局部最優(yōu)的缺點, 對粒子群算法進行了優(yōu)化和仿真； 文獻[14]對傳統(tǒng)的粒子群算法的慣性權(quán)重進行了優(yōu)化, 使權(quán)值能自適應改變, 從而提高了粒子群搜索的速度； 文獻[15]采用了一種時變學習因子的方式設計了一種自組織粒子群優(yōu)化器； 文獻[16]提出了一種基于改進粒子群算法的彈性阻抗非線性反演方法, 獲得了良好的應用效果。也有很多學者將其他優(yōu)化算法與粒子群算法相結(jié)合, 互相取長補短, 優(yōu)化其性能。如與神經(jīng)網(wǎng)絡、 模擬退火、 混沌、 量子等相結(jié)合[17-20]。筆者在上述文獻的基礎上, 考慮將粒子群算法與混沌和遺傳算法相結(jié)合。混沌序列具有遍歷性, 混沌優(yōu)化方法具有全局漸進收斂、 易跳出局部極小點和收斂速度快等特點[18], 因此考慮將混沌思想引入到粒子群中。但是混沌粒子群優(yōu)化算法又具有不徹底性, 盲目隨機性, 因此將遺傳算法嵌入到粒子群算法中, 可以增強種群之間的互相交流, 提高種群的多樣性[21]。將筆者算法應用于波阻抗反演中, 通過對模型和實際數(shù)據(jù)的反演, 結(jié)果表明筆者算法具有較快的收斂速度和較高的反演精度。

1 算法描述

1.1 粒子群算法基本原理

粒子群算法(PSO: Particle Swarm Optimization)與其他進化算法一樣, 也采用群體和進化的概念。它通過個體間的協(xié)作與進化, 實現(xiàn)在復雜空間中搜索最優(yōu)解的過程。算法將每個個體看作是n維搜索空間中的一個沒有體積和重量的粒子, 每個粒子的位置都為優(yōu)化問題的一個潛在解, 粒子在解空間中飛行, 由速度決定其運動方向和所在位置, 并由目標函數(shù)為之確定一個適應值, 適應值的大小衡量粒子的優(yōu)劣。在搜索過程中, 粒子通過跟蹤本身迄今找到的最優(yōu)位置和全種群迄今找到的最優(yōu)位置不斷尋優(yōu), 經(jīng)逐代搜索直到滿足結(jié)束條件。具體過程如下[11,13]。

設在一個d維的目標搜索空間中, 有m個粒子組成一個群體。其中在第t次迭代時粒子i的位置表示為xi(t)=(xi1(t),xi2(t),…,xid(t))； 相應的飛行速度表示為vi(t)=(vi1(t),vi2(t),…,vid(t))； 粒子本身迄今搜索到的最優(yōu)解(個體極值)表示為pi(t)=(pi1(t),pi2(t),…,pid(t))； 整個粒子群到目前為止找到的最優(yōu)解(全局極值)表示為pg(t)=(pg1(t),pg2(t),…,pgd(t))。在每次迭代中, 粒子通過跟蹤這兩個極值更新自己的速度和位置。在第t+1次迭代計算時, 粒子i的第k維分量根據(jù)下列規(guī)則更新自己的速度和位置

vik(t+1)=wvik(t)+c1r1(pik(t)-xik(t))+c2r2(pgk(t)-xik(t))(1)

xik(t+1)=xik(t)+vik(t+1)i=1,2,…,m,k=1,2,…,d(2)

式(1)中,w為慣性權(quán)重, 體現(xiàn)了前一次迭代的速度對下一次迭代速度的影響。w大全局搜索能力強,w小局部搜索能力強。c1和c2為兩個學習因子, 體現(xiàn)了粒子向個體最優(yōu)解和全局最優(yōu)解學習的步長。r1和r2為兩個均勻分布在(0,1)之間的隨機數(shù)。

1.2 遺傳算法基本原理

遺傳算法(GA: Genetic Algorithms)是一種以自然選擇和基因遺傳理論為基礎的隨機優(yōu)化搜索算法[22]。算法結(jié)合了生物進化過程中的適者生存規(guī)則與群體內(nèi)部染色體的隨機信息交換機制。具體過程為： 首先用一定長度和編排規(guī)則的數(shù)字或字符將每個可能的解編碼為一個染色體個體。在搜索之初, 從搜索空間中隨機地選擇一組染色體作為初始群體； 選定目標函數(shù)計算群體中所有個體的適應度值, 以度量個體的優(yōu)劣； 根據(jù)適應度值從中選取優(yōu)秀個體作為父本, 這體現(xiàn)了適者生存的進化理論； 由交換概率挑選的每兩個父本通過將相異的部分基因進行交換, 得到新的個體, 體現(xiàn)了信息交換的思想； 為了在群體中能夠引入新的個體, 用極小的概率對個體進行變異操作, 為新個體的產(chǎn)生提供機會。算法通過不斷地使用選擇、 交叉、 變異等遺傳操作, 使獲得的染色體一代更比一代好, 最終得到適應度最佳的解[23]。

1.3 混合粒子群算法

通過對粒子群算法和遺傳算法進行分析、 比較, 發(fā)現(xiàn)這兩種算法都是在自然特性的基礎上模擬個體種群的適應性, 然后采用一定的變換規(guī)則通過搜索空間求解的。主要不同之處如下。 1) PSO算法在計算適應度后, 所有個體都被更新, 形成下一代個體, 使算法具有了良好的記憶性, 但搜索速度受到影響。而遺傳算法是通過某種選擇機制, 選出若干個樣本, 兩兩雜交生成下一代的新個體, 以前的知識隨著種群的改變被破壞, 記憶性差, 但雜交操作體現(xiàn)了算法信息交換的思想。2) 遺傳算法有變異操作, 雖然變異概率很低, 但這為新個體的產(chǎn)生提供了機會, 更容易跳出局部最優(yōu)。而粒子群算法沒有變異操作?？梢妰煞N算法基本流程類似, 但各有利弊。筆者以PSO算法為主框架, 在原粒子群更新操作的基礎上, 加入了遺傳算法中的選擇、 交叉、 變異思想。采用這種結(jié)合方式, 粒子群優(yōu)化算法可以通過學習借鑒以前的搜索經(jīng)驗, 保存所有粒子好的解的知識, 減少無效迭代的次數(shù), 提高收斂速度。也可以發(fā)揮遺傳算法在保持種群多樣性和全局搜索的優(yōu)勢。此外, 由于粒子群算法的初值是隨機生成的, 不能保證種群的多樣性和粒子搜索的遍歷性。因此在初值選取時, 還引入了混沌序列的概念?；煦缡且环N看上去沒有規(guī)律的運動, 它是存在于確定性的非線性系統(tǒng)中的類似隨機的行為。它具有遍歷性、 隨機性和規(guī)律性等特點, 能根據(jù)其自身的規(guī)律, 在一定范圍內(nèi)不重復遍歷整個狀態(tài)。因此具有易跳出局部最優(yōu)、 搜索速度快、 全局漸進收斂等優(yōu)點[24]。所以筆者在初值選取時, 引入了混沌序列代替隨機初始化。粒子的位置初值和速度初值以混沌序列形式排列, 既具有隨機性的本質(zhì), 又具有混沌特性, 從而提高了種群的多樣性和粒子搜索的遍歷性。筆者混合粒子群算法具體步驟如下。

1) 參數(shù)初始化。確定粒子群算法的種群規(guī)模m、 種群進化次數(shù)以及學習因子等參數(shù)； 利用logistic混沌映射模型產(chǎn)生d個混沌變量

(3)

(4)

2) 計算粒子群的每個粒子的適應度值, 將粒子的個體極值記為pbest, 全局極值記為gbest。

3) 根據(jù)粒子群優(yōu)化算法式(1)和式(2)更新粒子的速度和位置, 重新計算粒子的適應度值, 并更新pbest和gbest。

4) 對適應度進行排序, 并擇優(yōu)選取樣本作為下一步遺傳操作的父本。

5) 執(zhí)行遺傳算法的交叉, 變異操作, 生成新一代種群。

6) 判斷是否滿足算法的停止條件, 若滿足, 算法結(jié)束, 輸出最優(yōu)解； 否則返回2)。

混合粒子群算法的具體流程如圖1所示。

圖1 混合粒子群算法流程圖Fig.1 Flow chart of hybrid particle swarm optimization

2 粒子群波阻抗反演

2.1 地震波阻抗反演

目前常用的波阻抗反演方法基本上都是建立在褶積模型的基礎上。地震記錄用褶積可表示為

S(t)=R(t)*ω(t)+N(t)(5)

2.2 混合粒子群波阻抗反演思路

由褶積模型(5)可以看出, 地震記錄與波阻抗為非線性關系, 因此求取波阻抗的過程就是一個求解非線性問題的過程。粒子群算法是一個典型的求解非線性問題的算法, 將其應用于地震道的波阻抗反演, 問題的關鍵在于如何構(gòu)建一個能夠?qū)Ψ囱萁Y(jié)果進行評價的目標函數(shù)[4]。筆者在簡化模型的基礎上, 不考慮平衡白噪聲, 直接利用具有一定抗噪能力的最小二乘原理建立目標函數(shù), 其數(shù)學表達式如下

(6)

其中S(t)為實際地震記錄,N為采樣點數(shù),ω(t)為提取的地震子波,R(t)為所求的反射系數(shù)。

標準PSO波阻抗反演過程中, 每個粒子的位置, 代表著每個波阻抗模型。每個粒子的速度, 則代表著波阻抗模型的修正量。PSO的任意一個粒子在波阻抗空間里進行搜索, 每搜索到一個新位置就相當于構(gòu)建了一個新的波阻抗模型。將波阻抗根據(jù)式(5)與子波褶積, 即可求出模型響應S。在每步迭代過程中依據(jù)粒子群速度和位置更新公式、 模型以及模型的修正量, 并最終達到目標函數(shù)最小。

筆者混合粒子群算法反演波阻抗過程具體描述如下。

1) 確定求解模型空間。在分析地質(zhì)、 測井等先驗信息的基礎上, 建立反演初始模型。以該模型為依據(jù), 隨機產(chǎn)生一組波阻抗合理解, 構(gòu)成模型求解空間。并利用混沌發(fā)生器Logistic映射使之成為一個混沌變量, 也即粒子。粒子的維數(shù)對應地下介質(zhì)的層數(shù)。

2) 設定算法參數(shù)。為每個粒子賦一個初始速度v, 作為波阻抗模型的修正量初值。同時為慣性權(quán)重、 學習因子、 交叉變異概率等參數(shù)設置初值。

3) 計算每個粒子的適應度。目標函數(shù)值越小, 模型響應越接近實際地震記錄。記錄每次迭代粒子個體最優(yōu)位置pbest和所有粒子的最優(yōu)位置gbest。

4) 迭代尋優(yōu)。采用粒子群位置和速度公式對粒子進行更新, 并對粒子進行選擇、 交叉、 變異這3個遺傳操作。粒子通過反復追隨自身歷史最優(yōu)位置和種群最優(yōu)位置, 并以一定的概率進行信息交換, 從而逐步趨近全局最優(yōu)解模型。

5) 波阻抗的輸出。當滿足迭代終止條件(達到理想精度或者最大迭代次數(shù)時)時, 停止迭代, 輸出種群最優(yōu)位置gbest即為最終反演的波阻抗模型值。否則返回3)。

3 算例分析

3.1 標準粒子群波阻抗反演

筆者設計了一個一維的五層層狀介質(zhì)波阻抗模型, 如圖2所示。該模型各層介質(zhì)波阻抗為1.5,1.6,2.0,1.7,1.8, 單位為106g·m-2·s-1。模型空間采樣間隔為2 ms, 有60個采樣點。選擇雷克(Ricker)子波合成地震記錄, 褶積建立合成地震記錄如圖3所示。雷克(Ricker)子波是最常用的一種子波。它有一個幅值較大的主瓣, 兩個幅值較小、 對稱的旁瓣, 非常接近真實地震子波。雷克子波的數(shù)學表達式為f(t)=[1-2(?ft)2]e-(πft)2。筆者選擇的Ricker子波頻率為50 Hz, 采樣率為0.002 s。

圖2 波阻抗模型 圖3 合成地震記錄 Fig.2 Wave impedance model Fig.3 Synthetic seismogram

將標準粒子群算法用于上述模型的波阻抗反演中。反演波阻抗結(jié)果如圖4所示, 反演地震記錄如圖5所示。從反演的波阻抗和地震記錄中, 可以看出反演的波阻抗模型趨勢與實際模型一致, 基本能反映出模型的各個層位波阻抗的變化。但是反演結(jié)果與實際模型相比還是存在一定的誤差, 誤差精度在104, 誤差曲線如圖6所示。

圖4 標準粒子群波阻抗反演 圖5 標準粒子群反演地震記錄 Fig.4 Inversive wave impedance with Fig.5 Inversive seismogram with ordinary ordinary particle swarm optimization particle swarm optimization

3.2 混合粒子群算法波阻抗反演

由于標準粒子群算法具有后期收斂速度較慢, 容易早熟且易陷于局部極小值的局限性, 因此采用筆者的混合粒子群算法對上述模型進行反演。設計學習因子大小為1.494 45, 種群為60, 慣性權(quán)重采用線性遞減慣性權(quán)重ω(t)=0.9-0.5t/T, 迭代次數(shù)T為1 000次, 交叉概率為0.5, 變異概率為0.05。反演結(jié)果如圖7～圖9所示。

圖6 標準粒子群波阻抗反演誤差曲線 圖7 混合粒子群波阻抗反演結(jié)果 Fig.6 Error curve with ordinary Fig.7 Inversive wave impedance with particle swarm optimization hybrid particle swarm optimization

圖8 混合粒子群反演地震記錄 圖9 混合粒子群反演誤差曲線 Fig.8 Inversive seismogram with hybrid Fig.9 Error curve with hybrid particle swarm optimization particle swarm optimization

從圖7～圖9可見, 在相同參數(shù)設置的前提下, 混合粒子群算法反演出的波阻抗、 地震記錄與模型的吻合度更高, 雖然在部分區(qū)域存在誤差, 但在誤差量級上有了大幅度的減小。這足以說明, 筆者的混合粒子群算法對波阻抗反演有較好的效果, 精度較高, 具有一定的可信度。

3.3 算法抗噪性分析

考慮到實際環(huán)境中噪聲的存在, 筆者在理論地震記錄中加入15%的隨機噪聲, 對算法進行了抗噪性試驗與分析。

圖10 反演波阻抗(15%噪聲) 圖11 反演地震記錄(15%噪聲) Fig.10 Inversive wave impedance with 15% noise Fig.11 Inversive seismogram with 15% noise

由圖10, 圖11可以看出, 當?shù)卣鹩涗浿屑尤朐肼暫? 反演結(jié)果仍然能跟蹤模型參數(shù), 可以較真實的反映有效信息。噪聲顯然對結(jié)果造成了影響, 使某些區(qū)域反演誤差有所增加, 但還在可接受的范圍之內(nèi)。這說明筆者混合粒子群算法的波阻抗反演具有一定的抗噪能力。

4 應用實例

筆者利用混合粒子群算法對位于松遼盆地某區(qū)塊的實際地震資料進行波阻抗反演。該數(shù)據(jù)采樣間隔0.5 ms, 地震道為60道, 每道有100個采樣點, 采用40 Hz的雷克子波作為提取的子波。反演結(jié)果如圖12所示。

圖12 實際地震資料反演波阻抗剖面Fig.12 Inversive wave impedance section with seismic data

從該剖面的儲層預測結(jié)果可以看出, 目標層段內(nèi)儲層縱向各薄層砂體在波阻抗反演剖面上顯示效果良好, 預測結(jié)果未受到地震資料宏觀模型影響, 井間薄層砂體連續(xù)性關系得到了清晰刻畫, 對于油田范圍內(nèi)薄層儲層預測具有很好的參考價值。

5 結(jié) 語

筆者針對傳統(tǒng)粒子群算法在波阻抗反演中效率低, 誤差偏大的缺點, 將混沌和遺傳算法引入到粒子群中, 并對地震波阻抗模型進行了反演。結(jié)果表明, 筆者采用的混合粒子群算法在收斂精度上, 明顯優(yōu)于傳統(tǒng)粒子群反演。通過在模型中加入噪聲, 驗證了筆者的混合粒子群算法也具有良好的抗噪聲效果。利用混合粒子群算法對實際地震資料進行波阻抗反演, 獲得了殘差較小的反演地震記錄, 說明筆者方法具有一定的應用價值。