0 引言

地震台站噪声水平是影响观测数据质量的重要因素（林彬华等，2020），功率谱密度（power spectral density，简称PSD）是研究台基噪声并对其进行量化的重要手段（McNamara et al，2009；葛洪魁等，2013；杨千里等，2019）。噪声的频段范围分布较广，长周期的自然界噪声和高频段的人类活动是地震观测中的主要噪声源（葛洪魁等，2013；谢江涛等，2021）；可靠的台基噪声功率谱估计方法可为地震观测环境等级划分、观测数据可用性评估、地震监测能力评估、噪声源分析提取等提供依据。

地震预警系统是基于P波初始记录进行计算的。预警台站数量的迅速增加对观测数据质量提出了更高的要求（Enferadi et al，2021）。台基噪声功率谱估计及RMS计算已成为地震台网的日常工作。Welch修正周期图法（Welch，1967）是计算噪声PSD的常用方法，但对其参数的选择目前尚无统一标准，这导致计算结果存在较大偏差（Li et al，2015；Xu et al，2019）。Li等（2015）对Welch方法中几个常用参数，如窗口函数、窗口长度、分段数据长度、重叠率等的选取进行了研究，发现不同参数组合对PSD的计算结果会产生较大影响。本研究在此基础上，进一步开展工作，以期对规范日常台基噪声计算中的参数选择提供参考。

1 方法 1.1 功率谱估计方法

通常用PSD来描述噪声的频域特性，评价一个地震台站的观测等级。利用Welch方法进行傅里叶变换的数据较短，在许多情况下Welch方法比其他方法计算量更少（Welch，1967），是目前进行功率谱估计的通用方法（Xu et al，2019）。该方法假设原始数据为一个二阶稳态随机序列，并将其分成若干长度相等的数据段，对分段数据进行加窗重叠，并通过修正周期图来进行谱的平滑和方差的减小，从而在一定程度上减小谱估计的误差。通常窗口长度与分段数据长度相等，文中统一称为窗口长度。Welch方法数据分段示意图如图 1所示。

Welch（1967）在进行功率谱估计方法研究时，将方差作为一个重要衡量指标，从理论角度推导出了通过分段和重叠可在一定程度上减小方差，并从计算耗时角度指出，当窗口长度小于原始数据长度的平方根时，耗时更短。本文对重叠率的选择进行了研究，研究结果可对该方法进行一定的补充。

1.2 最低重叠率、匹配重叠率

窗口长度和重叠率取值组合不同，功率谱估计结果亦不同。本文从分段数据匹配的角度，研究窗口长度与重叠率间的关系。

窗口长度占比（本文简称为窗口长度）R_cd = L_c/L_d，当R_cd＞50%且重叠率较小时，第2段数据末端会超出原始数据，从而导致该数据段不参与计算。当第2段数据末端与原始数据重合时（图 2），将对应的重叠率定义为最低重叠率R_om，计算公式为

$ R_{\mathrm{om}}=\frac{L_{\mathrm{c}}-\left(L_{\mathrm{d}}-L_{\mathrm{c}}\right)}{L_{\mathrm{c}}}=2-\frac{L_{\mathrm{d}}}{L_{\mathrm{c}}}=2-\frac{1}{R_{\mathrm{cd}}} $

(1)

当窗口长度固定时，仅当S_n末端与原始数据重合时，各段叠加之后的数据才是完整的（图 1），将此时的重叠率定义为匹配重叠率R_m；当R_o≠R_m时（图 3），L_n+1末端会超出范围，无法参与计算，使得叠加后无法完整覆盖原始数据，造成长度为L_s的数据缺失。

匹配重叠率计算过程如下：

能够参与计算的有效分段数为

$ n=\frac{L_{\mathrm{d}}\left(1-R_{\mathrm{cd}}\right)}{L_{\mathrm{c}}\left(1-R_{\mathrm{o}}\right)}+1=\frac{1-R_{\mathrm{cd}}}{R_{\mathrm{cd}}\left(1-R_{\mathrm{o}}\right)}+1 $

(2)

当n≥2且为整数时，对应的重叠率R_o即为匹配重叠率R_m，公式为

$ R_{\mathrm{m}}=\frac{n}{n-1}-\frac{1}{R_{\mathrm{cd}}(n-1)} \quad(n=2, 3, \ldots) $

(3)

当n = 2时，即最低重叠率；当n = 2，3，…，1 000时，匹配重叠率随窗口长度的变化曲线如图 4所示。当选择其下方区域的参数组合时，仅第1段数据参与计算，造成数据缺失。最低重叠率曲线可作为对Li等（2015）提出的K线的解释。

