基于高精度時頻瞬時相位譜的多尺度曲率及其應用

李 勇,張固瀾,何承杰,王曉凱,王婷一,殷 瑛,張建軍

(1.西南石油大學地球科學與技術學院,四川成都610500;2.西安交通大學電子與信息工程學院,陜西西安710049;3.中國石油集團東方地球物理公司新興物探開發(fā)處,河北涿州072751)

幾何屬性是最重要的地震屬性[1-2]之一,在地震資料解釋中發(fā)揮著重要作用。基于地震數(shù)據(jù)的層位信息,可獲得地層傾角及方位角[3-4]等幾何屬性,估算曲率[5-9]和相干屬性[10-13]等,用于地震資料解釋。在傾角估計方面,PICOU等[14]提出了非歸一化互相關掃描法;BARNES[15-16]提出了復數(shù)分析法;MARFURT等[17]發(fā)展了相干掃描方法;BAKKER等[18]提出了梯度結構張量法;LUO等[19]提出了加權結構法。為減少地震數(shù)據(jù)振幅橫向變化對傾角估算精度的影響,WANG等[20]深入分析了復偏導數(shù)與瞬時振幅及瞬時相位的偏導數(shù)間的關系,提出了基于瞬時相位的梯度結構張量法,并基于多窗口分析技術提出了基于瞬時相位的多窗口梯度結構張量法[21]。利用時頻分析的多尺度分辨能力對實際地震資料進行高精度時頻分析獲得高精度的時頻瞬時相位譜,抽取不同頻率的多尺度高精度時頻瞬時相位譜估算地層傾角和曲率屬性可以提高該方法的應用效果。

地層傾角及曲率屬性的估算精度依賴于時頻瞬時相位譜的精度,而時頻瞬時相位譜的精度受時頻分析方法及隨機噪聲的影響。短時傅里葉變換[22]、S變換[23-24]和廣義S變換[25-27]可獲得信號的瞬時振幅及初始相位隨時間和頻率的分布,即時頻瞬時振幅譜和時頻初始相位譜;小波變換[28-33]可獲得信號的瞬時振幅及瞬時相位[34-38]隨時間和尺度(或頻率)的分布,即時頻瞬時振幅譜和時頻瞬時相位譜。改進的短時傅里葉變換[39]通過對待分析的實信號解析信號加上對稱時窗函數(shù)并進行傅里葉變換,獲得待分析實信號的瞬時振幅、初始相位和瞬時相位隨時間和頻率的分布(分別稱為高精度時頻瞬時振幅譜、時頻初始相位譜和時頻瞬時相位譜)。本文利用改進的短時傅里葉變換獲得高精度時頻瞬時相位譜,然后抽取不同頻率的高精度時頻瞬時相位譜估算多尺度曲率屬性,并將該屬性應用于實際三維地震數(shù)據(jù)的不連續(xù)性檢測。

1 改進的短時傅里葉變換

1.1 基本原理

基于短時傅里葉變換、廣義S變換和小波變換,實信號x(t)改進的短時傅里葉變換[39]可表示為:

(1)

(2)

方程(1)可寫為:

(3)

式中,H0(τ,f)反映復信號z(t)p(t-τ)中頻率成分f的瞬時振幅和初始相位,H1(τ,f)反映復信號z(t)p(t-τ)中頻率成分f在t=τ時刻的瞬時振幅和瞬時相位。若定義

(4)

且

(5)

則稱A0(τ,f),θ0(τ,f),A0(τ,fτ),θ0(τ,fτ)和fτ分別為實信號x(t)在t=τ時刻的高精度時頻瞬時振幅譜、時頻初始相位譜、時頻峰值瞬時振幅譜、時頻峰值初始相位譜和時頻峰值頻率;A1(τ,f),θ1(τ,f),A1(τ,fτ)和θ1(τ,fτ)分別為實信號x(t)在t=τ時刻的高精度時頻瞬時振幅譜、時頻瞬時相位譜、時頻峰值瞬時振幅譜和時頻峰值瞬時相位譜。由于

(6)

因此,有

(7)

1.2 合成記錄分析

對改進的短時傅里葉變換、短時傅里葉變換及小波變換進行了時頻分析對比,如圖1、圖2所示。這里假設短時傅里葉變換可表示為:

(8)

且

(9)

式中,f為頻率,λ>0且β≥0。假設小波變換可表示為:

