引 言
结构抗震设计中，地震动、结构、结构响应是主要的研究对象．地震动模型是研究结构地震动响应的依据和基础．目前，地震动模型可以分为时程激励模型、反应谱激励模型和随机激励模型三种．时程激励模型被认为是确定性模型，是研究其他地震动模型的基础；反应谱激励模型是建立在反应谱理论上的一种地震动模型；从统计学角度上来看，地震动的发生具有显著的随机性，因此在过去的几十年内有许多学者对地震动随机激励的模型化进行研究，取得了不小的成就，有多种模型相继被提出和改进[1-2]．1947年，美国学者Housner首次将地震动看作随机过程，提出用白噪声模型来模拟地震动，开启了用随机方法研究地震动的历史[2]；20世纪50年代以后，日本学者Kanai[3]和Tajimi在白噪声模型的基础上提出了过滤白噪声激励模型，即Kanai-Tajimi模型(K-T谱)，但K-T谱假定基岩地震加速度为白噪声，不能反映基岩地震动的频谱特征．为此，很多学者在此基础上尝试使用不同形式的线性滤波器，获得改进地震动的模型，其中美国学者Clough和Penzien[4]提出了C-P谱，该激励模型为消除低频成分的影响，把基岩上覆盖土层当做2个关联的单自由度线性滤波器，改善了K-T谱过分夸大低频能量的缺点，因此，在工程上得到了广泛的应用．李创第等[5]研究了基于C-P谱的隔震结构随机地震动响应分析的复模态法，运用扩阶的解法将地震动协方差表述为2对共轭振子的线性组合，但所获结果比较复杂，且仅获得了结构的位移和速度方差．李春祥等[6]利用虚拟激励法研究了基于C-P谱的多重调谐质量阻尼器(MTMD)的地震特性，但所得响应表达式仅为功率谱，且表达式复杂，计算效率较差．
振动系统的阻尼特性及其动力响应是建筑结构动力分析领域中的一个重要问题．近年来，随着越来越多的高分子材料和高强度的聚合物被广泛应用于建筑工程中，传统的黏滞阻尼模型已经不能准确地描述结构复杂的阻尼耗能特性．而非黏滞阻尼模型将阻尼力表现为质点速度与某一核函数的卷积，这种模型比传统的黏滞阻尼更具有一般性，因此越来越多的学者将系统的阻尼建模为非黏滞阻尼来研究[7-8]．李春江[9]研究了卷积型非黏滞阻尼模型在建筑结构中的适用性,将利用振动台实验所得到的非黏滞阻尼模型参数与用电脑软件仿真计算所得到的模型参数识别结果进行对比，两者吻合度较高，表明卷积型非黏滞阻尼更适用于建筑结构．王禹[10]研究了指数型非黏滞阻尼模型参数的识别，提出了更为准确的松弛因子的迭代计算方法．段忠东等[11]用精细积分法研究了卷积型非黏滞阻尼模型系统的时程分析，所得结果精度较高，但该方法计算量大，效率低．Liu[12]利用卷积的微分性质和加速度的中心差分公式，得出了非黏滞性阻尼结构体系的动力响应的时程分析方法．以上研究均是对地震动时程激励下的非黏滞阻尼结构的动力响应分析，而实际地震动具有随机性，为此，研究非黏滞阻尼结构基于随机激励下的响应具有重要的工程意义．
结构在平稳随机激励下的地震动响应分析主要有时域法和频域法[13]．时域法又称之为模态分析法，其应用的前提是激励要有协方差．方同等[14]利用复模态法研究了结构基于地震动激励下的响应问题，用分析的形式给出了系统在激励下的稳态响应的谱矩，并将地震动激励的协方差表示成一对共轭振子的线性组合，但所得结论表达式复杂难懂，不利于应用到实际工程当中．频域法中具有代表性的方法有传递函数法和虚拟激励法[15]等．特别是林家浩教授提出的虚拟激励法已成功应用于各种振动领域[16]．然而工程中需要获得结构响应谱，还要进一步求系统的均方响应，或者可靠度分析时的谱矩等参数，这些都得通过相当繁杂的积分运算．传递函数法往往要用到复杂的、不容易理解的拉氏变换，而且往往还得求助于数值积分，对于大型复杂的结构来说，计算效率比较低．李创第等[17]研究了非平稳激励下六参数黏弹性阻尼结构的响应问题，运用虚拟激励法与精细积分形式相结合的方法得到了响应的方差，但并没有分析谱矩、可靠度这些重要性质，且计算效率比较低．
本文针对非黏滞阻尼结构卷积型的地震动方程分析复杂的问题，提出了一种简明封闭解的新方法．首先，依据非黏滞阻尼的积分型指数模型，推导出其等效的微分模型；其次利用C-P谱基于白噪声激励的滤波方程组，两者结合重构其地震动方程，将其转化为基于白噪声激励的微分表达式的地震动方程；然后，利用复模态法并结合白噪声激励的特点，获得结构响应的功率谱及0～2阶谱矩的简明封闭解；最后，分析了基于首次超越破坏准则以及Markov分布规则[18-19]的结构动力可靠度．
1 非黏滞阻尼结构地震动方程的重构
假设一个单自由度结构，质量、阻尼和刚度分别为m,c,k，其在地震动作用下的运动方程为
(1)
其中
分别为结构相对于地面的加速度、速度和位移；G(t)=μ1e-μ1t+μ2e-μ2t为阻尼核函数[9]，其中μ1,μ2为松弛因子．
1.1 非黏滞阻尼的精确等效微分形式
式(1)为积分微分混合型运动方程，不易获得简明解，为此，本文提出其精确等效的本构关系．由于核函数G(t)为2个指数函数的线型组合，则可令
(2)
对式(2)进行求导得
(3)
式(3)简化为
(4)
至此，式(1)的阻尼力项可改写为
(5)
1.2 地震动方程的重构
C-P谱把基岩上覆盖土层当作两个关联的单自由度线性滤波器，其运动方程为
(6)
其中，ξg,ωg为基岩以上场地土(第一滤波)的阻尼比和卓越频率；ξf,ωf分别为地面处场地土(第二滤波)的阻尼比和卓越频率；
分别为地面相对于基岩运动的速度和位移；
为基岩运动加速度，其为白噪声随机过程．
把式(5)、(6)代入式(1)，则式(1)可改写为
(7)
式中
(8)
r=[-1 -1 0 0 0 0 0 0]T，
式中，上标“T”表示转置．
2 非黏滞阻尼结构地震动响应的复模态解耦
根据复模态理论[13]，方程(8)存在对角阵p和右、左特征向量U,V使之解耦，即存在如下关系:
(9)
引入复模态变换:
y=Uz,
(10)
式中,z为复模态广义参数．把式(10)代入式(7)，并左乘VT，可推得
(11)
式中η=VTr/(VTMU)．再结合式(9)，则式(11)可以化简为
(12)
由于p为对角阵，则式(12)的分量形式为
(13)
其中，zi,ηi,pi分别为z，η，p的分量．
由式(12)、(13)可得各个响应的表达式，以结构的位移x和速度
为代表，可表示为
(14)
(15)
其中，ui为右特征向量的第i行向量;λn,i为结构响应的强度系数，可以表示为
λn,i=un,iηi．
(16)
3 结构响应的协方差及谱矩的简明解
3.1 结构响应的协方差
工程随机振动问题中经常涉及平稳随机过程，因而有必要进一步研究它的统计特性．二阶平稳随机过程最重要的矩特性是它的相关函数与协方差函数[13]．本文采用C-P谱激励模型，其基岩运动可以看做白噪声激励，均值为零，其相关函数等同于协方差函数，故下文将直接计算协方差．
为了行文方便、简单，可以将式(14)、(15)统一表示为
(17)
其中
(18)
由随机振动理论[16]及式(17)，结构响应的协方差可以表示为
(19)
其结构响应分量的协方差为
(20)
把白噪声随机过程协方差函数
代入式(20)，可得
(21)
其中,S0为基岩运动加速度功率谱强度系数，δ(τ)为Dirac delta函数．
计算式(21)的积分，可化为
(22)
由式(19)、(22),可以得到结构基于C-P谱的协方差函数为
(23)
令
则式(23)可以简化为
(24)
由式(22)可以将结构响应协方差函数化为结构振动特征值指数函数的线性组合，表达式简单明了，易于计算．由随机振动理论，容易理解，当τ=0时，结构响应的协方差即为结构响应的均方差：
(25)
式中
表示结构响应协方差．
3.2 谱矩
上述协方差函数可提供有关随机过程的时域统计信息，而协方差函数的Fourier变换可用于描述频域统计特性．平稳随机过程的功率谱密度定义为其协方差函数的Fourier变换，称为Wiener-Khinchin关系式[13]．由此关系式，可得结构响应的单边功率谱为
(26)
式中CSn(τ)为结构响应Sn的协方差函数，可由式(19)得到．下标“Sn”表示相对应的响应．
下面来推导谱矩的公式，由谱矩定义，可以得到谱矩αSn,k的表达式为
αSn,k=
SSn(ω)ωkdω,
(27)
k=0,1,2时分别表示响应的0～2阶谱矩．
结合式(22),把式(26)代入式(27)，可以得到
(28)
由随机振动理论[16]可知，k=0,k=2时结构的0阶谱矩、2阶谱矩与结构位移方差和速度方差相对应，可以表示为
(29)
(30)
当k=1时，为结构响应的一阶谱矩，即
(31)
对式(31)进行积分,可得
(32)
因为2阶谱矩存在,则1阶谱矩也必然存在,所以式(32)右边第一项存在:
(33)
即
(34)
故1阶谱矩可以表示为
(35)
式(29)、(35)、(30)为结构响应的0～2阶谱矩的解析解，与传统方法相比，本文方法不需要大量积分，且表达式比较简洁，计算效率大大提高．
4 动力可靠度分析
目前线性结构可靠度计算的最重要方法是基于首次超越破坏准则及Markov超越分布规则的动力可靠度分析．平稳Gauss激励下的线性结构响应，仍为Gauss分布，Crandall提出了能计算结构可靠度的方法，对于双侧相同超越界限的动力可靠度[18-19],其公式为
(36)
式中，Rb(n)为第n年界限值，是正数，按照缓变量方法进行计算；T为一次动荷载的激励时间；
分别为结构某响应S的方差和响应变化率的方差，它们均是动力方程中结构相对位移响应x、相对速度
的确定性函数．
qS为谱参数，其表达式为
(37)
式中，αSn,0,αSn,1,αSn,2表示响应的0～2阶谱矩，可分别按式(29)、(35)、(30)来计算．
5 算 例
建立在地震烈度为8度，二类场地上的单自由度结构，设计分组为2组，场地特征周期为Tg=0.4 s； 层高为4.5 m， 其梁尺寸为400 mm×550 mm， 柱子的截面尺寸为500 mm×500 mm，质量m=4×103 kg，刚度k=3.76×104 N/m，阻尼c=1.226 4×103 N·s/m，阻尼比ξ=0.05，松弛因子μ1=1 s-1,μ2=0.2 s-1 ，地震动采用C-P谱，基岩以上场地土(第一滤波)的阻尼比ξg=0.72，卓越频率ωg=15.71 rad/s；地面处场地土(第二滤波)的阻尼比ξf=ξg，卓越频率ωf=0.15ωg；基岩运动加速度功率谱强度系数S0=6.193×10-3 m2/s3．依据此算例，进行了如下分析．
5.1 结构振动特征的分析
用本文方法，结合MATLAB编程得到了其结构振动特征值与
模态强度系数表，如表1所示,其中等效阻尼比ξi=-Re(pi)/|pi|．
从表1可知，模态3和4的特征值为复共轭项，其对应的等效频率分别与地面处场地土(第二滤波)的卓越频率ωf以及基岩以上场地土(第一滤波)的卓越频率ωg相等，两个模态的等效阻尼比与基岩以上场地土(第一滤波)的阻尼比ξg相等，的模态强度系数不为零，而其余模态强度系数均为0，这是由与场地的地震动仅与场地有关而与结构无关的特点所决定的，也从侧面说明本文方法的正确性．
表1 结构振动特征值与
模态强度系数表
Table 1 Structural vibration eigenvalues and modal strength coefficients of 
mode numbereigenvalue piequivalent frequency |pi|equivalent damping ratio ξixg modal strength coefficient βi1-0.003 80.0041.0002-0.335 10.3351.0003-1.696 7±1.635 3i2.3560.72-1.700 1±0.107 0i4-11.311±10.902 3i15.7100.72-13.011 3±2.454 0i5-0.430 6±38.480 8i38.4830.010
5.2 功率谱对比分析
5.2.1 地震加速度
功率谱
本文方法中，基于C-P谱的地震加速度功率谱，可以由前文中“协方差分析”以及式(26)得
(38)
式中
传统方法中，基于C-P谱的地震加速度功率谱，其表达式为
(39)
式中
表示传统方法的地震加速度功率谱．
5.2.2 结构位移x的功率谱
本文方法中，基于C-P谱的单自由度结构位移功率谱，由前文中“协方差分析”以及式(26)得
(40)
式中λn,i为结构位移强度系数，可由式(16)求得．
传统方法中，其基于C-P谱的单自由度结构位移功率谱以及其系统的幅频特性表达式分别为
(41)
式中
由式(39)可得．
图1与图2给出了本文方法与传统方法中地震加速度与单自由度结构位移响应x的功率谱密度函数对比图，两者完全吻合，说明了本文方法计算结构响应功率谱密度函数的正确性．
图1 地震加速度功率谱对比
Fig. 1 Comparison of seismic acceleration power spectrums
图2 结构位移功率谱对比
Fig. 2 Comparison of structural displacement power spectrums
5.3 计算效率及精度对比分析
传统方法计算功率谱时有积分关系式，为此用本文方法与传统方法比较．根据传统方法，结构响应的0～2阶谱矩的计算公式如下：
(42)
式中i=0,1,2时分别表示结构位移的0～2阶谱矩，
见式(41)．
上述响应值，采用数值方法,ω在[0,∞) rad/s区间进行求解，是无法实现的．由图2可知，随着频率ω的增大，功率谱越来越小，同时在[0,50] rad/s区间内，功率谱有显著值．为此对式(42)进行积分，选[0,50] rad/s为积分区间，为了比较积分精度，积分步距分别取：① 为0.05 rad/s；② 为0.1 rad/s；③ 为0.5 rad/s．
本文方法得到的谱矩公式由式(29)、(35)、(30)可以得到．表2列出了ω在[0,50] rad/s积分区间下，本文方法和传统方法所得谱矩计算值和CPU耗时对比以及两种方法的误差分析．
表2 本文方法和传统方法的谱矩计算值和CPU耗时对比以及误差分析([0,50] rad/s)
Table 2 The comparison of the calculated values of spectral moments, CPU time costs and errors between this method and traditional method([0,50] rad/s)
αx,0αx,1αx,2CPU time tc/serror η/%αx,0αx,1αx,2this paper5.699 290 490 14E-62.121 897 659 596E-48.063 925 776E-30.000 921---①5.685 095 538 72E-62.113 675 579 527E-48.014 533 291E-30.037 073-0.249 06-0.387 48-0.612 51②5.685 095 167 17E-62.113 675 409 436E-48.014 532 519E-30.014 737-0.249 07-0.387 49-0.612 52③5.743 131 496 77E-62.131 934 649 841E-48.125 004 952E-30.006 8580.769 240.473 020.757 44
从表2可看出，本文方法的CPU耗时比传统方法的分别提高了40倍、16倍和7倍，同时随着积分步长的增大，误差越来越大，说明了本文方法所得封闭解的正确性．传统方法与本文方法的误差出现正负值与功率谱曲线的凹凸有关．从图2可以看出若功率谱曲线为凸，则用本文方法求得功率谱面积小于传统方法所得面积，有
若功率谱曲线为凹，则用本文方法求得的功率谱面积大于传统方法法所得的面积，故有
5.4 可靠度计算
设定结构的层间位移为H/500，层高H为5.0 m，按上下界限均取10 mm，动荷载的激励时间T为15 s，结合上文算例中给出的数据来计算可靠度．
根据式(33)，可以得到本算例中的
数值较小，可以忽略不记，说明1阶谱矩的计算是正确的．再利用式(36)、(37)，可以得到谱参数qS=0.999 0，可靠度Pf=0.999 9．所得结果较为理想，说明本文方法的正确性．
6 结 论
本文针对非黏滞阻尼结构基于C-P谱随机激励下目前解法比较复杂的问题，提出了简明的封闭解法，结论如下：
1) 非黏滞阻尼模型利用核函数的卷积形式表示结构的耗能部分，但其在微分型的动力方程中不易求解，本文推导出其精确等效的一阶微分型本构关系；C-P谱的滤波振动方程可将复杂的地面运动表示成基于白噪声的随机地震动；本文利用非黏滞阻尼模型微分型本构关系和C-P谱的滤波方程重构非黏滞阻尼结构的地震动方程，基于随机振动理论获得了结构响应的协方差函数、功率谱和0～2阶谱矩的统一简明封闭解．
2) 本文将得到的结构响应的功率谱表达式与传统方法的相对比，表明本文方法所得表达式的正确性与简洁性，且在随机响应方差和谱矩分析时能获得简明的封闭解，避免了复杂的积分运算，具有计算效率较高的特点．
3) 由于时域法在分析结构随机响应时需要简明的激励协方差，故本文利用C-P谱的滤波方程，将复杂的C-P谱表示为白噪声激励的线性微分方程组，从而达到简化分析，获得简明封闭解的目的．因此，本文方法可适用于诸如K-T谱、欧进萍谱等作用下的线型结构随机地震动响应分析．
4) 本文是针对单自由度非黏滞阻尼结构基于C-P谱的随机响应分析，由于采用C-P滤波方程和非黏滞阻尼的一阶微分方程，重构运动方程后即为8参数的微分方程组．因此，本文解法就是针对“多自由度”的求解．由于多自由度非黏滞阻尼结构可以表示为阻尼系数矩阵与指数型核函数的卷积的乘积形式，故本文方法也适用于多自由度非黏滞阻尼结构的随机响应分析．
致谢 本文作者衷心感谢广西科技大学创新团队支持计划项目对本文的资助．
[1] 欧进萍, 刘会仪. 基于随机地震动模型的结构随机地震反应谱及其应用[J]. 地震工程与工程振动, 1994, 14(1): 14-23.(OU Jinping, LIU Huiyi. Random earthquake response spectrum of structures based on random ground motion model and its application[J]. Earthquake Engineering and Engineering Vibration, 1994, 14(1): 14-23.(in Chinese))
[2] 欧进萍, 牛荻涛，杜修力. 设计用随机地震动的模型及其参数确定[J]. 地震工程与工程振动, 1991, 11(3): 45-53.(OU Jinping, NIU Ditao, DU Xiuli. Model and parameter determination of random ground motion for design[J]. Earthquake Engineering and Engineering Vibration, 1991, 11(3): 45-53.(in Chinese))
[3] KANAI K. An empirical formula for the spectrum of strong earthquake motions[J]. Bulletin of the Earthquake Research Institute, 1961, 39(1): 85-95.
[4] CLOUGH R W, PENZIEN J. Dynamics of Structures[M]. New York: McGraw-Hill Inc, 1975.
[5] 李创第, 丁晓华, 陈俊忠. 基础隔震结构基于Clough-Penzien谱随机地震响应分析的复模态法[J]. 振动与冲击, 2006, 25(5): 162-165,199.(LI Chuangdi, DING Xiaohua, CHEN Junzhong. Complex modal method for analysis of random earthquake response of structures with base isolation on the basis of Clough-Penzien spectrum[J]. Journal of Vibration and Shock, 2006, 25(5): 162-165,199.(in Chinese))
[6] 李春祥, 熊学玉. 基于Kanai-Tajimi/Clough-Penzien模型时MTMD的动力特性[J]. 振动与冲击, 2002, 21(4): 41-45,121.(LI Chunxiang, XIONG Xueyu. Dynamic characteristics of MTMD based on Kanai-Tajimi Clough-Penzien model[J]. Journal of Vibration and Shock, 2002, 21(4): 41-45,121.(in Chinese))
[7] 潘玉华, 王元丰. 基于微分求积法求解非粘滞阻尼系统的时程响应[J]. 计算力学学报, 2011, 28(4): 517-522,529.(PAN Yuhua, WANG Yuanfeng. Time history response of non viscous damping system based on differential quadrature method[J]. Chinese Journal of Computational Mechanics, 2011, 28(4): 517-522,529.(in Chinese))
[8] 沈洪宇, 段忠东, 李惠. 基于Volterra积分方程Taylor展开的卷积型非粘滞阻尼系统响应的时程分析方法[J]. 应用力学学报, 2018, 35(2): 261-266.(SHEN Hongyu, DUAN Zhongdong, LI Hui. Time history analysis method of response of convolution type non viscous damping system based on Taylor expansion of Volterra integral equation[J]. Journal of Applied Mechanics, 2018, 35(2): 261-266.(in Chinese))
[9] 李春江. 卷积型非粘滞阻尼模型在建筑结构中的适用性研究[D]. 硕士学位论文. 哈尔滨: 哈尔滨工业大学, 2017.(LI Chunjiang. Study on applicability of convolutional non viscous damping model in building structure[D]. Master Thesis. Harbin: Harbin University of Technology, 2017.(in Chinese))
[10] 王禹. 指数型非粘滞阻尼模型参数识别的理论与试验研究[D]. 硕士学位论文. 哈尔滨: 哈尔滨工业大学, 2017.(WANG Yu. Theoretical and experimental research on parameter identification of exponential non viscous damping model[D]. Master Thesis. Harbin: Harbin University of Technology, 2017.(in Chinese))
[11] 段忠东, 沈洪宇. 非粘滞阻尼系统时程响应分析的精细积分方法[J]. 计算力学学报, 2009, 26(5): 638-644.(DUAN Zhongdong, SHEN Hongyu. Precise integration method for time history response analysis of non-viscous damping system[J]. Journal of Computational Mechanics, 2009, 26(5): 638-644.(in Chinese))
[12] LIU Q M. Explicit computational method of dynamic response for non-viscously damped structure systems[J]. Mechanics Research Communications, 2016, 71: 48-55.
[13] 方同. 工程随机振动[M]. 北京: 国防工业出版社, 1995.(FANG Tong. Random Vibration of Engineering[M]. Beijing: National Defense Industry Press, 1995.(in Chinese))
[14] 方同, 张天舒. 计算响应谱矩的复模态方法[J]. 力学学报, 1989, 21(1): 65-74.(FANG Tong, ZHANG Tianshu. Complex mode method for calculating response spectral moment[J]. Acta Mechanica Sinica, 1989, 21(1) 65-74.(in Chinese))
[15] 林家浩. 随机振动的虚拟激励法[M]. 北京: 科学出版社, 2004.(LIN Jiahao. Virtual Excitation Method of Random Vibration[M]. Beijing: Science Press, 2004.(in Chinese))
[16] 赵雷, 刘宁国. 大跨度叠合梁斜拉桥多维多点非平稳随机地震响应分析[J]. 应用数学和力学, 2017, 38(1): 118-125.(ZHAO Lei, LIU Ningguo. Non-stationary random seismic analysis of large-span composite beam cable-stayed bridges under multi-support and multi-dimentional earthquake excitations[J]. Applied Mathematics and Mechanics, 2017, 38(1): 118-125.(in Chinese))
[17] 李创第, 昌明静, 柏大炼, 等. 六参数粘弹性阻尼耗能减震系统非平稳地震响应分析的精细积分法[J]. 广西科技大学学报, 2019, 30(4): 15-22,41.(LI Chuangdi, CHANG Mingjing, BAI Dalian, et al. Precise integration method for the analysis of non-stationary seismic response of six parameter viscoelastic damping energy dissipation system[J]. Journal of Guangxi University of Science and Technology, 2019, 30(4): 15-22,41.(in Chinese))
[18] LI L, HU Y, WANG X. Design sensitivity analysis of dynamic response of nonviscously damped systems[J]. Mechanical Systems and Signal Processing, 2013, 41(1/2): 613-638.
[19] 曹宏, 李桂青, 李秋胜. 结构动力可靠性理论及其应用[M]. 北京: 地震出版社, 1993.(CAO Hong, LI Guiqing, LI Qiusheng. Theory and Application of Structural Dynamic Reliability[M]. Beijing: Earthquake Press, 1993.(in Chinese))