CN103926616B - 一种基于叠前crp道集的多尺度各向异性扩散滤波方法 - Google Patents

Abstract

本发明涉及一种基于叠前CRP道集的多尺度各向异性扩散滤波方法，其包括步骤：1）对输入的叠前CRP道集进行正则化处理；2）采用二维Mallat算法进行多尺度分解，并初始化分解后的每个子剖面；3）经扩散系数法求得扩散门限值，进入步骤6）；进行基于扩散张量法的各向异性扩散求得扩散张量参数后，代入非线性各向异性方程进行各向异性扩散滤波，并进入步骤4）；4）对每次迭代后的子剖面计算SNR、MSE和PSNR，优选出最佳的地震子剖面；5）对各向异性扩散滤波处理后的各个子剖面进行步骤4）处理，进入步骤6）；6）采用二维Mallat重构算法对迭代后的各个子剖面信息进行无损重构叠前CRP道集，重构之后输出最佳的叠前CRP道集。本发明可以广泛用于各种石油勘探地震资料的处理过程中。

Description

一种基于叠前CRP道集的多尺度各向异性扩散滤波方法

技术领域

本发明涉及一种石油勘探过程中的地震数据处理方法，特别是关于一种基于叠前CRP(共反射点)道集的多尺度各向异性扩散滤波方法。

背景技术

在实际生产中，随着地震勘探精度的提高，人们对地震资料品质的要求越来越高。地震资料的“三高”（高分辨率、高信噪比、高保真度）处理技术一直是地球物理勘探领域研究的热点和难点，其中，信噪比制约着分辨率的提高，保真度直接影响资料解释成果的可靠性。随着勘探目标从常规油气藏逐渐向岩性、隐蔽型油气藏转变，叠前反演成为重要的手段，这种手段的主要影响因素是信噪比，因此提高隐蔽油气藏勘探的成功率，关键之一是提高叠前CRP道集资料的信噪比。有效的去噪滤波方法应在压制噪声的同时，更好地保留有效同相轴信息。

各向异性扩散滤波是一种噪音压制方法，可以等价于一个物理热传导的求解过程。偏微分方程在图像处理中的应用由经典的物理问题推导而来，偏微分方程模型首先应用于图像去噪。假设一幅图像是一块热量不同的二维平板，板上不同区域的温度不同。图像中的细节特征、噪声等频域内值较大的像素点可以看作是高温点，而同质、平滑区域的低频像素点可看作是低温点。图像去噪过程可看作是使整块板温度趋于平衡的过程：高温点向低温点输送能量，高温点降温而低温点接受能量升温。当能量传输的过程结束时，整块板中的温度实现平衡。地震剖面中噪声压制的过程就如同温度平衡过程，通过扩散滤波处理提高资料的信噪比。

Perona等人提出了P-M模型（非线性各向异性扩散方程），其保留了图像的边界信息，但它的适定性和稳定性存在问题，并没有考虑局部图像的结构方位。Weickert引入扩散张量代替扩散系数来提取图像结构的方向性信息，改进了P-M模型在边缘细节上损失有效信息的缺陷，后来又提出具有旋转不变性的构造扩散张量的方法，克服了时间步长的限制，但是构造模版较为复杂。

发明内容

针对上述问题，本发明的目的是提供一种能够有效地压制叠前CRP道集中随机噪声和提高资料信噪比品质的基于叠前CRP道集的多尺度各向异性扩散滤波方法。

为实现上述目的，本发明采取以下技术方案：一种基于叠前CRP道集的多尺度各向异性扩散滤波方法，其包括以下步骤：1）对输入的叠前CRP道集进行正则化处理，利用不同方差σ值的高斯核函数G_σ与地震剖面褶积；2）采用二维Mallat算法对正则化处理后的地震剖面做多尺度分解，把对称小波sym6作为基小波，对原始的图像进行N层小波分解；并对分解后的每个子剖面进行初始化；3）将初始化后多尺度分解的各个子剖面经扩散系数法求得扩散门限值，进入步骤5）；同时，将初始化后多尺度分解的各个子剖面进行基于扩散张量法的各向异性扩散求得扩散张量参数后，代入非线性各向异性方程进行各向异性扩散滤波，并进入步骤4）；4）对每一次迭代得到的子剖面，计算信噪比、均方差和峰值信噪比，最终确定最佳迭代次数，根据该迭代次数得到最佳的地震子剖面；5）采用二维Mallat重构算法对迭代后的各个子剖面信息进行无损重构叠前CRP道集，重构之后输出最佳的叠前CRP道集；其中最佳的叠前CRP道集为：

$C_j=H_r^{*}{H}_c^{*}{C}_{j+1}+H_r^{*}{H}_c^{*}{D}_{j+1}^{*}+H_r^{*}{H}_c^{*}{D}_{j+1}^2+H_r^{*}{H}_c^{*}{D}_{j+1}^3,$

式中，C_j+1为低频分量，为水平高频分量，为垂直高频分量，为对角高频分量，和分别代表用尺度滤波器系数对阵列的行和列作用的算子的共轭。

所述步骤2）中，所述多尺度分解包括以下步骤：（1）设为待分析的图像信号，其二维逼近图像为： $A_{j}f=A_{j+1}f+D_{j+1}^{1}f+D_{j+1}^{2}f+D_{j+1}^{3}f,---\left(1\right)$ 其中，

其中，A_j为分解后的第j层的结果，C_j+1(m,n)为低频分量，为水平高频分量，为垂直高频分量，为对角高频分量，为小波基函数，m为按列方向滤波器的平移量n为按行方向滤波器的平移量；利用尺度函数和小波函数的正交性，由式(1)和(2)得：

$C_{j+1}(m,n)=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)h(l-2n)C_j(k,l)$

