基于自回歸滑動(dòng)平均模型和粒子群算法的地震子波提取

戴永壽，牛 慧，彭 星，王少水

(中國(guó)石油大學(xué)信息與控制工程學(xué)院，山東東營(yíng) 257061)

基于自回歸滑動(dòng)平均模型和粒子群算法的地震子波提取

戴永壽，牛 慧，彭 星，王少水

(中國(guó)石油大學(xué)信息與控制工程學(xué)院，山東東營(yíng) 257061)

基于自回歸滑動(dòng)平均(ARMA)模型理論，對(duì)地震子波進(jìn)行參數(shù)化建模，采用累積量擬合法精確估計(jì)參數(shù)，使地震子波提取問(wèn)題最終歸結(jié)為一個(gè)多參數(shù)、多極值的非線性函數(shù)優(yōu)化問(wèn)題。對(duì)基本粒子群算法進(jìn)行改進(jìn)，通過(guò)自適應(yīng)參數(shù)調(diào)整和邊界約束，克服基本粒子群算法易陷入局部極值的缺陷，同時(shí)提高算法尋優(yōu)精度和計(jì)算效率。仿真數(shù)據(jù)試驗(yàn)結(jié)果驗(yàn)證了改進(jìn)的粒子群算法在地震子波提取方法中的有效性和穩(wěn)定性。

地震數(shù)據(jù)處理;自回歸滑動(dòng)平均模型;地震子波;系統(tǒng)辨識(shí);累積量擬合;粒子群算法

地震子波提取是地震資料反褶積處理、波阻抗反演和正演模型的基礎(chǔ)［1］。Lazear［2］對(duì)基于四階累積量的混合相位子波提取方法進(jìn)行研究，提出將子波的四階矩和地震數(shù)據(jù)的四階累積量進(jìn)行擬合優(yōu)化以提取子波，但由于參數(shù)模型的未知性，確定模型初始參數(shù)范圍較為困難，嚴(yán)重影響算法的尋優(yōu)效率。粒子群算法是基于群體智能的多峰值尋優(yōu)算法，概念簡(jiǎn)單且收斂速度快，但算法尋優(yōu)精度與參數(shù)選擇有著密切的關(guān)系。筆者基于自回歸滑動(dòng)平均(ARMA)模型理論，將累積量擬合法和參數(shù)化辨識(shí)方法相結(jié)合，將子波提取問(wèn)題最終歸結(jié)為一個(gè)多參數(shù)、多極值的非線性函數(shù)優(yōu)化問(wèn)題［3］。

1 地震子波ARMA模型和累積量擬合函數(shù)

通常地震激勵(lì)假設(shè)為一零均值平穩(wěn)隨機(jī)過(guò)程y(n)，用如下褶積模型［3］描述一道地震記錄：

式中，h(n)為地震子波;lc和lnc分別為子波因果、非因果部分的長(zhǎng)度;r(n)為發(fā)射系數(shù)序列;v(n)為加性噪聲。模型(1)一般有如下假設(shè)：①反射系數(shù)序列r(n)為零均值、獨(dú)立同分布的非高斯過(guò)程［4］，其方差＜∞，且至少存在一k階累積量滿足＜∞;②E(θ)≤4為環(huán)境噪聲、儀器噪聲及激發(fā)產(chǎn)生的多次反射噪聲等的合成噪聲信號(hào)，為一與r(n)統(tǒng)計(jì)獨(dú)立的隨機(jī)過(guò)程，一般為加性色噪聲［5-6］，且其高斯成分遠(yuǎn)大于非高斯部分;③h(n)為非因果、混合相位的地震子波，其非因果性表征了檢波器或信道引入的失真。

根據(jù)文獻(xiàn)［4］，將地震子波的褶積模型轉(zhuǎn)換為

其中，x(n)為平穩(wěn)隨機(jī)過(guò)程，為一ARMA模型在輸入為反射系數(shù)序列r(n)時(shí)的響應(yīng)，該ARMA模型滿足以下差分方程：

式中，p和q分別為ARMA模型的AR(自回歸)和MA(滑動(dòng)平均)部分的階數(shù)，模型對(duì)應(yīng)沖激響應(yīng)為式(1)中的地震子波h(n)。

