求助：利用三角级数和模拟人工地震波（MATLAB）
打算考虑二维相干性（Hao相干模型），采用杜修力-陈厚群功率谱模型，用三角级数和的方法生成空间多点的地震时程曲线，用MATLAB编，不熟练，程序中有错误，哪位高手帮忙看看，另外有类似程序能否学习一下。学位论文急用，谢谢！M文件和相关理论依据见附件。 用MATLAB编，不熟练，程序中有错误
呵呵，好多循环呀
是什么错误，是不能运行出来呢，还是算法有问题
前者的话，你看提示修改吧
[ 本帖最后由 jimin 于 2006-9-19 16:03 编辑 ] ll=chol(gama);
这一句有问题
用chol这个命令，其中gama必须是方阵
显然在你的循环过程中，大部分都不满足这个要求 将其放在第二层循环的外部就可以，这句话已经调通，程序好像可以通了，还在看结果如何:@L:@L
[ 本帖最后由 ChaChing 于 2009-12-30 18:34 编辑 ] 为什么看不了啊
回复 6楼 wuhongju 的帖子
可能时间太久了，附件不好用了。 麻烦楼主吧附件重新传一下！附件重新传
烦劳楼主将附件重新传一下 下不了但是很有用啊！ 请问 楼主 你的问题解决了没有？ 回复 1 # double_li 的帖子你好，这个程序不能下载啊，我正在做论文，也遇到这个问题，不知大侠能否给我发份程序，非常感激啊，pansidong311@163.com 这个帖子好久以前的，没想到还有人回
那个程序最后改好了，多交流。
%%%%%%%%%%%%%三角级数法模拟人造地震波程序%%%%%%%%%%%%%%%% % 本程序利用杜修力－陈厚群功率谱，考虑空间二维相干的Hao相干函数模型的沿桥纵向地震波 % 变量说明： % s0: 谱强度因子 % dt: 地震波采样时间间隔 % kao: 地震波总时间长度 % Ts: 平稳地震动的持续时间 % t1,t2,C: 控制主震段首末时间和衰减快慢的参数 % keceg,omigag: 地表覆盖层的阻尼比和卓越频率% d,omiga0:
与地震震级有关的参数% omiga1: 截断频率% gama: 相关系数% beita1,beita2,alf1,alf2,a,b,c,d,e,g: Hao相干函数中参数% ddl,ddt: 距震中点水平距和垂直震中距间距离% vapp: 视波速% ak、fan合成系数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear %******************************* nn=1000;
%频率采样点数N=6;
%支撑点数Ts=9.0;
%平稳地震动的持续时间amax=1.084;
%最大加速度s0=0.002797665;
%谱强度omigag=13.16;
%场地卓越频率keceg=0.97;
%场地阻尼比dd=0.01137;
%自功率谱中参数Domiga0=1.83;
%自功率谱中参数omiga1=25*pi;
%截断频率beita1=0.000225;
%相干函数中系数beita2=0.00051;
%相干函数中系数a=0.01066;
b=0.0000265;c=-0.0001;d=0.006655;e=0.000059;g=-0.00112;
%相干函数中系数ddl=[0 77.8 77.8 265.8 265.8 343.6
-77.8 0 0 188 188 265.8
-77.8 0 0 188 188 265.8
-265.8 -188 -188 0 0 77.8
-265.8 -188 -188 0 0 77.8
-343.6 -265.8 -265.8 -77.8 -77.8 0];
%顺桥向位置ddt=
-17.8 0 -35.6 0 -35.6 -17.8
17.8 35.6 0 35.6 0 17.8
-17.8 0 -35.6 0 -35.6 -17.8
17.8 35.6 0 35.6 0 17.8
0 17.8 -17.8 17.8 -17.8 0];
%横桥向位置%******************************* gama=
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1];
ss=zeros(nn,1);
xxg=zeros(N,nn);
xg=zeros(N,N,nn);
ggg=zeros(N,N,nn);
dt=0.02;
cc=0.5;
%控制主震段首末时间和衰减快慢的参数Domiga=25*pi/nn;
%频率变化值for ii=1:nn %计算功率谱
omiga(ii)=Domiga*(ii);
Gomiga1(ii)=4*(keceg*omiga(ii)/omigag)^2;
Gomiga2(ii)=(1-(omiga(ii)/omigag)^2)^2;
Gomiga3(ii)=(1+Gomiga1(ii))/(Gomiga2(ii)+Gomiga1(ii));
Gomiga4(ii)=s0/(1+(dd*omiga(ii))^2);
Gomiga5(ii)=omiga(ii)^4;
Gomiga6(ii)=(omiga0^2+omiga(ii)^2)^2;
Gomiga(ii)=Gomiga3(ii)*Gomiga4(ii)*Gomiga5(ii)/Gomiga6(ii);
qrand(ii)=rand(1)*2*pi;
ss(ii)=Gomiga(ii);
%随频率变化的谱强度计算
vapp(ii)=3344+1095*log(omiga(ii)/2/pi);
end%计算相干函数
for ii=1:nn
if omiga(ii)<=20*pi
alf1(ii)=2*pi*a/omiga(ii)+b*omiga(ii)/2/pi+c;
alf2(ii)=2*pi*d/omiga(ii)+e*omiga(ii)/2/pi+g;
else
alf1(ii)=2*pi*a/(20*pi)+b*(20*pi)/2/pi+c;
alf2(ii)=2*pi*d/(20*pi)+e*(20*pi)/2/pi+g;
end
for j=1:N
for k=1:N
gama1(j,k,ii)=exp(-(beita1*abs(ddl(j,k))+beita2*abs(ddt(j,k))));
gama2(j,k,ii)=exp(-(alf1(ii)*sqrt(abs(ddl(j,k)))+alf2(ii)*sqrt(abs(ddt(j,k))))*(omiga(ii)/2/pi)^2);
gama(j,k)=gama1(j,k,ii)*gama2(j,k,ii);
end
end
ll=chol(gama);
%乔列斯基分解
for j=1:N
for k=j:N
ggg(j,k,ii)=ll(j,k);
end
end
for j=1:N
for k=j+1:N
ggg(k,j,ii)=ll(j,k);
end
end%计算合成系数
for j=1:N
for k=1:N
%幅值
ak(j,k,ii)=ggg(j,k,ii)*sqrt(4*ss(ii)*Domiga);
end
endend
for j=1:N
for k=1:N
for ii=1:nn
%相位差
ff(j,k,ii)=omiga(ii)*ddl(j,k)/vapp(ii);
%行波效应
hh(j,k,ii)=sqrt(ss(ii))*ggg(j,k,ii)*(exp(-1*ff(j,k,ii)*i));
%相位角
fan(j,k,ii)=atan(imag(hh(j,k,ii))/real(hh(j,k,ii)));
end
end
end%随机相位角
for ii=1:N
for j=1:nn
qrand(ii,j)=rand(1)*2*pi;
end
end%平稳样本合成
for ii=1:N
for kk=1:nn
for k=1:N
for j=1:nn
xg(ii,k,kk)=xg(ii,k,kk)+ak(ii,k,j)*cos(omiga(j)*kk*dt+fan(ii,k,j)+qrand(k,j));
end
xxg(ii,kk)=xxg(ii,kk)+xg(ii,k,kk);
end
end
end%调制函数计算
for ii=2:N
for j=1:nn
t1(1,j)=1.2;
t2(1,j)=10.2;
t1(ii,j)=t1(1,j)+0.00034*ddl(1,ii);
t2(ii,j)=t1(1,j)+9+0.002*ddl(1,ii)+t1(ii,j);
end
end
for ii=1:N
for j=1:nn
if j*dt<t1(ii,j)
fai(ii,j)=(j*dt/t1(ii,j))^2;
elseif (j*dt>=t1(ii,j)&j*dt<t2(ii,j))
fai(ii,j)=1;
else j*dt>=t2(ii,j)
fai(ii,j)=exp(-cc*(j*dt-t2(ii,j)));
end
endend%非平稳时程计算for ii=1:N
for j=1:nn
xxg2(ii,j)=xxg(ii,j)*fai(ii,j);
x(j)=j*dt;
endend%采用比例法调整
for ii=1:N
xxg1(ii)=max(abs(xxg2(ii,:)));
amax1(ii)=amax/xxg1(ii);
end
for ii=1:N
for kk=1:nn
y(ii,kk)=xxg2(ii,kk)*amax1(ii);
end
endyy=y';subplot(2,3,1);plot(x,y(1,:))subplot(2,3,2);plot(x,y(2,:))subplot(2,3,3);plot(x,y(3,:))subplot(2,3,4);plot(x,y(4,:))subplot(2,3,5);plot(x,y(5,:))subplot(2,3,6);plot(x,y(6,:))
好东西啊，正需要呢，谢谢啊 本帖最后由 jack_chen 于 2012-3-15 15:53 编辑
回复 12 # seismic123 的帖子
能把这个模拟地震波的程序给我传一个吗？多谢了 我邮箱：chenxiangxh@126.com 本帖最后由 jack_chen 于 2012-3-15 15:53 编辑
回复 11 # double_li 的帖子
这个附件下不了啊，大侠能给我传个这个模拟地震波的程序吗？多谢了 我邮箱：chenxiangxh@126.com 收下了，十分感谢啊
页:
[1]
2