CN108646291B - 基于果蝇神经网络算法的子波整形反褶积处理方法及装置 - Google Patents

Abstract

本发明提供一种基于果蝇神经网络算法的子波整形反褶积处理方法及装置，涉及地震勘探技术。该方法包括：获取原始地震数据；所述原始地震数据包括各道地震记录；在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波；采用最小平方原理构建维纳滤波器；根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。本发明可以提高地震剖面的分辨率和成像精度，有利于对砂体纵横向展布规律的认识，同时可有效提高后续地震反演的精度，减少其多解性。

Description

基于果蝇神经网络算法的子波整形反褶积处理方法及装置

技术领域

本发明涉及地震勘探技术，尤其涉及一种基于果蝇神经网络算法的子波整形反褶积处理方法及装置。

背景技术

近年来，随着世界油气工业的发展趋势以及油气勘探开发的深入，油气勘探的主体目标已经由构造油气藏转为岩性油气藏，地震勘探从平原地区转入山区，沙漠，高原地区，从物性良好的厚层转向陆相碎屑岩薄互层。因此，就目前的研究形势，提高地震资料的分辨率和信噪比已成为地震勘探的关键和焦点。

在地震数据处理中，反褶积技术是提高地震资料分辨率的重要手段之一，而子波提取方法的深化与完善是地震信号反褶积处理的一个重要研究内容和关键技术，因此如何通过地震子波的提取及反褶积方法来发展一种好的反褶积方法提高地震资料的分辨率和计算效率就显得尤为重要。

发明内容

本发明实施例的主要目的在于提供一种基于果蝇神经网络算法的子波整形反褶积处理方法及装置，以获得高分辨率反褶积因子，提高地震资料处理精度。

为了实现上述目的，本发明实施例提供一种基于果蝇神经网络算法的子波整形反褶积处理方法，包括：

获取原始地震数据；所述原始地震数据包括各道地震记录；

在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波；

采用最小平方原理构建维纳滤波器；

根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；

将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

具体的，所述在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波，包括：

步骤1、建立待估计自回归滑动平均子波模型；

步骤2、根据待估计自回归滑动平均子波模型与满足反射系数序列假设的随机序列合成地震数据；

步骤3、根据合成的地震数据，运用奇异值分解法和基于奇异值分解的总体最小二乘法分别确定待估计自回归滑动平均子波模型的阶数p、q以及自回归参数AR和滑动平均参数MA；

步骤4、在阶数p、q以及自回归参数AR和滑动平均参数MA的基础上，确定模型参数向量θ的搜索范围；

步骤5、在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x；

步骤6、计算E(θ)的解x的评价函数值；

步骤7、在所述评价函数值的变化率小于预先设置的变化阈值时，对待估计自回归滑动平均子波模型的阶数p、q进行扰动，更新阶数p、q；在步骤7之后返回执行步骤3；

步骤8、在所述评价函数值的变化率大于等于预先设置的变化阈值时，将评价函数值对应的当次扰动前的阶数p、q以及E(θ)的解x作为待估计自回归滑动平均子波模型的最优解；

步骤9、根据待估计自回归滑动平均子波模型的最优解确定每道地震记录的子波。

具体的，所述在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x，包括：

随机产生N个果蝇个体，其中每个果蝇个体表示一组决定地震子波的自回归滑动平均参数，所述自回归滑动平均参数包括自回归参数AR和滑动平均参数MA；

根据味道浓度判定函数：确定N个果蝇个体对应的味道浓度值Smell，并选取味道浓度最大的值时对应的自回归滑动平均参数；其中η为放大因子，e_i表示第i组组决定地震子波的自回归滑动平均参数的味道浓度判定函数与累积量的拟合误差，N表示果蝇个体的样本数量。

具体的，所述采用最小平方原理构建维纳滤波器，包括：

利用最小均方误差获取最优维纳滤波器系数w，使得维纳滤波器的输出信号y(n)与期望信号s(n)误差的均方值最小；

所述最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)；其中，w_opt表示获取得到的最优维纳滤波器系数；E(xx^H)为输入信号的自相关；E(sx)为期望信号与输入信号的互相关。

具体的，所述根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数，包括：