(10)

式中,a為尺度算子,q(t)為小波函數(shù),*表示取共軛,且

(11)

方程(10)的小波變換可表示為:

(12)

圖1a為合成記錄,圖1b和圖1c分別為改進的短時傅里葉變換得到的高精度時頻瞬時振幅譜A1(τ,f)和A0(τ,f),圖1d和圖1e分別為小波變換和短時傅里葉變換得到的時頻瞬時振幅譜。對比圖1b～圖1e 可知,各種時頻瞬時振幅譜的變化趨勢完全相同,改進的短時傅里葉變換得到的高精度時頻瞬時振幅譜變化范圍為0～3,與合成記錄的瞬時振幅變化區(qū)間相同,而小波變換和短時傅里葉變換得到的時頻瞬時振幅譜在0～1.5之間變化。圖1f和圖1g 分別為改進的短時傅里葉變換得到的高精度時頻瞬時振幅譜與小波變換得到的時頻瞬時振幅譜的差值與比值。由圖1f可見,兩者的差值全部大于0;由圖1g可見,當頻率從4Hz變化到26Hz時,比值系數(shù)為2.0,較大比值主要分布在低頻和高頻端。

圖2a和圖2d分別為改進的短時傅里葉變換得到的高精度時頻瞬時相位譜和高精度時頻初始相位譜,兩者變化趨勢相同;圖2b和圖2e分別為小波變換和短時傅里葉變換得到的時頻瞬時相位譜和時頻初始相位譜,兩者變化趨勢也相同。圖2c為改進的短時傅里葉變換的高精度時頻瞬時相位譜與小波變換的時頻瞬時相位譜之差值,圖2f為改進的短時傅里葉變換的高精度時頻初始相位譜與短時傅里葉變換的時頻初始相位譜之差值,差值的變化區(qū)間為[-45°,45°]。當頻率從4Hz變化到26Hz時,差值幾乎為零,較大的差值主要分布在低頻和高頻端。

綜合分析可知,在相同時窗函數(shù)情況下,改進的短時傅里葉變換具有更好的時頻精度。

圖1 合成記錄時頻分析a 合成記錄; b 改進的短時傅里葉變換得到的高精度時頻瞬時振幅譜(n=1); c 改進的短時傅里葉變換得到的高精度時頻瞬時振幅譜(n=0); d 小波變換得到的時頻瞬時振幅譜; e 短時傅里葉變換得到的時頻瞬時振幅譜; f 圖1b與圖1d的差值; g 圖1b與圖1d的比值

圖2 合成記錄時頻分析a 改進的短時傅里葉變換得到的高精度時頻瞬時相位譜; b 小波變換得到的時頻瞬時相位譜; c 圖2a與圖2b的差值; d 改進的短時傅里葉變換得到的高精度時頻初始相位譜; e 短時傅里葉變換得到的時頻初始相位譜; f 圖2d與圖2e的差值

2 基于高精度時頻瞬時相位譜的多尺度曲率屬性

2.1 基于瞬時相位的曲率屬性

曲率屬性估算的核心是利用地震數(shù)據(jù)準確估算地層傾角,而利用地震信號的振幅信息估算地層傾角是常用的方法。由于振幅不僅反映地層傾角信息,還反映波阻抗橫向變化,因此基于振幅估算的曲率屬性受地層傾角及波阻抗橫向變化的綜合影響,在波阻抗橫向變化劇烈的地區(qū)應用效果不佳。瞬時相位僅反映地層傾角信息,與波阻抗橫向變化無關,因此可基于地震信號的瞬時相位估算地層傾角,最終獲得較高精度的曲率屬性。

以下通過二維合成地震記錄及相關屬性進行說明。圖3a為無噪合成地震記錄,共三套層位:第一套(0.05s處)和第三套(0.15s處)均為水平層;第二套層位(0.10s附近)共有7個起伏(中心分別在30,90,150,210,270,330和390道處),且起伏程度隨道號增大而增加,另外,當?shù)捞柸≈禐?40～470時,第二套層位僅有波阻抗變化而無傾角變化。對圖3a合成地震記錄數(shù)據(jù)進行希爾伯特變換,得到的瞬時相位如圖3b所示:當?shù)捞枮?40～470時,第二套層位對應的地震信號瞬時相位無任何變化。圖3c為利用圖3a 所示的地震記錄數(shù)據(jù)計算得到的地層傾角:當?shù)捞枮?40～470時,第二套層位的傾角信息受地震信號振幅的影響嚴重,與理論模型不匹配。圖3d為利用圖3b所示的瞬時相位計算得到的地層傾角,與理論模型的傾角信息完全匹配。

