CN115857013A - 一种利用改进welch法计算地震计自噪声的方法 - Google Patents

Abstract

本发明提供了一种利用改进welch法计算地震计自噪声的方法，包括：获取待计算自噪声的两台地震计的数据文件；对数据文件进行分段得到分段数据；对分段数据进行预处理得到预处理后的分段数据；对预处理后的分段数据依次进行加窗和内插补零得到补零后的数据段；计算补零后的数据段的平均自功率谱密度和平均互功率谱密度；获取两台地震计各自的传递函数，根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声。本发明通过对分段数据进行加窗和内插补零，可在保证频率分辨率的同时，提高了功率谱估计的平滑度，同时克服了频谱泄漏，从而提高了地震计自噪声的测试精度。

Description

一种利用改进welch法计算地震计自噪声的方法

技术领域

本发明涉及地震计技术领域，特别是涉及一种利用改进welch法计算地震计自噪声的方法。

背景技术

自噪声是宽频带地震计性能的重要指标。目前，在地震计行业标准中，地震计自噪声的测试方法主要使用两仪器法和三仪器法。这两种方法思路相似，都是利用两台或三台地震计的观测数据进行噪声功率谱的计算，由功率谱的相关性求得该地震计的自噪声。目前进行功率谱计算最常用的方法是加窗平均周期图法。

周期图法虽然是一种渐近无偏的估计，但其方差较大，无法收敛到真谱。造成这种计算结果性能不好的根本原因在于傅里叶变换的特点。傅里叶变换是把所分析的信号分解成无穷多正弦信号的叠加，这些信号的幅度、频率及相位都是固定不变的。但是随机信号的幅度、频率和相位是随机变化的，而周期图实际上是把随机信号视为确定性的信号，因此必然带来估计上的质量不高。

发明内容

为了克服现有技术的不足，本发明的目的是提供一种利用改进welch法计算地震计自噪声的方法以解决地震计自噪声测算精度低的问题。

为实现上述目的，本发明提供了如下方案：

一种利用改进welch法计算地震计自噪声的方法，包括：

获取待计算自噪声的两台地震计的数据文件；

对所述数据文件进行分段得到分段数据；

对所述分段数据进行预处理得到预处理后的分段数据；

对所述预处理后的分段数据依次进行加窗和内插补零得到补零后的数据段；

计算所述补零后的数据段的平均自功率谱密度和平均互功率谱密度；

获取两台地震计各自的传递函数，根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声。

优选地，所述对所述分段数据进行预处理得到预处理后的分段数据，包括：

对所述分段数据进行去均值处理得到去除均值后的分段数据；

对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据。

优选地，所述对所述分段数据进行去均值处理得到去除均值后的分段数据，包括：

采用公式：

/>

对所述分段数据进行去均值处理得到去除均值后的分段数据；其中，x¹ _i(n)为第i段数据，L为分段数据的段数，且1≤i≤L，N为分段数据的各段数据长度，且0≤n≤N，x² _i(n)为去除均值后的分段数据。

优选地，所述对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据，包括：

采用公式：

对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据；其中，x_i(n)为预处理后的分段数据。

优选地，计算所述补零后的数据段的平均自功率谱密度和平均互功率谱密度，包括：

对所述补零后的数据段进行傅里叶变换得到傅里叶变换后的补零数据段；

求取所述傅里叶变换后的补零数据段的自功率谱密度和互功率谱密度；

根据所述自功率谱密度和互功率谱密度得到所述补零后的数据段的平均自功率谱密度和平均互功率谱密度。

优选地，补零后的数据段的平均自功率谱密度计算方法为：

采用公式：

得到平均自功率谱密度；其中，

为补零后的数据段，X(k)为傅里叶变换后的补零数据段，Pⁱ _XX(ω)为自功率谱密度，/>

为平均自功率谱密度，d(n)为数据窗，j为虚数单位，ω为角频率，k为数据段内数据点的序号。

优选地，补零后的数据段的平均互功率谱密度计算方法为：

采用公式：

得到平均互功率谱密度；其中，Pⁱ _XY(ω)为互功率谱密度，

为平均互功率谱密度，X^*(k)为X(k)的转置，Y(k)为第二地震计经傅里叶变换后的补零数据段。

优选地，根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声，包括：

采用公式：

计算各自地震计的自噪声；其中，P₁为第一地震计的自噪声，P₂为第二地震计的自噪声，H₁为第一地震计的传递函数，H₂为第二地震计的传递函数。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供了一种利用改进welch法计算地震计自噪声的方法，与现有技术相比，本发明通过对分段数据进行加窗和内插补零，可在保证频率分辨率的同时，提高了功率谱估计的平滑度，同时克服了频谱泄漏，从而提高了地震计自噪声的测试精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的实施例中的一种利用改进welch法计算地震计自噪声的方法流程图；

图2为本发明提供的实施例中的一种利用改进welch法计算地震计自噪声的方法原理图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。

本申请的说明书和权利要求书及所述附图中的术语“第一”、“第二”、“第三”和“第四”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤、过程、方法等没有限定于已列出的步骤，而是可选地还包括没有列出的步骤，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤元。

