matlab海杂波统计建模及CFAR阈值分割

海洋状况不同，选取的海杂波分布模型不同。低海况均匀杂波分布一般服从较简单的统计模型，如高斯分布、瑞利分布、

对数正态分布、韦布尔分布等，而高海况不均匀杂波由于其较重的拖尾特性一般由K分布来描述。

SAR图像统计建模按其建模过程可大致分为两类：参量模型和非参量模型。参量模型的建模过程是首先给出图像数据大致

服从的已知模型分布，其模型参数未知，然后通过图像区域数据对未知的模型参数进行估计，最后按照一定的拟合优度

评价准则选取最优的概率分布，作为图像区域数据的模型分布。相对应的，非参量模型的建模过程是事先不假设图像区域

数据服从某种概率模型分布，而是通过图像区域的数据驱动方式获取最优的概率分布模型。由于非参量模型的建模过程

涉及到复杂的逼近操作，计算耗时大，并且需要大样本数据的支持，不能满足实际应用的需要，因此SAR图像的统计建模

一般采用参量模型的形式。

海杂波建模的一般过程

clear;clc;
fid=fopen('D:\Ship\Clutter\2','rb');
width=512;
height=512;
imgdata=fread(fid,[width,height],'uint16');
fclose(fid);

L=2; %视数
imgdata=double(imgdata);
imgdata2=imgdata.^2;
mean1=mean(mean(imgdata));
mean2=mean(mean(imgdata2));
vmuc=1/(mean2/mean1/mean1/(1+1/L)-1); %形状参数

y=imgdata(:);
[h,w]=hist(y,max(y));
bar(w,h/sum(h)); %直方图
hold on;

%高斯分布
[f,s]=normfit(y);
y1=normpdf(w,f,s);
[h1,p1,k1,cv1]=kstest(y,[y,normcdf(y,f,s)]); %ks检验
hold on;
%瑞利分布
[p,n]=raylfit(y,0.05);
y2=raylpdf(w,p);
[h2,p2,k2,cv2]=kstest(y,[y,raylcdf(y,p)]);  %ks检验
hold on;
%威布尔分布
a=wblfit(y);
y3=wblpdf(w,a(1),a(2));
[h3,p3,k3,cv3]=kstest(y,[y,wblcdf(y,a(1),a(2))]);  %ks检验
hold on;
%对数正态分布
par=lognfit(y);
y4=lognpdf(w,par(1),par(2));
[h4,p4,k4,cv4]=kstest(y,[y,logncdf(y,par(1),par(2))]);  %ks检验
hold on;
%Gamma分布
pha=gamfit(y);
y5=gampdf(w,pha(1),pha(2));
[h5,p5,k5,cv5]=kstest(y,[y,gamcdf(y,pha(1),pha(2))]);  %ks检验
hold on;
%K分布
if (vmuc<0)
    y6=chi2pdf(w,L); %形状因子为负数时，卡方分布
else
    %y6=2/w/gamma(vmuc)/gamma(L)*((L*vmuc*w)/mean(y)).^(L+vmuc)/2*besselk((vmuc-L),2*sqrt((L*vmuc*w)/mean(y)));
    alpha=sqrt(std(y).^2/(2*vmuc));
    y6=2*((w/(2*alpha)).^vmuc).*besselk((vmuc-1),w/alpha)./(alpha*gamma(vmuc)); %K分布
end
[h6,p6,f6]=kstest2(y,y6);
%画图
plot(w,y1,'r-',w,y2,'k:',w,y3,'m-.',w,y4,'y--',w,y5,'c',w,y6,'k');
legend('直方图','Normal','Rayleigh','Weibull','Lognormal','Gamma','K')
title('直方图拟合PDF');

xlabel('杂波像素值');ylabel('概率密度值');
disp('峰度系数：'); 
K= kurtosis(y)-3 
disp('偏度系数：'); 
S=skewness(y)
disp('方差/均值：'); 
A=std(y)/mean(y) 
disp('临界值：'); 
P=[cv1,cv2,cv3,cv4,cv5]
disp('检验统计量：'); 
F=[k1,k2,k3,k4,k5]

杂波拟合结果

function value = K_Func(x,L,u,v)
%K-distribution function
%   Detailed explanation goes here
% L 视数
% u 均值
% v 形状参数
if(x==0)
    value=0;
else
    t=sqrt(L*v*x/u);
    alpha=v-L;
    lamda=v+L-1;
    a=(lamda+alpha-1)/(lamda-alpha+1)/(lamda+alpha+1)*t^(lamda+1)*besselk(alpha,t)*GeneHypFunc(1,(lamda-alpha+3)/2,(lamda+alpha+1)/2,t*t/4);
    b=1/(lamda-alpha+1)/(lamda+alpha+1)*t^(lamda+2)*besselk(alpha-1,t)*GeneHypFunc(1,(lamda-alpha+3)/2,(lamda+alpha+3)/2,t*t/4);
    value=a+b;
end
end

function value = GeneHypFunc( a,b,c,z )
%GeneHypFunc 广义超几何函数
%   Detailed explanation goes here
eps=1e-7;
sum=0.0;
k=0;
t=gamma(b)*gamma(c)*gamma(a)/gamma(a)/gamma(b)/gamma(c);
while(t>eps)
    sum=sum+t;
    k=k+1;
    t=gamma(b)*gamma(c)*gamma(a+k)/gamma(a)/gamma(b+k)/gamma(c+k)*(z^k)/factorial(k);
end
value=sum;
end