牛双晨 程冰洁

(①成都理工大学油气藏地质及开发工程国家重点实验室，四川成都 610059；②成都理工大学地球勘探与信息技术教育部重点实验室，四川成都 610059；③四川中成煤田物探工程院有限公司，四川成都 610072)

0 引言

众所周知，河道砂体往往是油气藏的重要载体，随着勘探开发不断深化，埋深更深、宽度更窄、厚度更薄的河道逐渐成为油气勘探领域的研究热点。针对河道空间展布特征识别的难题，业界广泛运用振幅、频率、相位等地震属性和速度、密度、阻抗等反演参数开展研究。然而，受河道埋藏深度大、纵横向分布复杂和地震资料主频偏低、频带窄、空间分辨率不足等因素的制约，传统方法对河道的识别精度往往不高，识别效果欠佳，不利于气藏空间分布预测、储量估算、水平井钻井等勘探开发和工程作业。

近几年来，时频分析技术在地震勘探领域发挥着越来越重要的作用[1]，被广泛应用于主振幅和主频率等参数计算[2]、地震与测井信号的多分辨率分析[3]、储层预测[4]和吸收补偿[5]等。目前，地震信号时频分析技术也成为河道识别的重要手段，应用越来越广泛，针对地下河道沉积相的空间延展、厚度和宽度变化等特征，基于地震波在河道介质中的传播特性，利用地震信号开展时频分析，在时间域和频率域同时研究振幅、功率谱、能量、相位等地震属性。如，陈杰等[6]分析了短时傅里叶变换(STFT)的适用条件，并应用于河道砂体识别；杨道庆等[7]通过离散傅里叶变换，利用分频技术识别河道砂体；冯金义[8]应用时频分析技术提取了不同频率的瞬时振幅沿层切片，以检测河道砂体的空间分布。

1946年，Gabor[9]最早提出了利用声波信号频谱图等时频分析技术。根据核函数的不同，把时频分析法分为线性、双线性和非线性三大类[10]。其中，线性方法包括Gabor频谱图、STFT[11]、连续小波变换(CWT)[12]、S变换(ST)[13]和广义S变换(GST)[14-16]等，该类方法计算速度较快，但多具有窗口依赖性，因而无法同时提高时间分辨率与频率分辨率。双线性方法又称为二次型时频分布，包括WVD[17]、平滑伪Wigner-Ville分布(SPWVD)[18]和Choi-Williams分布(CWD)[19]等，以WVD为例，本质属于信号自相关系数的傅里叶变换，不受Heisenberg测不准原理[20]的制约，具有理论上优良的数学性质[21]。非线性方法包括基于数据驱动的自适应分解方法，如自适应STFT[22-23]，相比STFT窗口长度灵活可变；Huang等[24]提出的经验模态分解(EMD)、集合经验模态分解(EEMD)以及完全集合经验模态分解(CEEMD)[25-26]，时频分辨率高，但仍存在受信噪比影响大和计算效率低等问题；融合CWT与EMD的同步挤压小波变换[27]，以及融合STFT、ST的同步挤压变换(SST)[28-29]，分辨率和计算效率都很高，缺点是对低信噪比信号分析鲁棒性不足，且严格意义上同步挤压技术只是对时频变换结果的谱图重排处理。此外，还有基于信号稀疏表示的非线性方法，如基于Gabor原子的匹配追踪分解[30]等。