綜上所述,基于瞬時相位估算的地層傾角信息可以避免地震信號振幅的空間變化帶來的假象。基于時頻分析中低頻時頻譜具有較高頻率分辨率、高頻時頻譜具有較高時間分辨率的多尺度思想,可利用改進的短時傅里葉變換方法對地震資料進行時頻譜分析,獲得高精度的時頻瞬時相位譜,并抽取高精度的時頻瞬時相位譜進行曲率屬性計算,以提高曲率屬性的估算精度。

圖3 合成地震記錄及相關屬性分析a 合成地震記錄; b 瞬時相位; c 基于振幅的地層傾角信息; d 基于瞬時相位的地層傾角信息

2.2 基于高精度時頻瞬時相位譜的多尺度曲率屬性

為了減小隨機噪聲對時頻瞬時相位譜的影響,首先利用改進的短時傅里葉變換對地震信號進行時頻分析,并利用高精度的時頻峰值瞬時振幅譜和高精度的時頻峰值瞬時相位譜對輸入的地震信號進行重構;然后利用改進的短時傅里葉變換對重構信號進行時頻分析得到高精度的時頻瞬時相位譜;最后抽取共頻率的高精度時頻瞬時相位譜進行多尺度曲率屬性計算。

基于方程(4)和方程(5),輸入實信號x(t)的重構信號可表示為:

(13)

圖4驗證了高精度時頻峰值瞬時振幅譜和高精度時頻峰值瞬時相位譜重構地震信號的可行性。圖4a 為20Hz余弦信號;圖4b和圖4c分別為該余弦信號的瞬時振幅和瞬時相位;圖4d為圖4a所示余弦信號加入強隨機噪聲的結果;圖4e和圖4f分別為圖4d 所示20Hz含噪余弦信號的瞬時振幅和瞬時相位,其均被隨機噪聲污染。

利用改進的短時傅里葉變換對圖4d所示含噪信號進行時頻分析,得到高精度時頻瞬時振幅譜和高精度時頻瞬時相位譜,并基于高精度時頻峰值瞬時振幅譜及高精度時頻峰值瞬時相位譜進行信號重構。圖5a 為圖4a所示無噪信號的高精度時頻瞬時振幅譜;圖5b為圖4a所示無噪信號的高精度時頻瞬時相位譜;圖5c為圖4d所示含噪信號的高精度時頻瞬時振幅譜,其受隨機噪聲影響較小;圖5d為圖4d 所示含噪信號的高精度時頻瞬時相位譜,其受隨機噪聲干擾嚴重。圖6a為圖5c中各時刻的高精度時頻峰值瞬時振幅譜,圖6b為圖5d中各時刻的高精度時頻峰值瞬時相位譜;圖6c為利用圖6a和圖6b重構的信號,其與圖4a所示的20Hz余弦信號完全一致,說明可基于高精度時頻峰值瞬時振幅譜和高精度時頻峰值瞬時相位譜重構信號。

圖4 合成信號瞬時屬性分析a 無噪信號; b 無噪信號的瞬時振幅; c 無噪信號的瞬時相位; d 含噪信號; e 含噪信號瞬時振幅; f 含噪信號瞬時相位

圖5 無噪信號的高精度時頻瞬時振幅譜(a)和相位譜(b)與含噪信號的高精度時頻瞬時振幅譜(c)和相位譜(d)

圖6 含噪信號的高精度時頻峰值瞬時振幅譜(a)和瞬時相位譜(b)及其重構信號(c)

3 應用實例

