張亞南，戴永壽，丁進(jìn)杰，張漫漫，王蓉蓉

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

傳統(tǒng)的地震反褶積處理常常將地震子波假設(shè)為因果最小相位的，但在實(shí)際應(yīng)用過程中，地震波在地層中傳播受各種復(fù)雜地質(zhì)因素的影響，子波的最小相位假設(shè)很難滿足;且由于接收檢波器等儀器性能的限制將非因果成分引入到地震子波中，地震子波的因果性也不能嚴(yán)格保證。故傳統(tǒng)反褶積結(jié)果中往往殘留子波相位成分，從而導(dǎo)致分辨率無(wú)顯著提高或者地層假象的出現(xiàn)。

為解決最小相位反褶積后相位殘余問題，許多學(xué)者假設(shè)子波振幅已經(jīng)準(zhǔn)確提出，對(duì)反褶積后的剖面應(yīng)用相位校正方法。周興元[1]在最小相位反褶積處理后的剖面中對(duì)子波剩余相位進(jìn)行常相位校正，一定程度上提高了剖面的信噪比和分辨率。郭向宇等[2]提出了一種估算剩余混合相位子波相位的算法，對(duì)地震記錄中的子波進(jìn)行相位校正，可使子波接近或達(dá)到零相位，從而消除剩余相位的影響，提高了地震記錄的分辨率。王君等[3]提出了一種在最小相位反褶積后估算子波相位補(bǔ)償因子的方法，并對(duì)地震記錄進(jìn)行相位校正，可使剩余相位接近或達(dá)到零相位，從而提高資料的分辨率。但以上方法由于假設(shè)條件的苛刻性、模型描述的局限性和選取的補(bǔ)償因子的特殊性，其效果并不十分理想，且不能涵蓋所有因果和相位性的地震子波相位殘余。

我們采用自回歸滑動(dòng)平均(ARMA)模型和單位化振幅譜的方式來(lái)構(gòu)造純相位濾波器，以描述任意因果和任意特性的地震子波相位殘余，并在最大方差模準(zhǔn)則約束下，采用改進(jìn)的粒子群算法對(duì)子波相位殘余進(jìn)行非線性尋優(yōu)，最終實(shí)現(xiàn)子波相位的精確校正。通過正演模擬數(shù)據(jù)和實(shí)際地震資料處理結(jié)果驗(yàn)證了所提方法的有效性和實(shí)用性。

1 純相位濾波器的構(gòu)造

在子波提取的過程中往往假設(shè)子波為最小相位，故反褶積結(jié)果中子波殘余相位有可能是非因果的。為準(zhǔn)確消除子波殘余相位影響，我們采用單位化ARMA模型振幅譜構(gòu)造純相位濾波器，以描述任意特性的子波殘余相位。相對(duì)于其它參數(shù)化模型，ARMA模型具有參數(shù)吝嗇的特點(diǎn)，能夠采用較少的參數(shù)描述一個(gè)精確的濾波器。ARMA模型可以表示為

(1)

式中：x(n)為濾波器輸出;e(n)為濾波器輸入信號(hào);自回歸(AR)部分的參數(shù)ai和滑動(dòng)平均(MA)部分的參數(shù)bk均為實(shí)數(shù);p為AR階數(shù)，q為MA階數(shù)，不失一般性，令a0=b0=1。若將地震子波視為一個(gè)濾波器系統(tǒng)，Robinson褶積模型[4]也可由(1)式表示。

假設(shè)ARMA描述下濾波器系統(tǒng)z域變換的系統(tǒng)函數(shù)W(z)無(wú)零極點(diǎn)對(duì)消，則

(2)

式中：w(m)為濾波器響應(yīng)的時(shí)間域序列；W(z)為其z域表示。由于ai和bk均為實(shí)系數(shù)，因此系統(tǒng)零點(diǎn)ck與極點(diǎn)di均為實(shí)根或共軛復(fù)根。

(2)式所示的ARMA濾波器模型中，當(dāng)零點(diǎn)全部在單位圓內(nèi)時(shí)為零相位濾波器，內(nèi)外兼有為混合相位濾波器，全部在單位圓上及單位圓外時(shí)為最大相位濾波器;當(dāng)極點(diǎn)全部在單位圓內(nèi)時(shí)為因果濾波器，內(nèi)外兼有為混合因果濾波器，全部在單位圓上及單位圓外時(shí)為反因果濾波器。故ARMA模型可以描述任意因果特性、任意相位特性的濾波器系統(tǒng)。

令z=ejω，(2)式中濾波器的z域表示轉(zhuǎn)換為頻率響應(yīng)，即

(3)

