地震工程作业哈工大 MATLAB程序

Matlab在地震数据处理与分析中的应用指南

Matlab在地震数据处理与分析中的应用指南地震是一种自然灾害，对人们的生命和财产安全造成了巨大威胁。

了解地震的发生和传播规律，对于地震风险评估、灾害预警和防御措施的制定都具有重要意义。

然而，地震数据的处理和分析是一项复杂而繁琐的工作。

在这个过程中，Matlab作为一种功能强大、易于使用的数学建模软件，可以帮助地震学家和研究人员高效地进行地震数据的处理和分析。

本文将介绍Matlab在地震数据处理与分析中的应用指南，以帮助读者更好地运用Matlab进行相关工作。

一、地震数据的读取与可视化处理地震数据通常以数值形式存储在地震波形文件中，这些文件的格式各不相同。

Matlab提供了丰富的函数库，可以读取多种地震数据文件格式，并将其转换为方便处理的矩阵数据。

以SAC文件为例，可以使用sacread函数读取SAC文件，并将其转换为Matlab中的矩阵数据。

读取地震数据后，我们可以使用Matlab强大的图形绘制功能，对地震波形进行可视化处理，更直观地了解地震数据的特征。

Matlab的plot函数可以绘制地震波形的时间序列曲线，利用subplot函数可以将多个波形图像进行排列，方便对比不同地震事件。

二、地震波形的滤波与去噪处理地震数据中通常包含大量的噪声干扰，这些噪声对于地震数据的分析和解释会产生不利影响。

Matlab提供了一系列信号滤波函数，可以有效地去除地震数据中的噪声。

常用的滤波方法包括低通滤波、高通滤波和带通滤波等。

我们可以根据地震波形的频率特征选择适当的滤波方法，并利用Matlab的filter函数进行滤波处理。

此外，Matlab还提供了多种经典的去噪算法，如小波变换去噪、小波阈值去噪等，这些方法可以更精确地去除地震波形中的噪声成分。

三、地震数据的频率域分析地震波形的频率域分析是对地震数据进行深入研究和理解的重要手段。

Matlab 提供了丰富的频率域分析函数，可以计算地震波形的功率谱密度、相位谱、互相关等频域特征。

matlab地震反应谱

matlab地震反应谱摘要：I.引言- 介绍Matlab 在地震工程领域的应用- 介绍本文的主要内容II.地震反应谱的计算方法- 加载地震波数据- 预处理地震波数据- 计算地震波数据的频谱- 计算地震反应谱III.地震反应谱的分析方法- 评估结构的地震响应特性- 为地震工程设计提供依据IV.结论- 总结Matlab 在地震反应谱计算中的应用- 展望未来可能的研究方向正文：Matlab 是一种功能强大的数学软件，可以用于解决各种工程和科学问题。

在地震工程领域，Matlab 可以用于计算地震反应谱，从而评估结构的地震响应。

本文将介绍如何使用Matlab 计算地震反应谱，以及如何根据反应谱分析结构的地震响应。

首先，我们需要加载地震波数据。

可以使用Matlab 内置的函数fiducial() 从文件中读取地震波数据，也可以使用自定义函数从其他来源获取数据。

接下来，我们需要对地震波数据进行预处理，以去除噪声和异常值。

可以使用Matlab 内置的函数removeoutliers() 对数据进行去噪处理。

然后，我们可以使用Matlab 内置的函数freqz() 计算地震波数据的频谱。

通过频谱，我们可以了解地震波的频率分量以及振幅。

最后，我们可以使用Matlab 内置的函数responsespectrum() 计算地震反应谱。

该函数可以根据结构的阻尼比和自振周期计算反应谱。

通过反应谱，我们可以了解结构在不同地震波作用下的加速度响应。

综上所述，Matlab 可以用于计算地震反应谱，从而评估结构的地震响应。

matlab计算单自由度的地震反应的程序

