基于快速傅里叶变换的地震波加速度构成及其幅频特性研究*

付 荣 傅荣华 付安生

1) 中国山西晋中030620山西省第三地质工程勘察院 2) 中国成都610059成都理工大学环境与土木工程学院

基于快速傅里叶变换的地震波加速度构成及其幅频特性研究*

1) 中国山西晋中030620山西省第三地质工程勘察院 2) 中国成都610059成都理工大学环境与土木工程学院

根据2008年汶川MS8.0地震广元曾家台(GYZJ)的强震记录,运用Matlab软件对东西向55—75 s的加速度时程进行快速傅里叶变换,得到有限复傅里叶系数,从而将离散的加速度记录转化为加速度随时间变化的函数表达式． 该表达式中包含了50组不同频率的正弦波和余弦波,由此确定了地震波加速度的主要构成． 经回归分析可知, 函数计算值与原始记录的相关系数R为0.9,拟合程度良好； 利用SeismoSignal软件对地震波加速度时程进行谱分析,结果表明在原来的波中频率约为2.8 Hz的谐振分量振幅最大且0—10 Hz这一频段的波最为丰富. 这一结果与函数表达式中频率与振幅的分布基本吻合, 为深入研究汶川地震的动力特性及其对山体动力效应的分析提供了重要的理论依据,具有广阔的应用前景．

强震记录 加速度 快速傅里叶变换 幅频特性

引言

2008年汶川MS8.0地震诱发了数以万计的次生地质灾害,大量斜坡的变形破坏反映出许多动力特征, 如青川东河口特大型滑坡的临空抛射现象. 斜坡在地震作用下的动力响应特征与斜坡所遭受的地震动特性(振幅、 频谱、 持续时间)密切相关(祁生文,2006)． 因此,深入研究强震记录中地震波振幅与频率的组合特征及其动力输入,是研究地震作用下斜坡动力响应问题的一个重要方面．

汶川地震后,布设在龙门山断裂带附近的强震观测仪获得了大量的强震记录,这些记录为研究地震诱发地质灾害与强震记录地震动参数之间的关系提供了宝贵的数据资料． 研究人员利用强震记录进行了多方面研究(王秀英等,2009； 于海英等,2009； 刘晓,2010； 周荣军等,2010； 张冬丽等,2010； 郭晓云,2011； 罗永红,2011; 谢俊举等,2011)． 本文对广元曾家台(GYZJ)获得的汶川地震加速度记录进行快速傅里叶变换,研究该台记录到的地震波加速度的主要构成及其幅频特性,旨在为研究汶川地震的地震波特性及其对斜坡的影响提供理论依据．

1 强震记录

广元曾家台观测场地条件为土层,仪器型号为MR2002,位于汶川地震发震断层下盘(图1),震中距为311.38 km(中国地震局震害防御司,2008)．

图2给出了广元曾家台在汶川地震时记录的未校正加速度时程(东西向),数据记录点为1万9356个,采样时间间隔为0.01 s,最大峰值加速度为4.24 m/s2．

地震动能量与加速度、 速度或位移幅值的平方成正比(谢礼立,周雍年,1984)． 图3为广元曾家台接收到的汶川地震加速度能量在时域上的分布图． 从图3中可以看出55—75 s是该记录能量最强的时间段． 为避免表达式冗长,对地震波加速度时程进行了截取,选择55—75 s时间段内东西向的加速度记录作为有限傅里叶近似的加速度数据(图4),数据点为1001个,时间间隔为0.02 s．

2 地震波加速度时程的有限傅里叶近似

2.1 有限傅里叶近似

广元曾家台记录的地震波加速度是按等时间间隔读取的离散采样点,本文采用大崎顺彦(1980)的有限傅里叶近似方法,用离散的采样值求取原光滑曲线的表达式． 其近似表达式为

图1 地震台站分布与2008年汶川MS8.0地震地质灾害区划的对应关系(据王秀英等 (2010)修改)

图2 广元曾家台(GYZJ)记录到的汶川地震未校正加速度时程(东西向)

图3 广元曾家台接收到的汶川地震 加速度能量在时域上的分布图

图4 广元曾家台记录的55—75 s时段 汶川地震东西向加速度时程

(1)

其中，

(2)

式中,Ak和Bk为傅里叶系数,m为采样点序号(m=0, 1, 2, …,N-1),xm为采样值,N为采样数,Δt为采样点的时间间隔．