純相位濾波器的時(shí)域?yàn)榧兿辔恍蛄校媪⑷A[5]在其研究純相位序列的能量傳遞性質(zhì)的過程中，給出了純相位序列的驗(yàn)證方法：一個(gè)序列g(shù)(n)為純相位序列的充要條件是其z變換G(z)z=ejω能表示成ejφ(ω)的形式。則根據(jù)(3)式可構(gòu)造純相位濾波器：

(4)

由(4)式可以看出，任意濾波器均可以通過頻率譜單位化來(lái)構(gòu)造純相位濾波器。為避免計(jì)算錯(cuò)誤，當(dāng)(4)式分母中某些頻率成分為零值時(shí)，相除結(jié)果規(guī)定為單位值1。

相對(duì)于其它模型，由ARMA模型單位化振幅譜構(gòu)造出的純相位濾波器可以描述任意特性的子波相位殘余，包括最小相位、混合相位、最大相位，以及因果、混合因果、反因果等。

2 子波相位校正

反褶積結(jié)果中的相位殘余可表示為

(5)

式中：D(ejω)和R(ejω)分別為反褶積結(jié)果和反射系數(shù)序列的頻域表示;φR(ω)為子波殘余的相位譜。采用(4)式中構(gòu)造的純相位濾波器對(duì)(5)式進(jìn)行相位校正處理，

(6)

式中：y(n)和Y(ejω)分別為相位校正結(jié)果的時(shí)間域和頻率域表示;F-1為逆傅里葉變換算子。當(dāng)構(gòu)造的純相位濾波器和子波相位殘余具有相同相位譜時(shí)，(6)式可表示為

(7)

此時(shí)地震子波被校正為零相位，反褶積結(jié)果被校正為y0(n)，即反射系數(shù)序列r(n)。

若反射系數(shù)序列為稀疏序列，則相位校正結(jié)果具有最大的縱向分辨率，本文采用輸出結(jié)果的最大方差模準(zhǔn)則[6]來(lái)檢驗(yàn)相位校正效果。最大方差模準(zhǔn)則表示為

(8)

其中，L為相位校正輸出的數(shù)據(jù)長(zhǎng)度。采用最大方差模準(zhǔn)則Ψ(y)對(duì)相位校正結(jié)果y(n)進(jìn)行評(píng)價(jià)，當(dāng)Ψ(y)取最大值時(shí)，可以認(rèn)為構(gòu)造出的純相位濾波器的相位譜最接近于子波殘余的相位譜。

對(duì)目標(biāo)函數(shù)Ψ(y)求取最大值的過程是一個(gè)多參數(shù)、多極值的非線性優(yōu)化問題，對(duì)該目標(biāo)函數(shù)應(yīng)用的優(yōu)化算法需具有參數(shù)向量整體的全局隨機(jī)搜索能力和對(duì)單個(gè)參數(shù)的深度搜索能力。楊文采[7]、師學(xué)明等[8-9]針對(duì)地球物理學(xué)中的非線性尋優(yōu)問題引入了多種非線性優(yōu)化算法，并對(duì)此類算法進(jìn)行了分析和探討[10-11]。在此我們采用改進(jìn)的粒子群算法對(duì)參數(shù)進(jìn)行尋優(yōu)，改進(jìn)的粒子群算法繼承了基本粒子群算法的優(yōu)點(diǎn)，并具有全局尋優(yōu)能力強(qiáng)、收斂速度快、可避免陷入局部最優(yōu)值和自適應(yīng)學(xué)習(xí)的優(yōu)點(diǎn)[12]。

基于純相位濾波器的子波相位校正方法的實(shí)現(xiàn)步驟為：

1) 采用ARMA模型描述子波殘余，并且通過相關(guān)函數(shù)法確定子波殘余的AR階數(shù)p，采用信息量準(zhǔn)則和高階累積量相結(jié)合的方法確定子波殘余的MA階數(shù)q[13];

2) 初始化AR參數(shù)和MA參數(shù);

3) 通過(4)式構(gòu)造純相位濾波器，通過(6)式進(jìn)行相位校正;

4) 通過(8)式的最大方差模準(zhǔn)則對(duì)相位校正結(jié)果進(jìn)行評(píng)價(jià)，將評(píng)價(jià)值與歷史最優(yōu)進(jìn)行比較，若優(yōu)于歷史最優(yōu)，則將當(dāng)前參數(shù)和評(píng)價(jià)值設(shè)定為歷史最優(yōu);

5) 判斷是否達(dá)到最大迭代次數(shù)，若達(dá)到則轉(zhuǎn)向6);若未達(dá)到，則粒子飛翔，生成新的參數(shù)，轉(zhuǎn)向3);

