一种从叠后地震数据采集脚印衰减的方法

摘要

本发明提供一种从叠后地震数据采集脚印衰减的方法，所述方法包括：将叠后地震数据划分为多个时窗，对每个时窗执行以下步骤，对时窗内的叠后地震数据进行截断奇异值分解，得到分解后的地震数据；计算分解后的地震数据在不同秩值下的平均功率谱；对叠后地震数据进行多项式拟合，利用多项式拟合的结果合成地震纯波模型数据；计算地震纯波模型数据的平均功率谱；利用地震纯波模型数据的平均功率谱和分解后的地震数据在不同秩值下的平均功率谱，计算分解后的地震数据在不同秩值下的信噪比谱；将计算的信噪比谱之中的最大信噪比谱所对应的秩值作为地震数据重构的最佳秩值；利用最佳秩值对叠后地震数据进行重构，得到重构后的地震数据。

说明书

技术领域

本发明属于石油勘探地震资料处理领域，特别涉及一种从叠后地震数据采集脚印衰减的方法。

背景技术

随着岩性油气藏和隐蔽油气藏勘探对地震成像分辨率要求的提高，对地震资料的细节刻画要求也更高了，不仅不规则的采集和空间采样不足会产生采集脚印，规则的三维观测系统、规则的空间采样同样也会产生采集脚印，它们都影响地质目标的高精度、高质量地震成像，影响AVO属性的分析和速度分析，这种“采集脚印”又称“采集痕迹”，一般表现在三维叠加数据上。

目前针对叠后地震数据采集脚印衰减的常用方法有：道混合法、F-K滤波法、三维真振幅倾角滤波法、针对时间切片的F-KX-KY法、自适应滤波法等，虽然这些方法都能不同程度地衰减叠后地震数据中的采集脚印，但处理后的时间切片上往往会以模糊横向的振幅变化为代价。例如，F-KX-KY陷波滤波法利用采集脚印在时间切片上的X、Y方向上存在等间隔的现象，与时间方向的50Hz工业电干扰在频率域表现为一样的脉冲，采集脚印会在频率切片的KX、KY域也表现为脉冲，则可在F-KX-KY域用陷波滤波来消除采集脚印，但由于采集脚印在X、Y方向的频率切片上存在数个等间隔脉冲，需要一定的频率计算宽度和个数，由于无法准确计算出其间隔宽度，为了衰减掉采集脚印，就会给出较大的计算宽度和个数，从而造成部分有效信号被当作噪声脉冲衰减掉，使得采集脚印衰减后的时间切片横向振幅较模糊。

对叠后地震数据进行采集脚印衰减，获得高信噪比地震数据是高精度地震解释的关键，但是目前常用的采集脚印衰减方法在处理地震数据过程中无法对其信噪比进行衡量和质控，使得对地震资料衰减程度不同，其噪声压制效果不同，往往得不到最佳信噪比的地震数据。例如，常用的SVD截断奇异值分解压制采集脚印法利用分解后对角线上的元素，将其作为奇异值，并根据大的奇异值代表着有效信号，小的奇异值代表着噪声，对大的奇异值进行重构从而获得采集脚印衰减后的地震数据，但是该方法没有一个有效判断奇异值大小的根据，以及大奇异值的个数，使得不同个数的奇异值重构得到不同信噪比的地震数据，无法对处理过程中的地震数据进行有效衡量和质控。

发明内容

本发明的目的为利用截断奇异值分解重构衰减采集脚印技术获得具有最佳信噪比的地震数据，使重构后的地震数据采集脚印与有效信号得到最佳分离，且时间切片上振幅较清晰，较大程度上改善地震数据去噪效果。

