波流信號交叉譜估計方法比較*

張 瑞 類淑河, 管長龍 王 紅

(1. 中國海洋大學數學科學學院 青島 266100; 2. 中國海洋大學物理海洋實驗室 青島 266100)

海洋中波浪和海流幾乎總是同時存在。大到錢塘江的怒潮,小到近海浮游生物的“隨波逐流”都與波流相互作用有密不可分的關系(高山,2003)。理論分析、數值模擬、實驗研究和現場觀測等多種手段被廣泛用于研究波流相互作用。一般來說,波浪主要通過海面風應力、波致輻射應力和底應力等影響海流。海流則主要通過水深變化和隨時空變化的流場影響波浪。兩者相互作用的影響極為廣泛,例如,流使波浪及波浪譜變形、流浪場的相互作用、波和流影響結構物載荷、運動及被動標量運輸等(王濤等,1999)。通過波流相互作用的深入討論,可以更好地解釋海洋動力現象及其演變過程。

交叉譜(cross-spectrum)可視為互相關函數的傅里葉變換。作為一種統計工具,它可用于分析不同時間序列對應頻率分量之間的相互關系。由交叉譜得到的凝聚譜可刻畫兩個序列在頻域內的相關性。在物理海洋方面,這些統計工具盡管已被應用于海流的相互作用研究,但是用于探討波流相互作用的研究工作仍相對較少。袁耀初等(2002)基于南海東北部錨定測流站450m以淺長達77d的Long Ranger ADCP測流資料以及2000m與2300m處長達7個月左右的測流資料,通過加窗周期圖法進行交叉譜估計,發現2000m與2300m流速時間序列在周期為2個多月和1個多月的振動有很好的相關性,100m與2300m流速時間序列在周期為15.5d和2d等的振動有很好的相關性。McKone(2003)基于IAS(Intra-Americas Sea)不同地點的 NLOM模擬海流數據,采用 Welch方法和MTM 方法進行交叉譜估計,發現年周期振動均存在強烈的相關性,其中一部分地點之間在9、6、4、2.4個月周期振動也存在相關性,而采用帶通濾波器濾掉年周期信號后,4個月和2.4個月已不存在相關性,6個月的相關性程度也明顯減弱。

實際觀測中,得到的時間序列長度是有限且離散的。交叉譜分析把有限長度數據得到的交叉譜估計作為真實譜的近似。交叉譜估計方法大致可分為非參數方法和參數方法兩類。非參數估計方法包括周期圖法(Jenkinset al,1968)、Welch方法(Welch,1967;Carter,1987)、多窗譜方法(Multitaper Method,簡稱“MTM方法”)(Thomson,1982)、最小方差無失真響應方法(The Minimum Variance Distortionless Response Method,簡稱“MVDR 方法”) (Benestyet al,2005)等。參數估計方法主要針對AR模型、MA模型、ARMA模型等進行參數建模,具體的方法包括最大熵估計方法(Strand,1977)、協方差法(Kay,1988)等。參數交叉譜估計方法雖然具有一些優勢如提高了頻率分辨率等,但其統計性能很難從數學上給出解析式,相比之下非參數估計方法應用更加廣泛。一種好的譜估計方法不一定對任何場合都有良好的估計性能,應依據具體應用實例選擇恰當的方法(黃建國,1991)。采用多種方法進行交叉譜估計并對比選擇,可以避免對偽峰的誤識別,有望給出更合理的結果。

本文主要目標是尋找適合于波流信號的交叉譜估計方法。第一部分介紹三種主要的非參數交叉譜估計方法,即Welch方法、MTM方法和MVDR方法。第二部分利用數值模擬,從相關信號的檢測能力和分辨能力兩個方面對三種方法進行了比較。第三部分給出風浪實驗數據凝聚譜估計,并對三種方法的估計結果進行分析比較。第四部分給出相應的結論與討論。

1 交叉譜及其估計方法

1.1 交叉譜及相關概念

設二維零均值平穩時間序列{Xt,Yt},t= 0,±1,±2,???,則Xt與Yt的互相關函數定義為

其中E[]表示數學期望。

交叉譜一般是復值的,單位為(Xt的單位)(Yt的單位)/(頻率單位Hz)。寫成極坐標形式為

其中αXY(f)為交振幅譜(cross-amplitude spectrum),表示Xt、Yt在頻率f上的相關程度。