本发明的目的是提供一种利用改进welch法计算地震计自噪声的方法以解决地震计传递函数标定精度低的问题。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

请一并参阅图1-2，一种利用改进welch法计算地震计自噪声的方法，包括：

步骤1：获取待计算自噪声的两台地震计的数据文件；

步骤2：对所述数据文件进行分段得到分段数据；

在实际应用中，对于待计算自噪声的两台地震计，选取其中一台的数据文件，截取优选的数据段x(n)。此数据段应为地震计在稳定状态下连续运行20天后的采集数据。

将x(n)进行重叠分段，即前一段信号和后一段信号有一部分是重叠的，特别是每段一半的重叠率能大大降低谱估计的方差。各段数据长度为n，共分为L段，记各段数据为x¹ _i(n)(1≤i≤L，0≤n≤N)。

对另一台地震计截取同一开始时间和终止时间的数据段，用上述同一方法进行分段，记各段数据为y¹ _i(n)(1≤i≤L，0≤n≤N)。

步骤3：对所述分段数据进行预处理得到预处理后的分段数据；

进一步的，步骤3包括：

对所述分段数据进行去均值处理得到去除均值后的分段数据；

对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据。

在本发明实施例中，为了去除数据中不表示真实地面运动的直流偏移，需对观测数据分段进行去均值及去趋势值处理，预处理后的数据记为x_i(n)和y_i(n)：

对每段x¹ _i(n)和y¹ _i(n)，以x¹ _i(n)为例，取段内平均值，

对x¹ _i(n)的各数据点进行去均值处理

x² _i(n)＝x¹ _i(n)-x_a (2)

2)去趋势值

求x² _i(n)前1/3段、后1/3段数据的平均值，计算公式为：

其中，x¹ _i(n)为第i段数据，L为分段数据的段数，且1≤i≤L，N为分段数据的各段数据长度，且0≤n≤N，x² _i(n)为去除均值后的分段数据，x_i(n)为预处理后的分段数据。

步骤4：对所述预处理后的分段数据依次进行加窗和内插补零得到补零后的数据段；

需要说明的是，本发明的数据加窗步骤为：

对步骤3中处理后的x_i(n)加窗，每段的数据窗口可以不是矩形窗口，例如可使用汉宁窗或汉明窗，数据窗记为d(n)。

加窗后的数据段为：

x^d _i(n)＝x_i(n)d(n) (6)

y^d _i(n)＝y_i(n)d(n) (7)

本发明的内插补零步骤为：

在加窗后的数据段x^d _i(n)末尾补零，补零的长度与原数据段长度相同，即在x^d _i(n)末尾添加长度为N的零向量。

补零后得到数据段：

其中，0≤k≤2N-1。

步骤5：计算所述补零后的数据段的平均自功率谱密度和平均互功率谱密度；

进一步的，步骤5包括：

对所述补零后的数据段进行傅里叶变换得到傅里叶变换后的补零数据段；

求取所述傅里叶变换后的补零数据段的自功率谱密度和互功率谱密度；

根据所述自功率谱密度和互功率谱密度得到所述补零后的数据段的平均自功率谱密度和平均互功率谱密度。

下面本发明结合具体的实施例对上述计算步骤做详细说明：

先由

傅里叶变换得到X(k)：

由X(k)求得自功率谱的估计Pⁱ _XX(ω)：

式中

是归一化因子，使用它是为了保证所得到的谱是渐近无偏估计。

上式中的

是对x⁰ _i(k)进行DFT运算，与常用的welch法相比，e^-jωk中指数

的分母N增加，即在更小的频次刻度上(即更高的频率分辨率)观察原信号数列。换言之，补零前的功率谱计算得到的是/>

处的频域值，而补零后原/>

处的频域值不变，在/>

处也有了插值，故能够更加仔细的反映出功率谱的特征，可以克服“栅栏”效应，使谱的外观平滑。

将

代入式(8)和式(9)，可求得Y(k)和Pⁱ _YY(ω)。

1、求平均自功率谱密度

对上一步中求得的Pⁱ _XX(ω)对应相加，再取平均，得到平均周期图

将Pⁱ _YY(ω)代入式(11)可得

2、计算互功率谱密度：

计算

和/>

的互功率谱：

式中

3、求平均互功率谱密度

对上一步中求得的Pⁱ _XY(ω)对应相加，再取平均，得到平均周期图

上述各式中，

为补零后的数据段，X(k)为傅里叶变换后的补零数据段，Pⁱ _XX(ω)为自功率谱密度，/>

为平均自功率谱密度，d(n)为数据窗，j为虚数单位，ω为角频率，k为数据段内数据点的序号，Pⁱ _XY(ω)为互功率谱密度，/>

为平均互功率谱密度，X^*(k)由式(8)求得，为X(k)的转置，Y(k)为式(7)经过式(8)和式(9)计算出的结果，即第二地震计经傅里叶变换后的补零数据段。

