CN102768366A - 基于高阶统计量的线性与非线性融合的地震子波提取方法 - Google Patents

Abstract

本发明公开一种基于高阶统计量的线性与非线性融合的地震子波提取方法，其特点在于：对ARMA模型子波建立拟合目标函数后，利用线性与非线性相结合的方法求解地震子波参数。首先采用累积量矩阵方程法对子波模型进行初步估计，然后将得到的子波预估计值用于拟合优化算法初始参数搜索空间的确定，并在此基础上用累积量拟合误差调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到模型准确的阶数及参数值。本发明提出的子波提取方法抗噪性能强，在数据较短时也能提取出较高精度的地震子波，在实际地震数据处理中有很好的应用价值。

Description

基于高阶统计量的线性与非线性融合的地震子波提取方法

技术领域：

本发明属于地震信号处理领域。 

背景技术：

地震子波估计是地震信号处理中的关键问题之一，是地震正演、地震反演和反褶积等地震处理技术的基础。高精度的地震子波估计具有重要的实用价值，是油气地震勘探领域高分辨率处理技术亟待解决的重要研究课题。 

在信号处理和系统理论等领域使用高阶统计量的主要动机与出发点可以归结为：(1)抑制加性有色噪声(其功率谱未知)的影响；(2)辨识非因果、非最小相位系统或重构非最小相位信号；(3)抽取由于高斯性偏离引起的各种信息；(4)检测和表征信号中的非线性以及辨识非线性系统；(5)检测和表征信号中的循环平稳性以及分析和处理循环平稳信号。高阶统计量能够保留系统的相位信息、并且理论上可以完全抑制高斯噪声，在解决非线性、非高斯和非最小相位系统的分析和处理方面已逐渐显示出独特优势，为高精度的地震子波提取提供了可能。由高阶累积量缺陷的分析可知，对于高阶统计的运算来说，数据量越大越有利于提高估计质量，该累积量的缺陷本可用于增加地震记录数据(如以多道联合求解)的方式进行缩减或消除。然而不实用的地方在于：一般来说，对于实际的地震记录，地震子波是时空变的，即使同一道记录，浅层和深层的地震子波也是不相同的，这主要是因为大地的滤波作用，使地震子波频率随着地层深度增加而降低。为了能够准确估计出地震子波，在地震子波为时不变模型的假设下，要求参与运算的地震记录长度不能太长，甚至应将浅层和深层数据分开来计算，导致了所求得的地震数据累积量有很大的估计方差。 

累积量矩阵方程法及累积量拟合法均可从理论上提取出非高斯随机过程的模型参数。然而，这两种方法又各有自己的缺陷：相对于非线性参数估计方法(累积量拟合法)，线性化参数估计方法(累积量矩阵方程法)运算速度较快，但由于仅利用了地震记录累积量的特殊切片信息，在对高噪信比、短数据地震记录进行子波提取时，该方法可能存在较大的估计误差；累积量拟合优化方法充分利用了地震记录的累积量信息，子波提取具有更高的精确度和数值稳定性，但用累积量拟合法确定模型初始参数范围较为困难，严重影响了算法的寻优效率。虽然比较稳定，但该方法对模型参数的求取异常困难。这是由于累积量拟合目标函数多为一多维多峰值的非线性函数，通常的优化算法对该目标函数的求解无能为力。在对目标函数的仿真中会发现，该目标函数的拟合误差对模型参数的变化极为敏感，使得线性化优化算法的参数估计结果严重依赖于初始估值，尤其在ARMA描述下，这一特点将更为突出。 

针对累计量矩阵方程法和累积量拟合法的缺陷，将两者结合，提出一种基于高阶累积量的线性与非线性融合的地震子波提取方法，首先利用矩阵方程法确定子波模型的阶数及参数范围，然后根据累积量拟合公式精确估计模型参数，以弥补两种方法单独使用时的缺陷。 

考虑利用累积量拟合误差反馈矩阵方程法所得初始估计，以拟合误差的大小对累积量矩阵方程法的估计参数及阈值进行调整，并根据调整过的参数重新确定参数寻优空间，进而求得模型的全局最优解。 

发明内容：

本发明的目的在于提出一种基于高阶累积量的线性与非线性方法融合的高精度子波提取方法。 