本发明提供一种从叠后地震数据采集脚印衰减的方法，所述方法包括：将叠后地震数据划分为多个时窗，对每个时窗执行以下步骤，对时窗内的叠后地震数据进行截断奇异值分解，得到分解后的地震数据；计算分解后的地震数据在不同秩值下的平均功率谱；对叠后地震数据进行多项式拟合，利用多项式拟合的结果合成地震纯波模型数据；计算地震纯波模型数据的平均功率谱；利用地震纯波模型数据的平均功率谱和分解后的地震数据在不同秩值下的平均功率谱，计算分解后的地震数据在不同秩值下的信噪比谱；将计算的信噪比谱之中的最大信噪比谱所对应的秩值作为地震数据重构的最佳秩值；利用最佳秩值对叠后地震数据进行重构，得到重构后的地震数据。

此外，所述相邻两个时窗的移动间隔为半个时窗。

此外，计算分解后的地震数据在不同秩值下的信噪比谱的公式为，

$>{\mathrm{SNR}}_r=\frac{\overline{G_s\left(f\right)}}{\overline{{G_n}_r\left(f\right)}},$>

其中，SNR_r表示分解后的地震数据在秩值r下的信噪比谱，为地震纯波模型数据的平均功率谱，为分解后的地震数据在秩值r下的噪声平均功率谱，为分解后的地震数据在秩值r下的平均功率谱。

根据本发明的从叠后地震数据采集脚印衰减的方法，利用信噪比谱质控下求取的地震数据重构的最佳秩值，对叠后地震数据进行重构，使得重构后的地震数据信噪比更高，时间切片的振幅清晰度也更高，为精细地质解释和储层描述提供更加可靠的基础地震数据。

附图说明

通过下面结合附图进行的详细描述，本发明的上述和其它目的、特点和优点将会变得更加清楚，其中：

图1示出根据本发明实施例的从叠后地震数据采集脚印衰减方法的流程图；

图2是采集脚印衰减前主测线方向剖面图；

图3是根据本发明实施例的采集脚印衰减后主测线方向剖面图；

图4是示出采集脚印衰减前T=840秒时间切片示意图；

图5是示出根据本发明实施例的采集脚印衰减后T=840秒时间切片示意图。

具体实施方式

下面，将参照附图详细描述本发明的实施例。

提供参照附图的以下描述以帮助对由权利要求及其等同物限定的本发明的实施例的全面理解。包括各种特定细节以帮助理解，但这些细节仅被视为是示例性的。因此，本领域的普通技术人员将认识到在不脱离本发明的范围和精神的情况下，可对描述于此的实施例进行各种改变和修改。此外，为了清楚和简洁，省略对公知的功能和结构的描述。

本发明的从叠后地震数据采集脚印衰减的方法包括：将叠后地震数据划分为多个时窗，对每个时窗执行以下步骤，（1）对时窗内的叠后地震数据进行截断奇异值分解，得到分解后的地震数据；（2）计算分解后的地震数据在不同秩值下的平均功率谱；（3）对叠后地震数据进行多项式拟合，利用多项式拟合的结果合成地震纯波模型数据；（4）计算地震纯波模型数据的平均功率谱；（5）利用地震纯波模型数据的平均功率谱和分解后的地震数据在不同秩值下的平均功率谱，计算分解后的地震数据在不同秩值下的信噪比谱；（6）将计算的信噪比谱之中的最大信噪比谱所对应的秩值作为地震数据重构的最佳秩值；（7）利用最佳秩值对叠后地震数据进行重构，得到重构后的地震数据。

具体地讲，将叠后地震数据划分为多个时窗，依次对每个时窗进行采集脚印衰减处理，所述相邻两个时窗的移动间隔为半个时窗，对每个时窗的执行步骤为：

图1示出根据本发明实施例的从叠后地震数据采集脚印衰减方法的流程图。

在步骤101，对时窗内的地震数据进行截断奇异值分解，得到分解后的地震数据。

野外采集到的叠后地震数据通常较大，运算时将叠后地震数据划分为多个时窗，分时窗进行运算，依次对不同时窗内的叠后地震数据进行采集脚印衰减处理。

例如，叠后地震数据X可被表示为，

X={X_nm|n∈[1,N]，m∈[1,M]}           （1-1）

式中，x_nm为叠后地震数据，n为地震道的数量，m为采样点的数量，且N、M为大于零的整数。

对叠后地震数据X进行截断奇异值分解表示为，

