一种基于高斯平滑的时频域反褶积方法与流程
本发明是关于一种基于高斯平滑的时频域反褶积方法，属于石油地震勘探领域。
背景技术：
随着地震勘探对复杂构造储层精确描述的难度越来越大，地震资料高分辨率、高信噪比的要求也随之越来越高，通过反褶积来压缩子波进而提高地震分辨率的方法已经成为一种最常用的手段。反褶积模型最早是由robinson提出的，其假设地层反射系数满足白谱特性且地震子波是最小相位的，在此基础上建立了传统的反褶积模型，传统的反褶积模型都是基于地震波在地下传播的过程中能量没有发生损失且波形保持不变的前提，但是实际地震波在传播的过程中会受到地下介质的影响，产生能量耗散和速度频散，主频会随着传播时间的增长而向低频偏移，而传统的反褶积模型不能吻合这个动态过程。为更加稳合实际地震波的传播过程，考虑大地滤波效应，margrave将地震记录转换到gabor(伽柏)域，结合平滑函数，对每个时窗振幅谱进行平滑，直接估计出衰减子波振幅谱，提出一种时频域反褶积方法，这种时频域反褶积方法可以提高地震记录的分辨率，但是由于平滑函数自身的问题和gabor变换时窗的固定性，使得该反褶积方法具有一定的局限性。
传统基于傅里叶变换在频率域实现反褶积是提高地震分辨率的重要方法，但是傅里叶变换是对单道地震记录进行谱分析，只能在时域和频域相互映射，缺乏对单道地震记录的时间和频率同时定位的能力，不能突出单道地震记录的局部频谱特征，对单道地震记录在时频谱的分析方法有gabor变换、小波变换和s变换等，然而，gabor变换受时窗固定的缺陷不能自适应时频分辨率变化；小波变换需要对小波基进行合理选择，高频区分辨率差；s变换基本函数固定，在实际资料处理过程中缺乏灵活性。
反褶积方法由静态发展到动态，从频域分析发展到时频域分析，虽然改善了反褶积结果，但由于在谱模拟时使用平滑函数受数据体干扰大，当多项式拟合阶数低时，拟合误差大，阶数高时计算效率低且容易溢出，使用传统的平滑函数拟合会使子波振幅谱“发胖”，不符合子波振幅谱的真实形状，给反褶积结果带来误差，进而不能提高地震记录分辨率和恢复地震波深层被衰减的能量。
技术实现要素：
针对上述问题，本发明的目的是提供一种能够提高地震分辨率和恢复地震波深层被衰减能量的基于高斯平滑的时频域反褶积方法。
为实现上述目的，本发明采取以下技术方案：一种基于高斯平滑的时频域反褶积方法，其特征在于，包括以下步骤：步骤1)：获取包括有多道非平稳地震记录的地震数据；步骤2)：基于高斯函数与广义s变换，得到改进广义s变换的公式；步骤3)：选取地震数据中某道非平稳地震记录，设定该道非平稳地震记录的均衡因子；步骤4)：设定改进广义s变换的参数，并根据设定的参数和均衡因子对该道非平稳地震记录进行改进广义s变换，得到均衡后非平稳地震记录的时频谱；步骤5)：对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱均进行高斯平滑处理，得到高斯平滑后非平稳地震记录的动态子波振幅谱；步骤6)：计算该道非平稳地震记录的反褶积因子，并将该反褶积因子与均衡后非平稳地震记录的时频谱相乘，得到该道非平稳地震记录的反射系数时频谱；步骤7)：对该道非平稳地震记录的反射系数时频谱进行s反变换，得到该道非平稳地震记录的时间域反射系数；步骤8)：重复步骤3)～7)直至求得地震数据中所有道非平稳地震记录的时间域反射系数，完成地震数据的时频域反褶积。
进一步地，所述步骤2)中基于高斯函数与广义s变换，得到改进广义s变换的公式，具体过程为：将高斯函数g(t)扩展为另一种形式：
其中，f为频率，t为时间，r为窗函数窗口的高度，σ为窗函数的方差；将上述高斯函数与广义s变换结合，得到改进广义s变换的公式：
其中，s(τ,f)为选取的非平稳地震记录x(t)的原始时频谱，τ为每一时间点，f0为信号主频率。
进一步地，所述步骤3)中选取地震数据中某道非平稳地震记录，设定该道非平稳地震记录的均衡因子，具体过程为：步骤3.1)：设定改进广义s变换的时频谱最大值范围；步骤3.2)：选取地震数据中某道非平稳地震记录，引入均衡因子均衡其原始时频谱，得到均衡后的该道非平稳地震记录xr(t)：
xr(t)＝x(t)*ra
其中，x(t)为地震数据中某道非平稳地震记录，ra为该道非平稳地震记录的均衡因子；步骤3.3)：根据设定的时频谱最大值范围设定该道非平稳地震记录的均衡因子，使均衡后非平稳地震记录的时频谱最大值落入设定的时频谱最大值范围内。
进一步地，所述步骤4)中设定改进广义s变换的参数，并根据设定的参数和均衡因子对该道非平稳地震记录进行改进广义s变换，得到均衡后非平稳地震记录的时频谱，具体过程为：改进广义s变换的参数包括窗函数的方差和窗口的高度，根据设定的窗函数的方差、窗口的高度和均衡因子，对均衡后的非平稳地震记录进行改进广义s变换，得到均衡后非平稳地震记录的时频谱sr(τ,f)：
进一步地，所述步骤5)中对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱均进行高斯平滑处理，得到高斯平滑后非平稳地震记录的动态子波振幅谱，具体过程为：步骤5.1)：对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱使用高斯函数：
其中，|sg(τ,f)|est为高斯平滑后非平稳地震记录的时频谱即整个时间点的振幅谱，amax(i)为均衡后非平稳地震记录的时频谱中时间点i的振幅谱峰值，fi为均衡后非平稳地震记录的时频谱中时间点i的频率值，fmax为均衡后非平稳地震记录的时频谱中的峰值频率值，s(i)为均衡后非平稳地震记录的时频谱中时间点i的振幅谱半宽信息；步骤5.2)：采用最小二乘法求解上述高斯函数，得到高斯平滑后非平稳地震记录的时频谱；步骤5.3)：根据高斯平滑后非平稳地震记录的时频谱，得到高斯平滑后非平稳地震记录的动态子波振幅谱。
进一步地，所述步骤5.3)中根据高斯平滑后非平稳地震记录的时频谱，得到高斯平滑后非平稳地震记录的动态子波振幅谱，具体过程为：基于rosa理论，均衡后非平稳地震记录的时频谱近似等同于其静态子波频谱ωr(f)、衰减函数a(τ,f)和反射系数r(τ,f)的乘积，只考虑均衡后非平稳地震记录的振幅谱，得到：
|sr(τ,f)|≈|ωr(f)||a(τ,f)||r(τ,f)|
其中，|sr(τ,f)|为均衡后非平稳地震记录的时频谱中每一时间点的振幅谱，|ωr(f)|为均衡后非平稳地震记录的静态子波振幅谱，|a(τ,f)|为均衡后非平稳地震记录的衰减函数振幅谱，|r(τ,f)|为均衡后非平稳地震记录的反射系数时频谱中每一个时间点的振幅谱；由于均衡后非平稳地震记录的动态子波振幅谱|ωα(τ,f)|等于其静态子波振幅谱与衰减函数振幅谱的乘积，即：
|ωα(τ,f)|＝|ωr(f)||a(τ,f)|
因此，得到均衡后非平稳地震记录的时频谱中每一时间点的振幅谱：
|sr(τ,f)|≈|ωα(τ,f)||r(τ,f)|
使用高斯函数对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱进行高斯平滑处理，能够消除噪声和反射系数的影响，因此，得到高斯平滑后非平稳地震记录的动态子波振幅谱|ωα(τ,f)|est：
|ωα(τ,f)|est≈|sg(τ,f)|est。
进一步地，所述步骤6)中计算该道非平稳地震记录的反褶积因子，并将该反褶积因子与均衡后非平稳地震记录的时频谱相乘，得到该道非平稳地震记录的反射系数时频谱，具体过程为：步骤6.1)：根据高斯平滑后非平稳地震记录的动态子波振幅谱，计算该道非平稳地震记录的反褶积因子振幅谱|r(τ,f)|：
其中，μ为调节因子；步骤6.2)：将该道非平稳地震记录的反褶积因子振幅谱进行希尔伯特变换，得到该道非平稳地震记录反褶积因子的相位谱θ(τ,f)：
θ(τ,f)＝hilbert{ln(|ωα(τ,f)|est+μamax)}
步骤6.3)：计算该道非平稳地震记录的反褶积因子r(τ,f)：
r(τ,f)＝|r(τ,f)|*exp(-iθ(τ,f))
步骤6.4)：在s域中将该道非平稳地震记录的反褶积因子和均衡后非平稳地震记录的时频谱相乘，得到该道非平稳地震记录的反射系数时频谱r(τ,f)est：
进一步地，所述步骤7)中非平稳地震记录的时间域反射系数为：
其中，r(t)为非平稳地震记录的时间域反射系数。
本发明由于采取以上技术方案，其具有以下优点：1、本发明通过对地震记录进行改进广义s变换和高斯平滑处理使得平滑算法更加稳定，能够更准确地求出地震记录的动态子波时频谱，进而得到反褶积因子，求出时间域反射系数，算法实现简单，不仅能有效的分辨薄层，提高地震分辨率，还能在保持信噪比的前提下恢复地震波深层被衰减的能量。2、本发明由于引入均衡因子，能够基本消除多道地震记录的时频谱最大值的不同对反褶积因子的影响，保证高斯平滑的稳定性，同时也缩小调节因子的不确定性带来的误差，可以广泛应用于石油地震勘探领域中。
附图说明
图1是将某道地震记录进行本发明时频域反褶积前后的对比图，其中，图1(a)为建立的随机反射系数模型图，图1(b)为平稳地震记录，图1(c)为非平稳地震记录，图1(d)为进行本发明时频域反褶积后的非平稳地震记录；
图2是图1中非平稳地震记录的时频谱分析图，其中，图2(a)为平稳地震记录的时频谱，图2(b)为原始非平稳地震记录的时频谱，图2(c)为高斯平滑后非平稳地震记录的时频谱，图2(d)为该道非平稳地震记录的反射系数时频谱；
图3是对图1中非平稳地震记录进行本发明高斯平滑前后的振幅谱对比图，其中，
图4是对中国西部三维工区时间深度为2.4～3.2s的地震数据进行本发明时频域反褶积前后的cdp(道集)对比图，其中，图4(a)为原始地震数据的cdp图，图4(b)为时频域反褶积后地震数据的cdp图；
图5是对中国西部三维工区时间深度为3.2～4s的地震数据进行本发明时频域反褶积前后的cdp对比图，其中，图5(a)为原始地震数据的cdp图，图5(b)为时频域反褶积后地震数据的cdp图；
图6是对图4和图5的地震数据中某道非平稳地震记录进行本发明时频域反褶积前后的对比图，其中，图6(a)为原始非平稳地震记录，图6(b)为时频域反褶积后的该道非平稳地震记录；
图7是图6中非平稳地震记录的时频谱分析图，其中，图7(a)为该道非平稳地震记录的原始时频谱，图7(b)为该道非平稳地震记录时频域反褶积后的时频谱；
图8是图6中非平稳地震记录的振幅谱对比图，其中，
图9是对中国西部三维工区时间深度为2～3.8s的地震数据进行本发明时频域反褶积前后的cdp对比图，其中，图9(a)为原始地震数据的cdp图，图9(b)为时频域反褶积后地震数据的cdp图；
图10是对图9的地震数据中某道非平稳地震记录进行本发明时频域反褶积前后的对比图，其中，图10(a)为原始非平稳地震记录，图10(b)为时频域反褶积后的该道非平稳地震记录；
图11是图10中非平稳地震记录的时频谱分析图，其中，图11(a)为该道非平稳地震记录的原始时频谱，图11(b)为该道非平稳地震记录时频域反褶积后的时频谱；
图12是图10中非平稳地震记录的振幅谱对比图，其中，
具体实施方式
以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。
本发明提供的基于高斯平滑的时频域反褶积方法，包括以下步骤：
1、获取地震数据，其中，地震数据包括多道非平稳地震记录。
2、基于高斯函数与广义s变换，得到改进广义s变换公式，具体过程为：
s变换是一种时频域分析方法，它吸收并发展了短时傅里叶变换和连续小波变换，其基本小波函数是由简谐波和高斯函数乘积构成，基本形态固定，不能根据实际情况调节窗口的大小，因此将高斯函数g(t)扩展为另一种形式：
其中，f为频率，t为时间，r为窗函数窗口的高度，σ为窗函数的方差，在实际资料处理时根据实际情况调节r和σ的值来控制时窗的大小，进而达到最佳效果。
因此，将上述高斯函数g(t)与广义s变换结合，得到改进广义s变换的公式：
其中，s(τ,f)为给定信号即选取的非平稳地震记录x(t)的原始时频谱，τ为每一时间点，f0为信号主频率。
3、选取地震数据中某道非平稳地震记录，设定该道非平稳地震记录的均衡因子，具体过程为：
3.1)设定改进广义s变换的时频谱最大值范围，其中，改进广义s变换的时频谱最大值范围可以根据实际情况进行设定，在此不做赘述。
3.2)选取地震数据中某道非平稳地震记录x(t)，引入均衡因子ra均衡其原始时频谱s(τ,f)，得到均衡后的该道非平稳地震记录xr(t)：
xr(t)＝x(t)*ra(3)
3.3)根据设定的时频谱最大值范围设定该道非平稳地震记录的均衡因子ra，使均衡后非平稳地震记录xr(t)的时频谱最大值落入设定的时频谱最大值范围(例如10-4～10-5)内。
4、设定改进广义s变换的参数，并根据设定的参数和均衡因子对该道非平稳地震记录进行改进广义s变换，得到均衡后非平稳地震记录的时频谱，其中，改进广义s变换的参数包括窗函数的方差σ和窗口的高度r，可以根据实际情况以达到地震记录时频谱的最高分辨率为基准进行设定，在此不做赘述。
基于公式(2)，根据设定的窗函数的方差σ、窗口的高度r和均衡因子ra，对均衡后的非平稳地震记录xr(t)进行改进广义s变换得到其时频谱sr(τ,f)：
5、对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱均进行高斯平滑处理，得到高斯平滑后非平稳地震记录的动态子波振幅谱，具体过程为：
5.1)对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱|sr(τ,f)|使用高斯函数：
其中，|sg(τ,f)|est为高斯平滑后非平稳地震记录的时频谱即整个时间点的振幅谱，amax(i)为均衡后非平稳地震记录的时频谱中时间点i的振幅谱峰值，fi为均衡后非平稳地震记录的时频谱中时间点i的频率值，fmax为均衡后非平稳地震记录的时频谱中的峰值频率值，s(i)为均衡后非平稳地震记录的时频谱中时间点i的振幅谱半宽信息。
5.2)采用最小二乘法求解上述高斯函数，得到高斯平滑后非平稳地震记录的时频谱。
将公式(5)两边同时取对数得到：
令ln|sg(τ,f)|est＝z(i)，
5.3)根据高斯平滑后非平稳地震记录的时频谱，得到高斯平滑后非平稳地震记录的动态子波振幅谱。
基于rosa理论，均衡后非平稳地震记录的时频谱sr(τ,f)近似等同于其静态子波频谱ωr(f)、衰减函数a(τ,f)和反射系数r(τ,f)的乘积，只考虑均衡后非平稳地震记录xr(t)的振幅谱，可以得到：
|sr(τ,f)|≈|ωr(f)||a(τ,f)||r(τ,f)|(7)
其中，|sr(τ,f)|为均衡后非平稳地震记录的时频谱中每一时间点的振幅谱，|ωr(f)|为均衡后非平稳地震记录的静态子波振幅谱，|a(τ,f)|为均衡后非平稳地震记录的衰减函数振幅谱，|r(τ,f)|为均衡后非平稳地震记录的反射系数时频谱中每一个时间点的振幅谱。
由于均衡后非平稳地震记录的动态子波振幅谱|ωα(τ,f)|等于其静态子波振幅谱|ωr(f)|与衰减函数振幅谱|a(τ,f)|的乘积，即：
|ωα(τ,f)|＝|ωr(f)||a(τ,f)|(8)
因此，可以得到均衡后非平稳地震记录的时频谱中每一时间点的振幅谱|sr(τ,f)|：
|sr(τ,f)|≈|ωα(τ,f)||r(τ,f)|(9)
上述公式(9)中，均衡后非平稳地震记录的时频谱中每一时间点的振幅谱主要趋势是由动态子波振幅谱引起的，而反射系数影响均衡后非平稳地震记录的时频谱中每一时间点振幅谱的细节部分，因此，使用更加接近动态子波振幅谱的高斯函数对均衡后非平稳地震记录的时频谱中每一时间点的振幅谱进行高斯平滑处理，能够消除噪声和反射系数的影响，进而得到高斯平滑后非平稳地震记录中每一时间点的动态子波振幅谱。
因此，基于上述rosa理论，可以认为高斯平滑后非平稳地震记录的时频谱消除了反射系数的影响，即得到高斯平滑后非平稳地震记录的动态子波振幅谱|ωα(τ,f)|est：
|ωα(τ,f)|est≈|sg(τ,f)|est(10)
6、计算该道非平稳地震记录的反褶积因子振幅谱，并将该反褶积因子振幅谱与均衡后非平稳地震记录的时频谱相乘，得到该道非平稳地震记录的反射系数时频谱，具体过程为：
6.1)根据高斯平滑后非平稳地震记录的动态子波振幅谱，计算该道非平稳地震记录的反褶积因子振幅谱|r(τ,f)|：
其中，r(τ,f)为反褶积因子，为防止分母出现零值引入参数μamax，其中，μ为调节因子，是一个很小的值，可以根据实际情况进行设定，在此不做赘述。
6.2)将该道非平稳地震记录的反褶积因子振幅谱进行希尔伯特变换，得到该道非平稳地震记录反褶积因子的相位谱θ(τ,f)：
θ(τ,f)＝hilbert{ln(|ωα(τ,f)|est+μamax)}(12)
6.3)计算该道非平稳地震记录的反褶积因子r(τ,f)：
r(τ,f)＝|r(τ,f)|*exp(-iθ(τ,f))(13)
6.4)在s域中将该道非平稳地震记录的反褶积因子和均衡后非平稳地震记录的时频谱相乘得到该道非平稳地震记录的反射系数时频谱r(τ,f)est：
7、对该道非平稳地震记录的反射系数时频谱进行s反变换得到该道非平稳地震记录的时间域反射系数r(t)：
通过得到该道非平稳地震记录的时间域反射系数，能够消除该道非平稳地震记录子波的影响并移除衰减函数的效应，达到了压缩子波进而提高地震分辨率，恢复地震波深层能量衰减的作用。
8、重复步骤3～7求得地震数据中所有道非平稳地震记录的时间域反射系数，完成地震数据的时频域反褶积。
如图1～3所示，下面通过建立一加噪声的随机反射系数模型，对本发明基于高斯平滑的时频域反褶积方法的合理性和实用性进行验证：
将该随机反射系数模型与主频为40hz的雷克子波进行褶积得到平稳地震记录如图1(b)所示，实际中地震波的传播通常是伴随着能量衰减的，因此在频率域对平稳地震记录引入衰减函数，再通过傅里叶反变换得到非平稳地震记录如图1(c)所示，其中，在0～0.25s时品质因子q取值70，0.25～1s时品质因子q取值50。将如图1(c)所示中的非平稳地震记录进行改进广义s变换，得到非平稳地震记录的时频谱如图2(b)所示，对非平稳地震记录的时频谱进行本发明的高斯平滑处理得到高斯平滑后非平稳地震记录的时频谱即每一时间点的动态子波振幅谱如图2(c)所示，利用该动态子波振幅谱求取反褶积因子，将其与非平稳地震记录的时频谱相乘得到该道非平稳地震记录的反射系数时频谱如图2(d)所示，再对该反射系数时频谱进行s反变换得到该道非平稳地震记录的时间域反射系数如图1(d)所示，实现该道非平稳地震记录的时频域反褶积。将图1(c)和图1(d)比较可以发现将非平稳地震记录进行时频域反褶积后被衰减的地震记录能量得到了恢复，浅层薄层刻画更明显，同时如图1(d)所示的时间域反射系数与如图1(a)所示的加噪声的随机反射系数模型相匹配，说明模型论证成功。如图3所示，取图1(c)中非平稳地震记录的某一时间点振幅谱进行分析，该道非平稳地震记录的振幅谱明显受反射系数影响很大，将高斯平滑后整个时间点的振幅谱与其动态子波振幅谱进行对比，两者波形基本吻合，消除了时间域反射系数对动态子波细节的影响，这表明本发明高斯平滑处理提取动态子波振幅谱的准确性。
如图4～8所示，下面以中国西部a地区某实际三维工区深层叠加地震数据为具体实施例对本发明时频域反褶积方法的效果进行说明：
如图4所示，调节因子μ取0.08，设定窗函数的方差σ＝0.6，窗口的高度r＝2.1，如图4(a)所示，可以明显看出原始地震数据受地层的吸收作用深层地震波能量被衰减严重，箭头位置同相轴连续性差，将如图4(b)所示的经本发明时频域反褶积处理后的地震数据与如图4(a)所示的原始地震数据对比，可以明显看出衰减的能量得到了很好的恢复，同相轴更加清晰，连续性更好。如图5所示，可以看出如图5(b)所示的经本发明时频域反褶积处理后的地震数据与如图5(a)所示的原始地震数据对比，在没有降低信噪比的情况下，提升了整个地震数据能量，提高了地震记录纵向分辨率，箭头标记的薄层也能很明显的识别出来，达到了高分辨效果。如图6和图7所示，可以看出对抽取的图4和图5中某道非平稳地震记录进行本发明的时频域反褶积处理后，该道地震记录的能量得到了增强，如图8所示，可以看出地震记录的频带也得到了拓宽，整个频带能量得到了恢复，提高了地震记录的分辨率。
如图9～12所示，下面以中国南海某工区海上地震数据为具体实施例对本发明时频域反褶积方法的效果进行说明：
如图9所示，μ取0.04，设定窗函数的方差σ＝0.5，窗口的高度r＝2.1，可以明显看出如图9(b)所示的经本发明时频域反褶积处理后的地震数据与如图9(a)所示的原始地震数据对比，地震数据的能量得到了补偿，同相轴横向连续性增强，能量聚集性更好，椭圆圈出的位置同相轴更加清晰，箭头位置的反射层更加清楚，整个地震数据的分辨率均得到了提高。如图10所示，可以看出经本发明时频域反褶积处理后地震记录的动态子波得到了压缩，衰减部分的能量获得恢复。如图11所示，可以从时频谱上更直观地看出经本发明时频域反褶积处理前后，地震记录深层被衰减的能量得到了恢复。如图12所示，可以看出经本发明时频域反褶积处理后地震记录的振幅谱得到了拓宽，说明地震记录的动态子波得到了压缩，提高了地震记录的分辨率。
上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。
- 还没有人留言评论。精彩留言会获得点赞！