基于演化相位譜的脈動風速模擬

閻 啟，李 杰，2

(1.同濟大學 建筑工程系，上海 200092;2.同濟大學 土木工程防災國家重點實驗室，上海 200092)

結構在風荷載作用下的計算一般有頻域和時域兩種方法。頻域方法速度快，但是只能用于線性階段的計算;而對于輸電線路、大跨橋梁、高層建筑等具有明顯非線性特征的結構，時域分析能夠更為準確的把握結構響應特征。

時域分析采用的風速時程通常由數字模擬的方法生成［1]，其中譜表現法(Spectral representation method)是一種較為成熟的數字模擬方法。經歷了從1954年［2]至今的發展，譜表現法從最初只能模擬一維、單變量、平穩的隨機過程，已發展為可以廣泛的模擬多維、多變量、非平穩隨機過程的方法［3－5]，并在脈動風場模擬方面得到了廣泛的應用［1，6]。譜表現法的實質，是利用具有隨機相位的諧波疊加、對具有目標功率譜的隨機過程進行模擬。但仔細分析發現:譜表現法在基礎上至少存在以下三點問題。首先，譜表現法本質上是基于二階矩的模擬。在進行隨機過程模擬時，各諧波的幅值是利用功率譜密度進行計算，對目標隨機過程的模擬精度也是通過功率譜密度的符合程度來體現。由于功率譜密度在本質上屬于隨機過程的二階數值特征，因此，所模擬隨機過程的概率信息不會高于二階。其次，譜表現法存在無法重現樣本的困難。由于功率譜本質上屬于集合特征［7]，因此基于功率譜模擬的隨機過程無法重現現實中的隨機過程樣本。再次，譜表現方法中需要用到初始相位，通常是在不同頻率處取區間［0，2π)內均勻分布并且相互獨立的隨機變量作為初始相位［3－5]。這種做法使得在實際風場模擬中隨機變量的數目巨大(400～600個)，而不同頻率相位之間、相位與能量之間的可能的內在聯系并未被考慮。上述問題帶來的后果是:應用譜表現法得到的風速時程進行結構隨機響應分析，不僅計算工作量大，也很難得到響應的概率密度信息，進而，造成了結構可靠度分析的困難。

隨機Fourier函數模型的提出為動力激勵的建模與模擬提供了新的途徑［8]。對地震動、脈動風速以及海浪等動力激勵的時程，應用Fourier變換，可以將其分解為幅值部分和相位部分，稱為Fourier幅值譜和Fourier相位譜。幅值譜是對動力激勵能量的一種譜分解，基于大氣湍流的物理機制，可以得到脈動風速隨機Fourier幅值譜的物理模型［9]。相位譜則主要控制時程的波形、反映過程的非高斯性［10]。作為隨機Fourier譜中十分重要的一部分，建立合理的相位譜模型，對于脈動風速的模擬以及結構抗風可靠度的計算，具有至關重要的意義。

本文首先從概念上明確了相位譜和相位譜主值的區別，基于對風場湍流物理背景的分析，提出了相位演化速度的概念，并由此建立了相位的零點演化時間Te和相位譜一一對應的關系，依據大量實際風速測量數據，建立了Te的概率模型。結合脈動風速的Fourier幅值譜、應用Fourier逆變換進行了脈動風速的模擬。與實測結果的對比，表明了本文建議方法的可行性。

1 相位譜和相位譜主值

設脈動風速時程為u(t)，其離散Fourier變換為:

其中n為自然頻率，T為脈動風速的持時，d t為采樣間隔，與采樣頻率Fs互為倒數，F(n)稱為脈動風速時程u(t)的Fourier譜。乘以是為了取單邊譜。F(n)是一組復數，可以寫為幅值和相位的形式:

圖1 實測脈動風速時程及其Fourier幅值譜和相位譜Fig.1 Measured fluctuating wind speed and corresponding Fourier amplitude spectrum and phase spectrum