6) 讀取歷史最優(yōu)參數(shù)，算法結(jié)束。

3 正演模擬數(shù)據(jù)測(cè)試

為驗(yàn)證本文所提方法的有效性，首先對(duì)合成地震記錄進(jìn)行子波相位校正測(cè)試。合成地震記錄中反射系數(shù)序列為獨(dú)立同分布(IID)的隨機(jī)過程，且服從Bernoulli-Gaussian分布，序列采樣間隔為1ms，長(zhǎng)度為1000ms，如圖1a所示。為了體現(xiàn)ARMA模型描述下子波的一般性，選用混合因果混合相位的地震子波合成地震記錄。合成地震記錄所采用的原始地震子波在ARMA模型描述下差分方程形式為

x(t)-4.02x(t-1)+8.43x(t-2)-8.15x(t-3)+2.86x(t-4)=

r(t)-0.8r(t-1)+0.2r(t-2)-0.82r(t-3)

(11)

其z域系統(tǒng)函數(shù)形式為

(12)

子波的時(shí)間域波形、合成地震記錄和子波相位譜如圖1b,圖1c和圖1d所示，其中相位譜的取值范圍為(-π,π]。

針對(duì)合成地震記錄進(jìn)行地震子波提取，在模擬過程中構(gòu)造與原始合成地震子波具有相同振幅譜的因果最小相位地震子波進(jìn)行反褶積處理。反褶積處理所用最小相位地震子波時(shí)間域波形和相位譜如圖2a和圖2b所示。應(yīng)用構(gòu)造的子波進(jìn)行反褶積處理，反褶積結(jié)果如圖2c所示，反褶積結(jié)果中所含有的子波相位殘余如圖2d所示。

采用本文所提方法構(gòu)造純相位濾波器，并采用最大方差模準(zhǔn)則進(jìn)行約束，通過改進(jìn)的粒子群算法針對(duì)圖2c中的反褶積結(jié)果進(jìn)行相位校正，相位校正結(jié)果如圖3所示。

圖1 合成地震記錄a 反射系數(shù)序列; b 原始地震子波; c 合成的地震記錄; d 原始地震子波相位譜

圖2 子波估計(jì)與反褶積結(jié)果a 構(gòu)造的因果最小相位地震子波時(shí)域波形; b 構(gòu)造的因果最小相位地震子波相位譜; c 反褶積結(jié)果; d 子波相位殘余

圖3 子波相位校正結(jié)果a 提取的子波相位殘余; b 子波相位校正結(jié)果

對(duì)比圖2d與圖3a可以看出，提取的子波相位殘余和理論子波相位殘余略有差異，但整體相似度較高，可以認(rèn)為本文所提方法能夠較為準(zhǔn)確地提取出子波相位殘余。對(duì)比圖3b與圖1a可以看出，由于提取的子波相位殘余與理論子波相位殘余的差異，子波相位校正結(jié)果與原始反射系數(shù)序列不能完全吻合，但二者僅在細(xì)節(jié)方面存在差異，整體吻合度較高，子波相位校正結(jié)果能夠準(zhǔn)確反映出原始反射系數(shù)序列的各個(gè)反射界面。模擬數(shù)據(jù)測(cè)試結(jié)果表明本文所提方法能夠?qū)崿F(xiàn)子波相位的校正，可以較為準(zhǔn)確地恢復(fù)出原始的反射系數(shù)序列，證明該方法是有效可行的。

4 實(shí)際地震資料處理

為進(jìn)一步驗(yàn)證本文所提方法對(duì)子波相位校正和提高地震剖面分辨率的作用，應(yīng)用該方法對(duì)實(shí)際地震資料進(jìn)行相位校正處理。圖4為勝利油田某區(qū)塊經(jīng)反褶積處理后的地震剖面，圖5為應(yīng)用本文方法進(jìn)行相位校正后的反褶積剖面。

由圖4和圖5對(duì)比可以看出，經(jīng)相位校正后的剖面質(zhì)量得到明顯改善：同相軸連續(xù)性變好，一些弱反射變得清晰?？梢哉J(rèn)為，本文所提基于純相位濾波器的子波相位校正方法可以實(shí)現(xiàn)提高反褶積后剖面分辨率的作用，具有一定的實(shí)際應(yīng)用價(jià)值。

圖4 勝利油田某區(qū)塊經(jīng)反褶積處理后的剖面

圖5 勝利油田某區(qū)塊經(jīng)本文方法相位校正后的反褶積剖面

5 結(jié)束語(yǔ)

