基于譜修正法的多維隨機(jī)載荷時(shí)域模擬

豆碩 劉志明 毛立勇

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

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

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

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

1 多維隨機(jī)載荷譜表示方法

譜表示法是隨機(jī)載荷模擬中應(yīng)用最為廣泛的方法，也是載荷時(shí)域模擬的基礎(chǔ)算法?？紤]多點(diǎn)相干性的多維隨機(jī)載荷表示為[15]

(1)

(j=1,2,…,m)

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

(2)

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

(3)

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

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

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

(4)

(5)

(6)

(7)

由功率譜和協(xié)方差的關(guān)系可得

(8)

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

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

(9)

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

2 基于譜修正法的多維隨機(jī)載荷時(shí)域模擬

2.1 相位角對(duì)載荷模擬精度的影響

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

對(duì)于第1點(diǎn)載荷，在頻率點(diǎn)ωl處功率譜僅分解為1個(gè)余弦函數(shù)，不存在相干成分，因此合成時(shí)域載荷的功率譜等于目標(biāo)值。而其他點(diǎn)載荷由多個(gè)余弦函數(shù)的線(xiàn)性疊加而成，存在相干成分，使得合成時(shí)域載荷的功率譜相較于目標(biāo)值有較大誤差。

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

(10)

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

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

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

(11)

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

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

(12)

(n=1,2,…,N)

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

(13)

(n=1,2,…,N)

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

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

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

2.2 譜修正的多維隨機(jī)載荷時(shí)域模擬算法

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

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

圖1 基于譜修正法的多維隨機(jī)載荷時(shí)域模擬流程圖

首先根據(jù)目標(biāo)功率譜矩陣(CSDMT)初始化首次迭代合成載荷的功率譜矩陣CSDM0，并計(jì)算CSDM0對(duì)應(yīng)的Cholesky分解矩陣H0(ω):

CSDM0=CSDMT

(16)

(17)

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

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

修正合成時(shí)域載荷的Cholesky分解矩陣，這里沒(méi)有直接對(duì)功率譜矩陣進(jìn)行修正，是因?yàn)楣β首V之間是耦合的，例如S21(ω)是第1點(diǎn)和第2點(diǎn)載荷的互功率譜，不能相對(duì)于S11(ω)和S22(ω)進(jìn)行獨(dú)立的修正。另外修正后的功率譜矩陣可能是非正定的，導(dǎo)致功率譜矩陣不能進(jìn)行Cholesky分解。為保證每次迭代的功率譜矩陣為正定矩陣，直接修正上一步功率譜矩陣對(duì)應(yīng)的Cholesky分解矩陣，采用式(18)進(jìn)行修正:

(18)

計(jì)算合成時(shí)域載荷的功率譜矩陣CSDMk與目標(biāo)功率譜矩陣CSDMT的誤差值，如式(19)所示。如果滿(mǎn)足誤差要求，則迭代結(jié)束，否則重新進(jìn)行迭代:

(19)

3 數(shù)值分析

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

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

(20)

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

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

(21)

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

圖2 地震動(dòng)和波浪波動(dòng)功率譜密度

下面對(duì)一個(gè)三維隨機(jī)載荷模型進(jìn)行數(shù)值分析，以證明所提出算法的精確性。首先定義一個(gè)三維功率譜矩陣CSDM3×3，然后通過(guò)譜表示法和本研究提出的譜修正算法，生成對(duì)應(yīng)的三維時(shí)域載荷，并與目標(biāo)功率譜矩陣進(jìn)行對(duì)比分析。

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

(22)

式中：j表示第j點(diǎn)隨機(jī)載荷；ω為圓頻率；βj為功率譜強(qiáng)度系數(shù)，第1點(diǎn)載荷β1=1，第2點(diǎn)載荷β2=0.7，第3點(diǎn)載荷β3=0.5。

互功率譜采用自功率譜和相干函數(shù)結(jié)合的表達(dá)形式：

(23)

式中，Si(ω)、Sj(ω)為第i、j兩點(diǎn)的自功率譜；γij是對(duì)應(yīng)的相干函數(shù)。

相干函數(shù)采用如下形式

(24)

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

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

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

(25)

式中，SYM表示矩陣對(duì)稱(chēng)元素。

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

(26)

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

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

圖4 模擬時(shí)域載荷的互功率譜密度和相干函數(shù)對(duì)比

圖5 模擬時(shí)域載荷的自功率譜密度對(duì)比

譜修正法合成時(shí)域載荷的功率譜與目標(biāo)功率譜的相對(duì)誤差隨迭代次數(shù)的變化曲線(xiàn)如圖6所示。結(jié)果顯示，該算法的效率較高，在迭代的初始階段，功率譜的相對(duì)誤差值隨迭代次數(shù)的增加迅速減小，當(dāng)?shù)螖?shù)達(dá)到9次時(shí)，各點(diǎn)自功率譜和互功率譜的相對(duì)誤差均小于5%。當(dāng)相對(duì)誤差達(dá)到一定精度以后，進(jìn)入了緩慢減小階段，迭代次數(shù)達(dá)到30次時(shí)，相對(duì)誤差值如式(24)所示，均在1.5%以下。

圖6 譜修正法的功率譜密度相對(duì)誤差隨迭代次數(shù)的變化曲線(xiàn)

4 結(jié)論

(1)在嚴(yán)格的理論公式推導(dǎo)基礎(chǔ)上，通過(guò)對(duì)譜表示法進(jìn)行誤差分析，指出與功率譜相關(guān)的隨機(jī)相位角是合成時(shí)域載荷功率譜矩陣與目標(biāo)功率譜矩陣有較大偏差的原因，為提高多維隨機(jī)載荷時(shí)域模擬精度提供了理論依據(jù)。

(2)為消除相位角在載荷時(shí)域模擬中的影響，以提高多維隨機(jī)載荷時(shí)域模擬精度，提出了基于譜修正法的多維隨機(jī)載荷時(shí)域模擬算法。利用合成時(shí)域載荷功率譜矩陣與目標(biāo)功率譜矩陣的偏差，引入誤差反饋?lái)?xiàng)，并采用迭代技術(shù)修正合成時(shí)域載荷的功率譜矩陣，使合成時(shí)域載荷的功率譜矩陣逐漸趨近于目標(biāo)功率譜矩陣，進(jìn)而提高多維隨機(jī)載荷時(shí)域模擬的精度。

(3)通過(guò)對(duì)一個(gè)三維隨機(jī)載荷模型進(jìn)行數(shù)值分析，驗(yàn)證了所提出算法的精確性。在譜表示法合成的時(shí)域載荷中，第1點(diǎn)載荷自功率譜的相對(duì)誤差為6.64%，其他各點(diǎn)載荷的自功率譜以及互功率譜的相對(duì)誤差較大，最大值達(dá)到36.19%。而通過(guò)本研究提出的譜修正法合成的時(shí)域載荷的功率譜，相對(duì)誤差均在1.50%以下。另外該算法的效率較高，當(dāng)?shù)螖?shù)達(dá)到9次時(shí)，各點(diǎn)載荷自功率譜以及互功率譜相對(duì)誤差迅速減小到5%以下。

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