CN108957552B - 基于ss-pca的地震数据海浪噪声压制方法 - Google Patents

Abstract

本发明涉及一种基于SS‑PCA的地震数据海浪噪声压制方法，是利用同步挤压小波变换分辨率高的优点，估算海浪噪声的主要频率，并利用主成分分析法压制海浪噪声。经试验，不仅提高信号信噪比且有效保留了与噪声同频段的信号能量。与现有的带通滤波方法相比，本方法在压制海浪噪声的同时，能更好地保留低频信号特征，有利于改善受到强海浪噪声干扰的地震信号质量，即使在频谱混叠情况下也能有效压制噪声，有利于提高震相拾取的可靠性和准确性，有利于减小地震群速度和相速度的估计误差，对于天然地震到时拾取、震中定位和震级估计，特别是面波震级估计具有误差小、准确度高的特点。提高了天然地震定位精度和减小了震级的不准确度。

Description

基于SS-PCA的地震数据海浪噪声压制方法

技术领域：

本发明涉及一种基于SS-PCA的地震数据海浪噪声压制方法，特别是适用于受海浪噪声干扰严重的天然地震信号压噪。

背景技术：

天然地震事件三要素(发震时间、震中方位和震级)的估算常由于强噪声的存在而产生较大的误差，经过广泛分析遍布世界的地震台站发现，天然地震噪声在不同频段噪声分布不均匀，其主要能量一般集中于1Hz以下，这种噪声又被称为地脉动信号，它主要来源于海洋波浪与固体地球之间的相互作用。

目前传统去噪方式有两大类。第一类为频域去噪，代表方法有傅里叶频谱分析、功率谱分析、中值滤波和高斯滤波等方法，这些方法都建立在背景噪声为平稳随机信号的基础上，然而经过观察分析可以发现很多台站的背景噪声不完全遵守这个前提条件，因此其去噪同时常常损失很多信号细节，甚至导致信号能量大幅衰减，令震级估算结果准确度低。第二类方法为时-频联合去噪，这类方法将会更精细的分析信号与噪声的频率特征随时间变化的情况，以此为前提实现压制噪声的目的。代表方法有短时傅里叶变换(STFT)、Wigner-will分布和小波变换(WT)等，然而此类方法分辨率相对较低且无法解决在信号与噪声存在频谱混叠情况下的噪声压制问题。

目前由于傅里叶变换频谱分析方法由于具有方法简单、处理速度快且处理效果在一定程度上可接受的优点被广泛应用于天然地震信号噪声压制的工作中，但是由于瑞雷波、深度震相等震相其主频较低，上述方法处理后使得瑞雷波震级可信度降低、震源深度估算误差大，特别对于受海浪噪声严重影响的地震台站，误差更严重。

地震监测需要确定地震事件的三要素，包括发震时间、震中方位和震级，其精度依赖于地震数据的质量，因此当背景噪声较强时会引起较大的估算误差，对于中小型地震来说更是如此。常规的地震信号噪声常常被假定为平稳随机噪声，通常在低频段采用一个带通滤波器达到压制噪声的目的，这种方法简单快捷，并且当信号与噪声的频带不混叠时去噪效果理想，但事实上天然地震信号和噪声存在频带混叠的现象，特别是较强的海浪噪声频带更是落在地震信号的频带范围内。经带通滤波后面波能量被大幅削弱，而且震相到时出现较大误差，同时对面波能量的大幅衰减也导致面波震级估算误差。

发明内容：

本发明的目的就在于针对上述现有技术的不足，提供一种特别适用于受海浪噪声干扰严重的天然地震信号压噪的基于SS-PCA的地震数据海浪噪声压制方法。

本发明的发明思想是，根据地震数据和海水噪声的特点，首先利用同步挤压小波变换方法确定稳定海水噪声的频率个数；再利用海水噪声相关性强而地震信号相关性弱的特点结合主成分分析方法实现压制噪声的目的。

基于SS-PCA的地震数据海浪噪声压制方法，包括以下步骤：