将所述每道地震记录的子波作为最优维纳滤波器系数对应的维纳滤波器的输入信号，将所述宽带巴特沃斯子波作为最优维纳滤波器系数对应的维纳滤波器的期望信号，得到最优维纳滤波器系数对应的维纳滤波器的输出结果为子波整形反褶积系数。

一种基于果蝇神经网络算法的子波整形反褶积处理装置，包括：

原始地震数据获取单元，用于获取原始地震数据；所述原始地震数据包括各道地震记录；

子波提取单元，用于在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波；

维纳滤波器构建单元，用于采用最小平方原理构建维纳滤波器；

子波整形反褶积系数确定单元，用于根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；

子波整形结果获得单元，用于将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

另外，所述子波提取单元，具体用于执行：

步骤1、建立待估计自回归滑动平均子波模型；

步骤2、根据待估计自回归滑动平均子波模型与满足反射系数序列假设的随机序列合成地震数据；

步骤3、根据合成的地震数据，运用奇异值分解法和基于奇异值分解的总体最小二乘法分别确定待估计自回归滑动平均子波模型的阶数p、q以及自回归参数AR和滑动平均参数MA；

步骤4、在阶数p、q以及自回归参数AR和滑动平均参数MA的基础上，确定模型参数向量θ的搜索范围；

步骤5、在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x；

步骤6、计算E(θ)的解x的评价函数值；

步骤7、在所述评价函数值的变化率小于预先设置的变化阈值时，对待估计自回归滑动平均子波模型的阶数p、q进行扰动，更新阶数p、q；在步骤7之后返回执行步骤3；

步骤8、在所述评价函数值的变化率大于等于预先设置的变化阈值时，将评价函数值对应的当次扰动前的阶数p、q以及E(θ)的解x作为待估计自回归滑动平均子波模型的最优解；

步骤9、根据待估计自回归滑动平均子波模型的最优解确定每道地震记录的子波。

此外，所述子波提取单元，具体用于：

随机产生N个果蝇个体，其中每个果蝇个体表示一组决定地震子波的自回归滑动平均参数，所述自回归滑动平均参数包括自回归参数AR和滑动平均参数MA；

根据味道浓度判定函数：确定N个果蝇个体对应的味道浓度值Smell，并选取味道浓度最大的值时对应的自回归滑动平均参数；其中η为放大因子，e_i表示第i组组决定地震子波的自回归滑动平均参数的味道浓度判定函数与累积量的拟合误差，N表示果蝇个体的样本数量。

另外，所述维纳滤波器构建单元，具体用于：

利用最小均方误差获取最优维纳滤波器系数w，使得维纳滤波器的输出信号y(n)与期望信号s(n)误差的均方值最小；

所述最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)；其中，w_opt表示获取得到的最优维纳滤波器系数；E(xx^H)为输入信号的自相关；E(sx)为期望信号与输入信号的互相关。

另外，所述子波整形反褶积系数确定单元，具体用于：

将所述每道地震记录的子波作为最优维纳滤波器系数对应的维纳滤波器的输入信号，将所述宽带巴特沃斯子波作为最优维纳滤波器系数对应的维纳滤波器的期望信号，得到最优维纳滤波器系数对应的维纳滤波器的输出结果为子波整形反褶积系数。

本发明实施例提供的一种基于果蝇神经网络算法的子波整形反褶积处理方法及装置，可拓宽地震记录优势频率的带宽，并有效减少子波旁瓣，突出高频弱反射信息，达到改善地震图像的质量，提高地震资料分辨率的目的，其在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波；采用最小平方原理构建维纳滤波器；根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。本发明实施例提供的方法既能够提高地震资料的分辨率和计算效率，又能够提高弱层和薄层的信噪比。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于果蝇神经网络算法的子波整形反褶积处理方法的流程图一；

图2为本发明实施例提供的一种基于果蝇神经网络算法的子波整形反褶积处理方法的流程图二；

图3为本发明实施例中的在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波的具体实现步骤示意图；

图4为本发明实施例中的Ricker子波滤波器响应示意图；

图5为本发明实施例中的带通子波滤波器响应示意图；

图6为本发明实施例中的宽带Butterworth子波滤波器响应示意图；

图7为时间域的10-40Hz Butterworth子波的频谱示意图；

