帮我看看这个算法的matlab描述
这里有公式的内容[ 本帖最后由 wuqiong 于 2008-12-29 17:49 编辑 ]
帮我看看这个算法的matlab描述
这里是一个根据三份量地震波记录求解方位角的算法：震源方位角是运用垂直及水平分量记录的互相关函数(cross-correlations)来求取的（Nakamura ，1988）。垂直及水平分量的互相关函数，可以由以下公式所求得：
其中
为第i点垂直与东西分量的互相关函数
为第i-1点垂直与东西分量的互相关函数
为第i点垂直与南北分量的互相关函数
为第i-1点垂直与南北分量的互相关函数
为第i点垂直分量记录
为第i点东西分量记录
为第i点南北分量记录
为平滑常数(0.0-1.0)
由公式(1)及(2)求得互相关函数之后，再带入下式可求得地震波入射之方位角 。
---(3)
方位角 之所在象限和 与 之正负符号有关，其对应关系如下：
象限
I
负
负
II
正
负
III
正
正
IV
负
正
我的matlab函数为：
function az = azimath_naka(record,alfa,type)
%this function is used to compute the azimth of the seismic record by using
%the method of NAKAMURA. and the alfa is a moving average parameter range
%from 0.0 to 1.0. the parameter type is the type of the seismeter:
% 1 ------------------- vocerty
% 2 ------------------- displacement
% 3 ------------------- aceleration
if nargin ~= 3
error('the parameter type and alfa must be confirmed at first');
end
if type == 3
record(:,1) = cumtrapz(record(:,1));
record(:,2) = cumtrapz(record(:,2));
record(:,3) = cumtrapz(record(:,3));
end
ud = detrend(record(:,1));
ew = detrend(record(:,2));
ns = detrend(record(:,3));
% figure
% subplot(311),plot(ud);ylabel('ud')
% subplot(312),plot(ns);ylabel('ns')
% subplot(313),plot(ew);ylabel('ew')
%compute the crosscorr of the ud-ew and ud-ns
in = length(ud);
maxlag = in;samp_seg = 100;overlap = 0;flag = 'biased';
% Rue= cum2x(ud,ew,maxlag, samp_seg, overlap, flag);
Rue = xcorr(ud,ew,'unbiased');
Rue = Rue(maxlag+1:end);
% Run = cum2x(ud,ns,maxlag, samp_seg, overlap, flag);
Run = xcorr(ud,ns,'unbiased');
Run = Run(maxlag+1:end);
figure
subplot(211),plot(Rue)
subplot(212),plot(Run)
for i = 2:in
Rue(i) = Rue(i-1)*alfa + ud(i)*ew(i);
Run(i) = Run(i-1)*alfa + ud(i)*ns(i);
Tue = Rue(i);
Tun = Run(i);
temp = atan(Tue./Tun)*180/pi;
if (Tue<0)&(Tun<0)
az(i) = temp;
elseif (Tue>0)&(Tun<0)
az(i) = 180 - abs(temp);
elseif (Tue>0)&(Tun>0)
az(i) = 180 + abs(temp);
elseif (Tue<0)&(Tun>0)
az(i) = 360 - abs(temp);
end
end
不知道怎么了，做得结果总是不对，不知道我的程序是什么地方错了，请来看看。郁闷我好久了，不好意思啊
[ 本帖最后由 wuqiong 于 2008-12-29 17:47 编辑 ] 附件中可以看到公式 是出现运行问题 还是结果不正确 分量记录什么的是否已知 不是运行问题,是结果不对.
页:
[1]