CN105182418A - 一种基于双树复小波域的地震信号降噪方法及系统 - Google Patents

Abstract

本发明涉及一种基于双树复小波域的地震信号降噪方法及系统。所述方法包括：步骤a：建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；步骤b：对同一方向实部小波系数或虚部小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；步骤c：利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。本发明考虑了实部(或虚部)系数与对应模之间的相关性，对实部(或虚部)小波系数与对应的模建立双变量模型进行降噪，有效压制随机噪声，地震信号整体变得清晰，提高信号质量。

Description

一种基于双树复小波域的地震信号降噪方法及系统

技术领域

本发明属于地震信号降噪技术领域，尤其涉及一种基于双树复小波域的地震信号降噪方法及系统。

背景技术

随着地震勘探的不断深入以及一些先进技术在地震勘探中的应用，人们对地震资料信噪比的要求也日益提高，这使得对降噪方法的研究显得尤为重要。实际中各种地震噪声的特殊性、复杂性给降噪工作带来了一定的难度，为此人们不遗余力地研究更好的降噪方法。

以傅氏变换为基础的频率域滤波方法是地震资料处理中的经典方法，该方法根据有效波和干扰波在频带上的差异来压制干扰波突出有效波，其不足在于信号的降噪效果依赖于信号频带和噪声频带分离的程度，若混入噪声的频带和信号频带有重叠时，则降噪效果不是十分理想。一些学者利用小波变换这一时频分析工具解决非稳态问题，相继提出了基于小波变换的降噪方法。小波包f-x域前后向预测去噪方法，能较好地抑制叠后地震剖面中的随机噪声。用小波分析和反演方法联合去除随机噪声，是将地震信号分解成若干频段，对高频随机噪声所在的频段利用常规的时空域反演方法衰减噪声，比直接用反演方法取得更好的去噪效果。基于小波变换的最小二乘光滑去噪法先对信号进行小波多尺度分解，再在每个尺度上利用光滑滤波来去除地震资料中的随机噪声。基于贝叶斯估计的小波域双变量法利用了相邻尺度间小波系数的依赖关系，取得较好的降噪效果。

综上所述，现有的地震信号降噪方法至少存在以下技术问题：

(1)、通常小波变换降噪方法将地震信号进行多尺度分解得到各子带系数，然后对包含随机噪声的子带系数进行处理，以达到降噪目的，但并未考虑小波域统计模型；

(2)、软阈值法开创性地分析了小波系数的统计特征，SURE-LET算法的先验模型在尺度内较好地模拟小波系数的概率分布，但忽略了小波系数尺度间的相关性；

(3)、小波变换系数具有零均值和重尾性，小波系数的边缘分布近似服从拉普拉斯分布、广义高斯分布或混合高斯分布等这类非高斯统计分布，但是仅仅把这些边缘分布特征作为小波系数的模型是不够准确的，因为这类边缘统计分布特征完全忽略了系数的相关依赖性；

(4)、现有降噪方法只运用传统小波，缺少多方向性。

发明内容

本发明提供了一种基于双树复小波域的地震信号降噪方法及系统，旨在解决现有的地震信号降噪方法忽略了系数间的相关依赖性、未建立小波域统计模型、以及只运用传统小波，缺少多方向性的技术问题。

本发明是这样实现的，一种基于双树复小波域的地震信号降噪方法，包括：

步骤a：建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；

步骤b：对同一方向实部小波系数或虚部小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；

步骤c：利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。

本发明实施例采取的技术方案还包括：在所述步骤a中，所述含噪地震信号模型为：

g＝x+v

其中，x为原始地震信号，v是服从分布为的高斯白噪声，g为含噪地震信号。

本发明实施例采取的技术方案还包括：在所述步骤a中，所述对含噪地震信号进行双树复小波域分解变换还包括：在双树复小波域中，含噪地震信号模型为：

y＝w+n

其中，y为含噪小波系数，w为原始地震信号小波系数，n为噪声的小波系数。

