基于Matlab的遥感图像变化监测系统

目 录

1. 研究背景和文献综述 1
2. 实验流程 3
3. 数据预处理 4
   3.1. 大气纠正 4
   3.2. 几何纠正/图像配准 4
   3.2.1. 裁件拼接 4
   3.2.2. 辐射度匹配和归一化处理 4
   3.2.3. 灰度直方图匹配算法 5
   3.2.4. SCR图像匹配算法 5
4. 基于多时相数据图像灰度差值法的地表变化检测 6
   4.1. 灰度差值法的原理 6
   4.2. 图像差值运算 6
   4.3. 差值运算后处理 6
   4.4. 基于NDVI数据的差值法进行多时相数据检测 9
5. 基于多时相数据叠合方法的变化检测 10
   5.1. 基于不同时相的波段组合 10
   5.2. 基于NDVI数据的波段组合 12
6. 基于变化向量分析的地表变化检测 13
   6.1. 实验原理 13
   6.1.1. 变化向量分析（CVA） 13
   6.1.2. CVA的变化强度计算 13
   6.1.3. CVA的变化类型确定方法 14
   6.1.4. 基于对象分类的变化向量（CVA）分析法 14
   6.1.5. 变化阈值的确定 14
   6.2. 基于变化向量分析的地表变化检测（NDVI BI数据） 18
   6.3. 基于CVA和双窗口变步长阈值提取的变化监测 20
   6.3.1. 计算变化向量的变化强度 20
   6.3.2. 计算变化向量的变化方向 25
7. 总结与反思 26
   7.1. 变化检测方法的发展历程 26
   7.2. 对变化检测方法的评价 26
   7.3. 对于阈值提取方法的评价 27
8. 参考文献: 28
   2.实验流程
   经过文献调研，我发现变化检测需要经过如下流程：
   1.数据预处理，包括对于数据的可用性和质量的评价，对于缺失数据的填补，对于不同时相数据做大气校正、几何精校正或者是相对配准。保证卫星图像的清晰、没有云雾的遮挡，分析前的几何精纠正（正射纠正），使得总误差小于一个像元，某些情况下，对于高分辨率的影响，为了减少数据量，可以将不同时相的数据进行重采样，需要将最后的RMS配准结果控制在0.5像素一下。如果只有两个时相的数据可以只进行相对配准，而如果出现多幅图像，则需要统一配准到坐标系里面。
   2.变化检测，包括使用基于像元的分析方法和基于对象的分析方法，分析不同时相数据的变化情况。将结果进一步使用Arcgis和Envi进行分级操作。
   3.阈值确定，对于第二步所获取的结果，人工或者使用相关算法确定较为合理的阈值。
   4.检测结果的精度评价和精度验证，本文转载自http://www.biyezuopin.vip/onews.asp?id=14929一般情况下需要野外采点获取数据，也可以基于已有的矢量数据，进行评价。或者基于高光谱的数据验证分类之后的结果。
   3.数据预处理
   3.1.大气纠正
   可以采用FLAASH模型
   3.2.几何纠正/图像配准
   可以采用历史TM影像（具备准确的空间位置，经野外 GPS 验证）作为基准影像，配准方法为多项式变换，配准误差控制在 0.5 个像元以内。
   多时相图像的几何配准是指多时相图像的同名像点互相重叠，将一幅图像作为基准图像(晚期图像)，将一张图像作为待配准图像（早期图像），通过选取的控制点采用多项式的方法确定基准图像与待配准图像之间的对应关系，进行几何精纠正（正射纠正），基于三次卷积运算，实现图像的配准，配准的精度至少要小于一个像元。
   3.2.1.裁件拼接
   如果实验中需要涉及到不同的图幅共同构成变化区域，我们还将需要对于不同图幅的图像进行裁剪拼接。
   3.2.2.辐射度匹配和归一化处理
   利用多波段遥感数据进行生态环境和土地资源变化识别和监测时, 因多时相数据集受到季节性地物辐射变化、 太阳光照条件差异、 气象条件 ( 大气散射、 吸收和云量变化等) 波动大等因素影响需要进行辐射的相对归一化处理。[17]
   传统的归一化方法是采用基于全景或波谱稳定的子集的统计参数方法, 如最大- 最小 (MM) 归一化法、 平均标准方差法 (MS) 、全景简单线性回归法 (SR) 、 直方图匹配 (HM)。

%% 读取数据
% clc 
% close all
% clear all
path ='C:\Users\Ren\Desktop\变化检测作业';
userpath(path);
[image1,p1,t1]=freadenvinew('C:\Users\Ren\Desktop\变化检测作业\first.img');
[image3,p3,t3] = freadenvinew('C:\Users\Ren\Desktop\变化检测作业\second_cut.img');

%% 进行数据的重新的分割
%使用image无法分割三个影像，关键问题就是，读取进来的东西是BSQ格式。。。
first_image = reshape(image1,[608,1676,6]);
% second_image = reshape(image3,[608,1676,6]);
second_image = reshape(image3,[608,1676,6]);
figure;imshow(first_image(:,:,1),[]);
figure;imshow(second_image(:,:,1),[]);
%% 假设配准之后，像元严格对应。