$>X_r=\underset{i=1}{\overset{r}{\Sigma }}{\sigma }_i{u}_i{v}_i^T---(1-2)$>

其中，X_r为分解后的地震数据，r为叠后地震数据X的秩，上标T表示转置，σ_i为叠后地震数据X的第i个奇异值，u_i为矩阵XX^T的第i个特征向量，v_i是矩阵X^TX的第i个特征向量，i为叠后地震数据X的秩值，1≤i<r，i为整数。

从（1-2）式的左端到右端，是对叠后地震数据X的截断奇异值分解过程，当叠后地震数据X的秩值为r时，对（1-2）式进一步简化，得到分解后的地震数据为，

$>X_r=\underset{i=1}{\overset{r}{\Sigma }}{I}_i---(1-3)$>

式中，是由奇异值和本征图像构成的参数变量，是叠后地震数据X的第i个本征图像，由于从（1-2）式的左端到右端是一个分解过程，那么从右端到左端就是一个重构过程，则由SVD理论可知，分解后地震数据X_r的范数为，

$>\left|\right|X_r\left|\right|={\left(\underset{i=1}{\overset{r}{\Sigma }}{\sigma }_i^2\right)}^{\frac{1}{2}}---(1-4)$>

分解后地震数据X_r的总能量可以表示为，

$>{\left|\right|X_r\left|\right|}^2\underset{i=1}{\overset{r}{\Sigma }}{\sigma }_i^2---(1-5)$>

由式（1-5）可以看出叠后地震数据X经截断奇异值分解后，叠后地震数据X的总能量可以用奇异值平方和来表示，那么对地震数据进行截断奇异值分解，地震数据的总能量也可以用它的奇异值平方和表示。由于有效信号是有规律的，而时间切片上的采集脚印在地震剖面上表现为杂乱无章的随机噪声，因此认为大奇异值对应的本征图像主要代表有效信号，小奇异值对应的本征图像主要代表噪音。而分解后的地震数据的奇异值是按递减顺序排列的，需要找到地震数据重构的最佳秩值k，1≤k<r，k为整数，选用前k个奇异值对应的本征图像对地震数据进行重构，这样就可以有效地衰减采集脚印。但要实现有效信号与采集脚印和谐地进行分离，必须准确知道地震数据重构所需的最佳秩值k，因此还必须进行以下步骤的操作。这里，SVD理论为本领域公知的常识，本发明省略了此部分的描述。

在步骤102，计算分解后地震数据在不同秩值下的平均功率谱。

具体地讲，分解后的地震数据为当秩值r取值不同时，即有不同个奇异值对应的本征图像进行叠后地震数据重构。例如：当取r=3时，X₃=I₁+I₂+I₃，则有3个奇异值对应的本征图像进行叠后地震数据重构；当取r=5时，X₅=I₁+I₂+I₃+I₄+I₅，则有5个奇异值对应的本征图像进行地震数据重构，在不同的秩值r下，重构的地震数据也不同。

根据本发明实施例，利用式（1-3）中分解后的地震数据通过平均功率谱的计算公式，计算分解后地震数据在秩值r下的平均功率谱这里，平均功率谱的计算公式为本领域的公知常识，在本发明中省略了对此部分的描述。

在步骤103，对叠后地震数据进行多项式拟合，利用多项式拟合的结果合成地震纯波模型数据。这里，可以利用现有的各种利用多项式拟合的结果合成地震纯波模型数据的技术实现步骤103，下面给出一个合成地震纯波模型数据的示例。

首先，计算叠后地震数据X的时间多项式。例如，利用正交多项式代替一般多项式来计算时间多项式，这样可以利用独立扫描多项式系数的方法减少计算量。

根据现有技术可以得到叠后地震数据X的时间多项式，时间多项式的计算公式如下，

$>T\left(x\right)=c_0+\underset{j=1}{\overset{3}{\Sigma }}{c}_j{p}_j\left(x\right)---(1-6)$>

式中，T(x)为叠后地震数据X的时间多项式，c_j为拟合系数，c₀为拟合系数初始值，p_j(x)为正交多项式，j代表阶数，且j为整数。