采用EL-CENTRO地震波计算单自由体系的振动位移反应谱前言：本课程论文采用软件是MATLAB R2010a，所编程序应用的方法是线性加速度法。

程序采用的是两层循环上实现位移谱的求取。

内循环实现在地震波下由上一步所得到的位移、速度和加速度得到下一时间间隔后的位移速度和加速度；外循环实现在不同结构自振周期（频率）下遭遇地震波的响应。

本文在求出相对位移反应谱的同时，也求出了相对速度反应谱、绝对加速度反应谱。

一．线性加速度法的简述：线性加速度法是直接数值积分法求解地震反应的方法之一，本文所采用的线性加速度法参考大崎顺彦的《地震动的谱分析入门》第二版。

具体计算公式详见大崎顺彦的《地震动的谱分析入门》第二版P116-P118。

二．所编程序及编译：clear% ***********读入地震记录***********fid = fopen('ei.txt');[Accelerate,count] = fscanf(fid,'%g'); %count 读入的记录的量time=0:0.02:(count-1)*0.02;% ***********线性加速度法计算各反应***********%初始化各储存向量Displace=zeros(1,count); %相对位移Velocity=zeros(1,count); %相对速度AbsAcce=zeros(1,count);%绝对加速度Damp=0.05; %结构阻尼比取为0.05Tc=0.0:0.05:10; %结构自振周期Dt=0.02; %地震记录的步长%记录计算得到的反应，MDis为最大相对位移，MVel为最大相对速度%MAcc某阻尼时最大绝对加速度，用于画图MDis=zeros(1,length(Tc));MVel=zeros(1,length(Tc));MAcc=zeros(1,length(Tc));t=1; %在下一个循环中控制不同的结构自振周期for T=0.0:0.05:10Frcy=2*pi/T ; %结构自振频率DamFrcy=Frcy*sqrt(1-Damp*Damp);%计算公式化简e_t=exp(-Damp*Frcy*Dt);s=sin(DamFrcy*Dt);c=cos(DamFrcy*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DamFrcy;A(2,1)=-Frcy*e_t*s/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp^2-1)/(Frcy^2*Dt); %计算公式化简d_3t=Damp/(Frcy^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/Frcy)*s/DamFrcy+(2*d_3t+1/Frcy^2)*c)-2*d_3t;B(1,2)=-e_t*(d_f*s/DamFrcy+2*d_3t*c)-1/Frcy^2+2*d_3t;B(2,1)=-e_t*(((Damp/(Frcy*Dt)+1)*s/DamFrcy)+(1/(Frcy^2*Dt))*c)+1/(Frcy^2*Dt);B(2,2)=e_t*((Damp/(Frcy*Dt)*s/DamFrcy)+(1/(Frcy^2*Dt))*c)-1/(Frcy^2*Dt);for i=1:(count-1) %根据地震记录,计算不同的反应Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i +1);Velocity(i+1)=A(2,1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i+ 1);AbsAcce(i+1)=-2*Damp*Frcy*Velocity(i+1)-Frcy^2*Displace(i+1);endMDis(1,t)=max(abs(Displace));MVel(1,t)=max(abs(Velocity));if T==0.0MAcc(1,t)=max(abs(Accelerate)); %当结构的自振周期为0时，其绝对加速度应等于地面加速度elseMAcc(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,count);%初始化各储存向量避免下次计算时引用到前面结果Velocity=zeros(1,count);AbsAcce=zeros(1,count);t=t+1; %t=length(Tc)，即所求结构自振周期有多少个，对应就运行多少次。

地震工程学-反应谱和地震时程波的相互转化matlab编程#(精选.)

地震工程学作业课程名称:地震工程学______ 指导老师:_______翟永梅_________ 姓名:史先飞________ 学号:1232627________一、地震波生成反应谱1 所取的地震波为Elcentro地震波加速度曲线，如图1所示。

