引用本文
徐天吉, 程冰洁, 牛双晨, 秦正晔, 王贞贞. 基于最大熵准则的Wigner-Ville分布与微型古河道刻画方法及应用[J]. 石油勘探与开发, 2021,48(6): 1175-1186.
XU Tianji, CHENG Bingjie, NIU Shuangchen, QIN Zhengye, WANG Zhenzhen. A microscopic ancient river channel identification method based on maximum entropy principle and Wigner-Ville Distribution and its application[J]. Petroleum Exploration & Development, 2021,48(6): 1175-1186.  
Doi: 10.11698/PED.2021.06.09
Permissions
Copyright©2021, 《石油勘探与开发》编辑部
《石油勘探与开发》编辑部 所有
基于最大熵准则的Wigner-Ville分布与微型古河道刻画方法及应用
徐天吉1,2, 程冰洁3,4, 牛双晨5, 秦正晔1, 王贞贞1
1.电子科技大学资源与环境学院,成都 611731
2.电子科技大学长三角研究院（湖州）,浙江湖州 313000
3.成都理工大学“油气藏地质及开发工程”国家重点实验室,成都 610059
4.成都理工大学“地球探测与信息技术”教育部重点实验室,成都 610059
5.四川中成煤田物探工程院有限公司,成都 610072
联系作者简介：程冰洁（1977-）,女,湖北武汉人,博士,成都理工大学“油气藏地质及开发工程”国家重点实验室、“地球探测与信息技术”教育部重点实验室副教授,主要从事油气资源地震勘探方法研究与教学工作。地址：成都市成华区二仙桥东三路1号,邮政编码：610059。E-mail: chengbj@cdut.edu.cn

第一作者简介：徐天吉（1975-）,男,四川射洪人,博士,电子科技大学资源与环境学院研究员,主要从事多波地震勘探、储集层地质力学、人工智能等方法研究与教学。地址：成都市高新区（西区）西源大道2006号创新中心,邮政编码：611731。E-mail: xutianji@uestc.edu.cn

收稿日期: 2021-01-29 修回日期: 2021-11-09
基金项目: 国家自然科学基金面上项目“基于应力诱导各向异性与岩石力学性质的深层页岩储层可压裂性评价方法研究”（42074160）; 国家自然科学基金面上项目“基于各向异性介质弹性参数的页岩 TOC及脆性预测方法”（41574099）; 四川省科技计划“页岩储层智能化含气性预测方法”（2020JDRC0013）
摘要

针对窄细古河道精细刻画难题,利用最大熵准则增强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

关键词: 最大熵准则; Wigner-Ville分布; 频谱聚焦性; 高分辨; 地震; 河流相; 窄细古河道
中图分类号:TE122 文献标志码:A
A microscopic ancient river channel identification method based on maximum entropy principle and Wigner-Ville Distribution and its application
XU Tianji1,2, CHENG Bingjie3,4, NIU Shuangchen5, QIN Zhengye1, WANG Zhenzhen1
1. School of Resources and Environment, University of Electronic Science and Technology of China, Chengdu 611731, China
2. Yangtze Delta Region Institute of University of Electronic Science and Technology of China, Huzhou 313000
3. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Chengdu University of Technology, Chengdu 610059, China
4. College of Geophysics, Chengdu University of Technology, Chengdu 610059, China
5. Sichuan Zhongcheng Institute of Coalfield Geophysical Engineering, Chengdu 610072
Abstract

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.

Keyword: maximum entropy principle; Wigner-Ville Distribution; spectral focusing; high resolution; seismic; fluvial facies; narrow and thin ancient channel
0 引言

平直、蛇曲、辫状等类型多样的河流相, 多年前就已经成为油气领域的研究热点。在俄罗斯西西伯利亚、美国怀俄明州及中国的鄂尔多斯盆地、四川盆地等区域, 已经发现了许多河流相大型和中型油气田。但是随着中国油气勘探开发程度的不断深化, 目前河流相资源目标已经由大型、中型逐渐转向微型油气藏。微型古河道广泛分布在陆相与海陆过渡相的浅、中、深等各类地层, 油气储量丰富, 勘探潜力巨大, 但具有河道宽度窄、砂体薄、油气聚集复杂等不利特征[1, 2]。此外, 随着埋藏深度的增加, 地震反射表现出能量弱、低频为主、频宽窄等不利特征, 垂向与横向分辨率降低, 制约了窄古河道、薄储集层及其他油气信息的精确提取。微型古河道的识别, 已经成为增加井位部署、水平井钻采等勘探开发风险的突出难题。

