基于逆二分递推奇异值分解的地震资料高分辨率处理方法与流程

本发明涉及地震勘探领域，具体涉及基于逆二分递推奇异值分解的地震资料高分辨率处理方法。

背景技术：

在地震勘探中，提高地震信号的分辨率是数据采集与处理中十分重要的问题。因为地震信号分辨率是地震勘探工作中获取地层细节信息的关键因素，对研究薄层或小的地质体具有十分重要的意义，为此很多地球物理学家提出并发展了几类地震资料高分辨率处理的方法：(1)谱白化，它通过展宽振幅谱来提高信号分辨率，且不改变子波的相位谱，是一种“纯振幅”的滤波过程；(2)反q滤波，一种补偿大地吸收衰减效应的技术，它不仅可以补偿振幅衰减和频率损失，而且还可以改善记录的相位特性,从而改善同相轴的连续性，提高弱反射波的能量和地震资料的信噪比、分辨率；(3)多尺度联合分析方法，该方法常常利用测井资料、井间地震、vsp等地球物理技术对地下同一目标地质体进行不同尺度性质的反映，通过他们之间的联合作用提高地震资料分辨率；

(4)反褶积，通过假设地震子波为最小相位、反射系数为高斯白噪分布，采用地震记录的自相关代替子波的自相关，并在此基础上使用基于二阶统计量的维纳滤波来实现子波估计和反褶积。

上述方法对提高地震资料分辨率均有不错的效果，然而这些方法或难以保持地震数据振幅相对关系、或高度依赖品质因子q值求取、或需要特殊的井中资料，只能在各自的适用范围内取得较好的补偿结果。

svd分解是利用地震信号之间相关性较大的这一特点，按照能量分布关系，通过伸缩旋转对地震资料进行分解的一种方法。svd方法能将信号分解为一系列反映信号细部特征的细节信号和反映信号主体骨架的近似信号的组合。多分辨率奇异值分解算法(multi-resolutionsingularvaluedecomposition,mrsvd)是将矩阵二分递推结构原理和svd方法结合，通过伸缩旋转对信号逐步进行多尺度细化的一种自适应的时频分析方法。该方法不存在确定矩阵行、列数的难题，矩阵的结构简单，但是利用递推分解方式和这种简单的二分矩阵结构相结合，却能实现一种对信号逐次剥离的多层次分解，很好地把信号中微弱的细节信号和主体信号多层次体现出来，从而有利于提取其中隐含的信号特征。目前已经成功的应用于信号识别、信号恢复和消噪、机械故障诊断等领域。

为了便于理解，对mrsvd算法原理进行说明。

mrsvd分解过程：对于离散地震信号x＝(x1,x2,x3,…,xn)，用此信号构造一个行数为2的hankel矩阵，

对此矩阵进行svd处理，得到

h＝usv^t(2)

式中正交矩阵u＝(u1,u2)，u∈r^2×2,正交矩阵v＝(υ1,υ2,…,υ(n-1))，v∈r^(n-1)×(n-1),对角矩阵s＝(diag(σa,σd),o)，s∈r^2×(n-1),σa＜＜σd。将式(2)改写成用列向量ui和υi表示形式：

式中，ui∈r^2×1，υi∈r^(n-1)×1，i＝1,2。令则ha∈r^2×(n-1)，它对应的是大奇异值，反映信号的主体特征，称其为近似矩阵；hd∈r^2×(n-1)，它对应的是小奇异值，反映信号的细节特征，称其为细节矩阵。

第一次svd得到的近似信号a1和细节信号d1分别从矩阵ha、hd获得。以细节信号d1＝(d1,d2,…,dn)的求取为例来说明其获取的过程，细节矩阵hd是一个两行的向量

其中，u2,1，u2,2为列向量u2的第1、2个坐标。

如(5)式，设ld1和ld2是细节矩阵hd两个行向量的子向量，且分别代表各自行向量中的d2,d3，…，dn-1，但ld1≠ld2。例如d2在ld1中的值为σd1u2,1υ2,2，而在ld2中的值为σd1u2,2υ2,1，这两个值显然不相等。所以为了获得信息完全的近似信号，将ld1和ld2求平均，再利用这种平均值作为细节信号d1中相应的数据。因此，d1最终可表示为如下形式：

d＝(d1，(ld1+ld2)/2,dn)(6)

