一般均值服从正态分布,幅值服从威布尔分布,拟合求得正态分布均值muHat,标准差sigmaHat;三参数威布尔分布尺度参数A,形状参数B,位置参数C。然后根据极值概率,求得均值、幅值的最大值Xmax、Ymax,最后根据联合概率密度函数计算形成8×8二维载荷谱。
x为载荷均值,y为载荷幅值,,
,
,
,
。
N为8×8二维载荷谱:
N=nn(muHat,sigmaHat,A,B,C,Xmin,Xmax,Ymax);
主程序:
function NN = nn(muHat,sigmaHat,A,B,C,Xmin,Xmax,Ymax)
% 8×8二维载荷谱计算
muHat=muHat; % 均值
sigmaHat=sigmaHat; % 标准差
A = A; % 尺度参数
B = B; % 形状参数
C = C; % 阈值参数
Xmin=Xmin;
Xmax=Xmax;
Ymax=Ymax;
% 均值均分为8级
for i = 2:8
M(:,1) = Xmin;
M(:,i) = M(:,i-1) + (Xmax-Xmin)/7;
end
% 幅值分级(不等间距)分为8级
pc = [0.125 0.275 0.425 0.575 0.725 0.85 0.95 1];
R = Ymax.*pc;
n=1e6;% 极值概率
% 联合概率密度函数
fun = @(x,y) (1/(sqrt(2*pi)*sigmaHat)*exp(-(x-muHat).^2/(2*sigmaHat^2))).*(B/A*((y-C)/A).^(B-1).*exp(-((y-C)/A).^B));
%fun = @(x,y) 1/(sqrt(2*pi)*sigmaHat).*exp(-(x-muHat).^2/(2*(sigmaHat^2))) .* (B/A*(y/A).^(B-1).*exp(-(y/A).^B));
% 均值x的下、上限u_w,u_w+1; 幅值y的下、上限v_j,v_j+1;
% u_1=Xmin; u_9=Xmax; v_1=Ymin=0.125*Ymax; v_9=Ymax;
% 主程序:
for i=1:8
if i==1
ymin = 0.125*Ymax;
ymax = (R(1)+R(2))/2;
for j=1:8
if j==1
xmin = Xmin;
xmax = (M(1)+M(2))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
elseif j==8
xmin = (M(7)+M(8))/2;
xmax = Xmax;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
else
xmin = (M(j-1)+M(j))/2;
xmax = (M(j)+M(j+1))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
end
end
elseif i==8
ymin = (R(7)+R(8))/2;
ymax = Ymax;
for j=1:8
if j==1
xmin = Xmin;
xmax = (M(1)+M(2))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
elseif j==8
xmin = (M(7)+M(8))/2;
xmax = Xmax;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
else
xmin = (M(j-1)+M(j))/2;
xmax = (M(j)+M(j+1))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
end
end
else
for j=1:8
if j==1
xmin = Xmin;
xmax = (M(1)+M(2))/2;
ymin = (R(i-1)+R(i))/2;
ymax = (R(i)+R(i+1))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
elseif j==8
xmin = (M(7)+M(8))/2;
xmax = Xmax;
ymin = (R(i-1)+R(i))/2;
ymax = (R(i)+R(i+1))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
else
xmin = (M(j-1)+M(j))/2;
xmax = (M(j)+M(j+1))/2;
ymin = (R(i-1)+R(i))/2;
ymax = (R(i)+R(i+1))/2;
NN(i,j)=n*integral2(fun,xmin,xmax,ymin,ymax);
end
end
end
end
% 四舍五入
NN = round(NN);
end