%%
% 问题就是reshape之后就变质了，数据变成了uint8
% figure;imshow(first_image(:,:,1),[]);title('第一幅影像图first image');
% figure;imshow(second_image(:,:,1),[]);title('第二幅影像图 进行几何校正之后 second image');
% 然后就是matlab和envi里面的数据是相反的排列方式。。。
% first_image2 = reshape(image1,[1676,608,6]);
% figure;imshow(first_image2(:,:,1),[]);

% image4 =third_image(:,:,1) - first_image(:,:,1);
% 因为matlab中矩阵维度的不一致，所以我们没有一个好的办法实现图像的。。

% 两个方式：一个是用matlab的图像配准工具箱
% 一个是把图像裁剪到相同的大小，因为现在图像的区域是对应的，所以裁剪可行。
%但是我倾向于使用第一种方式也就是图像配准工具箱的方式。
% 但是图像配准工具箱要首先输出图像
%% 输出图像
weight = 50.0;
first_image50 = (uint16(first_image.*50.0));
second_image50  = (uint16(second_image.*50.0));
% imwrite(first_image50(:,:,1),[path,'firstB1','.tiff'],'tiff','WriteMode','overwrite','ColorSpace','icclab','Compression','none','Resolution',144);
%% 全路径输出文件
imwrite(first_image50(:,:,1),[path,'firstB1','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
imwrite(second_image50(:,:,1),[path,'secondB1','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(first_image50(:,:,2),[path,'firstB2','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(second_image50(:,:,2),[path,'secondB2','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(first_image50(:,:,3),[path,'firstB3','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(second_image50(:,:,3),[path,'secondB3','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(first_image50(:,:,4),[path,'firstB4','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(second_image50(:,:,4),[path,'secondB4','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(first_image50(:,:,5),[path,'firstB5','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(second_image50(:,:,5),[path,'secondB5','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(first_image50(:,:,6),[path,'firstB7','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);
% imwrite(second_image50(:,:,6),[path,'secondB7','.tiff'],'tiff','WriteMode','overwrite','Compression','none','Resolution',144);


%% 读取输出的文件
firstB1 =double(imread('firstB1.tif'))./weight;
secondB1=double(imread('secondB1.tif'))./weight;
% firstB2 =double(imread('firstB2.tif'))./weight;
% secondB2=double(imread('secondB2.tif'))./weight;
% firstB3 =double(imread('firstB3.tif'))./weight;
% secondB3=double(imread('secondB3.tif'))./weight;
% firstB4 =double(imread('firstB4.tif'))./weight;
% secondB4=double(imread('secondB4.tif'))./weight;
% firstB5 =double(imread('firstB5.tif'))./weight;
% secondB5=double(imread('secondB5.tif'))./weight;
% firstB7 =double(imread('firstB7.tif'))./weight;
% secondB7=double(imread('secondB7.tif'))./weight;

%% 进行赋值操作
fixed_image1  = firstB1;
% fixed_image2  = firstB2;
% fixed_image3  = firstB3;
% fixed_image4  = firstB4;
% fixed_image5  = firstB5;
% fixed_image7  = firstB7;
moving_image1 = secondB1;
% moving_image2 = secondB2;
% moving_image3 = secondB3;
% moving_image4 = secondB4;
% moving_image5 = secondB5;
% moving_image7 = secondB7;

%% 进行自动配准操作
[secondB1_estimtor_result] = registerImagesMULTIMODAL(moving_image1,fixed_image1);
% [secondB2_estimtor_result] = registerImagesMULTIMODAL(moving_image2,fixed_image2);
% [secondB3_estimtor_result] = registerImagesMULTIMODAL(moving_image3,fixed_image3);
% [secondB4_estimtor_result] = registerImagesMULTIMODAL(moving_image4,fixed_image4);
% [secondB5_estimtor_result] = registerImagesMULTIMODAL(moving_image5,fixed_image5);
% [secondB7_estimtor_result] = registerImagesMULTIMODAL(moving_image7,fixed_image7);
%% 输出图像
secondB1_fixed = secondB1_estimtor_result.RegisteredImage;
% secondB2_fixed = secondB2_estimtor_result.RegisteredImage;
% secondB3_fixed = secondB3_estimtor_result.RegisteredImage;
% secondB4_fixed = secondB4_estimtor_result.RegisteredImage;
% secondB5_fixed = secondB5_estimtor_result.RegisteredImage;
% secondB7_fixed = secondB7_estimtor_result.RegisteredImage;

%% 拼成一个矩阵
second_image(:,:,1) = secondB1_fixed;
% second_image(:,:,2) = secondB2_fixed;
% second_image(:,:,3) = secondB3_fixed;
% second_image(:,:,4) = secondB4_fixed;
% second_image(:,:,5) = secondB5_fixed;
% second_image(:,:,6) = secondB7_fixed;
%% 图像读取进来了，进行变化检测的部分了。



%% 参考文献
% BIL、BIP 和 BSQ 栅格文件
% http://desktop.arcgis.com/zh-cn/arcmap/10.3/manage-data/raster-and-images/bil-bip-and-bsq-raster-files.htm