臺(tái)基噪聲功率譜估計(jì)中的Welch參數(shù)選擇

尹康達(dá) 李小軍 張曉剛 丁 成

1）中國(guó)邢臺(tái) 054000 河北紅山巨厚沉積與地震災(zāi)害國(guó)家野外科學(xué)觀測(cè)研究站

2）中國(guó)石家莊 050021 河北省地震局

0 引言

地震臺(tái)站噪聲水平是影響觀測(cè)數(shù)據(jù)質(zhì)量的重要因素（林彬華等，2020），功率譜密度（power spectral density，簡(jiǎn)稱PSD）是研究臺(tái)基噪聲并對(duì)其進(jìn)行量化的重要手段（McNamara et al，2009；葛洪魁等，2013；楊千里等，2019）。噪聲的頻段范圍分布較廣，長(zhǎng)周期的自然界噪聲和高頻段的人類活動(dòng)是地震觀測(cè)中的主要噪聲源（葛洪魁等，2013；謝江濤等，2021）；可靠的臺(tái)基噪聲功率譜估計(jì)方法可為地震觀測(cè)環(huán)境等級(jí)劃分、觀測(cè)數(shù)據(jù)可用性評(píng)估、地震監(jiān)測(cè)能力評(píng)估、噪聲源分析提取等提供依據(jù)。

地震預(yù)警系統(tǒng)是基于P 波初始記錄進(jìn)行計(jì)算的。預(yù)警臺(tái)站數(shù)量的迅速增加對(duì)觀測(cè)數(shù)據(jù)質(zhì)量提出了更高的要求（Enferadi et al，2021）。臺(tái)基噪聲功率譜估計(jì)及RMS 計(jì)算已成為地震臺(tái)網(wǎng)的日常工作。Welch 修正周期圖法（Welch，1967）是計(jì)算噪聲PSD 的常用方法，但對(duì)其參數(shù)的選擇目前尚無(wú)統(tǒng)一標(biāo)準(zhǔn)，這導(dǎo)致計(jì)算結(jié)果存在較大偏差（Li et al，2015；Xu et al，2019）。Li 等（2015）對(duì)Welch 方法中幾個(gè)常用參數(shù)，如窗口函數(shù)、窗口長(zhǎng)度、分段數(shù)據(jù)長(zhǎng)度、重疊率等的選取進(jìn)行了研究，發(fā)現(xiàn)不同參數(shù)組合對(duì)PSD 的計(jì)算結(jié)果會(huì)產(chǎn)生較大影響。本研究在此基礎(chǔ)上，進(jìn)一步開(kāi)展工作，以期對(duì)規(guī)范日常臺(tái)基噪聲計(jì)算中的參數(shù)選擇提供參考。

1 方法

1.1 功率譜估計(jì)方法

通常用PSD 來(lái)描述噪聲的頻域特性，評(píng)價(jià)一個(gè)地震臺(tái)站的觀測(cè)等級(jí)。利用Welch 方法進(jìn)行傅里葉變換的數(shù)據(jù)較短，在許多情況下Welch 方法比其他方法計(jì)算量更少（Welch，1967），是目前進(jìn)行功率譜估計(jì)的通用方法（Xu et al，2019）。該方法假設(shè)原始數(shù)據(jù)為一個(gè)二階穩(wěn)態(tài)隨機(jī)序列，并將其分成若干長(zhǎng)度相等的數(shù)據(jù)段，對(duì)分段數(shù)據(jù)進(jìn)行加窗重疊，并通過(guò)修正周期圖來(lái)進(jìn)行譜的平滑和方差的減小，從而在一定程度上減小譜估計(jì)的誤差。通常窗口長(zhǎng)度與分段數(shù)據(jù)長(zhǎng)度相等，文中統(tǒng)一稱為窗口長(zhǎng)度。Welch方法數(shù)據(jù)分段示意圖如圖1 所示。

圖1 Welch 方法數(shù)據(jù)分段示意圖Ld 為原始數(shù)據(jù)長(zhǎng)度；S1，S2，…，Sn 為分段數(shù)據(jù)；窗口長(zhǎng)度為L(zhǎng)c；重疊部分長(zhǎng)度為L(zhǎng)o；重疊率Ro=Lo/LcFig.1 Record segmentation of Welch method

Welch（1967）在進(jìn)行功率譜估計(jì)方法研究時(shí)，將方差作為一個(gè)重要衡量指標(biāo)，從理論角度推導(dǎo)出了通過(guò)分段和重疊可在一定程度上減小方差，并從計(jì)算耗時(shí)角度指出，當(dāng)窗口長(zhǎng)度小于原始數(shù)據(jù)長(zhǎng)度的平方根時(shí)，耗時(shí)更短。本文對(duì)重疊率的選擇進(jìn)行了研究，研究結(jié)果可對(duì)該方法進(jìn)行一定的補(bǔ)充。

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

窗口長(zhǎng)度和重疊率取值組合不同，功率譜估計(jì)結(jié)果亦不同。本文從分段數(shù)據(jù)匹配的角度，研究窗口長(zhǎng)度與重疊率間的關(guān)系。

