基于l1-l2范數的高分辨率時頻分析方法及應用

邢文軍，曹思遠，陳思遠，馬敏瑤

中國石油大學(北京)，北京 102249

0 引言

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

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

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

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

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

1 理論

1.1 稀疏反演時頻分析方法

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

STFT首先將離散信號s∈N×1分解為N個長度為M(M為奇數，M

yi=Gsi，(1)

其中G∈M×M表示對角線為高斯窗函數g∈M ×1的矩陣.si是以原信號s的第i個點為中心，信號兩邊各取(M-1)/2個數據作為子信號.同時為了減弱邊緣效應，原始信號兩端需補長度為(M-1)/2個零.

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

yi≈SF-1xi，(2)

其中，S表示截斷矩陣：

S=[I|O]，(3)

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

(4)

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

(5)

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

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

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

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

(6)

其中，α為加權系數.式(6)表示l1-l2范數約束下的稀疏時頻分析優化方程.正則化參數λ和α共同參與調節時頻譜的稀疏度，兩者越大，獲得的時頻譜稀疏度越高.反之兩者越小，時頻譜稀疏度越低，值得說明，當α減小為α=0時，l1-l2范數退化為l1范數.

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

圖1 稀疏約束的相平面圖(a) l2范數； (b) l1范數； (c) lp范數(p=0.65)； (d) l1-l2范數.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，預定義迭代步長ρ，引入二次懲罰項，修改目標函數(6)為：

(7)

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

(8)

式(8)只涉及l2范數可使用梯度下降法求解：

(9)

Z0更新為：

(10)

需注意，當‖y‖∞=λ/ρ時，存在無窮多個解，但在實際處理中，因為數值的離散性，這種情況很難產生.

最后，對偶變量u使用對偶上升法進行更新，相應的更新迭代方程為：

(11)

1.3 地震頻散屬性

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

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

Smith和Gildlow近似方程為：

(12)

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

(13)

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

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

(14)

(15)

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

(16)

(17)

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

(18)

地震數據S在時頻域可寫為反射系數和子波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，為反射系數的頻散項，與反射系數縱橫波頻散項Ia和Ib相關.無反射系數頻散時：

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

(20)

式(19)和式(20)說明，與時頻域頻譜白化類似，時頻分析分辨率越高，相應的反射系數和地震數據的在時頻域的特征越清晰，越有利于計算頻散參數，將式(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)與縱橫波頻散項Ia和Ib相關，可通過求解式(21)獲得縱橫波頻散參數.

1.4 基于l1-l2范數的縱波頻散參數計算流程

根據上述縱波頻散屬性的推導，基于l1-l2范數的縱波頻散參數計算流程如下：

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

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

(3)工區內含測井數據時，使用井數據和各角度道集提取角度子波，選取子波主頻作為參考頻率f0，計算子波振幅譜W(θi,f).工區內不含測井數據時可直接使用地震數據振幅譜的包絡作為子波的振幅譜.

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

2 模型測試

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

模型一合成采樣頻率1024 Hz的調頻信號(圖2)，共1024個采樣點，調頻信號s(t)的表達式為：

(22)

圖2 調頻信號Fig.2 Original chirp signals

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

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

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

圖4 含噪調頻信號Fig.4 Noisy chirp signals

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

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

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

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

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

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

3 實際數據測試

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

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

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

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

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

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

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

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

4 結論

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

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

圖14 測井曲線和所計算的過井頻散屬性曲線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