2 演化相位譜

2.1 相位的演化速度

由于測量條件的限制，很難得到某一時刻空間中各點的湍流速度。Taylor凍結假定將以空間變量為背景的湍流假設為空間一點以時間為變量的湍流［12]。根據這一假定，同一空間點、不同時刻的脈動風速時程(自變量為時間)和同一時刻、不同空間點的脈動風場分布(自變量為位置)是相同的。這里，不同空間點指沿主風方向長度為L的范圍，L=UT，U為平均風速，T是測量脈動風速的時距。自Taylor凍結假定提出后，幾乎所有的湍流測量都以此為基礎［13]。

眾所周知，風速儀記錄到的是空間一點處不同氣體質點的流動速度，風速記錄減去平均風速即得到脈動風速時程。由Taylor假定可以推論:如果風速儀以和平均風速同樣的速度沿主風方向前進，那么將記錄到同一氣體質點在空間不同位置的速度，該速度與前述脈動風速時程相同。此時，脈動風速時程可視為是同一氣體質點在一以平均風速值大小前進的參考系內的縱向“振動”速度。

在物理上，這樣的空氣質點振動速度可視為是一系列不同尺度、不同頻率渦旋振動的疊加。渦旋的特征振動速度可以表示為［14]:

對所有頻率范圍內特征振動速度的平方求和，便得到脈動風速時程的能量總和即方差σ2:

特征速度為v(n)的渦旋從時間t0至t1行進的距離與其波長之比，表征了該時間間隔內渦旋變化的周期數，每個周期對應2π的相位變化，如圖2所示。不同尺度、即不同頻率渦旋在這一時間間隔τ內的相位改變Δφ(n)可寫為:

其中:τ=t1－t0;l(n)為渦旋的波長，表征了渦旋的尺度大小;k(n)為波數，與l(n)互為倒數;波數與自然頻率的關系為:

這里U為平均風速。

對式(5)關于時間t求導，便可得到不同頻率渦旋的相位演化速度:

圖2 渦旋相位變化示意圖Fig.2 Schematic diagram of the change of eddy’s phase angle

2.2 相位的零點演化時間

不同頻率渦旋振動速度不同，因而相位演化速度也不同。一般來說，低頻、大尺度的渦旋相位變化慢，高頻、小尺度的渦旋相位變化快。可以設想，真實脈動風速可視為是由一簇具有相同初始相位的渦旋(諧波)經過時間Te的演化后疊加而成的。最簡單的初始相位取值是零相位，如圖3所示。記錄到的某一段風速時程，可以看作是一批具有相同初始零相位的渦旋經過時間Te演化而來。我們稱Te為零點演化時間，單位為s。

圖3 相位演化時間Te的示意圖Fig.3 Schematic diagram of phase evolution time Te

應用實測的風速時程，可以識別出Te的具體數值。本文采用10 Hz采樣、10分鐘持時的脈動風速樣本進行識別，識別原理如下。

首先根據式(1)和(2)求取樣本風速時程的Fourier幅值譜和相位譜φ(n)，應用式(7)，可得到不同頻率渦旋的相位演化速度Δ·φ(n);以時程開始時刻為起點沿時間軸t向反方向推進(為計算方便，t仍取正值)，用下式求各頻率渦旋的相位:

求得相位φ(n，t)后需要換算到［0，2π)的主值區間中，得到相位主值φm(n，t)。若某一時刻所有渦旋相位主值均為零，則對應此時的t值即為所要識別的零點演化時間Te。

考慮理論抽象與真實問題背景之間的差距，本文推薦采取下述零點演化時間識別方法。

設最低頻渦旋的初始相位為φ(n1，0)，特征速度為v(n1)，以時程開始時刻為起點向反方向推進，用下式求得一系列時間點T0(T0仍為正值)，使得此時相位主值 φm(n0，T0)=0

