基于EVSEWT的下擊暴流時變平均風提取*

李春祥， 李 洲

(上海大學土木工程系 上海，200444)

引 言

極端風(熱帶氣旋、下擊暴流、龍卷風等)特性研究對于提高結構抗風設計水平、完善風荷載規范等意義重大。與大氣邊界層(atmospheric boundary layer，簡稱ABL)風不同，極端風具有明顯的非平穩性，即風信號的統計特征值具有時變性。為了便于分析，學者們通常將非平穩風視為時變平均風與零均值隨機脈動風的疊加[1]。在非平穩風特性分析中時變平均風提取的準確與否會影響時域中的相關分析或頻域中的功率譜分析，從而產生較大誤差，甚至使低頻譜完全失去真實性。另外，在非平穩脈動風模擬與預測、結構非平穩振動響應分析和結構參數識別等研究中都涉及到信號時變平均項的提取[2]。因此，合理時變平均風速的提取對于非平穩風研究具有重要意義。

在時變均值提取的各類方法中，平均斜率法、差分法、最小二乘擬合等通常需要預先知道信號趨勢項的類型，對處理具有復雜變化趨勢的信號并不適用[3]。線性擬合、高階濾波和滑動平均等缺乏統一的定量提取方法或準則，提取效果依賴于使用者的經驗以及計算參數的選擇，影響了模型的實用[4]。與此同時，由于經驗模態分解(empirical mode decomposition，簡稱EMD)具有自適應性，不需要先驗信息，使用方便；離散小波變換(discrete wavelet transform，簡稱DWT)在處理復雜信號上具有良好的時頻分辨率，故這兩種方法得到了學者們的普遍使用。陳雋等[4-5]提出使用EMD方法，并引入脈動風平穩度標準和截止頻率，來提取臺風風速記錄中的時變平均風。申建紅等[3]基于人工模擬風速數據，提出利用單尺度小波能量的突變確定時變均值分解的準確層次，較EMD分解方法更加準確有效。Su等[6]基于下擊暴流實測數據，針對Kernel回歸、DWT和集合經驗模態分解(ensemble empirical mode decomposition，簡稱EEMD)方法給出了選擇合適的方法和時間窗大小以獲得合理的時變平均值的建議，推薦采用高階的DWT方法用于下擊暴流時變非平穩風提取。Wang等[7]將DWT方法用于時變均值提取，用非平穩風速模型對蘇通大橋風場數據進行了分析，分析結果與早期結果吻合。但是EMD方法仍存在模態混疊、端部效應和缺乏數學理論等不足，DWT方法對信號分解層次和母小波類型的選擇有較嚴格的條件，相關參數的確定需要依據經驗。另一方面，目前的非平穩風研究雖然取得一定進展和突破，但仍存在以下一些挑戰：由于下擊暴流、龍卷風的發生具有偶然性，相關實測樣本仍不豐富；對極端風機理揭示不足，數學分析理論有待創新；傳統方法仍需要憑借豐富的經驗，新一代精細化時頻分析技術急需研究。

本研究的重點是提出一種基于能量的波谷尋找經驗小波變換方法(energy valley searching empirical wavelet transform，簡稱EVSEWT)，并給出了時變平均風的確定準則。首先，簡要介紹了時變平均風提取的常用方法，包括EMD，EEMD和DWT方法；其次，針對時變平均風的特征，對經驗小波變換(empirical wavelet transform，簡稱EWT)進行改進，提出EVSEWT，著重于提高對低頻信息的分辨率；最后，給出了時變平均風提取的結論和關于EVSEWT的改進和應用方面的展望。為避免以往研究只在信號層面進行效果驗證，筆者將提取的時變平均風和脈動風用于結構非平穩響應分析，通過2組下擊暴流風速試驗，對比EMD，EEMD，DWT(db10，db18)和EVSEWT的效果。

1 數據與方法

1.1 數據說明

2002年6月4日和15日，在圖1所示的沿直線均勻分布的7座觀測塔上，同時記錄兩組不同高度的下擊暴流風速時間序列。測量了超過40 m/s的陣風風速。一個記錄的下擊暴流被歸類為后翼下沉氣流(rear-flank downdrafts,簡稱RFD)，另一個歸類為Derecho。該實測數據具有強烈非平穩特性，其中采樣頻率為1 Hz，總時長為30 min，具體實測信息可參考文獻[8]。筆者選取塔4上高度15m處的RFD數據和Derecho數據，如圖2所示。

圖1 下擊暴流測點布置示意圖Fig.1 Diagram of measuring point arrangement of downburstL

圖2 兩組下擊暴流實測數據Fig.2 Two sets of downburst data

1.2 常用方法

