基于Welch方法凝聚譜的偏差與方差估計

孫豐霖

(中國海洋大學海洋與大氣學院，山東青島 266100)

1 引言

在譜分析中，凝聚譜(MSC)是衡量兩個序列在頻域上相關性的重要工具，有十分廣泛的應用[1-4]。在實際應用中，常常采用Welch方法對MSC進行估計[1-4]，通過加窗和重疊分段兩種方式來降低估計的偏差與方差，從而提高MSC估計量的準確度[1-5]。

估計量的偏差和方差是衡量一種估計方法優劣的兩個重要方面，因此了解Welch方法得到的MSC估計量的偏差與方差是十分必要的。非重疊情形下，MSC估計量的偏差與方差可以直接通過分布函數獲得[3]，然而，重疊情形下MSC分布尚未導出，致使MSC偏差與方差無法計算。

一些國外學者針對重疊情形下的MSC估計量的分布進行了研究。例如Bortel和Sovka[4]給出了75%重疊率下MSC估計量的漸進分布，該方法以非重疊情形下的MSC估計量的分布函數為基礎，對其中的自由度參數進行修正，并給出適用范圍，然而文章并沒有涉及實際中常用的50%重疊率。Gallet和Julien[2]針對50%和75%重疊率下，對MSC的顯著性閾值進行討論，針對MSC總體值γ2為0的情形，對MSC估計量的分布進行討論，從而得出閾值的適用范圍，不過文章沒有涉及γ2不為0的情形。

國內學者對Welch方法估計量的偏差與方差的討論大多集中在功率譜方面，包括窗函數、重疊率和分段數等因素的影響[5-10]，而對通過Welch方法估計MSC研究較少，因此本文基于這一現狀，對Welch方法下MSC估計量的偏差與方差進行研究，給出四種估計偏差和方差的方法，并通過隨機模擬，基于估計值和實際值之間的相對誤差，對四種方法進行比較并給出各個方法的適用范圍，最后對結論的穩定性進行分析。

2 Welch方法

假設x(t)和y(t)為兩個總長度為N的離散時間序列，分別進行可重疊的分段，即

xi(t)=x((i-1)D+t)

yi(t)=y((i-1)D+t)

(1)

其中i=1，2，…，Ns，t=1，2，…，L，Ns為總分段數，L為段長，D分段間的延遲量，則重疊率為p=100(1-D/L)%，通過式(2)估計MSC。

(2)

其中F為Fourier變換，*為復共軛，w(t)為數據窗。在非重疊情形下，MSC的分布已經由Carter 等[3]證得，此時MSC的偏差(記為Bias)與方差(記為Var)如下：

(3)

(4)

其中Nd是非重疊情形下的分段數

Nd=(Ns-1)(1-p)+1

(5)

若將式(3)和式(4)中的Nd替換為Ns來估計MSC的偏差與方差，這實際上默認了重疊并不會對MSC的分布產生影響，但事實上，重疊的做法讓部分數據多次使用，讓相鄰分段相關性增強。Gallet和Julien[2]針對MSC顯著性檢驗問題，給出了兩種修正Ns的方法。第一種方法(式(6))基于修正功率譜自由度的思想，即將窗函數和重疊率的影響考慮進來。

(6)

(7)

(8)

第二種方法是基于參數估計的思想，當γ2為0時，在非重疊情形下有E2=1/Nd，因此，若假定重疊前后MSC估計量分布未發生較大變化，則可以用Nc2(式(9))作為Nd的修正。

(9)

然而，這兩種修正方法僅討論當γ2為0時，修正后的MSC分布函數是否與隨機模擬得到的累積經驗分布函數(基于兩個獨立白噪聲序列計算得到)的一致。在本文，探究這兩種修正方法是否能應用于γ2≥0時的偏差與方差估計。

3 隨機模擬

采用隨機模擬生成不同MSC真實值γ2、分段數Ns和重疊率p下的MSC樣本，隨機模擬的參數見表1，模擬序列生成方法參考Carter等[3]，采用Hanning數據窗[2]，重復m=1×105次實驗，得到MSC估計量的偏差B(式(10))與方差V(式(11))的實際值，作為真實值的代表。

表1 隨機模擬參數

(10)

(11)

首先，探究分段數和重疊率對MSC估計量實際偏差與方差的影響。從圖1a可見，偏差會隨著γ2的增大而逐漸減小直至0(當γ2=1時)。當分段數Ns=5時，25%和50%重疊率下的偏差曲線比較接近，而75%重疊率下的偏差明顯增大，當提高分段數至15和30時，偏差出現明顯的減小。從圖1b可見，方差隨著γ2的增大出現先上升后下降的趨勢，并在γ2=1/3時達到最大值。當分段數Ns=5時，25%和50%重疊率下的實際方差曲線比較接近，而75%重疊率下的方差明顯增大，當提高分段數至15和30時，方差曲線出現明顯的降低，這與非重疊情形下的結論相吻合[3]。因此，重疊情形下，隨著分段數的增加，偏差和方差明顯降低。

