多光谱数据的贝叶斯分类

1.原理介绍

1.1最佳统计分类器

来自类 ωj 的特定模式概率表示为 p(ωj/x)。如果模式分类器判断 x 来自类 ωj，而实际它来自类 ωi，
那么分类器就会导致一次损失，表示为 Lij。由于模式 x 可能来自 W 中的任何一个类，故将模式 x 判
决为 ωj 的平均损失为：

根据贝叶斯公式， p(A/B) = [P(A)P(B/A)]/P(B)，所以上式可以写为：

经过讨论，我们知道 0-1 损失函数的贝叶斯分类器的决策函数的计算：

1.2高斯模式类的贝叶斯分类器

高斯分类器由样本均值向量和协方差矩阵决定，均值向量和协方差的估计如下：

和

根据前述的贝叶斯分类器决策函数，高斯模式下的贝叶斯分类器决策函数为：

如果所有的协方差矩阵都相等，则：

这是一个线性决策函数（超平面），其中 j = 1, 2, . . . , W。进一步如果 C = I，这里 I 为单位矩阵，且
P(ωj) = 1/W, j = 1, 2, . . . , W，则有：

2.多光谱图像贝叶斯分类

地面上的每个点都可以形如 x = (x1, x2, x3, x4)T 的一个四维模式向量表示，其中 x1 为蓝光图像，
x2 为绿光图像等等。如实验原理所述，高斯模式分类器要求估计每个类的均值向量和协方差矩阵。
根据前述高斯模式类贝叶斯分类器原理，完成该实验分为如下步骤：

2.1分类结果

2.2MATLAB代码

主函数：

clc;
clear;
close all;
%读入蓝光、绿光、红光与红外波段图像
imga=imread(['blue', '.tif']);
[M,N]=size(imga);
imgb=imread(['green', '.tif']);
imgc=imread(['red', '.tif']);
imgd=imread(['outside', '.tif']);
imge=imread(['model', '.tif']);
%样本模板图像
water=imread(['water', '.tif']);
city=imread(['city', '.tif']);
tree=imread(['tree', '.tif']);
%样本数
waterN=sum(sum(water));%967
cityN=sum(sum(city));%1865
treeN=sum(sum(tree));%965
%%%%%%%%训练过程 计算每个类均值方差
[waterm,waterC]=junzhifangc(M,N,imga,imgb,imgc,imgd,water,waterN);
[citym,cityC]=junzhifangc(M,N,imga,imgb,imgc,imgd,city,cityN);
[treem,treeC]=junzhifangc(M,N,imga,imgb,imgc,imgd,tree,treeN);

%%%%%%%%%分类为水体%%%%%%%%%
waterdiv=zeros(M,N);
for i=1:M
    for j=1:N
        x=[imga(i,j),imgb(i,j),imgc(i,j),imgd(i,j)]';
        a=panjue(x,waterm,waterC,citym,cityC,treem,treeC);
        if a==1 %对应水体
            waterdiv(i,j)=255;
        end
    end
end
% imshow(waterdiv);
% imwrite(waterdiv, sprintf('result/%s.jpg','waterdiv'));

%%%%%%%%%分类为市区%%%%%%%%
citydiv=zeros(M,N);
for i=1:M
    for j=1:N
        x=[imga(i,j),imgb(i,j),imgc(i,j),imgd(i,j)]';
        a=panjue(x,waterm,waterC,citym,cityC,treem,treeC);
        if a==2 %对应城市
            citydiv(i,j)=255;
        end
    end
end
% imshow(citydiv);
% imwrite(citydiv, sprintf('result/%s.jpg','citydiv'));

%%%%%%%%%分类为植被%%%%%%%%
treediv=zeros(M,N);
for i=1:M
    for j=1:N
        x=[imga(i,j),imgb(i,j),imgc(i,j),imgd(i,j)]';
        a=panjue(x,waterm,waterC,citym,cityC,treem,treeC);
        if a==3 %对应植被
            treediv(i,j)=255;
        end
    end
