高斯/非高斯混合隨機風壓場的模擬方法

羅俊杰，蘇 成，韓大建

（1.華南理工大學 土木與交通學院 土木工程系，廣州 510640，2.亞熱帶建筑科學國家重點實驗室，廣州 510640）

風荷載是空間屋蓋結構設計和分析的主要考慮因素。隨著結構跨度不斷增加，其剛度的幾何非線性特征也愈加明顯。傳統頻域分析方法，只適用于線性結構，無法滿足此類結構的風振響應分析、隨機分析和動力可靠度分析的要求，這就需要進行非線性時域分析。由此，對隨機風荷載時程的有效模擬是影響時域分析準確性的關鍵因素。

過去，作用于建筑結構的脈動風壓過程被視為具有空間相關性的平穩多點高斯隨機場。然而，由于屋蓋結構的外形不可能呈現完美的流線型，導致來流風在遇到結構的迎風面后產生分離現象；對于大跨度結構而言，當其順風向的結構尺度較大時，還有可能產生分離流重附著現象。經過大量實測數據的統計，Kumar和Stathopoulos等［1］提出利用脈動風壓樣本的偏斜度（Skewness，記為S）和峰態度（Kurtosis，記為K）這兩個高階統計量指標來劃分高斯區域和非高斯區域。據此標準，可認為處于屋蓋結構大部分相對較平緩區域的隨機風壓場樣本屬于高斯區；在氣流的分離和重附著區域，隨機風壓場樣本則屬于非高斯區［2-5］。由此，構成了一個空間相關的多點高斯/非高斯區域混合的隨機場。

目前，對這類隨機場的模擬方法主要有相關性變形法［6］，譜校正法［7］等。這些方法為了滿足實測樣本的統計特性（即功率譜密度函數和概率密度函數等）要求，在通過保持其中一個函數不變的前提下，不斷修正模擬樣本，使另外的一個函數逼近目標函數。由此，需要不斷生成新的隨機樣本。對于具有多個樣本點的風壓場而言，不但計算效率不高，而且有可能由于兩個函數的不相容而導致算法的不收斂現象［6］。

本文基于零記憶非線性轉化法［8］（ZMNL）理論，介紹這類隨機風壓場的具體模擬過程。針對隨機風壓場樣本服從對數正態分布和韋布爾分布的情況，推導了多點非高斯隨機過程與相應高斯隨機過程間的標準化協方差轉化函數；并且提出通過修正分解譜密度函數的方法，解決利用諧波合成法模擬多點高斯隨機過程時，功率譜密度函數矩陣在某些頻率點出現負定的問題。最后，通過算例表明本文所提出方法的可行性與準確性。

1 隨機風壓場模擬步驟

利用Kumar和Stathopoulos的劃分標準可以對隨機風壓場進行方便的劃分。但是，為保證隨機風壓場的整體相關性以及數值模擬的穩定性，本文假定整個風壓場的樣本均服從同一類的非高斯分布。即使這樣，從后述的算例表明：對常用于描述隨機風場的非高斯分布模型，如對數正態分布［2］或韋布爾分布［9］等，其偏斜度和峰態度也有可以符合Kumar等的劃分標準。要模擬這樣的隨機風壓場，其樣本需要同時滿足特定功率譜密度函數的和概率密度函數的要求，而且還要保持各樣本間的空間相關性。本節著重介紹利用ZMNL理論模擬隨機風壓場樣本的過程。

假設具有m個空間點的平穩隨機風壓場，可以視為m變量平穩隨機過程向量［2］，可將它記為（t）=，（注：本文的“風壓”是指無量綱的風壓系數）。首先對它進行標準化處理，使之成為均值為0，方差為1的隨機過程向量Z（t）={z1（t），z2（t），…，zm（t）}T。ZMNL法的模擬思路是：生成高斯隨機過程向量樣本X（t）={x1（t），x2（t），…，xm（t）}T，然后利用以下公式：

得到非高斯隨機過程向量樣本。其中，gi（·）為轉化函數。很明顯，由此所得到的非高斯隨機過程只能滿足概率密度函數的要求，不能達到目標功率譜密度函數的要求。此時，根據維納-辛欣定理，各點非高斯隨機過程構成的功率譜密度函數矩陣SNG，可通過傅里葉變換轉化為協方差函數矩陣ξ。由其定義可知，它與相應的高斯隨機過程的協方差函數矩陣ρ存在轉化關系。具體到矩陣中各元素間ξij和ρij的關系為：

