基于谱修正法的多维随机载荷时域模拟

豆硕 刘志明 毛立勇

(北京交通大学 机械与电子控制工程学院，北京 100044)

工程结构经常受到多维随机载荷的作用，如大跨度桥梁受到的多点地震激励，高层建筑受到的多点风压场，以及车辆在不平顺路面行驶过程中产生的轮心振动[1- 4]。结构在这些多维随机载荷的作用下容易发生失效，为保证结构的安全性，在结构评估中需要精确的多维随机载荷时间历程。由于载荷的随机性，人们通常采用功率谱矩阵描述这些多维随机载荷，如何从功率谱矩阵反演出精确的多维随机载荷时间历程对于结构的安全评估至关重要。

研究人员提出了多种多维随机载荷时域模拟方法。Smallwood等[5]利用Cholesky分解方法将功率谱矩阵分解为频响函数和自功率谱函数，然后在频域中生成一帧随机数据，利用傅里叶逆变换技术将数据转换到时域中，将连续的数据帧加窗、重叠并相加，以形成任意长度的时域数据，但需要大量计算成本。为提高大规模随机载荷时域模拟速度，学者们进行了相关研究。Li等[6]在多维随机载荷时域模拟中，基于样条插值快速傅里叶变换技术减少功率谱矩阵的Cholesky分解次数，随后采用快速傅里叶算法进一步提高了载荷时域模拟的计算速度。Gao等[7]在谱表示法基础上，提出了一种简化的遍历空间相关载荷模拟方法，给出了功率谱矩阵的Cholesky分解下三角矩阵元素的显式解，避免了互功率谱矩阵在每个频率点的重复Cholesky分解，极大地提高了计算效率。针对功率谱矩阵由于奇异性而导致矩阵分解受到限制的问题。Liu等[8]基于Fourier-Stieltjes积分公式，提出了降维谱表示法和降维本征正交分解的多维随机载荷时域模拟算法，能够用少量的基本随机变量表示多维平稳随机过程。Peng等[9]提出了一种基于时空随机场和基于本征正交分解插值的混合多维时域载荷模拟算法。该方法消除了传统的谱矩阵分解要求，应用二维快速傅里叶变换技术极大地提高了模拟效率。但是模拟样本数据的近似性不高，相关函数也不具有遍历性。Hao等[10]在研究空间相关多点地震动时提出了一种多维随机载荷时域模拟方法，称为HOP法。该方法首先生成一点随机载荷，在生成后面载荷时考虑已生成载荷的相关性。学者们以HOP法为基础也提出了一些改进算法。屈铁军等[11- 12]提出了生成任一点随机载荷时均考虑其他点载荷相关性的HOP改进算法，减小了合成各点载荷的幅值差异。董汝博等[13]将合成各点载荷的随机相位角进行统一，相位角只与频率分量有关，保证了随机载荷的局部收敛性。赵博等[14]对HOP法、屈铁军和董汝博提出的方法进行了验证和分析，表明HOP法生成各点载荷之间的相干性与目标值具有较好的吻合度。而其他的改进算法，由于多余的附加项，导致各点载荷之间的相干性紊乱，使合成载荷的相干性无法匹配目标值。总结可得，以上所述的多维随机载荷时域模拟方法主要以单向模拟为主，即从功率谱矩阵一次合成所需的时域载荷，合成过程中误差较大，并且难以控制所需的精度。如果能利用合成时域载荷功率谱矩阵与目标功率谱矩阵的偏差，引入误差反馈项，采用迭代技术修正合成时域载荷的功率谱矩阵，将能极大地提高合成时域载荷的模拟精度。

本研究首先介绍了多维随机载荷谱表示法的基本理论，从功率谱矩阵得到多维随机时域载荷；其次，详细分析了合成时域载荷功率谱矩阵与目标功率谱矩阵存在较大偏差的原因，为多维随机载荷时域模拟修正提供了理论依据，进而提出了基于谱修正法的多维随机载荷时域模拟算法；最后，对一个三维随机载荷模型进行数值分析，证明了所提出算法的精确性。

1 多维随机载荷谱表示方法

谱表示法是随机载荷模拟中应用最为广泛的方法，也是载荷时域模拟的基础算法。考虑多点相干性的多维随机载荷表示为[15]

(1)

(j=1,2,…,m)

对于一个多维平稳高斯随机过程Y(t)=[y1(t)y2(t) …ym(t)]T，其功率谱矩阵为

(2)

式中，Sij(ω)为随机过程yi(t)和yj(t)之间的互功率谱密度函数，即

(3)

式中:ui(t)和uj(t)为第i点和第j点载荷时间历程；T为载荷持续时间；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)是连续傅里叶变换形式下的帕斯瓦尔定理，