由于奇數(shù)階累積量恒為零，偶數(shù)階累積量不為零，故應(yīng)采用四階累積量進(jìn)行地震子波估計(jì)。對(duì)地震記錄求四階累積量，應(yīng)用Bartlett-Brillinger-Rosenblatt公式，由假設(shè)條件結(jié)合累積量定義和性質(zhì)，得

式中，C4y(t1，t2，t3) 為地震數(shù)據(jù)的四階累積量;m4b(t1，t2，t3)為估計(jì)的地震子波四階矩;r4r為反射系數(shù)序列峰度。r4r為一標(biāo)量，在不影響所求取子波形態(tài)的情況下，可作為系數(shù)被m4b(t1，t2，t3)吸收，即地震記錄的四階累積量與子波的四階矩函數(shù)之間僅差一個(gè)標(biāo)量因子，可通過(guò)計(jì)算地震記錄的四階累積量來(lái)估算子波。由于高階累積量關(guān)于其變?cè)菍?duì)稱的，因此在最小平方誤差意義下可建立目標(biāo)函數(shù)［5］

式中，L為地震子波的長(zhǎng)度。目標(biāo)函數(shù)φ為地震子波的函數(shù)，當(dāng)φ達(dá)到最小值時(shí)，所對(duì)應(yīng)的一組變量即為子波的估計(jì)值。此目標(biāo)函數(shù)的求解是一個(gè)非線性多參數(shù)多極值的優(yōu)化問(wèn)題，需要應(yīng)用非線性優(yōu)化算法進(jìn)行求解。通過(guò)對(duì)優(yōu)化算法進(jìn)行研究，提出了改進(jìn)的粒子群算法，并將其運(yùn)用到模型求解當(dāng)中，以改善累積量擬合優(yōu)化的計(jì)算效率和尋優(yōu)精度。

2 改進(jìn)的粒子群算法

粒子群算法［6］概念簡(jiǎn)明，實(shí)現(xiàn)方便，參數(shù)設(shè)置少和收斂速度快，但算法精度與參數(shù)選擇有著密切的關(guān)系，不恰當(dāng)?shù)膮?shù)選擇，易導(dǎo)致算法結(jié)果發(fā)散，尋優(yōu)精度降低［7-8］。

基本粒子群算法的數(shù)學(xué)表達(dá)式如下：

式中，ω為慣性權(quán)重，控制前一代對(duì)下一代速度的影響程度;c1和c2為學(xué)習(xí)因子;r1和r2為在［0，1］范圍內(nèi)服從均勻分布的隨機(jī)變量;t為迭代次數(shù)。

基本粒子群算法存在如下缺陷［9］：①算法精度同參數(shù)選擇有著密切的關(guān)系，若參數(shù)選擇不當(dāng)，算法精度大大降低;②在整個(gè)實(shí)現(xiàn)過(guò)程中，參數(shù)均為定值，而不考慮具體的優(yōu)化模型和迭代過(guò)程，影響了算法精度。為此，本文中提出一種改進(jìn)的粒子群算法：令慣性權(quán)重ω、學(xué)習(xí)因子c1和c2隨著迭代次數(shù)的增加線性遞減或線性遞加，從而保證算法初期個(gè)體能搜索整個(gè)空間而不陷入局部最優(yōu)值，后期能朝著全局最優(yōu)值收斂而找到全局最優(yōu)值;為避免粒子位置超出給定位置，對(duì)粒子邊界條件進(jìn)行設(shè)置，保證尋優(yōu)解在有效的范圍內(nèi)，從而應(yīng)用到地震子波提取中。

2.1 慣性權(quán)重

Shi和 Eberhart［10］提出隨著迭代次數(shù)增加，慣性權(quán)重線性遞減的策略，表達(dá)如下：

式中，ωmax為最大慣性權(quán)重;ωmin為最小慣性權(quán)重;ts為算法迭代總次數(shù);t*為算法當(dāng)前迭代次數(shù)。

2.2 學(xué)習(xí)因子