图1 Elcentro地震波加速度曲线2 所调用的Matlab程序为：% ***********读入地震记录***********ElCentro;Accelerate= ElCentro(:,1)*9.8067;%单位统一为m和sN=length(Accelerate);%N 读入的记录的量time=0:0.005:(N-1)*0.005; %单位 s%初始化各储存向量Displace=zeros(1,N); %相对位移Velocity=zeros(1,N); %相对速度AbsAcce=zeros(1,N); %绝对加速度% ***********A,B矩阵***********Damp=0.02; %阻尼比0.02TA=0.0:0.05:6; %TA=0.000001:0.02:6; %结构周期Dt=0.005; %地震记录的步长%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;for T=0.0:0.05:6NatualFrequency=2*pi/T ; %结构自振频率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %计算公式化简e_t=exp(-Damp*NatualFrequency*Dt);s=sin(DampFrequency*Dt);c=cos(DampFrequency*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DampFrequency;A(2,1)=-NatualFrequency*e_t*s/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp^2-1)/(NatualFrequency^2*Dt);d_3t=Damp/(NatualFrequency^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c)-2*d_3 t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFrequency^2 )*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(NatualF requency^2*Dt);for i=1:(N-1) %根据地震记录,计算不同的反应Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i +1);Velocity(i+1)=A(2,1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i +1);AbsAcce(i+1)=-2*Damp*NatualFrequency*Velocity(i+1)-NatualFrequency^2*Displace(i+1);endMaxD(1,t)=max(abs(Displace));MaxV(1,t)=max(abs(Velocity));if T==0.0MaxA(1,t)=max(abs(Accelerate));elseMaxA(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,N);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End% ***********PLOT***********close allfigure %绘制地震记录图plot(time(:),Accelerate(:))title('PEER STRONG MOTION DATABASE RECORD')xlabel('time(s)')ylabel('acceleration(g)')gridfigure %绘制位移反应谱plot(TA,MaxD(1,:),'-.b',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=0.02')Gridfigure %绘制速度反应谱plot(TA,MaxV(1,:),'-.b',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k') title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=0.02')Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k') title('Absolute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s^2)')legend('ζ=0.02')Grid3 运行的结果得到的反应谱图2 位移反应谱图3 速度反应谱图4 加速度反应谱一、反应谱生成地震波1所取的反应谱为上海市设计反应谱图5 上海市设计反应谱2反应谱取值程序为：%%规范反应谱取值程序参照01年抗震规范function rs_z=r_s_1(pl,zn,ld,cd,fz) %%%pl 圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz %%%%烈度选择if ld==6arfmax=0.11;endif ld==7arfmax=0.23;endif ld==8arfmax=0.45;endif ld==9arfmax=0.90;end%%%%场地类别，设计地震分组选择if cd==1if fz==1Tg=0.25;endif fz==2Tg=0.30;endif fz==3Tg=0.35;endendif cd==2if fz==1Tg=0.35;if fz==2Tg=0.40;endif fz==3Tg=0.45;endendif cd==3if fz==1Tg=0.45;endif fz==2Tg=0.55;endif fz==3Tg=0.65;endendif cd==4if fz==1Tg=0.65;endif fz==2Tg=0.75;endif fz==3Tg=0.90;endend%%%%%%%%%ceita=zn; %%%%%阻尼比lmt1=0.02+(0.05-ceita)/8;if lmt1<0lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita); if lmt2<0.55lmt2=0.55;endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita); %%%%%分段位置 T1 T2 T3T1=0.1;T2=Tg;T_jg=2*pi./pl;%%%% 第一段 0～T1if T_jg<=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg;end%%%% 第二段 T1～T2if T1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%%%% 第三段 T2~T3if T2<T_jg&T_jg<=T3arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax;end%%%% 第四段 T3～6.0if T3<T_jg&T_jg<=6.0arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%% 第五段 6.0～if 6.0<T_jgarf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax;end%%%%%%反应谱值拟加速度值rs_z=arf_jg*9.8;end3生成人造地震波主程序：%%%主程序%%%%%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M，反应谱控制容差rc Tyz=[0.04:0.016:0.1,0.15:0.05:3.0,3.2:0.05:5.0];rc=0.06;nTyz=length(Tyz);ceita=0.035;%%%阻尼比：0.035for i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %%%%8度，2类场地，第1地震分组end%%%%%% 变换的频率差：2*pi*0.005(可以保证长周期项5s附近有5项三角级数)；%%%%频率变化范围 N1=30, 30*0.005*2*pi ;N2=3000, 5000*0.005*2*piplc=2*pi*0.005;pl=30*0.005*2*pi:0.005*2*pi:10000*0.005*2*pi;npl=length(pl);P=0.9; %%%保证率%%%%%%人造地震动持续时间40s，时间间隔：0.02sTd=40;dt=0.02;t=0:0.02:40;nt=length(t);%%%%%%% 衰减包络函数t1=8; %%%%上升段t2=8+24; %%%%%平稳段; 下降段则为40－32＝8sc=0.6; %%%%衰减段参数for i=1:ntif t(i)<=t1f(i)=(t(i)/t1)^2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2));endend%%%%%%% 反应谱转换功率谱for i=1:nplSw(i)=(2*ceita/(pi*pl(i)))*r_s_1(pl(i),ceita,8,2,1)^2/(-2*log(-1*pi*log(P)/(pl(i)*Td))); Aw(i)=sqrt(4*Sw(i)*plc);end%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfai(i)=rand(1)*2*pi;for j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=1;while max(rcrsp(:))>rc%%%%%循环体函数blxs=Syz./rspa;for xsxs=1:nplif 2*pi/pl(xsxs)<Tyz(1)blxs1(xsxs)=blxs(1);endfor sxsx=1:nTyz-1if (2*pi/pl(xsxs)>=Tyz(sxsx)) & (2*pi/pl(xsxs)<=Tyz(sxsx+1))blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxsx));endendif 2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfor j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=jsnum+1max(rcrsp(:))end%%%%%%% 最终的反应谱与规范谱%%%%%%%%%%%% response spectra of callidar%%%%%%% parameter%% Tjs=0.05:0.01:6;%% nTjs=length(Tjs);g=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i)/g;rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g;endsubplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4生成的人造地震波如图所示。

