中石化勘探分公司勘探研究院，四川 成都

收稿日期：2023年5月5日；录用日期：2023年5月31日；发布日期：2023年6月13日

摘要

近年来，致密油气藏勘探形式发展迅猛，部分井上钻遇到了一些前期目标评价未有效识别的河道砂体，这类砂体被证实也具有良好勘探潜力，如何有效识别隐蔽类河道成为气田增储上产急需解决的突出问题。由于谱分解技术可从全频地震资料中剥离出反映河道地质特征的特定频段信息，因此该技术有助于河道识别。然而，传统时频分析方法分频结果分辨率较低，难以清晰刻画出隐蔽河道边界。为此，本文引入了LR算法同步挤压广义S变换进行地震时频谱分解。数值模拟表明，同步挤压算法极大地提高了广义S变换时频分辨率，较好地区分开了非平稳信号中不同频率信号分量，更适应于实际地震信号的时频谱分解；而LR算法则进一步消除了信号间存在的虚假交叉项。实际资料应用表明，本文提出的谱分解方法精确刻画出了川东南地区凉高山组一期“隐蔽”河道砂体的空间展布，证实了方法的有效性和优越性。

关键词

同步挤压，时频分析，致密砂岩，边界检测，“隐蔽”型河道识别

Hidden Channel Sand Bodies Identification Based on the Synchrosqueezing Generalized S-Transform with the Lucy-Richardson Algorithm

Yijiang Zhang, Qianshu Chen, Yuxi Zhang, Meng Yuan, Yuan Zhou, Yu Jiang

Exploration Research Institute of Sinopec Exploration Company, Chengdu Sichuan

Received: May 5^th, 2023; accepted: May 31^st, 2023; published: Jun. 13^th, 2023

ABSTRACT

In recent years, the exploration of tight oil and gas reservoirs has developed rapidly. Some wells have encountered a few channel sand bodies that were not effectively identified by early target evaluations, and these sand bodies have been proven to have good exploration potential. How to effectively identify hidden channels has become a prominent problem that urgently needs to be solved for gas field storage and production increase. Due to the fact that spectral decomposition technology can extract specific frequency band information reflecting the geological characteristics of channels from full frequency seismic data, it is helpful for channel sand bodies identification. However, traditional time-frequency analysis methods have low resolution in frequency decomposition results, making it difficult to clearly detect hidden channel boundaries. Therefore, this article introduces the synchrosqueezing generalized S-transform with the Lucy-Richardson algorithm for seismic time-frequency spectrum decomposition. Numerical simulation shows that the synchrosqueezing algorithm greatly improves the time-frequency resolution of the generalized S-transform, effectively distinguishes different frequency signal components in non-stationary signals, and it is more suitable for the time-frequency decomposition of actual seismic signals. The Lucy-Richardson algorithm further eliminates false cross terms between signals. The practical application of data shows that the spectral decomposition method proposed in this paper accurately characterizes the spatial distribution of a hidden channel sand body in the Lianggaoshan Formation in southeastern Sichuan, which confirming the effectiveness and superiority of the method.

Keywords:The Synchrosqueezing Transform, Time-Frequency Analysis, Tight Sandstone, Boundary Detection, Identification of “Hidden” Channel Sand Bodies

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

川东南凉高山组河道砂岩发育，其整体以湖泊–三角洲沉积体系为主。前期研究证实，凉高山组不同层段的沉积环境不同，河道砂的地震响应特征也不相同。目前主要探明的致密气藏以“亮点”型河道砂岩气为主，该类河道砂岩的纵波速度、密度明显高于周围泥岩，响应特征以中–低频强波峰波谷反射为主，地震上易于识别。除“亮点”型河道外，王荐等(2022)研究发现地震上也存在中弱振幅响应的“暗点”型河道砂体，其AVO特征以Ⅰ类和Ⅱ类为主，应用常规地球物理手段难以有效识别该类河道砂体，为此他们引入优势入射角部分叠加技术实现了隐蔽河道的有效识别 [1] 。郑公营等(2022)利用谱反演技术对优势道集叠加数据开展拓频处理，进一步提升了隐蔽河道砂的识别精度 [2] 。事实上，隐蔽型河道识别如此困难，不光因为砂体自身AVO特征形成的弱反射，还可能由于其响应特征遭到外界其他因素的屏蔽 [3] 。本文研究区凉二下亚段就沉积一期三角洲水下分流河道，但由于最大湖泛面影响，该期河道的识别受到严重制约。