為排除量綱的影響,對交振幅譜標準化并引入相干譜 (coherency spectrum)。記Xt的單譜為SXX(f),Yt的單譜為SYY(f),則兩信號的相干譜可表示為

相干譜表征了Xt、Yt在頻率f上的相關程度大小,在頻域內刻畫了兩個序列的相關性。粗略地講,相干譜為頻域上的相關系數。

實際應用中,多采用相干譜的平方表示頻域內信號之間的相關性,稱為凝聚譜(magnitude-squared coherence),其表達式為:

顯然,相干譜與凝聚譜都是無量綱的,在[0,1]區間內取值。該值越接近于 1,表明二者在相應頻率分量上的相關性越強。特別地,若wXY(f)=0或Msc(f)=0,說明Xt、Yt在頻率f上是不相關的; 若對所有頻率均有wXY(f)=1或 Msc(f)=1,說明Xt和Yt在頻域上是完全相關的。

1.2 交叉譜估計

類似于單個信號的譜估計,周期圖法是估計交叉譜的基本方法,但其估計方差往往較大,且容易造成頻譜泄漏。

Welch方法是經典的交叉譜估計方法。其基本思想是: 先對數據進行分段(允許其間有重疊),然后對每一段數據進行加窗交叉譜估計,最后對各段交叉譜求平均。這種方法對周期圖法進行了改進: 一方面通過分段估計再取平均,提高了交叉譜估計的自由度,從而減小了誤差; 另一方面通過加窗減少泄露,從而降低了估計偏差。

MTM方法的基本思想是以一簇數據窗代替單一數據窗,對每一數據窗所構成的時間序列進行傅里葉變換,再加權平均。數據窗可采用正交離散扁長的球體序列(Discrete Prolate Spheroidal Sequences)。與Welch方法相比,MTM方法不需要主觀定義分辨率,而且交叉譜估計的自由度大于 2,有效地減小了估計的誤差和偏差。

MVDR方法(Capon,1969)在陣列信號中有廣泛的應用,其基本思想是保證期望方向上的陣列響應,同時減少陣列輸出中的干擾,從而改善天線陣列的分辨率。MVDR方法亦可用于交叉譜和凝聚譜估計,且比Welch方法分辨率更高(Benestyet al,2005)。特別地,可選取傅里葉矩陣這一特殊的酉陣進行凝聚譜估計。

若已估計出交叉譜和單譜,可進一步估計交振幅譜、相干譜和凝聚譜。記交叉譜估計為?XY(f),Xt,Yt的單譜估計相應為?XX(f)、?YY(f),則交振幅譜估計可表示為

相干譜估計可表示為

凝聚譜估計可表示為

以往的研究(Koopmans,1974; Bortelet al,2007;Galletet al,2011)詳細討論了凝聚譜的臨界值: 采取Welch方法進行估計時,每段有 50%的重疊且使用Hanning窗,則漸近服從自由度為的χ2分布(Percivalet al,1993)。這一結果可用于凝聚譜的非零相關性檢驗: 在顯著性水平α下,凝聚譜估計的臨界值為 1–α2/(ν–2),其中ν為自由度,超過這一臨界值,可以認為兩序列在相應頻率上顯著相關。

2 模擬數據分析

為比較三種估計方法的相關信號的檢測能力和頻率分辨能力,同時驗證前面的臨界值結果,參考前人的研究(Benestyet al,2005),進行了一系列數值模擬實驗。

模擬時間序列Xt,Yt由下式產生

其中n,m,r為正整數,Ai,Bi,Ci,Di為有理數,ε1,t,ε2,t是獨立的高斯白噪聲過程且其期望為0。在Yt中,相位φ1,φ2,???,φn是隨機變量。Xt和Yt的相同頻率為fi,i=1,2,???,n,Xt的特有頻率為fjx,j=1,2,???,m,Yt的特有頻率為fjy,j=1,2,???,r。則理論凝聚譜為

這里給出三對模擬結果:

第一對序列長度取1024,采樣頻率為1Hz,Ai,Bi,Ci,Di均為1。選取兩者的相同頻率為f1=0.10,f2=0.20,f3=0.21,f4=0.22,f5=0.23,f6=0.24;f7=0.45,Xt的特有頻率為f1x=0.07,Yt的特有頻率為f1y=0.32。其中頻率f2、f3、f4、f5和f6非?？拷?這樣便于考察各方法的頻率分辨能力。分別采用Welch方法(Hanning窗,窗長度256)、MTM方法 (分辨率帶寬4/1024)、MVDR方法(窗寬 256,分辨率為 400) 進行凝聚譜估計,并做非零相關性檢驗。在顯著性水平0.05下,得到的臨界值為0.41,具體結果如圖1中a、b、c所示。