其中：為書寫方便，此處用zi、xi（i=1，2，…m）代替zi（t）、xi（t）（下同）；μi、μj為非高斯樣本的均值；Φ［xi，xj，ρij］為非高斯隨機過程樣本的邊緣分布函數。公式（2）說明，只要確定了Φ［xi，xj，ρij］，便可用協方差函數為橋梁，得到各點高斯隨機過程和非高斯隨機過程間功率譜密度函數的轉化關系。具體的模擬步驟見圖1。值得說明的是，為了數值求解的穩定性，有必要將目標功率譜密度函數SNG中的自譜密度函數進行歸一化處理，同時互譜密度函數也作相應處理。由此得到標準化協方差函數矩陣ξ，其中各元素的定義域為［-1，1］。

圖1 ZMNL法模擬步驟Fig.1 Simulation procedures of the ZMNL approach

2 標準化協方差函數的確定

由于不同分布類型的ξ與ρ的轉化關系式并不一樣，在此推導服從對數正態分布和服從韋布爾分布的關系式。

2.1 服從對數正態分布模型的標準化協方差函數

假設隨機過程服從正態分布N（bi，），bi，為其均值和方差。令=ln（zi-ai），則隨機過程zi服從具有三參數的對數正態分布，其概率密度函數為：

經過詳細的推導，可得標準化互協方差函數ξij與ρij的轉化關系為：

當i=j時，得到zi的標準化自協方差函數：

ξii和ξij構成了非高斯樣本的標準化協方差函數矩陣ξ，而ρii和ρij構成高斯樣本的標準化協方差函數矩陣ρ。需要指出的是，由于隨機變量經過標準化處理，故ξii和ξij的取值范圍為［-1，1］，而從式（5）和（6）可以看出，ρii和ρij的取值范圍不一定落于［-1，1］的區間內，這就要研究參數c對ξ和ρ之間轉化關系的影響。

圖2為不同ci下對應的標準化自/互協方差函數。其中，圖 2（a）表明，當參數ci小于 0.1 時，ξii與ρii比較接近。隨著ci的增大，ξii與ρii僅僅在ρii=0和1兩點相交，其它地方都產生了不同程度的偏離。尤其在ρii=-1 時，ξii卻大于 -1，由于ρii的定義域為［-1，1］，這說明在進行標準化自協方差函數轉換前，首先要確定ξii的下限值，以保證ρii處于定義域范圍內。圖2（b）為相同ci對應不同cj時的標準化互協方差函數ξij。從中可以看出，ξij不但具有下限值，而且還具有上限值的要求。因此，當根據ξ來求取ρ時，首先要對各個元素進行上、下限檢驗，以保證ρ的合理性。

2.2 服從韋布爾分布模型的標準化協方差函數

對服從韋布爾分布隨機過程zi的模擬要相對復雜一些。假設zi的三參數概率密度函數模型為：

其中：ai為位移參數；pi＞0為形狀參數；（qi）pi＞0為尺度參數。令wi=zi-ai，并假設ri是服從瑞利分布的隨機過程，根據文獻［10］有ri和wi之間的關系式為：

經過推導，wi和wj的標準化互協方差函數ξij為［11］：

其中：Γ（·）為 Gamma 函數；2F1（·）表示廣義超幾何分布函數［12］；τij為隨機過程與的標準化互協方差協方差函數。當求得τij后，還要知道它與服從瑞利分布隨機過程ri與rj的標準化互協方差協方差函數γij的關系。經過推導，這個關系表達為：

從式（10）可看出，該式不涉及形狀參數pi與pj。經過數值計算，該式在τij=1和τij=-1附近會出現一些奇異點。為使該式在定義域τij∈［-1，1］內有意義，同時也為了減少計算量，可以用多項式擬合的方法得到其近似表達式。最后，ri用以下關系式求得：

其中：xi和yi為服從高斯分布N（0，2）的兩個獨立隨機過程。由此可得到γij和服從高斯分布隨機過程xi與xj的標準化互協方差函數ρij的關系為［13］：

至此，服從韋布爾分布的隨機樣本的標準化協方差函數矩陣ξ，經過從ξ→τ→γ→ρ的轉化過程，就可以得到服從高斯分布的隨機樣本的標準化協方差函數矩陣ρ。需要說明的是，根據式（9）、（10）和（12）可以發現，參數p是ξ和ρ之間的轉化的關鍵因素。經過研究表明，與2.1節討論的情況相似，當根據ξ來求取ρ時，也需要對各個元素進行上、下限的檢驗，以保證ρ的合理性。

3 分解譜密度函數的修正

