针对窄细古河道精细刻画难题,利用最大熵准则增强Wigner-Ville分布的聚焦特性,在有效提升地震信号时频分辨能力的基础上,建立了一种微型古河道识别新方法。基于最大熵功率谱与自回归模型（AR）功率谱等效的原理,首先利用Burg算法和Levinson-Durbin递推规则,求取AR模型的预测误差、自回归系数等参数;然后,在自相关函数一阶导数为0的条件下,计算地震信号的Wigner-Ville分布,获取微型古河道最大熵准则约束下的Wigner-Ville时频功率谱（MEWVD）。通过仿真地震信号和窄薄模型数值模拟信号实验分析,发现MEWVD既能有效避免Wigner-Ville分布的交叉项干扰,还能获得比短时傅里叶变换（STFT）、连续小波变换（CWT）等信号分析方法更加精准的频谱特征;同时,还证实了利用不同频率的MEWVD,可以有效识别不同尺度的窄细古河道。将该方法应用于四川盆地中江气田侏罗系沙溪庙组（J2s33-2小层）气藏,准确地识别出宽度小于500 m、砂岩厚度小于35 m的窄细古河道的宽度、走向等空间信息,可为井位部署、水平井压裂选段等提供依据。图7表2参31
In view of the problem of fine characterization of narrow and thin channels, the maximum entropy criterion is used to enhance the focusing characteristics of Wigner-Ville Distribution. On the basis of effectively improving the time-frequency resolution of seismic signal, a new method of microscopic ancient river channel identification is established. Based on the principle of the equivalence between the maximum entropy power spectrum and the AR model power spectrum, the prediction error and the autoregression coefficient of AR model are obtained by using the Burg algorithm and Levinson-Durbin recurrence rule. Under the condition of the first derivative of autocorrelation function being 0, the Wigner-Ville Distribution of seismic signal is calculated, and the Wigner-Ville Distribution time-frequency power spectrum (MEWVD) is obtained under the maximum entropy criterion of the microscopic ancient river channel. Through analysis of emulational seismic signal and numerical simulation signal of narrow thin model, it is found that MEWVD can effectively avoid the interference of cross term of Wigner-Ville Distribution, and obtain more accurate spectral characteristics than STFT and CWT signal analysis methods. It is also proved that the narrow and thin river channels of different scales can be identified effectively by using MEWVD of different frequencies. The method is applied to the third member of Jurassic Shaximiao Formation (J2s33-2) gas reservoir of the Zhongjiang gas field in Sichuan Basin. The spatial information of width and direction of narrow and thin river channel with width less than 500 m and sandstone thickness less than 35 m is accurately identified, providing basis for well deployment and horizontal well fracturing section selection.
本文基于信号最大熵准则（Maximum Entropy Principle）, 在克服WVD交叉项干扰的前提下, 研究增强WVD频谱聚焦特性和提升地震信号时频分辨率的方法, 以突出窄细古河道响应特征。并采用仿真信号、正演数值模拟等理论数据开展试验与论证分析, 以形成微型古河道识别新方法, 为井位部署、水平钻井等提供支撑。
$H(x)=\int\limits_{-\infty }^{\infty }{p(x)\ln p(x)\text{d}x}$ （1）
此时, Burg定义信号$X$的功率谱熵为：
$H[P(f)]=\frac{1}{2\pi }\int\limits_{-\pi }^{\pi }{\ln P(f)df}$ （2）
$P(f)=\frac{{{E}_{M-1}}\Delta t}{{{\left| 1+\sum\limits_{j=0}^{M-1}{\alpha (j){{\text{e}}^{-2\text{i}\pi fj\Delta t}}} \right|}^{2}}}$ （3）
其中 $0\le f\le \frac{1}{2\Delta t}$
$W(t, f)=\int\limits_{-\infty }^{\infty }{z\left( t+\frac{\tau }{2} \right)}z'\left( t-\frac{\tau }{2} \right){{\text{e}}^{-2i\pi f\tau }}d\tau $ （4）
然而, 实测的信号并非连续信号, 实际应用时, 需要对$W(t, f)$进行离散化。设$x(t)$、$z(t)$的离散表达式分别为$x(n)$和$z(n)$（$0\le n\le N-1$）, 对$x(n)$作希尔伯特（Hilbert）变换, 在计算出$z(n)$后即可获得$W(t, f)$的离散表述：
${{W}_{\text{z}}}(n, f)=2\sum\limits_{l=-N}^{N}{z(n+l)z'(n-l)}{{\text{e}}^{-4i\pi lf}}$ （5）
其中 $0\le |l|\le N$
${{k}_{n}}(l)=\left\{ \begin{matrix} z(n-l)z'(n+l)\ \ \ \ \ \ \ |l|\le \min \left( n, N-n \right) \\ 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ |l|\min \left( n, N-n \right) \\ \end{matrix} \right.$. （6）
对于每一个样点n, ${{k}_{n}}(l)$将组成一个卷积核序列${{k}_{\text{z}}}(n)$, 如（7）式。当$n=0$和$n=N-1$时, 卷积核序列${{k}_{\text{z}}}(n)={{k}_{\text{z}}}(0)$和${{k}_{\text{z}}}(n)={{k}_{\text{z}}}(N-1)$最短; 当$n=\frac{N}{2}$时, ${{k}_{\text{z}}}(n)={{k}_{z}}\left( \frac{N}{2} \right)$最长。以${{k}_{n}}(0)$为中心, 对${{k}_{\text{z}}}(n)$左右补充0值, 可将${{k}_{\text{z}}}(n)$扩充为长度为$N$的卷积核序列。
$\left[ \begin{align} & {{k}_{\text{z}}}(0)=\{z(0)z'(0)\}=\{{{k}_{0}}(0)\} \\ & {{k}_{\text{z}}}(1)=\{z(2)z'(0), z(1)z'(1), z(0)z'(2)\}=\{{{k}_{1}}(-1), {{k}_{1}}(0), {{k}_{1}}(1)\} \\ & {{k}_{\text{z}}}(2)=\{z(4)z'(0), z(3)z'(1), z(2)z'(2), z(1)z'(3), z(0)z'(4)\} \\ & \ \ \ \ \ \ \ \ =\{{{k}_{2}}(-2), {{k}_{2}}(-1), {{k}_{2}}(0), {{k}_{2}}(1), {{k}_{2}}(2)\} \\ & \vdots \\ & {{k}_{\text{z}}}(N-1)=\{z(N-1)z'(n-1)\} \\ \end{align} \right.$ （7）
对第$n$个采样点的卷积核序列${{k}_{\text{z}}}(n)$, 进行离散傅立叶变换, 即可获得WVD瞬时功率谱, 如下：
${{W}_{\text{z}}}(n, f)=\left\{ {{w}_{n}}\left( -\frac{N-1}{2}, f \right), \cdots , {{w}_{n}}(0, f), \cdots , {{w}_{n}}\left( \frac{N-1}{2}, f \right) \right\}$ （8）
其中, ${{W}_{\text{z}}}(n, f)$的子项可表述为：
${{w}_{n}}(m, f)=\frac{1}{N}\sum\limits_{l=-\frac{N-1}{2}}^{\frac{N-1}{2}}{{{k}_{n}}(l){{\text{e}}^{-2i\pi \frac{ml}{N}}}}$ （9）
通过计算每一个采样点n对应的瞬时功率谱, 最终可以得到离散信号$x(n)$的Wigner-Ville功率谱, 即：
$P(f)=\sum\limits_{n=0}^{N-1}{{{W}_{z}}(n, f)}$ （10）
$\left\{ \begin{align} & {{E}_{\text{F}, m}}(n)={{E}_{\text{F}, m-1}}(n)+{{\alpha }_{m}}(m){{E}_{\text{B}, m-1}}(n-1) \\ & {{E}_{\text{B}, m}}(n)={{E}_{\text{B}, m-1}}(n-1)+{{\alpha }_{m}}(m){{E}_{\text{F}, m-1}}(n) \\ & {{\alpha }_{m}}(m)=\frac{-2\sum\limits_{n=m}^{N-1}{\left| {{E}_{\text{F}, m-1}}(n) \right|{{E}_{\text{B}, m-1}}(n-1)}}{\sum\limits_{n=m}^{N-1}{{{\left| {{E}_{\text{F}, m-1}}(n) \right|}^{2}}}+\sum\limits_{n=m}^{N-1}{{{\left| {{E}_{\text{B}, m-1}}(n-1) \right|}^{2}}}} \end{align} \right.$ （11）
通过离散信号$x(n)$及其解析信号$z(n)$, 可以计算预测误差的初始功率, 即
$P_{0}^{2}(f)=\frac{1}{N}\sum\limits_{n=0}^{N-1}{{{\left| z(n)z'(n) \right|}^{2}}}$ （12）
利用$m$阶前向预测误差和后向预测误差, 可以计算预测误差的平均功率, 即
$P_{m}^{2}(f)=\frac{1}{2}\sum\limits_{n=m}^{N-1}{\left\{ {{\left| {{E}_{\text{F}, m}}(n) \right|}^{2}}+{{\left| {{E}_{\text{B}, m}}(n) \right|}^{2}} \right\}}$ （13）
据Levinson-Durbin递推约束规则, 设${{\alpha }_{m}}(j)$为$M$阶AR模型在阶次序号为$m$时的第j个自回归系数, ${{E}_{m}}$为$m$阶次时的最小预测误差, 则其递推关系为
$\left\{ \begin{align} & {{\alpha }_{m}}(j)={{\alpha }_{m-1}}(j)+{{\alpha }_{m}}(m){{\alpha }_{m-1}}(m-j)\ \ \ \ (j=1, 2, \cdots , m-1) \\ & {{E}_{m}}=[1-{{\alpha }_{m}}{{(m)}^{2}}]{{E}_{m-1}}\ \ \ \ (m=1, 2, \cdots , M) \\ \end{align} \right.$ （14）
这样, 将${{\alpha }_{m}}(1)$, ${{\alpha }_{m}}(2)$, …, ${{\alpha }_{m}}(m-1)$和${{E}_{m}}$代入（3）式, 就能计算出WVD最大熵功率谱。WVD最大熵功率谱在时间域和频率域的分布, 即为Wigner-Ville最大熵时频谱（简称MEWVD）。
在时间域, MEWVD并未将瞬时自相关函数前后数据相乘, 而仅使用当前$n$点附近的数据, 避免了交叉项干扰; 在频率域, 受最大熵准则约束, 使每个$n$点对应的功率谱更精确, 能量更聚焦。同时, 在最大熵的前提下, 采用Levinson-Durbin约束递推未知瞬时自相关函数, 使计算的功率谱误差减小, 提升了信号时频分辨能力。
地震波在地下传播时, 受地层埋深、孔隙度、填充流体等因素影响, 将产生能量衰减, 引起动力学、运动学和几何学等特征发生变化。采用振幅、吸收系数、频率、相位和走时等参数, 可以仿真描述地震波传播过程, 如下：
$\left\{ \begin{align} & {{x}_{\text{F}}}(t)={{x}_{1}}(t)+{{x}_{2}}(t) \\ & {{x}_{1}}(t)={{A}_{1}}{{e}^{-Qt}}\sin (2\pi {{f}_{1}}t+{{\varphi }_{1}}) \\ & {{x}_{2}}(t)={{A}_{2}}\sin (2\pi {{f}_{2}}t+{{\varphi }_{2}}) \end{align} \right.$ （15）
上式中, 仿真地震信号${{x}_{\text{F}}}(t)$分别由衰减分量${{x}_{1}}(t)$和平稳分量${{x}_{2}}(t)$叠加而成。其中, ${{x}_{1}}(t)$包含了吸收系数、振幅、频率、相位和走时等参数, ${{x}_{2}}(t)$不包含吸收系数, 传播能量将不会产生衰减。
为了对比分析地震波传播过程中的时频特征, 利用${{x}_{\text{F}}}(t)$的参数, 构造出了关联特征较强的高、低频能量强弱不一样的仿真地震信号。如表1, 信号${{x}_{\text{F}1}}$和${{x}_{\text{F}2}}$的起始能量${{A}_{1}}$和传播时间t相同, ${{x}_{\text{F}1}}$中${{A}_{2}}$、${{f}_{2}}$、${{\varphi }_{2}}$为0, 表示${{x}_{\text{F}1}}$仅由衰减分量构成, 且其${{A}_{1}}$、${{f}_{1}}$和$Q$均较大, 类似强能量、高频和强衰减的地震信号。${{x}_{\text{F}2}}$包含了${{x}_{\text{F}1}}$, 仅通过改变${{f}_{2}}$, 在${{x}_{\text{F}2}}$的基础上增加了低频、弱振幅的平稳分量, 意味着${{x}_{\text{F}2}}$在t时间内的传播能量、振幅等将随频率变化而改变。如图1a和图1b中第1列, 分别显示了${{x}_{\text{F}1}}$和${{x}_{\text{F}2}}$的波形特征。可见, ${{x}_{\text{F}1}}$随着t的变化, 振幅减弱, 能量衰减很快; ${{x}_{\text{F}2}}$随着t的变化, 振幅减弱, 能量也在衰减, 但衰减分量的控制作用越来越弱, 在信号传播约0.25 s之后, 信号主要受平稳分量控制, 且由于${{f}_{2}}$差异, 使${{x}_{\text{F}1}}$振荡周期短、波形“ 瘦” , ${{x}_{\text{F}2}}$振荡周期长、波形“ 胖” 。
同时, 对比不同的时频分析方法, 发现MEWVD具有非常优良的时频聚焦特性和分辨率。如图1a和图1b, 第2列显示了仿真地震信号的快速傅里叶变换（FFT）频谱, 揭示信号${{x}_{\text{F}1}}$的能量峰值聚集在80 Hz, ${{x}_{\text{F}2}}$的能量峰值聚集在10 Hz和80 Hz。第3至10列, 分别显示了STFT、GST、CWT、WVD、PWVD、SPWVD、CWD和MEWVD等方法的时频效果。对比分析, 尽管${{x}_{\text{F}1}}$和${{x}_{\text{F}2}}$的最强能量均聚集在80 Hz高频位置, 但是, 不同的时频分析方法显示的能量聚焦程度与频率分辨能力存在差异, 且MEWVD聚焦效果与分辨率最优、WVD次之、STFT最弱。GST和CWT的时频分辨能力相似, 对80 Hz高频强振幅与10 Hz低频弱振幅均能有效识别, 但能量聚焦程度、时间分辨能力等较MEWVD差。此外, WVD时频特征还存在交叉项干扰, PWVD虽在一定程度上压制了交叉项干扰, 但能量聚焦性减弱, 降低了${{x}_{\text{F}2}}$的低频、弱振幅信息的分辨能力; SPWVD和CWD基本压制了交叉项干扰, 但也几乎无法有效识别${{x}_{\text{F}2}}$中低频、弱振幅响应, 较PWVD的分辨率损失更严重。
为了进一步证实Wigner-Ville最大熵时频谱（MEWVD）的频谱聚焦性和时频高分辨能力, 采用表2所示参数, 设计出二维数值模型, 包含了埋深、宽度和厚度各异的窄薄地质体, 其“ 透镜状” 的地质特征与微型古河道的横截面类似。基于二维地震波动方程正演数值模拟方法, 以35 Hz的Rick子波作为点震源激发, 以10 m采样间距、1 ms采样率接收, 获得了正演模拟记录。模型④— ⑨的砂体厚度逐渐变薄、宽度逐渐变窄; 模型④— ⑥较宽、厚度大于1/4波长（约30 m）, 顶底界面的反射同相轴明显, 而⑦— ⑨较窄、厚度小于1/4波长, 地震记录表现为谐波反射, 且⑦和⑧表现出强调谐反射, 极薄、极窄的⑨则表现为弱调谐反射（见图2）。
分别抽取代表模型④— ⑨的单道地震信号, 计算MEWVD后, 可分析模型的时频差异。如图2, 显示了模型④— ⑨的单道地震信号及MEWVD频谱特征, 可见, 窄薄各异的数值模型的反射信号, MEWVD均表现出了较强的频谱聚焦性和高分辨特征。其中, 模型④— ⑥较宽、较厚, 单道地震信号表现出明确的顶底界反射, MEWVD在频率轴方向表现出较窄的频谱分布。模型⑦— ⑨较窄、较薄, 单道地震信号仅显示出了顶界面的反射, 底界面受波场调谐作用而难以识别。MEWVD在频率轴方向表现出较宽的频谱分布, 且厚度越薄、频谱分布范围越宽。
在单频剖面上, MEWVD的频谱聚焦性非常显著, 不同频率的功率谱展现出了不同的时频分辨率。如图3, 显示了正演地震信号的瞬时振幅及25, 35, 45 Hz的MEWVD频谱特征。其中, 图3a显示瞬时振幅受时频分辨能力的局限, 仅较厚的模型④表现出了顶底界面清晰的“ 双轴” 强瞬时振幅特征, 而厚度逐渐减薄的模型⑤— ⑨却表现为“ 单轴” 强瞬时振幅分布, 连厚度大于1/4波长的模型⑤和⑥的顶、底界面也难以有效识别。图3b— 图3c分别显示了3种频率的MEWVD 频谱分布, 且均比瞬时振幅具有更高的分辨率, 但不同频率的MEWVD频谱响应各异。其中, 45 Hz的MEWVD频谱分辨率最高, 可以准确识别窄薄模型⑤和⑥的顶底界面及横向宽度。35 Hz的MEWVD频谱分辨率次之, 也能够有效识别⑤和⑥的顶底界面及横向宽度, 对极窄、极薄模型⑧和⑨的频谱聚焦响应优良, 且受不同宽度的影响, 二者的频谱分布厚度相近、宽度差异显著。进一步分析极窄、极薄模型⑧和⑨, 二者厚度相同、宽度差5倍, 瞬时振幅横向差异不明显, MEWVD频谱异常显著。尤其在图3d中, 45 Hz高频聚焦性极好, 显示出了厚度一样、宽度差5倍的MEWVD频谱异常特征。同时, 虽然25 Hz的MEWVD频谱分辨率总体不如35 Hz和45 Hz的频谱, 但是对极窄、极薄模型⑧和⑨的频谱聚焦响应却有明显差异。其中, 35 Hz和45 Hz的频谱表现为“ 团状” 现象, 而25 Hz的频谱却表现出“ 层状” 频谱聚焦现象, 表明Wigner-Ville最大熵谱与地震信号的强弱无关, 当调谐作用使模型顶底界面的功率谱分布不确定时, 熵最大的功率谱即为最合理的分布。
总之, 仿真地震信号和窄薄正演数值模型实验表明, MEWVD具有良好的频谱聚焦性和时频高分辨能力。MEWVD能够精确表征信号的时频分布, 且不同频率的MEWVD能够准确获取不同尺度窄薄模型的频谱响应; 同时, 由于Wigner-Ville最大熵谱仅与频谱分布概率相关, 当地质体较窄或薄度小于1/4波长而发生调谐作用时, 利用地震反射同相轴不能直接识别, 熵最大的功率谱即为地质体的最合理的分布空间。因此, 综合不同频率的MEWVD频谱特征, 可以识别类似窄薄模型的微型古河道的空间分布。
中江气田面积为2 350 km2, 区内侏罗系沙溪庙组以浅水三角洲沉积为主, 在三角洲平原亚相中发育分流河道、分流间湾、河口坝、决口扇等沉积微相。其中, 分流河道主要发育细— 中粒岩屑长石砂岩、长石岩屑砂岩、岩屑砂岩和岩屑石英砂岩等砂岩储集层, 平均孔隙度达8.94%, 平均渗透率达0.5× 10-3μm2。目前, 针对资源丰富的沙溪庙组气藏, 已投产工业气井约150口, 日产气超过220× 104m3。然而, 随着开发程度的深化, 微型分流河道气藏更加隐蔽, 储集层以条带状叠置薄层河道砂体为主, 横向宽度窄, 纵向厚度薄, 精细刻画难度极大。需要探索窄细微型古河道的精确识别方法, 为气藏描述、水平钻井、水力压裂等提供支撑。
中江气田侏罗系中统沙溪庙组划分为上（J2s1）、中（J2s2）、下（J2s3）共3段; 其中, J2s3段由J2s31、J2s32和J2s33组成。依据古河道和砂岩沉积特征, J2s33层被细分为J2s33-1和J2s33-2小层。目前, J2s33-2小层是主力产气层, 层内网状和条带状古河道主要沿北东流向南西, 延伸长度18~58 km, 宽度0.2~3.5 km, 砂厚6.5~45.2 m。储集层呈席状、箱状和钟型分布, 表现为低频、“ 亮点” 、“ 空白” 等地震反射特征, 沿古河道方向, 地震反射总体较连续、稳定, 垂直河道方向, 呈“ 透镜状” 或弱反射异常。
近几年, 随着中江气田J2s33-2小层气藏的深入开发, 宽度小于500 m、厚度小于35 m的窄薄古河道砂岩成为重点目标, 微型古河道的精确识别也成为气藏开发的关键环节。然而, 受埋藏深度、宽度、厚度、岩性、物性等多种因素的影响, 不同的微型古河道呈现出明显的地震反射差异。如图4a图显示J2s33-2小层的地震主频约30 Hz, 频带8~70 Hz; 图4b显示E井J2s33-2小层古河道出现砂体GR低值异常, 地震反射为局部双波峰“ 透镜状” 异常, 预测河道宽约354 m、厚约35 m; 图4c显示D井J2s33-2小层古河道出现GR低值异常, 预测河道砂岩宽约310 m、厚约30 m, 较E井河道更窄、更薄, “ 透镜状” 反射更弱; 图4d显示C井J2s33-2小层古河道出现GR低值异常, 预测河道砂岩宽约248 m、厚约30 m, 与D、E井相比较, “ 透镜状” 异常消失, 呈弱反射响应。显然, 通过对比图4中的C、D、E井微型古河道地震反射特征, 随着宽度变窄、厚度减薄, 河道隐藏性增强, 地震反射越弱, 识别难度越大。
采用Wigner-Ville最大熵时频谱计算方法, 提取了中江气田沙溪庙组的MEWVD频谱特征, 发现不同宽度和厚度的微型古河道呈现出了不同的MEWVD异常, 揭示了窄细古河道的空间展布。如图5a显示的瞬时振幅剖面上, E井J2s33-2小层古河道出现瞬时振幅强异常, 异常范围较大, 聚焦性较差, 分辨率较低; 图5b— 图5d显示的MEWVD频谱剖面上, E井J2s33-2小层古河道出现强异常, 异常范围随着频率的增加而减小, 聚焦能力增强, 分辨率提高。比较瞬时振幅与25 Hz、35 Hz的MEWVD频谱特征, 在45 Hz的MEWVD剖面上, 分辨能力更高, 纵向频谱被分开, 横向频谱收窄, 频谱异常由“ 团状” 变化为“ 双层状” , 与砂体GR低值异常更加吻合, 有效地揭示出了E井J2s33-2小层35 m薄古河道砂岩的顶底界面。可见, 利用MEWVD频谱异常, 能够有效刻画窄细古河道的纵向分布。
同时, 利用MEWVD频谱异常, 还能刻画微型古河道的横向展布。如图6a瞬时振幅高值异常刻画出了相对较宽的Ⅰ — Ⅳ 、Ⅵ 、Ⅷ 号主河流的沿层展布; 但在黑色框示区域, 河道很窄、河道砂岩很薄, 地震反射能量较强, 分辨率低, 利用瞬时振幅无法识别窄细古河道的空间位置、宽度和流向。在图6b中, 利用MEWVD方法获得的25 Hz频谱出现高值异常, 不仅刻画出了较宽的Ⅰ — Ⅳ 、Ⅵ 、Ⅷ 号主河流的沿层展布, 而且中低值异常刻画出了更窄的Ⅴ 和Ⅶ 号河道的空间展布; 尤其在左侧黑色框内, Ⅴ 和Ⅶ 号河流的宽度和流向非常清晰。与Ⅰ — Ⅷ 微型河道比, Ⅸ 号河道更窄, 且流向弯曲变化较大, 隐蔽性更强, 利用瞬时振幅和25 Hz频谱难以识别, 而在35 Hz和45 Hz的MEWVD频谱特征中, Ⅸ 号曲流河道的沿层特征却非常清晰。由图6可见, MEWVD较瞬时振幅分辨更高, 且不同频率的MEWVD频谱特征, 可以识别不同宽度的窄细河道。当然, 由于任何物质都有相应的敏感频段, 河道也不例外, Ⅴ 和Ⅶ 号河道对低频更敏感, 故在35 Hz和45 Hz的MEWVD频谱特征中反而不如25 Hz刻画更清晰（见图6黑色方框）, 但整体规律仍然是高频频谱刻画效果更佳。
采用RGB融合（红、绿、蓝混合生成其他色彩）显示方式, 将不同频率的MEWVD频谱进行融合, 能更清晰的展示不同宽度和厚度河道的空间分布。如图7a显示了采用RGB融合方式, 将25, 35, 45 Hz的MEWVD频谱进行融合, 有效刻画出J2s33-2小层Ⅰ — Ⅸ 号古河道的展布; 尤其是具有极窄特征的Ⅴ 、Ⅶ 和Ⅸ 号河道, 以及Ⅷ 号河道以东、Ⅵ 号河道以西的更多微型古河道（红色箭头所指）, 也精确地揭示出了其宽度和流向等空间信息。当然, 结合区内沉积构造背景, 利用岩心分析与测井解释资料, 可以证实采用MEWVD频谱异常识别微型古河道的可靠性。如图7b所示, 在纵向上, D井主要发育分流河道、分流间弯、河口坝等三角洲前缘沉积微相。其中, J2s33-2小层以分流河道沉积微相为主, 声波时差和伽马测井曲线表现出低值响应特征, 深侧向和浅侧向电阻率曲线表现出高值异常。在分流河道中部2 869.57~2 871.49 m井段的钻井取心显示, 古河道水平层理发育, 在较强的水动力作用下, 沉积形成了厚度约30 m的致密砂岩储集层。受岩性、物性、厚度等因素影响, 储集层与围岩的阻抗差异不大, 在图4c所示的“ 透镜状” 河道反射较弱, 但在图6和图7却具有较突出的MEWVD频谱异常。可见, 综合地震反射、钻井取心、测井响应等信息, 利用MEWVD频谱异常刻画微型古河道具有较高可靠性。
此外, MEWVD频谱异常还可以刻画断层分布特征, 但断层与古河道具有显著的MEWVD差异。如图6d和图7a显示, 中江气田沉积构造环境整体比较稳定, 在试验工区南部仅发育少量近东西向的断层, 裂缝不发育。断层的MEWVD频谱异常表现为上盘与下盘边界频谱较强、中间的频谱较弱。断层与古河道虽皆呈细长条带状, 但古河道边界频谱较弱、中间频谱较强。可见, 断层与古河道的MEWVD频谱异常虽有相似之处, 却也存在显著差异。本质上, 二者的共性与差异皆由地质环境决定, 这些地震响应被MEWVD异常精确的刻画了出来。
总之, MEWVD频谱聚焦性较强, 时频分辨率较高, 且不同频率的MEWVD对窄细不同的河道具有不同程度的敏感性, 有效的揭示了中江气田沙溪庙组微型古河道的宽度、砂岩厚度和流向等空间信息。
熵最大时的频谱特征, 最能代表地震信号的频谱分布状态, MEWVD不仅避免了WVD产生的交叉项干扰, 而且最大程度的增强了频谱聚焦性, 有效提升了地震信号的时频分辨率。
仿真地震信号和窄薄模型正演模拟信号试验揭示, 较STFT、GST、CWT、WVD、PWVD、SPWVD、CWD等方法, MEWVD能更加精确的提取信号频谱特征; 尤其是针对窄薄各异的数值模型, 由于最大熵的差异而表现出了特殊的频谱响应, 采用不同频率的MEWVD能够准确的识别尺度各异的窄薄地质体。
基于最大熵准则与Wigner-Ville分布的微型古河道识别方法, 充分发挥了MEWVD频谱聚焦性较强、时频分辨率较高等优势, 有效的识别出了中江气田不同尺度的微型古河道, 且不同频率的MEWVD对窄细各异的河道具有不同程度的敏感性, 准确地刻画了微型古河道的宽度、砂岩厚度和流向等空间信息。
符号注释：
$\alpha $— — 自回归系数, 又称为预测误差因子, 无因次; ${{\alpha }_{m}}$— — 序号为$m$的自回归系数, 无因次; ${{A}_{1}}$— — 衰减分量的振幅, dB; ${{A}_{2}}$— — 平稳分量的振幅, dB; E— — 预测误差, W/Hz; ${{E}_{\text{B}m}}$— — 序号为$m$的后向预测误差, W/Hz; $E_{\text{F}m}^{{}}$— — 序号为m的前向预测误差, W/Hz; f— — 信号的频率, Hz; ${{f}_{1}}$— — 衰减分量的频率, Hz; ${{f}_{2}}$— — 平稳分量的频率, Hz; GR— — 自然伽马, API; H— — 熵函数, 无因次; i— — 虚数单位; j— — 自回归系数${{\alpha }_{m}}$的序号, 无因次; ${{k}_{n}}$— — 自相关函数, 无因次; ${{k}_{z}}$— — ${{k}_{n}}$的卷积核, 无因次; l— — 离散信号扩展采样点序号, 无因次; m— — ${{\alpha }_{m}}$、$E_{\text{F}m}^{{}}$、${{E}_{\text{B}m}}$、${{w}_{n}}$等函数的序号, 无因次; M— — AR模型的阶数, 无因次; n— — 离散信号采样点序号, 无因次; N— — 离散信号长度, 无因次; p— — 分布密度函数, 无因次; P— — 功率谱, W/Hz; ${{P}_{0}}$— — 初始功率谱, W/Hz; ${{P}_{m}}$— — 平均功率谱, W/Hz; Q— — 吸收系数, 无因次; Rlld— — 深侧向电阻率, Ω · m; Rlls— — 浅侧向电阻率, Ω · m; t— — 时间, s; W— — 连续WVD函数, W/Hz; ${{W}_{z}}$— — 离散WVD函数, W/Hz; ${{w}_{n}}$— — ${{W}_{z}}$的子项, W/Hz; X— — 信号, 无因次; x— — 某时段信号, dB; ${{x}_{1}}$— — 衰减分量, dB; ${{x}_{2}}$— — 平稳分量, dB; Δ t— — 时间采样率, s; ${{x}_{F}}$— — 仿真信号, dB; ${{\varphi }_{1}}$— — 衰减分量的相位, rad; ${{\varphi }_{2}}$— — 平稳分量的相位, rad; τ — — 时延, s; $z$— — 离散信号, dB。