基于l1-l2范數(shù)的高分辨率時(shí)頻分析方法及應(yīng)用

邢文軍，曹思遠(yuǎn)，陳思遠(yuǎn)，馬敏瑤

中國(guó)石油大學(xué)(北京)，北京 102249

0 引言

時(shí)頻分析是地震資料處理中的重要分析方法，通過對(duì)地震數(shù)據(jù)進(jìn)行時(shí)頻分解，可獲得每一時(shí)刻的頻率成分；由于時(shí)頻分析方法存在海森堡測(cè)不準(zhǔn)原理的約束(Mallat，1999)，即時(shí)間、頻率的分辨率不能同時(shí)最優(yōu)，使得帶窗口類的時(shí)頻分析方法應(yīng)用受到了限制，包括短時(shí)傅里葉變換(Portnoff，1980；Margrave and Lamoureux，2001)、小波變換(Daubechies，1990；Chakraborty and Okaya，1995)、S變換(Stockwell et al., 1996; 高靜懷等, 2003)等.在此基礎(chǔ)上，受經(jīng)驗(yàn)?zāi)B(tài)分解算法啟發(fā)，基于帶窗口類的時(shí)頻分析可求取瞬時(shí)頻率、通過能量擠壓的方式提高頻率分辨率，這一類算法統(tǒng)稱同步擠壓變換，根據(jù)窗口類時(shí)頻分析方法的不同，具體可命名為同步擠壓短時(shí)傅里葉變換(Oberlin et al., 2014; Wang et al., 2014; Mahdavi et al., 2021)、同步擠壓小波變換(Chen et al., 2014; 潘曉等, 2020)等.

基于瞬時(shí)頻率的思想，認(rèn)為地震信號(hào)是由多個(gè)窄帶信號(hào)(模態(tài))相加而成，代表時(shí)頻分析方法包括經(jīng)驗(yàn)?zāi)B(tài)分解算法(Huang et al., 1998; Han and van der Baan, 2013)、集合經(jīng)驗(yàn)?zāi)B(tài)分解(Wu and Huang, 2009; Wang et al., 2012)、變分模態(tài)分解(Torres et al., 2011；龍丹等，2020；鄔蒙蒙等，2020)等，然后對(duì)分解后的固有模態(tài)分量進(jìn)行Hilbert變換，最后將變換后的瞬時(shí)振幅、瞬時(shí)頻率、瞬時(shí)相位等信息排列在時(shí)頻平面上.一般來說，除含加性噪聲的地震數(shù)據(jù)外，地震數(shù)據(jù)通常呈現(xiàn)不可分離的狀態(tài)，即并不是由多個(gè)分量組成，且該類算法是數(shù)據(jù)驅(qū)動(dòng)的自適應(yīng)分解，穩(wěn)定性較差.

同樣基于瞬時(shí)頻率的想法，稀疏反演時(shí)頻分析方法認(rèn)為在時(shí)頻平面上，每一時(shí)刻振幅譜能量的極大值點(diǎn)是瞬時(shí)頻率，然后通過稀疏約束反演的方法得到極大值點(diǎn)；代表方法包括以l1范數(shù)約束的短時(shí)傅里葉變換(Chen et al., 2020)、lp范數(shù)約束的短時(shí)傅里葉變換(Wang et al., 2020)等，這類算法的時(shí)頻分辨率較高，目前這一類算法仍屬于探索階段(田琳和胡津健，2021；楊子鵬等，2021).

對(duì)于稀疏反演類時(shí)頻分析方法，稀疏約束能力是決定時(shí)頻分辨率的關(guān)鍵.已知稀疏能力最好的范數(shù)是l0范數(shù)，但是其求解困難(Natarajan，1995)，在實(shí)際中一般使用l0范數(shù)的最優(yōu)凸近似——l1范數(shù)代替求解(Chen et al.，1998)，這將導(dǎo)致稀疏約束能力的下降，繼而降低了時(shí)頻分辨率；lp擬范數(shù)(Chartrand and Yin，2008)、l1-l2范數(shù)(Yin et al., 2015)、lp-l1范數(shù)(Zhao et al., 2020)等作為非凸范數(shù)同樣具備稀疏約束能力，其中l(wèi)1-l2范數(shù)的稀疏約束能力已被證明高于lp擬范數(shù)，已被廣泛應(yīng)用于地震數(shù)據(jù)處理中(Wang et al., 2018, 2019).