當得到標準化協方差函數矩陣ρ后，利用維納-辛欣定理，可求得相應的功率譜密度函數矩陣SG，然后利用諧波合成法來生成高斯隨機場樣本。然而，雖然非高斯互譜密度函數矩陣SNG具有非負定的特點，在從ξ→ρ的轉化過程中，對于某些頻率點，受到參數的影響，會使矩陣ρ中各元素間的相關性發生不同程度的改變；此外，根據第2節分析可知，ρ的元素還有超出上、下限的可能。這些因素都會破壞SG的非負定性，最終導致某些頻點上的SG不能進行Cholesky分解。因此，需要對SG在這些頻率點處的值進行修正，以確保SG在整個頻率段范圍內的非負定性。

具體的修正過程用以下例子來說明。選取具有四個空間點的隨機風壓場為模擬目標，其中前三個點處于非高斯區域，第四個點處于高斯區域，前四階的統計參數見表1。假定點3與其它點的風壓樣本呈現負相關性，其它點間為正相關，而且四個點上的隨機過程樣本均服從韋布爾分布。各個點的風壓譜函數取下述形式［14］：

其中：為各點隨機風壓過程的方差；m1i、m2i為位置參數；n1i、n2i為形狀參數，具體取值見表2。隨機風壓樣本間的互相干函數定義為指數形式：

式中：Δij為節點間的距離，在此假設為12 m；參考高度的風速（z）=24 m/s；λij為指數衰減系數，取值如表2所示。

采用Cholesky分解算法按行進行分解［15］，每完成一行分解就進行相應的修正，然后以修正后的結果為基礎，進行下一行的分解。假定前三行的分解譜密度函數已經過修正，為H'jk（j=1，2，3；k=1，…，j），然后得到未修正分解譜密度函數（i=1，…，4），即圖 4 中的UGH43、UGH44。其中，SG在某些低頻率點由于非線性的轉化關系而導致出現負定，在圖中以零值表示。需要說明的是，即使隨機過程樣本服從對數正態分布，其分解譜密度函數曲線同樣會出現上述情況，在此不再贅述。以下對第4行的元素進行修正，具體步驟為：

表1 樣本前四階統計參數目標值及模擬值Tab.1 The first forth order of statistical parameters of target and simulating samples

表2 節點風壓譜模型及相干函數模型參數Tab.2 Parameters of the power spectral function and coherence function

（2）第2列至第3列元素，要考慮到已經修正的分解譜函數的影響，故表達為：

（3）平移該列曲線。如果M=0，則該列曲線不用平移；如果M＞0，則將連續的曲線整體前移M位，覆蓋原來的奇異值，得到新的函數曲線。這時，的尾端會出現M個空置值，可以根據曲線的發展趨勢，通過插值的方法補充完整。上述做法的目的在于保持函數曲線的連續性和平順。同時，可以重新求得前三列的互譜密度函數：

（4）第4列元素按

（5）修正第4列元素。由于上述平移過程的實質是對該節點的自譜密度函數進行的一次校正。這時求得：

同時可求得：

為便于計算機程序的編制，將上述流程如圖3所示。

圖3 分解譜密度函數修正的流程圖Fig.3 Flow chart of the modification of the factorized power spectral density function

圖4 節點4的分解自/互譜密度函數比較Fig.4 Comparisons of the factorized self-/cross-power spectral density functions of Point 4

圖5 修正后節點4自/互功率譜密度函數Fig.5 The self-/cross-power spectral density functions of Point 4 after modification

利用上述步驟，可以得到第四行的修正分解譜密度函數曲線，如圖4的MGH43、MGH44所示。從中可以看出，該方法不但修正了原曲線的奇異點，函數曲線也保持了相對完整和平順。圖5為利用修正后的分解譜密度函數，求得的譜密度函數曲線MPSDG4i（i=1，…，4）?？梢钥闯?，函數曲線的形態得以很好地保持；而且自譜密度函數值比互譜密度函數值都要大，這保證了Cholesky分解的順利進行。值得說明的是，上述方法從物理意義上講，是為了調整作用在各空間點上隨機風壓過程的空間相關關系，以保持非高斯風壓場和由其轉化得到的高斯風壓場的整體空間相關性。在修正了整個風壓場的協方差函數矩陣后，利用諧波合成法可以獲得整個高斯隨機風壓場樣本X={x1，x2，x3，x4}T。再根據式（4）或式（8）與式（11），轉化成相應的服從對數正態分布或韋布爾分布的非高斯隨機風壓過程樣本。

4 數值算例

