一种基于变化因子的多时相SAR影像变化检测方法与流程

本发明涉及卫星遥感应用领域，尤指一种多时相合成孔径雷达影像变化检测方法。

背景技术：

遥感影像变化检测是以一个时相的影像为参考，检测出其他时相影像相对于参考影像的差异，广泛应用于国民经济和国防建设的诸多领域，例如，资源和环境监测中的土地利用和土地覆盖变化、森林和植被变化、湿地变化、城区建设变化等信息获取，农作物生长监测及自然灾害中地震、洪水、泥石流和森林砍伐、火灾等灾情监测与评估等。

与可见光和红外遥感相比，微波遥感成像不受时间(昼夜)及成像条件(云雾)等的影响，具有全天候、全天时的工作能力，为变化检测提供了重要的数据源。sar(syntheticapertureradar)影像变化检测是指利用覆盖同一地区的多时相sar影像发现和分析地表变化，提供地物空间分布的定性与定量信息，是遥感影像变化检测研究中的一个重要方向.随着radarsat-2、terrasar-x、cosmosky-med、sentinel-1等大幅宽、高分辨率、高重访、多极化卫星的发展，基于sar影像的变化检测被广泛应用于土地利用类型调查、变化检测等应用，特别是镇区域的土地利用变化检测。因此，开展sar影像变化检测方法的研究具有重要的现实意义和应用价值。

遥感影像的变化检测技术中，需要利用阈值分割手段提取变化信息.如何确定变化区域与非变化区域之间的最优阈值被公认为是变化检测中的挑战性问题，目前常用的阈值选择方法是通过人机交互方式选择数值，这要求操作人员有比较丰富的经验，并且工作量比较大，通常难以得到较好的检测结果。

现有的sar影像变化检测方法对于永久变化和随机变化的区分能力有限。基于两个时期的sar影像能够检测出变化区域，然而这些检测出的变化区域中含有大量的虚假变化信息，例如农田、植被的生长等季节性变化，这些随机变化会很大程度的升高虚警率，降低了变化检测精度。因此需要发展有效的多时相sar影响变换检测方法，有效区分地物散射单元在时间序列上的永久变化与随机变化，抑制随机变化而放大永久型变化，提高变化检测的准确性与稳定性。

技术实现要素：

针对现有技术存在的问题，本发明提供一种多时相sar影像变化检测方法，该方法通过对多时期sar影像数据辐射定标和滤波处理，基于时域变化因子构建地物变化差异图像，经过阈值分割后得到地物变化信息。

为实现上述目标，本发明的一种基于变化因子的多时相sar影像变化检测方法，具体步骤为：1)对多个时相的sar数据进行辐射校正，生成多时相后向散射系数影像；2)对多时相后向散射系数影像进行sar影像滤波处理，生成多时相滤波影像；3)对多时相滤波影像进行时序差分运算得到地物发生变化的时间；4)根据变化因子计算方法由多时相滤波影像生成地物变化差异图像；4)对地物变化差异图像运用阈值分割方法得到地物变化信息。

进一步，步骤1)中所述的辐射校正是一种将sar影像的像素值转变为后向散射系数值的处理方法，sar数据辐射校正是实现从sar图像dn值到地面目标的后向散射系数的转化，标定sar系统测量目标后向散射信号幅度和相位的能力。sar图像绝对定标的公式：

σ0＝k(dn²·sinθ-n)，

其中dn为sar图像的像素值，θ为像元入射角，n是像元噪声等效值；k为定标常数，与传感器全系统增益、天线辐射方向图、以及系统成像传递函数等因子有关，是通过sar系统外定标予以确认的待定常数。

进一步，步骤2)中所述的sar影像滤波是一种sar影像斑点噪声抑制方法，sar成像系统是基于相干成像原理的，在雷达回波信号中，相邻像素点的灰度会由于相干性而产生一些随机的变化，并且这种随机变化是围绕着某一均值而进行的，从而会影响从图像中提取后向散射系数值的精度.所述滤波采用增强lee滤波，滤波公式为：

其中，mij为滤波后的sar图像中第i行第j列的像元值，aij为滤波前的sar图像中第i行第j列的像元值，为滤波前的sar图像中以第i行第j列为中心的7×7窗口内像元值的均值，β＝1-cu²/cl²，cu＝0.5227，σij为滤波前的sar图像中以第i第j列为中心的7×7窗口内像元局部标准差.

进一步，步骤3)中所述的时序差分运算过程为：首先，将滤波后的sar影像按照成像时间排序形成时序影像s＝{mp|(p＝1，2，...，n)}，n为序列影像的个数；然后，对于任意小于n的正整数k，将时序影像s＝{mp|(p＝1，2，...，n)}分成前序列s¹k＝{mp|(p＝1，2，...，k)}和后序列s²k＝{mp|(p＝k+1，...，n)}，后序列s²k的均值减去前序列s¹k的均值为dk，对于k＝1，2，...，n-1，可以得到时序差d＝{dk|(k＝1，2，...，n-1)}；最后，找到时序差d中的差异最大值dmax，其对应的序号即为地物发生变化时间t。