本发明的特征在于，对ARMA模型子波建立拟合目标函数后，利用线性与非线性相结合的方法求解子波参数。首先采用矩阵方程法对子波模型进行初步估计，然后将得到的子波预估计值用于拟合优化算法初始参数搜索空间的确定，并在此基础上用累积量拟合误差调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到模型准确的阶数及参数值。该方法依次含有以下步骤： 

步骤1数据生成：用待估计ARMA子波模型与满足反射系数序列假设的随机序列合成地震数据记录y(n)； 

步骤2构建拟合目标函数E(θ)：基于地震褶积模型的假设，构建多峰多极值目标函数E(θ)：  $E\left(\mathrm{\θ}\right)=\underset{{\mathrm{\τ}}_1}{\mathrm{\Σ}}\underset{{\mathrm{\τ}}_2}{\mathrm{\Σ}}\underset{{\mathrm{\τ}}_3}{\mathrm{\Σ}}{[{\tilde{C}}_{{4y}_m}({\mathrm{\τ}}_1,{\mathrm{\τ}}_2,{\mathrm{\τ}}_3)-C_{4y}({\mathrm{\τ}}_1,{\mathrm{\τ}}_2,{\mathrm{\τ}}_3|\mathrm{\θ})]}^2$ 其中 C_4y(τ₁，τ₂，τ₃|θ)分别为地震记录四阶累积量和估计模型的四阶矩归一化后的值，令待估计子波的参数向量为θ′＝(a₁，…，a_p，b₁，…，b_q)，子波模型的四阶累积量可表示为：  $C_{4y}({\mathrm{\τ}}_1,{\mathrm{\τ}}_2,{\mathrm{\τ}}_3|\mathrm{\θ})=\underset{i=-n}{\overset{n}{\mathrm{\Σ}}}h(i,\mathrm{\θ})h(i+{\mathrm{\τ}}_1,\mathrm{\θ})h(i+{\mathrm{\τ}}_2,\mathrm{\θ})h(i+{\mathrm{\τ}}_3,\mathrm{\θ})$ 其中，h(i，θ)为模型参数为θ的ARMA模型的冲激响应；n≥max(l_c，l_nc)，l_c，l_nc分别为参数向量θ′对应子波模型因果、非因果部分的长度，时延τ_i的取值区间仅需包含：{(0≤τ₁≤q+3p)，max(0，q-p)≤τ₂≤min(q+3p，τ₁)，0≤τ₃≤min(2p，τ₂)}。 

步骤3阶数预估计：用累积量矩阵方程法初步确定待估计模型的阶数p，q； 

步骤4确定参数搜索空间：在给定的模型阶数值p，q下，用矩阵方程法初步估计模型参数，并以此确定模型参数向量θ的搜索范围； 

步骤5拟合优化：在给定参数向量的搜索范围内用优化算法对累积量拟合目标函数E(θ)进行参数精确估计，寻找最优解θ_e； 

步骤6评价函数分析：将拟合误差作为目标适应度函数，并计算最优解θ_e的目标适应度函数值，当其有明显增加或降低时，转至步骤8，否则，转至步骤7； 

步骤7阶数调整：根据累积量拟合误差对模型的阶数进行调整，生成一组新的p，q值，并根据拟合误差和生成的新阶数值调整累积量矩阵方程法中PODE或SVD的参数阈值，然后转至步骤4； 

步骤8结束：将当次阶数调整前所得阶数p，q及对应精确参数估计值θ_e视为模型最优解。 

本发明的特征之二在于基于地震褶积模型假设构建多峰多极值目标函数E(θ)，采用累积量矩阵方程法对子波模型进行初步估计，即利用观测信号累积量组成Hankel矩阵的扩展矩阵来估计AR反因果部分的秩和参数，用估计出的AR反因果部分参数滤波原观测信号，用基于样本高阶累积量估计因果系统AR参数的SVD-TLS算法来估计模型AR因果部分参数，利用AR模型预测残差矩阵有效秩和PODE(对角元素乘积)实验相结合确定MA阶次q，然后将得到的子波预估计值作为改进的遗传优化算法初始值并用于初始参数搜索空间的确定，在此基础上用累积量拟合误差即E(θ)调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到最佳模型参数向量θ，准确的阶数p，q和模型参数a₁，…，a_p，b₁，…b_q。 实验仿真说明： 

