遥感图像变化检测方法与流程

本发明属于图像处理技术领域，更进一步涉及图像分割与图像分类识别技术领域中的一种遥感图像变化检测方法。该方法可快速检测出两时相遥感图像的变化信息，可应用于对地震前后地物变化监测、农作物生长状态的动态监测。

背景技术：
变化检测是通过分析同一地区不同时刻的多幅遥感图像，检测出该地区地物随时间发生变化的信息。随着遥感技术和信息技术的发展，多时相遥感图像变化检测已经成为当前遥感图像分析研究的一个重要方向。在多时相遥感图像变化检测方法的研究中，常见的一种检测方法是先分类后比较法，即首先将遥感图像中不同地物分类，然后再针对不同地物分别构造差异图像，最后利用阈值法确定变化类和非变化类。先分类后比较法的优点在于能够得到针对不同地物的变化检测图像，对阈值选择的方法要求不高，但该方法的存在的不足是，存在分类误差累计问题，即图像质量会对分类结果会造成误差，该误差会影响变化检测的精度。因为光学遥感图像受制于天气，而雾霾的影响尤其明显，雾霾往往污染遥感影像的可见光波段，被摄物体表面反射的光线在大气中由于吸收和散射作用而衰减；同时由于散射作用，摄影镜头中会混入大气光，造成遥感图像的对比度以及分辨率降低，从而降低变化检测精度。所以，在变化检测之前有必要消除或减轻雾霾对遥感图像的影响，提高遥感图像的质量和变化图像的精确度。北京市遥感信息研究所和中国科学院自动化研究所在其申请的专利“高分辨率可见光遥感图像变化检测方法”(专利申请号：201210247866,公开号：CN102842044A)中提出了一种对高分辨率可见光遥感图像变化检测的方法。该方法首先获取图像目标，然后去除虚假目标，并且提取目标的变化特征，然后计算目标的初始变化置信度，选择变化置信度高的部分目标作为正样本，选择变化置信度低的部分目标作为负样本，得到变化监测模型，并且利用训练得到的变化监测模型对目标进行变化检测，最后利用目标的纹理特性对变化目标进行提炼，得到最终的变化检测结果。该方法虽然能够通过形状分析有效去除由于局部配准误差引起的虚假目标，但是仍然存在的不足是，该方法的变化检测精确度依赖于可见光遥感图像的高分辨率，由于雾霾天气会导致遥感图像的分辨率下降，而在获取图像目标之前未进行去雾操作，该方法在对受雾霾影响的遥感图像处理时，分割和识别的精确度较低，导致变化检测结果中存在较多漏检信息和虚警信息，降低了遥感图像的变化检测精度，减小了该方法适用的遥感图像范围。西安电子科技大学在其申请的专利“基于treelet特征融合的遥感图像变化检测方法”(专利申请号：2013104153995，公开号：CN103456020A)中提出了一种基于treelet特征融合的遥感图像变化监测方法。该方法的实现步骤为：(1)输入两时相遥感图像X1和X2；(2)计算两幅图像的差值图，对数比值图和均值比图；(3)获取由三幅不同差异图中像素的邻域向量组成的样本矩阵；(4)用treelet方法对样本矩阵进行特征提取，得到特征向量矩阵；(5)采用Kmeans方法将特征向量矩阵聚为两类，根据聚类结果得到变化检测结果。该方法虽然能够解决单一类型差异图检测精度低，适用范围窄的问题，但是仍然存在的不足是，该方法面向像素点进行操作，从而忽略了面向某一特定地物对象的变化结果，仅得到变化和未变化的两类结果，不能够对遥感图像进行分类，无法得到不同地物的变化检测图像。

