一种基于稀疏窗S变换的分频-重构波阻抗反演方法*

冉喜阳 周怀来 张益明 李雷豪 杨吉鑫

(1.成都理工大学地球物理学院 四川成都 610059； 2.成都理工大学油气藏地质及开发工程国家重点实验室 四川成都 610059；3.四川省地震局 四川成都 610041； 4.中海油研究总院有限责任公司 北京 100028)

在地震勘探中，不同频率的地震反射波特征反映了不同的地层岩性、物性及含流体性质[1],故相比常规地震资料，分频技术的应用能够从地震数据中获取更多的有效信息，有利于对目标体进行更全面的分析解释，进而提高储层预测精度,降低解释多解性。地震分频技术的基础是时频分析方法，目前常见的时频分析方法包括短时傅里叶变换、Gabor变换、Wigner Ville时频分布、小波变换、S变换等[2]。

为了从地震数据中获得更多的有效信息，研究人员对分频技术做了很多的研究工作[3-12]。前人利用各种时频分析方法对测井数据和地震数据进行分频处理，发现不同频率尺度地层特征在此基础上进行波阻抗反演获得了较高分辨率的反演结果，但在应用过程中也存在一些问题。王振卿 等[5]对塔中西部台内数据用广义S变换进行谱分解，使用带通滤波处理得到15～25 Hz优势频段，对其进行波阻抗反演，虽然提高了识别精度，应用效果良好，但是在滤波过程中也滤掉了大量的有效波，忽视了部分细节。

另外，在常规的时频分析方法中，由于短时Fourier变换的窗函数固定，不能适应非平稳信号中不同频率成分特性的分析；S变换通过改进时间窗函数的尺度因子，具备多尺度分析能力，但S变换中基本小波是固定的，使其在应用中受到了限制；高静怀 等[12]在S变换的基础上，增加4个可调属性参数得到广义S变换，有效地增大了其应用范围，方法效果尤为依赖对窗参数的调节。

稀疏窗S变换在保证高分辨率的同时还可以降低时间成本[13]。本文在前人研究基础上，提出一种基于稀疏窗S变换的分频-重构波阻抗反演方法，将地震数据分割成N个不同频带的数据体，并利用分频井曲线约束进行反演，以提高对薄层的识别精度。

1 基本原理

1.1 稀疏窗S变换原理

稀疏窗S变换基于S变换，通过对窗口参数进行稀疏优化，实现自适应地调整非平稳信号的频率突变，具有与非优化算法相同的计算复杂度。传统滤波方法减少了频谱带宽，虽然使道分析结果变得稳定，却无法获得多个薄层的高分辨图像。稀疏窗S变换作为一个强大的时频分析方法能实现高分辨率复杂道分析，可以解决上述冲突。S变换的频率域表达式为

(1)

用稀疏概念来解决这个问题，就是通过定义窗函数来控制它们的“影响区域”，在“影响区域”内定位强振幅频率分量，并在整个时频域内混叠处理低振幅分量。为了以有效的2D方式解决这个问题，Hamid Sattari[13]将式(1)以矩阵形式重写为如下自适应形式：

(2)

基于振幅谱构建As[l][l,m]以区分不同频率分量，需要满足以下2个前提条件：

一是单位分解化，也就是其所有列加起来等于一个列单位向量，其矩阵形式为

(3)

为了使时频分解可逆，要将其进行归一化处理。而单位分解在窗参数优化中不仅仅是强制性的前提条件，同时也将稀疏性引入了该方法。

二是窗函数对称性，表达式为

a(s[l],m)=a(s[-l],m)

(4)

该条件保持了S变换分析窗正负频率的对称性。

对于s[l]的求取，由于需要根据相应的振幅谱作调整以获得自适应窗，在其求取过程中应用了输入信号频谱平滑后的结果。Hamid Sattari[11]构建的窗参数优化公式如下：

(5)

(6)

(7)

(8)

