bp.m

利用bp神经网络对遥感图像进行分类

clc;
clear;



%处理农田样本
crI=imread('CR.jpg');
crr=crI(:,:,1);crg=crI(:,:,2);crb=crI(:,:,3);%提取图像三个波段
crrr=im2double(crr);crgg=im2double(crg);crbb=im2double(crb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
crro=[];
crgo=[];
crbo=[];
[crlen crlen2]=size(crrr);%计算rr矩阵的宽度和长度
crx=[];

%将矩阵转换成一维的矩阵
for i=1:crlen
    crro=[crro,crrr(i,:)];   %cat
end
for i=1:crlen
    crgo=[crgo,crgg(i,:)];
end
for i=1:crlen
    crbo=[crbo,crbb(i,:)];
end

crx=[crro;crgo;crbo];    %将三个波段的矩阵合并为一个矩阵


crlength=length(crx);

%处理海水样本
seaI=imread('sea.jpg');
sear=seaI(:,:,1);seag=seaI(:,:,2);seab=seaI(:,:,3);%提取图像三个波段
searr=im2double(sear);seagg=im2double(seag);seabb=im2double(seab);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
searo=[];
seago=[];
seabo=[];
[sealen sealen2]=size(searr);%计算rr矩阵的宽度和长度
seax=[];
seax1=[];
%将300*300的矩阵转换成一维的矩阵
for i=1:sealen
    searo=[searo,searr(i,:)];   %cat
end
for i=1:sealen
    seago=[seago,seagg(i,:)];
end
for i=1:sealen
    seabo=[seabo,seabb(i,:)];
end

seax=[searo;seago;seabo];    %将三个波段的矩阵合并为一个矩阵


sealength=length(seax);

%处理芦苇样本
reedI=imread('reed.jpg');
reedr=reedI(:,:,1);reedg=reedI(:,:,2);reedb=reedI(:,:,3);%提取图像三个波段
reedrr=im2double(reedr);reedgg=im2double(reedg);reedbb=im2double(reedb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
reedro=[];
reedgo=[];
reedbo=[];
[reedlen reedlen2]=size(reedrr);%计算rr矩阵的宽度和长度
reedx=[];
reedx1=[];
%将矩阵转换成一维的矩阵
for i=1:reedlen
    reedro=[reedro,reedrr(i,:)];   %cat
end
for i=1:reedlen
    reedgo=[reedgo,reedgg(i,:)];
end
for i=1:reedlen
    reedbo=[reedbo,reedbb(i,:)];
end

reedx=[reedro;reedgo;reedbo];    %将三个波段的矩阵合并为一个矩阵


reedlength=length(reedx);

%处理绿地样本
grI=imread('GR.jpg');
grr=grI(:,:,1);grg=grI(:,:,2);grb=grI(:,:,3);%提取图像三个波段
grrr=im2double(grr);grgg=im2double(grg);grbb=im2double(grb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
grro=[];
grgo=[];
grbo=[];
[grlen grlen2]=size(grrr);%计算rr矩阵的宽度和长度
grx=[];
grx1=[];
%将矩阵转换成一维的矩阵
for i=1:grlen
    grro=[grro,grrr(i,:)];   %cat
end
for i=1:grlen
    grgo=[grgo,grgg(i,:)];
end
for i=1:grlen
    grbo=[grbo,grbb(i,:)];
end

grx=[grro;grgo;grbo];    %将三个波段的矩阵合并为一个矩阵


grlength=length(grx);

%处理养殖场样本
culI=imread('cul.jpg');
culr=culI(:,:,1);culg=culI(:,:,2);culb=culI(:,:,3);%提取图像三个波段
culrr=im2double(culr);culgg=im2double(culg);culbb=im2double(culb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
culro=[];
culgo=[];
culbo=[];
[cullen cullen2]=size(culrr);%计算rr矩阵的宽度和长度
culx=[];
culx1=[];
%将矩阵转换成一维的矩阵
for i=1:cullen
    culro=[culro,culrr(i,:)];   %cat
end
for i=1:cullen
    culgo=[culgo,culgg(i,:)];
end
for i=1:cullen
    culbo=[culbo,culbb(i,:)];
end

culx=[culro;culgo;culbo];    %将三个波段的矩阵合并为一个矩阵


cullength=length(culx);

%处理滩涂样本
mudI=imread('mudflat.jpg');
mudr=mudI(:,:,1);mudg=mudI(:,:,2);mudb=mudI(:,:,3);%提取图像三个波段
mudrr=im2double(mudr);mudgg=im2double(mudg);mudbb=im2double(mudb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
mudro=[];
mudgo=[];
mudbo=[];
[mudlen mudlen2]=size(mudrr);%计算rr矩阵的宽度和长度
mudx=[];
mudx1=[];
%将矩阵转换成一维的矩阵
for i=1:mudlen
    mudro=[mudro,mudrr(i,:)];   %cat
end
for i=1:mudlen
    mudgo=[mudgo,mudgg(i,:)];
end
for i=1:mudlen
    mudbo=[mudbo,mudbb(i,:)];
end

mudx=[mudro;mudgo;mudbo];    %将三个波段的矩阵合并为一个矩阵


mudlength=length(mudx);

%处理互米花草样本
huI=imread('humi.jpg');
hur=huI(:,:,1);hug=huI(:,:,2);hub=huI(:,:,3);%提取图像三个波段
hurr=im2double(hur);hugg=im2double(hug);hubb=im2double(hub);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
huro=[];
hugo=[];
hubo=[];
[hulen hulen2]=size(hurr);%计算rr矩阵的宽度和长度
hulength=hulen*hulen2
hux=[];
hux1=[];
%将矩阵转换成一维的矩阵
for i=1:hulen
    huro=[huro,hurr(i,:)];   %cat
end
for i=1:hulen
    hugo=[hugo,hugg(i,:)];
end
for i=1:hulen
    hubo=[hubo,hubb(i,:)];
end

hux=[huro;hugo;hubo];    %将三个波段的矩阵合并为一个矩阵


hulength=length(hux);

%处理房屋样本
fwI=imread('fw.jpg');
fwr=fwI(:,:,1);fwg=fwI(:,:,2);fwb=fwI(:,:,3);%提取图像三个波段
fwrr=im2double(fwr);fwgg=im2double(fwg);fwbb=im2double(fwb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
fwro=[];
fwgo=[];
fwbo=[];
[fwlen fwlen2]=size(fwrr);%计算rr矩阵的宽度和长度
fwx=[];
fwx1=[];
%将矩阵转换成一维的矩阵
for i=1:fwlen
    fwro=[fwro,fwrr(i,:)];   %cat
end
for i=1:fwlen
    fwgo=[fwgo,fwgg(i,:)];
end
for i=1:fwlen
    fwbo=[fwbo,fwbb(i,:)];
end
fwx=[fwro;fwgo;fwbo];    %将三个波段的矩阵合并为一个矩阵


fwlength=length(fwx);

%处理海三棱藨草的样本
sanI=imread('san.jpg');
sanr=sanI(:,:,1);sang=sanI(:,:,2);sanb=sanI(:,:,3);%提取图像三个波段
sanrr=im2double(sanr);sangg=im2double(sang);sanbb=im2double(sanb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
sanro=[];
sango=[];
sanbo=[];
[sanlen sanlen2]=size(sanrr);%计算rr矩阵的宽度和长度
sanx=[];
sanx1=[];
%将矩阵转换成一维的矩阵
for i=1:sanlen
    sanro=[sanro,sanrr(i,:)];   %cat
end
for i=1:sanlen
    sango=[sango,sangg(i,:)];
end
for i=1:sanlen
    sanbo=[sanbo,sanbb(i,:)];
end

sanx=[sanro;sango;sanbo];    %将三个波段的矩阵合并为一个矩阵


sanlength=length(sanx);

%处理未利用土地的样本
unI=imread('un.jpg');
unr=unI(:,:,1);ung=unI(:,:,2);unb=unI(:,:,3);%提取图像三个波段
unrr=im2double(unr);ungg=im2double(ung);unbb=im2double(unb);%将黑白图像转换成双精度格式并灰度值归一化为[0,1]之间
unro=[];
ungo=[];
unbo=[];
[unlen unlen2]=size(unrr);%计算rr矩阵的宽度和长度
unx=[];
unx1=[];
%将矩阵转换成一维的矩阵
for i=1:unlen
    unro=[unro,unrr(i,:)];   %cat
end
for i=1:unlen
    ungo=[ungo,ungg(i,:)];
end
for i=1:unlen
    unbo=[unbo,unbb(i,:)];
end

unx=[unro;ungo;unbo];   %将三个波段的矩阵合并为一个矩阵

unlength=length(unx);

%将样本置于一个一维向量中
%
% P_G=[crgo,seago,reedgo,grgo,culgo,mudgo,hugo,fwgo,sango,ungo];
% P_B=[crbo,seabo,reedbo,grbo,culbo,mudbo,hubo,fwbo,sanbo,unbo];
% P=[P_RP_G,P_B];
 P=[crx,seax,reedx,grx,culx,mudx,hux,fwx,sanx,unx];



%设定目标向量，一维
T=[];
for i=1:100
    T(1,i)=0.1;
    T(2,i)=0.9;
    T(3,i)=0.1;
end
for i=101:200
    T(1,i)=0.2;
    T(2,i)=0.2;
    T(3,i)=0.9;
end
for i=201:300
    T(1,i)=0.9;
    T(2,i)=0.5;
    T(3,i)=0.2;
end
for i=301:400
    T(1,i)=0.6;
    T(2,i)=0.5;
    T(3,i)=0.7;
end
for i=401:464
    T(1,i)=0.3;
    T(2,i)=0.1;
    T(3,i)=0.1;
end
for i=465:564
    T(1,i)=0.8;
    T(2,i)=0.5;
    T(3,i)=0.1;
end
for i=565:664
    T(1,i)=0.4;
    T(2,i)=0.1;
    T(3,i)=0.7;
end
for i=665:764
    T(1,i)=0.5;
    T(2,i)=0.9;
    T(3,i)=0.9;
end
for i=765:864
    T(1,i)=0.25;
    T(2,i)=0.6;
    T(3,i)=0.9;
end
for i=865:964
    T(1,i)=0.7;
    T(2,i)=0.5;
    T(3,i)=0.1;
end


clc;
net=newcf(minmax(P),[15,3],{'tansig','purelin'},'traingdx');
%当前输入层权值和阈值
imputweights=net.iw{1,1}
inputbias=net.b{1}  %产生隐含层各神经元的阈值
%当前网路层权值和阈值 
layerweights=net.lw{2,1}
layerbias=net.b{2}

%设置训练参数
net.trainparam.show=50 %迭代50次显示一次
net.trainparam.lr=0.05 %学习率为0.05（一般大小在0-1之间）
net.trainparam.mc=0.9;  %动量因子0.9（0-1之间）
net.trainparam.epochs=1000;%迭代次数1000次
net.trainparam.goal=1e-3;%误差
clc
%调用traingdm算法训练bp网络
newnet=train(net,P,T);
clc
%仿真
I=imread('tm2000mask.jpg');

ro=I(:,:,1);
go=I(:,:,2); 
bo=I(:,:,3); 
rr=im2double(ro);%将r,g,b图像变为双精度格式，其灰度值归一到0，1之间
gg=im2double(go);
bb=im2double(bo);
[m,n]=size(rr)
%将图像变为一维向量形式
x=[];
for i=1:m
    x=[x,rr(i,:)];
end

y=[];
for i=1:m
    y=[y,gg(i,:)];
end

z=[];
for i=1:m
    z=[z,bb(i,:)];
end
tm=[x;y;z];
a=sim(newnet,tm);


rR=[];
R=[a(1,:)];
for i=0:779
   rR=[rR;R((i*1090+1):(i*1090)+1090)];
end



rG=[];
G=[a(2,:)];
for i=0:779
   rG=[rG;G((i*1090+1):(i*1090)+1090)];
end



rB=[];
B=[a(3,:)];
for i=0:779
   rB=[rB;B((i*1090+1):(i*1090)+1090)];
end
cr=0;
sea=0;
reed=0;
gr=0;
cul=0;
mud=0;
hu=0;
fw=0;
san=0;
un=0;
for i=1:780
 for j=1:1090
    if ((0.098<rR(i,j))&(rR(i,j)<0.12))
        cr=cr+1;
    else if ((0.19<rR(i,j))&(rR(i,j)<0.22))
            sea=sea+1;
        else if ((0.88<rR(i,j))&(rR(i,j)<0.92))
            reed=reed+1;
            else if ((0.58<rR(i,j))&(rR(i,j)<0.62))
                    gr=gr+1;
                else if ((0.28<rR(i,j))&(rR(i,j)<0.32))
                        cul=cul+1;
                    else if ((0.78<rR(i,j))&(rR(i,j)<0.82))
                            mud=mud+1;
                        else if ((0.38<rR(i,j))&(rR(i,j)<0.42))
                                hu=hu+1;
                            else if ((0.48<rR(i,j))&(rR(i,j)<0.52))
                                    fw=fw+1;
                                else if ((0.24<rR(i,j))&(rR(i,j)<0.26))
                                        san=san+1;
                                    else if ((0.68<rR(i,j))&(rR(i,j)<0.78))
                                            un=un+1;
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
 end
end
figure,pie([cr,sea,reed,gr,cul,mud,hu,fw,san,un]);
title('各类型土地利用比例图');
legend('农田','海水','芦苇','绿地','养殖场','滩涂','互米花草','房屋','海三棱藨草','未利用地');

%将数据转成0-255格式
R1=abs(rR)*255;
G1=abs(rG)*255;
B1=abs(rB)*255;
%转成彩图RGB格式显示
R2=uint8(R1);
G2=uint8(G1);
B2=uint8(B1);
RGB=cat(3,R2,G2,B2);
figure,imagesc(RGB);


disp('1--分类后红色波段图像');
disp('2--分类后绿色波段图像');
disp('3--分类后蓝色波段图像');
choice=input('请输入数字，对应要显示的图像：');
if choice==1
    figure,imshow(abs(rR))
else if choice==2
        figure,imshow(abs(rG))
    else if choice==3
            figure,imshow(abs(rB))
        end
    end
end