基于Welch方法凝聚譜的偏差與方差估計(jì)

孫豐霖

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

1 引言

在譜分析中，凝聚譜(MSC)是衡量?jī)蓚€(gè)序列在頻域上相關(guān)性的重要工具，有十分廣泛的應(yīng)用[1-4]。在實(shí)際應(yīng)用中，常常采用Welch方法對(duì)MSC進(jìn)行估計(jì)[1-4]，通過(guò)加窗和重疊分段兩種方式來(lái)降低估計(jì)的偏差與方差，從而提高M(jìn)SC估計(jì)量的準(zhǔn)確度[1-5]。

估計(jì)量的偏差和方差是衡量一種估計(jì)方法優(yōu)劣的兩個(gè)重要方面，因此了解Welch方法得到的MSC估計(jì)量的偏差與方差是十分必要的。非重疊情形下，MSC估計(jì)量的偏差與方差可以直接通過(guò)分布函數(shù)獲得[3]，然而，重疊情形下MSC分布尚未導(dǎo)出，致使MSC偏差與方差無(wú)法計(jì)算。

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

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

2 Welch方法

假設(shè)x(t)和y(t)為兩個(gè)總長(zhǎng)度為N的離散時(shí)間序列，分別進(jì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為總分段數(shù)，L為段長(zhǎng)，D分段間的延遲量，則重疊率為p=100(1-D/L)%，通過(guò)式(2)估計(jì)MSC。

(2)

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

(3)

(4)

其中Nd是非重疊情形下的分段數(shù)

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

(5)

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

(6)

(7)

(8)

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

(9)

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

3 隨機(jī)模擬

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

表1 隨機(jī)模擬參數(shù)

(10)

(11)

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

圖1 重疊率和分段數(shù)對(duì)MSC估計(jì)量偏差和方差的影響

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

因此，在實(shí)際應(yīng)用中，50%成為估計(jì)MSC常用的重疊率，原因包括，一是在有限的數(shù)據(jù)量下能夠讓MSC估計(jì)量的偏差和方差有效降低，二是計(jì)算簡(jiǎn)單，MSC估計(jì)量的統(tǒng)計(jì)分布變化小，可利用非重疊情形下的結(jié)論進(jìn)行擴(kuò)展。

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

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

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

圖2 偏差和方差的估計(jì)效果

從圖2可見(jiàn)，重疊情形下，由于Bd和Vd與真實(shí)值之間差距過(guò)大，因此使用Nd來(lái)估計(jì)偏差與方差是不合適的，這個(gè)結(jié)論在其它Ns同樣成立，因此在接下來(lái)評(píng)價(jià)估計(jì)效果時(shí)，不再對(duì)Nd進(jìn)行討論。

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

(12)

(13)

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

表2 偏差估計(jì)值與實(shí)際值之間相對(duì)誤差ζ/%

表3 方差估計(jì)值與實(shí)際值之間相對(duì)誤差η/%

從表2可見(jiàn)，在25%重疊率下，Ns、Nc1和Nc2三種修正方法得到的偏差估計(jì)值均能使相對(duì)誤差控制在5%以?xún)?nèi)，在50%重疊率下，從Ns>9開(kāi)始，ζs出現(xiàn)一定程度的增長(zhǎng)，超過(guò)了5%的閾值范圍，且隨Ns增大沒(méi)有表現(xiàn)出明顯規(guī)律，說(shuō)明50%重疊率下使用Ns估計(jì)偏差的效果不夠穩(wěn)定。ζc1和ζc2依然能夠控制在5%以?xún)?nèi)。在75%重疊率下，ζs在所有Ns下均超過(guò)了34%，且隨Ns增大呈現(xiàn)一定上升趨勢(shì)。在Ns<11時(shí)的ζc1超過(guò)了5%，在Ns≥11時(shí)的ζc1小于5%，而Nc2始終能夠控制ζc2在5%以?xún)?nèi)。

從表3可見(jiàn)，在25%重疊率下，Ns、Nc1和Nc2三種修正方法得到的方差估計(jì)值均能使相對(duì)誤差控制在0.6%以?xún)?nèi)，在50%重疊率下，從Ns≥15開(kāi)始，ηs出現(xiàn)一定程度的增長(zhǎng)，超過(guò)了5%的閾值范圍，且隨Ns增大呈現(xiàn)一定程度的上升，說(shuō)明50%重疊率下使用Ns計(jì)算方差效果不夠穩(wěn)定。ηc1和ηc2能夠控制在1.2%以?xún)?nèi)。在75%重疊率下，ηs在所有Ns下均超過(guò)了15%，且隨Ns增大呈現(xiàn)一定上升趨勢(shì)。在Ns<9時(shí)的ηc1超過(guò)了5%，在Ns≥9時(shí)ηc1控制在5%以?xún)?nèi)，且隨著Ns增大呈現(xiàn)一定程度的下降，Nc2在Ns≥5能夠?qū)ⅵ莄2控制在5%以?xún)?nèi)，也隨著Ns增大逐漸減小。

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

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

4 穩(wěn)定性分析

在第3節(jié)中，生成給定γ2的兩個(gè)時(shí)間序列是基于服從N(0，1)白噪聲序列(以下稱(chēng)為“基序列”)，下面將基序列分布進(jìn)行修改，首先提高方差，變?yōu)榉腘(0，9)的白噪聲序列，其次，弱化不同時(shí)間點(diǎn)間相互獨(dú)立的假設(shè)，變?yōu)閰?shù)為0.5的AR(1)序列，重復(fù)偏差和方差實(shí)際值的計(jì)算過(guò)程，驗(yàn)證本文結(jié)論的穩(wěn)定性。

當(dāng)基序列變?yōu)镹(0，9)的白噪聲序列(或AR(1))時(shí)，計(jì)算此時(shí)的實(shí)際偏差與方差和N(0，1) 情況下實(shí)際偏差與方差的相對(duì)誤差

(14)

(15)

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

圖3 穩(wěn)定性分析

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

5 結(jié)論

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

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

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

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