图8为频率域的10-40Hz Butterworth子波的频谱示意图；

图9为东北某区三维地震叠前处理剖面示意图；

图10为本发明实施例的方法提高分辨率后的地震剖面示意图；

图11为本发明实施例的方法子波整形反褶积处理前后的频谱对比示意图；

图12为本发明实施例提供的一种基于果蝇神经网络算法的子波整形反褶积处理装置的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明实施例提供一种基于果蝇神经网络算法的子波整形反褶积处理方法，包括：

步骤101、获取原始地震数据；所述原始地震数据包括各道地震记录。

步骤102、在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波。

步骤103、采用最小平方原理构建维纳滤波器。

步骤104、根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数。

步骤105、将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

本发明实施例提供的一种基于果蝇神经网络算法的子波整形反褶积处理方法，可拓宽地震记录优势频率的带宽，并有效减少子波旁瓣，突出高频弱反射信息，达到改善地震图像的质量，提高地震资料分辨率的目的，其在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波；采用最小平方原理构建维纳滤波器；根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。本发明实施例提供的方法既能够提高地震资料的分辨率和计算效率，又能够提高弱层和薄层的信噪比。

为了使本领域的技术人员更好的了解本发明，下面列举一个更为详细的实施例，如图2所示，本发明实施例提供一种基于果蝇神经网络算法的子波整形反褶积处理方法，包括：

步骤201、获取原始地震数据，所述原始地震数据包括各道地震记录。

步骤202、在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波。

此处，需要说明的是，基于果蝇神经网络算法的累积量拟合法提取地震子波，是在自回归滑动平均模型(Auto Regressive Moving Average Model，简称ARMA模型)下假设反射系数序列是一个独立的、同分布的、非高斯的随机过程，此时地震数据的高阶累计量与地震子波的高阶累计量相当，这样通过估算地震数据的高阶累计量达到估算地震子波高阶累计量，进而达到估算地震子波目的。

自回归滑动平均模型是模型参量法高分辨率谱分析方法之一。ARMA模型方法是研究平稳随机过程有理谱的典型方法，它比AR模型法与MA模型法有较精确的谱估计及较优良的谱分辨率性能，但其参数估算比较繁琐。

具体的ARMA模型参数估计的方法有很多：

如果模型的输入序列与输出序列均能被测量时，则可以用最小二乘法估计其模型参数，这种估计是线性估计，模型参数能以足够的精度估计出来；

许多谱估计中，仅能得到模型的输出序列，这时，参数估计是非线性的，难以求得ARMA模型参数的准确估值。从理论上推出了一些ARMA模型参数的最佳估计方法，但它们存在计算量大和不能保证收敛的缺点。因此工程上提出次最佳方法，即分别估计AR和MA参数，而与最佳参数估计不同，一般不同时估计AR和MA参数，从而使计算量大大减少。

例如：

地震子波ARMA模型和累积量拟合函数：

通常地震激励假设为一零均值平稳随机过程y(n)，用如下褶积模型描述一道地震记录：

式中，h(n)为地震子波；m和n分别为子波因果、非因果部分的长度；r(n)为反射系数序列；v(n)为加性噪声。该模型有如下假设：

(1)反射系数序列r(n)是零均值、独立同分布的非高斯过程，其方差且至少存在一个k阶累积量满足|γ_kr|＜∞；

(2)环境噪声v(n)是零均值、因果的加性高斯色噪声，其高斯成分远大于非高斯成分；且与r(n)统计独立，因此与y(n)也相互独立；

(3)h(n)为非因果混合相位的地震子波，其非因果性是由地震道和检波器引入的失真造成的。

将地震子波的褶积模型转换为

y(n)＝x(n)+v(n) (2)

其中，x(n)为平稳随机过程，为一ARMA模型在输入为反射系数序列r(n)时的响应，该ARMA模型满足以下差分方程：

式中，p和q分别为ARMA模型的AR(自回归)和MA(滑动平均)部分的阶数，其中a(i)为ARMA模型的AR(自回归)参数，i为AR参数序列编号，b(j)为ARMA模型的MA(滑动平均)参数，j为MA参数编号，对应冲激响应为式(1)中的地震子波h(n)。

