基于變時窗Gabor變換的高分辨率處理技術(shù)研究與應(yīng)用＊

王清振 姜秀娣 翁 斌 劉志鵬 印海燕 陳劍軍

（中海油研究總院 北京 100028）

王清振，姜秀娣，翁斌，等.基于變時窗Gabor變換的高分辨率處理技術(shù)研究與應(yīng)用［J］.中國海上油氣，2015，27（6）：19-26.

目前常用的提高地震資料縱向分辨率的方法，如最小平方反褶積、預(yù)測反褶積、變模反褶積、同態(tài)反褶積、最小熵反褶積、基于互信息率準(zhǔn)則的盲反褶積以及譜白化等，其理論基礎(chǔ)是傳統(tǒng)的褶積模型［1-3］，該模型的基本假設(shè)之一就是子波在地下傳播過程中不隨時間變化。然而，實(shí)際地震子波常常是時變的，這就使得以該模型為理論基礎(chǔ)的提高分辨率的方法在許多情況下難以取得好的效果。為此，一些學(xué)者提出了反射地震記錄的另一種模型［4-5］，認(rèn)為子波在地下傳播過程中隨時間發(fā)生變化，不同時刻到達(dá)檢波器的子波的波形是不同的，反射地震記錄是這些具有不同到達(dá)時間的子波的疊加，因此這種模型被稱為反射地震記錄變子波模型。

地震記錄的子波時變性主要是由波前擴(kuò)散和地層吸收衰減效應(yīng)引起的［6］，其中波前擴(kuò)散可以用幾何擴(kuò)散函數(shù)來校正，而地層吸收衰減引起的子波時變效應(yīng)（即所謂的Q效應(yīng)）常采用反Q濾波和時變譜白化方法進(jìn)行補(bǔ)償［7-8］。但是，地下介質(zhì)構(gòu)造復(fù)雜導(dǎo)致Q值往往難以求準(zhǔn)，從而影響反Q濾波的準(zhǔn)確程度，而時變譜白化方法的參數(shù)需要人為選擇，對補(bǔ)償效果影響較大，并且處理后的數(shù)據(jù)的振幅相對關(guān)系難以描述。2001年，G F Margrave等［9］提出了一種Gabor反褶積方法，將地層視為單Q值的均勻黏彈性介質(zhì)，且僅考慮吸收衰減條件下直接將Wiener反褶積算法擴(kuò)展到待分析信號為時變的情況。這種方法不需要預(yù)先估計(jì)Q值就可完成反Q濾波，但對平滑窗依賴性較強(qiáng)，而且經(jīng)這種方法處理后的地震記錄不能很好地刻畫反射系數(shù)的局部能量相對關(guān)系（直觀地講會產(chǎn)生等幅效應(yīng)，很像經(jīng)自動增益控制（AGC）處理后的結(jié)果），在某些情況下不利于地震資料解釋。2002年，G F Margrave等［10］對上述方法進(jìn)行了改進(jìn)，用雙曲型平滑代替矩形窗平滑，壓制了AGC效應(yīng)，并于2005年系統(tǒng)地闡述了Gabor反褶積方法的原理［11］。由于Gabor反褶積方法采用均勻Gabor標(biāo)架分析信號，分析窗的長度難以適應(yīng)待分析信號各段的情況，并且按照單Q值模型同時求解不同時刻的反褶積算子，使得該方法不能用于地層Q值隨深度變化的情況。2009年，高靜懷 等［12］提出了基于地震記錄變子波模型提高地震資料分辨率的方法，但該方法是在單時窗內(nèi)采用維納反褶積方法進(jìn)行拓頻，沒有考慮不同深度時窗由吸收引起的振幅衰減影響。

筆者在高靜懷 等［12］方法的基礎(chǔ)上，考慮地層吸收的振幅衰減影響，基于反射地震記錄變子波模型與Gabor變換，在地震數(shù)據(jù)地層結(jié)構(gòu)的約束下構(gòu)造一組自適應(yīng)于地震記錄的分析時窗，對地震記錄進(jìn)行自適應(yīng)非均勻劃分，在時頻域利用非線性壓縮映射提取時變子波，對非平穩(wěn)地震記錄做頻譜拓寬和振幅校正，估計(jì)補(bǔ)償因子，對不同時窗之間由吸收引起的能量衰減進(jìn)行補(bǔ)償，反變換回時間域得到相對保幅的高分辨率處理結(jié)果。本文方法在渤海某油田進(jìn)行了實(shí)際應(yīng)用，顯著提高了該油田地震資料的分辨率，并在波阻抗反演、層序檢測等方面發(fā)揮了重要作用。

