豆碩 劉志明 毛立勇

(北京交通大學 機械與電子控制工程學院，北京 100044)

工程結構經常受到多維隨機載荷的作用，如大跨度橋梁受到的多點地震激勵，高層建筑受到的多點風壓場，以及車輛在不平順路面行駛過程中產生的輪心振動[1- 4]。結構在這些多維隨機載荷的作用下容易發(fā)生失效，為保證結構的安全性，在結構評估中需要精確的多維隨機載荷時間歷程。由于載荷的隨機性，人們通常采用功率譜矩陣描述這些多維隨機載荷，如何從功率譜矩陣反演出精確的多維隨機載荷時間歷程對于結構的安全評估至關重要。

研究人員提出了多種多維隨機載荷時域模擬方法。Smallwood等[5]利用Cholesky分解方法將功率譜矩陣分解為頻響函數和自功率譜函數，然后在頻域中生成一幀隨機數據，利用傅里葉逆變換技術將數據轉換到時域中，將連續(xù)的數據幀加窗、重疊并相加，以形成任意長度的時域數據，但需要大量計算成本。為提高大規(guī)模隨機載荷時域模擬速度，學者們進行了相關研究。Li等[6]在多維隨機載荷時域模擬中，基于樣條插值快速傅里葉變換技術減少功率譜矩陣的Cholesky分解次數，隨后采用快速傅里葉算法進一步提高了載荷時域模擬的計算速度。Gao等[7]在譜表示法基礎上，提出了一種簡化的遍歷空間相關載荷模擬方法，給出了功率譜矩陣的Cholesky分解下三角矩陣元素的顯式解，避免了互功率譜矩陣在每個頻率點的重復Cholesky分解，極大地提高了計算效率。針對功率譜矩陣由于奇異性而導致矩陣分解受到限制的問題。Liu等[8]基于Fourier-Stieltjes積分公式，提出了降維譜表示法和降維本征正交分解的多維隨機載荷時域模擬算法，能夠用少量的基本隨機變量表示多維平穩(wěn)隨機過程。Peng等[9]提出了一種基于時空隨機場和基于本征正交分解插值的混合多維時域載荷模擬算法。該方法消除了傳統(tǒng)的譜矩陣分解要求，應用二維快速傅里葉變換技術極大地提高了模擬效率。但是模擬樣本數據的近似性不高，相關函數也不具有遍歷性。Hao等[10]在研究空間相關多點地震動時提出了一種多維隨機載荷時域模擬方法，稱為HOP法。該方法首先生成一點隨機載荷，在生成后面載荷時考慮已生成載荷的相關性。學者們以HOP法為基礎也提出了一些改進算法。屈鐵軍等[11- 12]提出了生成任一點隨機載荷時均考慮其他點載荷相關性的HOP改進算法，減小了合成各點載荷的幅值差異。董汝博等[13]將合成各點載荷的隨機相位角進行統(tǒng)一，相位角只與頻率分量有關，保證了隨機載荷的局部收斂性。趙博等[14]對HOP法、屈鐵軍和董汝博提出的方法進行了驗證和分析，表明HOP法生成各點載荷之間的相干性與目標值具有較好的吻合度。而其他的改進算法，由于多余的附加項，導致各點載荷之間的相干性紊亂，使合成載荷的相干性無法匹配目標值?？偨Y可得，以上所述的多維隨機載荷時域模擬方法主要以單向模擬為主，即從功率譜矩陣一次合成所需的時域載荷，合成過程中誤差較大，并且難以控制所需的精度。如果能利用合成時域載荷功率譜矩陣與目標功率譜矩陣的偏差，引入誤差反饋項，采用迭代技術修正合成時域載荷的功率譜矩陣，將能極大地提高合成時域載荷的模擬精度。

本研究首先介紹了多維隨機載荷譜表示法的基本理論，從功率譜矩陣得到多維隨機時域載荷；其次，詳細分析了合成時域載荷功率譜矩陣與目標功率譜矩陣存在較大偏差的原因，為多維隨機載荷時域模擬修正提供了理論依據，進而提出了基于譜修正法的多維隨機載荷時域模擬算法；最后，對一個三維隨機載荷模型進行數值分析，證明了所提出算法的精確性。

1 多維隨機載荷譜表示方法

譜表示法是隨機載荷模擬中應用最為廣泛的方法，也是載荷時域模擬的基礎算法?？紤]多點相干性的多維隨機載荷表示為[15]

(1)

(j=1,2,…,m)

對于一個多維平穩(wěn)高斯隨機過程Y(t)=[y1(t)y2(t) …ym(t)]T，其功率譜矩陣為

