基于MATLAB语言的地震波动力特性分析方法

基于MAT LAB语言的地震波动力特性分析方法

马乐为1,钟骁珺1,谢异同1,张同亿2

(11西安建筑科技大学土木工程学院,陕西西安　710055;21中元国际工程设计研究院,北京　1000)

摘　要:地震波频谱组成的研究主要是利用Fourier分析或小波理论,而地震波对建筑物动力

反应的研究则是通过拟加速度反应谱来进行.现行分析方法的数值化过程必须通过计算机编

程来实现,而采用常规的编程语言实现起来往往较为复杂.基于此,作者提出了基于MA T2

L AB平台的分析方法,可以大大提高地震波动力特性分析的效率,便于分析结果的图形化输

出,为工程结构抗震计算提供有益的参考.

关键词:MA TL AB语言;地震波频谱组成;拟加速度反应谱;Fourier分析;小波理论

中图分类号:P315.3+1;TP312　　　　　　文献标识码:A

0　引言

地震波的动力特性是指其频谱组成以及对不同建筑物的动力反应特点,地震波的这些特性对建筑结构抗震计算是非常重要的.一般而言,在地震波三要素振幅、持时和频率中以第三者最为复杂.目前,有关地震波频谱特性分析或动力反应研究的理论主要包括时域或频域内的分析方法,即时程分析法、Fourier 变换法、拉普拉斯变换法(Laplace t ransfer)以及小波分析法等;而地震波对建筑物动力反应的研究则往往是通过拟加速度反应谱来进行的.

结构地震作用计算中的反应谱理论是各国规范普遍采用的计算原则,其中以地震波的加速度反应谱最为有效.绘制加速度反应谱曲线的常规方法是时域内的Duhamel积分,其数值积分可利用线性加速度法或Willsion2θ法等.上述分析方法的数值化过程必须通过计算机编程来实现,而采用常规的编程语言如Fort ran或C语言实现起来往往较为复杂.其实,动力方程的求解也可以通过拉普拉斯变换来实现,而且利用MA TL AB提供的函数可使求解过程大大简化.

作为目前科技人员普遍使用的计算工具,MA TL AB最大的优点在于其极为丰富的函数库以及几乎能够满足各学科要求的工具箱(Toolbox),人们可以用最少的语句完成最复杂的计算[1].本文上述的Fou2 rier变换、拉普拉斯变换以及小波分析等方法在MA TL AB中均可以快速实现,从而大大提高地震波动力特性的分析效率,便于分析结果的图形化输出,为工程结构抗震计算提供有益的参考.

1　基于Matlab的地震拟加速度反应谱

单自由度线性系统的运动微分方程可表示为:

m¨x+c x+kx=f(t)(1)式中m、c、k及f(t)分别表示质点的质量、阻尼系数、刚度和受到的激励;x为质点的位移.则质点的位移x(t)、速度 x(t)、加速度¨x(t)及激励f(t)的拉普拉斯变换分别为[2]:

3收稿日期:2008-10-21

作者简介:马乐为(1971-),男,陕西省西安市人,副教授,博士,研究方向:结构工程抗震

基金项目:国家自然科学基金项目(50478045)

第1期马乐为等:基于MA TL AB 语言的地震波动力特性分析方法X (s )≡Ls (t )≡∫∞0e -st x (t )d t

L x (t )=sX (s )-x (0)

L ¨x (t )=s 2X (s )-sx (0)- x (0)

F (s )=∫∞0e -st f (t )d t

(2)对于反应谱问题,由于外部激励f (t )为质点质量m 与地震波加速度¨x 0(t )的乘积,因此质点的绝对加速度a (t )可近似表示为:

a (t )=¨x 0(t )+¨x (t )=-k m x (t )=-ω2x (t )(3)

此时,如设质点为单位质量,即m =1,则求解反应谱的绝对加速度极值的问题就变成了求解体系最大位移响应的问题,因为两者只差一个常数-k.再利用拉普拉斯变换的性质及系统的零初始条件x (0)= x (0)=0,有:

X (s )=F (s )D (s )=∫∞0

-e -st ¨x (t )d t s 2+cs +k

(4)式中1/D (s )在控制理论中称为传递函数,其物理意义为系统输入与输出的拉普拉斯变换之比.而反应谱中质点绝对最大加速度值|a (t )|max 则可以由拉普拉斯逆变换方便的求解如下:

|a (t )|max =2

πT |L -1X (s )|(5)

如前所述,Matlab 语言具有丰富的函数库,这为解决拉普拉斯变换提供了极为方便的途径.其中tf ()函数可以产生上述的传递函数1/D (s ),而lsim ()函数则可以对时域系统响应进行仿真,得到线性时不变连续系统的冲激响应.由此,编制求解地震波拟加速度反应谱的Matlab 函数pars ()如下所示,此函数的功能为已知一条地震波加速度记录et hw v 、阻尼比dm p r 及地震波记录的时间间隔tmitvl ,即可得到我国