a、读取待处理的受强海浪噪声干扰的地震信号x[l],l∈[1,L]，L为地震信号总采样点数，读取震前具有一定长度的背景噪声n[h],h∈[1,H]，H为噪声总采样点数，噪声长度取100～200s之间，信号采样率为FS；

b、计算n[h]的小波系数，小波变换时尺度因子采样率为fs，在尺度因子a_m,m∈[1,M]和平移因子τ_k,k∈[1,K]处定义小波系数W_n(a_m,τ_k)为

其中ψ(t)为给定母小波函数，ψ(t)满足平方可积且无直流分量的条件，“*”表示复共轭，M和K分别为尺度因子和平移因子的个数，由信号采样点数及信号采样率决定；

c、定义瞬时频率

其中i为虚数单位；

d、在ω_m＝m×fs/M处的同步挤压变换定义如下

其中Δω＝fs/M，

e、各频率同步挤压系数的均值

方差

f、定义

其中f_m＝ω_m/2π，分析满足g(f_m)≥0且小于1Hz的频率的个数，记为α；

g、对X利用主成分分析(PCA)方法进行噪声压制，首先构造Hankel矩阵,当L为偶数时，令

当L为奇数时，令

h、计算Hankel矩阵H_x的协方差矩阵Γ_x，并进行特征值分解：

其中特征值矩阵Λ＝Diag[λ₁,λ₂,...,λ_Q]，且特征根按着由大到小的顺序排列，“·”为矩阵乘法，当L为偶数时，Q＝L/2，当L为奇数时，Q＝(L+1)/2，R为特征矩阵，其与特征值一一对应；最终得到主成分矩阵

Φ＝H_x·R^T (11)

i.将Φ第1至2α行元素置零，得到重构矩阵Φ′，则重构后的Hankel矩阵H′_x为：

j、取H′_x的第一行以及最后一列的第二行至最后一行数据重构，抑制海浪噪声后的信号x′[l]＝[x′[1],x′[2],...,x′[L]]，x′[l]即为x[l]经SS-PCA方法压制海浪噪声后的结果。

有益效果：经试验，本发明能够有利于改善受到强海浪噪声干扰的地震信号质量，即使在频谱混叠情况下也能有效压制噪声，保留信号特征，有利于提高震相拾取的可靠性和准确性，有利于减小地震群速度和相速度的估计误差，对于天然地震到时拾取、震中定位和震级估计，特别是面波震级估计具有误差小、准确度高的特点。提高了天然地震定位精度和减小了震级的不准确度。

附图说明：

图1是SS-PCA与带通滤波去噪效果对比；

a.理想地震信号；b.含海浪噪声地震信号；

c.SS-PCA方法去噪效果；d.带通滤波去噪效果。

图2是某台站记录的1Hz以下的环境噪声频率特性分析；

具体实施方式：

下面结合附图和实施例对本发明提出的一种基于SS-PCA的地震数据海浪噪声压制方法做进一步的详细说明。

基于同步挤压小波变换与主成分分析方法的地震海水噪声压制方法，包括如下步骤：

a、某台站记录到的一次事件信号，由于其噪声小，具有很高的信噪比，故在本发明的实例中使用此次事件z分量s[l],l∈[1,L]作为理想信号,如图1(a)所示。截取在此事件发生与s[l]长度相同的环境噪声n[l]，与s[l]相加可得含海浪噪声地震信号x[l]，在本例中信号采样率FS＝100，L＝12000，如图1(b),此时信噪比SNR₀＝-3.3508，x[l]与s[l]相似性Cor₀＝0.5264。

x[l]＝s[l]+6*n[l] (1)

b、计算n[l]的小波系数，小波变换时尺度因子采样率为fs＝0.025，在尺度因子a_m,m∈[1,M]和平移因子τ_k,k∈[1,K]处定义小波系数W_n(a_m,τ_k)为

其中ψ(t)为给定母小波函数，ψ(t)满足平方可积且无直流分量的条件，“*”表示复共轭，M和K分别为尺度因子和平移因子的个数；

c、定义瞬时频率

其中i为虚数单位；