然后，计算叠后地震数据X的振幅多项式。振幅多项式的计算公式如下，

$>A\left(x\right)=\underset{j=0}{\overset{3}{\Sigma }}{c}_j{p}_j\left(x\right)---(1-7)$>

式中，A(x)为叠后地震数据X的振幅多项式，c_j为拟合系数，p_j(x)为正交多项式，j代表阶数，且j为整数。

最后，根据叠后地震数据X的时间多项式和振幅多项式，形成地震纯波模型数据。例如，将时间多项式确定的同一时窗内的n道地震数据波形沿拟合出的同相轴方向相加，并对相加结果缩放，使其均方根振幅归一。波形相加公式如下，

$>A_l=\underset{n=-N}{\overset{N}{\Sigma }}S(n,T_n+l-L-1),l\in 2L+1---(1-8)$>

式中，A_l为地震纯波模型数据的期望波形，S(n,T_n+l-L-1)为第n道上时窗内地震道的样点值，n为地震道数，这里-N≤n≤N，N为整数，l为时间方向上时窗大小，且l为大于零的整数，L为时间方向上相邻两个时窗移动的大小，T_n为第n道的时窗中点时间，由时间多项式确定。

公式（1-8）为地震纯波模型数据的期望波形，该期望波形在时间方向是准确的，即每个同相轴波形在对应的时间点是准确的，但是每个期望波形的振幅是一样的，并不是实际的振幅值（实际地震数据的振幅应有强有弱），所以，合成在地震纯波模型数据时，还需考虑到振幅对地震纯波模型数据的影响。

根据本发明实施例，利用每个时窗的时间多项式（1-6）、振幅多项式（1-7）、和地震纯波模型数据期望波形（1-8），来合成地震纯波模型数据，即任意一道地震数据按照其在时窗中的位置计算时间多项式、振幅多项式，将振幅与期望波形的乘积放到计算出的时间位置上就可以合成地震纯波模型数据。

在步骤104，计算地震纯波模型数据的平均功率谱。

如果按照现有技术中的方式计算有效信号的平均功率谱对于非水平界面和低信噪比资料的计算存在很大的误差，因为实际的地震数据中地震道的相关性会随着相离距离的增大相关性减弱，尤其是低信噪比的地震道相关性受噪声影响更大。

本发明对时窗内地震数据按照步骤103进行多项式拟合，利用多项式拟合的结果合成地震纯波模型数据，再通过平均功率谱的计算公式，计算地震纯波模型数据的平均功率谱把地震纯波模型数据的平均功率谱作为时窗内地震数据有效信号的平均功率谱，避免非水平界面和低信噪比资料中噪声对求取信噪比谱的影响。这里，平均功率谱的计算公式为本领域的公知常识，在本发明中省略了对此部分的描述。

在步骤105，利用地震纯波模型数据的平均功率谱和分解后的地震数据在不同秩值下的平均功率谱，计算分解后的地震数据在不同秩值下的信噪比谱。

根据步骤102中得到的分解后地震数据在秩值r下的平均功率谱与步骤104中得到的地震纯波模型数据的平均功率谱可得到分解后的地震数据在秩值r下的噪声平均功率谱

$>\overline{{G_n}_r\left(f\right)}=\overline{{G_x}_r\left(f\right)}-\overline{G_s\left(f\right)}---(1-9)$>

式中，为地震纯波模型数据的平均功率谱，为分解后的地震数据在秩值r下的平均功率谱。

然后，再根据信噪比谱的计算公式，计算分解后的地震数据在不同秩值下的信噪比谱SNR_r：

$>{\mathrm{SNR}}_r=\frac{\overline{G_s\left(f\right)}}{\overline{{G_n}_r\left(f\right)}}---(1-16)$>

其中，为地震纯波模型数据的平均功率谱，为分解后的地震数据在秩值r下的噪声平均功率谱。

根据本发明实施例，在计算分解后的地震数据在秩值r下的信噪比谱时，将地震纯波模型数据的平均功率谱作为有效信号的功率谱而地震纯波模型数据的求取通过步骤103的多项式拟合的结果得到。