$\left\{\begin{array}{c}{D}_{j+1}^1=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)g(l-2n)C_j(k,l)\\ D_{j+1}^2=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)h(l-2n)C_j(k,l)\\ D_{j+1}^3=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)g(l-2n)C_j(k,l)\end{array}\right.$

对低频分量C_j+1、水平高频分量垂直高频分量和对角高频分量进行二维Mallat分解：

$\left\{\begin{array}{c}{C}_{j+1}=H_r{H}_c{C}_j\\ D_{j+1}^1=H_r{G}_c{C}_j\\ D_{j+1}^2=G_r{H}_c{C}_j\\ D_{j+1}^3=G_r{G}_c{G}_j\end{array}\right.,j=\mathrm{0,1},...,J$

式中，H_r和H_c分别表示用尺度滤波器系数对阵列的行和列作用的算子，G_r和G_c分别表示用小波滤波器系数对行和列作用的算子；（2）分解后的地震剖面图像的有效信息主要集中在低频分量C_j+1中，直接对低频分量C_j+1进行各向异性扩散滤波处理；对于地震剖面图像的高频分量，对其进行阈值收缩处理；阈值选取公式为：

$\mathrm{\η}=\mathrm{\δ}\·2\sqrt{\log N}/\sqrt{\log (l+1)}\mathrm{\δ}=\mathrm{median}\left(u\right)/0.6745$

其中，η为设定的小波阈值，δ为高频分量的标准偏差估计，u是图像小波分量，median表示取中间数，l为分解层数，N为分量中像素点个数。

所述步骤3）中，所述扩散张量法包括以下步骤：（1）设定扩散滤波方程为：

$\frac{\∂u(x,y,t)}{\∂t}=\▿\·(D\left(\right|{\▿u}_{\mathrm{\σ}}\left|\right)\▿u)$

其中，为扩散张量，σ为尺度函数，为对进行正则化处理，其中

${\▿u}_{\mathrm{\σ}}=\▿(K_{\mathrm{\σ}}*u(x,y,t)),K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\σ}}^2})$

式中，K_σ表示尺度为σ的高斯函数；（2）对初始化后多尺度分解的各个子剖面沿空间和时间方向求导，确定构造张量，并求解构造张量的特征值及单位特征向量，以确定地震剖面局部图像结构：①构造张量T为：假设构造张量

$T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right),$

T为2×2的对称正定矩阵；②采用带加速因子的Gauss-Seidel方法对构造张量T进行求解，得到构造张量T矩阵的特征值为：