最后将得到a(s[l],m)分别置于As[l][l,m]的m列上即可得到自适应窗函数矩阵，再根据式(2)即可得到稀疏自适应S变换的时频谱。

1.2 基于稀疏窗S变换的分频原理

基于稀疏窗S变换的时频滤波器可表示为

(9)

基于稀疏窗S变换的分频,实际上是将通过稀疏窗S变换获得的时频谱与时频滤波器相乘的过程。利用原始信号的频谱图加以控制，通过分频处理使有用的信息被保留，不必要的信息被去除,获得地震道时频谱并进行滤波后得到所需的频率分量。然而，这些数据仍然只是时频数据，还需进行反变换才能获得最终时间域的分频结果[14]。使用反变换为

(10)

式(10)中：hf(t)是分频后的地震信号；FFT-1是傅里叶反变换。先将二维的时频谱乘以频率域窗函数，接着对时间轴积分，再进行傅里叶逆变换，最终获得了时域的分频结果。

1.3 基于稀疏窗S变换的分频-重构波阻抗反演流程

基于稀疏窗S变换的分频-重构波阻抗反演流程包含叠后地震资料频谱分析及分频处理、测井资料预处理、波阻抗反演、数据融合等主要环节，最终获得重构波阻抗反演数据体。利用上述分频方法，设计了基于稀疏窗S变换的分频-重构反演技术流程，如图1所示。

图1 分频-重构波阻抗反演流程图Fig.1 Frequency-divided impedance inversion flow chart

首先，对原始叠后地震数据进行频谱分析，再根据其频谱图选取适当的频率范围将原始地震数据分为3段，通过分频处理得到低、中、高频地震剖面。然后，利用分频井曲线约束分别对3个分频地震数据进行波阻抗反演，之后将反演得到的各数据体再次进行频谱分析，取3个反演结果中相对应的频率部分进行数据重构，最终得到重构波阻抗数据体。

2 模型试算

为了验证稀疏窗S变换的优越性，本文设计了2个模拟实验。图2为利用短时傅里叶变换、S变换、稀疏窗S变换3种时频分析方法分别对Chirp信号进行的频谱分析图，可以看到短时傅里叶变换和S变换由于窗参数的固定，对不同频率成分的能量聚集性较弱，随着频率变化，对于非平稳信号的时频刻画能力有所降低；而稀疏窗函数由于对窗参数进行了稀疏的优化，对强振幅的频率分量作了相应的窗函数调节，得到的时频谱分辨率较高、能量聚集性较强，且能很好的适应该信号。

图2 Chirp信号仿真模拟Fig.2 Chirp signal simulation

图3为合成记录仿真模拟图，其中第1列为设计的反射系数，第2列为反射系数与45 Hz雷克子波褶积的合成记录。由图3可以看出，稀疏窗S变换在时间分辨率变小的情况下依然能较准确的分离信号，而短时傅里叶变换在时间间隔较小的情况下信号的时频谱出现耦合，无法分辨两个相邻信号，对识别造成了影响。这表明，稀疏窗S变换确实能有效提高时频谱的分辨率。

图3 合成记录仿真模拟Fig.3 Synthetic record simulation

3 方法验证

不同频率的测井曲线对储层的敏感性不同，大致来说，低频曲线反映宏观趋势，高频曲线反映薄层细节。为了验证该分频方法的可行性，本文选取含气B井中的P-wave曲线进行试处理，图4为曲线时频谱，图5为曲线分解图。

图4精确地反映了P-wave曲线的时频分布规律，据此图可将井曲线分成低、中、高频3段；从图5中原始井曲线与3条曲线的对比可以看出，声波曲线对地下岩性变化较为敏感；低频曲线反映尺度较大，主要刻画地层组层速度变化；高频曲线反映尺度较小，重点突出薄层内物性变化。3条分频曲线与岩性剖面吻合较好，说明该分频方法具有可行性。