(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 谱修正的多维随机载荷时域模拟算法

上节分析表明，与功率谱相关的随机相位角是合成时域载荷的功率谱矩阵与目标功率谱矩阵偏差较大的原因。这也为提高多维随机载荷时域模拟精度提供了思路。如果能利用合成时域载荷功率谱矩阵与目标功率谱矩阵的偏差，引入误差反馈项，并采用迭代技术修正合成时域载荷的功率谱矩阵，使合成时域载荷的功率谱矩阵逐渐趋近于目标功率谱矩阵，将会极大地提高多维随机载荷时域模拟的精度。

为消除相位角在载荷模拟中的影响，本研究提出了一种基于谱修正的多维随机载荷时域模拟算法，算法流程如图1所示，下面进行详细地介绍。

图1 基于谱修正法的多维随机载荷时域模拟流程图

首先根据目标功率谱矩阵(CSDMT)初始化首次迭代合成载荷的功率谱矩阵CSDM0，并计算CSDM0对应的Cholesky分解矩阵H0(ω):

CSDM0=CSDMT

(16)

(17)

生成m×N维的随机相位角ψmN，需要注意的是在后续的迭代过程中相位角ψ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 数值分析

一般采用功率谱描述随机载荷，而目标功率谱的获取通常采用如下方法：①如果载荷具有普适性特征，可以总结出经验表达式，并且在相应领域形成行业规范，进而直接获取载荷功率谱。例如在地震领域广泛采用的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)优化速度指数α=0.3，通过30次迭代，谱修正法合成时域载荷的功率谱矩阵相对于目标功率谱矩阵的相对误差表示为

(26)

由式(25)可得谱表示法合成时域载荷的功率谱相对误差，在第1点载荷的自功率谱的相对误差较小，为6.64%，其他各点载荷的自功率谱和互功率谱的相对误差较大，最大值达到36.19%。结果与第2.1小节“谱表示法误差分析中第1点载荷功率谱的相对误差较小，而其他各点载荷的功率谱相对误差较大”的结论相对应。式(26)显示，本研究提出的谱修正法合成时域载荷的功率谱相对误差较小，第1点载荷自功率谱的相对误差仅为0.59%，其他各点载荷的自功率谱和互功率谱的相对误差也在1.50%以下，表明合成时域载荷的功率谱矩阵能够较好地匹配目标功率谱矩阵。

谱表示法以及谱修正法合成时域载荷的互功率谱和相干函数与目标值的对比结果如图4所示，合成时域载荷的自功率谱与目标值对比结果如图5所示。图中互功率谱12和相干函数12代表第1点和第2点载荷之间的互功率谱和相干函数，自功率谱1代表第1点载荷的自功率谱。可以看出，通过谱表示法合成时域载荷的功率谱围绕目标功率谱波动，这与第2.1小节分析相对应。当合成大量样本载荷时，合成时域载荷的功率谱数学期望等于目标功率谱。但是对于一个合成样本载荷而言，具有较大的误差，功率谱幅值最大误差可达一倍以上。而本研究提出的谱修正法载荷模拟算法由于引入了功率谱误差反馈项，在迭代的过程中逐步消除了功率谱误差，使合成时域载荷的功率谱与目标功率谱高度匹配，几乎完全趋近于目标功率谱。

图4 模拟时域载荷的互功率谱密度和相干函数对比

图5 模拟时域载荷的自功率谱密度对比

谱修正法合成时域载荷的功率谱与目标功率谱的相对误差随迭代次数的变化曲线如图6所示。结果显示，该算法的效率较高，在迭代的初始阶段，功率谱的相对误差值随迭代次数的增加迅速减小，当迭代次数达到9次时，各点自功率谱和互功率谱的相对误差均小于5%。当相对误差达到一定精度以后，进入了缓慢减小阶段，迭代次数达到30次时，相对误差值如式(24)所示，均在1.5%以下。

图6 谱修正法的功率谱密度相对误差随迭代次数的变化曲线

4 结论

(1)在严格的理论公式推导基础上，通过对谱表示法进行误差分析，指出与功率谱相关的随机相位角是合成时域载荷功率谱矩阵与目标功率谱矩阵有较大偏差的原因，为提高多维随机载荷时域模拟精度提供了理论依据。

(2)为消除相位角在载荷时域模拟中的影响，以提高多维随机载荷时域模拟精度，提出了基于谱修正法的多维随机载荷时域模拟算法。利用合成时域载荷功率谱矩阵与目标功率谱矩阵的偏差，引入误差反馈项，并采用迭代技术修正合成时域载荷的功率谱矩阵，使合成时域载荷的功率谱矩阵逐渐趋近于目标功率谱矩阵，进而提高多维随机载荷时域模拟的精度。

(3)通过对一个三维随机载荷模型进行数值分析，验证了所提出算法的精确性。在谱表示法合成的时域载荷中，第1点载荷自功率谱的相对误差为6.64%，其他各点载荷的自功率谱以及互功率谱的相对误差较大，最大值达到36.19%。而通过本研究提出的谱修正法合成的时域载荷的功率谱，相对误差均在1.50%以下。另外该算法的效率较高，当迭代次数达到9次时，各点载荷自功率谱以及互功率谱相对误差迅速减小到5%以下。

(4)以上分析结果表明，本研究建立的基于谱修正法的多维随机载荷时域模拟算法具有较高的精度，可以应用于结构在多维载荷作用下的时域模拟研究。