本研究提出基于l1-l2范數(shù)稀疏反演時(shí)頻分析(L12-STFT)方法.L12-STFT通過短時(shí)傅里葉變換的逆變換構(gòu)造目標(biāo)方程，以l1-l2范數(shù)作為稀疏約束，通過交替方向乘子法(ADMM)進(jìn)行求解，獲得高分辨率的時(shí)頻譜.模型部分證明了L12-STFT在時(shí)頻分析分辨能力上的優(yōu)越性，同時(shí)通過加噪數(shù)據(jù)的測(cè)試，證明L12-STFT具備一定抗噪性.實(shí)際數(shù)據(jù)部分基于L12-STFT進(jìn)行譜分解，計(jì)算縱波頻散屬性，精確刻畫儲(chǔ)層.

1 理論

1.1 稀疏反演時(shí)頻分析方法

本部分以短時(shí)傅里葉變換(STFT)為基礎(chǔ)，推導(dǎo)基于稀疏反演的時(shí)頻分析方法.

STFT首先將離散信號(hào)s∈N×1分解為N個(gè)長(zhǎng)度為M(M為奇數(shù)，M

yi=Gsi，(1)

其中G∈M×M表示對(duì)角線為高斯窗函數(shù)g∈M ×1的矩陣.si是以原信號(hào)s的第i個(gè)點(diǎn)為中心，信號(hào)兩邊各取(M-1)/2個(gè)數(shù)據(jù)作為子信號(hào).同時(shí)為了減弱邊緣效應(yīng)，原始信號(hào)兩端需補(bǔ)長(zhǎng)度為(M-1)/2個(gè)零.

假設(shè)加窗后子信號(hào)yi的傅里葉變換為xi∈N×1，則xi和yi可表示為：

yi≈SF-1xi，(2)

其中，S表示截?cái)嗑仃嚕?/p>

S=[I|O]，(3)

I∈M×M為單位矩陣，O∈M×(N-M)為零矩陣.S矩陣的作用為截取F-1xi的前M個(gè)點(diǎn)，即去除無效值，F(xiàn)表示傅里葉變換矩陣，其形式為：

(4)

令A(yù)=SF-1，即A表示部分傅里葉變換矩陣；則每一時(shí)刻的稀疏的反演方程為：

(5)

式中，λ表示正則化算子，用以調(diào)整時(shí)頻譜稀疏度，系數(shù)1/2作用在于簡(jiǎn)化求導(dǎo)過程.遍歷信號(hào)si(i=1,2,…,N)即可得到信號(hào)s∈N×1的時(shí)頻分析.

1.2 基于l1-l2范數(shù)的高分辨率時(shí)頻分析及求解策略

式(5)使用l0范數(shù)作為稀疏約束，其求解是NP-Hard問題，通常使用l1、lp等范數(shù)近似求解，本文使用l1-l2范數(shù)代替l0范數(shù)進(jìn)行求解，圖1為l2范數(shù)、l1范數(shù)、lp范數(shù)(p=0.65)和l1-l2范數(shù)的相平面圖，顯然，l1-l2與坐標(biāo)軸近似度更高，更易獲得稀疏解，相比于其他三種范數(shù)可以更好的近似l0范數(shù).

基于上述討論，修改式(5)為：

(6)

其中，α為加權(quán)系數(shù).式(6)表示l1-l2范數(shù)約束下的稀疏時(shí)頻分析優(yōu)化方程.正則化參數(shù)λ和α共同參與調(diào)節(jié)時(shí)頻譜的稀疏度，兩者越大，獲得的時(shí)頻譜稀疏度越高.反之兩者越小，時(shí)頻譜稀疏度越低，值得說明，當(dāng)α減小為α=0時(shí)，l1-l2范數(shù)退化為l1范數(shù).

