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

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

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

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

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

1 改進(jìn)的短時傅里葉變換

1.1 基本原理

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

(1)

(2)

方程(1)可寫為:

(3)

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

(4)

且

(5)

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

(6)

因此,有

(7)

1.2 合成記錄分析

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

(8)

且

(9)

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

(10)

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

(11)

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

(12)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(13)

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

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

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

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

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

3 應(yīng)用實(shí)例

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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