为了求取ARMA模型的阶数和参数，引入累积量拟合法。由于奇数阶累积量恒为零，偶数阶累积量不为零，故应采用四阶累积量进行地震子波估计。对地震记录求四阶累积量，应用Bartlett-Brillinger-Rosenblatt公式，由假设条件结合累积量定义和性质，得：

C_4y(t₁,t₂,t₃)＝r_4rM_4h(t₁,t₂,t₃) (4)

式中，C_4y(t₁,t₂,t₃)为地震数据的四阶累积量；M_4h(t₁,t₂,t₃)为估计的地震子波四阶矩；r_4r为反射系数序列峰度。r_4r为一标量，在不影响所求取子波形态的情况下，可作为系数被M_4h(t₁,t₂,t₃)吸收，即地震记录的四阶累积量与子波的四阶矩函数之间仅差一个标量因子，可通过计算地震记录的四阶累积量来估算子波。由于高阶累积量关于其变元是对称的，因此在最小平方误差意义下可建立目标函数：

式中，k为地震子波的长度。目标函数E为地震子波的函数，当E达到最小值时，所对应的一组变量即为子波的估计值。此目标函数的求解是一个非线性多参数多极值的优化问题，需要应用非线性优化算法进行求解。通过对优化算法进行研究，使用了果蝇全局寻优算法，并将其运用到模型求解当中，以改善累积量拟合优化的计算效率和寻优精度。

具体的，如图3所示，上述步骤202中，在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波，具体可以采用如下方式：

步骤1、建立待估计自回归滑动平均子波模型。

步骤2、根据待估计自回归滑动平均子波模型与满足反射系数序列假设的随机序列合成地震数据。

步骤3、根据合成的地震数据，运用奇异值分解法和基于奇异值分解的总体最小二乘法分别确定待估计自回归滑动平均子波模型的阶数p、q以及自回归参数AR和滑动平均参数MA。

步骤4、在阶数p、q以及自回归参数AR和滑动平均参数MA的基础上，确定模型参数向量θ的搜索范围。

步骤5、在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x。

此处，对果蝇全局优化算法介绍如下：

1、基本果蝇优化算法

果蝇优化算法(Fruit fly optimization algorithm，FOA)是一种基于果蝇觅食行为的群智能优化算法。果蝇有很好的嗅觉和视觉器官，能够依靠嗅觉感觉到40公里外的食物源，然后在临近食物源时，依靠敏锐的视觉发现食物的具体位置。果蝇优化算法模拟该过程，基于嗅觉和视觉行为进行迭代搜索，通过对果蝇种群位置中心的优化，最终获得问题的优化解。果蝇算法的基本步骤如下：

步骤①：初始化种群中个体的位置。

步骤②：嗅觉搜索。由个体的当前位置，随机选择方向和位置进行搜索。

步骤③：个体评价。对个体搜索到的新位置，计算其浓度判定值和味道浓度。

步骤④：视觉搜索。选择味道浓度最大的位置，个体根据视觉向该位置搜索。

步骤⑤：判断算法是否结束。是，则输出最优解；否则，转步骤②进行迭代。

2、果蝇神经网络优化算法

与其他优化算法相比，果蝇算法具有可调参数少、程序容易实现等优点。为了提高优化精度和效率，提出如下的改进策略：在原有嗅觉搜索的基础上，增加优秀个体引导的协作搜索策略，使个体有方向地被引导。在评价每一代果蝇个体的同时，保存种群中评价最高的前30％的个体作为引导集合。对于果蝇个体X_i，选择引导集合的种群中心X_r，根据公式1向新的食物源方向搜索，其中，ω表示[0，1]之间的随机数。在算法中增加协作概率Q，果蝇个体在选择新的位置之前，首先产生一个随机数与协作概率比较，大于Q则采取嗅觉搜索，反之则采取协作搜索。

X_i＝X_i+ω(X_r-X_i)

3、基于果蝇算法的累积量拟合法

(1)果蝇群体初始化

随机产生N个果蝇个体，每个个体X_i＝(x₁,x₂,...,x_t)表示一组决定地震子波的ARMA参数，其中包括AR参数和MA参数，其中t表示连接权重的个数。

(2)味道浓度判定函数