$u_1=\frac{1}{2}(T_{11}+T_{12}+\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

$u_2=\frac{1}{2}(T_{11}+T_{12}-\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

单位特征向量v₁，v₂平行于梯度方向，并且单位特征向量v₁满足

$\left(\begin{array}{c}\cos \mathrm{\α}\\ \sin \mathrm{\α}\end{array}\right)=v_1\left|\right|\left(\begin{array}{c}2T_{12}\\ T_{22}-T_{11}+\sqrt{{(T_{11}-T_{22})}^2+{{4T}_{12}}^2}\end{array}\right);$

（3）根据步骤（2）中求得的构造张量T矩阵的特征向量构建扩散张量D，扩散张量D也为对称正定矩阵且其特征向量与构造张量T的特征向量相同，扩散张量D的特征值为λ₁和λ₂，求解扩散张量D矩阵中的参数：扩散张量D为：

$D=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right);$

（4）求出扩散张量参数a、b、c、d之后，实现非线性各向异性方程的求解，扩散张量方程如下所示：

$\frac{\∂u}{\∂t}=\mathrm{div}(D\·{\▿u}^k)---\left(3\right)$

对于扩散张量方程（3）的左边求解，用差分代替微分得到扩散张量方程的离散形式：

$\frac{\∂u}{\∂t}=\frac{u^{k+1}-u^k}{\mathrm{\Δt}}$

式中u^k和u^k+1分别为剖面在k次和k+1次迭代后的结果，Δt为步长；对于扩散张量方程（3）右边求解有：

$\begin{array}{c}\mathrm{div}(D\·{\▿u}^k)=\left(\begin{array}{cc}{\∂}_x& {\∂}_y\end{array}\right)\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)\left(\begin{array}{c}{\∂}_x{u}^k\\ {\∂}_y{u}^k\end{array}\right)\\ ={\∂}_x\left({a\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left({b\∂}_x{u}^k\right)+{\∂}_y\left({c\∂}_y{u}^k\right)\end{array}$

将扩散张量方程（1）左右两边联立起来为：

$u^{k+1}=u^k+\mathrm{\Δt}({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left(b{\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right))---\left(4\right)$

其中，Δt∈(0,0.3]，通过公式（4）实现对分解后的各层剖面进行扩散滤波处理，压制叠前CRP道集中的随机噪声增强同相轴的连续性；（5）根据构造张量T和扩散张量D构建均方根突变算子Θ：

$\mathrm{\Θ}=\sqrt{\frac{\mathrm{\Σ}\left(\mathrm{diag}\left(T^{0}T\right)\right)}{\mathrm{\Σ}\left(\mathrm{diag}\left(T^0\right)\right)\·\mathrm{\Σ}\left(\mathrm{diag}\left(T\right)\right)}};$

其中，diag(·)为矩阵T⁰T的主对角线，T⁰初始构造张量，T为当前迭代构造张量，Θ∈[0,1]，在突变附近，Θ趋近于0；远离突变时，Θ趋近于1；将均方根突变算子Θ加到扩散张量方程(3)中，得到非线性各向异性方程为：

$u^{k+1}=u^k+\mathrm{\Δt\Θ}\·({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right)).$

所述步骤（3）中，所述扩散张量D矩阵中的四个参数a、b、c和d的求解如下：①求解扩散张量D的特征值λ₁和λ₂，根据实际的地震剖面要求采用突变增强系数方法和相干增强系数方法两种方法估计特征值λ₁和λ₂：突变增强系数方法：取特征值λ₂=1；特征向量v₁对应的特征值λ₁为构造张量T特征值u₁的单调减函数，且

λ₂＝1

相干增强系数方法：令特征值λ₁和λ₂取值如下：

λ₁＝c₁

其中，为类波阻抗因子；②特征值λ₁和λ₂求解扩散张量D的参数a、b、c和d：

$D={(v_1,v_2)}^T\·\left(\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right)\·(v_1,v_2)=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)$

$\left\{\begin{array}{c}a={\mathrm{\λ}}_1{\cos }^{2}a+{\mathrm{\λ}}_2{\sin }^{2}a\\ b=({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)\cos a\sin a\\ c={\mathrm{\λ}}_1{\sin }^{2}a+{\mathrm{\λ}}_2{\cos }^{2}a\end{array}\right..$

所述步骤3）中，所述扩散系数法包括以下步骤：（1）设定扩散系数：将Biweight范数作为扩散系数

式中，k为扩散门限值；（2）根据扩散系数采用邻域内梯度绝对偏差估计方法求解扩散门限值，扩散门限k(l)与分解层数l相关：

$k\left(l\right)=\frac{1}{0.6745}\·\mathrm{median}\left(\right||\▿u\left(l\right)-\mathrm{median}\left(\right||\▿u\left(l\right)|\left|\right)|\left|\right)\sqrt{2\log (l+1)/l}.$

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于采用软阈值小波处理高频分量，因此可以减少对剖面有效信息的损伤，提高了扩散过程中剖面上突变信息的准确性，即能够在强噪声环境中区分由噪声引起的灰度值变化和突变引起的灰度值变化。2、本发明根据待处理剖面的特点提出了具有特定响应特性的扩散方程，并寻求降噪和细节保持间的平衡，因此可以对分解后的各层剖面进行扩散滤波处理，压制叠前道集中的随机噪声，同时增强同相轴的连续性，占用的内存少，计算效率高。3、本发明通过提高扩散参数及收敛条件对不同统计特性剖面的适应性，即自适应地估计扩散门限、扩散迭代次数等参数，因此有效地提高了扩散结果的鲁棒性。4、应用本发明多尺度各向异性扩散滤波方法对地震剖面进行处理，不但使噪声得到有效压制，资料的信噪比得到改善；而且有效地增强了同相轴的连续性。5、本发明由于引入均方根突变算子，使突变信息得到保留，因此真正达到有效保留原始有效信息和压制噪声的效果。6、本发明在实际地震资料处理中，既能增强同相轴连续性，又可以保持叠前CRP道集中的突变信息，增强对突变信息的识别能力；采用本发明处理后的剖面地震波特征自然，信噪比得到提高，横向连续性得到改善，构造信息更加丰富。本发明可以广泛用于各种石油勘探地震资料的处理过程中。

附图说明

图1是本发明各向异性扩散滤波流程示意图

图2是采用已有技术Focus去噪模块滤除随机噪音的结果示意图

图3是采用本发明滤除随机噪音的结果示意图

图4是采用本发明滤除噪声后的剖面同相轴信息同其它方法对比结果示意图

图5是对图4中各CRP道集煤层顶界面振幅追踪的曲线示意图

图6是去噪前后的CRP道集不同角度的叠加剖面对比示意图

图7是某地区A井流程采用现有技术优化前、后叠前资料的VP/VS反演剖面示意图

图8是某地区A井流程采用本发明优化前、后叠前资料的VP/VS反演剖面示意图

图9是基于常规流程和本发明优化前、后资料得到的工区VP/VS反演时间切片对比示意图

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

如图1所示，本发明的基于叠前CRP道集的多尺度各向异性扩散滤波方法包括以下步骤：

1）对输入的叠前CRP道集进行正则化处理，即利用不同方差σ值的高斯核函数G_σ（σ控制平滑程度，σ值越大图像越平滑，实际中根据需求选取）与地震剖面褶积，这样可以消除部分奇异值，得到平滑剖面，使得后续处理更加稳定。

2）采用二维Mallat（马拉）算法对正则化处理后的地震剖面做多尺度分解，把对称小波sym6作为基小波，对原始的图像进行N层小波分解；并对分解后的每个子剖面进行初始化，如给定迭代次数及迭代步长。

其中，多尺度分解包括以下步骤：

（1）设为待分析的图像信号，其二维逼近图像为：

$A_{j}f=A_{j+1}f+D_{j+1}^{1}f+D_{j+1}^{2}f+D_{j+1}^{3}f,---\left(1\right)$

其中，

其中，A_j为分解后的第j层的结果，C_j+1(m,n)为低频分量，为水平高频分量，为垂直高频分量，为对角高频分量，为小波基函数，m为按列方向滤波器的平移量n为按行方向滤波器的平移量。

利用尺度函数和小波函数的正交性，由式(1)和(2)可得：

$C_{j+1}(m,n)=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)h(l-2n)C_j(k,l)---\left(3\right)$

$\left\{\begin{array}{c}{D}_{j+1}^1=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)g(l-2n)C_j(k,l)\\ D_{j+1}^2=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)h(l-2n)C_j(k,l)\\ D_{j+1}^3=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)g(l-2n)C_j(k,l)\end{array}\right.---\left(4\right)$

对低频分量C_j+1、水平高频分量垂直高频分量和对角高频分量进行二维Mallat分解：

$\left\{\begin{array}{c}{C}_{j+1}=H_r{H}_c{C}_j\\ D_{j+1}^1=H_r{G}_c{C}_j\\ D_{j+1}^2=G_r{H}_c{C}_j\\ D_{j+1}^3=G_r{G}_c{G}_j\end{array}\right.,j=\mathrm{0,1},...,J---\left(5\right)$

式中，H_r和H_c分别表示用尺度滤波器系数对阵列的行和列作用的算子，G_r和G_c分别表示用小波滤波器系数对行和列作用的算子。

（2）分解后的地震剖面图像的有效信息主要集中在低频分量C_j+1中，可直接对低频分量C_j+1进行各向异性扩散滤波处理（即扩散张量法）；对于地震剖面图像的高频分量，主要是一些细微的突变点或随机噪声，因此要对它进行阈值收缩处理（即扩散系数法）。其中阈值选取公式为：

$\mathrm{\η}=\mathrm{\δ}\·2\sqrt{\log N}/\sqrt{\log (l+1)},\mathrm{\δ}=\mathrm{median}\left(u\right)/0.6745---\left(6\right)$

其中，η为设定的小波阈值，δ为高频分量的标准偏差估计，u是图像小波分量，median表示取中间数，l为分解层数，N为分量中像素点个数。

选取的阈值收缩方法为软阈值小波处理，这样可以减少对剖面有效信息的损伤，则软阈值大小为：

$\hat{u}=\left\{\begin{array}{cc}\mathrm{sgn}\left(u_j\right)\left(\right|u_j-\mathrm{\η}\left|\right)& \left|u_j\right|\≥\mathrm{\η}\\ 0& \left|u_j\right|\≤\mathrm{\η}\end{array}\right..---\left(7\right)$

3）将初始化后多尺度分解的各个子剖面经扩散系数法求得扩散门限值，进入步骤5）；同时，将初始化后多尺度分解的各个子剖面进行基于扩散张量法的各向异性扩散求得扩散张量参数后，代入非线性各向异性方程进行各向异性扩散滤波，并进入步骤4）；其中非线性各向异性方程为：

$\frac{\∂u}{\∂t}=\mathrm{div}[c\left(\right||\▿(G_{\mathrm{\σ}}*u_k)|\left|\right)\·\▿u_k],---\left(8\right)$

式中，u_k为子剖面，为图像梯度，G_σ表示方差为σ的高斯核函数，对剖面具有平滑作用，div为散度。

4）对每一次迭代得到的子剖面，计算信噪比SNR、均方差MSE和峰值信噪比PSNR，最终确定最佳迭代次数，根据该迭代次数得到最佳的地震子剖面。

5）采用二维Mallat重构算法对迭代后的各个子剖面信息进行无损重构叠前CRP道集，重构之后可以输出最佳的叠前CRP道集；其中最佳的叠前CRP道集为：

$C_j=H_r^{*}{H}_c^{*}{C}_{j+1}+H_r^{*}{H}_c^{*}{D}_{j+1}^{*}+H_r^{*}{H}_c^{*}{D}_{j+1}^2+H_r^{*}{H}_c^{*}{D}_{j+1}^3---\left(9\right)$

式中，C_j+1为低频分量，为水平高频分量，为垂直高频分量，为对角高频分量，和分别代表用尺度滤波器系数对阵列的行和列作用的算子的共轭。

上述步骤3）中，基于扩散张量的滤波方法是根据局部剖面信息构建构造张量，再由构造张量得到扩散张量，扩散张量法包括以下步骤：

（1）设定扩散滤波方程为：

$\frac{\∂u(x,y,t)}{\∂t}=\▿\·(D\left(\right|{\▿u}_{\mathrm{\σ}}\left|\right)\▿u)---\left(10\right)$

其中，为扩散张量，σ为尺度函数，可以根据噪声的均方差自适应调节，为对进行正则化处理，相当于作用高斯滤波，转化为一个适定问题，如下：

${\▿u}_{\mathrm{\σ}}=\▿(K_{\mathrm{\σ}}*u(x,y,t)),K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\σ}}^2})---\left(11\right)$