本发明通过合成地震记录的子波估计实验来检验线性与非线性结合方法提取地震子波算法的效果。首先合成一长度为800(采样间隔为1ms)的短数据地震记录。反射系数序列服从伯努利分布，子波为非因果、非最小相位，其参数如表1所示，长度取100ms，并分别加入5％(噪声与信号的能量比)的高斯有色噪声进行实验。分别采用累积量矩阵方程法和线性与非线性相结合方法提取地震子波，所得子波参数如表1所示，两种方法所提取出的子波波形如图2所示。 

从表1和图2中可以看出，使用线性与非线性融合方法提地震子波可以在短数据情况下有效地估计出非因果、非最小相位子波。在这个例子中，由于合成地震记录的长度仅为800ms，若仅使用基于累积量的矩阵方程法所得的子波于原始波形相差较大，而采用基于线性和非线性的混合算法的估计误差并不十分明显，相对具有更高的估计精度。 

为了分析在不同数据长度下，线性与非线性方法提取的子波的稳定性，进行如下实验仿真：在25次试验中每次均随机生成长度分别为20000ms，5000ms，2000ms，500ms的合成 地震记录，采用线性与非线性融合方法提取子波，结果如图3所示。由图3可以看出，对于不同数据长度的地震记录，提取的地震子波与真实子波十分接近，因此该方法可提取高精度的地震子波，特别是在地震数据短的情况。 

为了分析在不同噪声强度和成份下，线性与非线性融合方法提取的子波的稳定性，进行如下实验仿真：为考察加性噪声对本发明方法的影响，每次实验均任意生成一段长度为20000ms的合成记录，对该记录分别加入不同强度的高斯、非高斯色噪声后应用本发明方法进行子波提取。 

由图4(a)看出，本发明方法在高斯色噪声成份较高的情况下仍能较为准确的辨识出实际的波形。由图4(b)展示了非高斯色噪声的加入后对子波提取的效果影响并不明显。 

为验证本发明方法的实际应用价值，对一实际地震剖面进行地震子波提取试验。如图5(a)所示，该地震剖面数据为201道，采样率为1ms，数据记录长度为2s。地震子波提取方式为每隔20道均随机选取该地震剖面中的一道应用本发明方法进行模型参数估计，其结果如图5(b)所示。应用本发明方法所提取出的10道地震子波，其波形在相位上保持了一定的连续性，从另一个方面表明了本发明方法的一致性。在实际地震资料的子波提取过程中发现，由于随机性地选择20道中的一道地震记录，对剖面中某些道数据进行子波提取时，子波提取结果并非完全一致，由图5(b)也可看出，各地震波波形均有一定的偏差，分析认为是实际地震资料中微曲多次反射以及反射系数序列分布偏离白噪假设时对地震子波波形的滤波效应所致。 

综上所述，本发明提出的方法可以有效的提高地震子波的估计精度，抗噪性能强，特别是在地震数据短的情况下，也能得到精度比较高的地震子波。 

附图说明：

图1、本发明所属方法的计算机流程图 

图2、使用矩阵方程法和线性与非线性融合方法提取的子波比较 

图3、不同数据长度下提取的子波模型参数 

图4、不同强度噪声成份下提取的子波波形 

图5、实际地震剖面(a)和真实地震资料中提取的子波波形(b) 

具体实施方式：

本发明提出了基于高阶累积量的线性与非线性优化算法相结合的ARMA模型参数估计方法。这一算法充分利用了线性和非线性模型参数辨识方法的优点，且具有以下特点： 

1.相对于MA模型，基于高阶累积量的ARMA模型地震子波提取方法对地震子波的描述还具有参数更为节俭的特点，理论仿真分析和实际数据处理结果表明了该方法的可行性 和高效性。 

2.运用累积量矩阵方程法求解ARMA模型初始值和初始搜索范围，提高了子波提取的运算效率。 

3.运用累积量拟合法构建拟合目标函数，进而求得基于高阶累积量的ARMA模型描述下的地震子波提取方法，充分地利用了地震记录的累积量信息，使方法有更好的稳定性，参数估计结果受数据量大小以及噪声强度的影响较小。 