与非线性方法等其他时频分析方法相比，WVD具有极高的时频聚焦性，然而WVD存在交叉项干扰的缺陷。目前，最常用的消除交叉项的方法是核函数平滑法，如伪Wigner-Ville分布(PWVD)、SPWVD以及CWD等，然而这些方法均以牺牲时频分辨率为代价，且对交叉项干扰抑制不彻底。有学者也提出了许多消除交叉项的新方法，刘希康等[31]提出联合STFT和WVD消除交叉项，结合STFT不存在交叉项干扰和WVD时频聚焦性高的优点，既抑制了交叉项干扰，又保留了高分辨率特性，但仍存在依赖窗函数和地震信噪比等问题。孔慧芳等[32]提出基于分数阶傅里叶变换(FRFT)抑制WVD的交叉项，在有效抑制交叉项干扰的同时具有一定的抗噪性能，但仅限于线性调频信号。随后，纪永祯等[33]提出了结合贝叶斯学习算法和WVD消除交叉项，该方法在分辨率和保真性能方面具有一定优势，但依赖雷克子波库。与上述消除交叉项的方法相比，Zoukaneri等[34]提出的WVD-MEM具有极高的时频分辨率，计算简便、抗噪性能好，不依赖雷克子波等各种子波库，特别是对交叉项干扰有良好的抑制效果，且适用于包括地震信号在内的非平稳随机信号。

本文在前人研究基础上，利用非线性调频(FM)信号、合成地震信号和实际地震数据研究WVD-MEM消除交叉项干扰的方法。在对比分析目前广泛应用的STFT、CWT、GST、SPWVD、CWD等时频分析方法的基础上，针对调频信号、仿真信号、单道实测地震信号和楔状正演模型，开展WVD-MEM薄层分辨率定量分析和抗噪性能测试，并成功应用于四川盆地侏罗系沙溪庙组气藏的河道识别，基于WVD-MEM方法实现了对实测地震信号功率谱、瞬时振幅等高分辨地震属性的提取，获取了河道空间展布、宽度、厚度等重要信息，为气藏勘探开发提供了方法支撑。

1 方法原理

1.1 Wigner-Ville分布(WVD)

Wigner-Ville给出的WVD定义为

(1)

式中：z(t)为信号x(t)的解析信号(解析信号包括原始信号和希尔伯特变换)；上标“*”表示共轭转置；t是时间；f是频率；τ是时延；j为虚数单位。

在时间域进行离散处理，信号x(n)、z(n)的关系被定义为

z(n)=x(n)+jH[x(n)]

(2)

式中：H[x(n)]为信号x(n)的希尔伯特变换；n为采样点整数，且n∈[0，N-1]；N为最大采样点数。

利用解析信号z(n)可以生成信号自相关系数协方差矩阵C

C=zzT

(3)

式中：C为N×N阶矩阵；上标“T”表示复共轭转置。

协方差矩阵C具有沿主对角线两端对称及复数实部相等、虚部相反的性质，求解C的一半即可完整求出矩阵。矩阵C沿每条交叉对角线的序列，构成了离散WVD的核函数，可以表示为

K(n)={kn(-l),…,kn(0),…,kn(l)}

(4)

式中：kn(0)表示矩阵C的主对角线元素；l∈[0，N-1]。定义kn(l)为

(5)

当l=0时，中心项kn(0)也是矩阵C的主对角线元素，有

kn(0)=z(n)z*(n)n=1，…，N

(6)

W(n)={Wn[-(N-1)],…,

Wn(0),…,Wn(N-1)}

(7)

1.2 基于最大熵(MEM)原理的功率谱计算方法

最大熵(Maximum Entropy Method，MEM)原理指出，信号熵值最大的概率分布是最合理的分布。这样，基于最大熵原理的功率谱求解与自回归(autoregressive,AR)模型功率谱的求解是等价的[35]，基本方程可以表示为

(8)

式中：P(f)为功率谱；p代表最大熵AR模型的阶数；m表示阶数p的序数，取值范围为[0，p-1]；cm是序数为m阶时的自回归系数，当m=0时c=1；Ep-1对应p-1阶时的预测误差能量；Δt表示采样间隔。

在其他变量已知的前提下，方程求解的关键在于自回归系数c和预测误差能量E的求取。利用式(6)可得到E0，利用Burg算法[35]可直接从z(n)中求出自回归系数c

(9)

c(m,i)=c(m-1,i)+c(m,m)c*(m-1,m-i)