技术实现要素：
本发明针对上述现有技术存在的不足，提出了一种遥感图像变化检测方法。本发明既可以较好地处理受雾霾影响的遥感图像，又可以分别得到不同类型地物的变化检测图像，具有较广泛的适用范围和较高的检测精度。本发明实现上述目的的思路是：首先构造输入遥感图像的暗原色通道图像，并判断输入的遥感图像是否需要去雾，若需要，则用暗原色先验方法去雾，然后对去雾后的遥感图像使用模糊C均值FCM聚类方法进行分割，并且使用K近邻KNN算法进行分类识别，最后分别对不同地物变化前后的遥感图像进行变化检测。本发明的步骤包括如下：(1)输入变化前后的遥感图像：输入同一地区，不同时刻获取的两幅变化前后的遥感图像；(2)判断输入的遥感图像是否为彩色遥感图像，若是，执行步骤(3)，否则，将输入的遥感图像作为去雾后的遥感图像，执行步骤(6)；(3)构造暗原色通道图像：(3a)选取彩色遥感图像的某一个像素点，从红色R、绿色G、蓝色B三个颜色通道中，选取其中亮度值最小的颜色通道，将该颜色通道的亮度值作为该像素点的灰度值；(3b)重复步骤(3a)，直至处理完彩色遥感图像中的全部像素点，得到所有像素点的灰度值，将所有像素点的灰度值构成一幅最暗灰度图像；(3c)以最暗灰度图像中的某一像素点为中心像素点，选取一个大小为N1×N1个像素的正方形窗口，其中N1的取值范围为{7,9,11,13,15}；(3d)将正方形窗口中全部像素点的灰度值，按照由小到大的顺序排列，组成一个灰度序列，选取灰度序列中最小的灰度值作为滤波值，用该滤波值替代步骤(3c)中的中心像素点的灰度值；(3e)重复步骤(3c)和步骤(3d)，直至处理完最暗灰度图像中的全部像素点，得到一幅暗原色通道图像；(3f)对暗原色通道图像，按照下式，计算暗像素概率：其中，P表示暗原色通道图像的暗像素概率,P的取值范围为[0,1]，表示将等式左边的暗像素概率P定义为等式右边公式的操作，X表示暗原色通道图像中灰度值小于等于灰度阈值D的像素点个数，灰度阈值D的取值范围为{25,26,27…,35}，X的取值范围为{0,1,2…,M×N}，M、N分别表示暗原色通道图像矩阵的行数和列数；(4)判断彩色遥感图像是否有雾，若彩色遥感图像的暗像素概率P小于先验阈值T1，则认为该彩色遥感图像为有雾彩色遥感图像，执行步骤(5)；否则，将彩色遥感图像作为去雾后的遥感图像，执行步骤(6)，其中，T1表示先验阈值，T1的取值范围为{0.75,0.8,0.85,0.9,0.95}；(5)对有雾彩色遥感图像进行去雾：(5a)将有雾彩色遥感图像对应的暗原色通道图像中的所有像素点，按其灰度值由大到小排列，得到像素点序列，从像素点序列中灰度值最大的像素点开始，依次提取占该序列0.1％比例的像素点，记录所提取像素点在暗原色通道图像中的坐标；(5b)从有雾彩色遥感图像中，提取所记录坐标上的像素点，得到像素点集合O；(5c)将像素点集合O中所有像素点的红色R、绿色G、蓝色B三个颜色通道中的亮度值，由大到小进行排列，得到亮度值序列，从该亮度值序列中选取最大的亮度值作为全局大气光；(5d)选取暗原色通道图像中的某一个像素点，按照下式，计算所选取像素点在彩色遥感图像的透射率图中对应的像素点的灰度值：其中，t(x,y)表示所选取像素点在彩色遥感图像的透射率图中对应的像素点的灰度值，t(x,y)取值范围为[0,1]，(x,y)表示所选取像素点在彩色遥感图像的透射率图中对应的像素点的坐标，E1(x,y)表示暗原色通道图像中所选取像素点的灰度值，E1(x,y)的取值范围为{0,1,2…,255}，(x,y)表示所选取像素点在暗原色通道图像中的坐标，x的取值范围为{0,1,2…,M-1}，y的取值范围为{0,1,2…,N-1}，M、N分别表示暗原色通道图像矩阵的行数和列数，A表示全局大气光；(5e)重复步骤(5d)，直至处理完暗原色通道图像中的所有像素点，得到暗原色通道图像中所有像素点在彩色遥感图像的透射率图中对应像素点的灰度值，将彩色遥感图像的透射率图中对应像素点的灰度值构成透射率图；(5f)按照下式，计算彩色遥感图像去雾后的无雾图像，得到去雾后的遥感图像：其中，L表示去雾后的遥感图像，F表示彩色遥感图像，A表示全局大气光，t表示彩色遥感图像的透射率图，H表示一个大小为M×N的全1矩阵，M、N分别表示彩色遥感图像矩阵的行数和列数；(6)利用模糊C均值FCM聚类算法，对去雾后的遥感图像进行分割，得到分割后各类地物的遥感图像；(7)利用K近邻KNN算法，对分割后各类地物的遥感图像进行分类识别，获得分类识别后的遥感图像；(8)利用作差法，对分类识别后的遥感图像进行变化检测，获得所有地物种类的变化检测图像。本发明与现有技术相比具有以下优点：第一，本发明通过计算暗原色通道图像的暗像素概率判断遥感图像是否有雾，克服了现有技术的变化检测精度依赖于可见光遥感图像的高分辨率的缺点，使得本发明不仅能够对高分辨率的无雾遥感图像进行处理，还能够对分辨率较低的有雾遥感图像进行有效处理。第二，本发明在变化检测前对有雾的遥感图像进行去雾，克服了现有技术在对受雾霾影响的遥感图像处理时变化检测结果精确度较低的缺点，使得本发明减轻了遥感图像受雾霾的影响，提高了遥感图像的分辨率，进而提高了分割和分类识别的准确率，减少了变化检测结果的漏检信息和虚检信息。第三，本发明通过对分类识别后不同地物变化前后的遥感图像进行变化检测，克服了现有技术面向像素点进行操作而无法得到不同地物的变化检测图像的缺点，使得本发明能够获得不同地物的变化检测结果图像，获得更全面的变化检测信息。附图说明图1为本发明的流程图；图2为本发明仿真采用的部分遥感图像；图3为本发明仿真的变化检测图像。具体实施方式下面结合附图对本发明做进一步的详细描述。结合附图1，本发明的步骤如下。步骤1，输入变化前后的遥感图像。输入同一地区，不同时刻获取的两幅变化前后的遥感图像。步骤2，判断输入的遥感图像是否为彩色遥感图像，若是，执行步骤3，否则，将输入的遥感图像作为去雾后的遥感图像，执行步骤6。步骤3，构造暗原色通道图像。(3a)选取彩色遥感图像的某一个像素点，从红色R、绿色G、蓝色B三个颜色通道中，选取其中亮度值最小的颜色通道，将该颜色通道的亮度值作为该像素点的灰度值。(3b)重复步骤(3a)，直至处理完彩色遥感图像中的全部像素点，得到所有像素点的灰度值，将所有像素点的灰度值构成一幅最暗灰度图像。(3c)以最暗灰度图像中的某一像素点为中心像素点，选取一个大小为N1×N1个像素的正方形窗口，其中N1的取值范围为{7,9,11,13,15}。(3d)将正方形窗口中全部像素点的灰度值，按照由小到大的顺序排列，组成一个灰度序列，选取灰度序列中最小的灰度值作为滤波值，用该滤波值替代步骤(3c)中的中心像素点的灰度值。(3e)重复步骤(3c)和步骤(3d)，直至处理完最暗灰度图像中的全部像素点，得到一幅暗原色通道图像。(3f)对暗原色通道图像，按照下式，计算暗像素概率：其中，P表示暗原色通道图像的暗像素概率,P的取值范围为[0,1]，表示将等式左边的暗像素概率P定义为等式右边公式的操作，X表示暗原色通道图像中灰度值小于等于灰度阈值D的像素点个数，灰度阈值D的取值范围为{25,26,27…,35}，X的取值范围为{0,1,2…,M×N}，M、N分别表示暗原色通道图像矩阵的行数和列数。步骤4，判断彩色遥感图像是否有雾。若彩色遥感图像的暗像素概率P小于先验阈值T1，则认为该彩色遥感图像为有雾彩色遥感图像，执行步骤5；否则，将彩色遥感图像作为去雾后的遥感图像，执行步骤6，其中，T1表示先验阈值，T1的取值范围为{0.75,0.8,0.85,0.9,0.95}。步骤5，对有雾彩色遥感图像进行去雾。(5a)将有雾彩色遥感图像对应的暗原色通道图像中的所有像素点，按其灰度值由大到小排列，得到像素点序列，从像素点序列中灰度值最大的像素点开始，依次提取占该序列0.1％比例的像素点，记录所提取像素点在暗原色通道图像中的坐标。(5b)从有雾彩色遥感图像中，提取所记录坐标上的像素点，得到像素点集合O。(5c)将像素点集合O中所有像素点的红色R、绿色G、蓝色B三个颜色通道中的亮度值，由大到小进行排列，得到亮度值序列，从该亮度值序列中选取最大的亮度值作为全局大气光。(5d)选取暗原色通道图像中的某一个像素点，按照下式，计算所选取像素点在彩色遥感图像的透射率图中对应的像素点的灰度值：其中，t(x,y)表示所选取像素点在彩色遥感图像的透射率图中对应的像素点的灰度值，t(x,y)取值范围为[0,1]，(x,y)表示所选取像素点在彩色遥感图像的透射率图中对应的像素点的坐标，E1(x,y)表示暗原色通道图像中所选取像素点的灰度值，E1(x,y)的取值范围为{0,1,2…,255}，(x,y)表示所选取像素点在暗原色通道图像中的坐标，x的取值范围为{0,1,2…,M-1}，y的取值范围为{0,1,2…,N-1}，M、N分别表示暗原色通道图像矩阵的行数和列数，A表示全局大气光。(5e)重复步骤(5d)，直至处理完暗原色通道图像中的所有像素点，得到暗原色通道图像中所有像素点在彩色遥感图像的透射率图中对应像素点的灰度值，将彩色遥感图像的透射率图中对应像素点的灰度值构成透射率图。(5f)按照下式，计算彩色遥感图像去雾后的无雾图像，得到去雾后的遥感图像：其中，L表示去雾后的遥感图像，F表示彩色遥感图像，A表示全局大气光，t表示彩色遥感图像的透射率图，H表示一个大小为M×N的全1矩阵，M、N分别表示彩色遥感图像矩阵的行数和列数。步骤6，利用模糊C均值FCM聚类算法，对去雾后的遥感图像进行分割，得到分割后各类地物的遥感图像。将去雾后的遥感图像中的每一个像素点看做一个样本，并对去雾后的遥感图像提取原始像素灰度矩阵、高斯平滑后的像素灰度值矩阵、边缘像素灰度值矩阵三种特征矩阵，组成特征值矩阵，利用特征值矩阵对原图像进行FCM分割，聚类数设为地物类型数。首先定义几个公式如下：聚类中心：隶属度：目标函数：其中，d2(xi,vj)表示相似度测量函数，其表达式为：d2(xi,vj)＝|xi-vj|其中，vj表示聚类中心，N表示总样本数，图像中也认为是像素点的个数，uji表示第i个像素点对第j类的隶属度，m表示控制模糊度的权重指数，一般设为2，xi表示样本点，c表示聚类数目，设为地物数目，dji表示第i个像素点到第j个聚类中心的欧氏距离，dki表示第i个像素点到第k个聚类中心的欧氏距离，d2(xi,vj)表示相似度测量函数。具体步骤为：(6a)初始化聚类中心ci,其中i＝1,…,c，典型的做法是从所有数据点中任取c个点。(6b)利用上述公式(2)确定隶属度矩阵U。(6c)利用上述公式(3)计算目标函数，如果它相对上次目标函数值的改变量小于目标函数最小变化量阈值ε，则算法停止。(6d)利用上述公式(1)修正聚类中心，返回步骤(6b)。算法停止后，得到聚类中心、隶属度矩阵和目标函数，根据最大隶属度原则将每个像素点分别归为所有地物种类中的某一类。步骤7，利用K近邻KNN算法，对分割后各类地物的遥感图像进行分类识别，获得分类识别后的遥感图像。(7a)分别从每种选定的已知地物类别信息的遥感图像中，随意截取20个大小为10×10的子图像矩阵，建立不同地物类别遥感图像的训练集。(7b)分别从分割后各类地物的遥感图像中，随意截取50个大小为10×10的子图像矩阵，建立不同地物的测试集。(7c)分别对不同地物类别的遥感图像的训练集中的每一个训练样本，使用4个尺度5个方向的伽柏滤波器Gaborfilter，提取20维的特征列向量，得到训练集的特征向量集：其中，Sw表示第w种地物的训练集的特征向量集，w表示训练集的特征向量集所属的地物类别标号，w的取值范围为{1,2…，Z}，Z表示遥感图像中的地物种类数，si表示第w种地物训练集中的第i个训练样本的特征列向量。(7d)分别对不同地物的测试集中的每一个测试样本，使用4个尺度5个方向的伽柏滤波器Gaborfilter，提取20维的特征列向量，得到测试集的特征向量集：其中，Ur表示第r种地物测试集的特征向量集，r表示测试集的特征向量集的标号，r的取值范围为{1,2,…,Z}，Z表示遥感图像中的地物种类数，uj表示第r种地物测试集中的第j个测试样本的特征列向量。(7e)选取第r种地物测试集的特征向量集Ur，按照下式，计算所选取的特征向量集Ur到第w种地物的训练集的特征向量集Sw的欧氏距离：其中，dr,w表示第r种地物测试集的特征向量集Ur到第w种地物的训练集的特征向量集Sw的欧氏距离，r表示测试集的特征向量集的标号，r的取值范围为{1,2,…,Z}，w表示训练集的特征向量集所属的地物类别标号，w的取值范围为{1,2…，Z}，Z表示遥感图像中的地物种类数，uj表示第r种地物测试集中的第j个测试样本的特征列向量，si表示第w种地物训练集中的第i个训练样本的特征列向量。(7f)重复第五步，直至计算完第r种地物测试集的特征向量集Ur到所有地物的训练集的特征向量集Sw的欧氏距离，得到第r种地物的欧氏距离集合{dr,1,dr,2…,dr,Z}。(7g)将第r种地物的欧氏距离集合中的所有元素由小到大进行排列，得到该集合中元素的序列，从该集合中元素的序列中选取一个最小的元素，将该最小元素的第2个下标，作为第r种地物的遥感图像所属的地物种类标号。(7h)重复第五步、第六步、第七步，直至处理完所有地物测试集的特征向量集，得到分割后各类地物的遥感图像所属的地物种类标号。步骤8，利用作差法，对分类识别后的遥感图像进行变化检测，获得所有地物种类的变化检测图像。(8a)对分类识别后的遥感图像，选取某一类地物在同一地区，不同时刻获取的变化前后的两幅遥感图像，将这两幅遥感图像矩阵相减后取绝对值，得到该种地物的差值图像。(8b)设定变化检测阈值Y，其中Y的取值范围为{80,90,100,110,120}，将该种地物的差值图像中灰度值小于变化检测阈值Y的像素点的灰度值置0，得到该种地物的变化检测图像。(8c)重复步骤(8a)和步骤(8b)，直至处理完所有地物种类的遥感图像，得到所有地物种类遥感图像的变化检测图像。下面结合仿真图对本发明的效果做进一步说明。1.仿真条件：本发明的仿真是在主频2.4GHZ的Intel双核CPUB940、内存4GB的硬件环境和MATLABR2013a的软件环境下进行的。本发明实施例中，构造暗原色通道图像步骤中采用了15×15个像素的正方形窗口，并且灰度阈值D＝35，判断遥感图像是否有雾步骤中选取先验阈值T1＝0.8，分割步骤中设置聚类数c＝3，控制模糊度的权重指数m＝2，初始化聚类中心ci时任意选取c个像素点，目标函数最小变化量阈值ε＝1e-5，变化检测步骤中选取的变化检测阈值Y＝100。2.仿真内容：本发明仿真实验所用数据为一组真实遥感数据集，数据来自于互联网上对汶川地震报道中的新闻图片，图像的大小均为300×238像素，它们之间发生的变化是由地震破坏了大面积的植被和山上的泥土和石块阻塞河道所致。3.仿真效果分析：本仿真首先对汶川地震新闻报道中的12幅地震前后的遥感图像进行实验，首先计算出每一幅输入遥感图像对应的暗像素概率P，然后根据先验阈值T1对输入遥感图像判断是否有雾。判断结果如表1所示，表中Y表示有雾，N表示无雾。表1.对12幅遥感图像判断是否有雾的结果在实验数据中，无雾图像的暗像素概率P均大于0.8，而无雾图像的暗像素概率P均小于0.8，12幅遥感图像均得到了正确的判断结果。根据汶川地震图像前后的特点及地震自然灾害的主要类型，着重对遥感图像中的三种地物，即草地、裸地和水域进行识别。去雾前后的识别正确率如表2所示。表2.去雾前后的识别正确率图像类型草地识别率(％)裸地识别率(％)水域识别率(％)总识别率(％)去雾前图像93758358去雾后图像100839275经过去雾后，各种地物的识别率和总识别率均有不同程度地提升。实验结果说明在遥感图像变化检测前进行去雾，能够提高图像识别的准确度和分类的精确度。图2为本发明仿真采用的部分遥感图像，其中图2(a)和图2(b)分别为2008年5月12日地震之前和2008年5月13日地震之后的的北川遥感图像。图3为本发明仿真的变化检测图像，其中图3(a)为草地的变化检测图像，图3(b)为水域的变化检测图像，图3(c)为裸地的变化检测图像。由图3(b)和图3(c)可以看出，水域有部分河道变化为裸地，所以有出现河道堵塞和河水泛滥的可能，应及时转移周边群众到安全地带。由图3(a)和图3(b)可以看出，裸地面积增加，草地面积减少，出现泥石流的可能性增大，应警告民众及车辆绕道行走。