c1和c2分別表征搜索局部極值和全局極值的能力。算法初期，較大的c1和較小的c2使粒子在整個(gè)可行空間進(jìn)行搜索，保證地震子波初始范圍的多樣性;算法后期，較小的c1和較大的c2使粒子快速收斂于全局最優(yōu)值，找到地震子波的準(zhǔn)確解。根據(jù)分析，數(shù)學(xué)表達(dá)式如下：

式中，c1start和c2start為學(xué)習(xí)因子的初始值;c1end和c2end表示兩個(gè)學(xué)習(xí)因子的終值。

2.3 粒子邊界條件

當(dāng)粒子的位置超出了給定位置時(shí)，重新設(shè)置粒子的當(dāng)前位置能夠使逃逸出去的粒子又回到種群中，而且其值也不再單一地確定為邊界值，確保了如果粒子過(guò)多發(fā)生逃逸時(shí)新種群的多樣性，保證地震子波提取的有效性。

3 地震子波提取的具體實(shí)現(xiàn)步驟

(1)數(shù)據(jù)生成。用待估計(jì)ARMA子波模型與滿足反射系數(shù)序列假設(shè)的隨機(jī)序列合成地震數(shù)據(jù)。

(2)參數(shù)預(yù)估計(jì)。運(yùn)用奇異值分解(SVD)法、Zhang算法、基于奇異值分解的總體最小二乘法(SVD-TLS)等模型辨識(shí)方法確定待估計(jì)模型的階數(shù)p、q和自回歸參數(shù)、滑動(dòng)平均參數(shù)。

(3)確定參數(shù)搜索空間。在初步估計(jì)模型參數(shù)的基礎(chǔ)上，確定模型參數(shù)向量θ的搜索范圍。

(4)擬合優(yōu)化。在給定參數(shù)向量的搜索范圍內(nèi)用改進(jìn)的粒子群算法對(duì)累積量擬合目標(biāo)函數(shù)E(θ)進(jìn)行參數(shù)精確估計(jì)，尋找最優(yōu)解x。

(5)評(píng)價(jià)函數(shù)分析。計(jì)算最優(yōu)解x的評(píng)價(jià)函數(shù)值，當(dāng)其明顯增加或降低時(shí)繼續(xù)，否則轉(zhuǎn)至步驟(7)。

(6)階數(shù)擾動(dòng)。根據(jù)擬合誤差，調(diào)整上述方法的閾值，對(duì)模型的階數(shù)進(jìn)行擾動(dòng)，生成一組新的p、q值，然后轉(zhuǎn)至步驟(3)。

(7)結(jié)束。將當(dāng)次擾動(dòng)前所得階數(shù)p、q及精確參數(shù)估計(jì)值x視為模型最優(yōu)解。

4 合成地震記錄處理

通過(guò)合成地震記錄進(jìn)行仿真試驗(yàn)來(lái)驗(yàn)證方法的有效性和實(shí)用性。合成地震道所需的地層反射系數(shù)序列為獨(dú)立同分布(IID)的隨機(jī)過(guò)程，且服從伯努力分布;地震子波為因果混合相位的，其ARMA模型的差分方程為

地震子波ARMA模型的零極點(diǎn)分布如圖1所示。圖中兩個(gè)極點(diǎn)全部在單位圓內(nèi)，根據(jù)子波分類(lèi)標(biāo)準(zhǔn)可知，子波是因果的。同理，圖中兩個(gè)零點(diǎn)，一個(gè)在單位圓內(nèi)，一個(gè)在單位圓外，子波為混合相位的，即子波為因果混合相位的(子波波形見(jiàn)圖2)。

地震記錄中的噪聲一般假設(shè)為高斯有色噪聲，若地震記錄表示為y(t)=s(t)+v(t)，其中s(t)為有效波，v(t) 為噪聲，定義噪信比 RNS［11］為

在頻域定義噪信比為

式中，ω為地震記錄在頻域中的頻率;Avv(w)和Ass(w)分別為噪聲和有效波的振幅譜;Rvv(w)和Rss(w)分別為噪聲和有效波的功率譜。本文中采用式(13)形式定義噪信比。