我們提出了一種基于純相位濾波器的地震子波相位校正方法，通過ARMA模型和振幅譜單位化的方式構(gòu)造純相位濾波器以描述地震子波相位殘余，在最大方差模準(zhǔn)則約束下采用改進(jìn)的粒子群算法針對(duì)殘余的子波相位進(jìn)行準(zhǔn)確尋優(yōu)，以實(shí)現(xiàn)反褶積結(jié)果中子波殘余相位的校正。正演模擬數(shù)據(jù)和實(shí)際地震資料的應(yīng)用結(jié)果表明該方法是有效的，并具有一定的實(shí)際應(yīng)用價(jià)值。

將反褶積中子波殘余相位校正問題轉(zhuǎn)化為信號(hào)的智能尋優(yōu)問題是提高地震記錄分辨率的一種新途徑，下一步研究重點(diǎn)是進(jìn)一步提高相位譜的尋優(yōu)精度。

參 考 文 獻(xiàn)

[1] 周興元.常相位校正[J].石油地球物理勘探,1989,24(2):119-129

Zhou X Y.Constant phase correction [J].Oil Geophysical Prospecting,1989,24(2):119-129

[2] 郭向宇,周興元,董敏煜.混合相位子波的相位估算及校正[J].石油地球物理勘探,1998,33(2):214-221

Guo X Y,Zhou X Y,Dong M Y.Estimation and correction of mixed-phase wavelet phase [J].Oil Geophysical Prospecting,1998,33(2):214-221

[3] 王君,周興元.子波相位校正及效果分析[J].石油地球物理勘探,2008,43(增刊2):142-145

Wang J.Zhou X Y.Wavelet phase correction and analysis of effect [J].Oil Geophysical Prospecting,2008,43 (S2):142-145

[4] Robinson E A.Predictive decomposition of time series with application to seismic exploration [J].Geophysics,1967,32(3):418-484

[5] 舒立華.純相位序列的能量傳遞性質(zhì)[J].數(shù)學(xué)學(xué)報(bào),1974,17(1):20-27

Shu L H.The energy transfer properties of phase-only sequence[J].Acta Mathematica Sinica,1974,17(1):20-27

[6] 王有新,周興元.最小熵反褶積和集約化準(zhǔn)則[J].石油地球物理勘探,1993,28(6):678-684

Wang Y X,Zhou X Y.Minimum entropy deconvolution and parsimonious criterion [J].Oil Geophysical Prospecting,1993,28(6):678-684

[7] 楊文采.地球物理反演的遺傳算法[J].石油物探,1995,34(1):116-122

Yang W C.Genetic algorithm for geophysical inversion [J].Geophysical Prospecting for Petroleum,1995,34(1):116-122

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

Shi X M,Wang J Y,Yi Y Y,et a1.A study on the simulated atomic transition algorithm for geophysical inversion [J].Chinese Journal of Geophysics,2007,50(1):305-312

[9] 師學(xué)明,王家映,張勝業(yè),等.多尺度逐次逼近遺傳算法反演大地電磁資料[J].地球物理學(xué)報(bào),2000,43(1):122-130

Shi X M,Wang J Y,Zhang S Y,et a1.Multi-scale Genetic Algorithm and its application in magnetotelluric sounding data inversion [J].Chinese Journal of Geophysics,2000,43(1):122-130

[10] 楊文采.非線性地球物理反演方法:回顧與展望[J].地球物理學(xué)進(jìn)展,2002,17(2):255-261

Yang W C.Non linear geophysical inversion methods:review and perspective [J].Process in Geophysics,2002,17(2):255-261

[11] 王家映.地球物理資料非線性反演方法講座(一):地球物理反演問題概述[J].工程地球物理學(xué)報(bào),2007,4(1):1-3

Wang J Y.Lecture on non-linear inverse methods in geophysics(1):introduction to geophysical inverse problems [J].Chinese Journal of Engineering Geophysics,2007,4(1):1-3

[12] 戴永壽,?；?彭星,等.基于自回歸滑動(dòng)平均模型和粒子群算法的地震子波提取[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,35(3):47-50

Dai Y S,Niu H,Peng X,et al.Seismic wavelet extraction based on auto-regressive and moving average model and particle swarm optimization [J].Journal Of China University of Petroleum (Edition Of Natural Science),2011,35(3):47-50

[13] 張亞南,戴永壽,王少水,等.高效ARMA模型高分辨率地震子波提取方法[J].石油地球物理勘探,2011,46(5):686-694

Zhang Y N,Dai Y S,Wang S S,et al.High resolution wavelet estimation by ARMA modeling [J].Oil Geophysical Prospecting,2011,46(5):686-694