式中，K_σ表示尺度为σ的高斯函数。

（2）对初始化后多尺度分解的各个子剖面沿空间和时间方向求导，进而确定构造张量，并求解构造张量的特征值及单位特征向量，以确定地震剖面局部图像结构：

①构造张量T为：

$T(\▿u)=\▿u_{\mathrm{\σ}}\⊗{\▿u}_{\mathrm{\σ}}---\left(12\right)$

假设构造张量 $T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right),$ T为2×2的对称正定矩阵；

②由于构造张量T为对称正定矩阵，因此采用带加速因子的Gauss-Seidel（高斯-赛德尔）方法对构造张量T进行求解，采用Gauss-Seidel方法的优势在于：只要构造张量T是对称正定的，必收敛；在求解过程中不断地用来及时替代之前的数，这样占用的内存少，计算效率高；其求解过程如下：

$x_i^{k+1}=(b_i-\underset{j=1}{\overset{i-1}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^{k+1}-\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^k)/a_{\mathrm{ii}},i=\mathrm{1,2},...n---\left(13\right)$

为了提高计算效率，在迭代过程中引入最佳加速因子α：

$x_i^{k+1}=x_i^k+\mathrm{\α}(x_i^{k+1}-x_i^k),\mathrm{\α}=\frac{2}{1+\sqrt{1-{\mathrm{\ρ}}^2\left(T\right)}}---\left(14\right)$

其中，ρ为迭代矩阵的谱半径。

通过求解可知，构造张量T矩阵的特征值为：

$\begin{array}{c}{u}_1=\frac{1}{2}(T_{11}+T_{12}+\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})\\ u_2=\frac{1}{2}(T_{11}+T_{12}-\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})\end{array}---\left(15\right)$

单位特征向量v₁，v₂平行于梯度方向，并且单位特征向量v₁满足

$\left(\begin{array}{c}\cos \mathrm{\α}\\ \sin \mathrm{\α}\end{array}\right)=v_1\left|\right|\left(\begin{array}{c}2T_{12}\\ T_{22}-T_{11}+\sqrt{{(T_{11}-T_{22})}^2+{{4T}_{12}}^2}\end{array}\right)---\left(16\right)$