式(6)可以使用凸差算法(DCA)和交替方向乘子法(ADMM)聯(lián)合求解(Ma et al., 2017)，這種求解方案在Wang等(2018)的研究中被得到了證實(shí).由于凸差算法迭代速度較慢，本文建議只使用交替方向乘子法(ADMM)進(jìn)行式(6)的求解(Lou and Yan, 2018)，具體求解方式如下：

圖1 稀疏約束的相平面圖(a) l2范數(shù)； (b) l1范數(shù)； (c) lp范數(shù)(p=0.65)； (d) l1-l2范數(shù).Fig.1 Phase plane of sparse constraint(a) l2-norm; (b) l1-norm; (c) lp-norm (p=0.65); (d) l1-l2 norm.

分裂變量xi=zi，預(yù)定義迭代步長(zhǎng)ρ，引入二次懲罰項(xiàng)，修改目標(biāo)函數(shù)(6)為：

(7)

式(7)需要兩個(gè)變量交替進(jìn)行求解，則xi更新為：

(8)

式(8)只涉及l(fā)2范數(shù)可使用梯度下降法求解：

(9)

Z0更新為：

(10)

需注意，當(dāng)‖y‖∞=λ/ρ時(shí)，存在無窮多個(gè)解，但在實(shí)際處理中，因?yàn)閿?shù)值的離散性，這種情況很難產(chǎn)生.

最后，對(duì)偶變量u使用對(duì)偶上升法進(jìn)行更新，相應(yīng)的更新迭代方程為：

(11)

1.3 地震頻散屬性

地震波在傳播過程中的衰減往往伴隨著頻散，地層的速度與頻率有關(guān)，則反射界面處的反射系數(shù)也與頻率相關(guān)，即地震縱波頻散屬性可用作地層含氣性的識(shí)別.

小波變換最早被應(yīng)用于縱波頻散的譜分解中，結(jié)合Smith和Gidlow方程有效描述了儲(chǔ)層流體特征(Wilson et al., 2009; Wilson, 2010)；近些年，隨著高分辨率時(shí)頻分析方法的不斷革新，Wigner-Ville分布(Wigner, 1932)、匹配追蹤(Mallat and Zhang, 1993)、VMD(Liu et al., 2016)、反演譜分解(黃廣譚等, 2017)方法被用于縱波頻散計(jì)算中，均不同程度提高了儲(chǔ)層流體的識(shí)別精度，本文中所提出的L12-STFT也將被用于頻散屬性的計(jì)算.

Smith和Gildlow近似方程為：

(12)

式中，θi為入射角，VP和VS分別為縱波速度和橫波速度，ΔVP和ΔVS分別為縱波速度變化量和橫波速度變化量.

(13)

式中，系數(shù)A(θi)和B(θi)與入射角和速度有關(guān)，可以通過射線追蹤計(jì)算.

根據(jù)Chapman等(2006)的頻散介質(zhì)理論，可以得到具有頻率依賴性的AVO近似方程，且對(duì)每個(gè)時(shí)刻t：

(14)

(15)

使用參數(shù)Ia和Ib表示反射率縱橫波頻散，即:

(16)

(17)

對(duì)于無頻散情況，(14)式可修改為：

(18)

地震數(shù)據(jù)S在時(shí)頻域可寫為反射系數(shù)和子波W(θi,f)的乘積，即：

S(t,θi,f)=W(θi,f)[R(t,θi,f0)+α(t,θi,f)]，(19)

其中，α(t,θi,f)=(f-f0)A(θi)Ia+(f-f0)B(θi)Ib，為反射系數(shù)的頻散項(xiàng)，與反射系數(shù)縱橫波頻散項(xiàng)Ia和Ib相關(guān).無反射系數(shù)頻散時(shí)：

S(t,θi,f0)=W(θi,f0)R(t,θi,f0).

(20)

式(19)和式(20)說明，與時(shí)頻域頻譜白化類似，時(shí)頻分析分辨率越高，相應(yīng)的反射系數(shù)和地震數(shù)據(jù)的在時(shí)頻域的特征越清晰，越有利于計(jì)算頻散參數(shù)，將式(20)代入式(19)中，消除R(t,θi,f0)，即：