果蝇个体代表决定地震子波的ARMA参数，味道浓度判定函数与累积量的拟合误差相关，误差越大，味道浓度越小，反之则越大。

上式中，η是放大因子，e_i表示一个样本的误差，N表示样本的数量。

因此，上述步骤步骤5可以采用如下方式：

随机产生N个果蝇个体，其中每个果蝇个体表示一组决定地震子波的自回归滑动平均参数，所述自回归滑动平均参数包括自回归参数AR和滑动平均参数MA。

根据味道浓度判定函数：确定N个果蝇个体对应的味道浓度值Smell，并选取味道浓度最大的值时对应的自回归滑动平均参数；其中η为放大因子，e_i表示第i组组决定地震子波的自回归滑动平均参数的味道浓度判定函数与累积量的拟合误差，N表示果蝇个体的样本数量。

步骤6、计算E(θ)的解x的评价函数值。

步骤7、在所述评价函数值的变化率小于预先设置的变化阈值时，对待估计自回归滑动平均子波模型的阶数p、q进行扰动，更新阶数p、q。

在步骤7之后返回执行步骤3。

步骤8、在所述评价函数值的变化率大于等于预先设置的变化阈值时，将评价函数值对应的当次扰动前的阶数p、q以及E(θ)的解x作为待估计自回归滑动平均子波模型的最优解。

步骤9、根据待估计自回归滑动平均子波模型的最优解确定每道地震记录的子波。

步骤203、利用最小均方误差获取最优维纳滤波器系数w，使得维纳滤波器的输出信号y(n)与期望信号s(n)误差的均方值最小。

所述最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)；其中，w_opt表示获取得到的最优维纳滤波器系数；E(xx^H)为输入信号的自相关；E(sx)为期望信号与输入信号的互相关。

值得说明的是，维纳滤波器通常用于提取被噪声污染的有用信号，它是以最小均方误差准则进行滤波的，设计维纳滤波器的任务，实际上就是选择滤波器系数w，使其输出信号y(n)与期望信号s(n)误差的均方值最小。

下面对最小均方误差准则下的最优滤波器系数w进行推导：

假设一个以自相关函数或者它的功率谱为前提的平稳随机过程，x(n)＝s(n)+v(n)，其中s(n)是有用信号，即期望得到的信号，称为期望信号，v(n)是噪声信号，x(n)是输入信号，即加噪声以后的信号。滤波器实际输出的观测信号。

则误差定义为：

均方误差为：J＝E(e²)；

此处的目的是推导J达到最小时的滤波器系数。

其中w表示维纳滤波器系数。

从上面式子可以看出，代价函数J是w的一个二次函数，在以w为坐标轴的坐标系中是一个凹函数，有唯一的一个极值点，只要求出该极值点，就可以得到使代价函数J最小的滤波器权值w。

令J对w求偏导：

令其为0，得到极值点：w_opt＝{E(xx^H)}^-1·E(sx)，即最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)。

步骤204、将所述每道地震记录的子波作为最优维纳滤波器系数对应的维纳滤波器的输入信号，将所述宽带巴特沃斯子波作为最优维纳滤波器系数对应的维纳滤波器的期望信号，得到最优维纳滤波器系数对应的维纳滤波器的输出结果为子波整形反褶积系数。

此处，对宽带巴特沃斯子波(简称宽带Butterworth子波)介绍如下：

地震数据处理的理想目标是得到反射系数，即子波是冲激函数。但由于受原始数据的频带和信噪比的限制，最终只能得到具有带限子波的记录。常用的子波是带通子波和瑞克子波(简称Ricker子波)。前者延续时间长，旁瓣波形复杂，后者旁瓣幅度比较大，都不是冲激函数的很好的近似。为此，本发明实施例中提出一种新的子波，它由不同宽度的Ricker子波合成，称为宽带Butterworth子波。其主瓣比较窄、旁瓣幅度小、波形简单，在同样的主瓣宽度的情况下，振幅谱的峰值频率比较低，因而在保真度和信噪比方面都优于常用的带通子波和Ricker子波。

1、Ricker子波

Ricker子波的表达式为

其振幅谱为

式中g为Ricker子波表征频率的一个参数，等于峰值频率。如图4所示为Ricker子波滤波器响应。

