elcentro地震波积分文档格式.docx
《elcentro地震波积分文档格式.docx》由会员分享，可在线阅读，更多相关《elcentro地震波积分文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
COMPS00E9
IMPERIALVALLEYEARTHQUAKEMAY18,1940-2037PST52
EPICENTER324400N,1152700W31
INSTRPERIOD=0.0990SECDAMPING=0.55242
NO.OFPOINTS=985DURATION=53.73SEC42
UNITSARESECANDG/10.23
RMSACCLNOFCOMPLETERECORD=0.4876G/10.43
ACCELEROGRAMISBAND-PASSFILTEREDBETWEEN0.070AND25.000CYC/SEC
2688INSTRUMENTANDBASELINECORRECTEDDATA
ATEQUALLY-SPACEDINTERVALSOF0.02SEC.
PEAKACCELERATION=341.69531CMS/SEC/SECAT2.1200SEC
PEAKVELOCITY=33.44914CMS/SECAT2.1800SEC
PEAKDISPLACEMENT=10.86678CMSAT8.5800SEC
INITIALVELOCITY=-4.66421CMS/SECINITIALDISP.=2.15852CMS
IMPERIALVALLEYEARTHQUAKEMAY18,1940-2037PST
读取excel中的ELCENTRO地震波南北向的数据，将时间数组命名为t，加速度数组命名为a；
plot（t,a,'
r-'
）
加速度波形图如下：
图1.ELcentro地震波（加速度时程）图（单位
2.积分方法
本次作业数据为离散的数据，下面将分别采用梯形公式以及辛普森公式进行数值积分。
2.1梯形公式法
将积分区间
分为
等份，记步长
，节点为
，由定积分性质，
（1）
上面求积分公式的右端值可以看成是由线段x=a,x=b,过点（a,f（a））,（b,f（b））的直线以及x轴围成的梯形公式。
Matlab编程命令如下：
v=zeros（length（t）,1）;
v
（1）=-46.6421;
fori=2:
length（v）
v（i）=v（i-1）+1/2.*0.02.*（a（i-1）+a（i））;
end
u=zeros（length（t）,1）;
u
（1）=21.5852;
length（u）
u（i）=u（i-1）+1/2.*0.02.*（v（i-1）+v（i））;
plot（t,v,'
）;
plot（t,u,'
图2.梯形公式下的速度时程曲线图（单位
图3.梯形公式下的位移时程曲线图（单位
2.2辛普森公式
若f（x）用通过节点
，
的二次插值多项式
代替，则
（2）
得到辛普森公式，或称抛物型公式
（3）
具体做法是先将已有的原始数据进行拉格朗日二次插值得到两个时间坐标之间的一个数据，这样可以保证积分过程中不会因为积分方法而减少数据数量，然后再进行辛普森公式的积分计算，过程如下:
x=zeros（length（t）-1,1）;
fori=1:
length（x）
x（i）=t（i+1）-0.01;
d=zeros（length（t）-1,1）;
d
（1）=（x
（1）-t
（2））.*（x
（1）-t（3））./（（t
（1）-t
（2））.*（t
（1）-t（3）））.*a
（1）+（x
（1）-t
（1））.*（x
（1）-t（3））./（（t
（2）-t
（1））.*（t
（2）-t（3）））.*a
（2）+（x
（1）-t
（1））.*（x
（1）-t
（2））./（（t（3）-t
（1））.*（t（3）-t
（2）））.*a（3）;
fori=2:
size（d）
d（i）=（x（i）-t（i））.*（x（i）-t（i+1））./（（t（i-1）-t（i））.*（t（i-1）-t（i+1）））.*a（i-1）+（x（i）-t（i-1））.*（x（i）-t（i+1））./（（t（i）-t（i-1））.*（t（i）-t（i+1）））.*a（i）+（x（i）-t（i-1））.*（x（i）-t（i））./（（t（i+1）-t（i-1））.*（t（i+1）-t（i）））.*a（i+1）;
vsim=zeros（size（t）,1）;
vsim
（1）=-46.6421;
length（vsim）
vsim（i）=vsim（i-1）+1/6.*0.02.*（a（i-1）+4.*d（i-1）+a（i））;
dd=zeros（length（t）-1,1）;
dd
（1）=（x
（1）-t
（2））.*（x
（1）-t（3））./（（t
（1）-t
（2））.*（t
（1）-t（3）））.*vsim
（1）+（x
（1）-t
（1））.*（x
（1）-t（3））./（（t
（2）-t
（1））.*（t
（2）-t（3）））.*vsim
（2）+（x
（1）-t
（1））.*（x
（1）-t
（2））./（（t（3）-t
（1））.*（t（3）-t
（2）））.*vsim（3）;
size（dd）
dd（i）=（x（i）-t（i））.*（x（i）-t（i+1））./（（t（i-1）-t（i））.*（t（i-1）-t（i+1）））.*vsim（i-1）+（x（i）-t（i-1））.*（x（i）-t（i+1））./（（t（i）-t（i-1））.*（t（i）-t（i+1）））.*vsim（i）+（x（i）-t（i-1））.*（x（i）-t（i））./（（t（i+1）-t（i-1））.*（t（i+1）-t（i）））.*vsim（i+1）;
usim=zeros（size（t）,1）;
usim
（1）=21.5852;
length（usim）
usim（i）=usim（i-1）+1/6.*0.02.*（vsim（i-1）+4.*dd（i-1）+vsim（i））;
plot（t,vsim,’r-‘）
plot（t,usim,’r-‘）
图4.辛普森公式下的速度时程曲线图（单位
图5.辛普森公式下的位移时程曲线图（单位
2.3两种积分比较
在matlab中将两种积分公式积出来的图形结果输出在一张图中，进而比较两者积分的相似程度。
holdon
plot（t,vsim,'
b-'
holdoff
图6.速度时程曲线对比图（单位
plot（t,usim,'
图7.位移时程曲线对比图（单位
注：
红色表示梯形法，蓝色表示辛普森法。
由此可见，两种积分方法在本作业中精度较为接近，在速度时程曲线对比图中两种积分方法基本没有区别，在位移时程曲线对比图中两种积分方法开始出现一定的区别，但仍在可接受范围内。
由数值分析原理可知，梯形公式和辛普森公式的误差估计如下：
；
（4）
（5）
通过对绝对误差的分析可知，辛普森公式积分精度要更高一点，所以接下来的讨论中采取辛普森公式积分的结果进行进一步的分析。
3.零线漂移调整
目前，地震反应分析所使用的地震波来源于真实地震动的数据采集和地震动的人工合成。
地震动采集的数据大都是以加速度时程的形式给出，而速度和位移时程通常由加速度积分获得。
如果不考虑地层断裂、裂缝等因素，在地震结束后，地面运动应该回归初始状态，即地面位移、速度、加速度应该归零。
但是，由于地震动采集过程中存在低频仪器噪声、低频环境噪声、加速初始值和速度初始值及人为操作误差等诸多原因，可能导致由积分得到的速度和位移时程在终点时刻非零，即零线漂移。
另一方面，在地震动的人工合成计算中，人们往往比较注重对加速度时程的研究，而速度和位移时程仅仅通过地震加速度时程简单的积分运算得到，这也导致速度和位移时程表现漂移现象。
同时，在结构地震作用下的时程分析中，若采用相同的数值积分方法和积分步长，位移输入模型可获得比加速度输入模型更高的计算精度。
因此，对于地震加速度时程曲线进行校正，以获得合理的位移时程曲线非常重要。
对地震加速度的校正方法有两类：
一是滤波，即将加速度记录中不合适的波频过滤掉；
二是更改加速度的初始值，或者调整加速度的记录，使加速度积分后的位移时程为零。
采用最小二乘法对加速度记录进行修正：
用二次多项式对a进行修正，设修正后的加速度为a1,修正后的速度为v1，修正后的位移为u1：
在地震结束后，速度和位移应均为0，根据公式（6）可以得到：
由式（7）和式（8）将c和b用d的代数式表示出来，再带入下式：
根据最小二乘法曲线拟合可以得到:
然后对加速度进行修正并重新进行积分计算，具体过程如下：
fori=1:
length（a1）
a1（i）=a（i）-（xb+xc.*t（i）+xz.*（t（i）.^2））;
d1=zeros（length（t）-1,1）;
d1
（1）=（x
（1）-t
（2））.*（x
（1）-t（3））./（（t
（1）-t
（2））.*（t
（1）-t（3）））.*a1
（1）+（x
（1）-t
（1））.*（x
（1）-t（3））./（（t
（2）-t
（1））.*（t
（2）-t（3）））.*a1
（2）+（x
（1）-t
（1））.*（x
（1）-t
（2））./（（t（3）-t
（1））.*（t（3）-t
（2）））.*a1（3）;
size（d1）
d1（i）=（x（i）-t（i））.*（x（i）-t（i+1））./（（t（i-1）-t（i））.*（t（i-1）-t（i+1）））.*a1（i-1）+（x（i）-t（i-1））.*（x（i）-t（i+1））./（（t（i）-t（i-1））.*（t（i）-t（i+1）））.*a1（i）+（x（i）-t（i-1））.*（x（i）-t（i））./（（t（i+1）-t（i-1））.*（t（i+1）-t（i）））.*a1（i+1）;
vsim1=zeros（size（t）,1）;
vsim1
（1）=-46.6421;
length（vsim1）
vsim1（i）=vsim1（i-1）+1/6.*0.02.*（a1（i-1）+4.*d1（i-1）+a1（i））;
dd1=zeros（length（t）-1,1）;
dd1
（1）=（x
（1）-t
（2））.*（x
（1）-t（3））./（（t
（1）-t
（2））.*（t
（1）-t（3）））.*vsim1
（1）+（x
（1）-t
（1））.*（x
（1）-t（3））./（（t
（2）-t
（1））.*（t
（2）-t（3）））.*vsim1
（2）+（x
（1）-t
（1））.*（x
（1）-t
（2））./（（t（3）-t
（1））.*（t（3）-t
（2）））.*vsim1（3）;
size（dd1）
dd1（i）=（x（i）-t（i））.*（x（i）-t（i+1））./（（t（i-1）-t（i））.*（t（i-1）-t（i+1）））.*vsim1（i-1）+（x（i）-t（i-1））.*（x（i）-t（i+1））./（（t（i）-t（i-1））.*（t（i）-t（i+1）））.*vsim1（i）+（x（i）-t（i-1））.*（x（i）-t（i））./（（t（i+1）-t（i-1））.*（t（i+1）-t（i）））.*vsim1（i+1）;
usim1=zeros（size（t）,1）;
usim1
（1）=21.5852;
length（usim1）
usim1（i）=usim1（i-1）+1/6.*0.02.*（vsim1（i-1）+4.*dd1（i-1）+vsim1（i））;
plot（t,a1,’r-’）
plot（t,vsim1,’r-’）
plot（t,usim1,’r-’）
图8.ELcentro地震波（修正后加速度时程）图（单位
图9.修正后速度时程曲线图（单位
图10.修正后位移时程曲线图（单位
4．总结
通过对比修正前后加速度时程曲线对比图（图11）可以发现，加速度峰值基本未发生较大变动，且修正前后的图形几乎没有差异，故此方法对于加速度修正可行。
通过对比修正前后速度时程曲线对比图（图12）可以发现，地震结束后位移基本归零，且修正前后的图形几乎没有差异，故此方法可行。
通过对比修正前后位移时程曲线对比图（图13）可以发现，地震波位移以及速度的零位漂移现象基本得到了解决，故此方法可行。
图11.加速度时程曲线对比图（单位
图12.速度时程曲线对比图（单位
红色线表示修正前，蓝色线表示修正后。
参考文献
[1].同济大学计算数学教研室,现代数值计算[M],人民邮电出版社,2012.
[2].李洁涛,杨庆山,地震波基线漂移的处理方法[J],北京交通大学学报,2010,34
（1）.
[3].冯仲仁，李文俊，地震波积分计算与基调方法，中国科技论文在线[2008-03-03]
升级会员 