4.利用累积量拟合误差反馈矩阵方程法所得初始估计，以拟合误差的大小对初始估计进行参数调整，并根据调整过的参数重新确定参数寻优空间，进而求得模型的全局最优解，可稳定而准确地从短数据高噪声污染的地震记录中提取出混合相位的地震子波。 

5.本发明所提出的技术也可以应用于系统辨识领域。 

本发明按以下步骤实施： 

1.用待估计ARMA子波模型与满足反射系数序列假设的随机序列合成地震数据记录y(n)。 

2.基于地震褶积模型假设，构建多峰多极值目标函数E(θ)。 

3.用累积量矩阵方程法初步确定待估计模型的阶数p，q，确定模型参数向量θ的搜索范围。 

4.将累积量矩阵方程法得到的模型阶数和参数作为预估计值，在给定参数向量的搜索范围内用优化算法对累积量拟合目标函数E(θ)进行参数精确估计，寻找最优解θ_e。 

5.将拟合误差作为目标适应度函数，当其有明显增加或降低时，转至6，否则，根据累积量拟合误差对模型的阶数进行调整，生成一组新的p，q值，转至4。 

6.将当次阶数调整前所得阶数p，q及对应精确参数估计值θ_e视为模型最优解，输出地震子波。 具体原理如下： 

本发明利用累积量拟合误差反馈矩阵方程法所得初始估计，以拟合误差的大小对累积量矩阵方程法的估计参数及阈值进行调整，并根据调整过的参数重新确定参数寻优空间，进而求得模型的全局最优解。 

1.应用累积量矩阵方程法对模型参数进行初步估计 

通常地震记录可以假设为一零均值平稳随机过程y(n)，构建一随机ARMA过程如下： y(n)＝x(n)+v(n)    (1) 其中x(n)为满足下列差分方程的非高斯信号，是ARMA(p，q)模型在输入为反射系数序列r(n) 时的响应：  $\underset{i=0}{\overset{p}{\mathrm{\Σ}}}a\left(i\right)x(n-i)=\underset{j=0}{\overset{q}{\mathrm{\Σ}}}b\left(j\right)r(n-j)---\left(2\right)$ v(n)是满足下式的有色高斯观测噪声：  $\underset{i=0}{\overset{p_v}{\mathrm{\Σ}}}c\left(i\right)v(n-i)=\underset{j=0}{\overset{q_v}{\mathrm{\Σ}}}d\left(j\right)e(n-j)---\left(3\right)$

对上述信号作如下假设： 

假设1：反射系数序列r(n)为一般可假设为零均值、独立同分布的非高斯过程。其方差 

有限，且存在k阶累积量满足|r_kr|＜∞。 

假设2：v(n)为环境噪声、仪器噪声及激发产生的多次反射噪声等的合成噪声信号，与r(n)统计独立，一般为加性色噪声，且其高斯成分远大于非高斯部分。 

假设3：与式(2)对应的ARMA(p，q)模型为非因果、非最小相位系统，其中a(0)＝1。其传递函数为：  $H\left(z\right)=\underset{j=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h\left(i\right)z^{-j}=\frac{B\left(z\right)}{A\left(z\right)}=\frac{\underset{j=0}{\overset{q}{\mathrm{\Σ}}}b\left(j\right)z^{-j}}{1+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}a\left(i\right)z^{-i}}---\left(4\right)$ 其中x(n)为稳定信号，所以当|z|＝1时，A(z)≠0，并且H(z)中无零极点对消。 

假设4：H(z)有p₁个圆内极点，有p₂＝p-p₁个圆外极点。令α_j为H(z)的极点，则有 A(z)＝A₁(z)A₂(z)g(z)    (5) 式中  $A_1\left(z\right)=\underset{j=1}{\overset{p_1}{\mathrm{\&Pi;}}}(1-{\mathrm{\&alpha;}}_j{z}^{-1})=1+\underset{j=1}{\overset{p_1}{\mathrm{\&Sigma;}}}{a}_1\left(j\right)z^{-j},\left|{\mathrm{\&alpha;}}_j\right|<1---\left(6\right)$ 称之为H(z)的因果AR部分；而  $A_2\left(z\right)=\underset{j=1}{\overset{p_2}{\mathrm{\&Pi;}}}(1-{\mathrm{\&beta;}}_{j}z)=1+\underset{j=1}{\overset{p_2}{\mathrm{\&Sigma;}}}{a}_2\left(j\right)z^j;$ ${\mathrm{\&beta;}}_j=\frac{1}{{\mathrm{\&alpha;}}_{j+p_1}},\left|{\mathrm{\&alpha;}}_{j+p_1}\right|>1;$ $g\left(z\right)=z^{-p_2}\underset{j=1}{\overset{p_2}{\mathrm{\&Pi;}}}(-{\mathrm{\&alpha;}}_{j+p_1})---\left(7\right)$ 称 

