楊培杰,隋風(fēng)貴

(1.中國(guó)石油化工股份有限公司勝利油田分公司勘探開(kāi)發(fā)研究院,山東東營(yíng)257015;2.中國(guó)石油化工股份有限公司勝利油田分公司,山東東營(yíng)257001)

時(shí)頻分析[1]在地球物理勘探中的應(yīng)用十分廣泛。目前比較常用的時(shí)頻分析方法主要有短時(shí)傅里葉變換[2]、小波變換[3]、S變換[4]、希爾伯特-黃變換[5]、二次型時(shí)頻分布[6]和匹配追蹤分解[7]等。短時(shí)傅里葉變換相當(dāng)于地震道在一個(gè)固定的時(shí)間窗上與正弦基的互相關(guān),一旦窗函數(shù)選定,時(shí)頻分辨率便可確定下來(lái),這使得短時(shí)傅里葉變換對(duì)非平穩(wěn)地震信號(hào)的分析存在局限性;小波變換是地震道與小波字典的互相關(guān),是對(duì)短時(shí)傅里葉變換的改進(jìn),小波變換受測(cè)不準(zhǔn)原理的影響,需要在時(shí)間分辨率和頻率分辨率之間進(jìn)行權(quán)衡,因此很難滿足精細(xì)儲(chǔ)層預(yù)測(cè)的要求;S變換是短時(shí)傅里葉變換和小波變換的繼承和發(fā)展,S變換具有良好的時(shí)頻特性,在時(shí)頻分辨率上有所改進(jìn);希爾伯特-黃變換將信號(hào)分解為有限的幾個(gè)固有模態(tài)函數(shù),進(jìn)而得到信號(hào)的時(shí)頻譜,其不足是對(duì)噪聲較敏感;二次型時(shí)頻分布在時(shí)頻平面存在非常嚴(yán)重的交叉項(xiàng),導(dǎo)致其在地震勘探中的應(yīng)用并不廣泛;匹配追蹤分解是將地震信號(hào)看作由許多被稱作為原子信號(hào)的加權(quán)總和,將每個(gè)原子信號(hào)的頻譜進(jìn)行疊加就可得到地震道的時(shí)頻譜。

最小二乘法在地震勘探中的應(yīng)用十分廣泛,約束最小二乘地震反演[8]使用最小二乘方法構(gòu)建目標(biāo)函數(shù),極大提髙了解的收斂速度,其優(yōu)點(diǎn)是反演過(guò)程穩(wěn)定、計(jì)算效率高;最小二乘偏移[9-10]是一種基于線性化反演理論的真振幅成像方法,它將成像問(wèn)題當(dāng)作一個(gè)反問(wèn)題來(lái)處理,以尋求更接近于真實(shí)的地下反射系數(shù);羅水亮等[11]提出了一種基于變形Pride模型和Lee模型(P-L模型)的矩陣方程迭代精細(xì)橫波預(yù)測(cè)方法,其核心是通過(guò)最小二乘將橫波估計(jì)目標(biāo)函數(shù)的最優(yōu)化問(wèn)題轉(zhuǎn)化為求解線性矩陣方程組問(wèn)題,從而極大地提高了橫波估計(jì)的效率和穩(wěn)定性。

本文在前人研究成果[12-16]的基礎(chǔ)上,提出了改進(jìn)的約束最小二乘時(shí)頻分析方法(improved constrained least squares time frequency analysis,ICLS)?；诜囱莸乃悸穼㈩l譜分析轉(zhuǎn)換成一個(gè)求解反問(wèn)題的過(guò)程,通過(guò)正則化約束最小二乘方法求解該問(wèn)題,進(jìn)而得到地震信號(hào)的時(shí)頻分析結(jié)果,理論模型和實(shí)際應(yīng)用驗(yàn)證了本文方法的有效性。

1 方法原理

定義如下由待求的頻譜和不同頻率的三角函數(shù)計(jì)算信號(hào)的正演過(guò)程:

d=Sh

(1)

式中:d表示待分析的信號(hào),在實(shí)際應(yīng)用中,d往往是實(shí)數(shù),為了提高時(shí)頻估計(jì)結(jié)果的穩(wěn)定性和分辨率,可以通過(guò)Hilbert變換將d轉(zhuǎn)換成復(fù)數(shù)的形式;S表示不同頻率的三角函數(shù)構(gòu)成的矩陣;h表示待求的頻譜向量。已知信號(hào)d和矩陣S求頻譜h,這是一個(gè)典型的反問(wèn)題,這里用正則化最小二乘的思路進(jìn)行求解。

S的具體形式如下:

S=

(2)

式中:Δt表示時(shí)間采樣間隔;β表示時(shí)間采樣點(diǎn)數(shù);Δf表示頻率采樣間隔;α表示頻率采樣點(diǎn)數(shù)。

對(duì)于精確的數(shù)據(jù),待求的頻譜向量為:

h=S-1d

(3)

但在實(shí)際計(jì)算過(guò)程中,正演方程((1)式)往往是超定的或是欠定的,所以沒(méi)有精確解,即:

d=She+e

(4)

式中:he為計(jì)算的頻譜結(jié)果;e表示實(shí)際數(shù)據(jù)中存在的誤差。

通過(guò)普通最小二乘法即可以得到he的解:

he=(S′S)-1S′d

(5)

式中:S′表示S的共軛轉(zhuǎn)置矩陣。

進(jìn)一步,為了得到信號(hào)的時(shí)頻譜,需要以待分析點(diǎn)為參考,對(duì)待分析信號(hào)加窗函數(shù),目的是提取待分析點(diǎn)處前后的若干點(diǎn),當(dāng)待分析的信號(hào)加上窗函數(shù)后,矩陣S的元素變得相關(guān),因此,為了獲得穩(wěn)定的解,需要加約束,引入對(duì)角矩陣Z:

(6)

式中:Nw=ceil(β/2),ceil表示將實(shí)數(shù)向無(wú)窮遠(yuǎn)的方向舍入到最接近的整數(shù);Diag表示求對(duì)角矩陣。此時(shí),(1)式可以進(jìn)一步寫成:

dw=Swh

(7)

式中:dw=Zd;Sw=ZS。

利用最小二乘法求解(7)式。為了進(jìn)一步提高求解的穩(wěn)定性和唯一性,需要對(duì)求解過(guò)程加約束[17],由于L2范數(shù)正則化可以有解析形式的解,更易于求解,因此選用L2范數(shù)正則化約束最小二乘法(附錄A),最終可以得到ICLS時(shí)頻分析結(jié)果:

he=S′w(SwS′w+μΩ)-1dw

(8)

式中:Ω表示約束對(duì)角矩陣[13];μ為權(quán)重系數(shù),用于調(diào)節(jié)約束的大小,實(shí)際求解過(guò)程中,可不斷調(diào)節(jié)μ的大小,直到獲得一個(gè)時(shí)間分辨率和頻率分辨率折中的時(shí)頻分析結(jié)果。加入μΩ約束分量后,對(duì)于信噪比較低、時(shí)間段較短的信號(hào)的頻譜分析會(huì)有更好的效果,可以有效地提高頻譜分析過(guò)程的穩(wěn)定性。逐點(diǎn)加窗計(jì)算(8)式即可得到待求信號(hào)段的時(shí)頻譜。

進(jìn)一步,可以通過(guò)公式(9)方便地實(shí)現(xiàn)ICLS反變換:

dinv=She

(9)

式中:dinv為經(jīng)過(guò)反變換后重構(gòu)的信號(hào)。

2 模型試算

設(shè)計(jì)一個(gè)由3個(gè)信號(hào)疊加而成的復(fù)合信號(hào),如圖1 所示,3個(gè)信號(hào)分別為平穩(wěn)單頻信號(hào)①、頻率隨時(shí)間變化的信號(hào)②和加高斯窗函數(shù)的單頻高頻信號(hào)③。