针对微型古河道的识别难题, 目前主要从地震采集、处理和解释等角度着手, 利用高分辨地震数据获取古河道的横向宽度、垂向厚度等关键参数和空间展布特征。就高分辨地震采集而言, 尽管在高密度、“ 两宽一高” （宽方位、宽频带、高密度）等先进观测系统下采集到的地震资料, 提升了浅中层油气藏预测与描述精度, 但是面对深层地质目标, 至今未能突破震源能量衰减快、检波器难以接收到高频有效信号等激发与接收方面的“ 瓶颈” 问题。就高分辨地震资料处理与解释而言, 尽管形成了反Q补偿[3]、谱整形[4]、频宽拓展[5]、薄层反射系数反演[6]、低频扩展深度学习[7]、高分辨阻抗反演[8, 9]等方法, 可以使地震数据“ 无中生有” 地出现更多的反射轴、波阻抗界面增加等现象, 但存在信噪比降低、波组有效性评估难、甚至出现假反射界面等特征, 必然降低地质预测的可靠性。同时, 埋藏越深、规模越小的微型古河道的识别难度越大; 尤其是极窄、流向复杂的微型古河道, 需要从地震波的传播能量、反射频率等角度出发, 不断探索突破地震分辨率极限的有效方法, 才能逐步解决微型古河道的识别难题。

识别微型古河道的关键, 是基于高分辨时频分析方法提取地震信号的高分辨时频响应特征, 利用地震波的能量、频率、相位等信息, 有效刻画微型古河道的宽度、流向等特征。目前, 反Q补偿、频宽拓展等许多高分辨处理与解释方法, 均以时频分析为基础, 却未实现高分辨时频响应特征提取, 致使类似微型古河道的小规模地质目标难以实现高清识别。其实, 自1946年Gabor首次提出了声波信号的频谱图后[10], 迅速形成了线性、双线性和非线性3大类核函数不同的时频分析方法[11]。其中, Gabor变换、短时傅里叶变换（STFT）、连续小波变换（CWT）、S变换（ST）和广义S变换（GST）等属于线性时频分析方法, 虽具有计算效率较高的优势, 但多数具有窗口依赖性, 且受海森堡（Heisenberg W）不确定性原理控制[12], 时间和频率的分辨能力不能同时提升。自适应STFT[13]、模态分解[14, 15]、同步挤压小波变换[16]、匹配追踪（MP）[17]等属于非线性时频分析方法, 虽具有分辨率较高的优势, 但计算量较大, 鲁棒性不足, 难以适应低信噪比信号时频分析。Wigner-Ville分布（WVD）[18]、伪Wigner-Ville分布（PWVD）、平滑伪Wigner-Ville分布（SPWVD）和Choi-Williams分布（CWD）[19]等属于双线性（又称二次型）时频分析方法, 不受海森堡不确定性原理约束[12], 优良的数学性质为高分辨时频分析奠定了良好的理论基础。

本文基于信号最大熵准则（Maximum Entropy Principle）, 在克服WVD交叉项干扰的前提下, 研究增强WVD频谱聚焦特性和提升地震信号时频分辨率的方法, 以突出窄细古河道响应特征。并采用仿真信号、正演数值模拟等理论数据开展试验与论证分析, 以形成微型古河道识别新方法, 为井位部署、水平钻井等提供支撑。

1 方法原理
1.1 最大熵准则

熵（Entropy）, 是指离散随机事件或特定信息的出现概率, 由Shannon C E在1948年首次定义, 主要用于表征随机变量或特征信息的不确定性[20]。1957年, Jaynes E T 提出了最大熵准则, 即在先验信息准确的前提下, 最能代表当前信息状态的概率分布, 是熵最大的概率分布, 换言之, 当符合先验信息的概率分布具有不确定性时, 熵最大的信息分布是合理的[21]。

1967年, Burg将最大熵准则引入信号频谱分析中, 形成了最大熵时频分析Burg算法[22]。Burg认为, 熵是信号$X$的分布函数, 仅与$X$的取值概率有关, 而与实际取值无关。当信号$X$的取值为$x$时, 分布密度函数为$p(x)$, 则其熵$H(x)$为：