利用上述分频井曲线约束井旁地震道做单道分频反演，以验证该方法的正确性(图6)。图6中，第一列曲线为原始P-wave曲线，曲线1为利用本文方法得到的分频反演曲线，曲线2为常规反演曲线。从原始井曲线与2条反演曲线的对比可以看出，2条反演曲线的大体趋势都能较好地吻合原始井曲线。由于分频反演较为充分的利用了地震数据的高频成分，故在细节反映上分频反演得出的结果明显优于常规反演方法得出的结果，说明该方法的正确性以及优越性。

图4 B井P-wave稀疏窗S变换时频谱Fig.4 P-wave Sattari time spectrum of Well B

图5 B井P-wave曲线分解图Fig.5 P-wave curve decomposition diagram of Well B

图6 B井旁道分频反演图Fig.6 Frequency-divided inversion of seismic traces of Well B

4 应用实例

为了验证分频-重构波阻抗反演的实用性，选取中国东南部L工区过A—F井的连井地震剖面进行分频-重构波阻抗反演测试。图7为原始地震剖面及其频谱分析图，可以看出该地区目的储层较发育，但地震资料分辨率较低制约了砂体接触关系的确定，储层连通性存在不确定性，局部砂体较薄，常规波阻抗反演方法受到了限制。根据分频-重构波阻抗反演方法进行反演，图8～10分别为分频后得到的低、中、高频地震剖面，图11为传统波阻抗反演剖面，图12为分频-重构波阻抗反演剖面。由图11与图12对比可见，2种反演方法均能识别出有效的气层，但由于常规反演方法更多受主频控制，得到的反演结果分辨率较低，难以区分气层和气-水过渡带；而分频-重构反演合理地利用了地震数据中相对低频与相对高频的成分，使得反演结果与井吻合度更高，有效地提高了分辨率。此外，区域中C井附近断层较发育，但传统波阻抗反演难以获得有效相关信息，对断层反映极其不敏感，而分频-重构反演的剖面可以在目标井附近显示出明显的地层错断，说明其对断层的识别更有效；并且传统的波阻抗反演难以反映出地层的薄厚变化，而分频重构波阻抗反演却能更细腻准确地描述地层。以上3个优势显示了该反演方法有较好的应用潜力。

图7 中国东南部L工区原始地震剖面及其频谱分析图Fig.7 Original seismic profile and its spectrum analysis in L work area,Southeast China

图8 中国东南部L工区低频地震剖面及其频谱分析图Fig.8 Low-frequency seismic profile and its spectrum analysis in L work area,Southeast China

图9 中国东南部L工区中频地震剖面及其频谱分析图Fig.9 Mid-frequency seismic profile and its spectrum analysis in L work area,Southeast China

图10 中国东南部L工区高频地震剖面及其频谱分析图Fig.10 High-frequency seismic profile and its spectrum analysis in L work area,Southeast China

图11 中国东南部L工区常规波阻抗反演剖面Fig.11 Conventional impedance inversion profile in L work area,Southeast China

图12 中国东南部L工区分频-重构波阻抗反演剖面Fig.12 Frequency-divided impedance inversion profile in L work area,Southeast China

5 结论

1) 稀疏窗S变换对窗参数进行了基于振幅谱的稀疏优化，令窗参数根据不同信号的振幅谱作不同适应性窗参数优化，提升了运算速度，该时频分析方法具有分辨率高和能量聚集性好的优势。分频-重构波阻抗反演可以充分利用地震资料中的低频和高频信息，获得的最终剖面具有与井吻合度高、断层反应清晰、地层刻画准确等特点。

2) 本文方法的关键之处在于如何调节参数控制分频过程和效果以及如何将低、中、高频波阻抗数据体进行重构，这都将直接影响到最终反演结果。重构结果高频成分过重，会产生一些假象；重构结果低频成分过重，将影响反演的分辨率。如何更精确控制分频反演中各频率成分占比，使反演结果更合理，也是作者之后的研究内容。