利用广元曾家台记录的加速度随时间变化的离散值,进行有限傅里叶近似． 令k=0,则

(3)

显然,A0代表的是全体采样值的平均值,即整个波形对零线的漂移． 本文对广元曾家台在截取时段55—75s记录的1001个加速度采样值求平均得0.14,所以东西向A0/2=0.14．

从式(1)可以看出,原来的波形表现为正弦与余弦的集合． 因此,可把原来的波形看成已被分解的正弦波和余弦波．

若令

(4)

则式(1)可写成

(5)

由式(5)可知,原来的振动被分解成(N-1)/2个不同频率的波,即f1,f2, …,f(N-1)/2. 这里称f1为1次频率(基本频率),f2为2次频率,一般项fk为k次频率,并把k次频率的振动波称为原来波形的k次谐振分量．

2.2 有限傅里叶系数

欲求波形的有限傅里叶系数,按照式(2)作傅里叶变换是最简单明了的． 但是,该方法计算过于费时,而且随着采样数N的增大,计算时间将以N2的比例增加． 为了克服这个缺点,Gooley与Tukey研究出一种非常快速的傅里叶分析方法,称为Gooley-Tukey法或快速傅里叶变换(fast Fourier transform,简写为FFT)(大崎顺彦,1980). 该方法是对等间隔采样点的N个采样值xm(m=0, 1, 2, …,N-1),求其有限复傅里叶系数Ck(k=0, 1, 2, …,N-1)． 采用这种傅里叶变换方法计算速度非常快,所需要的计算时间不再与数据个数的平方成比例,而是当N=2p时,与pN成比例(大崎顺彦,1980)． 这为将采样点极多的地震数据转化为函数表达式,从而研究地震波的构成提供了方便快捷的方法．

本文利用Matlab软件对截取时段的采样值进行快速傅里叶变换,得到有限复傅里叶系数Ck,根据下式(大崎顺彦,1980)求得Ak和Bk：

(6)

式中,Re(·)表示实部,Im(·)表示虚部．

按k次分量的振幅由大到小排序,得到前50个正弦波与余弦波的组合(表1)．

表1 广元曾家台加速度记录的有限复傅里叶系数

图5 广元曾家台在截取时段上的原始加速度 时程与拟合函数曲线的对比图 Fig.5 Comparison of uncorrected acceleration observation curve (bule) and fitted function curve (red) recorded at GYZJ station in the time period 55—75 s

为观察函数表达式的拟合效果,将表1中的系数Ak和Bk及其所对应的频率fk代入式(5),得到加速度随时间变化的拟合函数； 再将55—75 s(Δt=0.02 s)的时间值代入该拟合函数,即可得到计算后的加速度值. 加速度拟合曲线与原始曲线的对比可以看出拟合结果与原始数据基本吻合(图5)．

为进一步检验拟合效果,利用Matlab统计工具箱中的regress函数命令计算实测值x与公式计算值Y之间的线性回归方程,并将其残差的置信区间不包括零点的异常数据剔除,最终得Y=-0.4+0.8x,相关系数R=0.9,拟合效果良好． 其残差及置信区间见图6．

从图6可以看出,除少数数据(红线及圆圈所对应的数据)外,其余数据的残差值(蓝点)均较接近零点,且残差的置信区间(绿线)均包含零点. 这说明回归模型能较好地拟合原始数据,即计算值与实测值比较接近．

图6 广元曾家台记录的原始加速度值与拟合函数值的残差分布图

3 地震波的谱分析

谱分析为研究地震波频率的构成提供了一种有效的手段(刘晓,2010)． 利用SeismoSignal软件(SeismoSoft Ltd, 2013)即可得到傅里叶幅值谱和功率谱． 从图7的地震波傅里叶幅值谱中可以看出地震波的主频率值约为2.8 Hz． 从图8的傅里叶功率谱可以看出广元曾家台所获得的东西向强震记录的能量主要集中在10 Hz以内的低频范围内． 这一结果与表1中频率fk和振幅Xk的分布情况基本吻合,从另一方面检验了函数表达式

的拟合效果,其中,Ak,Bk和fk由表1给出．

图7 广元曾家台加速度记录的傅里叶幅值谱

图8 广元曾家台加速度记录的傅里叶功率谱