为H(z)的反因果AR部分，显然 

(常数)； 

假设5：H(z)有q₁个圆内零点，有q₂＝q-q₁个圆外零点。则与极点情况类似有 B(z)＝B₁(z)B₂(z)s(z)其中B₁(z)为最小相位部分，B₂(z)s(z)为非最小相位部分。 

应用累积量矩阵方程法求解子波初始解，步骤如下： 

(1)求功率谱等价意义下的因果最小相位模型的AR参数 

由式(1)～(3)可以得到：  $Y\left(z\right)=\frac{B\left(z\right)}{A\left(z\right)}R\left(z\right)+\frac{D\left(z\right)}{C\left(z\right)}E\left(z\right)---\left(8\right)$

其中Y(z)、R(z)、C(z)、D(z)、E(z)分别是y(n)、r(n)、c(i)、d(j)、e(n)的Z变换，根据式(5)和式(7)易知观测信号y(n)的功率谱可以表示为：  $P\left(z\right)=\frac{B\left(z\right)B\left(z^{-1}\right)}{A\left(z\right)A\left(z^{-1}\right)}{\mathrm{\&sigma;}}_r^2+\frac{D\left(z\right)D\left(z^{-1}\right)}{C\left(z\right)C\left(z^{-1}\right)}{\mathrm{\&sigma;}}_e^2=\frac{B\left(z\right)B\left(z^{-1}\right)C\left(z\right)C\left(z^{-1}\right){\mathrm{\&sigma;}}_r^2+A\left(z\right)A\left(z^{-1}\right)D\left(z\right)D\left(z^{-1}\right){\mathrm{\&sigma;}}_e^2}{A\left(z\right)A\left(z^{-1}\right)C\left(z\right)C\left(z^{-1}\right)}---\left(9\right)$

由于A₂(z)的根均在单位圆外，将A₂(z)所对应的ARMA模型的极点取共轭倒数，而反演到单位圆内得到 

则 

的根均在单位圆内。因 

与A₂(z)的功率谱等价，现构造 

则A₁(z)A₂(z)与 

功率谱等价。由功率谱等价原理，并利用Cadzow的自相关法估计AR参数的SVD-TLS(奇异值分解-总体最小二乘)算法可以估计出 这种基于观测信号自相关来估计因果信号AR参数的算法已经比较成熟，且仿真结果证明它有很高的分辨率和精度。 

(2)求H(z)的反因果部分AR参数 

利用估计出的 

去滤波观测信号y(n)得到y1(n)，其Z变换为：  $Y_1\left(z\right)=A_1\left(z\right){\stackrel{\&OverBar;}{A}}_2\left(z\right)C\left(z\right)Y\left(z\right)=\frac{B\left(z\right){\stackrel{\&OverBar;}{A}}_2\left(z\right)C\left(z\right)}{A_2\left(z\right)g\left(z\right)}\&CenterDot;R\left(z\right)+A_1\left(z\right){\stackrel{\&OverBar;}{A}}_2\left(z\right)D\left(z\right)E\left(z\right)---\left(10\right)$

可见y₁(n)是一反因果信号，上式等号右边第二项为一高斯信号，其高阶累积量为零。用基于观测信号样本高阶累积量估计反因果系统AR参数的SVD-TLS算法可以估计出反因果AR参数 

(3)求H(z)的因果部分AR参数及MA参数 