圖1 凝聚譜估計 (三對數值模擬)Fig.1 Estimation on magnitude-squared coherence (numerical simulation on three pairs of data)a、b、c分別為Welch、MTM、MVDR方法對第一對模擬數據的凝聚譜估計; d、e、f分別為Welch、MTM、MVDR方法對第二對模擬數據的凝聚譜估計;g、h、i分別為Welch、MTM、MVDR方法對第三對模擬數據的凝聚譜估計

第二對改變序列長度為4096。做非零相關性檢驗,在顯著性水平0.05下,臨界值為0.10,具體結果如圖1中d、e、f所示。

第三對序列增加特有頻率的個數,Xt的特有頻率為f1x=0.07;f1x=0.30;f1x=0.40,Yt的特有頻率為f1y=0.01;f1y=0.15;f1y=0.48。做非零相關性檢驗,在顯著性水平0.05下,臨界值為0.41。結果如圖1中g、h、i所示。

相關信號的檢測能力可通過觀察顯著相關的頻率是否能夠被檢測出來進行判斷。當凝聚譜估計值大于臨界值時,相應分量是顯著相關的。由圖1可以看出,三種方法均能檢測出存在相關性的頻率,但Welch方法和MTM方法的估計結果存在偽峰。

可通過考查相近且相關程度相同的頻率分量并觀察其凝聚譜估計曲線是否存在兩個對應的可分辨的尖峰來衡量分辨能力。由圖1可知,在相近的頻率f2、f3、f4、f5、f6上,MVDR方法的估計譜能夠清晰地分辨每一個頻率點,每個頻率點上的峰狹窄且互不交疊; Welch方法的估計譜可以勉強分辨出存在相關的頻率點,但是一些譜峰連在一起; MTM 方法的估計結果也能夠顯示出具有相關性的頻率點,但其峰較寬,分辨率介于MVDR和Welch兩種方法之間。

此外,在第一對數值模擬的基礎上,從改變模擬數據長度、變換共有頻率、增加特有頻率的個數、改變振幅系數、加大噪聲方差等方面進行了系統的模擬。結論與以上實驗類似: MVDR方法的分辨率最高,不易出現偽峰; Welch方法分辨率較差,容易出現偽峰; MTM方法的估計性能介于兩者之間。

3 實驗室波流數據分析

3.1 數據及其預處理

為了解以上三種方法在具體應用中的差異,我們在中國海洋大學物理海洋實驗室的大型風-浪-流水槽中進行了波流實驗。水槽長65m,寬1.2m,高1.5m,實驗時水深0.7m左右。分別選定距離出風口30m、34.2m、39.7m三個風區位置進行實驗。在每個風區位置,風速分別設定為6.35、7.0、7.5、8.0和9.0m/s,實測風速在設定水平上下 5%—10%范圍內波動。波面位移采用實驗室自主研發的 OUC-2型鉭絲測波儀(以下簡稱“測波儀”)測量,采樣頻率約為26.78Hz。水下的流速采用聲學多普勒流速測量儀(以下簡稱“ADV”)測量,可探測探頭下方 10cm 位置的三維流速,測量深度分別位于15、20、25、30、35和40cm,垂直方向有重復觀測,采樣頻率為25Hz。實驗時,先開動流機,機械造波 5分鐘,然后攪動水槽內的水體,使水中雜質均勻散布于水體中。為便于ADV測量,流機流速范圍設定為 0—40cm/s。之后停止機械造波,待水面平靜后,設定風速,并啟動風機造風。波面起伏平穩后,同步開始測波儀、ADV的數據采集,測量記錄波面位移與水下給定深度處的流速。受到測波儀處理記錄的限制,每次記錄時間長度為10分鐘(類淑河,2010)。