圖1 重疊率和分段數對MSC估計量偏差和方差的影響

從圖1可見，在給定分段數Ns=5情形下，25%和50%重疊率的偏差和方差曲線十分接近，不過25%重疊率所使用的數據是50%情形下的1.33倍，數據量提升了約33%，但方差和偏差卻沒有得到相應程度的減小。雖然75%重疊率進一步利用了有限的數據，但這導致了MSC的分布出現較大變化，使得MSC估計量的統計性質由非重疊下情形向75%重疊率推廣的過程中產生一定問題，這已經在75%重疊率閾值的適用范圍的結論中得到了驗證[2]。

因此，在實際應用中，50%成為估計MSC常用的重疊率，原因包括，一是在有限的數據量下能夠讓MSC估計量的偏差和方差有效降低，二是計算簡單，MSC估計量的統計分布變化小，可利用非重疊情形下的結論進行擴展。

下面，在不同重疊率下，比較四種計算偏差和方差的方法的估計效果。將Nd、Ns、Nc1和Nc2分別代入式(3)，計算偏差估計值，分別記為Bd、Bs、Bc1和Bc2，代入式(4)，計算方差估計值，分別記為Vd、Vs、Vc1和Vc2，以Ns=5情形為例，繪制偏差和方差估計值曲線(圖2)。

對偏差而言，在25%重疊率下，除Bd外，其它三條估計曲線與實際值曲線B十分接近，但隨著重疊率提高到50%，Bs出現一定程度的偏離，Bc1和Bc2接近實際值曲線B，當重疊率上升到75%，Bs的偏離程度更加明顯，同時Bc1也位于實際值曲線B下方，只有Bc2仍然接近實際值曲線。

對方差而言，結論與偏差相似。在25%重疊率下，Vs、Vc1和Vc2與實際值曲線非常吻合，50%的重疊率令Vs開始出現細微的偏離，Vc1和Vc2與實際值吻合度較好，然而75%重疊率下，Vd高估了方差，而Vs和Vc1估計值偏低，Vc2也在γ2<1/3時高估了實際方差，不過這種影響在γ2>1/3后逐漸減弱消失。

圖2 偏差和方差的估計效果

從圖2可見，重疊情形下，由于Bd和Vd與真實值之間差距過大，因此使用Nd來估計偏差與方差是不合適的，這個結論在其它Ns同樣成立，因此在接下來評價估計效果時，不再對Nd進行討論。

下面，針對不同的分段數，對Ns、Nc1和Nc2估計偏差和方差的整體效果進行討論。對方法Ns，定義偏差和方差的相對誤差為

(12)

(13)

其中Φ為隨機模擬中γ2的取值集合，mean為取平均，同理，Nc1和Nc2偏差和方差的相對誤差記為ζc1，ζc2，ηc1和ηc2。ζ和η能衡量在整個γ2的取值范圍內，估計方法得到的偏差和方差估計值與實際值之間的整體差距。三種估計方法的相對誤差結果見表2，選擇5%作為判斷該方法是否適用的閾值，即如果某估計方法在重疊率p和分段數Ns下的相對誤差小于5%，則認為該方法可以用來估計偏差和方差。選擇5%作為閾值的理由在于，Gallet和Julien[2]利用KS檢驗對γ2=0情形下對MSC估計量的分布進行了討論，并利用式(6)和式(9)對非重疊情形下的分布函數進行修正，給出50%和75%重疊率下修正后分布函數的適用范圍。在這些適用范圍中，當γ2=0時，本文計算的偏差與方差的相對誤差均處于5%以內，因此將5%作為衡量估計值與實際值之間是否吻合的閾值。

表2 偏差估計值與實際值之間相對誤差ζ/%

表3 方差估計值與實際值之間相對誤差η/%

從表2可見，在25%重疊率下，Ns、Nc1和Nc2三種修正方法得到的偏差估計值均能使相對誤差控制在5%以內，在50%重疊率下，從Ns>9開始，ζs出現一定程度的增長，超過了5%的閾值范圍，且隨Ns增大沒有表現出明顯規律，說明50%重疊率下使用Ns估計偏差的效果不夠穩定。ζc1和ζc2依然能夠控制在5%以內。在75%重疊率下，ζs在所有Ns下均超過了34%，且隨Ns增大呈現一定上升趨勢。在Ns<11時的ζc1超過了5%，在Ns≥11時的ζc1小于5%，而Nc2始終能夠控制ζc2在5%以內。