i=1,…,m-1

(10)

Em=Em-1[1-c(m,m)c*(m,m)]

铝股和芯棒在稳定状态下的应力水平仿真结果如表9所示。由图11和表9可看出：导线铝股自外向内应力逐渐减小，导线内层铝股与芯棒接触位置应力最大。仿真结果证明导线外层所施加的径向载荷能够很好地传递给芯棒。

(11)

式中kz(m)为输入信号自相关函数。联合式(9)～式(11)，进入Levinson递归求解完整的自回归系数c和预测误差能量E，再由式(8)求得AR模型的最大熵功率谱。

1.3 基于WVD-MEM的高分辨时频属性计算方法

离散WVD核函数的累计求和，是输入信号自相关函数的偶数项序列，表示为

(12)

Burg在Levinson递归方程中推导出自相关函数的计算公式

c(m,m)Em-1

(13)

联合式(12)和式(13)，得到WVD-MEN的核函数

c(m,m)Em-1

(14)

与WVD核函数求取方法类似，根据协方差矩阵的共轭性质，只需求取右半部分的值即可求出完整的WVD-MEM的核函数。之后对其进行快速傅里叶变换，可求出WVD的最大熵瞬时功率谱，即WVD-MEM时频分布。由于时频分布的峰值对应瞬时振幅的平方，且峰值在时频面的投影为瞬时频率，因而利用峰值检测法[36]可以获得瞬时频率和瞬时振幅等属性参数。

WVD-MEM高分辨时频属性的主要计算步骤如下：

(1)由式(2)对原始输入数据x(n)做希尔伯特变换，获得解析信号z(n)。

(2)为了抑制数据变换结果的端点效应(即时频表示平面上两端数据随窗口增大而消失的现象)，需对z(n)进行加窗口扩展，在z(n)首尾各添一段数据，窗口长度均为L，设定值均为0。

(15)

式中：L为窗口长度，取奇数；zL(n)为解析信号z(n)扩展后的新数据，以z(n)为中心。

(3)设定最大熵AR模型阶数。

(4)由式(6)计算最大熵模型的预测误差初值(即WVD-MEM核函数的中心项)。

(5)利用式(9)～式(11)求自回归系数c。

(6)利用式(14)及自回归系数c，求WVD-MEM的核函数knWM(m)。

(7)对knWM(m)做快速傅里叶变换，获得WVD-MEM高分辨时频属性。

需要说明的是，由于第(2)步扩展数据时首尾加上了窗口长度L的数据，因此需对求得的WVD-MEM时频属性结果矩阵边缘长度为L的数据做舍去处理。

采用上述步骤不仅可以求取地震信号的功率谱，还可以计算瞬时频率、瞬时振幅等地震属性。由于WVD-MEM时频属性计算方法结合了WVD的时频高聚焦特性和MEM的最大概率分布两项优势，因此获得的功率谱、瞬时频率、瞬时振幅等地震属性具有显著的高分辨特征，可为薄储层预测、窄河道识别及其他微型地质体研究提供方法支撑。

2 方法试验与分辨率定量分析

2.1 调频信号定量试验

利用非线性调频(FM)信号研究WVD-MEM的时频分辨能力。图1a显示了正弦调频信号x1、双曲线调频信号x2和二者线性叠加信号x的时间域波形特征。由图可见，x1的振幅随着时延增大,曲线变密，x2的振幅随着时延减小，x的振幅随着时延变化并不明显，但曲线前半部分较为稀疏，后半段较密。图1b显示x1和x2在频率域的变化特征，其中，x1的频率呈高低起伏的“波浪状”，x2的频率随着时延逐渐降低，且降低的幅度越来越小。

图1 时间域和频率域的非线性调频信号