不失一般性,選取其中三對實驗數據信號進行研究。第一對是風區長度39.7m、風速6.35m/s、ADV探頭深度 20cm情形下的波面位移和鉛直流速信號,記為pair I。第二對是風區長度34.2m、風速8.0m/s、ADV探頭深度 25cm情形下的波面位移和鉛直流速信號,記為 pair II。第三對是風區長度 30m、風速9.0m/s、ADV探頭深度40cm情形下的波面位移和鉛直流速信號,記為pair III。其中波面位移時間序列有15000個點,鉛直流速時間序列約有14000個點。為進行交叉譜分析,對鉛直流速序列進行線性插值,使其與波面位移的15000個點一一對應。鑒于流速數據不平穩,所以采用 Butterworth濾波器濾去頻率小于0.5Hz的信號,即去除低頻信號,然后對波面位移和預處理之后的鉛直流速序列進行零均值化處理。

3.2 凝聚譜估計

由于實驗室數據足夠長,研究中首先將數據均進行分段。每段包括1024個數據點,且相互無重疊,由此共分為14段。對每段數據分別用Welch方法、MTM方法、MVDR方法進行凝聚譜估計,參數選擇與模擬數據分析相同。最后對14段凝聚譜估計取平均,并進行非零相關性檢驗,在顯著性水平 0.05下,得到的臨界值為 0.054。三對信號的平均凝聚譜估計結果如圖2所示。

不難看出,Welch和MTM方法的估計結果相近,在整個可識別頻率范圍內凝聚譜的估計值均高于臨界值,是顯著相關的。MVDR方法給出的估計在譜峰頻率fp附近即主導波頻段內(Babanin,2009)與 Welch和MTM 方法的估計結果走勢一致,估計值略低,但在高頻段,大約在1.5fp之后,有一個明顯的指數衰減趨勢,而在3fp之后,凝聚譜估計值小于臨界值,相關性不顯著。

分別對三對實驗信號采用Welch、MTM、MVDR方法估計得到的14段凝聚譜對應譜值求方差后再取平均,得到各自的平均方差。第一對信號三種方法對應的平均方差分別為 0.0318、0.0272、0.0178; 第二對信號分別為 0.0260、0.0208、0.0139; 第三對信號則為 0.0319、0.0266、0.0194。顯然,MVDR方法得到的平均方差最小,此方法最穩定。

圖2 三對風浪實驗信號的凝聚譜估計Fig.2 Estimation on magnitude-squared coherence of three pairs of wind-wave experiment data a、b、c分別為 pair I、pair II、pair III 實驗信號的凝聚譜估計

三種方法給出的凝聚譜估計均表明波流信號在主導波頻段上存在顯著的相關性。然而高頻段的信號極容易被噪聲干擾甚至掩蓋。就此實驗數據而言,并不能判斷高頻段波流信號是否一直存在顯著的相關。另一方面,模擬數據分析發現,MVDR方法提供的估計有更強的頻率識別與分辨能力; 且實驗數據分析表明,MVDR方法得到的平均方差最小,估計結果更穩定。因而,在這三種方法提供的凝聚譜估計中,MVDR方法的估計結果更為可靠,即波流信號僅在主導波波段內存在顯著相關性。

圖 3中分別給出了三對數據的波面位移單譜估計、鉛直流速單譜估計和 MVDR方法給出的凝聚譜估計(為便于比較,估計結果均取對數)。由圖可知,波面譜只有一個尖峭的主峰,與鉛直流速譜譜形相似。凝聚譜與兩個單譜的譜形走勢大體一致,但譜峰較寬,在稍高于譜峰對應的頻率的較寬的頻段內,仍保持較高的凝聚譜譜值。這表明,在這一頻率范圍內,表面波動與水下鉛直流速存在很強的相關性。

4 討論與結語

圖3 三對風浪實驗信號的波面位移譜估計、鉛直流速譜估計、凝聚譜估計Fig.3 Estimation of wave spectrum,current spectrum,and magnitude-squared coherence of three pairs of wind-wave experiment data a、b、c分別為pair I、pair II、pair III的波面位移譜估計; d、e、f分別為pair I、pair II、pair III的鉛直流速譜估計; g、h、i分別為pair I、pair II、pair III的凝聚譜估計

本文討論了Welch、MTM和MVDR三種交叉譜估計方法,在波流信號相關性分析中的應用,并通過數值模擬,從檢測相關信號能力和分辨率能力兩個方面評估了三種估計方法的性能。我們發現,Welch方法能夠鑒別出相近的相關頻率,但分辨率不高,且容易出現偽峰; MTM方法出現的偽峰較Welch方法少,分辨率能力仍較差; 而 MVDR方法既能夠檢測到相關信號,又不容易出現偽峰,分辨率較高,在三種方法中估計性能最優。