时频分析技术是地震资料处理解释中常用方法，同时也被证明是识别砂岩河道的一种有效手段。该技术主要依据地层反射调谐原理，若对地震数据进行时频分解，便可产生一系列单频体，在不同单频体上就可以看到薄层的干涉特征，进而描述地层的形态和变化(赵秀红，2007) [4] 。肖高杰等(2008)在基于常规阻抗反演等方法效果不佳的情况下，应用Wigner-Ville分布(WVD)分频技术准确识别到了W区河道沉积体 [5] 。王红岩等(2015)指出传统河道识别方法是根据全频段地震数据进行研究，而全频地震数据包含的信息过于繁杂，不同频率信号相互压制或抵消，从而降低了对目标地质体的预测精度。他们通过短时傅里叶变换(STFT)生成振幅频率体和相位体来检测薄层及地质异常体横向不连续性，提高了储层预测的精度 [6] 。分频技术有助于河道的识别，但所选用的时频分析方法对分频结果的分辨率影响较大。牛双晨等(2021)分析并认为STFT、连续小波变换(CWT)、S变换(ST) [7] 和广义S变换(GST) [8] [9] [10] 等线性类时频分析方法，计算效率较高，但对窗函数依赖性较高，难以兼顾时间–频率分辨率。为此他们提出Wigner-Ville分布(WVD)与最大熵方法(MEM)相结合的非平稳信号分析方法(WVD-MEM)，该方法对河道细节的刻画具有一定改善 [11] 。WVD作为一种非线性时频分析方法，摆脱了测不准原理的制约，但对于复杂信号会产生交叉项干扰。Lu等(2009)提出反褶积短时傅里叶变换(DSTFT)，不仅提高了时–频分辨率，且减少了交叉项干扰 [12] 。同步挤压变换(SST)由Daubechies和Maes (1996)提出，该算法能够将信号的频带汇聚在其中心频率附近更小的一段范围，因此其时频分布更够突破1/8波长，极大地提高了信号识别的精度 [13] 。由于SST核心还是对STFT时频结果进行重排，分析时窗相对固化，张懿疆等(2020)提出LR算法同步挤压广义S变换(LR-SSGST)，其分析时窗更加灵活，能够获得极高分辨率的时频结果，同时消除信号间交叉项干扰，适用于目前窄细或隐蔽类河道的识别 [14] 。

本文针对地震资料存在砂体边界被特殊地质背景屏蔽的情况，提出了一种基于分频技术的隐蔽河道识别方法，利用LR算法同步挤压广义S变换将全频地震信息时频分解，然后对分解结果进行扫描，优选出河道特征最明显的单频体开展河道刻画。由于该方法结合了SST高时频分辨率和WVD消除交叉项的优势，使得其提取的单频体对河道边界辨识能力更强。通过仿真信号和正演模型验证了该分频方法的有效性。最后，应用该技术成功在川东南地区刻画出一期隐藏在最大湖泛面期的水下分流河道。

2. 方法原理

通过数学推导，对广义S变换结果进行“挤压”和重排，再使用Lucy-Richardson算法对同步挤压广义S变换谱进行二维反褶积，从而得到地震信号精确的时频分布。其具体实现方法如下：

首先求取信号 $x\left(t\right)$ 的广义S变换，其变换公式为

$GS{T}_x\left(\tau ,f\right)={\displaystyle {\int }_{-\infty }^{\infty }x\left(t\right)}\frac{\lambda {\left|f\right|}^p}{\sqrt{2\pi }}{\text{e}}^{\frac{-{\lambda }^2{f}^{2p}{\left(t-\tau \right)}^2}{2}}{\text{e}}^{-i2\pi ft}\text{d}t,\lambda >0,p>0$ . (1)

$\lambda ,p$ 为广义S变换窗函数调节参数。当取 $p=0$ ，上式即为短时傅里叶变换；当 $\lambda =1,p=1$ 时，上式即为标准的S变换。由此可见短时傅里叶变换、S变换与广义S变换之间有着内在的联系，彼此可相互推导得出。

令 $\delta \left(f\right)=\lambda {\left|f\right|}^p$ ，其中 $\delta$ 可看成是决定窗函数时间宽度的尺度因子，上式可写为