1 基于變時窗Gabor變換的高分辨率處理技術(shù)

本文提出的基于變時窗Gabor變換的高分辨率處理技術(shù)包括自適應(yīng)時窗分解、變時窗Gabor正變換、時變子波提取、高分辨率處理和變時窗Gabor反變換等5個步驟，其中最關(guān)鍵的3個步驟是自適應(yīng)時窗分解、時變子波提取和高分辨率處理。

1.1 自適應(yīng)時窗分解

1）Gabor變換。

如果一個函數(shù)集｛Ψj∶j∈Z｝的和為1，即

則此函數(shù)集構(gòu)成對單位1的分割［11］。令φ（t）滿足

則集合｛φm（t）∶m∈Z｝構(gòu)成一個單位1的分割，取分析時窗函數(shù)g（t）=φ（t）u，綜合時窗函數(shù)γ（t）=φ（t）1-u，其中0≤u≤1，可以證明由｛gm∶m∈Z｝和｛γm∶m∈Z｝能夠構(gòu)造一對Gabor標(biāo)架。取u=1，則g（t）=φ（t），γ（t）=1，對于待分析信號s（t），相應(yīng)于上述Gabor標(biāo)架的Gabor變換為

Gabor逆變換為

可以看出，這種Gabor變換方法是通過設(shè)置u=1，使得反變換所用的綜合時窗為單位1，不需要額外計(jì)算綜合時窗，并且式（3）、（4）可以通過快速Fourier變換實(shí)現(xiàn)，計(jì)算效率較高。為了表述直觀，這里用連續(xù)變量表示時間和頻率，用離散變量表示各時窗的中心在時間軸上的位置。對于滿足一定條件的和不為1的函數(shù)集，也可以通過歸一化的方法將其變?yōu)閷挝?的分割。

2）構(gòu)造分析時窗。

由于Gabor變換在每個時間采樣點(diǎn)都有一個分析時窗，冗余度很高，因此如果在Gabor域逐窗處理數(shù)據(jù)，則計(jì)算量會很大。而且如果選擇的時窗太窄，則頻率分辨率低，窗內(nèi)的時變子波無法提??；如果為了提高頻率分辨率，將窗口變寬，則窗內(nèi)平穩(wěn)假設(shè)的近似程度可能會變差。針對上述問題，提出了構(gòu)造自適應(yīng)于非平穩(wěn)地震記錄的分析時窗的方法，具體步驟為：

①提取包絡(luò)峰值。地震道的包絡(luò)峰值與反射界面有一定的對應(yīng)關(guān)系，在一定程度上可以大致反映地層的層序結(jié)構(gòu)，因此，如果使用地震道的包絡(luò)峰值約束分子窗的構(gòu)造，那么構(gòu)成的分子窗在橫向上將與地層結(jié)構(gòu)有關(guān)，有利于保持處理后資料的橫向連續(xù)性（圖1）。根據(jù)復(fù)道分析原理，設(shè)s＊（t）為s（t）的Hilbert變換，則s（t）的包絡(luò)為

圖1 實(shí)際地震剖面及其包絡(luò)峰值剖面Fig.1 Field seismic section and peak envelope section

②生成滿足單位分解的原子窗族。在每個采樣點(diǎn)建立的分析時窗被稱為原子窗，所有采樣點(diǎn)對應(yīng)的分析時窗的集合被稱為原子窗族。選擇基本原子窗函數(shù)G（t）為Lamoureux函數(shù)，N為地震道采樣點(diǎn)個數(shù)，用G（t）=G（t-jΔt）表示中心位于第j個采樣點(diǎn)上的原子窗，按式（6）對原子窗族｛Gj（t）∶1≤j≤N｝進(jìn)行歸一化［13］處理，由式（1）可得到一組滿足單位分解的原子窗族｛gj（t）∶1≤j≤N｝，即

③構(gòu)造初始分子窗。在包絡(luò)的控制下將地震記錄自適應(yīng)地分成多個時間段，通過對包絡(luò)求兩級差分可以得到這些時間段的分割點(diǎn)，然后將相鄰分割點(diǎn)間的小原子窗疊加起來就得到了與各個分段對應(yīng)的分析時窗（形象的被稱為分子窗）。設(shè)第k個分子窗的起點(diǎn)和終點(diǎn)分別為Mk-1和Mk，則第k個分子窗為