S(t,θi,f)W(θi,f0)-S(t,θi,f0)W(θi,f)=

W(θi,f)W(θi,f0)α(t,θi,f)，(21)

式中，α(t,θi,f)與縱橫波頻散項(xiàng)Ia和Ib相關(guān)，可通過求解式(21)獲得縱橫波頻散參數(shù).

1.4 基于l1-l2范數(shù)的縱波頻散參數(shù)計(jì)算流程

根據(jù)上述縱波頻散屬性的推導(dǎo)，基于l1-l2范數(shù)的縱波頻散參數(shù)計(jì)算流程如下：

(1)對(duì)疊前地震數(shù)據(jù)進(jìn)行幾何擴(kuò)散補(bǔ)償、地表一致性振幅補(bǔ)償、噪聲衰減、動(dòng)校正等處理，要求盡可能保幅處理.然后劃分角道集、疊加，獲得角道集疊加數(shù)據(jù).

(2)提取單道地震數(shù)據(jù)測(cè)試基于l1-l2范數(shù)的時(shí)頻分析方法的參數(shù)，然后對(duì)地震數(shù)據(jù)進(jìn)行分頻，獲得不同角度、頻率的分頻剖面S(t,θi,f)和S(t,θi,f0).

(3)工區(qū)內(nèi)含測(cè)井?dāng)?shù)據(jù)時(shí)，使用井?dāng)?shù)據(jù)和各角度道集提取角度子波，選取子波主頻作為參考頻率f0，計(jì)算子波振幅譜W(θi,f).工區(qū)內(nèi)不含測(cè)井?dāng)?shù)據(jù)時(shí)可直接使用地震數(shù)據(jù)振幅譜的包絡(luò)作為子波的振幅譜.

(4)基于式(21)進(jìn)行縱波頻散屬性計(jì)算，從而準(zhǔn)確指示流體.

2 模型測(cè)試

本部分測(cè)試基于l1-l2范數(shù)時(shí)頻分析方法的聚焦性，并測(cè)試其抗噪能力.

模型一合成采樣頻率1024 Hz的調(diào)頻信號(hào)(圖2)，共1024個(gè)采樣點(diǎn)，調(diào)頻信號(hào)s(t)的表達(dá)式為：

(22)

圖2 調(diào)頻信號(hào)Fig.2 Original chirp signals

分別使用短時(shí)傅里葉變換(STFT)、同步擠壓短時(shí)傅里葉變換(SST)、l1范數(shù)約束的反演時(shí)頻分析(L1-STFT)和L12-STFT對(duì)上述調(diào)頻信號(hào)進(jìn)行時(shí)頻分析，分析結(jié)果如圖3所示，對(duì)比四種時(shí)頻分析方法，L12-STFT的時(shí)頻聚焦性最好(圖3d)，可定性的反映瞬時(shí)頻率、瞬時(shí)振幅等信息；L1-STFT由于l1范數(shù)稀疏約束能力欠佳，時(shí)頻譜的聚焦性(圖3c)也弱于L12-STFT(圖3d).

圖3 調(diào)頻信號(hào)的時(shí)頻分析(a) STFT； (b) SST； (c) L1-STFT； (d) L12-STFT.Fig.3 Time-frequency spectrum of original chirp signals

模型二與模型一中的信號(hào)解析式相同，加噪后信噪比為6.488 dB，模型二信號(hào)如圖4所示，其波形已被大量高斯隨機(jī)噪聲破壞；圖5為四種方法對(duì)含噪調(diào)頻信號(hào)的時(shí)頻分析.綜合四種方法分析，圖5d中的L12-STFT有較高的抗噪性，依然可以良好的展示含噪調(diào)頻信號(hào)的時(shí)頻特征，主要原因是稀疏反演類時(shí)頻分析方法使用l2范數(shù)作為擬合項(xiàng)，并使用稀疏范數(shù)作為正則化項(xiàng)，可壓制高斯分布的隨機(jī)噪聲(Sun et al., 2021)，且L12-STFT的稀疏約束能力較強(qiáng).與之原理相近，由于l1范數(shù)稀疏能力較弱，L1-STFT在時(shí)頻平面上出現(xiàn)很多噪聲能量(圖5c)；如圖5b所示，SST得到時(shí)頻譜的原理是求取瞬時(shí)頻率后的能量重分配，在瞬時(shí)頻率準(zhǔn)確求取的前提下，可有較好的抗噪性，但當(dāng)噪聲嚴(yán)重時(shí)，噪聲被SST認(rèn)為是有效能量，繼而給予能量分配.故而相比于L12-STFT，SST抗噪性有限，在時(shí)頻平面上也有噪聲能量出現(xiàn).圖5a是STFT的時(shí)頻譜，其不具備抗噪性，受噪聲影響較嚴(yán)重.

