专利名称：基于高阶统计量的线性与非线性融合的地震子波提取方法
技术领域：
本发明属于地震信号处理领域。
背景技术：
地震子波估计是地震信号处理中的关键问题之一，是地震正演、地震反演和反褶积等地震处理技术的基础。高精度的地震子波估计具有重要的实用价值，是油气地震勘探领域高分辨率处理技术亟待解决的重要研究课题。在信号处理和系统理论等领域使用高阶统计量的主要动机与出发点可以归结为(I)抑制加性有色噪声(其功率谱未知)的影响；(2)辨识非因果、非最小相位系统或重构 非最小相位信号；(3)抽取由于高斯性偏离引起的各种信息；(4)检测和表征信号中的非线性以及辨识非线性系统；(5)检测和表征信号中的循环平稳性以及分析和处理循环平稳信号。高阶统计量能够保留系统的相位信息、并且理论上可以完全抑制高斯噪声，在解决非线性、非高斯和非最小相位系统的分析和处理方面已逐渐显示出独特优势，为高精度的地震子波提取提供了可能。由高阶累积量缺陷的分析可知，对于高阶统计的运算来说，数据量越大越有利于提高估计质量，该累积量的缺陷本可用于增加地震记录数据(如以多道联合求解)的方式进行缩减或消除。然而不实用的地方在于一般来说，对于实际的地震记录，地震子波是时空变的，即使同一道记录，浅层和深层的地震子波也是不相同的，这主要是因为大地的滤波作用，使地震子波频率随着地层深度增加而降低。为了能够准确估计出地震子波，在地震子波为时不变模型的假设下，要求参与运算的地震记录长度不能太长，甚至应将浅层和深层数据分开来计算，导致了所求得的地震数据累积量有很大的估计方差。累积量矩阵方程法及累积量拟合法均可从理论上提取出非高斯随机过程的模型参数。然而，这两种方法又各有自己的缺陷相对于非线性参数估计方法(累积量拟合法)，线性化参数估计方法(累积量矩阵方程法)运算速度较快，但由于仅利用了地震记录累积量的特殊切片信息，在对高噪信比、短数据地震记录进行子波提取时，该方法可能存在较大的估计误差；累积量拟合优化方法充分利用了地震记录的累积量信息，子波提取具有更高的精确度和数值稳定性，但用累积量拟合法确定模型初始参数范围较为困难，严重影响了算法的寻优效率。虽然比较稳定，但该方法对模型参数的求取异常困难。这是由于累积量拟合目标函数多为一多维多峰值的非线性函数，通常的优化算法对该目标函数的求解无能为力。在对目标函数的仿真中会发现，该目标函数的拟合误差对模型参数的变化极为敏感，使得线性化优化算法的参数估计结果严重依赖于初始估值，尤其在ARMA描述下，这一特点将更为关出。针对累计量矩阵方程法和累积量拟合法的缺陷，将两者结合，提出一种基于高阶累积量的线性与非线性融合的地震子波提取方法，首先利用矩阵方程法确定子波模型的阶数及参数范围，然后根据累积量拟合公式精确估计模型参数，以弥补两种方法单独使用时的缺陷。考虑利用累积量拟合误差反馈矩阵方程法所得初始估计，以拟合误差的大小对累积量矩阵方程法的估计参数及阈值进行调整，并根据调整过的参数重新确定参数寻优空间，进而求得模型的全局最优解。
发明内容
本发明的目的在于提出一种基于高阶累积量的线性与非线性方法融合的高精度子波提取方法。本发明的特征在于，对ARMA模型子波建立拟合目标函数后，利用线性与非线性相结合的方法求解子波参数。首先采用矩阵方程法对子波模型进行初步估计，然后将得到的子波预估计值用于拟合优化算法初始参数搜索空间的确定，并在此基础上用累积量拟合误差调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到模型准确的阶数及参数值。该方法依次含有以下步骤步骤I数据生成用待估计ARMA子波模型与满足反射系数序列假设的随机序列合 成地震数据记录y(n)；步骤2构建拟合目标函数E( 0 ):基于地震褶积模型的假设，构建多峰多极值目标函数e(0):五(沒)=1沒)]其中 r (,
T1 T2 T3^4ymC3J ^
C4y (T1, T2, T3| 0)分别为地震记录四阶累积量和估计模型的四阶矩归一化后的值，令待估计子波的参数向量为e' = ( , -, ap, b:,…，b,)，子波模型的四阶累积量可表示
n
为CJrw#)=刀冲+,吵(/ +巧,吵(/ + :2，吵(/ + r3,0)其中，h(i，9 )为模型参数
i=-n
为0的ARMA模型的冲激响应；n Smax(1。，Inc), lc, In。分别为参数向量0 !对应子波模型因果、非因果部分的长度，时延Ti的取值区间仅需包含{(0彡T1^ q+3p), max(0,q-p) ^ T 2 < min (q+3p,x 3 : min (2p, x 2)} 0步骤3阶数预估计用累积量矩阵方程法初步确定待估计模型的阶数P，q ;步骤4确定参数搜索空间在给定的模型阶数值p，q下，用矩阵方程法初步估计模型参数，并以此确定模型参数向量9的搜索范围；步骤5拟合优化在给定参数向量的搜索范围内用优化算法对累积量拟合目标函数E ( 0 )进行参数精确估计，寻找最优解0 e ；步骤6评价函数分析将拟合误差作为目标适应度函数，并计算最优解叭的目标适应度函数值，当其有明显增加或降低时，转至步骤8，否则，转至步骤7 ；步骤7阶数调整根据累积量拟合误差对模型的阶数进行调整，生成一组新的P，Q值，并根据拟合误差和生成的新阶数值调整累积量矩阵方程法中PODE或SVD的参数阈值，然后转至步骤4 ；步骤8结束将当次阶数调整前所得阶数P，q及对应精确参数估计值0 e视为模型最优解。本发明的特征之二在于基于地震褶积模型假设构建多峰多极值目标函数E(0)，采用累积量矩阵方程法对子波模型进行初步估计，即利用观测信号累积量组成Hankel矩阵的扩展矩阵来估计AR反因果部分的秩和参数，用估计出的AR反因果部分参数滤波原观测信号，用基于样本高阶累积量估计因果系统AR参数的SVD-TLS算法来估计模型AR因果部分参数，利用AR模型预测残差矩阵有效秩和PODE (对角元素乘积)实验相结合确定MA阶次q，然后将得到的子波预估计值作为改进的遗传优化算法初始值并用于初始参数搜索空间的确定，在此基础上用累积量拟合误差即E( 0 )调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到最佳模型参数向量9，准确的阶数P，q和模型参数a, ap, b1； ...bq0 实验仿真说明本发明通过合成地震记录的子波估计实验来检验线性与非线性结合方法提取地震子波算法的效果。首先合成一长度为800 (采样间隔为Ims)的短数据地震记录。反射系数序列服从伯努利分布，子波为非因果、非最小相位，其参数如表I所示，长度取100ms，并、分别加入5% (噪声与信号的能量比)的高斯有色噪声进行实验。分别采用累积量矩阵方程法和线性与非线性相结合方法提取地震子波，所得子波参数如表I所示，两种方法所提取出的子波波形如图2所示。从表I和图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)也可看出，各地震波波形均有一定的偏差，分析认为是实际地震资料中微曲多次反射以及反射系数序列分布偏离白噪假设时对地震子波波形的滤波效应所致。综上所述，本发明提出的方法可以有效的提高地震子波的估计精度，抗噪性能强，特别是在地震数据短的情况下，也能得到精度比较高的地震子波。
图I、本发明所属方法的计算机流程2、使用矩阵方程法和线性与非线性融合方法提取的子波比较图3、不同数据长度下提取的子波模型参数图4、不同强度噪声成份下提取的子波波形图5、实际地震剖面(a)和真实地震资料中提取的子波波形(b)
具体实施例方式
本发明提出了基于高阶累积量的线性与非线性优化算法相结合的ARMA模型参数估计方法。这一算法充分利用了线性和非线性模型参数辨识方法的优点，且具有以下特
占-
^ \\\ I.相对于MA模型，基于高阶累积量的ARMA模型地震子波提取方法对地震子波的描述还具有参数更为节俭的特点，理论仿真分析和实际数据处理结果表明了该方法的可行性和高效性。2.运用累积量矩阵方程法求解ARMA模型初始值和初始搜索范围，提高了子波提取的运算效率。3.运用累积量拟合法构建拟合目标函数，进而求得基于高阶累积量的ARMA模型描述下的地震子波提取方法，充分地利用了地震记录的累积量信息，使方法有更好的稳定性，参数估计结果受数据量大小以及噪声强度的影响较小。4.利用累积量拟合误差反馈矩阵方程法所得初始估计，以拟合误差的大小对初始估计进行参数调整，并根据调整过的参数重新确定参数寻优空间，进而求得模型的全局最优解，可稳定而准确地从短数据高噪声污染的地震记录中提取出混合相位的地震子波。5.本发明所提出的技术也可以应用于系统辨识领域。本发明按以下步骤实施I.用待估计ARMA子波模型与满足反射系数序列假设的随机序列合成地震数据记录 y (n) o2.基于地震褶积模型假设，构建多峰多极值目标函数E( 0 )。3.用累积量矩阵方程法初步确定待估计模型的阶数P，q，确定模型参数向量9的搜索范围。4.将累积量矩阵方程法得到的模型阶数和参数作为预估计值，在给定参数向量的搜索范围内用优化算法对累积量拟合目标函数E( 0)进行参数精确估计，寻找最优解e eD5.将拟合误差作为目标适应度函数，当其有明显增加或降低时，转至6，否则，根据累积量拟合误差对模型的阶数进行调整，生成一组新的P，q值，转至4。6.将当次阶数调整前所得阶数P，q及对应精确参数估计值0 e视为模型最优解，输出地震子波。具体原理如下本发明利用累积量拟合误差反馈矩阵方程法所得初始估计，以拟合误差的大小对累积量矩阵方程法的估计参数及阈值进行调整，并根据调整过的参数重新确定参数寻优空间，进而求得模型的全局最优解。I.应用累积量矩阵方程法对模型参数进行初步估计
通常地震记录可以假设为一零均值平稳随机过程y(n)，构建一随机ARMA过程如下y(n) =x(n)+v(n)(l)其中x(n)为满足下列差分方程的非高斯信号，是ARMA(p，q)模
型在输入为反射系数序列Hn)时的响应
权利要求
1.基于高阶累积量的线性与非线性融合的地震子波提取方法的思想，其特征在于，该方依次含有以下步骤 步骤(I)初始数据生成y(n):用待估计ARMA子波模型与满足反射系数序列假设的随机序列合成地震记录y(n)； 步骤(2)构建拟合目标函数E(0):基于地震褶积模型的假设，构建多峰多极值目标函数
2.基于高阶累积量的线性与非线性相结合的地震子波提取方法，其特征在于基于地震褶积模型假设构建多峰多极值目标函数E(0)，采用累积量矩阵方程法(线性化方法)对子波模型进行初步估计，即利用观测信号累积量组成Hankel矩阵的扩展矩阵来估计AR反因果部分的秩和参数，用估计出的AR反因果部分参数滤波原观测信号，用基于样本高阶累积量估计因果系统AR参数的SVD-TLS算法来估计模型AR因果部分参数，利用AR模型预测残差矩阵有效秩和PODE (对角元素乘积)实验相结合确定MA阶次q，然后将得到的子波预估计值作为改进的遗传优化算法初始值并用于初始参数搜索空间的确定，在此基础上用累积量拟合误差即E( 0 )作为适应度函数调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到最佳模型参数向量匕、准确的阶数P，q和模型参数(a1, ap b1,…，bq)。
全文摘要
本发明公开一种基于高阶统计量的线性与非线性融合的地震子波提取方法，其特点在于对ARMA模型子波建立拟合目标函数后，利用线性与非线性相结合的方法求解地震子波参数。首先采用累积量矩阵方程法对子波模型进行初步估计，然后将得到的子波预估计值用于拟合优化算法初始参数搜索空间的确定，并在此基础上用累积量拟合误差调整矩阵方程法的阈值，寻找精确寻优区间，进而由累积量拟合方法得到模型准确的阶数及参数值。本发明提出的子波提取方法抗噪性能强，在数据较短时也能提取出较高精度的地震子波，在实际地震数据处理中有很好的应用价值。
文档编号G01V1/28GK102768366SQ201110114850
公开日2012年11月7日 申请日期2011年5月4日 优先权日2011年5月4日
发明者张亚南, 彭星, 戴永寿, 王俊岭, 魏磊 申请人:张亚南, 彭星, 戴永寿