图2是分别采用STFT、CWT、GST、WVD、SPWVD、CWD和WVD-MEM等7种时频分析方法获得的叠加信号x的时频谱。图2a显示叠加信号x的实际时频谱是x1和x2两个分量的叠加，在时频域被分离成两部分，同时具有x1和x2的频谱特征。图2b～图2d为采用线性时频分析方法STFT、CWT和GST处理获得的时频谱，虽不存在交叉项干扰，对分量x2的频谱聚焦性也较好(黄色箭头所示)，但对分量x1的频谱聚焦性相对较差(红色箭头所示)。WVD方法虽然时频聚焦性好，但存在严重的交叉项干扰(图2e中绿色箭头所示)。SPWVD(图2f)和CWD(图2g)处理虽然消除了WVD产生的交叉项干扰，但降低了频率分辨率，尤其是CWD，存在时窗平滑效应产生的带状干扰(图2g白色箭头所示)。相比之下，WVD-MEM方法(图2h)不仅完整消除了WVD存在的交叉项干扰，而且保持了与实际频率基本一致的时频分布，清楚地显示出x1和x2频谱分别由低到高和由高向低的变化特征。

图2 采用不同方法获得的叠加信号x时频分布特征

总之，试验表明这7种方法均能分辨叠加信号x中的x1、x2分量，且线性时频分析方法中，GST的聚焦性优于STFT和CWT；二次型时频分析方法中SPWVD、CWD和WVD-MEM的聚焦性整体优于线性时频分析法。其中WVD-MEM方法具有更显著的高分辨率优势。

2.2 单道仿真地震信号——非平稳信号定量试验

利用单道仿真地震信号定量检测WVD-MEM方法对非平稳信号的高分辨时频分析能力，并为实际地震信号功率谱、瞬时振幅等高分辨时频属性的提取和应用奠定基础。图3a为使用5个雷克子波分量合成的单道仿真地震信号，子波分量的主频分别为：①20Hz；②30Hz和60Hz；③40Hz、30Hz和60Hz；④60Hz、40Hz和60Hz；⑤60Hz、60Hz和60Hz。最大振幅分别为：①30dB；②20dB和-10dB；③15dB、-10dB和20dB；④30dB、20dB和-30dB；⑤20dB、40dB和20dB。信号长度(采样总数)分别为①100；②100；③100；④100；⑤101。具体参数选取和计算方程如式(16)所示。

(16)

式中：r(f,t)为雷克子波表达式；f1～f5为5个雷克子波分量频率。

图3b～图3h分别为采用STFT、CWT、GST、WVD、SPWVD、CWD和WVD-MEM等7种方法获得的仿真信号时频分布。由图可见，7种方法均在0.1、0.3、0.5、0.7及0.9s附近检测到相应的时频能量响应，但能量分布具有明显的差异。三种线性时频分析方法中，GST聚焦性和分辨率最好，CWT次之，STFT最差；相比之下，二次型时频分布的聚焦性均优于GST，但WVD法在0.2、0.4、0.6和0.8s附近出现了严重的干扰项(图3e中绿色箭头所示)；SPWVD和CWD虽然消除了图3e中的交叉项干扰，但对于0.5、0.7和0.9s附近的交叉项干扰却无法完全抑制(黄色箭头所示)。而WVD-MEM则完全消除了WVD中存在的所有交叉项干扰，表现出最优的时频聚焦性和频率定位精度，清楚地显示了20、30、40、60Hz等频率位置的功率谱强弱变化。例如，在0.9s附近的仿真信号由3个60Hz主频、中间振幅较高(40dB)、两侧振幅较低且相同(20dB)的雷克子波合成，WVD-MEM时频分析方法不仅准确展示出子波能量变化特征，而且精确定位在60Hz频率位置，其他方法均未能达到如此高的分辨率。

图3 单道仿真地震信号不同方法时频分析效果对比

2.3 单道实测地震信号——非平稳信号试验