同理，可得到近似信号a1。这样就使用mrsvd方法得到了第1次分解的结果d1和a1，细节信号d1对应的是小奇异值σd1，反映的是信号的细节特征。近似信号a1对应的是大奇异值σa1，反映的是信号的主体特征。接着利用a1构造(1)式所示的矩阵，并进行同样的处理，可得两个分量信号d2和a2,如此逐层进行分解，最终将原信号分解为一系列的细节信号和近似信号。

如图2所示，发明人对mrsvd分解过程中得到的近似信号的振幅谱进行研究，发现随着分解次数不断地增加，原信号的高频部分以细节信号的形式不断地被分解出来，mrsvd实质上是不断的分解出信号的高频分量。

因此，发明人研究发现mrsvd的重构过程就是把细节信号和近似信号逐层叠加的过程，即将第m层的细节信号dm与近似信号am叠加构建第(m-1)层的近似信号am-1，然后近似信号am-1再和细节信号dm-1叠加构建第(m-2)层的近似信号am-2，如此逐层进行，可得到原信号x的重构公式为：

式中，m表示总的分解层数。

技术实现要素：

本发明目的在于建立一种基于逆二分递推奇异值分解(imrsvd)的地震资料自适应高分辨率处理方法，能够将检测到的地震信号的缺失高频部分逆推出来，从而叠加获得高分辨率的地震信号。

为了实现上述技术效果，本发明采用以下技术方案：

基于逆二分递推奇异值分解的地震资料高分辨率处理方法，包括以下步骤：

步骤1：获取单道地震信号x；

步骤2：采用mrsvd算法将地震信号分解，然后利用得到的细节奇异值逐层逆向递推获得新的细节信号和近似信号；

步骤3：将新的细节信号逐次叠加到原始信号上，补偿地震信号的高频部分，从而得到高分辨率的地震信号，采用的公式为：

式中，x表示原信号，a′i表示第i次高频补偿的结果，g表示总的逆向递推次数，d′i为细节信号。

作为一种优选技术方案，通过修正方差模对总的逆向递推次数进行控制，修正方差模计算公式为：

其中，a′i(t)表示第i次高频补偿的结果，t为时间，n为信号的长度，a为常数。对每次高频补偿的信号a'1，a'2,...,a'(g-1),a'g计算其修正方差模为v1,v2,...,v(g-1),vg，如果v(g-6)≈v(g-3)≈vg，即修正方差模收敛并达到最大值，这时总的逆向递推次数g得以确定，且a'g为最终得到的高分辨率地震信号。

作为一种优选技术方案，上述步骤2中利用mrsvd分解得到的细节奇异值σd1,σd2,…,σdm，通过拟合函数逆向递推出新的细节奇异值σ′di(i＝1,2,…)，再通过细节奇异值得到对应的细节信号d′i，拟合指数函数为：

其中，j表示mrsvd的分解次数；an代表多项式的系数；k是一个正数，常常小于3；n是多项式的阶数，在最小平方意义下使f(j)逼近已知的细节奇异值，求得k和多项式的系数。

作为一种优选技术方案，mrsvd正向分解级数通过以下公式得到：

ej＝∑|aj-1-aj|²/∑|aj-1|²,(j＝1,…,m)

j表示mrsvd正向分解第j层，当ej≤10^-6时，循环分解结束,m是mrsvd正向分解总层数；aj-1和aj分别为第j-1层和第j层分解得到的近似信号。

作为一种优选技术方案，利用新的细节奇异值构建细节矩阵从而得到对应的细节信号。

本发明的有益效果在于：

本发明基于原始信号的特征通过逆向递推地震信号的高频部分建立了逆二分递推奇异值分解(imrsvd)，该算法的主要思想：为了恢复大地滤波引起的地震信号缺失的高频，我们通过mrsvd得到的细节信号的特征逆向递推第1个细节信号，即第1次外推原信号的高频部分，将细节信号叠加到原信号上得到第1次高频补偿的结果，然后逆向递推第2个细节信号，即第2次外推原信号的高频部分，将细节信号叠加到原信号上得到第2次高频补偿的结果，如此逐次逆向递推进行，不断补偿地震信号高频部分，拓展地震信号频宽，从而实现地震资料的高分辨处理。

附图说明

图1为本发明提出的imrsvd分解示意图。

图2为mrsvd分解过程中得到的近似信号的振幅谱图，其中图2中线条由上到下分别代表原始、第10～50次分解。