本发明实施例采取的技术方案还包括：在所述步骤a中，所述对含噪地震信号进行双树复小波域分解变换具体为：采用双树算法对地震信号进行二维双树复小波分解变换，二维双树复小波变换通过一维双树复小波的张量积得到；变换过程中，各尺度分别生成2个低频子带系数和方向分别为±15°、±45°和±75°的6个不同高频子带系数，且每个方向有实部小波系数与虚部小波系数两个部分。

本发明实施例采取的技术方案还包括：在所述步骤b中，所述双变量模型具体为：

$\left\{\begin{array}{c}{y}_1=w_1+n_1\\ m=u+n_0\end{array}\right.$

其中，含噪地震信号小波系数的模y₁是实部小波系数，y₂为对应的虚部小波系数，u为原信号小波系数的模，n₀为噪声小波系数的模。

本发明实施例采取的技术方案还包括：所述阈值函数的计算方式为：通过估计值计算公式计算得到双变量模型中的参数估计，利用参数估计得到阈值函数中各参数；所述参数估计包括待定参数t、噪声小波系数方差和估计信号的方差σ²；所述噪声小波系数方差的估计值计算公式为

${\widehat{\mathrm{\σ}}}_n^2=\frac{median\left(\right|y_i\left|\right)}{0.6745}$

其中y_i为高频子带系数；

所述估计信号的方差σ²的估计值计算公式为：

${\widehat{\mathrm{\σ}}}^2=max({\widehat{\mathrm{\σ}}}_y^2-{\widehat{\mathrm{\σ}}}_n^2,0)$

其中，为含噪观测值的方差，含噪观测值的方差的估计值计算公式为：

${\widehat{\mathrm{\σ}}}_y^2=\frac{1}{N^2}{\mathrm{\Σ}}_{y_j\∈U}{y}_j$

其中，N是邻域U的大小；

所述参数t的最优范围为[2,4]。

本发明实施例采取的另一技术方案为：一种基于双树复小波域的地震信号降噪系统，包括双树复小波分解模块、阈值处理模块和系数重构模块；

所述双树复小波分解模块用于建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；

所述阈值处理模块用于对同一方向实部小波系数或虚部小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；

所述系数重构模块用于利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。

本发明实施例采取的技术方案还包括：所述双树复小波分解模块对含噪地震信号进行双树复小波域分解变换具体为：采用双树算法对地震信号进行二维双树复小波分解变换，二维双树复小波变换通过一维双树复小波的张量积得到；变换过程中，各尺度分别生成2个低频子带系数和方向分别为±15°、±45°和±75°的6个不同高频子带系数，且每个方向有实部小波系数与虚部小波系数两个部分。

本发明实施例采取的技术方案还包括：所述双变量模型具体为：

$\left\{\begin{array}{c}{y}_1=w_1+n_1\\ m=u+n_0\end{array}\right.$

其中，含噪地震信号小波系数的模y₁是实部小波系数，y₂为对应的虚部小波系数，u为原信号小波系数的模，n₀为噪声小波系数的模。

本发明实施例采取的技术方案还包括：所述阈值处理模块计算阈值函数的计算方式为：通过估计值计算公式计算得到双变量模型中的参数估计，利用参数估计得到阈值函数中各参数；所述参数估计包括待定参数t、噪声小波系数方差和估计信号的方差σ²；所述噪声小波系数方差的估计值计算公式为

${\widehat{\mathrm{\σ}}}_n^2=\frac{median\left(\right|y_i\left|\right)}{0.6745}$

其中y_i为高频子带系数；

所述估计信号的方差σ²的估计值计算公式为：

${\widehat{\mathrm{\σ}}}^2=max({\widehat{\mathrm{\σ}}}_y^2-{\widehat{\mathrm{\σ}}}_n^2,0)$

其中，为含噪观测值的方差，含噪观测值的方差的估计值计算公式为：

${\widehat{\mathrm{\σ}}}_y^2=\frac{1}{N^2}{\mathrm{\Σ}}_{y_j\∈U}{y}_j$

其中，N是邻域U的大小；

所述参数t的最优范围为[2,4]。

本发明实施例的基于双树复小波域的地震信号降噪方法及系统通过对含噪信号进行双树复小波变换，考虑实部(或虚部)系数与对应的模之间的相关性，对实部(或虚部)小波系数与对应的模建立双变量模型进行降噪，本发明有效地压制了随机噪声的同时较好地保持了有效信号的波形特征，地震信号整体变得清晰，提高了信号质量。

附图说明

图1是本发明实施例的基于双树复小波域的地震信号降噪方法的流程图；

图2是一维双树复小波分解变换示意图；

图3是参数t对信噪比(SNR)的影响效果图；其中，图3(a)是噪声方差为10时，参数t对SNR的影响效果图；图3(b)是噪声方差为20时，参数t对SNR的影响效果图；图3(c)是噪声方差为30时，参数t对SNR的影响效果图；

图4是降噪效果图及残差图；其中，图4(a)为原始的合成地震记录；图4(b)为加入方差为20的高斯白噪声后的地震记录；图4(c)为采用本发明技术降噪结果；图4(d)为降噪结果与原始地震记录的残差图；

图5是实际地震信号降噪效果图；其中，图5(a)为实际地震信号；图5(b)为本发明降噪后的地震信号；图5(c)为压制的噪声；

图6是本发明实施例的基于双树复小波域的地震信号降噪系统的结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请参阅图1，是本发明实施例的基于双树复小波域的地震信号降噪方法的流程图。本发明实施例的基于双树复小波域的地震信号降噪方法包括以下步骤：

步骤100：建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；

在步骤100中，假定x为原始地震信号，v是服从分布为的高斯白噪声，g为含噪地震信号，则含噪地震信号模型可表示为：

g＝x+v

在给定含噪地震信号g的前提下，地震信号降噪的目的是根据一定的准则尽可能地恢复原始地震信号x；

在双树复小波域中，含噪地震信号模型可以表示为：

y＝w+n

其中，y为含噪小波系数，w为原始地震信号小波系数，n为噪声的小波系数；

本发明采用双树算法对地震信号进行二维双树复小波分解变换；双树复小波分解变换包含两个平行的分解树，需用到两个独立的实小波来分别计算实部小波系数与虚部小波系数，实部和虚部滤波器需满足希尔伯特变换；双树结构和特殊滤波器的选择使变换具有近似平移不变性和频率无偏性。

二维双树复小波变换可通过一维双树复小波的张量积得到，分解时相当于分别对行和列进行一维双树复小波变换。具体请一并参阅图2，是一维双树复小波分解变换示意图。其中g0(n)和h0(n)分别是树A分解过程的低频和高频滤波器对，g1(n)和h1(n)分别是树B分解过程的低频和高频滤波器对，↓2表示下采样操作。变换过程中，各尺度分别生成2个低频子带系数和方向分别为±15°、±45°和±75°的6个不同高频子带系数，且每个方向有实部小波系数与虚部小波系数两个部分，对上一尺度的低频子带系数再进行分解操作可达到多尺度分解的目的，最后得到不同尺度不同方向的含噪小波系数。

步骤200：对同一方向实部(或虚部)小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；

在步骤200中，考虑同一方向实部(或虚部)小波系数与其对应的模之间的依赖关系，对其建立双变量模型。假设：

$\left\{\begin{array}{c}{y}_1=w_1+n_1\\ m=u+n_0\end{array}\right.$

其中，含噪地震信号小波系数的模y₁是实部小波系数，y₂为对应的虚部小波系数，u为原信号小波系数的模，n₀为噪声小波系数的模；

假设原地震信号实部系数与模的双变量联合概率密度函数为：

$p(w_1,u)=\frac{t^2}{2{\mathrm{\π\σ}}^2}\·\exp (-\frac{t}{\mathrm{\σ}}\sqrt{w_1^2+u^2})$

假设噪声向量(n₁,n₀)是服从高斯分布，得出实部小波系数w₁的最大后验概率估计值为：

${\hat{w}}_1=\frac{{(\sqrt{y_1^2+m^2}-\frac{{\mathrm{t\σ}}_n^2}{\mathrm{\σ}})}_{+}}{\sqrt{y_1^2+m^2}}\·y_1$

同理虚部小波系数w₂最大后验概率估计值为：

${\hat{w}}_2=\frac{{(\sqrt{y_2^2+m^2}-\frac{{\mathrm{t\σ}}_n^2}{\mathrm{\σ}})}_{+}}{\sqrt{y_2^2+m^2}}\·y_2$

其中t是待定参数。

在以上计算推导过程中，为了计算得到估计值及必须先通过计算得到模型中的未知待定参数t、噪声小波系数方差和估计信号的方差即参数估计。利用参数估计得到阈值函数中各参数，再利用阈值函数对小波域中的含噪小波系数进行阈值处理，得到不含噪小波系数的估计值。

由于噪声方差在各个高频子带中不相同且无法准确估计，采用中值估计求出噪声小波系数方差估计值计算公式为

${\widehat{\mathrm{\σ}}}_n^2=\frac{median\left(\right|y_i\left|\right)}{0.6745}$

其中y_i为高频子带系数；

经过变换处理后，假设含噪观测值的方差不含噪观测值的方差σ²和噪声小波系数方差三者之间满足

${\mathrm{\σ}}_y^2={\mathrm{\σ}}_n^2+{\mathrm{\σ}}^2$

方差的估计值计算公式为

${\widehat{\mathrm{\σ}}}_y^2=\frac{1}{N^2}{\mathrm{\Σ}}_{y_j\∈U}{y}_j$

其中N是邻域U的大小；接着可得σ²的估计值计算公式为：

${\widehat{\mathrm{\σ}}}^2=max({\widehat{\mathrm{\σ}}}_y^2-{\widehat{\mathrm{\σ}}}_n^2,0)$

在本发明实施例中，参数t的确定通过仿真实验，确定其最优范围为[2,4]，实际处理中可酌情选择。具体请一并参阅图3，是参数t对信噪比(SNR)的影响效果图；其中，图3(a)是噪声方差为10时，参数t对SNR的影响效果图；图3(b)是噪声方差为20时，参数t对SNR的影响效果图；图3(c)是噪声方差为30时，参数t对SNR的影响效果图。图中的5条曲线分别表示了5个不同地震数据降噪后SNR随着参数t取值变化而变化。降噪效果先是随着参数t取值的增大而逐渐上升，在[2,4]间达到最大，即在此处取得较好的降噪效果，然后为参数t取值增大而SNR在减小；对于加入了不同大小方差噪声的不同地震数据皆有这一特征。其中，SNR的公式定义如下：

$SNR=10log\frac{{\mathrm{\Σ\Σx}}_{ij}^2}{\mathrm{\Σ}\mathrm{\Σ}{(x_{ij}-y_{ij})}^2}$

其中，x＝(x_ij)_n×m、y＝(y_ij)_n×m分别为原始信号和降噪后的信号。SNR的数值越大，就表明降噪得到的信号与原信号具有更小的差异。

步骤300：利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。

为了进一步说明本发明的有效性和准确性，具体请一并参阅图4和图5，图4是降噪效果图及残差图；其中，图4(a)为原始的合成地震记录；图4(b)为加入方差为20的高斯白噪声后的地震记录；图4(c)为采用本发明技术降噪结果；图4(d)为降噪结果与原始地震记录的残差图；观察图4(c)知，本发明得到的降噪地震信号，有效压制随机噪声的同时较好地保持了有效信号的波形特征，同时视觉效果也比较好。

图5是实际地震信号降噪效果图，其中，图5(a)为实际地震信号；图5(b)为本发明降噪后的地震信号；图5(c)为压制的噪声；观察图5(a)和图5(b)知，本发明有效地压制了随机噪声，地震信号整体变得清晰，提高了信号质量；观察图5(c)得出所压制噪声多为随机噪声。

请参阅图6，是本发明实施例的基于双树复小波域的地震信号降噪系统的结构示意图。本发明实施例的基于双树复小波域的地震信号降噪系统包括双树复小波分解模块、阈值处理模块和系数重构模块；其中，

双树复小波分解模块用于建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；其中，假定x为原始地震信号，v是服从分布为的高斯白噪声，g为含噪地震信号，则含噪地震信号模型可表示为：

g＝x+v

在给定含噪地震信号g的前提下，地震信号降噪的目的是根据一定的准则尽可能地恢复原始地震信号x；

在双树复小波域中，含噪地震信号模型可以表示为：

y＝w+n

其中，y为含噪小波系数，w为原始地震信号小波系数，n为噪声的小波系数；

本发明采用双树算法对地震信号进行二维双树复小波分解变换；双树复小波分解变换包含两个平行的分解树，需用到两个独立的实小波来分别计算实部小波系数与虚部小波系数，实部和虚部滤波器需满足希尔伯特变换；双树结构和特殊滤波器的选择使变换具有近似平移不变性和频率无偏性。

二维双树复小波变换可通过一维双树复小波的张量积得到，分解时相当于分别对行和列进行一维双树复小波变换。具体请一并参阅图2，是一维双树复小波分解变换示意图。其中g0(n)和h0(n)分别是树A分解过程的低频和高频滤波器对，g1(n)和h1(n)分别是树B分解过程的低频和高频滤波器对，↓2表示下采样操作。变换过程中，各尺度分别生成2个低频子带系数和方向分别为±15°、±45°和±75°的6个不同高频子带系数，且每个方向有实部小波系数与虚部小波系数两个部分，对上一尺度的低频子带系数再进行分解操作可达到多尺度分解的目的，最后得到不同尺度不同方向的含噪小波系数。

阈值处理模块用于对同一方向实部(或虚部)小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；其中，考虑同一方向实部(或虚部)小波系数与其对应的模之间的依赖关系，对其建立双变量模型。假设：

$\left\{\begin{array}{c}{y}_1=w_1+n_1\\ m=u+n_0\end{array}\right.$

其中，含噪地震信号小波系数的模y₁是实部小波系数，y₂为对应的虚部小波系数，u为原信号小波系数的模，n₀为噪声小波系数的模。

假设原地震信号实部系数与模的双变量联合概率密度函数为：

$p(w_1,u)=\frac{t^2}{2{\mathrm{\π\σ}}^2}\·\exp (-\frac{t}{\mathrm{\σ}}\sqrt{w_1^2+u^2})$

假设噪声向量(n₁,n₀)是服从高斯分布，得出实部小波系数w₁的最大后验概率估计值为：

${\hat{w}}_1=\frac{{(\sqrt{y_1^2+m^2}-\frac{{\mathrm{t\σ}}_n^2}{\mathrm{\σ}})}_{+}}{\sqrt{y_1^2+m^2}}\·y_1$

同理虚部小波系数w₂最大后验概率估计值为：

${\hat{w}}_2=\frac{{(\sqrt{y_2^2+m^2}-\frac{{\mathrm{t\σ}}_n^2}{\mathrm{\σ}})}_{+}}{\sqrt{y_2^2+m^2}}\·y_2$

其中t是待定参数。

在以上计算推导过程中，为了计算得到估计值及必须先通过计算得到模型中的未知待定参数t、噪声小波系数方差和估计信号的方差σ²，即参数估计。利用参数估计得到阈值函数中各参数，再利用阈值函数对小波域中的含噪小波系数进行阈值处理，得到不含噪小波系数的估计值。

由于噪声方差在各个高频子带中不相同且无法准确估计，采用中值估计求出噪声小波系数方差估计值计算公式为

${\widehat{\mathrm{\σ}}}_n^2=\frac{median\left(\right|y_i\left|\right)}{0.6745}$

其中y_i为高频子带系数；

经过变换处理后，假设含噪观测值的方差不含噪观测值的方差σ²和噪声小波系数方差三者之间满足

${\mathrm{\σ}}_y^2={\mathrm{\σ}}_n^2+{\mathrm{\σ}}^2$

方差的估计值计算公式为

${\widehat{\mathrm{\σ}}}_y^2=\frac{1}{N^2}{\mathrm{\Σ}}_{y_j\∈U}{y}_j$

其中N是邻域U的大小；接着可得σ²的估计值计算公式为：

${\widehat{\mathrm{\σ}}}^2=max({\widehat{\mathrm{\σ}}}_y^2-{\widehat{\mathrm{\σ}}}_n^2,0)$

在本发明实施例中，参数t的确定通过仿真实验，确定其最优范围为[2,4]，实际处理中可酌情选择。具体请一并参阅图3，是参数t对信噪比(SNR)的影响效果图；其中，图3(a)是噪声方差为10时，参数t对SNR的影响效果图，图3(b)是噪声方差为20时，参数t对SNR的影响效果图，图3(c)是噪声方差为30时，参数t对SNR的影响效果图。图中的5条曲线分别表示了5个不同地震数据降噪后SNR随着参数t取值变化而变化。降噪效果先是随着参数t取值的增大而逐渐上升，在[2,4]间达到最大，即在此处取得较好的降噪效果，然后为参数t取值增大而SNR在减小；对于加入了不同大小方差噪声的不同地震数据皆有这一特征。其中，SNR的公式定义如下：

$SNR=10log\frac{{\mathrm{\Σ\Σx}}_{ij}^2}{\mathrm{\Σ}\mathrm{\Σ}{(x_{ij}-y_{ij})}^2}$

其中，x＝(x_ij)_n×m、y＝(y_ij)_n×m分别为原始信号和降噪后的信号。SNR的数值越大，就表明降噪得到的信号与原信号具有更小的差异。

系数重构模块用于利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。

本发明实施例的基于双树复小波域的地震信号降噪方法及系统通过对含噪信号进行双树复小波变换，考虑实部(或虚部)系数与对应的模之间的相关性，对实部(或虚部)小波系数与对应的模建立双变量模型进行降噪，本发明有效地压制了随机噪声的同时较好地保持了有效信号的波形特征，地震信号整体变得清晰，提高了信号质量。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于双树复小波域的地震信号降噪方法，包括：

步骤a：建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；

步骤b：对同一方向实部小波系数或虚部小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；

步骤c：利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。

2.根据权利要求1所述的基于双树复小波域的地震信号降噪方法，其特征在于，在所述步骤a中，所述含噪地震信号模型为：

g＝x+v

其中，x为原始地震信号，v是服从分布为的高斯白噪声，g为含噪地震信号。

3.根据权利要求1所述的基于双树复小波域的地震信号降噪方法，其特征在于，在所述步骤a中，所述对含噪地震信号进行双树复小波域分解变换还包括：在双树复小波域中，含噪地震信号模型为：

y＝w+n

其中，y为含噪小波系数，w为原始地震信号小波系数，n为噪声的小波系数。

4.根据权利要求1所述的基于双树复小波域的地震信号降噪方法，其特征在于，在所述步骤a中，所述对含噪地震信号进行双树复小波域分解变换具体为：采用双树算法对地震信号进行二维双树复小波分解变换，二维双树复小波变换通过一维双树复小波的张量积得到；变换过程中，各尺度分别生成2个低频子带系数和方向分别为±15°、±45°和±75°的6个不同高频子带系数，且每个方向有实部小波系数与虚部小波系数两个部分。

5.根据权利要求1所述的基于双树复小波域的地震信号降噪方法，其特征在于，在所述步骤b中，所述双变量模型具体为：

$\left\{\begin{array}{c}{y}_1=w_1+n_1\\ m=u+n_0\end{array}\right.$

其中，含噪地震信号小波系数的模y₁是实部小波系数，y₂为对应的虚部小波系数，u为原信号小波系数的模，n₀为噪声小波系数的模。

6.根据权利要求1所述的基于双树复小波域的地震信号降噪方法，其特征在于，所述阈值函数的计算方式为：通过估计值计算公式计算得到双变量模型中的参数估计，利用参数估计得到阈值函数中各参数；所述参数估计包括待定参数t、噪声小波系数方差和估计信号的方差σ²；所述噪声小波系数方差的估计值计算公式为

${\widehat{\mathrm{\σ}}}_n^2=\frac{median\left(\right|y_i\left|\right)}{0.6745}$

其中y_i为高频子带系数；

所述估计信号的方差σ²的估计值计算公式为：

${\widehat{\mathrm{\σ}}}^2=max({\widehat{\mathrm{\σ}}}_y^2-{\widehat{\mathrm{\σ}}}_n^2,0)$

其中，为含噪观测值的方差，含噪观测值的方差的估计值计算公式为：

${\widehat{\mathrm{\σ}}}_y^2=\frac{1}{N^2}{\mathrm{\Σ}}_{y_j\∈U}{y}_j$

其中，N是邻域U的大小；

所述参数t的最优范围为[2,4]。

7.一种基于双树复小波域的地震信号降噪系统，其特征在于，包括双树复小波分解模块、阈值处理模块和系数重构模块；

所述双树复小波分解模块用于建立含噪地震信号模型，对含噪地震信号进行双树复小波域分解变换，得到每一尺度每个方向的含噪小波系数；

所述阈值处理模块用于对同一方向实部小波系数或虚部小波系数与其对应的模建立双变量模型得到阈值函数，利用阈值函数对双树复小波域中的含噪小波系数进行双变量阈值处理，得到不含噪小波系数的估计值；

所述系数重构模块用于利用双树复小波逆变换对不含噪小波系数的估计值进行重构得到降噪后的原始地震信号。

8.根据权利要求7所述的基于双树复小波域的地震信号降噪系统，其特征在于，所述双树复小波分解模块对含噪地震信号进行双树复小波域分解变换具体为：采用双树算法对地震信号进行二维双树复小波分解变换，二维双树复小波变换通过一维双树复小波的张量积得到；变换过程中，各尺度分别生成2个低频子带系数和方向分别为±15°、±45°和±75°的6个不同高频子带系数，且每个方向有实部小波系数与虚部小波系数两个部分。

9.根据权利要求7所述的基于双树复小波域的地震信号降噪系统，其特征在于，所述双变量模型具体为：

$\left\{\begin{array}{c}{y}_1=w_1+n_1\\ m=u+n_0\end{array}\right.$

其中，含噪地震信号小波系数的模y₁是实部小波系数，y₂为对应的虚部小波系数，u为原信号小波系数的模，n₀为噪声小波系数的模。

10.根据权利要求7所述的基于双树复小波域的地震信号降噪方法，其特征在于，所述阈值处理模块计算阈值函数的计算方式为：通过估计值计算公式计算得到双变量模型中的参数估计，利用参数估计得到阈值函数中各参数；所述参数估计包括待定参数t、噪声小波系数方差和估计信号的方差σ²；所述噪声小波系数方差的估计值计算公式为

${\widehat{\mathrm{\σ}}}_n^2=\frac{median\left(\right|y_i\left|\right)}{0.6745}$

其中y_i为高频子带系数；

所述估计信号的方差σ²的估计值计算公式为：

${\widehat{\mathrm{\σ}}}^2=max({\widehat{\mathrm{\σ}}}_y^2-{\widehat{\mathrm{\σ}}}_n^2,0)$

其中，为含噪观测值的方差，含噪观测值的方差 ${\mathrm{\σ}}_y^2$ 的估计值计算公式为：

${\widehat{\mathrm{\σ}}}_y^2=\frac{1}{N^2}{\mathrm{\Σ}}_{y_j\∈U}{y}_j$

其中，N是邻域U的大小；

所述参数t的最优范围为[2,4]。