從表3可見，在25%重疊率下，Ns、Nc1和Nc2三種修正方法得到的方差估計值均能使相對誤差控制在0.6%以內，在50%重疊率下，從Ns≥15開始，ηs出現一定程度的增長，超過了5%的閾值范圍，且隨Ns增大呈現一定程度的上升，說明50%重疊率下使用Ns計算方差效果不夠穩定。ηc1和ηc2能夠控制在1.2%以內。在75%重疊率下，ηs在所有Ns下均超過了15%，且隨Ns增大呈現一定上升趨勢。在Ns<9時的ηc1超過了5%，在Ns≥9時ηc1控制在5%以內，且隨著Ns增大呈現一定程度的下降，Nc2在Ns≥5能夠將ηc2控制在5%以內，也隨著Ns增大逐漸減小。

表4 不同重疊率下各方法適用范圍

結合表2和表3結果，給出不同重疊率和分段數下，給出三種方法估計MSC的適用范圍(表4)。在25%重疊率下，三種方法計算的偏差和方差均與實際值比較吻合，其中直接使用Ns計算是比較準確和快速的方法。在50%重疊率下，Nc1和Nc2計算偏差和方差相對準確，Nc1只需根據重疊率和數據窗進行計算，而Nc2需通過隨機模擬生成MSC樣本估計，因此前者是50%重疊率時計算偏差與方差比較好的選擇。在75%重疊率下，Nc1的適用范圍小于Nc2的適用范圍，即使在Ns≥11時，使用Nc1計算的偏差與方差的相對誤差高于Nc2，不過差距隨著Ns增大而逐漸減小。

4 穩定性分析

在第3節中，生成給定γ2的兩個時間序列是基于服從N(0，1)白噪聲序列(以下稱為“基序列”)，下面將基序列分布進行修改，首先提高方差，變為服從N(0，9)的白噪聲序列，其次，弱化不同時間點間相互獨立的假設，變為參數為0.5的AR(1)序列，重復偏差和方差實際值的計算過程，驗證本文結論的穩定性。

當基序列變為N(0，9)的白噪聲序列(或AR(1))時，計算此時的實際偏差與方差和N(0，1) 情況下實際偏差與方差的相對誤差

(14)

(15)

其中，i=1代表N(0，9)，i=2代表AR(1)，繪制ζi和ηi的分布直方圖(圖3)。當模擬基序列改為N(0，9)時，MSC偏差的實際值變化幅度處于[-13%，13%]，約92%的相對誤差集中在[-5%，5%]，而方差實際值變化幅度處于[-2.75%，2.75%]，約93%的相對誤差集中在[-1.25%，1.25%]。當模擬基序列改為AR(1)時，偏差與方差實際值的相對誤差分布狀況與N(0，9)情形下的結果十分相似，約91%的偏差相對誤差集中在[-5%，5%]，約94%的方差相對誤差集中在[-1.25%，1.25%]，使用KS test檢驗圖3a、圖3b中樣本是否來自同一分布，p值分別為0.86和0.98，不拒絕來自同一分布的假設。

圖3 穩定性分析

綜上，隨機模擬中基序列改變并不會對偏差和方差的實際值產生較大影響，同時，Ns、Nc1和Nc2三種方法的適用范圍(表4)并未隨著基序列改變發生變化，因此，本文給出的Ns、Nc1和Nc2三種方法的適用范圍具有一定的穩定性。

5 結論

本文研究了基于Welch方法，不同重疊率下MSC估計量偏差和方差的估計問題，給出四種估計偏差和方差方法，并通過隨機模擬對這四種方法的估計效果進行比較，給出各個方法的適用范圍，并通過改變隨機模擬基序列的方式驗證本文結論的穩定性。主要結論包括：

1)在25%重疊率下，使用Ns、Nc1和Nc2三種方法得到的偏差與方差估計值與實際值比較吻合。由于Ns不需通過計算即可得到，是一種比較合適的估計方法。

2)在50%重疊率下，使用Nc1和Nc2得到的偏差與方差與實際值均比較吻合。前者計算需要窗函數和重疊率，后者需要隨機模擬樣本對參數進行估計，同時需要足夠樣本量保證參數的穩定性，因此前者是相對比較合適的方法。

3)在75%重疊率下，Nc1的適用范圍需要保證Ns≥11，而Nc2的適用范圍需要保證Ns≥5，且后者的相對誤差要優于前者，不過這種差距隨著Ns的增大逐步減小。