编程实现基于 NDVI 的植被变化检测

说明

编程实现植被变化检测的四个关键环节：

影像的显示（432 波段组合， 16bit 变为 8bit）

影像的相对配准

NDVI 计算及专题图生成

植被变化区域 显示及变化面积统计

---

关键

1：envi格式影像读取显示，灰度拉伸以增强显示效果，网上查找到的enviread函数读取，稍加修改，依据辅助文件来读取影像信息。

2.相对配准同名点获取：

自动提取：sift或者orb等特征提取，优质影像对可以提取提取很多同名点，质量不好提取很少。

IMG1=im2double(img1_8bit);
IMG2=im2double(img2_8bit);
A1=detectSIFTFeatures(img1_8bit(:,:,4));
A2=detectSIFTFeatures(img2_8bit(:,:,4));
[B1,A1]=extractFeatures(img1_8bit(:,:,4),A1);
[B2,A2]=extractFeatures(img2_8bit(:,:,4),A2);
boxPairs = matchFeatures(B1, B2);
movingPoints=A1(boxPairs(:,1),:);
fixedPoints=A2(boxPairs(:,2),:);

手动提取，人工刺点cpselect函数

代码

clear all;
clc;
%% 数据读取
[img1,info1]=enviread('Q1','Q1.HDR');
[img2,info2]=enviread('SV1_5','SV1_5.HDR');
%% 图像显示
%16位图像转8位图
%quickbird
img1_guiyi=mat2gray(img1);
img1_8bit=im2uint8(img1_guiyi);
%SV1
img2_guiyi=mat2gray(img2);
img2_8bit=im2uint8(img2_guiyi);
%直方图均衡化，便于展示
img1_junhenghua=histeq(img1_guiyi);
img1_show=im2uint8(img1_junhenghua);
img2_junhenghua=histeq(img2_guiyi);
img2_show=im2uint8(img2_junhenghua);
%quickbird
figure(1)
imgmul1=cat(3,img1_show(:,:,4),img1_show(:,:,3),img1_show(:,:,2));%合成3维矩阵
imshow(imgmul1);
title('Quickbird后三个波段组合图像');
%SV1
figure(2)
imgmul2=cat(3,img2_show(:,:,4),img2_show(:,:,3),img2_show(:,:,2));%合成3维矩阵
imshow(imgmul2);
title('SV1后三个波段组合图像');
%% 相对配准
% 选择匹配点对
[fixedPoints, movingPoints] = cpselect(img1_8bit(:,:,4), img2_8bit(:,:,4), 'Wait', true);
% 使用fitgeotrans函数得到二次多项式变换矩阵Tpoly
Tpoly = fitgeotrans(movingPoints, fixedPoints, 'polynomial', 2);
% 利用变换矩阵Tpoly进行图像配准
movingRegisteredPoly = imwarp(img2_8bit, Tpoly, 'OutputView', imref2d(size(img1_8bit(:,:,1))));
% 显示结果
%直方图均衡化，便于展示
movingRegisteredPoly_guiyi=mat2gray(movingRegisteredPoly);
movingRegisteredPoly_junhenghua=histeq(movingRegisteredPoly_guiyi);
movingRegisteredPoly_show=im2uint8(movingRegisteredPoly_junhenghua);
figure(3)
imshowpair(img1_show(:,:,4:-1:2), movingRegisteredPoly_show(:,:,4:-1:2))
title('二次多项式纠正叠加图')
figure(4)
imshow(movingRegisteredPoly_show(:,:,4:-1:2))
title('二次多项式纠正结果图')
%精度评定
%内精度
f1=transformPointsInverse(Tpoly,fixedPoints);
Vx=f1(:,1)-movingPoints(:,1);
Vy=f1(:,2)-movingPoints(:,2);
dx1=sqrt(sum(Vx.^2)/(7-6));
dy1=sqrt(sum(Vy.^2)/(7-6));
%外精度
% 刺新的点做精度评定
[fixedPoints_jc, movingPoints_jc] = cpselect(img1_8bit(:,:,4), img2_8bit(:,:,4), 'Wait', true);
f2=transformPointsInverse(Tpoly,fixedPoints_jc);
dx=movingPoints_jc(:,1)-f2(:,1);
dy=movingPoints_jc(:,2)-f2(:,2);
%% NDVI 专题图
img1_double=im2double(img1_8bit);
img2_double=im2double(movingRegisteredPoly);
for i=1:size(img1_8bit(:,:,1),1)
    for j=1:size(img1_8bit(:,:,1),2)
        NDVI1(i,j)=(img1_double(i,j,4)-img1_double(i,j,3))/(img1_double(i,j,4)+img1_double(i,j,3));
        NDVI2(i,j)=(img2_double(i,j,4)-img2_double(i,j,3))/(img2_double(i,j,4)+img2_double(i,j,3));
        %quickbird
        if NDVI1(i,j)<=0.1%非植被覆盖区,黑色表示
            I1(i,j,1)=0;I1(i,j,2)=0;I1(i,j,3)=0;
        end
        if NDVI1(i,j)>0.1&&NDVI1(i,j)<=0.25%低植被覆盖区，红色
            I1(i,j,1)=255;I1(i,j,2)=0;I1(i,j,3)=0;
        end
        if NDVI1(i,j)>0.25&&NDVI1(i,j)<=0.4%中植被覆盖区，蓝色
            I1(i,j,1)=0;I1(i,j,2)=0;I1(i,j,3)=255;
        end
        if NDVI1(i,j)>0.4%高植被覆盖区，绿色
            I1(i,j,1)=0;I1(i,j,2)=255;I1(i,j,3)=0;
        end
        %SV1
        if NDVI2(i,j)<=0.1%非植被覆盖区,黑色表示
            I2(i,j,1)=0;I2(i,j,2)=0;I2(i,j,3)=0;
        end
        if NDVI2(i,j)>0.1&&NDVI2(i,j)<=0.25%低植被覆盖区，红色
            I2(i,j,1)=255;I2(i,j,2)=0;I2(i,j,3)=0;
        end
        if NDVI2(i,j)>0.25&&NDVI2(i,j)<=0.4%中植被覆盖区，蓝色
            I2(i,j,1)=0;I2(i,j,2)=0;I2(i,j,3)=255;
        end
        if NDVI2(i,j)>0.4%高植被覆盖区，绿色
            I2(i,j,1)=0;I2(i,j,2)=255;I2(i,j,3)=0;
        end
    end