d、在ω_m＝m×fs/M处的同步挤压变换定义如下

其中Δω＝fs/M，

e、各频率同步挤压系数的均值

方差

f、定义

其中f_m＝ω_m/2π，分析满足g(f_m)≥0且小于1Hz的频率的个数，记为α，在此例中α＝6；

g、对X利用主成分分析(PCA)方法进行噪声压制，首先构造Hankel矩阵,当L为偶数时，令

h、计算Hankel矩阵H_x的协方差矩阵Γ_x，并进行特征值分解：

其中特征值矩阵Λ＝Diag[λ₁,λ₂,...,λ_Q]，且特征根按着由大到小的顺序排列，“·”为矩阵乘法，Q＝L/2，R为特征矩阵，其与特征值一一对应；最终得到主成分矩阵

Φ＝H_x·R^T (11)

i.将Φ第1至12行元素置零，得到重构矩阵Φ′，则重构后的Hankel矩阵H′_x为：

j、取H′_x的第一行以及最后一列的第二行至最后一行数据重构，抑制海浪噪声后的信号x′[l]＝[x′[1],x′[2],...,x′[L]]，x′[l]即为x[l]经SS-PCA方法压制海浪噪声后的结果,如图1(c)所示。

此时经本发明提出方法处理后的信号信噪比SNR′＝1.3433，x′[l]与s[l]相似性Cor′＝0.6446。图1(d)为采用传统带通滤波的信号x″[l]，带通范围为1～10Hz，由于滤波器的相移特性使得处理后信号信噪比及相似性都大大降低，其中SNR″＝-1.6401，Cor″＝0.0510，并且对于信号的瑞雷波等震相在带通滤波后基本不可见，而本发明提出的方法则较好的保留这些震相特征，有利于后期进一步的面波震级确定和震源深度确定。

Claims (1)

1.一种基于SS-PCA的地震海浪噪声压制方法，其特征在于，包括以下步骤：

(1)读取待处理的受强海浪噪声干扰的地震信号x[l],l∈[1,L]，L为地震信号总采样点数，读取震前具有一定长度的背景噪声n[h],h∈[1,H]，H为噪声总采样点数，噪声长度取100～200s之间，信号采样率为FS；

(2)计算n[h]的小波系数，小波变换时尺度因子采样率为fs，在尺度因子a_m,m∈[1,M]和平移因子τ_k,k∈[1,K]处定义小波系数W_n(a_m,τ_k)为

其中ψ(t)为给定母小波函数，ψ(t)满足平方可积且无直流分量的条件，“*”表示复共轭，M和K分别为尺度因子和平移因子的个数，由信号采样点数及信号采样率决定；

(3)定义瞬时频率

其中i为虚数单位；

(4)在ω_m＝m×fs/M处的同步挤压变换定义如下

其中Δω＝fs/M，

(5)各频率同步挤压系数的均值

方差

(6)定义

其中f_m＝ω_m/2π，分析满足g(f_m)≥0且小于1Hz的频率的个数，记为α；

(7)对x[l]利用主成分分析(PCA)方法进行噪声压制，首先构造Hankel矩阵,当L为偶数时，令

当L为奇数时，令

(8)计算Hankel矩阵H_x的协方差矩阵Γ_x，并进行特征值分解：

其中特征值矩阵Λ＝Diag[λ₁,λ₂,...,λ_Q]，且特征根按着由大到小的顺序排列，“·”为矩阵乘法，当L为偶数时，Q＝L/2，当L为奇数时，Q＝(L+1)/2，R为特征矩阵，其与特征值一一对应；最终得到主成分矩阵

Φ＝H_x·R^T (11)

(9)将Φ第1至2α行元素置零，得到重构矩阵Φ′，则重构后的Hankel矩阵H′_x为：

(10)取H′_x的第一行以及最后一列的第二行至最后一行数据重构，抑制海浪噪声后的信号x′[l]＝[x′[1],x′[2],...,x′[L]]，x′[l]即为x[l]经SS-PCA方法压制海浪噪声后的结果。

Images