上述仿真信号是由与实际地震信号近似的雷克子波组成，已知子波频率、相位、振幅等信息，而实测地震信号中这些信息是未知的，且更加复杂，二者存在明显差异。因此，有必要针对实测地震信号开展试验，进一步分析该方法对实际地震信号的适用性，为后续利用功率谱、瞬时频率、瞬时振幅等地震属性进行高分辨河道识别提供参考。图4为单道实测地震信号及采用7种方法的时频分析结果对比。同样，在STFT、CWT和GST等3种线性时频分析方法中，GST的时频分辨率和聚焦性最优；WVD方法受到交叉项的干扰，有效的时频特征被完全淹没，无法进行信号分析；SPWVD和CWD方法的时频聚焦性优于线性时频分析方法，但交叉项干扰抑制不彻底(黄色箭头所指)。WVD-MEM则彻底消除了WVD存在的严重交叉项干扰，在保持高度时频聚焦性的同时，具有明显更高的时频分辨特性，对能量随频率和时间变化的细节刻画得更精确。可见，针对更复杂的实际地震信号，WVD-MEM方法仍然能获得高聚焦、高时频分辨率的处理效果。

图4 单道实测地震信号时频分析效果对比

2.4 多道正演数值模拟数据的WVD-MEM时频分辨能力定量分析

河道识别本质上是地震极限分辨率条件下的薄层识别问题。通过提高瞬时功率谱、瞬时频率、瞬时振幅等地震属性的分辨率，可以提高地震信号的薄层分辨率，从而实现河道识别。设计一个楔状模型正演数值模拟地震数据(图5a)，震源为主频35Hz的雷克子波，检波器的采样间距为5m，接收长度为1.5s，共201道。在地层变薄位置(红色与蓝色圈内)，可以观察到薄层调谐效应产生的明显复合地震反射现象。

图5b～图5d分别显示了利用GST、SPWVD和WVD-MEM三种方法提取的35Hz单频剖面。由图可见，在红色圈标识的薄层调谐位置，GST、SPWVD及WVD-MEM方法的最大识别位置分别为第66、54、41道；在蓝色圈标识的薄层调谐位置，GST、SPWVD及WVD-MEM方法的最大识别位置分别为第132、144、161道。显然，WVD-MEM表现出更高的能谱聚焦性，频谱被高度细化，增强了时频分辨能力，能更好地消除薄层调谐效应的影响，更有利于薄层识别。

图5 正演数值模拟记录及采用不同方法获得的35Hz单频剖面

2.5 WVD-MEM的抗噪性能测试分析

实际地震信号通常含有一定程度的噪声，因而，有必要在不同信噪比(SNR)条件下测试WVD-MEM方法的抗噪性能，以进一步明确其实用价值。图6显示单道实测地震信号(图6a)及SNR分别为30、20和10噪声条件下，应用WVD-MEM获得的瞬时振幅谱特征。可见，随着SNR由高向低转变，相应的WVD-MEM处理效果与原始地震信号处理结果相比较并未发生明显改变。由此可以证实，WVD-MEM具有良好的抗噪性能，能够适用于实测地震信号的时频处理。

图6 不同SNR实测地震信号的WVD-MEM处理结果对比

综上所述，在非线性调频信号、合成地震信号、单道实测地震信号及多波正演数值模拟地震记录的定量测试分析中，WVD-MEM时频分析方法都表现出较强的能谱聚焦和谱细化特性，时频分辨能力显著，且抗噪性能良好。因此，当地震波在地下介质中传播，遇到目标地质体时可能出现功率谱、瞬时振幅、瞬时频率等时频属性异常，借助WVD-MEM方法可以更精确地突出时频属性异常特征，进而有效识别薄储层、窄河道及其他微型地质体。

3 应用实例

图7 过JS205井(左)和JS301井(右)河道的地震反射特征