图3为二维理论模型正演地震剖面图。

图4为imrsvd高分辨处理后的理论模型地震剖面图。

图5为二维实际地震剖面图。

图6为imrsvd高分辨处理后的二维实际地震剖面图。

图7为imrsvd高分辨处理前后第134道数据的振幅谱对比图。

具体实施方式

下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。以下本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

以mrsvd技术为基础，本发明提供了imrsvd算法，用于使得地震信号恢复高分辨率处，此方法的核心在于为了恢复大地滤波引起的地震信号缺失的高频，通过对mrsvd正向分解得到的所有细节奇异值进行拟合外推，得到第1个新的细节信号，即第1次外推原信号的高频部分d′1，将细节信号d′1叠加到原信号x上得到第1次高频补偿的结果x′1，然后逆向递推第2个细节信号d′2，即第2次外推原信号的高频部分d′2，将细节信号d′2叠加到原信号x′1上得到第2次高频补偿的结果x′2，如此逐次逆向递推，不断补偿地震信号高频部分，拓展地震信号频宽，从而实现地震资料的高分辨处理。

因此，本发明中，所述高分辨率处理办法包括以下过程：

基于逆二分递推奇异值分解的地震资料高分辨率处理方法，包括以下步骤：

步骤1：获取单道地震信号x；

步骤2：采用mrsvd算法将地震信号分解，然后利用得到的细节奇异值逐层逆向递推获得新的细节信号和近似信号；

具体的，mrsvd正向分解级数通过以下公式得到：

ej＝∑|aj-1-aj|²/∑|aj-1|²,(j＝1,…,m)

j表示mrsvd正向分解第j层，当ej≤10^-6时，循环分解结束,m是mrsvd正向分解总层数；aj-1和aj分别为第j-1层和第j层分解得到的近似信号。

利用mrsvd分解得到的细节奇异值σd1,σd2,…,σdm，通过拟合函数拟合这些细节奇异值，从而逆向递推出新的细节奇异值σ′di(i＝1,2,…)，再通过细节奇异值得到对应的细节信号d′i，拟合函数为：

其中，j表示mrsvd的分解次数；an代表多项式的系数；k是一个正数，常常小于3；n是多项式的阶数，在最小平方意义下使f(j)逼近已知的细节奇异值，求得k和多项式的系数。

利用新的细节奇异值构建细节矩阵，为从而得到对应的细节信号。

步骤3：将新的细节信号逐次叠加到原始信号上，补偿地震信号的高频部分，从而得到高分辨率的地震信号，采用的公式为：

式中，x表示原信号，a′i表示第i次高频补偿的结果，g表示总的逆向递推次数，d′i为细节信号。

通过修正方差模对总的逆向递推次数进行控制，修正方差模为：

其中，a′i(t)表示第i次高频补偿的结果，t为时间，n为信号的长度，a为常数。对每次高频补偿的信号a'1,a'2,...,a'(g-1),a'g计算其修正方差模为v1,v2,...,v(g-1),vg，如果v(g-6)≈v(g-3)≈vg，即修正方差模收敛并达到最大值，这时总的逆向递推次数g得以确定，且a'g为最终得到的高分辨率地震信号。

本发明中，图3为理论模型数据的原始图件，图4为经过本发明所述方法处理后的理论模型高分辨结果，对比图4、图3可以看到，在图3中的从上往下第二层薄层无法区分、第三层区分度不好，并且楔状模型的可以分辨到第23道，经过imrsvd处理后，如图4所示，从上往下第二层能够有一定程度的区分，第三层已经完全可以区分，楔状模型也由原来只能分辨第23道提高到了可以分辨到第18道。

图5、图6分别为imrsvd处理前后的实际地震数据，对比图5、图6可以看到，经过imrsvd处理后，地震分辨率明显增大，地震同相轴的连续性增强，特别是在1.0秒左右的主要目的层效果尤为明显，我们抽取了imrsvd处理前后的第134道数据进行了振幅谱分析，如图7所示，可以看到经过imrsvd方法处理后，低频部分能够得到有效保持，高频部分得到有效提升，对地震分辨率提升有着重要的意义。

按照上述实施例，便可很好地实现本发明。值得说明的是，基于上述结构设计的前提下，为解决同样的技术问题，即使在本发明上做出的一些无实质性的改动或润色，所采用的技术方案的实质仍然与本发明一样，故其也应当在本发明的保护范围内。