分別采用三種方法估計了實驗室三對波面位移與鉛直流速信號的凝聚譜。結果表明,在以譜峰頻率為中心的主導波頻段,三種方法給出的凝聚譜估計結果相近,均高于臨界值,表明波流信號在此頻率區間內顯著相關。在高頻段,估計結果存在差異。在高于1.5fp之后的頻段,MVDR給出的估計有衰減趨勢,在 3fp之后變得不顯著,而 Welch方法和 MTM 方法給出的凝聚譜卻在整個高頻段均大于臨界值,即使在很高的頻段上,波流信號仍顯著相關。本文研究表明,MVDR的估計結果更加合理。

在進行交叉譜估計和單譜估計時,各種估計方法都面臨諸多參數選擇問題,如Welch方法中窗寬、重疊比例與窗類型的選取。選擇的參數不同,估計結果也會有差異。Priestley(1981)指出,相對于窗類型和重疊比例的選擇,截斷點即窗寬的選取對估計結果影響更大。我們數值模擬的結果也表明,選擇不同的窗、變換重疊比例等,估計結果變化不大,但改變窗寬,得到的估計結果會有明顯不同。另外,在進行非零相關性檢驗時,臨界值是通過Welch方法估計的自由度進行檢驗的,而對于MTM方法和MVDR方法來說,這一臨界值是否適合,還有待進一步驗證。

在實驗室波流信號相關性分析的討論中,僅采用了具有代表性的三對數據,得到的一些物理海洋方面的結論可能不具有一般性,需要在未來的工作中深入研究。

王 濤,李家春,1999. 波流相互作用研究進展. 力學進展,29(3): 331—343

類淑河,2010. 風浪破碎的隨機性、非線性和能量耗散特征研究. 青島: 中國海洋大學博士學位論文

袁耀初,趙進平,王惠群等,2002. 南海東北部450m以淺水層與深層海流觀測結果及其譜分析. 中國科學(D輯),32(2):163—176

高 山,2003. 二維數值波浪水槽模式的建立和應用及浪流相互作用研究. 青島: 中國海洋大學博士學位論文

黃建國,1991. 現代譜分析的發展及應用(二). 數據采集與處理,6(1): 33—45

Babanin A V,2009. Breaking of ocean surface waves. Acta Physica Slovaca,59(4): 305—535

Benesty J,Chen J D,Huang Y A,2005. A generalized MVDR spectrum. IEEE Signal Processing Letters,12(12): 827—830

Bortel R,Sovka P,2007. Approximation of statistical distribution of magnitude squared coherence estimated with segment overlapping. Signal Processing,87(5): 1100—1117

Capon J,1969. High-resolution frequency-wave number spectrum analysis. Proceedings of the IEEE,57(8):1408—1418

Carter G C,1987. Coherence and time delay estimation.Proceedings of the IEEE,75(2): 236—255

Gallet C,Julien C,2011. The significance threshold for coherence when using the Welch’s periodogram method:Effect of overlapping segments. Biomedical Signal Processing and Control,6(4): 405—409

Jenkins G M,Watts D G,1968. Spectral Analysis and its Applications. San Francisco: Holden Day,256—300

Kay S M,1988. Modern spectral estimation. New Jersey:Prentice Hall,446—471

Koopmans L H,1974. The Spectral Analysis of Time Series.New York: Academic Press,281—285

McKone K P,2003. Multiple Methods of Spectral Analysis with Applications to the Florida Current. Ph.D. Thesis,Dept. of Marine Science,University of Southern Mississippi,Hattiesburg,MS. pp 120

Percival D B,Walden A T,1993. Spectral Analysis for Physical Applications. Cambridge: Cambridge University Press,289—295

Priestley M B,1981. Spectral Analysis and Time Series. London New York: Academic Press,654—726

Strand O N,1977. Multichannel complex maximum entropy(autoregressive) spectral analysis. IEEE Transactions on Automatic Control,22(4): 634—640

Thomson D J,1982. Spectrum estimation and harmonic analysis.Proceedings of the IEEE,70(9): 1055—1096

Welch P D,1967. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short,modified periodograms. IEEE Transactions on Audio and Electroacoustics,15(2): 70—73