$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）

当利用（2）式求取功率谱时, 需要在$H[P(f)]$最大的条件下, 递推求解信号$X$的自相关函数, 使功率谱误差最小, 分辨率获得提升。其实, Burg算法等价于线性自回归AR（Auto Regressive）预测模型, 最大熵功率谱$P(f)$的计算公式[22]为：

$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}$

1.2 离散Wigner-Ville分布与功率谱

1992年Boashash定义[23]连续时序信号$x(t)$的解析信号为$z(t)$时, Wigner-Ville分布（WVD）为：

$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$

据WVD的性质可知, 离散WVD是对瞬时自相关函数进行离散傅立叶变换。离散信号的瞬时自相关函数为：

${{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）

当然, 受WVD双线性的影响, 在利用瞬时自相关函数${{k}_{n}}(l)$前后数据相乘时, 将产生交叉项干扰, 不利于高精度功率谱特征分析。

1.3 最大熵准则约束下的Wigner-Ville功率谱计算方法

当瞬时自相关函数${{k}_{n}}(l)$的一阶导数为0时, 求得的Wigner-Ville分布即为最大熵谱。由于一维最大熵谱与自回归模型AR谱是等效的, 故只需要求得AR模型参数、计算出AR模型谱, 就能获得最大熵谱[24]。这样, 可以避免利用（10）式计算功率谱, 绕开瞬时自相关函数前后数据相乘时产生的交叉项干扰。

基于Burg算法求取最大熵功率谱$P(f)$时, （3）式与（10）式是等效的, 只需基于Levinson-Durbin递推规则[25], 在求取预测误差$E$、自回归系数$\alpha $等AR模型的参数后, 代入（3）式就能获得最大熵准则约束下的WVD功率谱。Burg算法通过前向预测误差和后向预测误差, 递推AR模型的自回归系数[26, 27]。将${{E}_{\text{F}, 0}}(n)={{E}_{\text{B}, m-1}}(n)=z(n)z'(n)$设为初始条件, 递推公式如下：

$\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）。

概括起来, Burg算法的WVD最大熵功率谱的计算步骤[28]主要包括：①设$m=1$, 计算${{P}_{0}}(f)$、${{E}_{\text{F}, 0}}(n)$与${{E}_{\text{B}, 0}}(n)$; ②计算${{\alpha }_{m}}(m)$; ③计算阶次为$m$时的第j个自回归系数${{\alpha }_{m}}(j)$, $j=1\ \ 2\ \ \cdots m-1$; ④计算${{P}_{m}}(f)$; ⑤计算阶次为$m$时的${{E}_{\text{F}, m}}(n)$和${{E}_{\text{B}, m}}(n)$; ⑥设$m=m+1$, 重复②— ⑤, 直到${{P}_{m}}(f)$不再明显减小时, 利用Levinson-Durbin递推关系, 可以获得AR模型的参数, 最终计算出WVD最大熵功率谱。

在时间域, MEWVD并未将瞬时自相关函数前后数据相乘, 而仅使用当前$n$点附近的数据, 避免了交叉项干扰; 在频率域, 受最大熵准则约束, 使每个$n$点对应的功率谱更精确, 能量更聚焦。同时, 在最大熵的前提下, 采用Levinson-Durbin约束递推未知瞬时自相关函数, 使计算的功率谱误差减小, 提升了信号时频分辨能力。

需要指出的是, 最大熵准则在数学、经济、信息等领域被广泛应用, 主要优势在于其不仅保留了信息分布的不确性时, 同时将预测风险降到了最低。在油气领域, 常常面对横向较窄、纵向较薄的古河道、砂坝、点礁等微型地质体, 地震波衰减、调谐、噪声等因素增加了微型地质体的预测难度, 计算高精度、高聚焦性的地震功率谱, 是识别微型地质体的重要基础。MEWVD在熵最大时计算最能代表微型地质体的地震频谱响应特征, 不仅能实现时频特征的分类取优, 而且可以在噪声环境下提高信号分析的鲁棒性[29], 有利于提升微型地质体识别的可靠性。

2 方法实验
2.1 仿真地震信号的MEWVD实验及定量分析

地震波在地下传播时, 受地层埋深、孔隙度、填充流体等因素影响, 将产生能量衰减, 引起动力学、运动学和几何学等特征发生变化。采用振幅、吸收系数、频率、相位和走时等参数, 可以仿真描述地震波传播过程, 如下：

$\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}}$振荡周期长、波形“ 胖” 。