地震数据频谱分析对提取瞬时振幅及分频属性具有指导作用。首先对研究区地震数据进行频谱分析(图8)，可以看出，工区地震数据主频约为35Hz，频带范围为5～80Hz，优势频带范围为7～65Hz。结合前文方法与试验研究，采用WVD-MEM方法处理工区地震数据并提取瞬时频率、瞬时振幅和不同频率功率谱等地震属性，从中发现时频属性异常，以多视角分析出河道的异常响应特征，实现河道空间展布特征识别。

图8 研究区频谱分析

图9 研究区WVD-MEM提取的不同频段频谱剖面

图10 利用WVD-MEM方法提取的时频属性识别层河道分布

采用分频融合技术对WVD-MEM计算的低、中、高三个频段的功率谱特征开展多频融合增强显示，进一步增强河道识别效果。图11a为RGB融合后识别效果，图11b为CMY融合后的河道识别效果。由图可见，经过多频融合后，不同尺度、不同流向的河道得到进一步清晰展示，对一些较宽河道刻画效果更加明显(红色箭头所示)，对于原始振幅切片上无法识别的两条不明显窄河道(绿色箭头所示)刻画效果明显改善。因此，利用WVD-MEM方法提取低、中、高三个频段的功率谱属性并融合显示，能够突出共性、弱化差异，有效改善单频段频谱在河道识别方面的不足，增强对工区河道的高分辨识别效果。

图11 研究区WVD-MEM不同频率功率谱融合层位切片特征

这里需要进一步指出的是，基于时频分析方法进行河道识别，是将地震数据从时间域转换到时频域，即将地震数据从一个时间采样点对应一个原始振幅值转换到一个时间采样点对应不同频率下的多个值(无论振幅谱值还是功率谱值)，本质上不同于从时间域的角度分析地震信号，属于在时频域对原始地震信号的直接应用，不会破坏原始振幅的连续性，因而识别结果具有可靠性。而影响河道高分辨识别效果的因素取决于时频分析方法的能谱聚焦性和分辨率。由于WVD-MEM方法结合了Wigner-Ville分布和最大熵原理的非平稳信号分析优势，不仅能够有效消除WVD信号分析中存在的交叉项，而且可以保持与瞬时频率基本一致的高度时频聚焦性，使计算的瞬时振幅、功率谱等地震属性具有高分辨率特性。

该方法在川西地区侏罗系沙溪庙组气藏河道识别的应用实例表明，WVD-MEM提取的功率谱、瞬时振幅等地震属性，能从垂向和横向较好地刻画河道的沉积厚度、延展宽度和流动方向等地质细节，可以在其他类似的油气探区开展推广应用。

4 结论与认识

通过WVD-MEM方法研究和对非线性调频信号、仿真地震信号、单道实测地震(非平稳)信号、楔状模型正演模拟信号等定量试验，以及在川西侏罗系沙溪庙组气藏河道识别的实际应用，得到以下认识：

(1)相比于STFT、CWT、GST、WVD、SPWVD、CWD等时频分析方法，WVD-MEM时频分析方法结合了WVD的高度时频聚焦性以及MEM的最大概率分布优势，不仅完整消除了WVD中存在的严重交叉项，且具有显著的频谱聚焦性和高分辨率优势，可以应用于薄储层、窄河道及其他地下微型地质体的识别研究；

(2)针对埋藏较深、宽度较窄、厚度较薄、流向变化的河道高分辨识别难题，WVD-MEM方法能够在一定程度上解决地震资料主频不够高、频带不够宽等不利因素导致的空间分辨率不足的问题，利用WVD-MEM获得的瞬时振幅、分频功率谱等高分辨地震属性以及多频融合等手段，能够有效突出河道的厚度、宽度、流向等空间展布细节，进一步提高河道的识别精度；

(3)在对WVD-MEM实测地震信号分析过程中，就整个工区而言，阶数p一般选择定值。如何在每道实测地震信号计算前自适应调整p的取值、优化WVD-MEM计算效率和程序设计以及方法的应用推广，将是下一步研究的方向。

感谢中国石化西南油气分公司为此研究提供实验数据。