窗口長(zhǎng)度占比（本文簡(jiǎn)稱為窗口長(zhǎng)度）Rcd=Lc/Ld，當(dāng)Rcd＞50%且重疊率較小時(shí)，第2 段數(shù)據(jù)末端會(huì)超出原始數(shù)據(jù)，從而導(dǎo)致該數(shù)據(jù)段不參與計(jì)算。當(dāng)?shù)? 段數(shù)據(jù)末端與原始數(shù)據(jù)重合時(shí)（圖2），將對(duì)應(yīng)的重疊率定義為最低重疊率Rom，計(jì)算公式為

圖2 最低重疊率示意圖Ld 為原始數(shù)據(jù)長(zhǎng)度；S1，S2，…，Sn 為分段數(shù)據(jù)；窗口長(zhǎng)度為L(zhǎng)c；重疊部分長(zhǎng)度為L(zhǎng)o；重疊率Ro=Lo/LcFig.2 Illustration of minimum overlap rate

當(dāng)窗口長(zhǎng)度固定時(shí)，僅當(dāng)Sn末端與原始數(shù)據(jù)重合時(shí)，各段疊加之后的數(shù)據(jù)才是完整的（圖1），將此時(shí)的重疊率定義為匹配重疊率Rm；當(dāng)Ro≠Rm時(shí)（圖3），Ln+1末端會(huì)超出范圍，無(wú)法參與計(jì)算，使得疊加后無(wú)法完整覆蓋原始數(shù)據(jù)，造成長(zhǎng)度為L(zhǎng)s的數(shù)據(jù)缺失。

圖3 非匹配重疊示意圖Ld 為原始數(shù)據(jù)長(zhǎng)度；S1，S2，…，Sn 為分段數(shù)據(jù)；窗口長(zhǎng)度為L(zhǎng)c；重疊部分長(zhǎng)度為L(zhǎng)o；重疊率Ro=Lo/LcFig.3 Illustration of unmatching overlap

匹配重疊率計(jì)算過(guò)程如下：

能夠參與計(jì)算的有效分段數(shù)為

當(dāng)n≥2 且為整數(shù)時(shí)，對(duì)應(yīng)的重疊率Ro即為匹配重疊率Rm，公式為

當(dāng)n=2 時(shí)，即最低重疊率；當(dāng)n=2，3，…，1 000 時(shí)，匹配重疊率隨窗口長(zhǎng)度的變化曲線如圖4 所示。當(dāng)選擇其下方區(qū)域的參數(shù)組合時(shí)，僅第1 段數(shù)據(jù)參與計(jì)算，造成數(shù)據(jù)缺失。最低重疊率曲線可作為對(duì)Li 等（2015）提出的K 線的解釋。

圖4 匹配重疊率隨窗口長(zhǎng)度的變化紅色曲線即最低重疊率Fig.4 Curves of matching overlap rate with window length

1.3 平滑度評(píng)價(jià)指標(biāo)

通過(guò)周期圖法計(jì)算PSD 時(shí)，對(duì)有限長(zhǎng)數(shù)據(jù)的周期性拓展會(huì)產(chǎn)生較大的隨機(jī)偏差，功率譜密度曲線波動(dòng)較大，平滑度較低，Welch 方法可在一定程度上改善功率譜密度曲線的平滑度。參考Welch（1967）的評(píng)價(jià)指標(biāo)，通過(guò)方差來(lái)表征譜曲線的平滑度，方差計(jì)算公式如下

其中，var 為PSD 曲線方差；D為原始數(shù)據(jù)采樣點(diǎn)數(shù)；pi為各頻點(diǎn)下的PSD 值；pa為各頻點(diǎn)PSD 均值。

2 數(shù)據(jù)處理和分析

寬頻帶地震計(jì)具有較小的自噪聲和較大的動(dòng)態(tài)范圍（Sleeman et al，2006；李小軍等，2019），是目前地震臺(tái)網(wǎng)常用觀測(cè)設(shè)備之一。選取CTS-1 型甚寬頻帶地震計(jì)進(jìn)行觀測(cè)實(shí)驗(yàn)，由于甚寬頻帶地震計(jì)對(duì)溫度變化較敏感（田文德等，2013），傳感器安裝在山洞內(nèi)。選取1 h觀測(cè)數(shù)據(jù)（采樣率為100Hz）進(jìn)行功率譜估計(jì)實(shí)驗(yàn)。通過(guò)改變窗口長(zhǎng)度和重疊率，來(lái)研究其對(duì)PSD 平滑度的影響，從而確定窗口長(zhǎng)度和重疊率的合理取值區(qū)間。

數(shù)據(jù)處理基于MATLAB 軟件平臺(tái)的Pwelch 函數(shù)對(duì)記錄的垂直分量進(jìn)行計(jì)算。由于線性趨勢(shì)會(huì)抵消低頻段的頻譜能量，使譜估計(jì)產(chǎn)生較大失真（Mcnamara et al，2004），因此，在進(jìn)行功率譜計(jì)算之前，對(duì)所有記錄數(shù)據(jù)進(jìn)行去均值和線性趨勢(shì)處理。