end
%显示
figure(5)
imshow(I1)
title('quickbird的NDVI图')
figure(6)
imshow(I2)
title('SV1的NDVI图')
%% 植被变化结果图
sum_up=0;
sum_down=0;
sum_un=0;
sum_no=0;
for i=1:size(img1_8bit(:,:,1),1)
    for j=1:size(img1_8bit(:,:,1),2)
        if NDVI1(i,j)>0.1||NDVI2(i,j)>0.1
              NDVI_change(i,j)=NDVI2(i,j)-NDVI1(i,j);
              if NDVI_change(i,j)<-0.1%减少,红
                  I_change(i,j,1)=1;I_change(i,j,2)=0;I_change(i,j,3)=0;
                  sum_down=sum_down+1;
              end
              if NDVI_change(i,j)<=0.1&&NDVI_change(i,j)>=-0.1%不变，深绿
                  I_change(i,j,1)=0;I_change(i,j,2)=0.5;I_change(i,j,3)=0;
                  sum_un=sum_un+1;
              end
              if NDVI_change(i,j)>0.1%增加,浅绿
                  I_change(i,j,1)=0.2;I_change(i,j,2)=0.8;I_change(i,j,3)=0.2;
                  sum_up=sum_up+1;
              end
        else %非植被区，灰色
            I_change(i,j,1)=0.37;I_change(i,j,2)=0.37;I_change(i,j,3)=0.37;
            sum_no=sum_no+1;
        end
    end
end
%面积
area_up=0.25*sum_up;
area_down=0.25*sum_down;
area_un=0.25*sum_un;
area_no=0.25*sum_no;
% 画图
figure(7)
Img_change=im2uint8(I_change);
imshow(Img_change)
title('变化图')