2、带通子波

带通子波的表达式为

B(t)＝sinπ(f_h-f_l)tcosπ(f_h+f_l)t

其振幅谱为

式中，f_l和f_h分别为低截频和高截频。如图5所示为带通子波滤波器响应。

3、宽带Butterworth子波

时间域宽带Butterworth子波表达式为：

频率域宽带Butterworth子波表达式为：

高分辨率宽带Butterworth子波期望输出算子中，P_f，Q_f为其表征频率参数，改变P_f影响谱的左支(低频)，改变Q_f影响谱的右支(高频)，如图6所示为宽带Butterworth子波滤波器响应。

图5和图6分别为带通滤波算子及宽带Butterworth子波算子的频率响应，可以确定相比带通滤波器，宽带Butterworth滤波器低频响应较陡，有利于地震记录中面波的衰减，高频响应较缓且有较宽的通放带，有利于地震记录的高分辨率处理。

图7和图8分别为时间域的10-40Hz Butterworth子波的频谱，以及频率域的10-40Hz Butterworth子波的频谱。

Butterworth子波最大的特点是最小相位和物理上可实现的。Butterworth子波在时间为零处开始，其在通频带内最平直，需要4个参数限定，从通带的上限频率和下限频率开始。然后是所需的上限和下限的截止陡度，即通带外的振幅响应的陡度。

步骤205、将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

为了说明上述步骤201至步骤205的具体实施过程和技术效果，下面列举一个具体的应用实例：

图9为东北某区三维地震叠前处理剖面，图10为本发明实施例的方法提高分辨率后的地震剖面。地震单炮资料的地质位置一半在河漫滩，一半在河床，左边含砂吸收严重，右边比较湿润，耦合性好，分辨率较高，使用本发明实施例的方法处理后可以看出剖面的纵向分辨率有显著的提高，同时子波一致性变好。

在图10的剖面中上部可以看分辨出相对于图9呈现更多的地质层位，并且各波组特征更加明显，更利于各地层层位横向连续追踪对比。在图10的剖面下部所表现的绕射波特征更加明显，分辨率更高。图11分别为本发明实施例的方法子波整形反褶积处理前后的频谱对比，处理前的地震剖面频带较窄，主频在8-10hz，处理后剖面频带变宽，主频提高，在18-21hz之间，对比分析的结果可以看出，本发明实施例的方法在该区可以将原始资料，主频提高8-10Hz频带拓宽20Hz以上。使用发明实施例中的方法可以提高地震剖面的分辨率和成像精度，有利于对砂体纵横向展布规律的认识，同时可有效提高后续地震反演的精度，减少其多解性。

对应于上述的方法实施例，如图12所示，本发明实施例还提供一种基于果蝇神经网络算法的子波整形反褶积处理装置，包括：

原始地震数据获取单元31，用于获取原始地震数据；所述原始地震数据包括各道地震记录。

子波提取单元32，用于在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波。

维纳滤波器构建单元33，用于采用最小平方原理构建维纳滤波器。

子波整形反褶积系数确定单元34，用于根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数。

子波整形结果获得单元35，用于将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

另外，所述子波提取单元32，具体用于执行：

步骤1、建立待估计自回归滑动平均子波模型。

步骤2、根据待估计自回归滑动平均子波模型与满足反射系数序列假设的随机序列合成地震数据。

步骤3、根据合成的地震数据，运用奇异值分解法和基于奇异值分解的总体最小二乘法分别确定待估计自回归滑动平均子波模型的阶数p、q以及自回归参数AR和滑动平均参数MA。

步骤4、在阶数p、q以及自回归参数AR和滑动平均参数MA的基础上，确定模型参数向量θ的搜索范围。

步骤5、在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x。

步骤6、计算E(θ)的解x的评价函数值。

步骤7、在所述评价函数值的变化率小于预先设置的变化阈值时，对待估计自回归滑动平均子波模型的阶数p、q进行扰动，更新阶数p、q；在步骤7之后返回执行步骤3。

步骤8、在所述评价函数值的变化率大于等于预先设置的变化阈值时，将评价函数值对应的当次扰动前的阶数p、q以及E(θ)的解x作为待估计自回归滑动平均子波模型的最优解。