$GST\left(\delta ,f,\tau \right)={\displaystyle {\int }_{-\infty }^{\infty }x\left(t\right)\frac{\delta }{\sqrt{2\pi }}{\text{e}}^{-\frac{{\delta }^2{\left(t-\tau \right)}^2}{2}}{\text{e}}^{-i2\pi f\left(t-\tau \right)}{\text{e}}^{-i2\pi f\tau }\text{d}t}$ . (2)

通过对广义S变换谱求导，估计瞬时频率 $\hat{f}$ 。信号的瞬时频率为

$\hat{f}\left(f,\tau \right)=f+{\left[i2\pi GST\left(f,\tau \right)\right]}^{-1}\frac{\partial GST\left(f,\tau \right)}{\partial \tau }$ . (3)

某一点信号的中心频率有 $\hat{f}\left(f,\tau \right)=f_s$ 。把中心频率 $f_s$ 附近区间 $\left[f_s-\frac{1}{2}{f}_m,f_s+\frac{1}{2}{f}_m\right]$ 内的频谱进行叠

加，并将叠加的值放到中心频率处，以获得同步挤压广义S变换，此举极大地提高了广义S变换的频率分辨率。其表达式为

$SSGS{T}_x\left(f_s,\tau \right)=\left(f_m^{-1}\right){\displaystyle {\sum }_{f_k:\left|f_s\left(f_m,\tau \right)-f_s\right|\le f_m/2}\left|GST\left(f,\tau \right)\right|f_k\Delta f_k}$ , (4)

$f_k$ 为GST的离散频率， $\Delta f_k=f_k-f_{k-1}$ ， $f_s$ 和 $\Delta f_s$ 分别为中心频率和挤压频带， $\Delta f_s=f_s-f_{s-1}$ 。这时我们可以得到同步挤压广义S变换(SSGST)谱，

$SPE{C}_x\left(f_s,\tau \right)={\left|SSGS{T}_x\left(f_s,\tau \right)\right|}^2$ . (5)

在时频分析中，我们知道受测不准原理制约，信号在时间和频率上的分辨率具有一定的矛盾性，不同频率信号之间存在交叉项，这使得正确识别不同频率分量的时频分布变得困难。为此我们在SSGST谱上应用了迭代的二维反卷积算法，如Lucy-Richardson算法，以消除交叉项，同时增强谱的能量均衡，避免挤压算法将谱的能量过于挤压到一个点上。同步挤压广义S变换谱可以写成以下二维褶积形式：

$SPE{C}_x={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }WV{D}_h\left(\tau ,\eta \right)\times WV{D}_x\left(t-\tau ,f-\eta \right)\text{d}\tau \text{d}\eta }}$ , (6)

其 中 $WV{D}_x$ 和 $WV{D}_h$ 分别为信号 $x\left(t\right)$ 和窗函数 $h\left(u\right)$ 的魏格纳分布。式(6)也可简写成：

$S_x=W_x\ast \ast W_h$ . (7)

其中 $S_x$ 是信号SSGST谱， $W_x$ 和 $W_h$ 分别 是 $WV{D}_x$ 和 $WV{D}_h$ 的简写， $\ast \ast$ 代表二维褶积。窗函数 $h\left(u\right)$ 我们选择高斯窗，则 $W_h$ 便可作为已知量，在已知 $S_x$ 和 $W_h$ 的情况下，我们可以反向计算出 $W_x$ 。

Lucy-Richardson算法可以表示为：

$W_x^{-}\left(k+1\right)=W_x^{-}\left(k\right)\left(W_h\ast \frac{S_x}{W_h\otimes W_x^{-}\left(k\right)}\right)$ . (8)

$k+1$ 是现在的迭代次数， $W_x^{-}\left(0\right)=S_x$ 。其中 $\ast$ 为相关算子， $\otimes$ 为褶积算子，反褶积运算抑制了交叉项的产生。

LR-SSGST高分辨时频属性分解的主要步骤如下：

1) 由式(1)对原始地震信号 $x\left(t\right)$ 做广义S变换，获得广义S变换谱 $GS{T}_x$ (为了抑制数据变换结果的端点效应，可对输入信号进行加窗处理，本文选择汉宁窗)；

2) 通过同步挤压算法将广义S变换谱 $GS{T}_x$ 转化为同步挤压广义S变换谱 $SSGS{T}_x$ ；

3) 将同步挤压广义S变换谱 $SSGS{T}_x$ 与高斯窗 $h\left(u\right)$ 代入式(8) Lucy-Richardson算法中进行反卷积迭代，消除信号间交叉项；

