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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

基本粒子群算法存在如下缺陷［9］：①算法精度同參數(shù)選擇有著密切的關(guān)系，若參數(shù)選擇不當(dāng)，算法精度大大降低;②在整個實現(xiàn)過程中，參數(shù)均為定值，而不考慮具體的優(yōu)化模型和迭代過程，影響了算法精度。為此，本文中提出一種改進(jìn)的粒子群算法：令慣性權(quán)重ω、學(xué)習(xí)因子c1和c2隨著迭代次數(shù)的增加線性遞減或線性遞加，從而保證算法初期個體能搜索整個空間而不陷入局部最優(yōu)值，后期能朝著全局最優(yōu)值收斂而找到全局最優(yōu)值;為避免粒子位置超出給定位置，對粒子邊界條件進(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使粒子在整個可行空間進(jìn)行搜索，保證地震子波初始范圍的多樣性;算法后期，較小的c1和較大的c2使粒子快速收斂于全局最優(yōu)值，找到地震子波的準(zhǔn)確解。根據(jù)分析，數(shù)學(xué)表達(dá)式如下：

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

2.3 粒子邊界條件

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

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

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

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

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

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

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

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

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

4 合成地震記錄處理

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

地震子波ARMA模型的零極點分布如圖1所示。圖中兩個極點全部在單位圓內(nèi)，根據(jù)子波分類標(biāo)準(zhǔn)可知，子波是因果的。同理，圖中兩個零點，一個在單位圓內(nèi)，一個在單位圓外，子波為混合相位的，即子波為因果混合相位的(子波波形見圖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算法估計AR參數(shù)，同時采用累積量法對MA參數(shù)進(jìn)行估計，建立初始子波尋優(yōu)范圍，并用改進(jìn)的粒子群算法進(jìn)行優(yōu)化，在不同數(shù)據(jù)長度、噪聲的情況下，對該記錄運用本文方法進(jìn)行子波提取，驗證方法的有效性和實用性。

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

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

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

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

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

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

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

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

5 結(jié)論與建議

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

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

(3)實際地震記錄是非因果、混合相位的，因此下一步的工作目標(biāo)是采用改進(jìn)的粒子群算法對非因果混合相位子波和實際地震記錄進(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é)報，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é)報，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

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

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

(編輯 修榮榮)