表1 仿真地震信号的参数描述

同时, 对比不同的时频分析方法, 发现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的分辨率损失更严重。

图1 仿真地震信号及不同时频分析方法呈现出的频谱响应差异

2.2 窄薄数值模型MEWVD实验及效果分析

为了进一步证实Wigner-Ville最大熵时频谱（MEWVD）的频谱聚焦性和时频高分辨能力, 采用表2所示参数, 设计出二维数值模型, 包含了埋深、宽度和厚度各异的窄薄地质体, 其“ 透镜状” 的地质特征与微型古河道的横截面类似。基于二维地震波动方程正演数值模拟方法, 以35 Hz的Rick子波作为点震源激发, 以10 m采样间距、1 ms采样率接收, 获得了正演模拟记录。模型④— ⑨的砂体厚度逐渐变薄、宽度逐渐变窄; 模型④— ⑥较宽、厚度大于1/4波长（约30 m）, 顶底界面的反射同相轴明显, 而⑦— ⑨较窄、厚度小于1/4波长, 地震记录表现为谐波反射, 且⑦和⑧表现出强调谐反射, 极薄、极窄的⑨则表现为弱调谐反射（见图2）。

表2 数值模型参数描述

图2 窄薄各异模型的单道模拟地震信号及MEWVD等频谱特征

分别抽取代表模型④— ⑨的单道地震信号, 计算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最大熵谱与地震信号的强弱无关, 当调谐作用使模型顶底界面的功率谱分布不确定时, 熵最大的功率谱即为最合理的分布。

图3 正演地震信号的瞬时振幅及25, 35, 45 Hz的MEWVD频谱特征

总之, 仿真地震信号和窄薄正演数值模型实验表明, MEWVD具有良好的频谱聚焦性和时频高分辨能力。MEWVD能够精确表征信号的时频分布, 且不同频率的MEWVD能够准确获取不同尺度窄薄模型的频谱响应; 同时, 由于Wigner-Ville最大熵谱仅与频谱分布概率相关, 当地质体较窄或薄度小于1/4波长而发生调谐作用时, 利用地震反射同相轴不能直接识别, 熵最大的功率谱即为地质体的最合理的分布空间。因此, 综合不同频率的MEWVD频谱特征, 可以识别类似窄薄模型的微型古河道的空间分布。

3 应用案例

四川盆地是一个多旋回富气沉积盆地, 截至2020年, 在震旦系、寒武系、志留系、石炭系、三叠系、侏罗系等27套层系中, 已经发现了构造、构造-岩性、岩性-构造和岩性4类气藏, 形成了四川盆地东北部、西部、南部和北部4大工业气区[30, 31]。其中, 川西气区的侏罗系沙溪庙组气藏, 主要分布在西邻龙门山推覆构造带、东接川中隆起区的川西坳陷。在川西东斜坡地区, 发现了埋深小于3 100 m的大型岩性-构造气田— — 中江气田。

中江气田面积为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井微型古河道地震反射特征, 随着宽度变窄、厚度减薄, 河道隐藏性增强, 地震反射越弱, 识别难度越大。

图4 目标地层地震振幅谱及典型井微型分流河道地震反射剖面

采用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频谱异常, 能够有效刻画窄细古河道的纵向分布。

图5 E井瞬时振幅及25, 35, 45 Hz的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黑色方框）, 但整体规律仍然是高频频谱刻画效果更佳。

图6 瞬时振幅及25, 35, 45 Hz的MEWVD分别刻画的J2s33-2小层微型古河道的展布

采用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频谱异常刻画微型古河道具有较高可靠性。

图7 中江气田J2s33-2小层RGB融合MEWVD刻画河道分布及D井沉积相测井解释