利用某工區(qū)三維地震資料檢驗了基于高精度時頻瞬時相位譜的多尺度曲率屬性估算方法的應用效果。處理流程如下:首先利用方程(9)所示的時窗函數(shù)對原始地震數(shù)據(jù)進行改進的短時傅里葉變換,得到高精度時頻瞬時振幅譜和高精度時頻瞬時相位譜;然后利用高精度時頻峰值瞬時振幅譜及高精度時頻峰值瞬時相位譜進行地震數(shù)據(jù)重構;再利用方程(9)所示的時窗函數(shù)對重構地震數(shù)據(jù)進行改進的短時傅里葉變換,得到重構數(shù)據(jù)的高精度時頻瞬時振幅譜和高精度時頻瞬時相位譜,并抽取共頻率的高精度時頻瞬時振幅譜和高精度時頻瞬時相位譜,分別估算最負曲率;最后沿地震解釋層位做切片對比(圖7～圖10)。圖7a所示原始數(shù)據(jù)幾乎不含隨機噪聲,主要檢驗是否可以利用高精度時頻峰值瞬時振幅譜及高精度時頻峰值瞬時相位譜重構原始地震數(shù)據(jù)。圖7b為重構數(shù)據(jù)振幅切片,對比可見,除藍色箭頭標注處外,重構數(shù)據(jù)和原始數(shù)據(jù)幾乎沒有差異。

圖8a和圖8b分別為重構數(shù)據(jù)30Hz和50Hz的高精度時頻瞬時振幅譜切片;圖8c和圖8d分別為基于重構數(shù)據(jù)30Hz和50Hz的高精度時頻瞬時振幅譜估算的最負曲率切片。對比圖8c和圖8d可見:圖8c主要反映大尺度地質(zhì)體結構,圖8d主要反映小尺度地質(zhì)體結構,兩者互補。對比圖8a和圖8c可見:

圖7 原始數(shù)據(jù)振幅切片(a)與重構數(shù)據(jù)振幅切片(b)

高精度時頻瞬時振幅譜相對穩(wěn)定區(qū)域的曲率精度較高,高精度時頻瞬時振幅譜變化劇烈區(qū)域的曲率精度較低。對比圖8b和圖8d得到的結論相同。

圖9a和圖9b分別為重構數(shù)據(jù)30Hz和50Hz的高精度時頻瞬時相位譜切片,圖10a和圖10b分別為基于重構數(shù)據(jù)30Hz和50Hz的高精度時頻瞬時相位譜估算的最負曲率切片。對比圖10a和圖10b 可見:相對于基于低頻成分的高精度時頻瞬時相位譜估算的曲率屬性而言,基于高頻成分的高精度時頻瞬時相位譜估算的曲率屬性具有更高的分辨率。圖10a 比圖8c具有更高的分辨率,圖10b比圖8d具有更高的分辨率,說明基于高精度時頻瞬時振幅譜估算的曲率屬性容易受地震數(shù)據(jù)振幅空間變化的影響,不能完全反映地層的真實結構;基于高精度時頻瞬時相位譜估算的曲率屬性較基于高精度時頻瞬時振幅譜估算的曲率屬性具有更高的精度。

圖8 沿層切片a 重構數(shù)據(jù)30Hz的高精度時頻瞬時振幅譜; b 重構數(shù)據(jù)50Hz的高精度時頻瞬時振幅譜; c 重構數(shù)據(jù)30Hz高精度時頻瞬時振幅譜估算的最負曲率; d 重構數(shù)據(jù)50Hz高精度時頻瞬時振幅譜估算的最負曲率

圖9 沿層切片a 重構數(shù)據(jù)30Hz的高精度時頻瞬時相位譜; b 重構數(shù)據(jù)50Hz的高精度時頻瞬時相位譜

圖10 沿層切片a 重構數(shù)據(jù)30Hz高精度時頻瞬時相位譜估算的最負曲率; b 重構數(shù)據(jù)50Hz高精度時頻瞬時相位譜估算的最負曲率

4 結論

1) 改進的短時傅里葉變換可獲得信號的高精度時頻瞬時振幅譜、時頻瞬時相位譜和時頻初始相位譜。

2) 瞬時相位反映了地下構造的變化,基于瞬時相位可精確估算地層傾角,提高曲率屬性的估算精度。但瞬時相位對噪聲極其敏感,需充分壓制隨機噪聲并提高瞬時相位的精度。

3) 基于改進的短時傅里葉變換獲得的高精度時頻瞬時相位譜估算的多尺度曲率屬性可提高宏觀及微觀地下構造的識別精度。為降低隨機噪聲對時頻瞬時相位譜的影響,進一步提高基于高精度時頻瞬時相位譜的多尺度曲率屬性的估算精度,可基于改進的短時傅里葉變換獲得的高精度時頻峰值瞬時振幅譜和高精度時頻峰值瞬時相位譜對原始信號進行重構,再利用重構信號的高精度時頻瞬時相位譜進行多尺度的曲率分析。