利用估计出的A₂(z)g(z)滤波原观测信号y(n)得到y₂(n)，则其Z变换为：  $Y_2\left(z\right)=A_2\left(z\right)g\left(z\right)Y\left(z\right)=\frac{B\left(z\right)}{A_1\left(z\right)}R\left(z\right)+\frac{A_2\left(z\right)g\left(z\right)D\left(z\right)}{C\left(z\right)}E\left(z\right)---\left(11\right)$

由上式可见y₂(n)是一因果信号和一高斯有色噪声的叠加。利用基于样本高阶累积量估计因果系统AR参数的SVD-TLS算法可估计出因果部分的AR参数a(1)，a(2)，...，a(p₁)。然后利用AR模型预测残差矩阵有效秩和PODE(对角元素乘积)实验相结合可确定MA阶次q。用求得的因果部分AR参数a(0)，a(1)，...，a(p₁)对y₂(n)进行滤波可得到一纯MA过程y₃(n)，其 Z变换为：  $Y_3\left(z\right)=A_1\left(z\right)Y_2\left(z\right)=B\left(z\right)R\left(z\right)+\frac{A_1\left(z\right)A_2\left(z\right)D\left(z\right)}{C\left(z\right)}E\left(z\right)---\left(12\right)$ 最后采用只用高阶累积量确定MA参数b(0)，b(1)，...，b(q)的方法即可得到MA参数。 

2.累积量拟合法优化参数 

基于地震褶积模型的假设，构建目标函数：  $E\left(\mathrm{\&theta;}\right)=\underset{{\mathrm{\&tau;}}_1}{\mathrm{\&Sigma;}}\underset{{\mathrm{\&tau;}}_2}{\mathrm{\&Sigma;}}\underset{{\mathrm{\&tau;}}_3}{\mathrm{\&Sigma;}}{[{\tilde{C}}_{{4y}_m}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)-C_{4y}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3|\mathrm{\&theta;})]}^2---\left(13\right)$ 其中， 

C_4y(τ₁，τ₂，τ₃|θ)分别为地震记录四阶累积量和估计模型的四阶矩归一化后的值。修改参数向量θ，令待估计子波的参数向量为θ′＝(a₁，...，a_p，b₁，...，b_q)，子波模型的四阶累积量C_4y(τ₁，τ₂，τ₃|θ)可表示为：  $C_{4y}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3|\mathrm{\&theta;})=\underset{i=-n}{\overset{n}{\mathrm{\&Sigma;}}}h(i,\mathrm{\&theta;})h(i+{\mathrm{\&tau;}}_1,\mathrm{\&theta;})h(i+{\mathrm{\&tau;}}_2,\mathrm{\&theta;})h(i+{\mathrm{\&tau;}}_3,\mathrm{\&theta;})---\left(14\right)$ 其中，h(i，θ)为模型参数为θ，且描述符合式(3)的ARMA模型的冲激响应；n≥max(l_c，l_nc)，l_c，l_nc分别为参数向量θ′对应子波模型因果、非因果部分的长度。由累积量的对称性，时延τ_i的取值区间仅需包含：{(0≤τ₁≤q+3p)，max(0，q-p)≤τ₂≤min(q+3p，τ₁)，0≤τ₃≤min(2p，τ₂)}即可。 

目标函数E(θ)多峰多极值函数，采用改进的遗传算法求解参数精确解。不仅各参数间的依赖性较大，且函数值模型对某些特定参数尤其是AR部分参数的变化较为敏感。因此，求解该目标函数所用寻优算法须兼有对参数向量整体的全局随机搜索能力和对单个参数的深度搜索能力，遗传算法在搜索空间的深度搜索和广度搜索之间维持了很好的平衡性，结合本课题目标函数的特点，本发明采用改进的遗传算法对其进行寻优，将遍历均匀性好、不容易陷入小循环或者不动点的猫映射(cat map)引入到混沌遗传算法中并进行了改进，提出了较遗传算法性能更为稳定有效的、基于猫映射的混沌遗传算法。计算最优解θ_e的目标函数适应度即拟合误差E(θ)(确定在3.5左右)，当其有明显增加或降低时，继续直到接近目标函数适应度，否则，结束。根据累积量拟合误差对模型的阶数进行调整，生成一组新的p，q值，并根据拟合误差和生成的新阶数值调整累积量矩阵方程法中PODE或SVD的参数阈值。 