此外, MEWVD频谱异常还可以刻画断层分布特征, 但断层与古河道具有显著的MEWVD差异。如图6d和图7a显示, 中江气田沉积构造环境整体比较稳定, 在试验工区南部仅发育少量近东西向的断层, 裂缝不发育。断层的MEWVD频谱异常表现为上盘与下盘边界频谱较强、中间的频谱较弱。断层与古河道虽皆呈细长条带状, 但古河道边界频谱较弱、中间频谱较强。可见, 断层与古河道的MEWVD频谱异常虽有相似之处, 却也存在显著差异。本质上, 二者的共性与差异皆由地质环境决定, 这些地震响应被MEWVD异常精确的刻画了出来。

总之, MEWVD频谱聚焦性较强, 时频分辨率较高, 且不同频率的MEWVD对窄细不同的河道具有不同程度的敏感性, 有效的揭示了中江气田沙溪庙组微型古河道的宽度、砂岩厚度和流向等空间信息。

4 结论

熵最大时的频谱特征, 最能代表地震信号的频谱分布状态, 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。

（编辑 魏玮）

参考文献
[1] 刘成川, 曹廷宽, 卜淘. 薄窄河道致密砂岩气藏高效开发技术对策[J]. 大庆石油地质与开发, 2020, 39(2): 157-165.
LIU Chengchuan, CAO Tingkuan, BU Tao. High-efficiency developing techniques for the tight sand stone gas reservoir with thin and narrow channels[J]. Petroleum Geology & Oilfield Development in Daqing, 2020, 39(2): 157-165. [本文引用:1]
[2] 梁宏刚, 先伟, 李文平. 塔河油田超深层、薄储层、窄河道砂体识别与描述技术[J]. 地质科技情报, 2015, 34(5): 180-183.
LIANG Honggang, XIAN Wei, LI Wenping. Identification and description technique for ultra-deep, thin reservoir, narrow, channel sand , Tahe oil field[J]. Geological Science and Technology Information, 2015, 34(5): 180-183. [本文引用:1]
[3] WANG Y. Inverse Q-filter for seismic resolution enhancement[J]. Geophysics, 2006, 71(3): 51-60. [本文引用:1]
[4] 李胜强, 刘振东, 严加永, 等. 高分辨深反射地震探测采集处理关键技术综述[J]. 地球物理学进展, 2020, 35(4): 1400-1409.
LI Shengqiang, LIU Zhendong, YAN Jiayong, et al. Review on the key techniques for high resolution deep reflection seismic acquisition and processing[J]. Progress in Geophysics, 2020, 35(4): 1400-1409. [本文引用:1]
[5] MICHAL S, GARY P, JAIME S, et al. Extending seismic band width using the continuous wavelet transform[J]. First Brieak, 2008, 26(6): 97-102. [本文引用:1]
[6] CHOPRA S, CASTAGNA J P, XU Y. Thin-bed reflectivity inversion and some applications[J]. First Break, 2009, 27(5): 27-34. [本文引用:1]
[7] SUN H Y, DEMANET L. Low-frequency extrapolation with deep learning[R]. Anaheim: 88th Society of Exploration Geophysicists International Exposition and Annual Meeting, 2018. [本文引用:1]
[8] 李祺鑫, 罗亚能, 张生, 等. 高分辨率波阻抗贝叶斯序贯随机反演[J]. 石油地球物理勘探, 2020, 55(2): 389-398.
LI Qixin, LUO Yaneng, ZHANG Sheng, et al. High-resolution Bayesian sequential stochastic inversion[J]. Oil Geophysical Prospection, 2020, 55(2): 389-398. [本文引用:1]
[9] 陈彦虎, 毕建军, 邱小斌, 等. 地震波形指示反演方法及其应用[J]. 石油勘探与开发, 2020, 47(6): 1149-1158.
CHEN Yanhu, BI Jianjun, QIU Xiaobin, et al. A method of seismic meme inversion and its application[J]. Petroleum Exploration and Development, 2020, 47(6): 1149-1158. [本文引用:1]
[10] GABOR D. Theory of communication[J]. Journal of the Institute of Electrical Engineers of Japan, 1946, 93(26): 429-457. [本文引用:1]
[11] 黄昱丞, 郑晓东, 栾奕, 等. 地震信号线性与非线性时频分析方法对比[J]. 石油地球物理勘探, 2018, 53(5): 975-989.
HUANG Yucheng, ZHENG Xiaodong, LUAN Yi, et al. Comparison of linear and nonlinear time-frequecy analysis on seismic signals[J]. Oil Geophysical Prospection, 2018, 53(5): 975-989. [本文引用:1]
[12] WAERDEN B L van der. Sources of quantum mechanics[M]. Amsterdam: North-Holland Publishing Co, 1967. [本文引用:2]
[13] PEI S C, HUANG S G. STFT with adaptive window width based on the chirp rate[J]. IEEE Transactions on Signal Processing, 2012, 60(8): 4065-4080. [本文引用:1]
[14] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings Mathematical Physical & Engineering Sciences, 1998, 454(1971): 903-995. [本文引用:1]
[15] TORRES M E, COLOMINAS M A, SCHLOTTHAUER G, et al. A complete ensemble empirical mode decomposition with adaptive noise[R]. Prague: IEEE International Conference on Acoustics, 2011. [本文引用:1]
[16] DAUBECHIES I, LU J, WU H T. Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[J]. Applied and Computational Harmonic Analysis, 2011, 30(2): 243-261. [本文引用:1]
[17] MALLAT S G, ZHANG Z F. Matching pursuits with time-frequency dictionaries[J]. IEEE Transaction on Signal Processing, 1993, 41(12): 3397-3415. [本文引用:1]
[18] WIGNER E P. On the quantum correction for thermodynamic equilibrium[J]. Physical Review, 1932, 40(5): 749-759. [本文引用:1]
[19] CHOI H I, WILLIAMS W J. Improved time-frequency representation of multicomponent signals using exponential kernels[J]. IEEE Transactions on Acoustics Speech & Signal Processing, 2002, 37(6): 862-871. [本文引用:1]
[20] SHANNON C E. A mathematical theory of communication[J]. Mobile Computing and Communications Review, 1948, 5(1): 3-55. [本文引用:1]
[21] JAYNES E T. Information theory and statistical mechanicsⅡ[J]. Physical Review, 1957, 108(2): 171-190. [本文引用:1]
[22] BURG J P. A new analysis technique for time series data[R]. Enschede, the Netherland s: NATO Advanced Study Institute on Signal Processing with Emphasis on Underwater Acoustics, 1968. [本文引用:2]
[23] BOASHASH B. Estimating and interpreting the instantaneous frequency of a signal: Part I: Fundamentals[J]. Proceedings of the IEEE, 1992, 80(4): 520-538. [本文引用:1]
[24] BOS A V D. Alternative interpretation of maximum entropy spectral analysis[J]. IEEE Trans on Information Theory, 1971, 17(4): 493-494. [本文引用:1]
[25] LEVINSON N. The Wiener RMS (root mean square) criterion in filter design and prediction[J]. Journal of Mathematical Physics, 1946, 25(1/2/3/4): 261-278. [本文引用:1]
[26] KAY S M, MARPLE S L. Spectrum analysis a modern perspective[J]. Proceedings of the IEEE, 1981, 68(11): 1380-1419. [本文引用:1]
[27] MARPLE S L. A new autoregressive spectrum analysis algorithm[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1980, 28(4): 441-454. [本文引用:1]
[28] BURG J P. Maximum entropy spectrum analysis[D]. Stanford, California: Stanford University, 1975. [本文引用:1]
[29] 谢林江, 尹东. 基于最大相关熵准则的鲁棒度量学习算法[J]. 计算机系统应用, 2018, 27(10): 146-153.
XIE Linjiang, YIN Dong. Robust metric learning based on maximum correntropy criterion[J]. Computer Systems & Applications, 2018, 27(10): 146-153. [本文引用:1]
[30] 马永生, 蔡勋育, 赵培荣, 等. 四川盆地大中型天然气田分布特征与勘探方向[J]. 石油学报, 2010, 31(1): 347-354.
MA Yongsheng, CAI Xunyu, ZHAO Peirong, et al. Distribution and further exploration of the large-medium sized gas fields in Sichuan Basin[J]. Acta Petrolei Sinica, 2010, 31(1): 347-354. [本文引用:1]
[31] 魏国齐, 杨威, 刘满仓, 等. 四川盆地大气田分布､主控因素与勘探方向[J]. 天然气工业, 2019, 39(6): 1-12.
WEI Guoqi, YANG Wei, LIU Mancang, et al. Distribution rules, main controlling factors and exploration directions of giant gas fields in the Sichuan Basin[J]. Natural Gas Industry, 2019, 39(6): 1-12. [本文引用:1]