④分子窗能量歸一化。第k個分子窗的能量為

則能量歸一化后第k個分子窗為Ψk（t）／Ek。

圖2給出了自適應(yīng)時窗的分解過程。

1.2 時變子波提取

在很多情況下（如爆炸震源等），地震記錄震源子波的振幅譜是頻率的單峰函數(shù)。對于子波振幅譜非單峰函數(shù)的地震記錄，也可以將其轉(zhuǎn)化為單峰函數(shù)的情況。假定地震記錄中地震子波的振幅譜是頻率的單峰、光滑函數(shù)，可以構(gòu)造一個壓縮算子，把提取子波振幅譜的問題轉(zhuǎn)化為求解該算子的不動點(diǎn)問題，從而根據(jù)不動點(diǎn)迭代算法從地震記錄中提取子波振幅譜。

設(shè)x0∈（a，b），取δ0＞0足夠小，對于任意給定的0＜δ?1，若

則Ωx0，δ0，δ是L2［a，b］空間中的子集，由于區(qū)間［a，b］是有限的，可以推導(dǎo)出L2［a，b］?Lp［a，b］，0＜p＜1，

圖2 自適應(yīng)時窗的分解過程Fig.2 Adaptive time-window decomposition process

顯然，0≤F（?；f）≤1。

設(shè)c＞0，α＞1，β＞1，任 取?∈Ωx0，δ0，δ，定 義Ωx0，δ0，δ上的算子P為

顯然，δ≤P（?；f）≤c＋δ，可以推導(dǎo)出P（?；f）∈L2［a，b］，因此，P是Ωx0，δ0，δ?L2［a，b］到L2［a，b］的非線性算子。由Banach空間算子不動點(diǎn)定理可以證明，P是Ωx0，δ0，δ到Ωx0，δ0，δ的壓縮算子，并且在Ωx0，δ0，δ中存在唯一不動點(diǎn)。

設(shè)?（f）為地震記錄的振幅譜，f為頻率，fcut為振幅譜的截止頻率，0＜p＜1，設(shè)迭代初值為

因此，Ωx0，δ0，δ?Lp［a，b］。

任取?∈Ωx0，δ0，δ，令

其中壓縮算子P的參數(shù)由最小二乘得到。因此，通過建立如下迭代

最終得到不動點(diǎn)Faim，則估計(jì)的子波振幅譜｜?w（f）｜為F1／paim（下文稱為F函數(shù)）。

圖3給出了圖2中合成記錄時變子波振幅譜的估計(jì)效果，其中反射系數(shù)采樣點(diǎn)數(shù)為500，采樣間隔2 ms，子波為50 Hz雷克子波。從圖3可以看出，估計(jì)的子波振幅譜和實(shí)際子波振幅譜吻合很好。

圖4展示了在圖2中各分子窗內(nèi)提取的時變子波的振幅譜，可以看出，受地層吸收的影響，隨著深度增加，子波振幅譜主頻逐漸降低。因此，下一步的高分辨率處理就是要把這些傳播過程中損失的信息補(bǔ)償回來。

1.3 高分辨率處理

各分子窗對應(yīng)的振幅譜既含有等效子波的信息，也含有該段反射系數(shù)序列的信息。對于待分析的非平穩(wěn)地震記錄來說，分子窗都是自適應(yīng)于地震記錄構(gòu)造的，因此可認(rèn)為各分子窗內(nèi)片段的F函數(shù)近似等于該段的等效子波的振幅譜乘以一個常數(shù)。等效子波的振幅譜引起該段地震記錄頻帶變窄，使分辨率降低，因此從該段地震記錄中剔除等效子波的效應(yīng)就可提高其分辨率。

圖3 圖2中合成記錄時變子波振幅譜估計(jì)效果Fig.3 Time-varying wavelet amplitude spectrum estimation of the synthetic signal in Fig.2

圖4 圖2中合成記錄各分子窗內(nèi)提取的時變子波的振幅譜Fig.4 Amplitude spectrum extracted in each molecules window of the synthetic signal in Fig.2

1）不帶衰減補(bǔ)償?shù)耐仡l。

以第k段為例，為了提高該段的分辨率，采用零相位Wiener反褶積算法拓寬該段的頻帶，即

式（14）中：?k（f）為原始數(shù)據(jù)的頻譜；?Bk（f）為數(shù)據(jù)拓頻后的頻譜；｜?wk（f）｜為第k段對應(yīng)的F函數(shù)，它與第k段的等效子波振幅譜＾wk（f）相似；Amax為｜?wk（f）｜的最大值；ε是一個小于1的無量綱常數(shù)。