步骤9、根据待估计自回归滑动平均子波模型的最优解确定每道地震记录的子波。

此外，所述子波提取单元32，具体用于(上述步骤5的具体实现过程)：

随机产生N个果蝇个体，其中每个果蝇个体表示一组决定地震子波的自回归滑动平均参数，所述自回归滑动平均参数包括自回归参数AR和滑动平均参数MA。

根据味道浓度判定函数：确定N个果蝇个体对应的味道浓度值Smell，并选取味道浓度最大的值时对应的自回归滑动平均参数；其中η为放大因子，e_i表示第i组组决定地震子波的自回归滑动平均参数的味道浓度判定函数与累积量的拟合误差，N表示果蝇个体的样本数量。

另外，所述维纳滤波器构建单元33，具体用于：

利用最小均方误差获取最优维纳滤波器系数w，使得维纳滤波器的输出信号y(n)与期望信号s(n)误差的均方值最小。

所述最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)；其中，w_opt表示获取得到的最优维纳滤波器系数；E(xxH)为输入信号的自相关；E(sx)为期望信号与输入信号的互相关。

另外，所述子波整形反褶积系数确定单元34，具体用于：

将所述每道地震记录的子波作为最优维纳滤波器系数对应的维纳滤波器的输入信号，将所述宽带巴特沃斯子波作为最优维纳滤波器系数对应的维纳滤波器的期望信号，得到最优维纳滤波器系数对应的维纳滤波器的输出结果为子波整形反褶积系数。

本发明实施例提供的一种基于果蝇神经网络算法的子波整形反褶积处理装置，可拓宽地震记录优势频率的带宽，并有效减少子波旁瓣，突出高频弱反射信息，达到改善地震图像的质量，提高地震资料分辨率的目的，其在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波；采用最小平方原理构建维纳滤波器；根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。本发明实施例提供的方法既能够提高地震资料的分辨率和计算效率，又能够提高弱层和薄层的信噪比。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读取存储介质中，比如ROM/RAM、磁碟、光盘等。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于果蝇神经网络算法的子波整形反褶积处理方法，其特征在于，包括：

获取原始地震数据；所述原始地震数据包括各道地震记录；

在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波，具体包括：

步骤1、建立待估计自回归滑动平均子波模型；

步骤2、根据待估计自回归滑动平均子波模型与满足反射系数序列假设的随机序列合成地震数据；

步骤3、根据合成的地震数据，运用奇异值分解法和基于奇异值分解的总体最小二乘法分别确定待估计自回归滑动平均子波模型的阶数p、q以及自回归参数AR和滑动平均参数MA；

步骤4、在阶数p、q以及自回归参数AR和滑动平均参数MA的基础上，确定模型参数向量θ的搜索范围；

步骤5、在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x，具体包括：

随机产生N个果蝇个体，其中每个果蝇个体表示一组决定地震子波的自回归滑动平均参数，所述自回归滑动平均参数包括自回归参数AR和滑动平均参数MA；

根据味道浓度判定函数：确定N个果蝇个体对应的味道浓度值Smell，并选取味道浓度最大的值时对应的自回归滑动平均参数为E(θ)的解x；其中η为放大因子，e_i表示第i组决定地震子波的自回归滑动平均参数的味道浓度判定函数与累积量的拟合误差，N表示果蝇个体的样本数量；

步骤6、计算E(θ)的解x的评价函数值；

步骤7、在所述评价函数值的变化率小于预先设置的变化阈值时，对待估计自回归滑动平均子波模型的阶数p、q进行扰动，更新阶数p、q；在步骤7之后返回执行步骤3；

步骤8、在所述评价函数值的变化率大于等于预先设置的变化阈值时，将评价函数值对应的当次扰动前的阶数p、q以及E(θ)的解x作为待估计自回归滑动平均子波模型的最优解；

步骤9、根据待估计自回归滑动平均子波模型的最优解确定每道地震记录的子波；

采用最小平方原理构建维纳滤波器；

根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；

将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

2.根据权利要求1所述的方法，其特征在于，所述采用最小平方原理构建维纳滤波器，包括：