4 讨论与结论

为研究地震波加速度的主要构成,本文对广元曾家台(GYZJ)记录到的汶川MS8.0地震的强震记录进行快速傅里叶变换,得到有限复傅里叶系数,从而将离散的加速度记录转化为加速度随时间变化的函数表达式,并利用SeismoSignal软件对地震波加速度时程进行了谱分析,得到以下结论：

1) 广元曾家台在汶川地震时监测到的东西向地震波(55—75 s)主要由表1中所列的50项不同频率的正弦波和余弦波构成． 通过分析可知,所得结果与原始数据拟合良好,表明利用快速傅里叶变换从强震记录中获取地震波加速度振幅和频率的构成特征是一个很好的方法,为深入研究汶川地震的动力特性及其对山体动力效应的分析提供了理论依据,具有广阔的应用前景．

2) 通过谱分析可得,广元曾家台记录的东西向地震波中频率约为2.8 Hz的谐振分量振幅最大且0—10 Hz这一频段的波最为丰富． 这一结果与表1中频率fk和振幅Xk的分布情况基本吻合,从另一方面检验了函数表达式的适用性．

中国强震动观测台网为本研究提供了强震记录数据,审稿专家为本文的完善提出了许多宝贵的意见,在此一并表示感谢．

大崎顺彦(著). 1980. 吕敏申,谢礼立(译). 2008. 地震动的谱分析入门[M]. 北京: 地震出版社: 1--297.

Osaki Y. 1980. Lü M S, Xie L L (trans). 2008.IntroductiontoSpectrumAnalysisofGroundMotion[M]. Beijing: Seismological Press: 1--297 (in Chinese)．

郭晓云. 2011. 汶川地震反应谱研究[D]. 哈尔滨： 中国地震局工程力学研究所: 1--135.