采用基于相關(guān)函數(shù)的SVD-TLS算法估計(jì)AR參數(shù)，同時(shí)采用累積量法對(duì)MA參數(shù)進(jìn)行估計(jì)，建立初始子波尋優(yōu)范圍，并用改進(jìn)的粒子群算法進(jìn)行優(yōu)化，在不同數(shù)據(jù)長(zhǎng)度、噪聲的情況下，對(duì)該記錄運(yùn)用本文方法進(jìn)行子波提取，驗(yàn)證方法的有效性和實(shí)用性。

4.1 數(shù)據(jù)長(zhǎng)度對(duì)子波ARMA模型參數(shù)估計(jì)的影響

每次試驗(yàn)均在噪信能量比為50%的情況下，隨機(jī)生成長(zhǎng)度分別為10、8、5和3 s的合成地震記錄。

由表1可知，AR和MA參數(shù)在長(zhǎng)數(shù)據(jù)情況下均能得到較好的估計(jì)結(jié)果，E(θ)也滿足E(θ)≤4的標(biāo)準(zhǔn)，而在短數(shù)據(jù)(5 s和3 s)時(shí)估計(jì)結(jié)果與實(shí)際值有一定偏差，不過(guò)在不同數(shù)據(jù)長(zhǎng)度情況下提取的子波波形具有較好的連續(xù)性，如圖3(a)所示。

4.2 環(huán)境噪聲對(duì)子波ARMA模型參數(shù)估計(jì)的影響

為考察環(huán)境噪聲對(duì)模型參數(shù)估計(jì)方法的影響，每次試驗(yàn)均生成長(zhǎng)度為5 s的合成地震記錄，對(duì)該記錄分別加入噪信能量比為20%、50%、100%、200%的高斯色噪聲后，利用本文方法對(duì)子波模型進(jìn)行參數(shù)估計(jì)，所得結(jié)果見(jiàn)表2，提取的子波波形見(jiàn)圖3(b)。由表2和圖3(b)可知，隨著噪聲強(qiáng)度的增加，模型參數(shù)估計(jì)誤差略有增加，即使噪聲與信號(hào)的能量比達(dá)到1∶1，模型參數(shù)估計(jì)仍得到很好的結(jié)果，驗(yàn)證了模型參數(shù)估計(jì)和粒子群算法既有很好的抗噪容噪能力，又能在適宜的噪聲環(huán)境中應(yīng)用，也說(shuō)明了高階累積量的適用性。

表1 不同數(shù)據(jù)長(zhǎng)度下地震子波參數(shù)估計(jì)結(jié)果Table 1 Seismic wavelet parameter estimation results with different length

圖3 數(shù)據(jù)長(zhǎng)度和高斯色噪聲對(duì)提取的子波波形的影響Fig.3 Effects of trace length and Gaussian noise on waveforms of extracted wavelets

表2 不同強(qiáng)度高斯色噪聲環(huán)境下地震子波參數(shù)估計(jì)結(jié)果Table 2 Seismic wavelet parameter estimation results with different noise-signal-ratio

5 結(jié)論與建議

(1)用改進(jìn)的粒子群算法進(jìn)行尋優(yōu)，假設(shè)地震子波為因果混合相位時(shí)的仿真結(jié)果表明，所提方法能夠有效地確定模型參數(shù)，在不同數(shù)據(jù)長(zhǎng)度和適宜噪聲環(huán)境下仍能得到很好的結(jié)果。

(2)改進(jìn)的粒子群算法提高了算法收斂的穩(wěn)定性和收斂質(zhì)量，但是迭代次數(shù)較高，一般迭代到80～100次時(shí)才會(huì)滿足目標(biāo)函數(shù)E(θ)≤4的條件。因此，如何在保證地震子波提取精度的同時(shí)，提高粒子群算法尋優(yōu)的效率，將是以后研究的重點(diǎn)。

(3)實(shí)際地震記錄是非因果、混合相位的，因此下一步的工作目標(biāo)是采用改進(jìn)的粒子群算法對(duì)非因果混合相位子波和實(shí)際地震記錄進(jìn)行優(yōu)化，并與前期的優(yōu)化算法，如遺傳算法等優(yōu)化結(jié)果進(jìn)行比較，以得到較好的優(yōu)化結(jié)果。

