[发明专利]一种基于CEEMDAN与Savitzky-Golay滤波的地震信号噪声抑制方法

权利要求书

1.一种基于CEEMDAN与Savitzky-Golay滤波的地震信号噪声抑制方法,其特征在于，该方法包括：

第1步：提供待处理的含噪地震信号，作为整个系统输入的原始含噪地震信号s(t)(t＝1,...,L),L是地震信号的采样点数；

第2步：对输入的含噪地震信号s(t)进行CEEMDAN算法分解，得到多个IMF分量和残差分量r_N(t)；

第3步：将IMF分量从高频到低频依次排列，通过计算与原始含噪地震信号s(t)的相关系数确定噪声与信号的分界M，其中M为第一个互相关系数取局部极小值时对应的IMF分量的下标加1；

第4步：Savitzky-Golay滤波：找到IMF分量中主要包含噪声信号的IMF分量将含噪声分量分别作为输入信号，送入SG滤波器中做进一步的平滑，得到滤波后的信号

第5步：数据重构得到降噪后的地震信号

2.按照权利要求1所述的方法,其特征在于，其中第2步中包括：CEEMDAN算法通过向待分解的含噪地震信号s(t)中添加高斯白噪声生成多个新的信号，再分别对新的信号应用EMD分解求取各信号的第一个IMF分量，并求多个IMF分量的平均值，将其作为CEEMDAN分解的第一个IMF分量，将含噪地震信号s(t)减去IMF分量得到的值作为新的待分解信号重复上述过程，计算出所有IMF分量。

3.按照权利要求2所述的方法,其特征在于，具体包括：将E_j()定义为EMD分解的第j个IMF分量，令n_i(t)为零均值单位方差白噪声；如果原始含噪地震信号s(t)是目标信号，则CEEMDAN算法步骤如下：

a、对于所有i＝1,...,I,通过EMD分解每个s_i(t)＝s(t)+ε₀n_i(t)得到CEEMDAN分解的第一个IMF分量，并将第一个IMF分量定义为:

b、在j＝1时，计算第一个残差:

c、对于每个i＝1,...,I，通过EMD分解每个r₁(t)+ε₁E₁(n_i(t))，直到找到它的第一个IMF分量为止，并且定义CEEMDAN分解的第二个IMF分量为:

d、对于j＝2,3,...,N，计算第j个残差:

e、对于每个i＝1,...,I，通过EMD分解每个r_j(t)+ε_jE_j(n_i(t))，直到找到它的第一个IMF分量为止，并且定义CEEMDAN分解第(j+1)个IMF分量为:

f、对于下一个j跳转到步骤d；重复步骤d到f，直到获得的残差不能再由EMD进一步分解，满足IMF标准或者小于三个局部极值；最终的残差满足:其中N是模式的总数；输入信号表示为下式：

4.按照权利要求1所述的方法,其特征在于，其中第3步包括：

相关系数可表示为：

其中s表示原始信号，表示第j个IMF分量，L为采样点个数，且代表第j个IMF分量与原始输入信号的相关程度。

5.按照权利要求1所述的方法,其特征在于，第4步中，Savitzky-Golay滤波具体过程：

SG滤波器采用的平滑公式式中每一个测量值乘平滑系数h_i是为了减小平滑对有用信号的影响，h_i/H应用最小二乘原理与多项式计算得到；

设n＝2m+1为SG滤波器窗口，测量值为x＝(-m,-m+2,…,0,1,…,m-1,m)，应用k-1次多项式将窗口内的数据做拟合处理：y＝a₀+a₁x+a₂x²+…+a_k-1x^k-1，n个方程构成线性方程组，若想要方程组有解则应满足n≥k，通常选择n＞k，并应用最小二乘法拟合进一步确定拟合参数A，有：

转化成矩阵形式为：Y_(2m+1)×1＝X_(2m+1)×k·A_K×1+E_(2m+1)×1，A的最小二乘积为Y的滤波值其中D＝X·(X^T·X)^-1·X^T；将含噪声分量分别作为输入信号代入到上述的过程中替换测量值x进行滤波。