u₁、u₂为构造张量T矩阵的两个特征值，且有特征值u₁大于等于特征值u₂，因此特征值u₁的方向就指向数据波动较大的方向，即垂直叠前CRP道集同相轴的方向，指示道集中的突变信息；特征值u₂指向数据波动较小的方向，也就是平行于同相轴的方向。v₁、v₂为构造张量T矩阵的两个特征向量，两者之间相互正交，特征向量v₁与特征值u₁相对应，它的方向也与同相轴方向垂直，特征向量v₂与特征值u₂相对应，它平行于同相轴方向。

（3）根据步骤（2）中求得的构造张量T矩阵的特征向量构建扩散张量D，扩散张量D也为对称正定矩阵且其特征向量与构造张量T的特征向量相同，扩散张量D的特征值为λ₁和λ₂，并求解扩散张量D矩阵中的参数。其中，扩散张量D为：

$D=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right).$

扩散张量D矩阵中的四个参数a、b、c和d的求解如下：

1○求解扩散张量D的特征值λ₁和λ₂，根据实际的地震剖面要求采用突变增强系数方法和相干增强系数方法两种方法估计特征值λ₁和λ₂，突变增强方法可以保留更多的细节突变信息，相干增强方法可以使叠前CRP道集同相轴的连续性增强。

突变增强系数方法：为了达到增强突变细节信息的目的，取特征值λ₂=1；特征向量v₁对应的特征值λ₁为构造张量T特征值u₁的单调减函数，且

λ₂＝1

相干增强系数方法：为了达到增强同相轴连续性的目的，令特征值λ₁和λ₂取值如下所示：

λ₁＝c₁

其中，为类波阻抗因子，表示二维数据在某个方向上具有旋转不变的性质，避开了寻找最优旋转不变矩阵的过程，使扩散过程大大简化；

②特征值λ₁和λ₂求解扩散张量D的参数a、b、c和d：

$D={(v_1,v_2)}^T\·\left(\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right)\·(v_1,v_2)=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)---\left(19\right)$