《抗震规范》[3]要求的6s 内的拟加速度反应谱曲线:

f unction pars (ethwv ,tmitvl ,dmpr )

%function pars (ethwv ,tmitvl ,dmpr )

%Purpose :To draw p seudo 2acceleration response spectrum figure in 6s

%variable names :

%ethwv acceleration of earthquake wave record in one column form

%tmitvl time interval of the earthquake wave record

%dmpr damping ratio

t =0:tmitvl :(length (ethwv )-1)3tmitvl ;

Prd =0:0.01:6;Prd =Prd ′;

omgar =23pi./Prd ;

for i =1:length (omgar )

　　SYS =tf (1,[1,23dmpr 3omgar (i ),omgar (i )3omgar (i )]);

　　disp (i )=max (abs (lsim (SYS ,ethwv ,t )));

　　spd (i )=23pi.3disp (i )/Prd (i );

　　acc (i )=23pi.3spd (i )/Prd (i );

end

plot (Prd ,acc );

在此输入Elcent ro 、taft 及Nort hbridge 地震波加速度记录,并取阻尼比为0.05,3条波的时间间隔均为0.02s ,计算结果如图1所示.与一些权威结果[4]相比较,证明图1曲线正确,说明加速度反应谱的拉普・

131・

陕西科技大学学报第27

卷

图1　Elcent ro 、taft 及Nort hbridge

地震波加速度反应谱曲线

拉斯解法在理论上是可行的,上述Matlab 源程序除函数说明

语句外,有效程序语句仅10行,可见Matlab 程序平台的高效

性.

2　基于MAT LAB 的地震波谱密度分析

形象地说,功率谱就是2个坐标轴,横轴是频率,竖轴是功

率,功率谱密度就是在整个频谱上功率大小分布随频率变化的

一个量.对于地震动来说,功率谱密度函数可以反映地震波在

各个频率成分上振动能量(及振幅)的大小,是直接提供地震动

激励有关信息的一种有用形式.对于结构抗震计算、拟动力和振动台试验中的地震波选取以及构造人工地震波等问题,功率谱密度均可以作为一种合理性的判别依据.

地震波谱密度分析的数学基础是Fourier 变换.当样本总量满足2的n 次幂时,可直接采用快速Fou 2rier 变换,不满足此条件时,也可用末位补零的方法进行快速Fourier 变换,而地震波加速度记录中某一频率成分所占的比重则可用功率谱密度来反映.设地震波加速度记录为一个离散函数样本x (n ),其Fourier 变换X (k )及功率谱密度S (x )可按下式计算:

X (k )=

∑N

n =1x (n )e -2in (k-1)(n-1)/n S (x )=∑N -1k =-N +1R n e

-2i πf k

(6)

式中,N 为地震波加速度记录总数,取N =2m ,k =1,2,…,N/2,R 为自相关函数.

由式(6)可知,对离散的地震动加速度信号进行快速傅里叶变换即可得到其线性谱X (k ),而离散信号的自功率谱则可以表示为其线性谱乘其共轭线性谱.在Matlab 语言中fft ()是离散序列的Fourier 变换函数,conj ()是复数的共轭函数,直接利用这两个函数便可迅速求解给定地震动加速度信号的功率谱密度.具体函数如下:

f unction psdd (ethwv )

n =2^(nextpow2(length (ethwv )))

y =fft (ethwv ,n );

p =y.3conj (y )/n ;

Elcentro 波及3条上海人工波SHW1、SHW2、SHW3的加速度记录如图2所示,采用上述p sdd ()函数分别求解,即可得到如图3所示的功率谱密度图.

3　基于MAT LAB 的小波分解

由于Fourier 变换是一种整个频域内的变换,我们只能通过Fourier 变换了解到某条地震波在整个波形内的主要频率组成,而无法了解某个频率成分在地震波的哪一时刻所占的比重最大,因此无法完成对于拟动力试验所需要的波的某一时段的截取.

有幸的是小波(Wavelet )分析这种新型工具提供了一种可以将频域和时域分析联系起来的方法,它能够表述地震波信号的时频局部性质,从而可以使我们了解到一条地震波中哪些时刻的频率成分与结构的基本周期相近,为拟动动力试验中地震波截取进行定性的判断,进而可以避免拟动力试验中地震波选择时仅以峰值作为主要条件的弊端.

小波分析的数学基础十分复杂,本文仅以MA TL AB 提供的常用小波函数Daubechies (dbN )小波进行地震波分解,Daubechies 函数是由世界著名的小波分析学者Inrid Daubechies 构造的小波函数,除了・231・

第1期马乐为等:基于MA TL AB

语言的地震波动力特性分析方法　　　图2　地震波原始记录　　　　　　　　　　　图3　