圖4 含噪調(diào)頻信號(hào)Fig.4 Noisy chirp signals

圖5 含噪調(diào)頻信號(hào)的時(shí)頻分析(a) STFT； (b) SST； (c) L1-STFT； (d) L12-STFT.Fig.5 Time-frequency spectrum of noisy chirp signals

模型三使用時(shí)間采樣間隔為1 ms，主頻50 Hz的Ricker合成地震記錄測(cè)試時(shí)頻分析的性能(圖6)，信號(hào)含隨機(jī)噪聲，信噪比為13.53 dB.對(duì)比方法包括基于lp范數(shù)的稀疏時(shí)頻分析，即Lp-STFT(p=0.65)，其中，短時(shí)傅里葉變換為窗長(zhǎng)21ms的高斯窗，其余四種時(shí)頻分析方法高斯窗的窗長(zhǎng)為11 ms.測(cè)試結(jié)果如圖7所示：STFT的時(shí)頻分辨率仍然最低(圖7a)；地震子波的同步擠壓變換(圖7b)表現(xiàn)為“線”的形式，而反演類稀疏時(shí)頻分析方法則表現(xiàn)為“抖動(dòng)的能量團(tuán)”(圖7c—e)，且能量團(tuán)的時(shí)頻聚焦性與正則化項(xiàng)的稀疏約束能力成正比，即L12-STFT時(shí)頻聚焦性 >Lp-STFT>L1-STFT.

圖6 反射系數(shù)(虛線)及合成地震記錄(實(shí)線)Fig.6 Reflectivity (dotted line) and synthetic seismogram (solid line)

圖7 合成地震記錄的時(shí)頻分析(a) STFT； (b) SST； (c) L1-STFT； (d) Lp-STFT(p=0.65)； (e) L12-STFT.Fig.7 Time-frequency spectrum of synthetic seismogram

為了測(cè)試算法的抗噪性，使用時(shí)頻聚焦性的表征參數(shù)Renyi熵測(cè)試不同信噪比(SNR)的算法性能，同一信噪比重復(fù)測(cè)試20次取平均值.Renyi值越小，表示時(shí)頻聚焦性越好.測(cè)試結(jié)果如圖8所示.測(cè)試結(jié)果表明，隨信噪比降低，參與測(cè)試的時(shí)頻分析方法的時(shí)頻聚焦性均變差.但SST和稀疏反演類時(shí)頻分析方法(L1-STFT、Lp-STFT、L12-STFT)的變化差異較小，且本文建議的L12-STFT的Renyi熵始終最小，這說明建議的方法在不同信噪比的數(shù)據(jù)測(cè)試中，均有較高的時(shí)頻聚焦性.

圖8 Renyi熵隨信噪比的變化曲線Fig.8 Variation curves of Renyi entropy with SNR

3 實(shí)際數(shù)據(jù)測(cè)試

本部分采用含氣層的實(shí)際數(shù)據(jù)測(cè)試計(jì)算縱波頻散屬性，數(shù)據(jù)采樣間隔為1 ms，目標(biāo)儲(chǔ)層約為2000 ms處，主要巖性為砂泥巖.圖9為實(shí)際數(shù)據(jù)的近(圖9a)、中(圖9b)、遠(yuǎn)(圖9c)三個(gè)角度疊加剖面，其中第293道存在橫波測(cè)井?dāng)?shù)據(jù)，如圖9所示，黑色曲線所示為縱橫波速度比曲線.