相对于单独使用累积量矩阵法和累积量拟合法进行ARMA模型参数估计，上述方法可以在数据量较短且噪声较大的情况下得到较高精度的地震子波，且具有较高的运算效率。 表1子波估计结果比较 

Claims (2)

1.基于高阶累积量的线性与非线性融合的地震子波提取方法的思想，其特征在于，该方依次含有以下步骤：

步骤(1)初始数据生成y(n)：用待估计ARMA子波模型与满足反射系数序列假设的随机序列合成地震记录y(n)；

步骤(2)构建拟合目标函数E(θ)：基于地震褶积模型的假设，构建多峰多极值目标函数：

$E\left(\mathrm{\&theta;}\right)=\underset{{\mathrm{\&tau;}}_1}{\mathrm{\&Sigma;}}\underset{{\mathrm{\&tau;}}_2}{\mathrm{\&Sigma;}}\underset{{\mathrm{\&tau;}}_3}{\mathrm{\&Sigma;}}{[{\tilde{C}}_{{4y}_m}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3)-C_{4y}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3|\mathrm{\&theta;})]}^2$

其中

C_4y(τ₁，τ₂，τ₃|θ)分别为地震记录四阶累积量和估计模型的四阶矩归一化后的值，令待估计子波的参数向量为θ′＝(a₁，…，a_p，b₁，…，b_q)，子波模型的四阶累积量可表示为：

$C_{4y}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,{\mathrm{\&tau;}}_3|\mathrm{\&theta;})=\underset{i=-n}{\overset{n}{\mathrm{\&Sigma;}}}h(i,\mathrm{\&theta;})h(i+{\mathrm{\&tau;}}_1,\mathrm{\&theta;})h(i+{\mathrm{\&tau;}}_2,\mathrm{\&theta;})h(i+{\mathrm{\&tau;}}_3,\mathrm{\&theta;})$

其中，h(i，θ)为模型参数为θ的ARMA模型的冲激响应；n≥max(l_c，l_nc)，l_c，l_nc分别为参数向量θ′对应子波模型因果、非因果部分的长度，时延τ_i的取值区间仅需包含：{(0≤τ₁≤q+3p)，max(0，q-p)≤τ₂≤min(q+3p，τ₁)，0≤τ₃≤min(2p，τ₂)}；

步骤(3)阶数预估计：用累积量矩阵方程法初步确定待估计模型的阶数p，q；

步骤(4)确定参数搜索空间：对给定阶数p，q值，用矩阵方程法对该阶数进行参数预估计，并由此确定解空间(a，b)，确定模型参数向量θ的搜索范围；

步骤(5)拟合优化：在给定解空间内用拟合优化算法对累积量拟合目标函数E(θ)进行参数精确估计，寻找最优解θ_e；

步骤(6)评价函数分析：将拟合误差作为目标适应度函数，并计算最优解θ_e的目标适应度函数值，当其有明显增加或降低时，转至(8)，否则，转至步骤(7)；

步骤(7)阶数调整：根据累积量拟合误差对模型的阶数进行调整，生成一组新的p，q值，并根据拟合误差和生成的新阶数值调整累积量矩阵方程法中PODE或SVD的参数阈值，然后转至步骤(4)；

步骤(8)结束：将当次扰动前所得阶数p，q及参数估计值θ_e视为模型最优解。

2.基于高阶累积量的线性与非线性相结合的地震子波提取方法，其特征在于基于地震褶积模型假设构建多峰多极值目标函数E(θ)，采用累积量矩阵方程法(线性化方法)对子波模型进行初步估计，即利用观测信号累积量组成Hankel矩阵的扩展矩阵来估计AR反因果部分的秩和参数，用估计出的AR反因果部分参数滤波原观测信号，用基于样本高阶累积量估计因果系统AR参数的SVD-TLS算法来估计模型AR因果部分参数，利用AR模型预测残差矩阵有效秩和PODE(对角元素乘积)实验相结合确定MA阶次q，然后将得到的子波预估计值作为改进的遗传优化算法初始值并用于初始参数搜索空间的确定，在此基础上用累积量拟合误差即E(θ)作为适应度函数调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到最佳模型参数向量θ_e、准确的阶数p，q和模型参数(a₁，…，a_p，b₁， …，b_q)。

Images