在時變均值提取的諸多方法中，本章主要介紹離散小波變換、經驗模態分解和經驗小波變換方法。

離散小波分析法是用小波基函數逼近信號的時頻分析方法。小波基函數是一組由快速衰減、有限長的母小波經縮放和平移而生成的函數序列，常用的母小波函數有db族小波、Harr小波、sym小波、Morlet小波和Meyer小波等。但離散小波分析法存在分解層數和小波母函數選擇的困難，對不同特征的數據，母函數的效果也有較大差距，常需要按經驗選擇。

EMD方法是一種信號分析方法，依據數據的時間尺度特征進行信號分解，無須預先設定基函數。它將一個信號分離為多個本征模態函數(intrinsic mode function,簡稱IMF)與一個余量R之和。本征模函數視為信號中的振動成分，反映了信號中不同頻率的成分。從原始信號中先分離出來的本征模函數較后分離出來的頻率高，余項的頻率最低，為趨勢成分，這樣通過選擇性地疊加低頻分量就可得到原始信號的低頻趨勢項。EMD方法使用簡單，適用于大多數非平穩信號分析，但EMD分解的結果會出現端部效應和模態混疊的問題，同時存在缺乏數學基礎理論、分解的子信號物理意義不明確等不足。為了改進EMD的端部效應和模態混疊，Wu等[9]提出集成經驗模態分解，通過在原信號中加入合適比例的白噪聲，來補充一些缺失的尺度，在機械、土木、航空航天等領域的信號分解中具有良好表現。

在時變均值提取的應用中，Holmes等[10]通過兩個準則確定移動平均法中的窗口大小，即時變平均風速需要反映原始風速低頻部分的趨勢、峰值以及脈動風的均值應趨近為0。Su等[6]根據如下3個準則確定合理的小波基函數和分解層數：①時變平均風速應反映原始風速的變化趨勢；②對應脈動風的演化功率譜密度(evolutionary power spectral density,簡稱EPSD)估計結果應具有物理意義的解釋；③EPSD作用下的高層建筑順風向結構響應是否相對保守。為了簡便，本研究中DWT方法按照尺度小波能量的突變確定時變均值的準確層次[3,11]；EMD和EEMD方法取分解結果中小于截止頻率fc的子信號為時變均值項，通過間歇檢測準則來完成，截止頻率確定取決于高通信號的平穩性，參考文獻[12]。

1.3 經驗小波變換

為了解決小波變換參數選擇困難、經驗模態分解缺乏數學理論以及模態混疊等問題，Gilles[13]提出了一種經驗小波變換。這種方法可以看作是構造1組帶通濾波器和1個低通濾波器，用于將原信號分解成多個頻率均一、穩定的子信號。具體實現流程如下：

1) 如圖3所示，對原信號進行快速傅里葉變換，得到原信號標準化傅里葉譜(橫軸標準化至0～π)；

圖3 標準化傅里葉譜的劃分Fig.3 Classification of standardized Fourier spectra

2) 在傅里葉譜上確定各帶通濾波器的邊界，并基于這些邊界，利用Meyer小波構造濾波器組；

3) 原信號通過低通濾波器和這些帶通濾波器，被分解為頻率從高到低的若干子信號集。

1.4 EVSEWT介紹

盡管EWT已經給出了幾種傅里葉邊界劃分的策略，如“locmaxmin”，“adaptive”和“scalespace”等，但由于應用的目的和場景不同，難以結合時變平均風的提取。不過，EWT方法有完整的數學理論，筆者基于EWT的理論框架，在傅里葉譜邊界劃分上更注重低頻頻段的提取，提出EVSEWT，圖4所示為其提取時變趨勢項的流程，具體如下：

圖4 基于能量的波谷尋找經驗小波變換原理圖Fig.4 Schematic diagram of energy-based valley searching empirical wavelet transform

1) 對原始信號計算進行快速傅里葉變換，獲得標準化傅里葉譜；

2) 在傅里葉譜的0～X之間，對每一處波谷設立劃分邊界，邊界數量自適應確定；

3) 計算相鄰兩邊界與傅里葉譜圍成的面積，作為子信號能量指標，并按照能量從大到小對子信號進行排序；

4) 依據步驟2確定的邊界，利用Meyer小波構造濾波器，按能量從大到小逐次提取子信號放入時變趨勢項中，剩余部分子信號作為脈動分量；

5) 利用Runtest法確定脈動風的平穩性；

6) 重復步驟4和5，直到脈動風為平穩數據且均值趨近于零，若對0～X的信號處理完后，脈動風的平穩性要求未能滿足，則只取0～X之間頻段作為時變平均風。

