楊培杰 羅紅梅 王金鐸
(中國石化勝利油田分公司勘探開發研究院,山東東營 257015)
利用地震數據時頻分析[1-3]進行儲層預測的應用非常廣泛,這些方法通過分析地震數據頻譜的振幅譜,利用含流體儲層高頻衰減、低頻陰影的特征輔助進行流體檢測以及薄互層識別,并取得了較好的應用效果。相位譜也是地震信號頻譜分析中非常重要的組成部分。目前,研究人員主要用相位信息輔助進行地層結構的解釋,如通過瞬時相位信息識別超剝線、計算薄層厚度等[4-5],但對于相位譜的進一步應用顯得不足,特別是針對地震數據中相位信息分解與重構的研究幾乎是空白,主要原因是相位譜不太直觀,并且對相位譜的有效應用有一定的難度。因此,如何從地震數據中有效地提取相位信息并進行分析是地震數據頻率域處理的難點,也是今后研究的熱點。
以高分辨率譜估計方法[6]為基礎,選擇合適的窗函數,逐點提取地震道的時頻信息,針對相位譜信息,在[-180°,180°]區間設定待重構的相位值及相位容許誤差,并將設定值以外的相位信息置為零,再進行傅里葉反變換,即可得到分相位重構后的地震道。該地震道只包含設定相位的信息,實現了地震數據的分相位重構。該方法能用于地震數據的處理和解釋,可以根據需要在不同的相位內分析地震數據,實現薄層或特殊反射體的有效識別,具有廣泛的應用前景。
譜估計的目的是得到地震數據的振幅譜和相位譜,這些振幅譜和相位譜包含了豐富的信息,進而可以進行儲層預測。目前常用的譜估計方法包括短時傅里葉變換[7]、小波變換[8]、S變換[9]、匹配追蹤分解[10]等,這些方法在石油勘探中得到了廣泛的應用。
本文采用基于反演理論[11-13]的譜估計方法(Inversion based Spectral Analysis,ISA)實現地震信號的譜估計,與其他的譜估計方法相比,ISA方法具有更高的分辨率和穩定性。
考慮下面的正演方程
Fm=d
(1)
式中:m表示待求的頻率系數向量;d表示加窗的地震數據;F表示不同頻率的三角函數構成的矩陣
F(t,f)=cos(2πft)+i sin(2πft)
(2)
式中:t表示時間;f表示頻率; i為虛數單位。
為了提高時頻估計結果的穩定性和分辨率,可以通過Hilbert[14-15]變換將d轉換成復數的形式
d=dr+idi
(3)
式中:dr表示d的實部;di表示d的虛部。
從式(1)中可以看出,已知觀測數據d和正演矩陣F,求m,這是一個典型的反問題,實際應用中,式(1)往往是超定的,一般可以通過最小二乘法求解
me=(FTF)-1FTd
(4)
式中:me為最小二乘反演結果;FT表示F的共軛轉置矩陣。
與其他地震反演過程[16-19]類似,為了提高式(4)求解的穩定性,需要加入約束信息[20],進一步得到
me=(FTF+αI)-1FTd
(5)
式中:α為常數;I表示單位對角矩陣。此時的me即為ISA結果。在加入αI約束分量后,對于信噪比較低、時間段較短的地震數據的譜估計會有更好的效果,可以提高求解過程的穩定性和分辨率。進一步將ISA的分析結果表示為
Fd(f,t)=ISA(d)
(6)
ISA只是一種譜估計的方法,為了得到時頻譜,需要逐點對地震數據進行ISA分析。以待分析點為中心,對待分析信號加窗函數,目的是提取待分析點前后的若干點,窗函數可以有效減小頻譜能量泄漏,窗函數可以加在時域上,也可以加在頻域上,但在時域上加窗更為普遍。實際應用的窗函數主要包括:①冪窗(如矩形、三角形、梯形等);②三角函數窗(如漢寧窗、海明窗等);③指數窗(如高斯窗等)。文中使用的是漢寧窗,漢寧窗適用于非周期性的連續信號,地震數據實際上就是一種非周期性的連續信號。一般設定窗函數的時間長度為主頻的倒數,如主頻為25Hz,則窗長為40ms,具體可根據分相位重構的效果進行判斷。
在通過ISA方法得到地震數據的時頻譜Fd(f,t)后,時變振幅譜和相位譜的計算公式如下
(7)
式中:A(f,t)表示地震數據的振幅譜;P(θ,t)表示地震數據的相位譜;θ表示相位; imag(·)表示求虛部; real(·)表示求實部。
進一步將振幅譜和相位譜統一表示為
B(f,θ,t)={A(f,t),P(θ,t)}
(8)
圖1為某一地震數據及其ISA時頻分析結果。從圖中可以看出,ISA時頻分析結果具有較高的頻率分辨率,并且相位譜在頻率方向也比較平穩,說明該方法具有較高的相位分辨率和穩定性。