由于各片段對應(yīng)的F函數(shù)和等效子波的振幅譜之間存在一個比例系數(shù)，這個比例系數(shù)和該段內(nèi)反射系數(shù)的能量有關(guān)。也就是說，｜?wk（f）｜中既包含第k段等效子波的能量，也包含該段內(nèi)反射系數(shù)的能量，因此，式（14）的結(jié)果會產(chǎn)生類似于AGC的等幅效應(yīng)［10］。為了讓得到的高分辨結(jié)果的振幅有意義，采用式（15）校正?Bk（f）的能量，即

2）帶衰減補(bǔ)償?shù)耐仡l。

對于第k個分子窗截出的地震記錄片段，把從參考子波（震源子波）到第k個分子窗之間的介質(zhì)視為均勻黏彈性介質(zhì)，介質(zhì)的等效品質(zhì)因子記為Qk。令參考子波從震源傳播到中心為k的分子窗處所用的時間為Tk，則平面波在頻率域滿足因果律的傳播算子可表示為［14-16］

式（16）中：H為在某個時刻t對頻率f的Hilbert變換。

用（f）表示參考子波的Fourier頻譜，用（f）表示第k段上等效子波的Fourier頻譜，根據(jù)波的傳播理論有

將式（16）代入式（17）中，且僅考慮振幅部分，可得

由于本文構(gòu)造的分子窗所覆蓋的信號段近似平穩(wěn)，所以每個窗截取的信號段有一個近似不變的等效子波，該段的F函數(shù)可表示為該等效子波的振幅譜乘以一個常數(shù)。記（f）的F函數(shù)為｜（f）｜，Ck為第k段的常數(shù)比例因子，則

將式（18）代入式（19）中，可得

對式（20）兩邊取對數(shù)，采用對數(shù)譜比值法可計(jì)算出Qk，從而得到｜h（f，Tk）｜。令

用βk乘以?k（f），得到補(bǔ)償后的片段為

記（f）的F函數(shù)為｜（f）｜，將式（14）和式（15）分別作用于（f），可得拓頻結(jié)果為

第k段對應(yīng)的帶衰減補(bǔ)償?shù)耐仡l結(jié)果為

圖5 圖2中合成記錄高分辨率處理效果Fig.5 High resolution processing result of the synthetic signal in Fig.2

對比圖5可以看出，反射系數(shù)①②③④（圖5a）在合成地震記錄（圖5b）上不清晰。反射系數(shù)①②在譜白化處理結(jié)果（圖5c）中比在合成地震記錄上清晰，而反射系數(shù)③④在譜白化處理結(jié)果中仍不清晰。在圖5d、e中，反射系數(shù)①②③④都得到了很好的刻畫，圖5d中的高分辨地震記錄保持了圖5b原始地震記錄的相對能量關(guān)系；圖5e中的高分辨地震記錄恢復(fù)了衰減的能量，與圖5a中的反射系數(shù)序列有較好的對應(yīng)關(guān)系。這表明，該技術(shù)不僅能夠提高地震資料分辨率，而且能夠補(bǔ)償由地層吸收引起的能量衰減。

2 應(yīng)用效果

圖6、7是渤海A油田提高分辨率前后地震剖面及頻譜對比圖，可以看出，提高分辨率后剖面波組特征保持得較好，頻譜的主頻和帶寬都得到了提升。

圖8是渤海A油田提高分辨率處理前后波阻抗反演結(jié)果對比圖，可以看出，提高分辨率后反演結(jié)果和井的吻合度更高，波阻抗資料的分辨率也有了較大提高，反演精度更高。

圖6 渤海A油田提高分辨率前后地震剖面對比Fig.6 Comparison between field data section and high resolution section in A oilfield，Bohai sea

圖7 渤海A油田提高分辨率前后地震資料頻譜對比Fig.7 Comparison of amplitude spectrum between field data and high resolution data in A oilfield，Bohai sea

圖8 渤海A油田提高分辨率前后波阻抗反演結(jié)果對比Fig.8 Comparison of impedance inversion result between field data and high resolution data in A oilfield，Bohai sea

圖9是渤海A油田提高分辨率前后層序檢測（瞬時頻率屬性）對比圖，可以看出，提高分辨率后紅圈范圍內(nèi)的幾套三角洲砂體的疊置關(guān)系更加明顯，砂體邊界也更加清晰。