圖9 部分角道集疊加剖面Fig.9 Partial stacked sections of angle gather

首先對(duì)過井地震道進(jìn)行時(shí)頻分析測(cè)試(圖10).圖11為圖10的時(shí)頻譜，四種方法都可以在時(shí)頻域上體現(xiàn)地震道的響應(yīng)特征，如2000 ms的強(qiáng)振幅特征等，并且圖11d所示的L12-STFT的時(shí)頻聚焦性優(yōu)于其余三種方法(圖11a、b、c).

圖10 過井地震道Fig.10 Cross-well seismic trace

圖11 過井地震道的時(shí)頻分析(a) STFT； (b) SST； (c) L1-STFT； (d) L12-STFT.Fig.11 Time-frequency spectrum of cross-well seismic trace

基于縱波速度(VP)，橫波速度(VS)和密度并利用Zoeppritz方程正演角道集，正演所使用的Ricker子波主頻與地震數(shù)據(jù)主頻相同(圖12).通過合成的角道集可以測(cè)試該數(shù)據(jù)的頻散屬性響應(yīng)特征.其中在1900～1950 ms和1950～2000 ms處呈現(xiàn)儲(chǔ)層AVO響應(yīng)(黑色虛線框)，結(jié)合測(cè)井曲線分析，1900～1950 ms處為砂泥巖互層，泥巖含量偏高，儲(chǔ)層質(zhì)量差.1975～2000 ms處為砂巖氣藏，上覆1975 ms為高GR的泥巖蓋層，滿足油氣的儲(chǔ)存條件；分別基于STFT、L1-STFT和L12-STFT對(duì)圖12的角道集進(jìn)行頻散參數(shù)的計(jì)算，如圖13所示，不同方法計(jì)算得到的頻散曲線均能反應(yīng)儲(chǔ)層，且L12-STFT的分辨率較高(黑色虛線框).

圖12 Zoeppritz方程正演的角道集Fig.12 Synthetic angle gather using Zoeppritz equation

對(duì)實(shí)際角道集疊加數(shù)據(jù)(圖9)進(jìn)行頻散參數(shù)測(cè)試，三種方法的頻散屬性曲線上，在1975～2000 ms處的位置均產(chǎn)生異常，并且L1-STFT和建議的L12-STFT的分辨率最高，SST次之，STFT的頻散屬性分辨率最差(圖14黑色虛線框).圖15是這三種方法形成頻散屬性的剖面，基于L12-STFT計(jì)算的縱波頻散屬性分辨率最高，相比于L1-STFT，其頻散屬性剖面更整潔，能清晰描繪儲(chǔ)層邊界，刻畫儲(chǔ)層特征(圖15黑色虛線框).需要指出，實(shí)際數(shù)據(jù)的頻散屬性與正演模擬的屬性存在差異，原因在于實(shí)際數(shù)據(jù)中存在噪聲、角道集劃分準(zhǔn)確度差、地震子波精度低等情況.

4 結(jié)論

本研究提出一種基于l1-l2范數(shù)約束的時(shí)頻分析方法，該方法受益于l1-l2范數(shù)強(qiáng)大的稀疏約束能力，聯(lián)合短時(shí)傅里葉變換，可獲得高分辨率的時(shí)頻譜.縱波頻散屬性可以較好的指示流體，我們將基于l1-l2范數(shù)的高精度時(shí)頻分析方法與頻散屬性相結(jié)合，提高了現(xiàn)有方法的縱波頻散的分析精度.模型和實(shí)際數(shù)據(jù)均表明，本文提出的時(shí)頻分析方法具有較高的時(shí)頻分辨率，并且可以適用于地震數(shù)據(jù)分析.

圖13 圖11所示角道集計(jì)算的頻散屬性Fig.13 Calculated dispersion property of angle gather in Fig.11

圖14 測(cè)井曲線和所計(jì)算的過井頻散屬性曲線Fig.14 Well-log curves and calculated dispersion attribute curves

圖15 縱波頻散屬性剖面(a) STFT； (b) SST； (c) L1-STFT； (d) L12-STFT.Fig.15 P-wave dispersion attribute