2024年2月5日发(作者：)

MATLAB雨流计数法程序

一、仿真随机过程

sw=2;

A=sqrt(2*sw*detaw)

（1）中心频率为10pi，detaw=0.01，带宽为2pi的低频窄带随机过程

for n=9:0.01:11;

xn=A*sin(n*pi*t+rand(1,1)*2*pi);

x=x+xn;

end

（2）中心频率为100pi，detaw=0.01，带宽为2pi的高频窄带随机过程

for n=99:0.01:101;

yn=A*sin(n*pi*t+rand(1,1)*2*pi);

y=y+yn;

end

（3）低频与高频组合成的宽带随机过程

z=x+y

低频窄带随机过程1050-5-16高频窄带随机过程1050-5-16组合成的宽带随机过程1050-5-16

问题一：关于带宽。带宽为2pi似乎太宽了，若带宽改为0.2pi，detaw依然不变，得到图如下。

低频窄带随机过程10.50-0.5-1高频窄带随机过程0.40.20-0.2-0.4组合成的宽带随机过程10.50-0.5-1

二、用雨流计数法截取循环

（1）得到随机过程x的极点值（储存在矩阵S中）。

n=length(x);

s=x(1);

for i=2:n-1;

if (x(i)>x(i-1))&&(x(i)>x(i+1))||((x(i)

s=[s,x(i)];

end

end

S=[s,x(n)];

（2）对极值点间变程进行比较分析，截取循环

while i+3

i=i+1;

detaS1=abs(S(i+1)-S(i)); %获得四个极值点间的三段变程

detaS2=abs(S(i+2)-S(i+1));

detaS3=abs(S(i+3)-S(i+2));

if (detaS2<=detaS1)&&(detaS2<=detaS3) %将变程满足截取条件的循环截取出

Ba=[Ba,(S(i+2)-S(i+1))/2]; %截取出循环的幅值

Bm=[Bm,(S(i+2)+S(i+1))/2]; %截取出循环的均值

Q=[Q;S(i+1),S(i+2)]; %构成截取出循环的极值点存于Q中，即原随机过程抛弃的点

S(i+1)=[];

S(i+2)=[];

i=i-1;

end

end

Sr=S; %剩余的极值点存于Sr中

（3）分别对低频窄带、高频窄带、宽带进行雨流处理

低频窄带

处理前载荷1050-5-1处理后的载荷1050-5-100

高频窄带

处理前载荷6420-2-4-6010001200处理后的载荷6420-2-4-6

宽带

处理前载荷1050-5-1处理后的载荷1050-5-1

问题二：由于观察到用雨流法处理这三种随机过程的效果差别比较大，就编程统计计算了雨流处理这三种随机过程抛弃的极值点数（即截取的循环）与原载荷极值点数之比rate。因为我理解为rate越大，雨流处理效果越好（不知道对不对？）得到的统计结果如下：

首先，每种随机过程都对一个信号雨流处理的结果比较

雨流 低频窄带随机过程 高频窄带随机过程 组合宽带随机过程

采样点数N1 1001 1001 1001

极值点数N2

剩余点数N3

截取点数N4

943

783

160

82

38

44

0.5366

941

277

624

0.6992 截取率rate=N4/N2 0.1697

其次，每种随机过程都对很多个信号雨流处理的统计结果（对100组信号进行统计）:

rate

低频窄带随机过程 0.0814

频窄带随机过程 0.3006

组合宽带随机过程 0.3178

这样看的话这种雨流计数法对这三种随机过程的适用性：宽带>高频窄带>低频窄带（自己这么想的，不知是否行的通）

三、对剩余载荷的处理

re=[ZSr,ZSr];

[reside,Ba,Bm,Q]=Rainflow(re);

residue2150200residue+residue25100-100residue+cycles515

问题三：我的理解是对于剩余载荷residue，在复制一个一样的residue，让它们首尾相接作为一个新的过程再进行雨流处理，得到{residue2+cycles}。我对文章的理解是residue2要和原本的residue一致，但这又如何保证呢？总之，对这里理解的不是很清晰。

四、载荷的重构（将循环插入剩余载荷中）

（1）筛选出要插入循环Q中高于门槛值threshold的循环存于S中

[c,r]=size(Q);

S=[];

Abs=[];

for i=1:c

if abs(Q(i,2)-Q(i,1))>threshold

S=[S;Q(i,1),Q(i,2)];

Abs=[Abs,abs(Q(i,2)-Q(i,1))];

end

end

（2）对S中循环按幅值大小降序排列

[Abs,index]=sort(Abs);

index=fliplr(index);

k1=S(:,1);S(:,1)=k1(index);

k2=S(:,2);S(:,2)=k2(index);

（3）按幅值由大到小插入载荷序列中（这部分编出来程序运行没有问题，但结果显示似乎并没有插进去，还在寻找问题中…）

[c,r]=size(S);

no1=[];no2=[];

while m<=c

if S(m,1)

n=length(reside);

rde=[];indexrde=[];

for i=1:n-1

if (reside(i)>=S(m,2))&&(reside(i+1)<=S(m,1)) %筛选出第m个循环能够插入的序列

rde=[rde;reside(i),reside(i+1)];

indexrde=[indexrde;i,i+1];

end

end

[k,g]=size(rde);

if k==0

m=m+1;

else

chu=randi(k,1,1); %在所有能够插入的序列中随机选取一个插入

p=indexrde(chu,1);

reside=[reside(1:p),S(m,1),S(m,2),reside(p+1:length(reside))];

m=m+1;

end

else

n=length(reside);

rin=[];indexrin=[];

for i=1:n-1

if (reside(i)<=S(m,2))&&(reside(i+1)>=S(m,1))

rin=[rin;reside(i),reside(i+1)];

indexrin=[indexrin;i,i+1];

end

end

[k,g]=size(rin);

if k==0

m=m+1;

else

chu=randi(k,1,1);

p=indexrin(chu,1);

reside=[reside(1:p),S(m,1),S(m,2),reside(p+1:length(reside))];

m=m+1;

end

end

end