這里p為零點個數搜索范圍，是自然數的子集，從1開始直到一個計算量允許的數，每個p值對應一個T0。在各個時間點T0，用下式計算前10個頻率點(對應頻率為1/600 Hz到1/60 Hz)的相位值:

受計算量的限制，搜索范圍設置為p≤106。在所搜索范圍內，如果T0滿足最大偏差條件:

以及加權約束條件

那么就認為該T0值為所要識別的Te。其中φmc(ni，T0)是為了進行識別而轉化到(－π，π] 區間的相位譜主值，轉化方式如下:

得到Te后，令所有諧波的初相位為零，從該時刻出發，利用公式(7)的相位演化速度，可以重建經過時間Te后的真實時程起點演化相位譜Φ(n，Te):

進而，結合實測的Fourier幅值譜，利用逆Fourier變換即可得到重建后的時程。

理論上，重建時程應當和原時程一致，但是由于采用了松弛識別準則，重建時程和原時程在細節上會有一定出入，但基本信息保持一致。圖4為一段重建時程和原時程的比較。可以看到，重建時程與原始時程的波形基本一致，但在高頻部分有一定差別。兩段時程的相關系數為0.81(完全相同時為1)，可以認為基本達到了重建的目標。

圖4 脈動風速重建時程和原始時程的比較Fig.4 Comparison of rebuilt time history and measured time history of fluctuating wind speed

2.3 零點演化時間的統計建模

由上節可知，任意一個風速時程樣本都可以看作是由一簇諧波自零相位出發經過時間Te演化而來。對于脈動風速隨機過程，Te顯然是一個隨機變量。若能夠統計得到Te的概率分布，結合Fourier幅值譜，便可以得到模型脈動風速時程的樣本集合，從而，可以用隨機Fourier函數反映脈動風速過程。

2006年，本研究小組于江蘇某地建立了國內第一座風場觀測臺陣，經近5年觀測，得到了大量的風速數據［15，16]。采用這一臺陣實測的 10 m、20 m、28 m 和 43 m高度處各200條10分鐘風速時程，根據上節的識別方法，識別了800條樣本脈動風速時程的零點演化時間Te。對800個Te值做統計并用不同概率統計模型進行擬合，發現Gamma分布可以很好的對Te的分布進行描述，如圖5(a)所示。

圖5 相位演化時間Te的概率分布擬合Fig.5 Probability distribution fitting of evolution time Te

Gamma分布是一個兩參數連續分布函數，一般用來模擬隨機事件的等候時間，其概率密度函數為［17]:

其中，θ為尺度參數、k為形狀參數，兩者都必須為正;Γ為Gamma函數。

應用實測數據做Gamma分布的參數估計時，形狀參數 k可用下式估算［17]:

其中:

即s為實測變量均值的對數值與實測變量對數值的均值之差。Gamma分布的均值為kθ，因此容易得到尺度參數θ的計算公式:

根據前述800個樣本統計，得到的模型參數值為:k≈1.1，θ≈8.2 ×108。進而，應用其他1 200 條 10 分鐘風速記錄對該分布參數進行了檢驗，所得效果仍然非常理想，見圖5(b)。

3 脈動風速模擬

3.1 隨機Fourier函數

前已指出:譜表現法是基于功率譜密度的脈動風速模擬方法。雖然許多學者提出了不同的脈動風速功率譜密度模型［18]，但從本質上來說，由于功率譜密度是二階矩，譜表現方法僅適用于平穩隨機過程［7]。針對譜表現方法的固有局限性，我們提出應采用隨機函數模型進行風速時程模擬的觀點［9]。隨機Fourier函數的基本表達式如下式:

式中，η和γ為影響脈動風速隨機過程的可測物理因素。γ為確定性變量，η為隨機變量。當需要計算具體樣本時，可取η的具體實現值。