地震工程学-反应谱和地震时程波的相互转化matlab编程

地震工程学作业课程名称:地震工程学______指导老师:_______翟永梅_________姓名:史先飞________学号:1232627________一、地震波生成反应谱1 所取的地震波为Elcentro地震波加速度曲线，如图1所示。

图1 Elcentro地震波加速度曲线2 所调用的Matlab程序为：% ***********读入地震记录***********ElCentro;Accelerate= ElCentro(:,1)*9.8067;%单位统一为m和sN=length(Accelerate);%N 读入的记录的量time=0:0.005:(N-1)*0.005; %单位 s%初始化各储存向量Displace=zeros(1,N); %相对位移Velocity=zeros(1,N); %相对速度AbsAcce=zeros(1,N); %绝对加速度% ***********A,B矩阵***********Damp=0.02; %阻尼比0.02TA=0.0:0.05:6; %TA=0.000001:0.02:6; %结构周期Dt=0.005; %地震记录的步长%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;for T=0.0:0.05:6NatualFrequency=2*pi/T ; %结构自振频率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %计算公式化简e_t=exp(-Damp*NatualFrequency*Dt);s=sin(DampFrequency*Dt);c=cos(DampFrequency*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DampFrequency;A(2,1)=-NatualFrequency*e_t*s/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp^2-1)/(NatualFrequency^2*Dt);d_3t=Damp/(NatualFrequency^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c)-2*d_3 t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFrequency^2 )*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(NatualF requency^2*Dt);for i=1:(N-1) %根据地震记录,计算不同的反应Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i +1);Velocity(i+1)=A(2,1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i +1);AbsAcce(i+1)=-2*Damp*NatualFrequency*Velocity(i+1)-NatualFrequency^2*Displace(i+1);endMaxD(1,t)=max(abs(Displace));MaxV(1,t)=max(abs(Velocity));if T==0.0MaxA(1,t)=max(abs(Accelerate));elseMaxA(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,N);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End% ***********PLOT***********close allfigure %绘制地震记录图plot(time(:),Accelerate(:))title('PEER STRONG MOTION DATABASE RECORD')xlabel('time(s)')ylabel('acceleration(g)')gridfigure %绘制位移反应谱plot(TA,MaxD(1,:),'-.b',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=0.02')Gridfigure %绘制速度反应谱plot(TA,MaxV(1,:),'-.b',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k') title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=0.02')Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k') title('Absolute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s^2)')legend('ζ=0.02')Grid3 运行的结果得到的反应谱图2 位移反应谱图3 速度反应谱图4 加速度反应谱一、反应谱生成地震波1所取的反应谱为上海市设计反应谱图5 上海市设计反应谱2反应谱取值程序为：%%规范反应谱取值程序参照01年抗震规范function rs_z=r_s_1(pl,zn,ld,cd,fz) %%%pl 圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz %%%%烈度选择if ld==6arfmax=0.11;endif ld==7arfmax=0.23;endif ld==8arfmax=0.45;endif ld==9arfmax=0.90;end%%%%场地类别，设计地震分组选择if cd==1if fz==1Tg=0.25;endif fz==2Tg=0.30;endif fz==3Tg=0.35;endendif cd==2if fz==1Tg=0.35;if fz==2Tg=0.40;endif fz==3Tg=0.45;endendif cd==3if fz==1Tg=0.45;endif fz==2Tg=0.55;endif fz==3Tg=0.65;endendif cd==4if fz==1Tg=0.65;endif fz==2Tg=0.75;endif fz==3Tg=0.90;endend%%%%%%%%%ceita=zn; %%%%%阻尼比lmt1=0.02+(0.05-ceita)/8;if lmt1<0lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita); if lmt2<0.55lmt2=0.55;endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita); %%%%%分段位置 T1 T2 T3T1=0.1;T2=Tg;T_jg=2*pi./pl;%%%% 第一段 0～T1if T_jg<=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg;end%%%% 第二段 T1～T2if T1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%%%% 第三段 T2~T3if T2<T_jg&T_jg<=T3arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax;end%%%% 第四段 T3～6.0if T3<T_jg&T_jg<=6.0arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%% 第五段 6.0～if 6.0<T_jgarf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax;end%%%%%%反应谱值拟加速度值rs_z=arf_jg*9.8;end3生成人造地震波主程序：%%%主程序%%%%%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M，反应谱控制容差rc Tyz=[0.04:0.016:0.1,0.15:0.05:3.0,3.2:0.05:5.0];rc=0.06;nTyz=length(Tyz);ceita=0.035;%%%阻尼比：0.035for i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %%%%8度，2类场地，第1地震分组end%%%%%% 变换的频率差：2*pi*0.005(可以保证长周期项5s附近有5项三角级数)；%%%%频率变化范围 N1=30, 30*0.005*2*pi ;N2=3000, 5000*0.005*2*piplc=2*pi*0.005;pl=30*0.005*2*pi:0.005*2*pi:10000*0.005*2*pi;npl=length(pl);P=0.9; %%%保证率%%%%%%人造地震动持续时间40s，时间间隔：0.02sTd=40;dt=0.02;t=0:0.02:40;nt=length(t);%%%%%%% 衰减包络函数t1=8; %%%%上升段t2=8+24; %%%%%平稳段; 下降段则为40－32＝8sc=0.6; %%%%衰减段参数for i=1:ntif t(i)<=t1f(i)=(t(i)/t1)^2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2));endend%%%%%%% 反应谱转换功率谱for i=1:nplSw(i)=(2*ceita/(pi*pl(i)))*r_s_1(pl(i),ceita,8,2,1)^2/(-2*log(-1*pi*log(P)/(pl(i)*Td))); Aw(i)=sqrt(4*Sw(i)*plc);end%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfai(i)=rand(1)*2*pi;for j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=1;while max(rcrsp(:))>rc%%%%%循环体函数blxs=Syz./rspa;for xsxs=1:nplif 2*pi/pl(xsxs)<Tyz(1)blxs1(xsxs)=blxs(1);endfor sxsx=1:nTyz-1if (2*pi/pl(xsxs)>=Tyz(sxsx)) & (2*pi/pl(xsxs)<=Tyz(sxsx+1))blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxsx));endendif 2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfor j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=jsnum+1max(rcrsp(:))end%%%%%%% 最终的反应谱与规范谱%%%%%%%%%%%% response spectra of callidar%%%%%%% parameter%% Tjs=0.05:0.01:6;%% nTjs=length(Tjs);g=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i)/g;rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g;endsubplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4生成的人造地震波如图所示。