(2)

式中，Sij(ω)為隨機過程yi(t)和yj(t)之間的互功率譜密度函數，即

(3)

式中:ui(t)和uj(t)為第i點和第j點載荷時間歷程；T為載荷持續(xù)時間；e為自然常數；當i=j時，Sii(ω)表示自功率譜密度函數。

由于自功率譜Sii(ω)是實數，而互功率譜Sij(ω)為復數，具有如下性質：Reij(ω)=Reji(ω)，Imij(ω)=-Imji(ω)。因此對于一個多維隨機過程Y(t)，其功率譜矩陣S(ω)是一個正定的Hermitian矩陣。

對功率譜矩陣S(ω)進行Cholesky分解:

(4)

(5)

(6)

(7)

由功率譜和協方差的關系可得

(8)

將式(8)帶入式(1)中，第j點隨機載荷可以表示為

θjk(ωl)+ψkl],j=1,2,…,m

(9)

式中:Hjk(ωl)是功率譜矩陣S(ω)的Cholesky分解矩陣元素；θjk(ωl)是Hjk(ωl)對應的相位角，

2 基于譜修正法的多維隨機載荷時域模擬

2.1 相位角對載荷模擬精度的影響

對于載荷之間相干的情況，式(9)譜表示法合成的多維時域載荷，除了第1點載荷的自功率譜外，其他載荷的自功率譜以及互功率譜與目標值均有較大誤差。下面分析出現這種誤差的原因。在式(9)中，功率譜矩陣分解為多個余弦函數的幅值，進行余弦函數的線性疊加時，各窄帶信號的隨機相位角保持不變，使得生成隨機信號的功率譜與相位角相關，進而導致合成載荷的功率譜圍繞真實值波動。

對于第1點載荷，在頻率點ωl處功率譜僅分解為1個余弦函數，不存在相干成分，因此合成時域載荷的功率譜等于目標值。而其他點載荷由多個余弦函數的線性疊加而成，存在相干成分，使得合成時域載荷的功率譜相較于目標值有較大誤差。

計算隨機載荷ui(t)和uj(t)的內積〈ui(t)uj(t)〉，可以得到合成時域載荷的功率譜與目標功率譜之間的誤差關系。由帕斯瓦爾定理[16]可得，周期信號可以等效為各次諧波的疊加，因此傅里葉系數的平方和與原信號的功率相等。式(10)是連續(xù)傅里葉變換形式下的帕斯瓦爾定理，

(10)

式中，*表示復共軛，Xi(f)和Xj(f)分別是ui(t)和uj(t)的傅里葉譜。

由帕斯瓦爾定理可得ui(t)和uj(t)時域平均值與功率譜之間的關系為

E(〈uin(t)ujn(t)〉)=Sij(ωn)Δω

(11)

ui(t)和uj(t)采用式(1)進行表示，帶入式(11)中可得

E(〈uin(t)ujn(t)〉)=

(12)

(n=1,2,…,N)

E(〈uin(t)ujn(t)〉)=

(13)

(n=1,2,…,N)

式中，Re表示實部，lim(ωn)=Aim(ωn)eiθim(ωn)，*表示復共軛。

隨機載荷ui(t)和uj(t)的時域平均值采用功率譜表示為

E(〈uin(t)ujn(t)〉)=Sij(ωn)Δω=

2ΔωRe(Sij(ωn)),n=1,2,…,N

(14)

將式(7)和(8)帶入式(13)中得

E(〈uin(t)ujn(t)〉)=

2ΔωRe(Sij(ωn)),n=1,2,…,N

(15)

由式(14)和(15)可得，通過譜表示法合成的各點時域載荷的自功率譜及互功率譜的期望數值等于目標值。所以，當生成大量樣本載荷時，合成時域載荷的平均功率譜能匹配目標功率譜。而式(13)成立的前提條件是兩個隨機變量的合成,對于一次合成樣本載荷而言，式(12)中隨機相位角ψmn和ψrn是與時間無關的，縮減為了兩個樣本的合成，所以當m>j和m≠r時，〈uin(t)ujn(t)〉≠0，進而導致式(14)和(15)不再相等。因此，由于相位角的影響，除了合成第1點載荷的自功率譜，其他各點載荷的自功率譜及互功率譜會圍繞真實值波動，始終存在偏差。

2.2 譜修正的多維隨機載荷時域模擬算法

