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

Ｏ ２００ ４００ ６００ ８００
。
１０００
１２∞
即１２．５～２５ Ｈｚ、５０／２３即６．２５～１２．５ Ｈｚ、５０／２４即
３．１２５～６．２５ Ｈｚ、５０／２５即１．５６～３．１２５ Ｈｚ、５０／２６即
０．７８，－－．１．５６
Ｈｚ．按以上方法编制的ＭＡＴＬＡＢ程序如
。
下：
ｌｏａｄ ｓｈｗｌ；
∞ｏ ∞Ｏ ∞Ｏ ∞ｏ ∞ｏ ∞Ｏ ∞０
№。１ ・１ ３０・
陕西科技大学学报
ＪＯＵＲＮＡｌ。ＯＦ ＳＨＡＡＮＸＩ ＵＮＩＶＥＲＳＩＴＹ ＯＦ ＳＣＩＥＮＣＥ ８Ｌ ＴＥＣＨＮＯＬＯＧＹ
Ｆｅｂ．２００９ Ｖ０１．２７
。文章编号：１０００—５８１１（２００９１０１—０１３０－０５
基于ＭＡＴＬＡＢ语言的地震波动力特性分析方法
马乐为１，钟骁瑶１，谢异同１，张同亿２
口（ｆ）一；ｏ（￡）＋三（￡）一一与（ｚ）＝一ｃｃＪ２ｘ（ｔ）
பைடு நூலகம்
（３）
此时，如设质点为单位质量，即ｍ一１，则求解反应谱的绝对加速度极值的问题就变成了求解体系最 大位移响应的问题，因为两者只差一个常数一志．再利用拉普拉斯变换的性质及系统的零初始条件ｚ（ｏ）＝ ｚ（Ｏ）一ｏ，有：
矾）：器一了：－－丽ｅ－＇ｘ（ｔ）ｄｔ
ｆｕｎｃｔｉｏｎ
ｐａｒｓ（ｅｔｈｗｖ，ｔｍｉｔｖｌ，ｄｍｐｒ）
％ｆｕｎｃｔｉｏｎ ｐａｒｓ（ｅｔｈｗｖ，ｔｍｉｔｖｌ，ｄｍｐｒ） ％Ｐｕｒｐｏｓｅ：Ｔｏ ｄｒａｗ ｐｓｅｕｄｏ－ａｃｃｅｌｅｒａｔｉｏｎ ％ｖａｒｉａｂｌｅ
Ｂｆｌｌｎｅｓ：
ｒｅｓｐｏｎｓｅ
ｓｐｅｃｔｒｕｍ ｆｉｇｕｒｅ ｉｎ ６
Ｓ
％ｅｔｈｗｖ ａｃｃｅｌｅｒａｔｉｏｎ ｏｆ ｅａｒｔｈｑｕａｋｅ
ｅｎｄ ｐｌｏｔ（Ｐｒｄ，ａｃｅ）；
在此输入Ｅｌｃｅｎｔｒｏ、ｔａｆｔ及Ｎｏｒｔｈｂｒｉｄｇｅ地震波加速度记录，并取阻尼比为０．０５，３条波的时间间隔均 为０．０２ ｓ，计算结果如图１所示．与一些权威结果［４］相比较，证明图１曲线正确，说明加速度反应谱的拉普
万方数据
－＿ｌｌ＿＿＿－●●ｌＩ●ｌ●＿●－＿＿Ｉｌ＿＿●●－ｌＩＩ－Ｉ ＩＩｌ＿●－●－●■ＩＩ＿ｌ＿●Ｉ＿－＿＿＿ｌ●＿●＿●●●●●ｌｌ－＿■＿－ｌ－●●＿ｌ＿●＿●Ｉｌｌ●●ＩＩ＿＿＿Ｉｌ●＿＿－
ｆｏｒｉ－－－－１
ｌ
６
ｄｅｃｍｐ＝ｗｒｃｏｅｆ（‘ａ’，ｃ。ｌ，‘ｄｂ３’，７一ｉ）； ｓｕｂｐｌｏｔ（７，１，ｉ＋１）； ｐｌｏｔ（ｄｅｃｍｐ）；
万方数据
・１３４． ！！！！！！！！！！！Ｉ＝
陕西科技大学学报
。
第２７卷 ＩＩ！！！！！！！！！！竺
Ｙｌａｂｅｌ（［‘ａ’，ｎｕｍ２ｓｔｒ（７－－ｉ）］）；
０ ２００ ４００ ６００ ８００ １０００
１２∞
样本总量
ｓ＝ｓｈｗｌ；ｉｓ＝ｌｅｎｇｔｈ（ｓｈｗｌ）；
［ｃ，１］－－－－－ｗａｖｅｄｅｃ（ｓ，６，‘ｄｂ３’）；
ｓｕｂｐｌｏｔ（７，１，１），ｐｌｏｔ（ｓ）｝
图４
ＳＨＷｌ波小渡分解
ｔｉｔｌｅ（‘ＳＨＷＩ“ＳＨＷｌ各层小波重构’）｝ Ｙｌａｂｅｌ（‘ＳＨＷｌ’）：
ｗａｖｅ ｒｅｃｏｒｄ ｉｎ ｗａｖｅ
ｏｎｅ
ｃｏｌｕｍｎ ｆｏｒｍ
％ｔｍｉｔｖｌ ｔｉｍｅ ｉｎｔｅｒｖａｌ ｏｆ ｔｈｅ ｃａｒｔｈｑｕａｋｅ ％ｄｍｐｒ ｄａｍｐｉｎｇ
ｒａｔｉｏ
ｒｅｃｏｒｄ
ｔ＝０：ｔｍｉｔｖｌ：（１ｅｎｇｔｈ（ｅｔｈｗｖ）一１）＊ｔｍｉｔｖｌ； Ｐｒｄ—Ｏ：０．０１：６；Ｐｒｄ—Ｐｒｄ’； ｏｍｇａｒ一２。ｐｉ．／Ｐｒｄ； ｆｏｒ
ｆｕｎｃｔｉｏｎ ｐｓｄｄ（ｅｔｈｗｖ） ｎ＝２‘（ｎｅｘｔｐｏｗ２（１ｅｎｇｔｈ（ｅｔｈｗｖ））） ｙ＝ｆｆｔ（ｅｔｈｗｖ，ｎ）；
ｔ晶１
ｐ＝ｙ－＊ｃｏｎｊ（ｙ）／ｎ；
Ｅｌｃｅｎｔｒｏ波及３条上海人工波ＳＨＷｌ、ＳＨＷ２、ＳＨＷ３的加速度记录如图２所示，采用上述ｐｓｄｄ（）函 数分别求解，即可得到如图３所示的功率谱密度图． ３基于ＭＡＴＬＡＢ的小波分解 由于Ｆｏｕｒｉｅｒ变换是一种整个频域内的变换，我们只能通过Ｆｏｕｒｉｅｒ变换了解到某条地震波在整个波 形内的主要频率组成，而无法了解某个频率成分在地震波的哪一时刻所占的比重最大，因此无法完成对于 拟动力试验所需要的波的某一时段的截取． 有幸的是小波（Ｗａｖｅｌｅｔ）分析这种新型工具提供了一种可以将频域和时域分析联系起来的方法，它能 够表述地震波信号的时频局部性质，从而可以使我们了解到一条地震波中哪些时刻的频率成分与结构的 基本周期相近，为拟动动力试验中地震波截取进行定性的判断，进而可以避免拟动力试验中地震波选择时 仅以峰值作为主要条件的弊端． 小波分析的数学基础十分复杂，本文仅以ＭＡＴＬＡＢ提供的常用小波函数Ｄａｕｂｅｅｈｉｅｓ（ｄｂＮ）Ｌ［＇＊波进 行地震波分解，Ｄａｕｂｅｅｈｉｅｓ函数是由世界著名的小波分析学者Ｉｎｒｉｄ Ｄａｕｂｅｅｈｉｅｓ构造的小波函数，除了
万方数据
第１期
ＩＩ
马乐为等：基于ＭＡＴＩ。ＡＢ语言的地震波动力特性分析方法
・１３３・
ｋ山山■…■．｛２●：．
１ｌｆｒ＂”…”唧一ｊ…。。；’
） １０ ２０ ３０ ４０ ５０ ６Ｉ
爵 霄 哀 霉
频率／Ｈｚ
图２
地震波原始记录
；
图３地震波的功率谱密度 长
．
ｄｂｌ（１ｌｐ ｈａａｒ小波）外，其他小波没有明确的表达式，但
Ｉ口（￡）Ｉ一＝－ｚ－亍兀ｌ Ｌ－１Ｘ（５）Ｉ
（４）
式中１／Ｄ（ｓ）在控制理论中称为传递函数，其物理意义为系统输入与输出的拉普拉斯变换之比．而反 应谱中质点绝对最大加速度值Ｉ口（￡）ｌ一则可以由拉普拉斯逆变换方便的求解如下：
（５）
如前所述，Ｍａｔｌａｂ语言具有丰富的函数库，这为解决拉普拉斯变换提供了极为方便的途径．其中玎（） 函数可以产生上述的传递函数１／Ｄ（ｓ），而ｌｓｉｍ（）函数则可以对时域系统响应进行仿真，得到线性时不变 连续系统的冲激响应．由此，编制求解地震波拟加速度反应谱的Ｍａｔｌａｂ函数ｐａｒｓ（）如下所示，此函数的 功能为已知一条地震波加速度记录ｅｔｈｗｖ、阻尼比＆ｎｐｒ及地震波记录的时间间隔ｔｍｉｔｖｌ，即可得到我国 《抗震规范》［３］要求的６ Ｓ内的拟加速度反应谱曲线：
转换函数＾（）的平方模是很明确的，它具有以下基本性 质：正交性、双正交性、紧支撑性、连续小波变换、离散 小波变换，支撑长度为２Ｎ—ｌ，滤波器长度为２Ｎ，并呈 近似对称性［洲． 本文采用ｄｂ３并利用ＭＡＴＬＡＢ对ＳＨＷＩ进行６
．
层分解，如图４所示．由于ＳＨＷＩ的采样频率为０．Ｏｌ孽． ｓ，故小波分析时根据采样定埋，其频率上限即Ｎｙｑｕｉｓｔ鞭。 频率为１／（２）４０．０１）＝５０ Ｈｚ，６层分解按２的幂次（从 １到６）的频率范围依次为５０／２１即２５～５０ Ｈｚ、５０／２２
关键词：ＭＡＴＬＡＢ语言；地震波频谱组成；拟加速度反应谱；Ｆｏｕｒｉｅｒ分析；小波理论
中图分类号：Ｐ３１５．３＋１；ＴＰ３１２ 文献标识码：Ａ
０引言 地震波的动力特性是指其频谱组成以及对不同建筑物的动力反应特点，地震波的这些特性对建筑结 构抗震计算是非常重要的．一般而言，在地震波三要素振幅、持时和频率中以第三者最为复杂．目前，有关 地震波频谱特性分析或动力反应研究的理论主要包括时域或频域内的分析方法，即时程分析法、Ｆｏｕｒｉｅｒ 变换法、拉普拉斯变换法（Ｌａｐｌａｃｅ ｔｒａｎｓｆｅｒ）以及小波分析法等；而地震波对建筑物动力反应的研究则往 往是通过拟加速度反应谱来进行的． 结构地震作用计算中的反应谱理论是各国规范普遍采用的计算原则，其中以地震波的加速度反应谱 最为有效．绘制加速度反应谱曲线的常规方法是时域内的Ｄｕｈａｍｅｌ积分，其数值积分可利用线性加速度 法或Ｗｉｌｌｓｉｏｎ－Ｏ法等．上述分析方法的数值化过程必须通过计算机编程来实现，而采用常规的编程语言如 Ｆｏｒｔｒａｎ或Ｃ语言实现起来往往较为复杂．其实，动力方程的求解也可以通过拉普拉斯变换来实现，而且 利用ＭＡＴＬＡＢ提供的函数可使求解过程大大简化． 作为目前科技人员普遍使用的计算工具，ＭＡＴＬＡＢ最大的优点在于其极为丰富的函数库以及几乎 能够满足各学科要求的工具箱（Ｔｏｏｌｂｏｘ），人们可以用最少的语句完成最复杂的计算［１］．本文上述的Ｆｏｕ— ｒｉｅｒ变换、拉普拉斯变换以及小波分析等方法在ＭＡＴＬＡＢ中均可以快速实现，从而大大提高地震波动力 特性的分析效率，便于分析结果的图形化输出，为工程结构抗震计算提供有益的参考． １基于Ｍａｔｌａｂ的地震拟加速度反应谱 单自由度线性系统的运动微分方程可表示为：
拉斯解法在理论上是可行的，上述Ｍａｔｌａｂ源程序除函数说明 语句外，有效程序语句仅１０行，可见Ｍａｔｌａｂ程序平台的高效
＼
能 馏 真
∞
．１３２・
陕西科技大学学报
第２７卷
性． ２基于ＭＡＴＬＡＢ的地震波谱密度分析 形象地说，功率谱就是２个坐标轴，横轴是频率，竖轴是功 率，功率谱密度就是在整个频谱上功率大小分布随频率变化的 一个量．对于地震动来说，功率谱密度函数可以反映地震波在 各个频率成分上振动能量（及振幅）的大小，是直接提供地震动 图１ Ｅｌｃｅｎｔｒｏ、ｔａｆｔ及Ｎｏｒｔｈｂｒｉｄｇｅ 地震波加速度反应谱曲线 激励有关信息的一种有用形式．对于结构抗震计算、拟动力和 振动台试验中的地震波选取以及构造人工地震波等问题，功率
赫＋＆＋妇＝，（￡）
移ｚ（ｆ）、速度５ｃ（ｔ）、加速度王（ｆ）及激励厂（￡）的拉普拉斯变换分别为‘２１：
（１）
式中ｍ、Ｃ、尼及厂（￡）分别表示质点的质量、阻尼系数、刚度和受到的激励；ｚ为质点的位移．则质点的位
・收稿日期：２００８—１０一２１ 作者简介：马乐为（１９７１一）。男，陕西省西安市入，副教授，博士，研究方向：结构工程抗震 基金项目：国家自然科学基金项目（５０４７８０４５）
（１．西安建筑科技大学土木工程学院，陕西西安７１００５５；２．中元国际工程设计研究院，北京１０００８９）
摘要：地震波频谱组成的研究主要是利用Ｆｏｕｒｉｅｒ分析或小波理论，而地震波对建筑物动力 反应的研究则是通过拟加速度反应谱来进行．现行分析方法的数值化过程必须通过计算机编 程来实现，而采用常规的编程语言实现起来往往较为复杂。基于此，作者提出了基于ＭＡＴ— ＬＡＢ平台的分析方法，可以大大提高地震波动力特性分析的效率，便于分析结果的图形化输 出，为工程结构抗震计算提供有益的参考．
万方数据
第１期
马乐为等：基于ＭＡＴＩ，ＡＢ语言的地震波动力特性分析方法
・１３１・
Ｘ（ｓ）兰厶（￡）兰ｌ
ｅ『｜ｚ（￡）ｄｔ
ｋ（￡）一ｓＸ（ｓ）一ｚ（０） 疋（￡）一ｓ２ｘ（５）一耽（０）一主（Ｏ） Ｆ（ｓ）＝ｌ一厂（￡）ｄｔ
（２）
对于反应谱问题，由于外部激励，（￡）为质点质量ｍ与地震波加速度互。（￡）的乘积，因此质点的绝对加 速度口（￡）可近似表示为：
Ｅｎｄ
图中横坐标为时域的样本总量，乘以采样频率０．０１ Ｓ后即为ＳＨＷＩ的总持时１０．２５ ｓ；纵坐标为频率 分解的振幅，如果已知结构物的基本周期为０．２４ ｓ，频率为４．３７ Ｈｚ，则模型基频落在ａ４频带范围内，而图 中纵坐标表示了这一频率的幅值．由此可见，ＳＨＷｌ地震波在第２～６ Ｓ时段内的频率成分与结构模型相 近，并且这一时段内包含了地震波加速度峰值最大的时段，故在拟动力试验截取地震波时应首先考虑这一 时段． ４结束语 拟加速度反应谱在结构地震作用计算中普遍采用，绘制拟加速度反应谱曲线的常规方法是时域内的 Ｄｕｈａｍｅｌ积分。本文的研究表明，根据拉普拉斯变换原理，并利用ＭＡＴＬＡＢ提供的传递函数ｔｆ（）及时域 系统响应仿真函数ｌｓｉｍ（）可使反应谱曲线的求解过程大大简化． 参考文献
Ｎ
ｘ（忌）＝至：ｚ（咒）ｅ－ｚ州卜”（””／”
例
Ｎ－－Ｉ
（６）
ｓ（ｚ）一芝：Ｒ。ｅ－ｚｉ，ｏ“
式中，Ｎ为地震波加速度记录总数，取Ｎ＝２”，ｋ＝１，２，…，Ｎ／２，Ｒ为自相关函数． 由式（６）可知，对离散的地震动加速度信号进行快速傅里叶变换即可得到其线性谱Ｘ（愚），而离散信号 的自功率谱则可以表示为其线性谱乘其共轭线性谱．在Ｍａｔｌａｂ语言中ｆｉｔ（）是离散序列的Ｆｏｕｒｉｅｒ变换函 数，ｃｏｎｊ（）是复数的共轭函数，直接利用这两个函数便可迅速求解给定地震动加速度信号的功率谱密度． 具体函数如下：
谱密度均可以作为一种合理性的判别依据． 地震波谱密度分析的数学基础是Ｆｏｕｒｉｅｒ变换．当样本总量满足２的砣次幂时，可直接采用快速Ｆｏｕ— ｒｉｅｒ变换，不满足此条件时，也可用末位补零的方法进行快速Ｆｏｕｒｉｅｒ变换，而地震波加速度记录中某一频 率成分所占的比重则可用功率谱密度来反映．设地震波加速度记录为一个离散函数样本ｚ（行），其Ｆｏｕｒｉｅｒ 变换Ｘ（是）及功率谱密度Ｓ（ｚ）可按下式计算：
ｉ＝１：ｌｅｎｇｔｈ（ｏｍｇａｒｊ
ＳＹＳ＝ｔｆ（１，［１，２＊ｄｍｐｒ＊ｏｍｇａｒ（ｉ），ｏｍｇａｒ（ｉ）＊ｏｍｇａｒ（ｉ）］）；
ｄｉｓｐ（ｉ）＝ｍａｘ（ａｂｓ（１ｓｉｍ（ＳＹＳ，ｅｔｈｗｖ，ｔ）））； ｓｐｄ（ｉ）＝２＋ｐｉ．＊ｄｉｓｐ（ｉ）／Ｐｒｄ（ｉ）； ａｃｅ（ｉ）＝２＊ｐＬ＊ｓｐｄ（ｉ）／Ｐｒｄ（ｉ）；