$\left\{\begin{array}{c}a={\mathrm{\λ}}_1{\cos }^{2}a+{\mathrm{\λ}}_2{\sin }^{2}a\\ b=({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)\cos a\sin a\\ c={\mathrm{\λ}}_1{\sin }^{2}a+{\mathrm{\λ}}_2{\cos }^{2}a\end{array}\right.---\left(22\right)$

特征值λ₁、λ₂是用来控制扩散程度的，λ₁为一较小的正实数，表示扩散程度很小，λ₂的数值大小与迭代矩阵特征值u₁、u₂关系密切，当迭代矩阵特征值u₁＝u₂时，表示扩散程度同样很小，但是当u₁≠u₂，就意味着叠前CRP道集同相轴存在突变特征，扩散需沿着突变方向进行，特征值λ₂越大，意味着扩散程度的增强；特征向量v₁、v₂用来控制扩散方向，通过增强沿叠前CRP道集同相轴方向上的扩散程度、减小垂直同相轴方向上的扩散程度，达到对增强叠前CRP道集同相轴的目的。

（4）求出扩散张量参数之后，实现非线性各向异性方程的求解，扩散张量方程如下所示：

$\frac{\∂u}{\∂t}=\mathrm{div}(D\·\▿u^k)---\left(20\right)$

对于扩散张量方程(20)的左边求解，可以用差分代替微分，得到扩散张量方程的离散形式：

$\frac{\∂u}{\∂t}=\frac{u^{k+1}-u^k}{\mathrm{\Δt}}---\left(21\right)$

式中u^k和u^k+1分别为剖面在k次和k+1次迭代后的结果，Δt为步长。

对于扩散张量方程(20)右边求解有：

$\begin{array}{c}\mathrm{div}(D\·{\▿u}^k)=\left(\begin{array}{cc}{\∂}_x& {\∂}_y\end{array}\right)\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)\left(\begin{array}{c}{\∂}_x{u}^k\\ {\∂}_y{u}^k\end{array}\right)\\ ={\∂}_x\left({a\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left({b\∂}_x{u}^k\right)+{\∂}_y\left({c\∂}_y{u}^k\right)\end{array}---\left(22\right)$

将扩散张量方程(20)左右两边联立起来为：

$u^{k+1}=u^k+\mathrm{\Δt}({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left(b{\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right))---\left(23\right)$

其中，Δt∈(0,0.3]，通过公式（23）就可实现对分解后的各层剖面进行扩散滤波处理，压制叠前CRP道集中的随机噪声增强同相轴的连续性。

（5）根据构造张量T和扩散张量D构建均方根突变算子Θ，均方根突变算子Θ既能增强同相轴连续性又可以保持叠前CRP道集中的突变信息，增强对突变信息的识别能力：

$\mathrm{\Θ}=\sqrt{\frac{\mathrm{\Σ}\left(\mathrm{diag}\left(T^{0}T\right)\right)}{\mathrm{\Σ}\left(\mathrm{diag}\left(T^0\right)\right)\·\mathrm{\Σ}\left(\mathrm{diag}\left(T\right)\right)}}---\left(24\right)$

其中，diag(·)为矩阵T⁰T的主对角线，T⁰初始构造张量，T为当前迭代构造张量，Θ∈[0,1]，在突变附近，Θ趋近于0；远离突变时，Θ趋近于1。

将均方根突变算子Θ加到扩散张量方程(20)中，得到非线性各向异性方程为：

$u^{k+1}=u^k+\mathrm{\Δt\Θ}\·({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right)).---\left(25\right)$

上述步骤3）中，扩散系数法包括以下步骤：

（1）设定扩散系数：为使扩散过程对参数设置更具鲁棒性,将Biweight（比瓦特）范数作为扩散系数

式中，k为扩散门限值，当剖面邻域梯度小于扩散门限k时，当前像素与其邻域像素属于同一光滑区域，则邻域像素对当前像素平滑的贡献较大；当邻域梯度大于扩散门限k并小于门限时，当前像素与其邻域像素属于同一光滑区域的概率较小，则邻域像素对当前像素平滑的贡献减小；当邻域梯度大于扩散门限时，当前像素与其邻域像素不属于同一光滑区域，则当前像素的平滑不受其邻域像素的影响。

（2）根据扩散系数采用邻域内梯度绝对偏差估计（MAD）方法求解扩散门限值：为了控制在不同分辨率级图像上进行尺度不同的滤波，采用MAD方法，扩散门限k(l)与分解层数l相关:

$k\left(l\right)=\frac{1}{0.6745}\·\mathrm{median}\left(\right||\▿u\left(l\right)-\mathrm{median}\left(\right||\▿u\left(l\right)|\left|\right)|\left|\right)\sqrt{2\log (l+1)/l}---\left(27\right)$

在上述步骤4）中，迭代次数定量分析包括以下三种方法：

（1）均方差（MSE）方法：

$\mathrm{MSE}=\frac{1}{M*N}\underset{(i,j)=1}{\overset{M,N}{\mathrm{\Σ}}}{(u*(i,j)-u(i,j))}^2---\left(28\right)$

式中，MSE为均方差值,u^*和u分别表示滤波后的图像和原始图像。MSE值越小，降噪效果越好。

（2）信噪比（SNR）方法：

$\mathrm{SNR}\left(\mathrm{dB}\right)=10{\log }_{10}\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}u{(i,j)}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(u^{*}(i,j)-u(i,j))}^2}---\left(29\right)$

u^*和u分别表示滤波后的图像和原始图像，处理后SNR越大说明处理效果越好。

（3）峰值信噪比（PSNR）方法：

$\mathrm{PSNR}=10{\log }_{10}\left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(30\right)$

其中，MSE为均方差值，峰值信噪比PSNR越大表示滤波的效果越明显。

根据上述三种定量分析标准，优选出三者达到统计最优的迭代次数，使达到平滑同相轴的同时保持边缘细节信息，保持较好的保真度。

如图2所示，其中图2（a）为原始道集，图2（b）为采用Focus去噪模块滤除随机噪音后的叠前CRP道集，图2（c）为滤出的噪声结果，从图中可以得出：Focus去噪模块对随机噪音的滤除效果不甚理想，从滤除的噪音剖面上，该方法较大地损害了目的层的有效反射信号以及浅部的同相轴信息。

如图3所示，其中图3（a）为原始道集，图3（b）为采用本发明滤除随机噪音后的叠前CRP道集，图3（c）为滤出的噪声结果；从图中可以得出：本发明不仅较好地滤除了随机噪音、保护有效同相轴信息，而且目的层上下的同相轴特征清晰，横向连续性得到较好改善，同相轴有强有弱，波组特征更加明显，远近偏移距的能量一致性加强，振幅的保真度较好。从滤除的噪音剖面上可以看到，本发明对有效信息的损伤较小，克服了现有技术Focus模块对同相轴较大的损伤的缺陷。

如图4所示，其中，a为测井资料，b为根据井上反射系数利用射线模型生成的正演道集，c是采用本发明压噪后道集，d是井控振幅恢复后道集，e是原始道集；从图中可以看出：原始道集e的信噪比较低，对同相轴的追踪困难，井控振幅恢复后道集d较好地恢复了能量，但同相轴的连续性、信噪比明显不如改进的扩散滤波压噪法。正演道集b和采用本发明压噪后道集不仅在标志层附近具有良好的吻合度，且其余深度的同相轴信息也有较好的吻合度，对真实地下反射层进行有效的识别；然而在井控振幅恢复后道集d和原始道集e上几乎无法识别这些有效的同相轴信息。因此本发明改善了剖面信噪比，且其分辨率是真实可靠的。

如图5所示，其中，A线为原始道集振幅曲线，B线为射线模型曲线，C线为井控振幅恢复后道集振幅曲线，D线为本发明压噪曲线；通过对比各条曲线可知，射线模型曲线B的特征是随着偏移距的增加，振幅逐渐增加，趋势较为平缓；原始道集振幅曲线A的趋势与模型基本不符，近偏移距处振幅较为稳定，随着偏移距的增加有一个较大幅度的下降，之后有一定程度的回升，远偏移距也有波动；井控振幅恢复曲线C在近偏移距的变化趋势与原始振幅曲线相近，中、远偏移距的变化趋势与模型曲线一致，但局部偏离幅度振荡较大；本发明压噪曲线D在近偏移距与射线模型曲线B吻合度较高，在中偏移距略小于射线模型曲线B，在远偏移距高于射线模型曲线B，但整体趋势走向与射线模型曲线B基本吻合，相似系数较高，且光滑性较好，振荡较少。通过对比处理后的各振幅曲线，可以看出本发明压噪曲线与射线模型曲线B吻合度很高，具有较好的保幅特性，有利于提取叠前CRP道集中的AVO属性和反演。

如图6所示，其中，分别选取A线：0-15°，B线：15°-30°，C线：30°-48°的CRP道集进行叠加，图6(a)为去噪前叠加剖面示意图，图6(b)为采用本发明去噪后叠加剖面示意图；通过比较去噪前后的叠加结果发现，去噪前，远、中、近偏移距的能量、相位在横向上的一致性较差，波形存在较大的差异；采用本发明去噪后，远、中、近偏移能量、相位的连续性得到有效的改善，不同信噪比的有效信号能量得到有效的恢复，且去噪后振幅的保真度较好。

如图7、图8所示，对比图7和图8可知，基于本发明的一体化流程得到的VP/VS反演剖面，其储层的砂泥岩横向连续性、清晰度比采用现有技术去噪反演剖面要高，结合测井资料对比，本发明处理后反演出的同相轴信息是真实有效的，同时砂泥岩分布特征也得到改善。

如图9所示，其中图9(a)采用常规流程得到的反演切片，显示6井在某时间深度可能是泥岩，5井可能是砂岩，而图9(b)采用本发明去噪后反演切片上显示的结果恰好相反，6井（实线圈）目标层是砂岩的可能性更大，5井（虚线圈）是砂岩的可能性较小。真实钻井资料显示，6井目的层有油气显示，5井目的层则无油气显示。因此可以说明采用本发明对叠前资料的分辨率、信噪比和保真度均取得了良好的效果。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (5)

1.一种基于叠前CRP道集的多尺度各向异性扩散滤波方法，其包括以下步骤：

1）对输入的叠前CRP道集进行正则化处理，利用不同方差σ值的高斯核函数G_σ与地震剖面褶积；

2）采用二维Mallat算法对正则化处理后的地震剖面做多尺度分解，把对称小波sym6作为基小波，对原始的图像进行N层小波分解；并对分解后的每个子剖面进行初始化；

3）将初始化后多尺度分解的各个子剖面经扩散系数法求得扩散门限值，进入步骤5）；同时，将初始化后多尺度分解的各个子剖面进行基于扩散张量法的各向异性扩散求得扩散张量参数后，代入非线性各向异性方程进行各向异性扩散滤波，并进入步骤4）；

4）对每一次迭代得到的子剖面，计算信噪比、均方差和峰值信噪比，最终确定最佳迭代次数，根据该迭代次数得到最佳的地震子剖面；

5）采用二维Mallat重构算法对迭代后的各个子剖面信息进行无损重构叠前CRP道集，重构之后输出最佳的叠前CRP道集；其中最佳的叠前CRP道集为：

$C_j=H_r^{*}{H}_c^{*}{C}_{j+1}+H_r^{*}{H}_c^{*}{D}_{j+1}^{*}+H_r^{*}{H}_c^{*}{D}_{j+1}^2+H_r^{*}{H}_c^{*}{D}_{j+1}^3,$

式中，C_j+1为低频分量，为水平高频分量，为垂直高频分量，为对角高频分量，和分别代表用尺度滤波器系数对阵列的行和列作用的算子的共轭。

2.如权利要求1所述的基于叠前CRP道集的多尺度各向异性扩散滤波方法，其特征在于：所述步骤2）中，所述多尺度分解包括以下步骤：

（1）设为待分析的图像信号，其二维逼近图像为：

$A_{j}f=A_{j+1}f+D_{j+1}^{1}f+D_{j+1}^{2}f+D_{j+1}^{3}f,---\left(1\right)$

其中，

其中，A_j为分解后的第j层的结果，C_j+1(m,n)为低频分量，为水平高频分量，为垂直高频分量，为对角高频分量，为小波基函数，m为按列方向滤波器的平移量n为按行方向滤波器的平移量；

利用尺度函数和小波函数的正交性，由式(1)和(2)得：

$C_{j+1}(m,n)=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)h(l-2n)C_j(k,l)$

$\left\{\begin{array}{c}{D}_{j+1}^1=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h(k-2m)g(l-2n)C_j(k,l)\\ D_{j+1}^2=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)h(l-2n)C_j(k,l)\\ D_{j+1}^3=\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{l=-\∞}{\overset{\∞}{\mathrm{\Σ}}}g(k-2m)g(l-2n)C_j(k,l)\end{array}\right.$

对低频分量C_j+1、水平高频分量垂直高频分量和对角高频分量进行二维Mallat分解：

$\left\{\begin{array}{c}{C}_{j+1}=H_r{H}_c{C}_j\\ D_{j+1}^1=H_r{G}_c{C}_j\\ D_{j+1}^2=G_r{H}_c{C}_j\\ D_{j+1}^3=G_r{G}_c{G}_j\end{array}\right.,j=\mathrm{0,1},...,J$

式中，H_r和H_c分别表示用尺度滤波器系数对阵列的行和列作用的算子，G_r和G_c分别表示用小波滤波器系数对行和列作用的算子；

（2）分解后的地震剖面图像的有效信息主要集中在低频分量C_j+1中，直接对低频分量C_j+1进行各向异性扩散滤波处理；对于地震剖面图像的高频分量，对其进行阈值收缩处理；阈值选取公式为：

$\mathrm{\η}=\mathrm{\δ}\·2\sqrt{\log N}/\sqrt{\log (l+1)},\mathrm{\δ}=\mathrm{median}\left(u\right)/0.6745$

其中，η为设定的小波阈值，δ为高频分量的标准偏差估计，u是图像小波分量，median表示取中间数，l为分解层数，N为分量中像素点个数。

3.如权利要求1或2所述的基于叠前CRP道集的多尺度各向异性扩散滤波方法，其特征在于：所述步骤3）中，所述扩散张量法包括以下步骤：

（1）设定扩散滤波方程为：

$\frac{\∂u(x,y,t)}{\∂t}=\▿\·(D\left(\right|{\▿u}_{\mathrm{\σ}}\left|\right)\▿u)$

其中，为扩散张量，σ为尺度函数，为对进行正则化处理，其中

${\▿u}_{\mathrm{\σ}}=\▿(K_{\mathrm{\σ}}*u(x,y,t)),K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\σ}}^2})$