地震工程学大作业代码3求人工波

w=linspace(0,100,10000); %生成周期矩阵kesai=[0.72 0.8]; %输入阻尼矩阵w0=[20.94 13.96]; %输入周期矩阵w1=0.15*w0;kesai1=kesai;t=linspace(0,40,2001); %生成时间矩阵，便于后期画图t1=[0.8 1.2]; %输入强度函数相关参数t2=[7 9];c=[0.35 0.25];file_id='D:\人造波\';mkdir(file_id); %创建文件夹保存生成的所有图和地震动数据for i=1:2 %利用循环以此计算两类地震动数据for j=1:10000 %利用循环计算功率谱a=(1+4*kesai(1)^2*(w(j)/w0(i))^2)*(w(j)/w1(i))^4;b=(1-(w(j)/w0(i))^2)^2+4*kesai(i)^2*(w(j)/w0(i))^2;d=(1-(w(j)/w1(i))^2)^2+4*kesai1(i)^2*(w(j)/w1(i))^2;S(j)=a*0.0252/(b*d);A(j)=sqrt(2*S(j)*(w(2)-w(1))/pi); %计算对应振幅end %功率谱计算结束fai=repmat(2*pi*rand(10000,1)-pi,1,2001);%生成相位矩阵for n=1:2001 %计算强度函数值gt(n)=g(t(n),t1(i),t2(i),c(i));endacceleration=A*cos(w'*t+fai).*gt; %计算地震动数据N= find(abs(acceleration)>=0.02*max(abs(acceleration)));%寻找截断点figure(2*i-1); %利用figure函数控制图形窗口plot(w,S); %得到功率谱图name=strcat(num2str(i+1),'类场地人工波功率谱');title(name);set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看F=getframe(gcf);imwrite(F.cdata,[file_id,strcat(name,'.png')]); %存储功率谱图形figure(2*i);plot(t(1:N(end)),acceleration(1:N(end))); %画地震动时程图name=strcat(num2str(i+1),'类场地人工波加速度时程');title(name);set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看F=getframe(gcf);imwrite(F.cdata,[file_id,strcat(name,'.png')]); %存储地震动时程图xlswrite(strcat(file_id,'第',num2str(i+1),'类场地人工波.xls'),[t(1:N(end));acceleration(1:N(end))]');%存储地震动数据end%%function gt=g(t,t1,t2,c)%强度包络函数计算if t<t1gt=(t/t1)^2;elseif t>t2gt=exp(-c*(t-t2));elsegt=1;endend。

地震射线弯曲法程序matlab

在地球科学领域中，地震射线弯曲法程序matlab是一个非常重要的工具。

地震是地球内部的一种自然现象，而地震射线弯曲法则是利用地震波在地球内部的传播路径来研究地球结构的一种方法。

在这篇文章中，我将深入探讨地震射线弯曲法程序matlab的原理、应用和意义，帮助你更全面地理解这一重要领域。

1. 地震射线弯曲法程序matlab是什么？地震射线弯曲法程序matlab是一种针对地震波传播路径和速度结构进行研究的数值模拟工具。

它利用了地震波在地球内部传播的特点，通过对地震波传播路径的模拟和反演，可以获得地球内部结构的信息。

在地震勘探、地质勘探和地球物理学研究中，地震射线弯曲法程序matlab扮演着至关重要的角色。

2. 地震射线弯曲法程序matlab的原理和方法地震射线弯曲法程序matlab的原理基于地震波在地球内部传播的轨迹。

当地震波穿过地球内部的不同介质时，由于介质的不均匀性，地震波路径会产生弯曲和偏折。

地震射线弯曲法程序matlab通过模拟地震波传播路径的弯曲情况，来推断地球内部的速度结构和地质构造。

它基于数学和物理原理，利用计算机编程进行地震波路径的数值模拟和反演，从而揭示地球内部结构的特征。

3. 地震射线弯曲法程序matlab的应用领域地震射线弯曲法程序matlab在地球科学领域有着广泛的应用。

它可以用于地震勘探，帮助地震学家和地球物理学家了解地球内部的速度结构和地质构造，从而寻找地下资源和研究地壳运动。

地震射线弯曲法程序matlab还可以用于地质勘探，帮助地质学家识别地下的地质体和构造，为工程建设和地质灾害防治提供重要信息。

4. 地震射线弯曲法程序matlab的意义和价值地震射线弯曲法程序matlab的意义和价值在于它可以帮助我们更好地理解地球内部的结构和演化过程。

通过对地震波传播路径的模拟和反演，我们可以揭示地球内部的奥秘，为地质学、地球物理学和地球科学的发展提供重要支持。

地震射线弯曲法程序matlab还可以为地球资源勘探和自然灾害预警提供科学依据，对于人类的生存和发展具有重要意义。

单自由度地震能量计算的程序matlab

funct ion responsenew%基于戏线性滞冋模型的单自由度体系的地震能虽分析程序%made by Yanan Li Faculty of Infrastructure Enginoering弔Dalian University of Technology%质量57041kg,阻尼36612 N s/叫初始刚度2350000N/m,刚度折减系数0. 2,屈服位移0・01m, 采用ELCENTRO波%参数替换直接在下面修改，然后运行clcformat long;m=57041;% 质量ug= importdataC* ELCENTRO. txt') ;%地震波txt 文件ug=ug/100;P=-m*ug;num=size(P, 1);c 二36612; % 阻尼k 1=2350000; %初始刚度k2二kl*0・2;a=zeros(l t num) ;v=zeros(l, num) ;x=zeros(l, num);%加速度速度位移Ei=zeros(l, num);Ek二zeros(1, num);Ed二zeros(1, num); Eh二zeros (1, num)输入能动能阻尼耗能滞回耗能EI=zeros(l t num); EK二zeros (1, num) ;ED=zeros(l, num) ;EH=zeros(l, num) ;%累积的各种能量t imc=zeros(l, num);a(l)=P (l)/m;hfl=zeros (1, num);h=0. 02;%地震波采样间隔Xy=0.01;%屈服位移pxmax=0;nxmax-0;pd=0;%双线型滞冋模型折线的标识，0农示弹性，1表示正向弹犁性，2表示反向弹性，-1 表示反向弹塑性，-2表示正向弹性for ii=l:numif pd==0 %弹性阶段k二kl;if x(ii)>Xypd=l;b=[(a(ii)-a(ii-l))/6/h a(ii~l)/2 v(ii-l) x(iiT)-Xy];% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==ldth=d(j);endendhp二h-dth;vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth*2/h/2);ap=a(i i-l) + (a (i i)-a(ii-l)) *dth/h;pp=m*ap+c*vp+kl*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(i i +1) -pp+m* (6*vp/hp+3*ap) +c* (3*vp^-hp*ap/2); dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=Xy+dtx;v(ii)=vp^dtv;dtf=k2*dtx;hfl(ii)=kl*Xy+dtf;a(ii) = (P(ii)-hfl (i i+1)-c*v(i i))/m;elseif x(ii)<-Xypd=-l;b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(iiT) x(iiT)+Xy];% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==ldth=d(j);endendhp=h-dth;vp=v(i i 1)・a(i i 1)*dth+ ((a(ii)-a(ii-l))*dth"2/h/2); ap a(ii-l) + (a(ii)-a(ii-l))*dth/h;pp=m*ap+c*vp_k1*Xy; kd=k2+3*c/hp+6*m/hp/hp;dtpd=P (i i + 1)-pp+m*(6*vp/hp+3*ap)(3*vp^hp*ap/2);dtx=dtpd/kd;dtv-3*dtx/hp-3*vp-hp*ap/2;x(ii)=-Xy+dtx： v(i i)=vp^dtv;dtf=k2*dtx; hfl(ii)=-kl*Xy+dtf;a(i i) = (P(i i)-hf 1 (i i)-c*v(ii))/m;endendif pd==l %正向弹塑性k=k2;if v(ii)<0pd=2;b=[(a(ii)-a(ii-l))/2/h a(ii-l) v(ii-l)];% 拐点处理d二rools(b);for j=l:length(d)if isreal (d(j))==ldth=d (j);endendhp=h-dth;xp=x(iiT)+v(i i-1) *dth+a (i i-1) *dth*2/2+ (a (i i) -a (i i-l))*dth*3/h/6; ap=a (i i-1)+ (a (i i) -a (i i-1)) *dth/h;pp-m*ap+c*O+k1*Xy+k2* (xp-Xy);kd=kl+3*c/hp+6*m/hp/hp:dtpd=P(ii+l)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(i i)=xp+dtx;v(ii)=O+dtv;dtf=kl*dtx;hfl(ii)=kl*Xy+k2*(xp-Xy)+dtf;a(ii) = (P(ii)-hfl (ii)-c*v(ii))/m;pxmax=xp;endif pd==2 %反向弹性k=kl;if x (i i)>pxmaxpd=l; ■b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(ii-l) x(ii-l)-pxmax]:% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==I dth=d(j);endendhp=h-dth; vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth"2/h/2); ap=a(ii-l)+(a(ii)-a(ii-l))*dth/h;pp-m*ap+c*vp+k1*Xy+k2*(pxmax-Xy); kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd;d t v=3*d1x/hp-3*vp-hp*ap/2;x(i i)二pxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl (i i)=kl*Xy+k2*(pxmax-Xy) +dtf;a(ii) = (P(ii)-hfl (ii)-c*v(ii))/m;el seif x (i i)<(pxmax~2*Xy)pd=-l;―.・pu og 、(( =)A *3—(二二Jq —(二)d )H(二)P,±P +AX盖〒Ax*上—殳艾舌 X F (二二jqEp*znp 2p+dA:=(=)A-xlp+AX*zlxpulxdu(二)x,z 、dB*dqld 家gldv、xlp^HA4p&<p d "P H X 4p-(z 、dp*dfdA*g)芯+(du*g+dq、dA*9)*u+dd —(7二)dupdpp -dq'dq'lu 关 9+dqb苗+0上HP上-(xx*^IAX*>l —z>l*xeulxd) +d A *o +d e *u "dd2>-p*(〒二)E —(=)B)+U丄一)fde二z 'q 'Fe p 決((TT OF(U)P ))+±P*II )P +UI )AH D A召p —qudqp u a)pUO」(「)p*pTH((「)P 二p冒岛(P)£8U 上二丄joj-(q)s_oof p三P+O A=)Awp+dxH(二)x-z 、dydq —0^—dq、w'pdup,(z 、d¥dq+o為)*¥(dp*e+dq/0*9)*lu+dd —(I +二)dupdlp‘ d<dq 、ul*9+d<。