el-centro地震波积分.doc

EL centro 地震波的积分及相关问题1.地震波的选取本作业的EL centro地震波来源于（NISEEthe national information service for earthquake engineering University of California, Berkeley）。具体链接为：/data/strong_motion/caltech/volume2.d/el_centro_1940_s00e这是一条南北向的地震波波，已经过一定的修正，网站中给出了地震波的初试速度和初始位移，原文具体说明如下：CORRECTED ACCELEROGRAM IIA001 40.001.0 COMP S00E FILE 1 CORRESPONDING TOFILE 1 OF UNCORRECTED ACCELEROGRAM DATA OF VOL. I-A, EERL 70-20 IMPERIAL VALLEY EARTHQUAKE MAY 18, 1940 - 2037 PST IA001 40.001.0 S 18 STATION NO. 117 32 47 43N,115 32 55W 38 EL CENTRO SITE IMPERIAL VALLEY IRRIGATION DISTRICT 50 COMP S00E 9 IMPERIAL VALLEY EARTHQUAKE MAY 18, 1940 - 2037 PST 52 EPICENTER 32 44 00N,115 27 00W 31 INSTR PERIOD = 0.0990 SEC DAMPING = 0.552 42 NO. OF POINTS = 985 DURATION = 53.73 SEC 42 UNITS ARE SEC AND G/10. 23 RMS ACCLN OF COMPLETE RECORD = 0.4876 G/10. 43 ACCELEROGRAM IS BAND-PASS FILTERED BETWEEN 0.070 AND 25.000 CYC/SEC 2688 INSTRUMENT AND BASELINE CORRECTED DATA AT EQUALLY-SPACED INTERVALS OF 0.02 SEC. PEAK ACCELERATION = 341.69531 CMS/SEC/SEC AT 2.1200 SEC PEAK VELOCITY = 33.44914 CMS/SEC AT 2.1800 SEC PEAK DISPLACEMENT = 10.86678 CMS AT 8.5800 SEC INITIAL VELOCITY = -4.66421 CMS/SEC INITIAL DISP. = 2.15852 CMS IMPERIAL VALLEY EARTHQUAKE MAY 18, 1940 - 2037 PST 读取excel中的EL CENTRO地震波南北向的数据，将时间数组命名为t，加速度数组命名为a；plot(t,a,r-)加速度波形图如下：图1. EL centro 地震波（加速度时程）图（单位）2.积分方法本次作业数据为离散的数据，下面将分别采用梯形公式以及辛普森公式进行数值积分。2.1梯形公式法将积分区间分为等份，记步长，节点为，由定积分性质，abf(x)abL1xdx=abx-ba-bfa+x-ab-a fbdx=12b-afa+fb (1)上面求积分公式的右端值可以看成是由线段x=a,x=b,过点(a,f(a),(b,f(b)的直线以及x轴围成的梯形公式。Matlab编程命令如下：v=zeros(length(t),1);v(1)=-46.6421;for i=2:length(v)v(i)=v(i-1)+1/2.*0.02.*(a(i-1)+a(i);endu=zeros(length(t),1);u(1)=21.5852;for i=2:length(u)u(i)=u(i-1)+1/2.*0.02.*(v(i-1)+v(i);endplot(t,v,r-);plot(t,u,r-);图2.梯形公式下的速度时程曲线图（单位）图3.梯形公式下的位移时程曲线图（单位）2.2辛普森公式若f(x)用通过节点x0=a，x1=a+b2，x2=b的二次插值多项式L2(x)代替，则fxL2x=x-x1x-x2x0-x1x0-x2fx0+x-x0x-x2x1-x0x1-x2fx1+x-x0x-x1x2-x0x2-x1fx2 (2)得到辛普森公式，或称抛物型公式 abfxdxabL2xdx=16(b-a)fa+4fa+b2+fb (3)具体做法是先将已有的原始数据进行拉格朗日二次插值得到两个时间坐标之间的一个数据，这样可以保证积分过程中不会因为积分方法而减少数据数量，然后再进行辛普森公式的积分计算，过程如下:x=zeros(length(t)-1,1); for i=1:length(x)x(i)=t(i+1)-0.01;endd=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); for i=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);endvsim=zeros(size(t),1);vsim(1)=-46.6421;for i=2:length(vsim)vsim(i)=vsim(i-1)+1/6.*0.02.*(a(i-1)+4.*d(i-1)+a(i);enddd=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);for i=2: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);endusim=zeros(size(t),1);usim(1)=21.5852;for i=2:length(usim)usim(i)=usim(i-1)+1/6.*0.02.*(vsim(i-1)+4.*dd(i-1)+vsim(i);endplot(t,vsim,r-)plot(t,usim,r-)图4.辛普森公式下的速度时程曲线图（单位）图5.辛普森公式下的位移时程曲线图（单位）2.3两种积分比较在matlab中将两种积分公式积出来的图形结果输出在一张图中，进而比较两者积分的相似程度。hold onplot(t,v,r-)plot(t,vsim,b-)hold off图6.速度时程曲线对比图（单位）hold onplot(t,u,r-)plot(t,usim,b-)hold off图7.位移时程曲线对比图（单位）注：红色表示梯形法，蓝色表示辛普森法。由此可见，两种积分方法在本作业中精度较为接近，在速度时程曲线对比图中两种积分方法基本没有区别，在位移时程曲线对比图中两种积分方法开始出现一定的区别，但仍在可接受范围内。由数值分析原理可知，梯形公式和辛普森公式的误差估计如下： ； (4) (5)通过对绝对误差的分析可知，辛普森公式积分精度要更高一点，所以接下来的讨论中采取辛普森公式积分的结果进行进一步的分析。3.零线漂移调整目前，地震反应分析所使用的地震波来源于真实地震动的数据采集和地震动的人工合成。地震动采集的数据大都是以加速度时程的形式给出，而速度和位移时程通常由加速度积分获得。如果不考虑地层断裂、裂缝等因素，在地震结束后，地面运动应该回归初始状态，即地面位移、速度、加速度应该归零。但是，由于地震动采集过程中存在低频仪器噪声、低频环境噪声、加速初始值和速度初始值及人为操作误差等诸多原因，可能导致由积分得到的速度和位移时程在终点时刻非零，即零线漂移。另一方面，在地震动的人工合成计算中，人们往往比较注重对加速度时程的研究，而速度和位移时程仅仅通过地震加速度时程简单的积分运算得到，这也导致速度和位移时程表现漂移现象。同时，在结构地震作用下的时程分析中，若采用相同的数值积分方法和积分步长，位移输入模型可获得比加速度输入模型更高的计算精度。因此，对于地震加速度时程曲线进行校正，以获得合理的位移时程曲线非常重要。对地震加速度的校正方法有两类：一是滤波，即将加速度记录中不合适的波频过滤掉；二是更改加速度的初始值，或者调整加速度的记录，使加速度积分后的位移时程为零。采用最小二乘法对加速度记录进行修正：用二次多项式对a进行修正，设修正后的加速度为a1,修正后的速度为v1，修正后的位移为u1： a1(t)=a(t)-(b+ct+dt2)v1t=vt-(bt+12ct2+13dt3)u1t=ut-(12bt2+16ct3+112dt4) (6)在地震结束后，速度和位移应均为0，根据公式（6）可以得到： vT-bT+12cT2+13dT3=0 (7) uT-12bT2+16cT3+112dT4=0 (8)由式（7）和式（8）将c和b用d的代数式表示出来，再带入下式： min=0Ts(t)-(12bt2+16ct3+112dt4)2dt (9)根据最小二乘法曲线拟合可以得到:b=0.249497492684421c=-0.010571559564482d=-6.58801118748905210-5然后对加速度进行修正并重新进行积分计算，具体过程如下：for i=1:length(a1)a1(i)=a(i)-(xb+xc.*t(i)+xz.*(t(i).2);endd1=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);for i=2: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);endvsim1=zeros(size(t),1);vsim1(1)=-46.6421;for i=2:length(vsim1)vsim1(i)=vsim1(i-1)+1/6.*0.02.*(a1(i-1)+4.*d1(i-1)+a1(i);enddd1=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);for i=2: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);endusim1=zeros(size(t),1); usim1(1)=21.5852;for i=2:length(usim1)usim1(