步骤6：获取两台地震计各自的传递函数，根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声。

进一步的，步骤6包括：

读入两台地震计的传递函数参数，计算各频率点的传递函数模值，两台地震计的传递函数分别记为H₁，H₂。

根据两通道法地震计自噪声测试原理，两台地震计的自噪声分别为：

通常选用两台传递函数完全相同的地震计，则上式可简化为：

其中，P₁为第一地震计的自噪声，P₂为第二地震计的自噪声，H₁为第一地震计的传递函数，H₂为第二地震计的传递函数。

在实际应用中，本发明可选用Matlab作为改进的加窗平均周期图法计算地震计自噪声的实现工具。将两台待测地震计置于监测环境优越的测试地点，采样频率设置为100Hz，选取地震计连续工作30天后的采集数据，取其中3600s进行自噪声计算。

将一台地震计中所选取的数据x(n)进行重叠分段，各段数据长度为N＝1024，共分为L＝352段，记各段数据为x¹ _i(n)(1≤i≤L，0≤n≤N)。

对另一台地震计截取同一开始时间和终止时间的数据段，用上述同一方法进行分段，记各段数据为y^l _i(n)(1≤i≤L，0≤n≤N)。最后根据对分段数据依次进行预处理、加窗和内插补零得到补零后的数据段；其中，加窗后的数据段x^d _i(n)末尾补零，补零的长度为N＝1024.

补零后得到数据段

其中0≤k≤2N-1，0≤n≤N-1。

计算所述补零后的数据段的平均自功率谱密度和平均互功率谱密度；根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明通过引入零值填充的方法，在保证频率分辨率的同时，能够对功率谱进行平滑处理。另外，加窗平均周期图法中的窗函数的引入，也对谱估计造成了“频谱泄漏”的影响。零值填充方法对于窗函数引起的“频谱泄漏”，造成频谱中出现难以确认的谱峰等也有消除作用。

综上，本发明改进的加窗平均周期图法，在保证频率分辨率的同时，提高了功率谱估计的平滑度，同时克服了频谱泄漏，从而提高了地震计自噪声的测试精度。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的方法而言，由于其与实施例公开的装置相对应，所以描述的比较简单，相关之处参见装置部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种利用改进welch法计算地震计自噪声的方法，其特征在于，包括：

获取待计算自噪声的两台地震计的数据文件；

对所述数据文件进行分段得到分段数据；

对所述分段数据进行预处理得到预处理后的分段数据；

对所述预处理后的分段数据依次进行加窗和内插补零得到补零后的数据段；

计算所述补零后的数据段的平均自功率谱密度和平均互功率谱密度；

获取两台地震计各自的传递函数，根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声。

2.根据权利要求1所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，所述对所述分段数据进行预处理得到预处理后的分段数据，包括：

对所述分段数据进行去均值处理得到去除均值后的分段数据；

对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据。

3.根据权利要求2所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，所述对所述分段数据进行去均值处理得到去除均值后的分段数据，包括：

采用公式：

对所述分段数据进行去均值处理得到去除均值后的分段数据；其中，x¹ _i(n)为第i段数据，L为分段数据的段数，且1≤i≤L，N为分段数据的各段数据长度，且0≤n≤N，x² _i(n)为去除均值后的分段数据。

4.根据权利要求3所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，所述对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据，包括：

采用公式：

对所述去除均值后的分段数据进行去趋势值处理得到预处理后的分段数据；其中，x_i(n)为预处理后的分段数据。

5.根据权利要求4所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，计算所述补零后的数据段的平均自功率谱密度和平均互功率谱密度，包括：

对所述补零后的数据段进行傅里叶变换得到傅里叶变换后的补零数据段；

求取所述傅里叶变换后的补零数据段的自功率谱密度和互功率谱密度；

根据所述自功率谱密度和互功率谱密度得到所述补零后的数据段的平均自功率谱密度和平均互功率谱密度。

6.根据权利要求5所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，补零后的数据段的平均自功率谱密度计算方法为：

采用公式：

得到平均自功率谱密度；其中，

为补零后的数据段，X(k)为傅里叶变换后的补零数据段，Pⁱ _XX(ω)为自功率谱密度，/>

为平均自功率谱密度，d(n)为数据窗，j为虚数单位，ω为角频率，k为数据段中的数据序号。

7.根据权利要求6所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，补零后的数据段的平均互功率谱密度计算方法为：

采用公式：

得到平均互功率谱密度；其中，Pⁱ _XY(ω)为互功率谱密度，

为平均互功率谱密度，X^*(k)为X(k)的转置，Y(k)为第二地震计经傅里叶变换后的补零数据段。

8.根据权利要求7所述的一种利用改进welch法计算地震计自噪声的方法，其特征在于，根据各自地震计的传递函数、平均自功率谱密度和平均互功率谱密度得到各自地震计的自噪声，包括：

采用公式：

计算各自地震计的自噪声；其中，P₁为第一地震计的自噪声，P₂为第二地震计的自噪声，H₁为第一地震计的传递函数，H₂为第二地震计的传递函数。

Images