地震波的功率谱密度

图4　SHW1波小波分解

db1(即haar 小波)外,其他小波没有明确的表达式,但

转换函数h ()的平方模是很明确的,它具有以下基本性

质:正交性、双正交性、紧支撑性、连续小波变换、离散

小波变换,支撑长度为2N -1,滤波器长度为2N ,并呈

近似对称性[528].

本文采用db3并利用MA TLAB 对SHW1进行6

层分解,如图4所示.由于SHW1的采样频率为0.01

s ,故小波分析时根据采样定埋,其频率上限即Nyquist

频率为1/(2×0.01)=50Hz ,6层分解按2的幂次(从

1到6)的频率范围依次为50/21即25～50Hz 、50/22

即12.5～25Hz 、50/23即6.25～12.5Hz 、50/24即

3.125～6.25Hz 、50/25即1.56～3.125Hz 、50/26即

0.78～1.56Hz.按以上方法编制的MA TL AB 程序如

下:

load shw1;

s =shw1;ls =length (shw1);

[c ,l ]=wavedec (s ,6,‘db3’

);subplot (7,1,1);plot (s );

title (‘SHW12‘SHW1各层小波重构’

);Y label (‘

SHW1’);for i =1∶6

　　decmp =wrcoef (‘a ’,c ,l ,‘db3’,7-i );　　subplot (7,1,i +1);

　　plot (decmp );

・

331・

陕西科技大学学报第27卷

　　Ylabel ([‘a ’,num2str (7-i )]);

End 图中横坐标为时域的样本总量,乘以采样频率0.01s 后即为SHW1的总持时10.25s ;纵坐标为频率分解的振幅,如果已知结构物的基本周期为0.24s ,频率为4.37Hz ,则模型基频落在a 4频带范围内,而图中纵坐标表示了这一频率的幅值.由此可见,SHW1地震波在第2～6s 时段内的频率成分与结构模型相近,并且这一时段内包含了地震波加速度峰值最大的时段,故在拟动力试验截取地震波时应首先考虑这一时段.

4　结束语

拟加速度反应谱在结构地震作用计算中普遍采用,绘制拟加速度反应谱曲线的常规方法是时域内的Duhamel 积分.本文的研究表明,根据拉普拉斯变换原理,并利用MA TLAB 提供的传递函数tf ()及时域系统响应仿真函数lsim ()可使反应谱曲线的求解过程大大简化.

参考文献

[1]陈怀琛.MA TL AB 及其在理工课程中的应用指南[M ].西安:西安电子科技大学出版社,2000:125.

[2]方　同,薛　璞.振动理论及应用[M ].西安:西北工业大学出版社,1998:84294.

[3]G B5001122001.建筑抗震设计规范[S].北京:中国建筑工业出版社,2001.

[4]武藤清[日].结构物动力设计[M ].北京:中国建筑工业出版社,1984:3382350.

[5]胡昌华.基于MA TLAB 的系统分析与设计—小波分析[M ].西安:西安电子科技大学出版社,1999.

[6]李　健,梁　琨.基于小波多分辨率分析的快速图像匹配[J ].陕西科技大学学报,2006,24(6):81284.

[7]马乐为,吴敏哲,谢异同.基于脉冲频响函数的MA TLAB 动力方程求解方法[J ].世界地震工程,2003,19(2):1682171.

[8]孙　凌,吴知丰.结构随机响应的一种改进的概率谱叠加方法[J ].世界地震工程,2003,19(2):1562162.

DY NAMIC CHARACTERISTICS ANALYSIS METH OD OF THE EARTHQUAKE WAVES ON THE BASIS OF MAT LAB LANGUAGE

MA Le 2wei 1,ZHON G Xiao 2jun 1,XIE Y i 2tong 1,ZHAN G Tong 2yi 2

(11School of Civil Eng.,Xi ′an Univ.of Arch.&Tech.,Xi ′an 　710055,China ;21IPPR International Engi 2neering ,Beijing 　1000,China )

Abstract :At present ,t he main t heoretic met hod of eart hquake wave spect ral composition main are by means of Fo urier or wavelet analysis.For st uding of t he dynamic response of buildings on t he impact of eart hquake waves ,t he p seudo 2acceleration response spect rum will usually be int roduced.The digital p rocessing of above t heoretic met hods should be used of comp utting p rogramme.It is difficult for t he usuall comp uter programme to finish t hose ana 2lytic p rocesses.The article point s out t hat t he analyzing efficiency of dynamic characteristics of t he eart hquake wave will be improved greatly and t he grap hic result will be simply got on basis of Matlab language ,t he result of t his article can be benefited for calculation of anti 2eart hquake of buildings ′st ruct ure.

K ey w ords :Matlab language ;spect ral compo sition of eart hquke waves ;p seudo 2acceleration response spect rum ;Fourier ′s analysis ;wavelet t heory

・

431・