經過上述三方面的驗證說明，所模擬的非高斯隨機風壓時程既能滿足功率譜密度函數的要求，又能符合對概率密度函數的要求，同時還保持了隨機風壓場整體的相關性。因此，本文提出方法適用于空間多點相關的隨機風壓場模擬。

需要指出的是，在實際結構上需要模擬隨機風壓過程的點有成百上千個，模擬所耗費的時間主要取決于高斯樣本的標準化協方差函數矩陣的生成和利用諧波合成法生成高斯樣本這兩個方面。對于前者，從圖2可以看出，ξ隨ρ的增加而單調增加，這對服從對數正態分布和服從韋布爾分布模型的情況均成立。所以，對于相同參數下ξ與ρ的轉化過程，可以先求取若干個ρ值，再利用插值方法求取其它ρ值，這樣可以大大減少求取ρ值的時間。對于后者，可以采用文獻［15-16］中的方法來提高高斯樣本的生成效率，在此不在贅述。

5 結論

作用于屋蓋結構上的隨機風壓場由呈現高斯和非高斯統計特性的區域共同構成。采用非線性無記憶轉換法可以一次生成所需要的整體隨機場樣本，適合于多點非高斯隨機過程的模擬。本文詳細介紹了此法的模擬過程，并解決了其中兩個關鍵問題：

（1）以對數正態分布和韋布爾分布模型為例，推導了非高斯隨機過程轉化為高斯隨機過程的標準化協方差函數矩陣的轉化公式，并指出高斯隨機過程的標準化協方差函數矩陣ρ中，各元素均要進行上、下限的檢驗。

（2）提出了分解譜密度函數修正法，保證了模擬高斯隨機過程時，功率譜密度函數矩陣的非負定性。經過算例驗證說明，利用本文方法所模擬的多點相關隨機風壓場樣本，能同時滿足特定功率譜密度函數和目標概率密度函數的要求。

［1］Kumar K S，Stathopoulos T.Wind loads on low building roofs：a stochastic perspective［J］.Journal of Structural Engineering，ASCE2000，126（8）：944-956.

［2］Gioffrè M，Gusella V，Grigoriu M.Non-Gaussian wind pressure on prismatic buildings I：stochastic field［J］.Journal of Structural Engineering，ASCE2001，127（9）：981-989.

［3］孫 瑛.大跨屋蓋結構風荷載特性研究［D］.哈爾濱：哈爾濱工業大學，2007.

［4］葉繼紅，侯信真.大跨屋蓋脈動風壓的非高斯特性研究［J］.振動與沖擊，2010，29（7）：9-15.

［5］李壽科，李壽英，陳政清，等.大跨開合式屋蓋峰值風壓的試驗研究［J］.振動與沖擊，2010，29（11）：66-72.

［6］Deodatis G，Micaletti R C.Simulation of highly skewed non-Gaussian stochastic processes［J］.Journal of Engineering Mechanics，ASCE2001，127（12）：1284-1295.

［7］ Gurley K R，Kareem A.A conditional simulation of nonnormal velocity/pressure fields［J］.J.Wind.Eng.Ind.Aerodyn.，1998，77-78：39-51.

［8］Grigoriu M.Applied non-Gaussian processes：example，theory，simulation，linear vibration and MATLAB solutions［M］.Englewood Cliffs，N.J.：Prentice Hall，1995.

［9］黃 鵬，顧 明.懸臂弧形挑篷脈動風壓的概率特征［J］.同濟大學學報（自然科學版），2007，35（9）：1153-1157.

［10］ Szajnowski W J.The generation of correlated weibull clutter for signal detection problems［J］.IEEE Transactions on Aerospace and Electronic Systems，1977，13（5）：536-540.

［11］羅俊杰.大跨度屋蓋風荷載模擬方法研究及張拉式索-膜結構風振響應分析［D］.廣州：華南理工大學，2008.

［12］ Abramowitz M，Stegun I A.Handbook of mathematical functions［M］.New York：Dover，1972.

［13］ Ertel R B，Reed J H.Generation of two equal power correlated Rayleigh fading envelopes［J］. IEEE Communications Letters，1998，2（10）：276-278.

［14］ Kumar K S，Stathopoulos T.Power spectra of wind pressures on low building roofs［J］.J.Wind.Eng.Ind.Aerodyn.，1998，74-76：665-674.

［15］羅俊杰，韓大建.諧波合成法模擬隨機風場的優化算法［J］.華南理工大學學報，2007，35（7）：105-109.

［16］羅俊杰，韓大建.大跨度結構隨機脈動風場的快速模擬方法［J］.工程力學，2008，25（5）：96-101.