end
% imshow(treediv);
% imwrite(treediv, sprintf('result/%s.jpg','treediv'));

%%%%%%%%%正确分类的结果%%%%%%%%
% waterdiv=waterdiv.*water;%水体正确分类结果
% citydiv=citydiv.*city;%市区正确分类结果
% treediv=treediv.*tree;%植被正确分类结果
divresult=zeros(M,N);%正确分类的结果
divresult=waterdiv.*water+citydiv.*city+treediv.*tree;
% imshow(divresult);
% imwrite(divresult, sprintf('result/%s.jpg','divresult'));

%%%%%%%%%%输出结果%%%%%%%%%%%
subplot(331); imshow(imga); title('可见蓝光'); axis on
subplot(332); imshow(imgb); title('可见绿光'); axis on
subplot(333); imshow(imgc); title('可见红光'); axis on
subplot(334); imshow(imgd); title('近红外波长'); axis on
subplot(335); imshow(imge); title('模板图像'); axis on
subplot(336); imshow(divresult); title('分类结果'); axis on
subplot(337); imshow(waterdiv); title('水体分类'); axis on
subplot(338); imshow(citydiv); title('市区分类'); axis on
subplot(339); imshow(treediv); title('植被分类'); axis on
% imwrite(imga, sprintf('result/%s.jpg','a'));
% imwrite(imgb, sprintf('result/%s.jpg','b'));
% imwrite(imgc, sprintf('result/%s.jpg','c'));
% imwrite(imgd, sprintf('result/%s.jpg','d'));
% imwrite(imge, sprintf('result/%s.jpg','e'));

计算均值方差函数：

function [m,waterC]=junzhifangc(M,N,imga,imgb,imgc,imgd,water,waterN)
%返回类的均值与方差
waterm=zeros(4,1);
% count=0;
for i=1:M
    for j=1:N
        if water(i,j)==1
            %模式向量分别求均值
            %count=count+1;
            %if count<=waterN %一半用于训练
                waterm(1)= waterm(1)+imga(i,j)/waterN;
                waterm(2)= waterm(2)+imgb(i,j)/waterN;
                waterm(3)= waterm(3)+imgc(i,j)/waterN;
                waterm(4)= waterm(4)+imgd(i,j)/waterN;
            %end
        end
    end
end
%求协方差矩阵
waterC=zeros(4,4);
% count=0;
for i=1:M
    for j=1:N
        if water(i,j)==1
            %count=count+1;
            %if count<=waterN %一半用于训练
                x=[imga(i,j),imgb(i,j),imgc(i,j),imgd(i,j)]';%均值向量
                x=im2double(x);
                waterC=waterC+x*(x')/waterN;
            %end
        end
    end
end
m=[waterm(1),waterm(2),waterm(3),waterm(4)]';
waterC=waterC-m*m';
end

判决函数：

function a=panjue(x,waterm,waterC,citym,cityC,treem,treeC)
%返回判决结果，1，2，3对应water,city,tree
m=[waterm(1),waterm(2),waterm(3),waterm(4)]';
m=im2double(m);
x=im2double(x);
%判决函数
d1=-log(det(waterC))-0.5*(x-m)'*inv(waterC)*(x-m);
m=[citym(1),citym(2),citym(3),citym(4)]';
d2=-log(det(cityC))-0.5*(x-m)'*inv(cityC)*(x-m);
m=[treem(1),treem(2),treem(3),treem(4)]';
d3=-log(det(treeC))-0.5*(x-m)'*inv(treeC)*(x-m);

if d1>=d2&&d1>=d3%判为水体
    a=1;
end
if d2>=d1&&d2>=d3%判为市区
    a=2;
end
if d3>=d1&&d3>=d2%判为植被
    a=3;
end
end