上節(jié)分析表明，與功率譜相關的隨機相位角是合成時域載荷的功率譜矩陣與目標功率譜矩陣偏差較大的原因。這也為提高多維隨機載荷時域模擬精度提供了思路。如果能利用合成時域載荷功率譜矩陣與目標功率譜矩陣的偏差，引入誤差反饋項，并采用迭代技術修正合成時域載荷的功率譜矩陣，使合成時域載荷的功率譜矩陣逐漸趨近于目標功率譜矩陣，將會極大地提高多維隨機載荷時域模擬的精度。

為消除相位角在載荷模擬中的影響，本研究提出了一種基于譜修正的多維隨機載荷時域模擬算法，算法流程如圖1所示，下面進行詳細地介紹。

圖1 基于譜修正法的多維隨機載荷時域模擬流程圖

首先根據目標功率譜矩陣(CSDMT)初始化首次迭代合成載荷的功率譜矩陣CSDM0，并計算CSDM0對應的Cholesky分解矩陣H0(ω):

CSDM0=CSDMT

(16)

(17)

生成m×N維的隨機相位角ψmN，需要注意的是在后續(xù)的迭代過程中相位角ψmN保持不變，使相位角在載荷模擬中的影響保持不變，進而保證每次迭代合成時域載荷之間具有相同的相位關系。

通過式(9)所示的多維隨機載荷譜表示方法，得到功率譜矩陣CSDMk對應的多維時域載荷yj(t),j=1,2,…,m。并計算yj(t)的功率譜矩陣CSDMk，及其對應的Cholesky分解矩陣Hk(ω)。

修正合成時域載荷的Cholesky分解矩陣，這里沒有直接對功率譜矩陣進行修正，是因為功率譜之間是耦合的，例如S21(ω)是第1點和第2點載荷的互功率譜，不能相對于S11(ω)和S22(ω)進行獨立的修正。另外修正后的功率譜矩陣可能是非正定的，導致功率譜矩陣不能進行Cholesky分解。為保證每次迭代的功率譜矩陣為正定矩陣，直接修正上一步功率譜矩陣對應的Cholesky分解矩陣，采用式(18)進行修正:

(18)

計算合成時域載荷的功率譜矩陣CSDMk與目標功率譜矩陣CSDMT的誤差值，如式(19)所示。如果滿足誤差要求，則迭代結束，否則重新進行迭代:

(19)

3 數值分析

一般采用功率譜描述隨機載荷，而目標功率譜的獲取通常采用如下方法：①如果載荷具有普適性特征，可以總結出經驗表達式，并且在相應領域形成行業(yè)規(guī)范，進而直接獲取載荷功率譜。例如在地震領域廣泛采用的Clough-Penzien地面加速度功率譜模型[17]和在海洋工程中描述波浪高度廣泛采用的P-M波高功率譜模型[18]。②對于另一些結構，其受到的載荷沒有通用性特征，例如車輛在路面行駛中產生的隨機振動，那么可以通過線路實測數據獲取目標載荷功率譜。

Clough-Penzien地面加速度功率譜模型為

(20)

式中，S0為譜強度因子，ζg和ωg分別為阻尼比和卓越頻率，ζf和ωf用于模擬地震動低頻能量。譜強度因子S0=7.123 cm2/s3，譜參數ζg=0.8，ωg=6.98 rad/s，ζf=0.8，ωf=0.1ωg，對應的地震動加速度譜如圖2(a)所示。

P-M波浪高度功率譜模型為

(21)

式中：g為重力加速度，g=9.81 m/s2；v為海面風速。海面以上19.5 m處的風速v=16.24 m/s，對應的波高功率譜如圖2(b)所示。

圖2 地震動和波浪波動功率譜密度

下面對一個三維隨機載荷模型進行數值分析，以證明所提出算法的精確性。首先定義一個三維功率譜矩陣CSDM3×3，然后通過譜表示法和本研究提出的譜修正算法，生成對應的三維時域載荷，并與目標功率譜矩陣進行對比分析。

由圖2可以看出，地震動和波浪載荷為低頻的單峰衰減功率譜模型，這也是工程領域廣泛應用的一類模型。為驗證本研究提出的算法在中高頻的寬頻帶范圍內同樣具有較高精度，本研究定義了一個帶寬為0～100 Hz的單峰衰減功率譜模型。自功率譜密度函數為

(22)

式中：j表示第j點隨機載荷；ω為圓頻率；βj為功率譜強度系數，第1點載荷β1=1，第2點載荷β2=0.7，第3點載荷β3=0.5。

互功率譜采用自功率譜和相干函數結合的表達形式：

(23)

式中，Si(ω)、Sj(ω)為第i、j兩點的自功率譜；γij是對應的相干函數。

相干函數采用如下形式

(24)

