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

(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]等,用于地震資料解釋。在傾角估計(jì)方面,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ù)振幅橫向變化對(duì)傾角估算精度的影響,WANG等[20]深入分析了復(fù)偏導(dǎo)數(shù)與瞬時(shí)振幅及瞬時(shí)相位的偏導(dǎo)數(shù)間的關(guān)系,提出了基于瞬時(shí)相位的梯度結(jié)構(gòu)張量法,并基于多窗口分析技術(shù)提出了基于瞬時(shí)相位的多窗口梯度結(jié)構(gòu)張量法[21]。利用時(shí)頻分析的多尺度分辨能力對(duì)實(shí)際地震資料進(jìn)行高精度時(shí)頻分析獲得高精度的時(shí)頻瞬時(shí)相位譜,抽取不同頻率的多尺度高精度時(shí)頻瞬時(shí)相位譜估算地層傾角和曲率屬性可以提高該方法的應(yīng)用效果。

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

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

1.1 基本原理

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

(1)

(2)

方程(1)可寫為:

(3)

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

(4)

且

(5)

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

(6)

因此,有

(7)

1.2 合成記錄分析

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

(8)

且

(9)

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

(10)

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

(11)

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

(12)

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

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

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

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

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

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

2.1 基于瞬時(shí)相位的曲率屬性

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

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

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

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

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

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

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

(13)

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

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

圖4 合成信號(hào)瞬時(shí)屬性分析a 無噪信號(hào); b 無噪信號(hào)的瞬時(shí)振幅; c 無噪信號(hào)的瞬時(shí)相位; d 含噪信號(hào); e 含噪信號(hào)瞬時(shí)振幅; f 含噪信號(hào)瞬時(shí)相位

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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