圖1 地震信號及其ISA分析結果
分相位數據的重構是方法的核心,設定相位值
Pr∈[-180°,180°]
(9)
式中Pr表示待重構的相位值。
由于實際地震數據中含有各種噪聲,其瞬時相位譜復雜多變,很難得到較穩定的相位信息,往往呈現斜線、折線、震蕩等特點,無法準確地提取所設定的相位值的相位信息。因此,為了保證相位重構后地震數據的穩定性,除了先要設定待重構的相位值外,還需要給定相位容許誤差
|P(θ,t)-Pr|<εb
(10)
式中εb表示相位容許誤差,即當某一時刻的瞬時相位信息與設定的相位值之差的絕對值小于相位容許誤差時,就保留該時刻的頻譜信息;否則就將該時刻時頻譜的實部和虛部都置為零,即只保留設定的待重構相位值的頻率信息,其他相位所對應的頻率信息都去除。
相位容許誤差的大小可以通過試算法得到,給定若干相位容許誤差,分別對數據進行分相位重構,并對重構結果進行分析,選擇合適的相位容許誤差。
用待重構的相位值Pr及相位容許誤差εb對相位信息進行處理,然后在頻譜的一定頻率范圍內進行傅里葉反變換,即可得到分相位重構后的地震數據,公式如下

(11)
式中:D′(θPr,t)表示重構后只包含指定相位信息的地震數據;θPr表示某一指定的相位值;f1、f2表示反變換所用的頻率范圍,一般令f1為0、f2為奈奎斯特頻率(fN)[21],當f1>0、f2 將[-180°,180°]區間所有分相位重構后地震數據集合,就得到了相位道集。相位道集是相位分解與分重構中一個非常重要的概念,為了進一步說明,假設θ只取[-180°,180°]區間的整數,定義下面的數據集合 D′(θ,t)=[D′(θ-180°,t),…,D′(θ0°,t),…, D′(θ180°,t)] (12) 式中:D′(θ,t)為相位道集;D′(θ-180°,t)為-180°相位的地震數據,依此類推。 圖2 方法流程圖 相位道集的縱坐標表示時間,橫坐標表示相位,包含了[-180°,180°]區間的相位信息。相位道集的一個重要性質是將相位道集的數據疊加,可無損恢復原始地震數據,即 (13) 式中:D(t)表示地震數據;D′(t)表示將相位道集疊加后的地震數據。 圖3為一個地震道的分相位重構效果。相位道集也稱為分相位體,與分頻體的概念類似,所有頻率的分頻體疊加,可無損恢復地震數據;同理,分相位體疊加,亦可無損恢復地震數據。圖3c中兩條曲線幾乎重合,說明了本文提出的分相位重構方法的準確性。 圖3 地震數據分相位重構 圖4a為一個合成信號,由5個不同相位的信號組成,從左到右分別為-90°、-40°、0°、40°、90°相位,并加入信噪比為4∶1的噪聲。用ISA方法對其進行時頻分解,得到其時頻振幅譜(圖4b)和時頻相位譜(圖4c),可以看出這5個子信號的振幅譜一樣,主頻都是25Hz左右,帶寬也基本一樣,但是相位譜卻明顯不同。 對該信號進行指定相位的重構,重構-90°、0°、90°相位分量的信號,如圖5所示。-90°相位分量只包含了合成信號中-90°相位分量的成分,波形特征為從波谷到波峰的變化; 0°相位分量只包含了合成信號中0°相位分量的成分,波形特征為近似對稱,一個主瓣外加兩個旁瓣; 90°相位分量只包含了合成信號中90°相位分量的成分,波形特征為從波峰到波谷的變化。通過對比可以看出,這3個相位分量信號與原始合成信號相應分量信號的波形基本一致,時間位置相同,說明了本文提出的分相位重構方法的有效性和準確性。 圖4 合成信號及其ISA分析結果 圖5 合成信號分相位重構 實際應用中地震數據的長度應大于一個波長,為了保證分相位重構過程的穩定性和準確性,要求地震資料的信噪比大于4∶1。在完成一維地震數據的分相位重構后,對地震數據進行道循環,即可完成三維地震數據的分相位重構。對于三維地震數據的分相位重構,相位容許誤差的設定非常重要,相位容許誤差過大,會包含過多的待重構相位以外的信息,相位容許誤差過小,相位重構后地震數據的連續性會很差。一般來說,相位容許誤差的設置原則是,在保證指定相位重構后剖面連續性的情況下,盡量減小相位容許誤差。 圖6為二維地震數據分相位重構結果,該數據摘自文獻[6],通過圖像處理技術轉為SEGY格式。從圖中可見,通過應用本文方法,地震亮點反射被放大聚焦,同時分辨率也得到了提高,并且不同相位分量地震剖面所表現出的地震同相軸以及亮點的形態不一樣。對于零相位的地震資料,當儲層波阻抗大于上下圍巖時,一般對應著90°相位分量地震數據; 當儲層波阻抗小于上下圍巖時,一般對應著-90°相位分量地震數據。實際過程中需結合研究區儲層的速度分析結果對相位重構結果進行綜合應用。 圖6 地震數據分相位重構 如何從地震數據中有效地提取相位信息,并進一步分析處理,是目前地震數據頻率域分析的熱點和難點之一。本文實現了地震數據的分相位重構,進而可以得到不同相位分量的地震數據,這些不同相位分量的數據可以從不同的角度反映地質體的特征,從而更好地凸顯地質異常體。下一步將開展分相位重構方法的實際應用研究。相信相位分析方法在巖性以及地層類油氣藏儲層預測中將發揮重要的作用。

2 一維信號驗證


3 地震數據的分相位重構

4 結束語