［1］ROBINSON E A，TREITEL Sven.Geophysical signal analysis［M］.Englewood Cliffs：Prentice-Hall，1980.

［2］LAZEAR G D.Mixed-phase wavelet estimation using fourthorder cumulant［J］.Geophysics，1993，58：1042-1051.

［3］師學(xué)明，王家映，易遠(yuǎn)元，等.一種新的地球物理反演方法——模擬原子躍遷反演法［J］.地球物理學(xué)報(bào)，2007，50(1)：305-312.

SHI Xue-ming，WANG Jia-ying，YI Yuan-yuan，et al.A study on the simulated atomic transition algorithm for geophysical inversion［J］.Chinese Journal Geophysics，2007，50(1)：305-312.

［4］ 戴永壽，王俊嶺，王偉偉，等.基于高階累積量ARMA模型線性非線性結(jié)合的地震子波提取方法研究［J］.地球物理學(xué)報(bào)，2008，51(6)：1851-1859.

DAI Yong-Shou，WANG Jun-ling，WANG Wei-wei，et al.Seismic wavelet extraction via cumulant-based ARMA model approach with linear and nonlinear combination［J］.Chinese Journal Geophysics，2008，51(6)：1851-1859.

［5］ADNAN Al-Smadi.Cumulant-based order selection of non-Gaussian autoregressive moving average models：the corner method［J］.Signal Processing，2005，85(3)：449-456.

［6］KENNEDY J，EBERHART C.Particle swarm optimization：proceeding of IEEE International Conference on Neural Networks［C］.Piscataway，NJ：IEEE，c1995：1942-1948.

［7］KENNEDY J，EBERHART R C.A new optimizer using particle swarm theory：proceedings of the Sixth International Symposium on Micro Machine and Human Science［C］.Nagoya：IEEE，c1995：1045-1049.

［8］NIU Hui，DAI Yongshou，PENG Xing.Particle swarm optimization with adaptive parameters and boundary constraints：proceedings of the 3rd International Conference on Networks Security，Wireless Communications and Trusted Computing［C］.Wuhan：IEEE，c2011：1082-1085.

［9］CLERC M，KENNEDY J.The particle swarm-explosion，stability，and convergence in a multidimensional complex space［J］.IEEE Transactions on Evolutionary Computation，2002，6(1)：58-73

［10］KANNAN S，SLOCHANAL S M R，SUBBARAJ P，et al.Application of particle swarm optimization technique and its variants to generation expansion planning problem［J］.Electric Power Systems Research，2004，70(3)：203-210.

［11］YIN P Y.A discrete particle swarm algorithm for optimal polygonal approximation of digital curves［J］.Journal of Visual Communication and Image Representation，2004，15(2)：241-260.

Seismic wavelet extraction based on auto-regressive and moving average model and particle swarm optimization

DAI Yong-shou，NIU Hui，PENG Xing，WANG Shao-shui
(College of Information and Control Engineering in China University of Petroleum，Dongying 257061，China)

A seismic wavelet parametric model was developed based on auto-regressive and moving average(ARMA)model theory.The model parameters were accurately determined based on cumulant fitting method.So the seismic wavelet can be a multi-parameters，multi-extremes nonlinear functional optimization problem.An improved particle swarm optimization with adaptive parameters and boundary constraints was proposed for the local extreme value defects of elementary particle swarm optimization.The optimization accuracy and computation efficiency are also improved.Simulation results show that the method has good applicability and stability in seismic wavelet extraction.

seismic data processing;auto-regressive and moving average(ARMA)model;seismic wavelet;system identification;cumulant fitting;particle swarm optimization

P 631

A

10.3969/j.issn.1673-5005.2011.03.009

1673-5005(2011)03-0047-04

2011－04－02

國(guó)家自然科學(xué)基金項(xiàng)目(40974072);山東省自然科學(xué)基金項(xiàng)目(ZR2010DM14)

戴永壽(1963－)，男(漢族)，安徽巢湖人，教授，博士，博士生導(dǎo)師，從事信號(hào)與信息處理和計(jì)算機(jī)測(cè)控技術(shù)等教學(xué)研究工作。

(編輯 修榮榮)