进一步，步骤4)中变化因子计算方法为：首先，利用地物发生变化时间t将时序影像s＝{mp|(p＝1，2，...，n)}分成变化前序列影像s¹t＝{mp|(p＝1，2，...，t)}和变化后序列影像s²t＝{mp|(p＝t+1，...，n)}；然后，将变化后序列影像s²t＝{mp|(p＝t+1，...，n)}减去差异最大值dmax得到新的变化后序列影像将变化前序列影像s¹t＝{mp|(p＝1，2，...，t)}和新的变换后序列影像组合成归一化序列影像最后，计算归一化序列影像的标准差std，利用计算公式cf＝dmax/std得到变化差异图像icf。

进一步，变化差异图像icf是l×c的矩阵，l是变化差异图像的高度，c为变化差异图像的宽度，矩阵中每个元素为对应的变化因子值。

进一步，步骤5)中地物变化信息提取包含增强变化提取和减弱变化提取。增强变化是在地物变化过程中，后向散射系数值增大；减弱变化是在地物变化过程中，后向散射系数值减小.

进一步，增强变化提取的过程为：对变化差异图像icf进行阈值分割，大于增强变化阈值th1的像素为增强变化，生成增强变化掩模图像r⁺。

进一步，增强变化阈值th1的计算方法为th1＝m1+δ1，m1为变化差异图像icf中大于0的像素值的均值，δ1为变化差异图像icf中大于0的像素值的标准差；

进一步，增强变化掩模图像r⁺的生成规则为：sij≥th1，r⁺ij＝+1；sij＜th1，r⁺ij＝0。sij为变化差异图像icf中第i行第j列的像素值，r⁺ij为增强变化掩模图像中第i行第j列的像素值，值为1则表示该像素发生了增强变化，值为0表示该像素未发生增强变化，i＝1，2，…，j＝1，2，....

进一步，减弱变化提取的过程为：对变化差异图像icf进行阈值分割，小于减弱变化阈值th2的像素为减弱变化，生成减弱变化掩模图像r^-.

进一步，减弱变化阈值th2的计算方法为th2＝m2-δ2，m2为变化差异图像icf中小于0的像素值的均值，δ2为变化差异图像icf中小于0的像素值的标准差。

进一步，减弱变化掩模图像r^-的生成规则为：sij≤th2，r^-ij＝+1；sij＞th2，r^-ij＝0。sij为变化差异图像icf中第i行第j列的像素值，r^-ij为减弱变化掩模图像中第i行第j列的像素值，值为1则表示该像素发生了减弱变化，值为0表示该像素未发生减弱变化，i＝1，2，…，j＝1，2，....

本发明针对sar影像变换检测，提出了一种基于变化因子的多时相sar影像变化检测方法，可以有效区分地物散射单元在时间序列上的永久变化与随机变化，抑制自然地物季节变化导致的随机变化放大与人类活动相关的永久型变化，提高变化检测的准确性与稳定性。

附图说明

图1为本发明的一种基于变化因子的多时相sar影像变化检测方法的流程图；

图2为8个时相sar影像数据；

图3为辐射校正后的sar影像；

图4为影像滤波处理后的sar影像；

图5为地物发生变化时间的图像；

图6为变化差异图像；

图7为增强变化的掩模图像；

图8为减弱变化的掩模图像。

具体实施方式

本发明的多时相sar影像变化监测方法是基于多时相sar影像，通过计算变化因子构建地物变化差异图像，经阈值分割后得到地物变化信息。变化检测流程图如图1所示，首先对多时相sar影像进行辐射校正，生成后向散射系数影像，然后进行sar影像滤波处理抑制斑点噪声，根据多时相滤波影像时序差分运算得到地物发生变化的时间，按照地物变化时间将多时相滤波影像分成变化前后序列影像，通过变换因子计算方法得到变化差异图像，最后利用阈值分割方法得到变化信息.

以南京市8个时相sar卫星遥感影像来说明本专利的实施方法。sar数据的具体参数为：radarsat-2卫星，超宽精细模式，vv极化，入射角为28.7°，像元大小为2.7米×2.9米(斜距向×方位向)，成像时间分别为2016年10月01日、2016年10月25日、2016年12月12日、2017年01月29日、2017年02月22日、2017年03月18日、2017年04月11日、2017年05月05日。图2为8个时相sar影像数据。

对本发明的变化检测方法具体介绍如下：

1.sar影像辐射校正方法