式中：i,j=1,2,3(i≠j)；aij、bij為第i、j兩點載荷的相干函數系數；第1點和第2點載荷的系數為a12=0.7，b12=0.4；第1點和第3點載荷的系數為a13=0.8，b13=0.2；第2點和第3點載荷的系數為a23=0.9，b23=0.3。

圖3 基于譜修正法的模擬載荷時間歷程

基于式(19)，譜表示法合成時域載荷的功率譜矩陣相對于目標功率譜矩陣的相對誤差表示為

(25)

式中，SYM表示矩陣對稱元素。

在式(18)優(yōu)化速度指數α=0.3，通過30次迭代，譜修正法合成時域載荷的功率譜矩陣相對于目標功率譜矩陣的相對誤差表示為

(26)

由式(25)可得譜表示法合成時域載荷的功率譜相對誤差，在第1點載荷的自功率譜的相對誤差較小，為6.64%，其他各點載荷的自功率譜和互功率譜的相對誤差較大，最大值達到36.19%。結果與第2.1小節(jié)“譜表示法誤差分析中第1點載荷功率譜的相對誤差較小，而其他各點載荷的功率譜相對誤差較大”的結論相對應。式(26)顯示，本研究提出的譜修正法合成時域載荷的功率譜相對誤差較小，第1點載荷自功率譜的相對誤差僅為0.59%，其他各點載荷的自功率譜和互功率譜的相對誤差也在1.50%以下，表明合成時域載荷的功率譜矩陣能夠較好地匹配目標功率譜矩陣。

譜表示法以及譜修正法合成時域載荷的互功率譜和相干函數與目標值的對比結果如圖4所示，合成時域載荷的自功率譜與目標值對比結果如圖5所示。圖中互功率譜12和相干函數12代表第1點和第2點載荷之間的互功率譜和相干函數，自功率譜1代表第1點載荷的自功率譜?？梢钥闯?，通過譜表示法合成時域載荷的功率譜圍繞目標功率譜波動，這與第2.1小節(jié)分析相對應。當合成大量樣本載荷時，合成時域載荷的功率譜數學期望等于目標功率譜。但是對于一個合成樣本載荷而言，具有較大的誤差，功率譜幅值最大誤差可達一倍以上。而本研究提出的譜修正法載荷模擬算法由于引入了功率譜誤差反饋項，在迭代的過程中逐步消除了功率譜誤差，使合成時域載荷的功率譜與目標功率譜高度匹配，幾乎完全趨近于目標功率譜。

圖4 模擬時域載荷的互功率譜密度和相干函數對比

圖5 模擬時域載荷的自功率譜密度對比

譜修正法合成時域載荷的功率譜與目標功率譜的相對誤差隨迭代次數的變化曲線如圖6所示。結果顯示，該算法的效率較高，在迭代的初始階段，功率譜的相對誤差值隨迭代次數的增加迅速減小，當迭代次數達到9次時，各點自功率譜和互功率譜的相對誤差均小于5%。當相對誤差達到一定精度以后，進入了緩慢減小階段，迭代次數達到30次時，相對誤差值如式(24)所示，均在1.5%以下。

圖6 譜修正法的功率譜密度相對誤差隨迭代次數的變化曲線

4 結論

(1)在嚴格的理論公式推導基礎上，通過對譜表示法進行誤差分析，指出與功率譜相關的隨機相位角是合成時域載荷功率譜矩陣與目標功率譜矩陣有較大偏差的原因，為提高多維隨機載荷時域模擬精度提供了理論依據。

(2)為消除相位角在載荷時域模擬中的影響，以提高多維隨機載荷時域模擬精度，提出了基于譜修正法的多維隨機載荷時域模擬算法。利用合成時域載荷功率譜矩陣與目標功率譜矩陣的偏差，引入誤差反饋項，并采用迭代技術修正合成時域載荷的功率譜矩陣，使合成時域載荷的功率譜矩陣逐漸趨近于目標功率譜矩陣，進而提高多維隨機載荷時域模擬的精度。

(3)通過對一個三維隨機載荷模型進行數值分析，驗證了所提出算法的精確性。在譜表示法合成的時域載荷中，第1點載荷自功率譜的相對誤差為6.64%，其他各點載荷的自功率譜以及互功率譜的相對誤差較大，最大值達到36.19%。而通過本研究提出的譜修正法合成的時域載荷的功率譜，相對誤差均在1.50%以下。另外該算法的效率較高，當迭代次數達到9次時，各點載荷自功率譜以及互功率譜相對誤差迅速減小到5%以下。

(4)以上分析結果表明，本研究建立的基于譜修正法的多維隨機載荷時域模擬算法具有較高的精度，可以應用于結構在多維載荷作用下的時域模擬研究。