前已述及，Fourier譜可以分解為幅值譜和相位譜，其中Fourier幅值譜是對脈動風速能量的一種譜分解，無論對于平穩隨機過程或是非平穩隨機過程都適用［7]。根據大氣湍流的物理圖景和經典湍流理論，我們已經建立了隨機Fourier幅值譜的雙線性模型［9]，其中基本的隨機變量為地面粗糙度z0和10分鐘平均風速。本文的脈動風速時程模擬中的幅值譜即基于該雙線性模型，模型的具體細節請參閱文獻［9] 。

考慮到本文重點在于Fourier相位譜，因此進行脈動風速模擬時可給定幅值譜，如圖6所示。其中，在剪切含能區滿足“－1”冪次規律，慣性子區滿足“－5/3”冪次規律，兩個子區交點在頻率0.12 Hz。平均風速取10 m/s、湍流強度取 0.15。

圖6 Fourier幅值譜雙線性模型Fig.6 Bilinear model of Fourier amplitude spectrum

3.2 演化相位譜檢驗

由于零點演化時間Te的概率模型是由800個時程樣本統計得到，為進行集合樣本檢驗，也等概率選取800個Te樣本值，分別生成演化相位譜。這里等概率取樣的含義是:等間距分割概率分布函數，并取分割位置概率分布函數所對應的樣本值。圖7為Gamma分布的概率分布函數與等概率取值示意。

圖7 Gamma分布概率分布函數與等概率取值示意Fig.7 CDF of Gamma distribution and schematic diagram of equally probability sampling

對所選取的800個Te的樣本值，應用式(14)生成演化相位譜并換算到［0)的主值區間。圖8比較了不同頻率位置實測的800條樣本相位譜和所生成的800條演化相位譜主值的概率分布直方圖，圖中直線為均勻分布的概率密度函數。可以看到，在各個頻率點處，實測和演化所得相位譜主值均基本符合［0，2π)內均勻分布。這與已有相關研究結果［10]是相符的。

3.3 脈動風速模擬

圖8 不同頻率處相位譜主值概率分布直方圖比較Fig.8 Histogram of the distribution of phase-spectrum main value at different frequency

其中，Re表示取實部，Fs為采樣頻率，d n是頻率間隔，與采樣持時T互為倒數。除以是將能量調整至與原單邊譜相同。

圖9為任意選取的四個Te值計算出的演化相位譜主值及其分布，圖10為用其生成的脈動風速時程。兩圖從上到下，Te均依次取值為 1.131e8 s、4.036e8 s、7.129e8 s、12.997e8 s。可見，不僅相位譜主值具有［0，2π)內均勻分布的特征，所生成的脈動風速時程也與常見的風速觀測結果有很好的相似性。

圖9 不同Te所生成的相位譜主值及其分布Fig.9 Main value of phase spectrum generated by different Te and its distribution

3.4 模擬風速的檢驗

一般認為，良態風主風向脈動速度的概率分布可用正態分布近似描述［19]。比較兩個概率密度函數相近程度的一個有效指標是相對熵［20]。實際風速樣本的概率密度p(u)和其對應正態分布概率密度函數f(u)的相對熵為:

圖10 不同Te所生成脈動風速時程Fig.10 Fluctuating wind speed history generated by different Te

實際風速樣本的概率密度p(u)用直方圖代替。由于每個時間截口脈動速度u的均值為零，因此正態分布的概率密度函數為:

其中σ為截口樣本的標準差。當兩個概率密度函數完全相同時，相對熵為零。當兩個概率密度函數相差越大時，相對熵絕對值越大。

計算得到的相對熵如圖11所示。可以看到，相對熵的值較小但不為零，說明各時間截口的脈動風速分布近似符合正態分布但又略有區別，這也與已有實測結果相吻合［19]。由此可以判斷，本文提出的脈動風速模擬方法是合理的。

圖11 脈動風速樣本概率密度與正態分布的相對熵Fig.11 Relative entropy of sampled fluctuating wind speed and normal distribution