sar数据辐射校正是实现从sar图像dn值到地面目标的后向散射系数的转化，标定sar系统测量目标后向散射信号幅度和相位的能力。sar图像绝对定标的公式：

σ0＝k(dn².sinθ-n)，

其中dn为sar图像的像素值，θ为像元入射角，n是像元噪声等效值；k为定标常数，与传感器全系统增益、天线辐射方向图、以及系统成像传递函数等因子有关，是通过sar系统外定标予以确认的待定常数.图3为辐射校正后的sar影像.

2.sar影像滤波方法

sar成像系统是基于相干成像原理的，在雷达回波信号中，相邻像素点的灰度会由于相干性而产生一些随机的变化，并且这种随机变化是围绕着某一均值而进行的，从而会影响从图像中提取后向散射系数值的精度。所述滤波采用增强lee滤波，滤波公式为：

其中，mij为滤波后的sar图像中第i行第j列的像元值，aij为滤波前的sar图像中第i行第j列的像元值，为滤波前的sar图像中以第i行第j列为中心的7×7窗口内像元值的均值，β＝1-cu²/cl²，cu＝0.5227，σij为滤波前的sar图像中以第i第j列为中心的7×7窗口内像元局部标准差.图4为影像滤波处理后的sar影像。

3.变化时间的确定方法

首先，将滤波后的sar影像按照成像时间排序形成时序影像s＝{mp|(p＝1，2，...，n)}，n为序列影像的个数，n＝8；然后，对于任意小于n的正整数k，将时序影像s＝{mp|(p＝1，2，...，n)}分成前序列s¹k＝{mp|(p＝1，2，...，k)}和后序列s²k＝{mp|(p＝k+1，...，n)}，后序列s²k的均值减去前序列s¹k的均值为dk，对于k＝1，2，...，n-1，可以形成时序差d＝{dk|(k＝1，2，...，n-1)}；最后，找到时序差d中的差异最大值dmax，其对应的序号即为地物发生变化时间t。图5为地物发生变化时间的图像.

4.变化差异图像生成方法

首先，利用地物发生变化时间t将时序影像s＝{mp|(p＝1，2，...，n)}分成变化前序列影像s¹t＝{mp|(p＝1，2，...，t)}和变化后序列影像s²t＝{mp|(p＝t+1，...，n)}；然后，将变化后序列影像s²t＝{mp|(p＝t+1，...，n)}减去差异最大值dmax生成新的变化后序列影像然后将变化前序列影像s¹t＝{mp|(p＝1，2，...，t)}和新的变换后序列影像组成归一化序列影像最后，利用计算公式cf＝dmax/std得到变化差异图像。变化差异图像icf是l×c的矩阵，l是变化差异图像的高度，c为变化差异图像的宽度，矩阵中每个元素为对应的变化因子。图6为变化差异图像.

5.地物变化信息提取方法

地物变化信息提取包含增强变化提取和减弱变化提取。增强变化是在地物变化过程中，后向散射系数值增大；减弱变化是在地物变化过程中，后向散射系数值减小。

增强变化提取的过程为：首先计算增强变化阈值th1，其计算方法为th1＝m1+δ1，m1为变化差异图像icf中大于0的像素值的均值，δ1为变化差异图像icf中大于0的像素值的标准差，本示例中m1＝1.89，δ1＝0.56；然后，对变化差异图像进行阈值分割，大于增强变化阈值th1的像素为增强变化，生成增强变化掩模图像r⁺，生成规则为：sij≥th1，r⁺ij＝+1；sij＜th1，r⁺ij＝0.sij为变化差异图像icf中第i行第j列的像素值，r⁺ij为增强变化掩模图像中第i行第j列的像素值，值为1则表示该像素发生了增强变化，值为0表示该像素未发生增强变化，i＝1，2，…，j＝1，2，....图7为增强变化的掩模图像，其中白色为增强变化提取结果。

减弱变化提取的过程为：首先计算减弱变化阈值th2，其计算方法th2＝m2-δ2，m2为变化差异图像icf中小于0的像素值的均值，δ2为变化差异图像icf中小于0的像素值的标准差，本示例中m2＝2.14，δ2＝0.82；。然后，对变化差异图像icf进行阈值分割，小于减弱变化阈值th2的像素为减弱变化，生成减弱变化掩模图像r^-，生成规则为：sij≤th2，r^-ij＝+1；sij＞th2，r^-ij＝0·sij为变化差异图像中第i行第j列的像素值，r^-ij为减弱变化掩模图像中第i行第j列的像素值，值为1则表示该像素发生了减弱变化，值为0表示该像素未发生减弱变化，i＝1，2，…，j＝1，2，....图8为减弱变化的掩模图像，其中白色为减弱变化提取结果.

需要指出的是根据本发明的具体实施方式所作出的任何变形，均不脱离本发明的精神以及权利要求记载的范围。