Guo X Y. 2011.StudyonResponseSpectrumofWenchuanEarthquake[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration: 1--135 (in Chinese)．

刘晓. 2010. 汶川地震区斜坡动力反应研究[D]. 武汉： 中国地质大学: 1--155.

Liu X. 2010.ResearchonDynamicResponseofSlopeinWenchuanEarthquakeArea[D]. Wuhan: China University of Geosciences: 1--155 (in Chinese)．

罗永红. 2011. 地震作用下复杂斜坡响应规律研究[D]. 成都： 成都理工大学: 1--211.

Luo Y H. 2011.StudyonComplexSlopesResponseLawUnderEarthquakeAction[D]. Chengdu: Chengdu University of Technology: 1--211 (in Chinese)．

祁生文. 2006. 单面边坡的两种动力反应形式及其临界高度[J]. 地球物理学报,49(2)： 518--523.

Qi S W. 2006. Two patterns of dynamic responses of single-free-surface slopes and their threshold height[J].ChineseJournalofGeophysics,49(2)： 518--523 (in Chinese).

王秀英,聂高众,王登伟. 2009. 利用强震记录分析汶川地震诱发滑坡[J]. 岩石力学与工程学报,28(11)： 2369--2376.

Wang X Y,Nie G Z,Wang D W. 2009. Analysis of landslide induced by Wenchuan earthquake by strong motion records[J].ChineseJournalofRockMechanicsandEngineering,28(11)： 2369--2376 (in Chinese)．

王秀英,聂高众,王登伟. 2010. 汶川地震诱发滑坡与地震动峰值加速度对应关系研究[J]. 岩石力学与工程学报,29(1)： 82--89.

Wang X Y,Nie G Z,Wang D W. 2010. Research on relationship between landslides and peak ground accelerations induced by Wenchuan earthquake[J].ChineseJournalofRockMechanicsandEngineering,29(1)： 82--89 (in Chinese).

谢俊举,温增平,高孟潭,袁美巧,何少林. 2011. 2008年汶川地震近断层地震动的非平稳特征[J]. 地球物理学报,54(3)： 728--736.

Xie J J,Wen Z P,Gao M T,Yuan M Q,He S L. 2011. Non-stationary characteristics of near-fault strong motions during the 2008 Wenchuan earthquake[J].ChineseJournalofGeophysics,54(3)： 728--736 (in Chinese)．

谢礼立,周雍年. 1984. 一个新的地震动持续时间定义[J]. 地震工程与工程振动,4(2)： 27--35.

Xie L L,Zhou Y N. 1984. A new definition of strong ground motion duration[J].JournalofEarthquakeEngineeringandEngineeringVibration,4(2)： 27--35 (in Chinese).

于海英,王栋,杨永强,解全才,江汶乡,周宝峰. 2009. 汶川8.0级地震强震动加速度记录的初步分析[J]. 地震工程与工程振动,29(1)： 1--13.

Yu H Y,Wang D,Yang Y Q,Xie Q C,Jiang W X,Zhou B F. 2009. The preliminary analysis of strong ground motion records from theMS8.0 Wenchuan earthquake[J].JournalofEarthquakeEngineeringandEngineeringVibration,29(1)： 1--13 (in Chinese)．

张冬丽,徐锡伟,周正华,于贵华. 2010. 汶川MS8.0级地震强震记录所揭示的地震断层特征分析[J]. 地学前缘,17(5)： 33--42.

Zhang D L, Xu X W, Zhou Z H, Yu G H. 2010. The causative fault properties from strong-motion records ofMS8.0 Wenchuan earthquake[J].EarthScienceFrontiers, 17(5): 33--42 (in Chinese).

中国地震局震害防御司. 2008. 汶川8.0级地震未校正加速度记录[M]. 北京： 地震出版社: 1--302.

Department of Earthquake Disaster Prevention, China Earthquake Administration. 2008.UncorrectedAccelerationRecordsoftheWenchuanMS8.0Earthquake[M]. Beijing： Seismological Press: 1--302 (in Chinese).

周荣军,赖敏,余桦,龙承厚,王世元,亢川川,梁明剑. 2010. 汶川MS8.0地震四川及邻区数字强震台网记录[J]. 岩石力学与工程学报,29(9)： 1850--1858.

Zhou R J,Lai M,Yu H,Long C H,Wang S Y,Kang C C,Liang M J. 2010. Strong motion records of WenchuanMS8.0 earthquake from digital strong earthquake networks in Sichuan and its neighbouring regions[J].ChineseJournalofRockMechanicsandEngineering, 29(9): 1850--1858 (in Chinese)．

SeismoSoft Ltd. 2013. SeismoSignal_v5.1.0[CP/OL]. [2013-09-20]. http:∥www.seismosoft.com/en/download_details.aspx?ID_Download=1.

Composition and amplitude-frequency characteristics of ground motion acceleration based on fast Fourier transform analysis

1)TheThirdInstituteofGeologicalEngineeringReconnaissance,ShanxiJinzhong030620,China2)CollegeofEnvironmentandCivilEngineering,ChengduUniversityofTechnology,Chengdu610059,China

This paper analyzed the strong ground motion records of the 2008 WenchuanMS8.0 earthquake at the station GYZJ, and used Matlab software to carry on the fast Fourier transform on the acceleration records of E--W component in 55—75 s, and then obtained a function expression of the acceleration variation with time. The expression contains 50 groups of sine and cosine waves with different frequencies, which are the main composition of the seismic wave. A correlation coefficient of 0.9 deduced from regression analysis shows well fitting between the calculated values and observed records. In addition, we used the SeismoSignal software to analyze the spectrum of the ground motion acceleration. Our results show that the harmonic component with frequency of 2.8 Hz in the original wave has the maximum amplitude and a dominant range of seismic waves has a frequency band of 0—10 Hz. These results generally agree well with the frequency and amplitude distribution of the functional expression. This paper provides an important theoretical basis for the research of wave characteristics and the analysis on the dynamic effect of slope in Wenchuan earthquake.

strong ground motion record； acceleration； fast Fourier transform； amplitude-frequency characteristic

10.3969/j.issn.0253-3782.2014.03.007.

国家自然科学基金(41072228)资助.

2013-04-18收到初稿,2013-09-02决定采用修改稿.

e-mail: fu-rong20080808@163.com

10.3969/j.issn.0253-3782.2014.03.007

P315.3+1

A

付荣, 傅荣华, 付安生. 2014. 基于快速傅里叶变换的地震波加速度构成及其幅频特性研究. 地震学报, 36(3): 417--424.

Fu R, Fu R H, Fu A S. 2014. Composition and amplitude-frequency characteristics of ground motion acceleration based on fast Fourier transform analysis.ActaSeismologicaSinica, 36(3): 417--424. doi:10.3969/j.issn.0253-3782.2014.03.007.