4 結論

依據分析，本文提出任何一段脈動風速時程都可以認為是由一簇具有零相位的諧波經過時間Te演化而來，根據800條實測風速記錄統計得到了Te的概率模型。這一概率模型經1 200條實測風速記錄檢驗正確。結合Fourier幅值譜，由Te的樣本值可以得到相位譜的樣本，進而進行脈動風速模擬。

與傳統方法相比，相位零點演化時間Te的概念將原本無規律的各諧波相位聯系在一起，且具有物理意義。演化相位譜的變量僅為Te，不僅大幅度減少了變量數目，也使得采用精細化的概率密度分析方法進行結構抗風可靠度的計算成為可能。

［1] Kareem A. Numerical simulation of wind effects:A probabilistic perspective［J] .Journal of Wind Engineering and Industrial Aerodynamics，2008，96:1472－1497.

［2] Rice S O.Mathematical analysis of random noise［C] .In Selected Papers on Noise and Stochastic Processes，edited by N.Wax.Dover.1954:133－294.

［3] Shinozuka M.Simulation of multivariate and multidimensional random process［J] .Journal of the Acoustical Society of America，1971，49(1):357－367.

［4] Shinozuka M，Jan C M.Digital simulation of random process and its application ［J] .Journal of Sound and Vibration，1972，25(1):111－128.

［5] Shinozuka M，Deodatis G.Simulation of stochastic processes and fields［J] .Probabilistic Engineering Mechanics，1997，12(4):203－207.

［6] Shinozuka M，Yun C B，Seya H.Stochastic methods in wind engineering［J] .Journal of Wind Engineering and Industrial Aerodynamics，1990，36:829 －843.

［7] Li J，Chen JB.Stochastic dynamics of structures［M] .John Wiley＆ Sons(Asia)Pte Ltd.2009.

［8] 李 杰.隨機動力系統的物理逼近［J] .中國科技論文在線，2006，1(2):95 －104.

［9] 李 杰，閻 啟.結構隨機動力激勵的物理模型——以脈動風速為例［J] .工程力學，2009，26 Sup II，175－183.

［10] Seong SH，Peterka J A.Experiments on Fourier phases for synthesis of non-Gaussian spikes in turbulence time series［J] . Journal of Wind Engineering and Industrial Aerodynamics，2001，89:421 －443.

［11] 金 星，廖振鵬.地震動相位特性的研究［J] .地震工程與工程振動，1993，13(1):7－13.

［12] Taylor G I.The spectrum of turbulence［J] .Proceedings of the Royal Society of London，Series A，1938，164:476.

［13] Bahraminasab A，Niry M D，Davoudi J，et al.Taylor's frozen-flow hypothesis in Burgers turbulence［J] .Physical Review，2008，E77.065302(R).

［14] Hinze J Q. Turbulence［M] . New York:McGraw-Hill，1975.

［15] 閻 啟，謝 強，李 杰.風場長期觀測與數據分析［J] .建筑科學與工程學報，2009，26(1):37－42.

［16] 李 杰，閻 啟，謝 強，等.臺風“韋帕”風場實測及風致輸電塔振動響應［J] .建筑科學與工程學報，2009，26(2):1－8.

［17] Choi S C，Wette R.Maximum likelihood estimation of the parameters of the Gamma distribution and their bias［J] .Technometrics，1969，11(4):683 －690.

［18] Dyrbye C，Hansen SO.Wind loads on structures［M] .New York:John Wiley＆ Sons，1997.

［19] Yim J Z，Chou C R，Huang W P.A study on the distributions of the measured fluctuating wind velocity components［J] . Atmospheric Environment，2000，34:1583－1590.

［20] Sobezyk K， Trebicki J. Maximum entropy principle in stochastic dynamics［J] . Probabilistic Engineering Mechanics，1990，5(3):102 －110.