在結構動力學分析中，一般認為平均風速最高頻率為結構基頻的1/5～1/10時，可以忽略其動力放大效應，按擬靜力分析[6]。因下面分析的結構基頻為0.23 Hz，故平均風速最高頻率為0.023 Hz。對采樣頻率為1 Hz的信號來說，0.023 Hz在標準化傅里葉譜上對應為0.144 Hz。筆者確定X為0.144 Hz，時變趨勢項的頻率控制在0～X范圍中，且X可以視使用情況調整。總體來講，EVSEWT注重低頻部分的分解，即0～1/10結構基頻處，這樣提取的時變均值項滿足按擬靜力分析的要求，也能較好地刻畫實測數據的趨勢。

2 非平穩隨機響應分析

2.1 理論介紹

將實測風速分為時變平均風和脈動風

(1)

(2)

其中：M1為結構第1階振型的廣義質量；ω1為結構第1階振型的角頻率。

采用精確且高效的虛擬激勵法對結構的隨機動力響應進行分析，結構的虛擬激勵力可以表示為

(5)

(6)

虛擬位移y(ω,t)可以由Newmark-β逐步精細積分法計算，參見文獻[15]，則廣義位移脈動分量的EPSD矩陣Sq′(ω,t)=y(ω,t)y(ω,t)T。當非平穩風荷載的EPSD時間變化率大于結構自然周期時，則可以采用擬平穩假設計算Sq′(ω,t)，結果更偏保守，計算快速

Sq′(ω,t)=|H(ω)|2SQ′(ω,t)

(7)

其中：H(ω)為基本振型下的頻率響應函數[16]，H(ω)=1/[M1(-ω2+2iξ1ω1ω+ω2)]。

(10)

其中：ρ為空氣密度，取1.225 kg/m3；CD為阻力系數，取1.2；B為建筑物長度；χD(ω)為氣動導納函數；g(z,ω,t)為調制函數；dΘ(z,ω)表示復值零均值正交增量隨機過程。

筆者選取Davenport氣動導納函數

|χD(ω)|2=2/λy(1-1/λy+1/λye-λy)

其中：λy=kyfD/Umax；ky=8為衰減因子；f為頻率；Umax為時變平均風速的最大值。

當只考慮結構的第1階振型時，廣義力的時變均值和脈動分量計算如下

(11)

(12)

其中：φ(z)=(z/H)β為第1階振型的振型函數；β為振型指數，當振型為線性時取為1。

將式(9)代入式(11)可得

(13)

廣義力的脈動分量Q′(t)的EPSD可表示為

當脈動風的EPSD假定為沿著建筑高度不變時，即Su′(z1,ω,t)=Su′(ω,t)，式(14)可表示為

(15)

若時變平均風剖面進一步簡化為

(16)

此時，時變平均風速引起的廣義力可表示為

(17)

相應的，式(15)可以簡化為

(18)

其中：|Jz(ω)|2為聯合接收函數。

(19)

由于實測樣本少以及在數學上處理的困難，在大多數研究中會采用時不變指數函數或Davenport提出的相干函數。筆者假設z1和z2高度處的相干函數是時不變的

(20)

其中：kz為衰減因子，取值為8。

本研究中，下擊暴流風速的豎向風剖面選取Wood[17]提出的時不變經驗模型

(21)

以上數值風剖面經驗公式已經通過實測風速數據驗證[18]，在其他研究中也被廣泛采用[6,16]。由于目前尚缺乏深入的下擊暴流風速知識以及足夠數目的觀測樣本，雖然上述假定在研究中會產生不確定性，但是這些簡化處理不會影響一般性的規律與結論。

除此之外，采用基于穿越率的極值分布理論可以得到響應的平均極值和標準差[19-20]。設R(t)為時間周期T內的非平穩響應，根據穿越率的泊松近似關系，在時刻t和水平r的平均穿越率可由下式計算

(22)

由文獻[21]可知，非平穩響應R(t)在某時間周期T內的極值的累積分布函數如下

(23)

(24)

2.2 結果分析

為了驗證EVSEWT在時變均值項提取中的有效性，將EVSEWT用于RFD和Derecho兩組實測下擊暴流風速的時變均值提取試驗，并選取EMD，EEMD和DWT做對比分析。根據文獻[6]的建議，母小波函數選擇db10和db18，最后將提取的結果用于某幢高層鋼結構建筑的隨機動力響應分析。

已有的研究表明[20-21]，下擊暴流風速會隨地面高度增加而增加，在50～200 m高度處達到峰值，再往上風速下降緩慢。因此，下擊暴流可能對懸臂結構產生顯著的風荷載影響[18]。下面以一座高H=200 m，長B和寬D均為40 m的結構為例。結構密度設為192 kg/m3，結構的基頻f1=46/H=0.23 Hz，振型阻尼比ξ1=1%。對于高層建筑，風振響應一般以基本振型為主，高層振型的貢獻可以忽略，筆者只考慮第1振型的影響，且振型形狀假定為線性函數。具體地，以結構頂端均值位移、基于擬平穩和非平穩的RMS位移等結構響應值以及結構非平穩響應的平均極值作為評價指標。

