求助：利用三角级数和模拟人工地震波（MATLAB）

double_li

打算考虑二维相干性（Hao相干模型），采用杜修力-陈厚群功率谱模型，用三角级数和的方法生成空间多点的地震时程曲线，用MATLAB编，不熟练，程序中有错误，哪位高手帮忙看看，另外有类似程序能否学习一下。学位论文急用，谢谢！
M文件和相关理论依据见附件。

jimin

用MATLAB编，不熟练，程序中有错误
呵呵，好多循环呀
是什么错误，是不能运行出来呢，还是算法有问题
前者的话，你看提示修改吧

[ 本帖最后由 jimin 于 2006-9-19 16:03 编辑 ]

AaronSpark

1. ll=chol(gama);

复制代码

这一句有问题
用chol这个命令，其中gama必须是方阵
显然在你的循环过程中，大部分都不满足这个要求

double_li

将其放在第二层循环的外部就可以，这句话已经调通，程序好像可以通了，还在看结果如何:@L:@L

[ 本帖最后由 ChaChing 于 2009-12-30 18:34 编辑 ]

wuhongju

为什么看不了啊

friendchj

可能时间太久了，附件不好用了。

xuelei

麻烦楼主吧附件重新传一下！

liujiashun

烦劳楼主将附件重新传一下 下不了但是很有用啊！

pansidong311

请问 楼主 你的问题解决了没有？

pansidong311

回复 1 # double_li 的帖子

你好，这个程序不能下载啊，我正在做论文，也遇到这个问题，不知大侠能否给我发份程序，非常感激啊，pansidong311@163.com

double_li

这个帖子好久以前的，没想到还有人回
那个程序最后改好了，多交流。

%%%%%%%%%%%%%三角级数法模拟人造地震波程序%%%%%%%%%%%%%%%%

% 本程序利用杜修力－陈厚群功率谱，考虑空间二维相干的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;
%自功率谱中参数D

omiga0=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=[0 17.8 -17.8 17.8 -17.8 0

-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

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

end

end

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

end

end

%非平稳时程计算

for ii=1:N

for j=1:nn

xxg2(ii,j)=xxg(ii,j)*fai(ii,j);

x(j)=j*dt;

end

end

%采用比例法调整

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

end

yy=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,:))

seismic123

好东西啊，正需要呢，谢谢啊

jack_chen

回复 12 # seismic123 的帖子

能把这个模拟地震波的程序给我传一个吗？多谢了 我邮箱：chenxiangxh@126.com

jack_chen

回复 11 # double_li 的帖子

这个附件下不了啊，大侠能给我传个这个模拟地震波的程序吗？多谢了 我邮箱：chenxiangxh@126.com

无一涯

收下了，十分感谢啊