1.3 平滑度评价指标

通过周期图法计算PSD时，对有限长数据的周期性拓展会产生较大的随机偏差，功率谱密度曲线波动较大，平滑度较低，Welch方法可在一定程度上改善功率谱密度曲线的平滑度。参考Welch（1967）的评价指标，通过方差来表征谱曲线的平滑度，方差计算公式如下

$ \operatorname{var}=\frac{1}{D / 2-1} \sum\nolimits_{i=1}^{D / 2}\left(p_i-p_{\mathrm{a}}\right)^2 $

(4)

其中，var为PSD曲线方差；D为原始数据采样点数；p_i为各频点下的PSD值；p_a为各频点PSD均值。

2 数据处理和分析

宽频带地震计具有较小的自噪声和较大的动态范围（Sleeman et al，2006；李小军等，2019），是目前地震台网常用观测设备之一。选取CTS-1型甚宽频带地震计进行观测实验，由于甚宽频带地震计对温度变化较敏感（田文德等，2013），传感器安装在山洞内。选取1 h观测数据（采样率为100Hz）进行功率谱估计实验。通过改变窗口长度和重叠率，来研究其对PSD平滑度的影响，从而确定窗口长度和重叠率的合理取值区间。

数据处理基于MATLAB软件平台的Pwelch函数对记录的垂直分量进行计算。由于线性趋势会抵消低频段的频谱能量，使谱估计产生较大失真（Mcnamara et al，2004），因此，在进行功率谱计算之前，对所有记录数据进行去均值和线性趋势处理。

通过方差对功率谱密度曲线的平滑度进行定量分析。窗口长度为1%—100%，步长1%；重叠率为0—99%，步长1%，对10 000种窗口长度和重叠率的组合，进行功率谱估计实验。

按式（4）计算10 000条速度功率谱密度曲线的方差，研究不同参数组合对平滑度的影响，计算结果如图 5所示。由图 5可见，方差随重叠率和窗口长度呈阶跃式变化，突变的位置可以形成一系列曲线与图 4匹配重叠率曲线相对应。

随着窗口长度的增加，方差呈趋势性增大，这主要是因为未互相重叠的分段数减少，减弱了利用分段平均进行平滑的效果，虽然通过较大的重叠率，可使分段数增加，但重叠部分的平均对平滑度的改善效果不显著。

当窗口长度大于50%时，在最低重叠率下方区域方差整体偏大。由于此时没有通过取均值进行平滑的过程，且与原始数据偏差较大，此区域内方差不随重叠率而改变，为不可靠区域。

当固定窗口长度时，在2个匹配重叠率间隔内，随着重叠率的增大，方差逐渐增大，这主要是因为窗口叠加之后不能完全覆盖原始数据。对于可近似为随机信号的台基噪声，随着缺失数据的增多，方差变大，谱估计的误差增大；当取匹配重叠率时，方差突然减小，因为此时窗口经过叠加对原始数据形成了全覆盖，功率谱估计相对准确，这也验证了选取匹配重叠率的重要性。同时，对于同一窗口长度取不同的匹配重叠率时，方差变化不大，因为重叠的主要作用是补偿由加窗所造成的信号两端衰减。在同一窗口长度下，增大重叠率，虽然在保证分辨率不降低的前提下，增加了分段数，但多为重复数据，取平均对功率谱密度曲线平滑度的改善效果不明显。

3 结论与讨论

针对台基噪声功率谱估计中的参数选择问题，通过对Welch方法不同参数组合下功率谱密度曲线的平滑度进行分析，并以方差指标来进行定量研究，可为地震台站观测环境综合评价和观测数据质量分析提供参考。

在本研究中定义了最低重叠率，当窗口长度大于原始数据的50%，重叠率取值小于最低重叠率时，功率谱估计结果不可靠；定义了匹配重叠率，特定窗口长度对应一系列匹配重叠率，且在匹配重叠率下的台基噪声功率谱估计结果更准确。

Welch方法将原始数据进行分段时，分段长度和重叠率取值都是单一的，这导致重叠率的选取受到限制。在2个匹配重叠率之间，随着重叠率的增大，缺失数据也会增加，而且当分段数较少时，2个匹配重叠率的间隔较大，此时若想通过增加重叠率对功率谱估计结果有所改善，必须进行跨越式增大，这将会使计算量大大增加。

因此，如何根据研究目的和原始数据特征，选择更合理的分段、重叠方式，并通过分段加权进行平均，需要开展进一步研究。