簡潔起見，只給出了RFD的情況作圖示說明。圖5是由EMD，EEMD，DWT(db10)，DWT(db18)和EVSEWT方法，依據各自的準則提取的時變平均風速和脈動風速。圖6是脈動風速的EPSD圖，EPSD圖刻畫了信號在不同時刻、不同頻率上的能量分布，一般認為大多數嵌入在脈動風中的能量分布在0～0.1 Hz的低頻范圍內。同時，在一定時間內，隨著頻率的增加，大多數EPSD都呈現出衰減的趨勢。脈動風EPSD也是后續結構響應分析中需要用到的物理量之一。從圖2中RFD的實測數據可以看出，信號的能量主要集中在150～700s之間，這與圖6結果相吻合，且能量集中在0～0.1 Hz頻段，與經驗相符。另外，從圖6的結果來看，EMD和EEMD的低頻部分能量最高，可能導致更大的結構動力響應，EVSEWT次之，DWT的最小。對應圖5中提取的時變平均風，DWT方法的波動比較明顯，EMD和EEMD的波動平緩，EVSEWT介于兩者之間。

圖6 不同分解方法對應脈動風的演化功率譜密度Fig.6 Evolutionary power spectral density of fluctuating wind speeds extracted by different methods

圖5 不同分解方法獲得的平均風速和脈動風速Fig.5 Time-varying mean and fluctuating wind speeds obtained by different decomposition methods

與ABL風相似，非平穩風激發的結構響應一般分為靜風響應和動力響應兩部分。前者由時變平均風速引起，通過擬靜力方法簡單計算；后者由平穩/非平穩脈動風引起，通過動態分析進行計算[22]。圖7是不同信號分解下，RFD風速在結構頂端產生的均值位移(Mean)、基于擬平穩(Pseudo-stationary)和非平穩的(RMS)位移。可以發現，均值位移與提取的時變平均風的趨勢類似，DWT方法的波動較大，EMD和EEMD方法的均值位移相對平緩。需要指出，基于非平穩假設計算的RMS位移的最大值略小于基于擬平穩假設的方法，且在時間上，非平穩響應最大值要滯后于擬平穩響應和均值位移的最大值，這與文獻[6,15]的結論一致。

表1和表2分別是RFD和Derecho兩組下擊暴流風速在不同方法下得到的結構200 m高度處的均值位移最大值、RMS位移最大值和平均極值位移。綜合圖7、表1和表2，可以得出以下結論：①在RFD風速中，最大相差發生在RMS位移，為11.24%(EMD法和EVSEWT法之間)，均值位移相差最大達9%(EMD法和EVSEWT法之間)，相對來看，Derecho風速中各方法相差較小，為RMS位移的7.81%(EMD法和EVSEWT法之間)，說明EVSEWT方法對結構響應分析相對保守；②從3項位移指標來看，2組試驗均表現出EVSEWT方法最佳、DWT方法次之、EEMD和EMD略差的結果，這可能是由于EVSEWT方法是基于信號傅里葉譜進行的分解，在信號低頻部分具有更好的分辨率，配合相應的準則能比較合理地提取時變平均風速。

圖7 頂端位移均值和RMS位移Fig.7 Mean and RMS responses of tip displacement

表1 試驗1(RFD)結構頂端位移的最大均值、最大RMS和平均極值位移

表2 試驗2(Derecho)結構頂端位移的最大均值、最大RMS和平均極值位移

3 結 論

1) EWT提取的時變均值能兼顧實測信號趨勢和高頻波動，對低頻信息的分辨率更優秀。

2) 在時變平均風極大值和脈動風的能量都不突出的情況下，能計算得到偏于安全的結構響應位移，說明EVSEWT的確定準則更加合理。

3) 3個響應位移之間相互關系符合經驗[19,23]，即非平穩RMS位移極大值略小于擬平穩RMS位移極大值，且滯后于擬平穩RMS位移極大值和均值位移最大值。

4) 總體來說，EVSEWT具有非常好的信號低頻分辨率，對2組試驗數據的時變均值提取試驗效果較好。后續的研究可以考慮用于橋梁和高層建筑結構健康監測信號時頻分析[24]、機械故障診斷和生物醫學信號處理等領域。另一方面，受到多元經驗模態分解的啟發，將EVSEWT拓展到空間信號處理問題，對于注重空間相關性的風場分析也很有意義。

致謝本研究風速數據來源于德州理工大學風工程研究中心，在此表示衷心感謝。