利用最小均方误差获取最优维纳滤波器系数w，使得维纳滤波器的输出信号y(n)与期望信号s(n)误差的均方值最小；

所述最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)；其中，w_opt表示获取得到的最优维纳滤波器系数；E(xx^H)为输入信号的自相关；E(sx)为期望信号与输入信号的互相关。

3.根据权利要求1所述的方法，其特征在于，所述根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数，包括：

将所述每道地震记录的子波作为最优维纳滤波器系数对应的维纳滤波器的输入信号，将所述宽带巴特沃斯子波作为最优维纳滤波器系数对应的维纳滤波器的期望信号，得到最优维纳滤波器系数对应的维纳滤波器的输出结果为子波整形反褶积系数。

4.一种基于果蝇神经网络算法的子波整形反褶积处理装置，其特征在于，包括：

原始地震数据获取单元，用于获取原始地震数据；所述原始地震数据包括各道地震记录；

子波提取单元，用于在自回归滑动平均模型下基于果蝇神经网络算法提取每道地震记录的子波，具体用于执行：

步骤1、建立待估计自回归滑动平均子波模型；

步骤2、根据待估计自回归滑动平均子波模型与满足反射系数序列假设的随机序列合成地震数据；

步骤3、根据合成的地震数据，运用奇异值分解法和基于奇异值分解的总体最小二乘法分别确定待估计自回归滑动平均子波模型的阶数p、q以及自回归参数AR和滑动平均参数MA；

步骤4、在阶数p、q以及自回归参数AR和滑动平均参数MA的基础上，确定模型参数向量θ的搜索范围；

步骤5、在所述模型参数向量θ的搜索范围内采用果蝇全局优化算法对地震子波相关的累积量拟合目标函数E(θ)进行参数估计，确定E(θ)的解x，具体包括：

随机产生N个果蝇个体，其中每个果蝇个体表示一组决定地震子波的自回归滑动平均参数，所述自回归滑动平均参数包括自回归参数AR和滑动平均参数MA；

根据味道浓度判定函数：确定N个果蝇个体对应的味道浓度值Smell，并选取味道浓度最大的值时对应的自回归滑动平均参数为E(θ)的解x；其中η为放大因子，e_i表示第i组决定地震子波的自回归滑动平均参数的味道浓度判定函数与累积量的拟合误差，N表示果蝇个体的样本数量；

步骤6、计算E(θ)的解x的评价函数值；

步骤7、在所述评价函数值的变化率小于预先设置的变化阈值时，对待估计自回归滑动平均子波模型的阶数p、q进行扰动，更新阶数p、q；在步骤7之后返回执行步骤3；

步骤8、在所述评价函数值的变化率大于等于预先设置的变化阈值时，将评价函数值对应的当次扰动前的阶数p、q以及E(θ)的解x作为待估计自回归滑动平均子波模型的最优解；

步骤9、根据待估计自回归滑动平均子波模型的最优解确定每道地震记录的子波；

维纳滤波器构建单元，用于采用最小平方原理构建维纳滤波器；

子波整形反褶积系数确定单元，用于根据所述维纳滤波器、所述每道地震记录的子波和预先设置的宽带巴特沃斯子波确定子波整形反褶积系数；

子波整形结果获得单元，用于将所述子波整形反褶积系数与所述原始地震数据进行褶积，获得子波整形结果。

5.根据权利要求4所述的装置，其特征在于，所述维纳滤波器构建单元，具体用于：

利用最小均方误差获取最优维纳滤波器系数w，使得维纳滤波器的输出信号y(n)与期望信号s(n)误差的均方值最小；

所述最优维纳滤波器系数w为w_opt＝{E(xx^H)}^-1·E(sx)；其中，w_opt表示获取得到的最优维纳滤波器系数；E(xx^H)为输入信号的自相关；E(sx)为期望信号与输入信号的互相关。

6.根据权利要求4所述的装置，其特征在于，所述子波整形反褶积系数确定单元，具体用于：

将所述每道地震记录的子波作为最优维纳滤波器系数对应的维纳滤波器的输入信号，将所述宽带巴特沃斯子波作为最优维纳滤波器系数对应的维纳滤波器的期望信号，得到最优维纳滤波器系数对应的维纳滤波器的输出结果为子波整形反褶积系数。