基于阈值优化的SAR图像变化检测方法与流程

本发明属于图像处理
技术领域：
，更进一步涉及图像分割与图像分类识别
技术领域：
中的一种基于阈值优化的合成孔径雷达sar(syntheticapertureradar)图像变化检测方法。该方法可快速检测出两时相合成孔径雷达图像的变化信息，可应用于对地震前后地物变化监测、农作物生长状态的动态监测。
背景技术：
随着遥感技术和信息技术的发展，基于合成孔径雷达sar图像的变化检测由于传感器具有不受时段、天气条件影响等优良特性而在近年内受到了广泛的关注，而阈值技术是图像变化检测中的关键技术之一，阈值技术在合成孔径雷达sar图像变化检测领域中至关重要。目前关于阈值技术的研究多需要事先建立模型对未变化像素和变化像素进行假设分布，并且只能获得非整数的阈值解，从而降低了变化检测的准确率。西安电子科技大学在其申请的专利文献“用于变化检测的高斯对数模型单边曲率拟合阈值方法”(专利申请号：201010548359.4，申请公布号：cn102005050a)中提出了一种用于合成孔径雷达sar图像变化检测的高斯对数模型单边曲率拟合阈值方法。该方法的实现过程为：首先对两幅合成孔径雷达sar图像构造差异图并求出差异图像的直方图，在直方图上确定单边拟合区域，应用高斯对数模型对单边直方图区域的直方图曲线进行曲率拟合，求出单边拟合区域中曲率拟合误差最小的阈值，确定无变化区域的直方图概率分布函数和初始阈值，根据初始阈值确定基于高斯模型变化区域直方图概率分布函数，最后用最大后验概率方法确定最终阈值，通过该阈值构造合成孔径雷达sar图像变化检测结果图。该方法虽然能够较好拟合差异图直方图分布，降低变化检测的错误率，但是，该方法仍然存在的不足之处是，需要使用高斯对数模型对单边直方图区域的直方图曲线进行曲率拟合，并且高斯对数模型不能完全拟合差异图直方图的直方图曲线，影响了合成孔径雷达图像变化检测精度。浙江环球星云遥感科技有限公司在其申请的专利文献“一种基于变化因子的多时相sar影像变化检测方法”(专利申请号：201710804479.8，申请公布号：cn107689051a)中提出了一种基于变化因子的多时相sar影像变化检测方法。该方法的实现过程为：(1)对多个时相的合成孔径雷达数据进行辐射校正，生成多时相后向散射系数影像；(2)对多时相后向散射系数影像进行滤波处理，生成多时相滤波影像；(3)对多时相滤波影像进行时序差分运算得到地物发生变化的时间；(4)根据变化因子计算方法生成变化差异图像；(5)对变化差异图像运用阈值分割方法得到地物变化信息。该方法虽然能够降低自然地物季节变化对图像造成的影响，但是，该方法仍然存在的不足之处是，在用阈值分割方法处理差异图时，由于所使用的阈值是整数，导致用整数阈值分割差异图中变化区域与未变化区域的划分界限不够准确的缺点。技术实现要素：本发明针对上述现有技术存在的不足，提出了一种基于阈值优化的合成孔径雷达sar图像变化检测方法。实现本发明目的的具体思路是，首先预处理输入的合成孔径雷达图像，构造输入图像的差异图，接着产生初始阈值，并根据阈值生成变化检测图，利用kappa系数计算出阈值的适应度，然后利用径向基函数计算每个候选阈值的替代估计值，利用每个候选阈值的替代估计值和距离度量值计算候选阈值的采样权重值，根据采样权重产生当前更新阈值并得到当前更新阈值的适应度，最后根于阈值适应度大小，产生最优阈值，得到最终的变化检测图。该方法的具体步骤包括如下：本发明的具体实现步骤如下：(1)输入合成孔径雷达sar图像和变化参考图：输入同一地区，不同时刻获取的两幅合成孔径雷达图像以及该地区的变化参考图；(2)预处理输入图像：(2a)判断输入的合成孔径雷达图像和变化参考图是否为彩色图像，若是，执行步骤(2b)，否则，执行步骤(2c)(2b)将输入的合成孔径雷达图像和变化参考图转化成灰度图像；(2c)选取3×3中值滤波器分别对变化前、变化后的合成孔径雷达图像进行中值滤波，得到中值滤波后的两幅变换前、变化后的合成孔径雷达图像；(3)按照下式，构造差异图：其中，i表示构造的差异图，log表示以2为底的对数操作，i1表示中值滤波后的变化前的图像，i2表示中值滤波后的变化后图像；(4)产生初始阈值：(4a)在区间[min，max]产生一个随机数，将该随机数作为初始阈值，其中，min表示差异图中所有像素的最小值，max表示差异图中所有像素的最大值；(4b)将初始阈值加入阈值集合中；(5)生成变化检测图：(5a)从差异图中任意选取的一个像素点；(5b)判断所选取像素点的像素值是否小于等于初始阈值，若是，将所选像素点的像素值设置为0，否则，将所选像素点的像素值设置为255；(5c)判断是否选取完差异图中的所有像素点，若是，执行步骤(6)，否则，执行步骤(5a)；(6)计算初始阈值的适应度；(6a)利用系数kappa系数计算法，计算变化检测图与变化参考图之间的kappa系数；(6b)用1减去kappa系数，将差值作为初始阈值的适应度值；(7)产生候选阈值：在[0,1]区间随机产生100个不同的随机数，利用候选解采样法，分别生成与所选100个不同的随机数对应的100个候选阈值；(8)计算每个候选阈值的替代估计值：(8a)按照下式，计算每个候选阈值的径向基插值：其中，sk(zk)表示第k个候选阈值zk的径向基插值，n表示阈值集合中阈值的个数，∑表示求和操作，i表示阈值xi在阈值集合中的序号，m表示径向基的序号，λm表示由线性方程组得到的第m个径向基的系数，||·||表示求欧式距离操作，zk表示第k个候选阈值，xi表示阈值集合中的第i个阈值，b表示由线性方程组得到的权重值，a表示由线性方程组得到的偏移量值；(8b)将100个候选阈值的径向基插值从大到小排序，得出径向基插值的最大值和最小值；(8c)利用径向基函数采用替代估计值计算公式计算每个候选阈值的替代估计值；(9)计算每个候选阈值的距离度量值：(9a)采用欧氏距离公式，分别计算每个候选阈值与阈值集合每个阈值的欧氏距离，得到每个候选阈值的欧氏距离集合，将每个候选阈值的欧氏距离集合中的最小值作为该候选阈值的最小欧氏距离；(9b)将100个候选阈值的最小欧氏距离从大到小排序，得出最大值和最小值；(9c)利用距离度量值公式，计算每个候选解的距离度量值；(10)产生当前更新阈值：(10a)利用每个候选阈值的替代估计值和距离度量值，采用采样权重公式，计算每个候选阈值的采样权重值；(10b)对100个候选阈值的采样权重从大到小排序，将最大的采样权重对应的候选阈值作为当前更新阈值；(10c)将当前更新阈值加入阈值集合中；(11)生成变化检测图：(11a)从差异图中任意选取的一个像素点；(11b)判断所选取像素点的像素值是否小于等于当前更新阈值，若是，将所选像素点的像素值设置为0，否则，将所选像素点的像素值设置为255；(11c)判断是否选取完差异图中的所有像素点，若是，执行步骤(12)，否则，执行步骤(11a)；(12)计算当前更新阈值的适应度；(12a)利用系数kappa系数计算法，计算变化检测图与变化参考图之间的kappa系数；(12b)用1减去kappa系数，将差值作为当前更新阈值的适应度值；(13)判断当前更新阈值与前一次阈值的差的绝对值是否小于终止精度，若是，执行步骤(14)，否则，执行步骤(7)；(14)产生最优阈值：对阈值集合中所有阈值的适应度从大到小排序，从中选取最小适应度的阈值，将该阈值作为最优阈值；(15)生成变化检测图：(15a)从差异图中任意选取的一个像素点；(15b)判断所选取像素点的像素值是否小于等于最优阈值，若是，将所选像素点的像素值设置为0，否则，将所选像素点的像素值设置为255；(15c)判断是否选取完差异图中的所有像素点，若是，执行步骤(16)，否则，执行步骤(15a)；(16)得到变化检测图。本发明与现有技术相比具有以下优点：第一，由于本发明利用径向基函数采用替代估计值计算公式计算每个候选阈值的替代估计值，用替代估计值拟合阈值，克服了现有技术使用高斯对数模型，对单边直方图区域的直方图曲线进行曲率拟合的问题，并且高斯对数模型不能完全拟合差异图直方图的直方图曲线，从而影响了合成孔径雷达图像的变化检测精度的问题，使得本发明不需要事先建立模型对差异图直方图中变化像素和未变化像素进行假设分布，提高了变化检测的检测正确率。第二，由于本发明利用每个候选阈值的替代估计值和距离度量值，采用采样权重公式，计算每个候选阈值的采样权重值，从候选阈值集合中选取采样权重值最大的候选阈值作为当前更新阈值，克服了现有技术在用阈值分割方法处理差异图时，由于所使用的阈值是整数，导致了用整数阈值分割差异图中变化区域与未变化区域的划分界限不够准确的问题，使得本发明获得的非整数阈值具有更高的分割精度的优点。附图说明图1为本发明的流程图；图2为本发明仿真实验使用的渥太华地区的合成孔径雷达图像；图3为本发明方法和三种对比方法对渥太华地区的仿真结果图。具体实施方式下面结合附图对本发明做进一步的详细描述。结合附图1，对本发明的步骤做进一步的详细描述。步骤1，输入合成孔径雷达sar图像和变化参考图。输入同一地区，不同时刻获取的两幅合成孔径雷达图像以及该地区的变化参考图。步骤2，预处理输入图像。a.判断输入的合成孔径雷达图像是否为彩色图像，若是，则执行本步骤的a，否则，执行本步骤的c。b.将输入的两幅合成孔径雷达图像和变化参考图转化成灰度图像。c.选取3×3中值滤波器分别对变化前、变化后的合成孔径雷达图像进行中值滤波，得到中值滤波后的两幅变换前、变化后的合成孔径雷达图像。步骤b中所述的将合成孔径雷达图像和变化参考图转化成灰度图像的步骤如下：第1步，从合成孔径雷达图像中任意选取的一个像素点。第2步，按照下式，计算选取像素点的灰度值。i灰＝r×0.3+g×0.59+b×0.11其中，i灰表示所选取像素点的灰度值，r表示所选取像素点红色通道的亮度，g表示所选取像素点绿色通道的亮度，b表示所选取像素点蓝色通道的亮度。第3步，对两幅合成孔径雷达图像和变化参考图重复本步骤的第1步和第2步，直至处理完图像中的全部像素点，得出变化前合成孔径雷达图像的灰度图像、变化后合成孔径雷达图像的灰度图像和变化参考图的灰度图像。步骤3，按照下式，构造差异图：其中，i表示构造的差异图，log表示以2为底的对数操作，i1表示中值滤波后的变换前的图像，i2表示中值滤波后的变换后图像。步骤4，产生初始阈值。在区间[min，max]内产生一个随机数，将该随机数作为初始阈值，其中，min表示差异图中所有像素的最小值，max表示差异图中所有像素的最大值。将初始阈值加入阈值集合中。步骤5，生成变化检测图。a.从差异图中任意选取的一个像素点。b.判断所选取像素点的像素值是否小于等于初始阈值，若是，将所选像素点的像素值设置为0，否则，将所选像素点的像素值设置为255。c.判断是否选取完差异图中的所有像素点，若是，执行步骤6，否则，执行本步骤的a。步骤6，计算初始阈值的适应度。利用kappa系数计算法，计算变化检测图与变化参考图之间的kappa系数。所述的kappa系数计算法的具体步骤如下：第1步，将变化检测图与变化参考图做差值处理，得出差值图，差值图中像素点值为0，255或-255。第2步，按照下式，计算变化检测图与变化参考图之间的kappa系数：其中，k表示变化检测图与变化参考图之间的kappa系数，tp表示由变化检测图中像素值为0与变化参考图中像素值为0做差得到的差值图中所有像素点数目，tn表示由变化检测图中像素值为255与变化参考图中像素值为255做差得到的差值图中所有像素点数目，fn表示差值图中像素值为-255的所有像素点数目，fn表示差值图中像素值为255的所有像素点数目。用1减去kappa系数，将差值作为初始阈值的适应度值。步骤7，产生候选阈值。在[0,1]区间随机产生100个不同的随机数，利用候选解采样法，分别生成与所选100个不同的随机数对应的100个候选阈值。所述候选解采样法的步骤如下：第1步，对阈值集合中所有阈值的适应度从大到小排序，从中选取最小适应度的阈值。第2步，用最大适应度的阈值加上一个[0,1]区间的随机数，将其和作为临时阈值。第3步，将差异图中最大像素值与临时阈值二者中的最小值作为中间阈值。第4步，将差异图中最小像素值与中间阈值二者中的最大值作为候选阈值。步骤8，计算每个候选阈值的替代估计值。按照下式，计算每个候选阈值的径向基插值：其中，sk(zk)表示第k个候选阈值zk的径向基插值，n表示阈值集合中阈值的个数，∑表示求和操作，i表示阈值xi在阈值集合中的序号，m表示径向基的序号，λm表示由线性方程组得到的第m个径向基的系数，||||表示求欧式距离，zk表示第k个候选阈值，xi表示阈值集合中的第i个阈值，b表示由线性方程组得到的权重值，a表示由线性方程组得到的偏移量。所述的线性方程组公式为：其中，φij＝(||xi-xj||)3表示三样条径向基函数，||·||表示欧几里得二范式，xi表示阈值集合中的第i个阈值，xj表示阈值集合中的第j个阈值，f(xi)表示阈值xi的适应度。将100个候选阈值的径向基插值从大到小排序，得出径向基插值的最大值和最小值。利用替代估计值计算公式，计算每个候选阈值的替代估计值。所述替代估计值计算公式如下：其中，vk(zk)表示第k个候选阈值zk的替代估计值，sk(zk)表示第k个候选阈值zk的径向基插值，zk表示第k个候选阈值，smax表示径向基插值的最大值，smin表示径向基插值的最小值。步骤9，计算每个候选阈值的距离度量值。采用欧氏距离公式，分别计算每个候选阈值与阈值集合每个阈值的欧氏距离，将每个候选阈值与阈值集合中所有阈值的欧氏距离组成集合，将每个候选阈值的欧氏距离集合中的最小值作为该候选阈值的最小欧氏距离。将100个候选阈值的最小欧氏距离从大到小排序，得出最大值和最小值。利用距离度量值的计算公式，计算每个候选解的距离度量值。所述距离度量值的计算公式如下：其中，dk(zk)表示第k个候选阈值zk的距离度量值，δmax表示100个最小欧氏距离中的最大值，dmin(zk)表示第k个候选阈值zk的欧氏距离集合中的欧氏距离的最小值，zk表示第k个候选阈值，δmin表示100个最小欧氏距离中的最小值。步骤10，产生当前更新阈值。按照下式，计算每个候选阈值的采样权重：w(zk)＝wvk(zk)+(1-w)dk(zk)其中，w(zk)表示第k个候选阈值zk的采样权重，w＝0.3表示平衡权重，zk表示第k个候选阈值，表示第k个候选阈值zk的替代估计值，dk(zk)表示第k个候选阈值zk的距离度量值。对100个候选阈值的采样权重从大到小排序，将最大的采样权重的对应的候选阈值作为当前更新阈值。将当前更新阈值加入阈值集合中。步骤11，生成变化检测图.a.从差异图中任意选取的一个像素点。b.判断所选取像素点的像素值是否小于等于当前更新阈值，若是，将所选像素点的像素值设置为0，否则，将所选像素点的像素值设置为255。c.判断是否选取完差异图中的所有像素点，若是，执行步骤12，否则，执行本步骤的a。步骤12，计算当前更新阈值的适应度。利用kappa系数计算法，计算变化检测图与变化参考图之间的kappa系数；用1减去kappa系数，将差值作为当前更新阈值的适应度值；步骤13，判断当前更新阈值与前一次阈值的差的绝对值是否小于终止精度10e-4，若是，则执行步骤14，否则，执行步骤7。步骤14，产生最优阈值。对阈值集合中所有阈值的适应度从大到小排序，从中选取最小适应度的阈值，将该阈值作为最优阈值。步骤15，生成变化检测图。a.从差异图中任意选取的一个像素点。b.判断所选取像素点的像素值是否小于等于最优阈值，若是，将所选像素点的像素值设置为0，否则，将所选像素点的像素值设置为255。c.判断是否选取完差异图中的所有像素点，若是，执行步骤16，否则，执行本步骤的a。步骤16，得到变化检测图。下面结合仿真实验对本发明的效果做进一步的说明。1.仿真条件：本发明仿真实验的运行环境是：windows10，64位操作系统，本发明的仿真实验是在主频2.60ghz的intelcore(tm)i7-4720cpu，内存8gb的硬件环境和matlabr2014b软件环境下完成的。2.仿真内容：本发明仿真实验是采用本发明和三种现有技术(nr-gkit阈值法、ln-gkit阈值法和wr-gkit阈值法)，分别对渥太华地区的合成孔径雷达图像进行仿真，得到合成孔径雷达图像变化检测的最优阈值，并根据最优阈值分割得到变化检测图。本发明仿真实验所用数据为一组真实的合成孔径雷达数据集，图2为加拿大国防研究与发展部提供的由radarsatsar传感器获得渥太华地区的数据集(290×350像素)，该图像的拍摄时间分别为1997年5月和1997年8月。图2(a)为渥太华地区变化前图像，图2(b)为渥太华地区变化后图像，图2(c)为渥太华地区真实的变化参考图。图3(a)为采用现有技术nr-gkit阈值法对渥太华地区合成孔径雷达图像进行仿真的结果图。图3(b)为采用现有技术ln-gkit阈值法对渥太华地区合成孔径雷达图像进行仿真的结果图。图3(c)为采用现有技术wr-gkit阈值法对渥太华地区合成孔径雷达图像进行仿真的结果图。图3(d)为采用本发明的阈值方法对渥太华地区合成孔径雷达图像进行仿真的结果图。图3(a)中的白色区域为渥太华地区合成孔径雷达图像中水域变化为裸地的区域。比较图3(a)和图2(c)可以看出，图3(a)中的白色区域减少，且白色区域轮廓边缘不清晰，说明现有技术nr-gkit阈值法不能准确检测出渥太华地区合成孔径雷达图像的变化信息。比较3(b)和图2(c)可以看出，图3(b)中杂斑较多，白色区域减少，说明现有技术ln-gkit阈值法对渥太华地区合成孔径雷达图像中水域变为裸地的信息没有检测出来。比较图3(d)和图2(c)可以看出，图3(d)中白色区域接近于真实的变化参考图2(c)中白色区域，说明采用本发明的阈值方法对渥太华地区合成孔径雷达图像中发生变化的区域的检测较为准确。3.仿真效果分析：在本发明与三种现有技术(nr-gkit阈值法、ln-gkit阈值法和wr-gkit阈值法)分别对渥太华地区的合成孔径雷达图像进行仿真的过程中，需要分别计算用于渥太华地区合成孔径雷达图像分割差异图的最优阈值以及渥太华地区合成孔径雷达图像仿真的结果图与真实变化参考图之间kappa系数值。kappa系数越大，仿真实验的结果图与真实变化参考图之间越接近，变化检结果能更好。本发明仿真实验的四种方法所计算的最优阈值和kappa系数值的比较结果相见表1。其中，nr-gkit表示现有技术的nr-gkit阈值法的计算结果，ln-gkit表示现有技术的ln-gkit阈值法的计算结果，wr-gkit表示现有技术的wr-gkit阈值法的计算结果。表1四种方法的计算结果比较表算法最优阈值kappa(％)nr-gkit1.088.46wr-gkit2.028.23ln-gkit1.088.46本发明0.844989.66从表1中可以看出，处理渥太华地区合成孔径雷达图像时，本发明达到了89.66％的kappa系数值，分别比三种现有技术(nr-gkit阈值法、ln-gkit阈值法和wr-gkit阈值法)的kappa系数值提高了1.20％、61.43％和1.20％。而且，三种现有技术方法得到的阈值只能为整数，忽略了最优阈值为小数的情况，本发明能够获得小数的最优阈值，提高了阈值分割差异图时的阈值分割精度。综上所述，相较于三种现有技术(nr-gkit阈值法、ln-gkit阈值法和wr-gkit阈值法)，本发明能够得到kappa系数更高的变化检测结果，并且能够获得非整数的最优阈值。当前第1页12