式中，K_σ表示尺度为σ的高斯函数；

（2）对初始化后多尺度分解的各个子剖面沿空间和时间方向求导，确定构造张量，并求解构造张量的特征值及单位特征向量，以确定地震剖面局部图像结构：

①构造张量T为：

$T(\▿u)=\▿u_{\mathrm{\σ}}\⊗{\▿u}_{\mathrm{\σ}}$

假设构造张量 $T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right),$ T为2×2的对称正定矩阵；

②采用带加速因子的Gauss-Seidel方法对构造张量T进行求解，得到构造张量T矩阵的特征值为：

$u_1=\frac{1}{2}(T_{11}+T_{12}+\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

$u_2=\frac{1}{2}(T_{11}+T_{12}-\sqrt{{(T_{11}-T_{22})}^2+{4T}_{12}^2})$

单位特征向量v₁，v₂平行于梯度方向，并且单位特征向量v1满足

$\left(\begin{array}{c}\cos \mathrm{\α}\\ \sin \mathrm{\α}\end{array}\right)=v_1\left|\right|\left(\begin{array}{c}2T_{12}\\ T_{22}-T_{11}+\sqrt{{(T_{11}-T_{22})}^2+{{4T}_{12}}^2}\end{array}\right);$

（3）根据步骤（2）中求得的构造张量T矩阵的特征向量构建扩散张量D，扩散张量D也为对称正定矩阵且其特征向量与构造张量T的特征向量相同，扩散张量D的特征值为λ₁和λ₂，求解扩散张量D矩阵中的参数：扩散张量D为：

$D=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right);$