在步骤106，将计算的信噪比谱之中的最大信噪比谱所对应的秩值作为地震数据重构的最佳秩值。

具体地讲，当秩值r取值不同时，分解后的地震数据的组成也不同，则分解后的地震数据在秩值r下的平均功率谱也不同，所以分解后的地震数据在秩值r下的信噪比谱SNRr的值就不同。在计算的信噪比谱SNR_r中找到SNR_r值最大的一个，将最大信噪比谱所对应的秩值作为地震数据重构的最佳秩值。

在步骤107，利用最佳秩值对叠后地震数据进行重构，得到重构后的地震数据。

具体地讲，地震数据重构的最佳秩值为k，则重构后的地震数据为由于大奇异值对应的本征图像主要代表有效信号，小奇异值对应的本征图像主要代表噪音，且分解后的地震数据的奇异值是按递减顺序排列的，因此，选用前k个奇异值对应的本征图像来进行地震数据的重构，这样就可以有效地衰减采集脚印，实现有效信号与采集脚印和谐地进行分离。

一个时窗内的地震数据处理完成后，将移动至下一个时窗，按照步骤101至步骤107的计算过程，对下一个时窗内的地震数据进行重构，待所有时窗计算完成后，获取多个时窗内的重构后的地震数据，得到完整的采集脚印衰减后的最终地震数据。

图2是采集脚印衰减前主测线方向剖面图，横坐标为InLine（主测线）方向测线号，纵坐标为Time（时间）方向，图3是根据本发明实施例的采集脚印衰减后主测线方向剖面图，横坐标为InLine（主测线）方向测线号，纵坐标为Time（时间）方向，图2和图3为三维数据体重在InLine（主测线）方向上的剖面图。

本发明的实施例来源于国内某地区地震勘探资料，从该地区地震勘探资料的图2剖面上看，杂乱无章的随机噪声随处可见，严重影响了地震勘探资料的信噪比，从比较示意图可以看出图3比图2中剖面的浅、中、深层信噪比大大提高、同相轴连续性更好，而且被噪声淹没影响的有效信号更加突出。

图4是示出采集脚印衰减前T=840秒时间切片示意图，横坐标为InLine（主测线）方向测线号，纵坐标为CrossLine（横测线）方向测线号，图5是示出根据本发明实施例的采集脚印衰减后T=840秒时间切片示意图，横坐标为InLine（主测线）方向测线号，纵坐标为CrossLine（横测线）方向测线号。

从图4所示T=840秒的时间切片示意图可知，剖面上的随机噪声实际上是由时间切片上表现为条带状的采集脚印造成的，这些采集脚印的存在给后续的精细解释带来很大困扰，从图5与图4中T=840秒的时间切片对比分析知，经本发明的从叠后地震数据采集脚印衰减的方法处理后，采集脚印衰减后的条带状采集脚印基本被衰减掉，突出了有效信号的能量，而且该技术并没有让切片的振幅模糊，保持了有效信号原有能量的精度，细节刻画清晰可见，为后续的精细解释提供信噪比更高，精度更好的剖面。

根据本发明的从叠后地震数据采集脚印衰减的方法，是在运用截断奇异值分解重构衰减采集脚印技术时，对时窗内的叠后地震数据利用多项式拟合求取地震纯波模型数据，进而计算地震纯波模型数据的平均功率谱，将其作为有效信号的平均功率谱，通过信噪比谱计算公式计算时窗内不同秩值下分解后的地震数据的信噪比谱，将信噪比谱最大值时的秩值作为地震数据重构的最佳秩值，从而使得重构后的地震数据采集脚印和有效信号得到最佳的分离，从实际地震勘探资料应用可知，该技术取得了不错的应用效果，算法稳定，去噪效果较好，应用前景较为广阔。

尽管已经参照其示例性实施例具体显示和描述了本发明，但是本领域的技术人员应该理解，在不脱离权利要求所限定的本发明的精神和范围的情况下，可以对其进行形式和细节上的各种改变。