圖1 復(fù)合信號(hào)

采用本文方法(ICLS)對(duì)該復(fù)合信號(hào)進(jìn)行時(shí)頻分析,如圖2所示,從圖2中可以看出,時(shí)頻分析結(jié)果客觀、準(zhǔn)確地描述了該復(fù)合信號(hào)頻率成分的變化,分辨率高。圖3為復(fù)合信號(hào)的瞬時(shí)振幅及本文方法時(shí)頻譜的局部放大顯示,可以看出,時(shí)頻譜為零的地方對(duì)應(yīng)著復(fù)合信號(hào)瞬時(shí)振幅零點(diǎn)的位置,時(shí)頻譜非零值對(duì)應(yīng)著復(fù)合信號(hào)的瞬時(shí)振幅非零值的位置,也就是說(shuō),本文方法時(shí)頻譜不僅可以描述頻譜隨時(shí)間的變化,同時(shí)也描述了原始信號(hào)瞬時(shí)振幅的大小。圖4為本文方法反變換復(fù)合信號(hào)重構(gòu)結(jié)果,可以看出,重構(gòu)后復(fù)合信號(hào)與原始復(fù)合信號(hào)幾乎重合,說(shuō)明了本文方法的準(zhǔn)確性。

圖2 復(fù)合信號(hào)的ICLS時(shí)頻譜

圖3 復(fù)合信號(hào)瞬時(shí)振幅及其ICLS時(shí)頻譜

圖4 ICLS反變換

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

實(shí)際資料來(lái)自勝利油田某工區(qū),原始剖面如圖5所示,在1.7s左右、CDP100～CDP200處有一背斜油氣藏。分別采用本文方法(ICLS)和廣義S變換(STFT)[18]對(duì)該剖面進(jìn)行時(shí)頻分析,結(jié)果如圖6至圖8 所示。從圖6至圖8可以看出,ICLS和STFT的時(shí)頻分析結(jié)果具有較好的一致性,STFT的時(shí)頻分析結(jié)果在頻率較高的頻譜分量上分辨率也較高,但在低頻的頻譜分量上分辨率較低,這已經(jīng)無(wú)法滿足目前隱蔽油氣藏精細(xì)勘探的要求,相比而言,ICLS時(shí)頻分析結(jié)果的分辨率更高、準(zhǔn)確性更好,如對(duì)10Hz頻譜分量依然具有很好的聚焦能力。

圖8 STFT和ICLS的時(shí)頻譜分析結(jié)果(50Hz)

圖6 STFT和ICLS的時(shí)頻譜分析結(jié)果(10Hz)

圖5 原始地震剖面

圖9為ICLS時(shí)頻分析結(jié)果的10Hz頻譜分量除以50Hz頻譜分量的頻譜比剖面,箭頭處數(shù)值較大,可以明顯看出該背斜存在低頻增加、高頻衰減的現(xiàn)象,指示了一套氣藏的存在。

圖7 STFT和ICLS的時(shí)頻譜分析結(jié)果(30Hz)

圖9 頻譜比剖面(10Hz/50Hz)

4 結(jié)論

基于正則化約束最小二乘估計(jì)的思路,提出了一種改進(jìn)的最小二乘時(shí)頻分析方法。模型試算結(jié)果表明,采用本文方法可以得到準(zhǔn)確的、高分辨率的時(shí)頻分析結(jié)果;實(shí)際應(yīng)用結(jié)果表明,本文方法對(duì)于地震信號(hào)的高中低頻部分都有很高的分辨率和穩(wěn)定性,可以有效地提高儲(chǔ)層預(yù)測(cè)的準(zhǔn)確度和可靠性,在隱蔽油氣藏精細(xì)勘探中具有廣闊的應(yīng)用前景。