通過(guò)方差對(duì)功率譜密度曲線的平滑度進(jìn)行定量分析。窗口長(zhǎng)度為1%—100%，步長(zhǎng)1%；重疊率為0—99%，步長(zhǎng)1%，對(duì)10 000 種窗口長(zhǎng)度和重疊率的組合，進(jìn)行功率譜估計(jì)實(shí)驗(yàn)。

按式（4）計(jì)算10 000 條速度功率譜密度曲線的方差，研究不同參數(shù)組合對(duì)平滑度的影響，計(jì)算結(jié)果如圖5 所示。由圖5 可見(jiàn)，方差隨重疊率和窗口長(zhǎng)度呈階躍式變化，突變的位置可以形成一系列曲線與圖4匹配重疊率曲線相對(duì)應(yīng)。

圖5 速度PSD 方差隨Welch 參數(shù)的變化Fig.5 Variance of velocity PSD at different Welch parameters

隨著窗口長(zhǎng)度的增加，方差呈趨勢(shì)性增大，這主要是因?yàn)槲椿ハ嘀丿B的分段數(shù)減少，減弱了利用分段平均進(jìn)行平滑的效果，雖然通過(guò)較大的重疊率，可使分段數(shù)增加，但重疊部分的平均對(duì)平滑度的改善效果不顯著。

當(dāng)窗口長(zhǎng)度大于50%時(shí)，在最低重疊率下方區(qū)域方差整體偏大。由于此時(shí)沒(méi)有通過(guò)取均值進(jìn)行平滑的過(guò)程，且與原始數(shù)據(jù)偏差較大，此區(qū)域內(nèi)方差不隨重疊率而改變，為不可靠區(qū)域。

當(dāng)固定窗口長(zhǎng)度時(shí)，在2 個(gè)匹配重疊率間隔內(nèi)，隨著重疊率的增大，方差逐漸增大，這主要是因?yàn)榇翱诏B加之后不能完全覆蓋原始數(shù)據(jù)。對(duì)于可近似為隨機(jī)信號(hào)的臺(tái)基噪聲，隨著缺失數(shù)據(jù)的增多，方差變大，譜估計(jì)的誤差增大；當(dāng)取匹配重疊率時(shí)，方差突然減小，因?yàn)榇藭r(shí)窗口經(jīng)過(guò)疊加對(duì)原始數(shù)據(jù)形成了全覆蓋，功率譜估計(jì)相對(duì)準(zhǔn)確，這也驗(yàn)證了選取匹配重疊率的重要性。同時(shí)，對(duì)于同一窗口長(zhǎng)度取不同的匹配重疊率時(shí)，方差變化不大，因?yàn)橹丿B的主要作用是補(bǔ)償由加窗所造成的信號(hào)兩端衰減。在同一窗口長(zhǎng)度下，增大重疊率，雖然在保證分辨率不降低的前提下，增加了分段數(shù)，但多為重復(fù)數(shù)據(jù)，取平均對(duì)功率譜密度曲線平滑度的改善效果不明顯。

3 結(jié)論與討論

針對(duì)臺(tái)基噪聲功率譜估計(jì)中的參數(shù)選擇問(wèn)題，通過(guò)對(duì)Welch 方法不同參數(shù)組合下功率譜密度曲線的平滑度進(jìn)行分析，并以方差指標(biāo)來(lái)進(jìn)行定量研究，可為地震臺(tái)站觀測(cè)環(huán)境綜合評(píng)價(jià)和觀測(cè)數(shù)據(jù)質(zhì)量分析提供參考。

在本研究中定義了最低重疊率，當(dāng)窗口長(zhǎng)度大于原始數(shù)據(jù)的50%，重疊率取值小于最低重疊率時(shí)，功率譜估計(jì)結(jié)果不可靠；定義了匹配重疊率，特定窗口長(zhǎng)度對(duì)應(yīng)一系列匹配重疊率，且在匹配重疊率下的臺(tái)基噪聲功率譜估計(jì)結(jié)果更準(zhǔn)確。

Welch 方法將原始數(shù)據(jù)進(jìn)行分段時(shí)，分段長(zhǎng)度和重疊率取值都是單一的，這導(dǎo)致重疊率的選取受到限制。在2 個(gè)匹配重疊率之間，隨著重疊率的增大，缺失數(shù)據(jù)也會(huì)增加，而且當(dāng)分段數(shù)較少時(shí)，2 個(gè)匹配重疊率的間隔較大，此時(shí)若想通過(guò)增加重疊率對(duì)功率譜估計(jì)結(jié)果有所改善，必須進(jìn)行跨越式增大，這將會(huì)使計(jì)算量大大增加。

因此，如何根據(jù)研究目的和原始數(shù)據(jù)特征，選擇更合理的分段、重疊方式，并通過(guò)分段加權(quán)進(jìn)行平均，需要開(kāi)展進(jìn)一步研究。