圖9 渤海A油田提高分辨率前后層序檢測（瞬時頻率屬性）對比圖Fig.9 Comparison of sequence detection（instantaneous frequency）between field data and high resolution data in A oilfield，Bohai sea

3 結(jié)束語

本文提出的基于變時窗Gabor變換的高分辨處理技術(shù)克服了時變譜白化等方法中窗長難以確定和等幅效應(yīng)的問題，以及反Q濾波等方法中Q值難以求準(zhǔn)的問題，理論分析表明其適用條件更接近地下實(shí)際情況，是一種相對保幅的提高地震資料分辨率的新方法。實(shí)際應(yīng)用表明，該方法顯著提高了地震資料的分辨率，能夠在波阻抗反演、層序檢測等領(lǐng)域發(fā)揮重要作用。

［1］YILMAZ O.Seismic data processing［M］.Tulsa：Society of Exploration Geophysicists，1987：49-61，85-94.

［2］孫雷鳴，萬歡，陳輝，等.基于廣義S變換地震高分辨率處理方法的改進(jìn)及在流花11-1油田的應(yīng)用［J］.中國海上油氣，2011，23（4）：234-237.Sun Leiming，Wan Huan，Chen Hui，et al.An improved method of seismic high-resolution processing based on generalized S transform and its application in LH11-1 oilfield［J］.China Offshore Oil and Gas，2011，23（4）：234-237.

［3］陳寶書，汪小將，李松康，等.海上地震數(shù)據(jù)高分辨率相對保幅處理關(guān)鍵技術(shù)研究與應(yīng)用［J］.中國海上油氣，2008，20（3）：162-166.Chen Baoshu，Wang Xiaojiang，Li Songkang，et al.Key technique research and application in relative amplitude preservation processing with high resolution for offshore seismic data［J］.China Offshore Oil and Gas，2008，20（3）：162-166.

［4］RICKER N.The form and nature of seismic wavelets and the structure of seismograms［J］.Geophysics，1940，5（4）：348-366.

［5］DOBRIN MB，SAVIT C H.Introduction to geophysical prospecting［M］.New York：McGraw Hill Book Co.，1988：200-220.

［6］ZIOLKOWSKI A.Why don′t we measure seismic signatures？［J］.Geophysics，1991，56（2）：190-201.

［7］HALE I D.Q adaptive deconvolution［R］.Palo Alto：Stanford Exploration Project Report Number 30，1982：133-158.

［8］GIBSON B，LARNER K L.Comparison of spectral flattening techniques［R］.Western Geophysical Company，1982：10-30.

［9］MARGRAVE G F，LAMOUREUX MP.Gabor deconvolution［R］.Calgary：CREWES Research Report，2001：241-276.

［10］MARGRAVE G F，HENLEY D C，LAMOUREUX MP，et al.An update on Gabor deconvolution［R］.Calgary：CREWES Research Report，2002：11-27.

［11］MARGRAVE G F，GIBSON P C，GROSSMAN J P，et al.The Gabor transform，pseudo differential operators，and seismic deconvolution［J］.Integrated Computer-Aided Engineering，2005，12（3）：43-45.

［12］高靜懷，汪玲玲，趙偉.基于反射地震記錄變子波模型提高地震記錄分辨率［J］.地球物理學(xué)報(bào)，2009，52（5）：1289-1300.Gao Jinghuai，Wang Lingling，Zhao Wei.Enhancing resolution of seismic traces based on the changing wavelet model of the seismogram［J］.Chinese Journal of Geophysics，2009，52（5）：1289-1300.

［13］GROSSMAN J P，MARGRAVE G F，LAMOUREUX MP.Constructing adaptive，nonuniform Gabor frames from partitions of unity［R］.Calgary：CREWES Research Report，2002：1-10.

［14］WANG Yanghua.A stable and efficient approach of inverse Q filtering［J］.Geophysics，2002，67（2）：657-664.

［15］WANG Yanghua.Inverse Q-filter for seismic resolution enhancement［J］.Geophysics，2006，71（3）：51-61.

［16］余連勇，胡光義，趙巖，等.穩(wěn)定的反Q濾波統(tǒng)一算法及其在地震資料高分辨率處理中的應(yīng)用［J］.中國海上油氣，2014，26（4）：29-33.Yu Lianyong，Hu Guangyi，Zhao Yan，et al.A unified algorithm of stable inverse Q filtering and its application to highresolution processing of seismic data［J］.China Offshore Oil and Gas，2014，26（4）：29-33.