4) 在时频属性分解时，需对单道地震信号进行时频分析，确定目标地质体优势频率，进而获得LR-SSGST高分辨单频体。

需要注意的是，调节广义S变换双参数时会对广义S变换分析时窗造成影响，而调节反褶积运算中高斯时窗的宽度亦会对信号分析造成影响。随着时窗长度的增加，信号频率分辨率提高，时间分辨率降低；反之则频率分辨率降低，时间分辨率提高。因此，反褶积时窗长度选取不应过大或过小，应与实际储层厚度的采样点相一致。

3. 仿真实验

为对比LR-SSGST变换的时频分辨率与时频聚集性，采用一个合成的复杂信号进行验证。该信号包含了低、中、高三个频段，其中301~510 ms包含了20 Hz和100 Hz两个频率分量的信息。详细频率信息如表1；该信号的表达式为：

$x\left(t\right)=\{\begin{array}{ll}\cos \left(40\pi t/1024\right),\hfill & 0\le t\le 160\hfill \\ \cos \left(80\pi t/1024\right),\hfill & 161\le t\le 300\hfill \\ \cos \left(200\pi t/1024\right),\hfill & 301\le t\le 359\hfill \\ \cos \left(200\pi t/1024\right)+\cos \left(40\pi t/1024\right),\hfill & 360\le t\le 450\hfill \\ \cos \left(200\pi t/1024\right),\hfill & 451\le t\le 510\hfill \\ \cos \left(160\pi t/1024\right),\hfill & 511\le t\le 720\hfill \\ \cos \left(100\pi t/1024\right),\hfill & 721\le t\le 1024\hfill \end{array}$ (9)

图1为合成信号，图2为合成信号的短时傅里叶变换时频谱，图3、图4分别为合成信号的广义S变换时频谱和LR-SSGST变换时频谱。比较图2、图3、图4可看出：首先，图2、图3中箭头处端点效应明显，出现了部分假频；其次，两种频率成分交界的位置时频分辨率较差(图中虚线椭圆所示)。相对而言，LR-SSGST变换时频谱的整体分辨率较高，不同频率成分交界处边界清晰(图4)，频率聚集性相对优于STFT变换和广义S变换。因此，LR-SSGST变换对非平稳信号中不同分量的区分能力更强。

表1. 合成信号的时频分布

图1. 合成信号

图2. STFT变换时频谱

图3. 广义S时频谱

图4. LR-SSGST时频谱

本文根据研究区实际地质背景，设计了一个包含三种河道类型的地质模型，并利用该模型验证了不同时频分析方法分频效果的差异。该模型的道数(CDP)和持续时间分别为100和0.34秒(s)，时间采样间隔设置为2 ms，如图5(a)所示，黄色区域代表河道砂体的位置，其余区域代表泥岩位置。图5(b)为相应地质模型正演生成的合成地震剖面，Ⅰ类孤立型河道在地震记录上表现为断续低频强振幅响应。Ⅱ类复合型河道由两期并列河道组成，地震上表现为强振幅响应，但因为两期河道相距较近，河道间的过渡带振幅也较强。Ⅲ类隐蔽型河道由于其沉积位置接近不同沉积期分界面，河道地震响应被地质分界面的连续强振幅所屏蔽。图5(c)是对合成地震记录做GST得到的时频结果，三类河道的响应大体被区分开，但Ⅱ类复合河道中间的过渡带区分效果不明显(图5(c)中红色箭头所指)，这说明广义S变换频率分辨率较有限。图5(d)是对信号做LR-SSGST得到的时频分布，该结果具有较高的时频分辨率与频率汇聚性，三类河道的特征均能被有效识别出来。

图5. 三类河道模型及不同时频分析方法获得单频剖面。(a) 砂岩透镜体模型；(b) 合成地震记录；(c) 广义S变换单频剖面；(d) LR-SSGST单频剖面

图6. F1井凉高山组沉积层序综合柱状图

4. 实例分析

研究区位于川东南弧形高陡构造带，目的层为侏罗系中统凉高山组，按岩性凉高山组可分为凉一、凉二和凉三段。钻井及地球物理预测显示，该区凉二中、下亚段发育两期河道砂组，河道规模相对较大，凉二下亚段烃源泥页岩也较为发育。按照“贴近烃源岩，寻找河道砂”的思路，近期部署实施了专探凉二中亚段③砂组的F1井，该砂组紧贴烃源岩上部沉积(图6)，源储关系较好，经测试该井获产稳定工业气流，取得了该区凉高山组致密气勘探的突破。从源储配置来看，②砂组夹在两套烃源岩之间(图6红色方框所示)，源储配置比③砂组更为有利，但受最大湖泛面地质环境影响，平面上难以有效刻画②砂组的展布。为解决该问题本文利用LR算法同步挤压广义S变换分频技术，来去除强反射背景的屏蔽，提高河道砂岩辨识度，从而实现对“隐蔽”型河道砂岩的识别。

本文首先分析了不同实际地震道的时频特征(图7)。从原始地震剖面中可以看出，J₂l²²层地震同相轴振幅能量整体较强，横向变化不大，而在时频谱上，第1435道和第1605道在J₂l²²层振幅能量较弱，第1475道和第1555道在J₂l²²层振幅能量较强(图中红色方框段对应地震剖面J₂l²²层)，这说明在时间域地质体的地震响应差异不明显时，在频率域其响应特征较突出，利用该特性可帮助研究人员在隐蔽条件下识别特殊目标。

图7. 不同实际地震道时频分析

接下来对研究区地震资料振幅–频率谱进行分析，明确资料主频及频宽后以便提取有效频带范围内的地震属性或分频数据体，可以看出工区地震数据主频约40 Hz，有效频带范围为10~70 Hz (图8)。F1井③砂组为一套高速、高波阻抗砂岩，与上下覆围岩形成泥包砂结构，因此在地震剖面上③砂组表现为典型的断续“亮点”响应特征，②砂组具有与③砂组相似的岩性组合，但受最大湖泛面影响，其河道响应特征被连续强能量同相轴所屏蔽，在常规地震剖面难以识别其河道边界(图9，剖面位置如图10中蓝线所示)。在部分入射角叠加剖面上，河道的响应特征有所增强，但精细识别河道边界还较为困难(图9(b)、图9(c))。相比原始及部分角度叠加地震剖面，基于LR-SSGST分解的50 Hz单频剖面(图9(d))，仅保留了②砂组的砂体地震响应，可清晰识别出该期河道的边界。

图8. 研究区目的层地震资料频谱分析

图9. 过F1井凉二段全叠加、部分叠加与分频地震剖面对比。(a) 原始全叠加地震剖面；(b) 近偏移距叠加剖；(c) 远偏移距叠加剖；(d) LR-SSGST单频剖面

图10. 研究区目的层分频前后地震资料河道刻画对比图。(a) 原始属性；(b) 分频属性

分别利用J₂l²²层提取②砂组原始和分频地震数据振幅平面属性，由图10(a)可以看出，基于原始地震数据的平面属性上未有效识别出河道的上游展布(图中黑色箭头所指)；而分频地震数据的平面属性图对上下游河道刻画更完整且清晰(图10(b)中黑色箭头及椭圆框所指)。

5. 结论

1) 本文研究区凉二下亚段沉积一期三角洲水下分流河道，其地震响应特征应为断续强振幅的“亮点”，但由于砂体位置靠近最大湖泛面，导致河道边界特征被淹没在连续强振幅地震同相轴中，在常规叠后或叠前部分叠加数据上难以准确识别这期河道边界及展布。谱分解技术可以从全频地震数据中剥离出特定频段信息，而单频数据可突出砂体地震响应，有助于识别“隐蔽”类河道。

2) 利用谱分解技术进行地震分频时，所选用的时频分析方法对分频结果的分辨率影响较大，为提高河道识别精度，本文采用Lucy-Richardson算法同步挤压广义S变换进行地震数据分解。仿真实验证明，常规时频分析方法的时频谱能量不够集中，产生了一定的虚假频带宽度，本文所选的LR算法同步挤压广义S变换具有时频分辨率高、频率汇聚性强的优点，对区分信号中的不同频率成分效果显著。将该技术应用于川东南某研究区凉高山②号砂组的描述，提高了剖面上砂体边界的识别度，更好地刻画了河道平面展布特征，为下步勘探部署提供了有力依据。

文章引用

张懿疆,陈芊澍,张玉玺,袁 萌,周 远,江 馀. 基于LR算法同步挤压广义S变换的隐蔽河道识别
Hidden Channel Sand Bodies Identification Based on the Synchrosqueezing Generalized S-Transform with the Lucy-Richardson Algorithm[J]. 石油天然气学报, 2023, 45(02): 119-129. https://doi.org/10.12677/JOGT.2023.452015

参考文献