（4）求出扩散张量参数a、b、c、d之后，实现非线性各向异性方程的求解，扩散张量方程如下所示：

$\frac{\∂u}{\∂t}=\mathrm{div}(D\·{\▿u}^k)---\left(3\right)$

对于扩散张量方程（3）的左边求解，用差分代替微分得到扩散张量方程的离散形式：

$\frac{\∂u}{\∂t}=\frac{u^{k+1}-u^k}{\mathrm{\Δt}}$

式中u^k和u^k+1分别为剖面在k次和k+1次迭代后的结果，Δt为步长；

对于扩散张量方程（3）右边求解有：

$\begin{array}{c}\mathrm{div}(D\·{\▿u}^k)=\left(\begin{array}{cc}{\∂}_x& {\∂}_y\end{array}\right)\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)\left(\begin{array}{c}{\∂}_x{u}^k\\ {\∂}_y{u}^k\end{array}\right)\\ ={\∂}_x\left({a\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left({b\∂}_x{u}^k\right)+{\∂}_y\left({c\∂}_y{u}^k\right)\end{array}$

将扩散张量方程（1）左右两边联立起来为：

$u^{k+1}=u^k+\mathrm{\Δt}({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left(b{\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right))---\left(4\right)$

其中，Δt∈(0，0.3]，通过公式（4）实现对分解后的各层剖面进行扩散滤波处理，压制叠前CRP道集中的随机噪声增强同相轴的连续性；

（5）根据构造张量T和扩散张量D构建均方根突变算子Θ：

$\mathrm{\Θ}=\sqrt{\frac{\mathrm{\Σ}\left(\mathrm{diag}\left(T^{0}T\right)\right)}{\mathrm{\Σ}\left(\mathrm{diag}\left(T^0\right)\right)\·\mathrm{\Σ}\left(\mathrm{diag}\left(T\right)\right)}};$

其中，diag(·)为矩阵T⁰T的主对角线，T⁰初始构造张量，T为当前迭代构造张量，Θ∈[0,1]，在突变附近，Θ趋近于0；远离突变时，Θ趋近于1；

将均方根突变算子Θ加到扩散张量方程(3)中，得到非线性各向异性方程为：

$u^{k+1}=u^k+\mathrm{\Δt\Θ}\·({\∂}_x\left(a{\∂}_x{u}^k\right)+{\∂}_x\left({b\∂}_y{u}^k\right)+{\∂}_y\left(b{\∂}_x{u}^k\right)+{\∂}_y\left(c{\∂}_y{u}^k\right)).$

4.如权利要求3所述的基于叠前CRP道集的多尺度各向异性扩散滤波方法，其特征在于：所述步骤（3）中，所述扩散张量D矩阵中的四个参数a、b、c和d的求解如下：

①求解扩散张量D的特征值λ₁和λ₂，根据实际的地震剖面要求采用突变增强系数方法和相干增强系数方法两种方法估计特征值λ₁和λ₂：

突变增强系数方法：取特征值λ₂=1；特征向量v₁对应的特征值λ₁为构造张量T特征值u₁的单调减函数，且

λ₂＝1

相干增强系数方法：令特征值λ₁和λ₂取值如下：

λ₁＝c₁

其中，为类波阻抗因子；

②特征值λ₁和λ₂求解扩散张量D的参数a、b、c和d：

$D={(v_1,v_2)}^T\·\left(\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\end{array}\right)\·(v_1,v_2)=\left(\begin{array}{cc}a& b\\ b& c\end{array}\right)$

$\left\{\begin{array}{c}a={\mathrm{\λ}}_1{\cos }^{2}a+{\mathrm{\λ}}_2{\sin }^{2}a\\ b=({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)\cos a\sin a\\ c={\mathrm{\λ}}_1{\sin }^{2}a+{\mathrm{\λ}}_2{\cos }^{2}a\end{array}\right..$

5.如权利要求1或2所述的基于叠前CRP道集的多尺度各向异性扩散滤波方法，其特征在于：所述步骤3）中，所述扩散系数法包括以下步骤：

（1）设定扩散系数：将Biweight范数作为扩散系数：

式中，k为扩散门限值；

（2）根据扩散系数采用邻域内梯度绝对偏差估计方法求解扩散门限值，扩散门限k(l)与分解层数l相关：

$k\left(l\right)=\frac{1}{0.6745}\·\mathrm{median}\left(\right||\▿u\left(l\right)-\mathrm{median}\left(\right||\▿u\left(l\right)|\left|\right)|\left|\right)\sqrt{2\log (l+1)/l}.$