一种顾及统计特性的时序PolSAR影像无监督变化检测方法与流程

本发明属于遥感图像处理技术领域，是一种时序影像变化检测算法，具体地说是一种对时间序列全极化合成孔径雷达(polarimetricsyntheticapertureradar，polsar)影像变化检测的新方法。

背景技术：

在遥感影像处理中，变化检测是对同一位置进行多次观测并快速识别地物变化信息的一种技术手段。sar传感器不受光照和云雾天气影响，能够对地物进行长时间变化检测。无监督变化检测方法具有设计简单、处理速度快等特点，广泛应用于雷达影像变化检测的研究中。无监督变化检测方法主要通过预处理、差异影像提取和最优阈值选取，获得研究区域内地物变化检测的结果。其中，差异图影像的质量和最优阈值的选取与无监督变化检测的结果紧密相关。迄今为止，生成差异影像的方法大多集中在单和双极化sar数据中，如：对数比、邻域比、特征融合、基于马尔可夫随机场模型、主成分分析、kullback-leibler散度和googleearthengine时间序列分析法等。但这些方法由于信息利用不充分，导致差异影像提取精度不高。现阶段自动阈值提取方法能够有效地从差异影像中选取阈值。最小误差阈值选取方法基于贝叶斯理论，通过最小化分类误差确定最佳阈值，在变化检测阈值提取中被广泛应用。传统的最小误差阈值选取方法主要利用高斯分布对差异影像进行估计。由于生成差异影像的数据源不同，差异影像的统计特性往往不同。传统的方法并没有考虑polsar差异影像的概率统计特性，导致差异影像概率分布欠拟合，影响变化检测结果。因此，为了获得更为准确的变化检测结果，应采用更合适的分布模型对利用polsar数据生成的差异影像进行估计。

技术实现要素：

针对当前多时相sar影像变化检测中出现的极化信息利用不充分、差异影像概率分布欠拟合等问题，本发明利用多时相polsar影像的优势，提出了一种高效的时序polsar变化检测流程。相较于传统方法，利用时序polsar影像提供的丰富时空信息并结合综合检验统计方法，能够更为高效和准确地获取相同位置不同时刻的地物相似程度，通过拟合更为充分的概率密度函数改进最小误差阈值选取方法，提高阈值选择的精度。该方法适用于对不同性质的变化进行监测，不仅适合长时间地物变换分析，而且能够及时获取突变信息，为分析地物变化规律和应急响应提供了条件。

本发明的目的在于从根本上克服了sar影像无监督变化检测中存在的信息利用不充分和差异影像概率分布欠拟合的问题，提出利用综合检验统计和改进最小误差阈值选取方法对时序polsar影像进行变化检测。利用长时间序列polsar数据能够提供更为丰富的时间-空间信息，结合综合检验统计方法，能够更好地获取不同时相、相同区域地物的相似度。同时利用基于weibull和gamma分布改进的最小误差阈值选取方法，能够自适应地拟合差异影像统计分布并获得最优阈值，为高效获取精准的变化检测提供了基础。

本发明提供的一种顾及统计特性的时序polsar影像无监督变化检测方法，包括如下步骤：

步骤1，对时序polsar影像进行预处理；

步骤2，利用综合检验统计方法计算不同时相、同一区域地物的相似度，获得时序polsar差异影像k代表起始影像位置，n表示参与计算的时序polsar影像的个数；

步骤3，利用不同的概率分布函数对差异影像统计分布特性进行拟合，并选择weibull和gamma函数概率分布模型对最小误差阈值选取方法进行优化，获得优化后的判别函数j；

步骤4，计算判别函数j的最小值，并获得对应的最优阈值t^*；

步骤5，通过步骤2和步骤4获得的差异影像和最优阈值t^*，进行不同策略选择：若且k＝1，说明整个时间序列polsar影像相同位置的地物没有发生变化，即：x1＝x2＝…＝xn，索引为0并记录检测结果，进入步骤7；否则，说明整个时序polsar影像中至少存在一个变化，进入步骤6；

步骤6，利用综合检验统计方法计算不同时间内的差异影像并利用改进的最小误差阈值选取计算最优阈值s表示参与计算的影像间步长；

6a)若则令s＝s+1进入步骤6b；否则，说明在时间区间[k,k+s-1]内，相同位置的地物没有发生改变，但在k+s时刻，地物出现变化：xk+s≠xk+s-1，则cn＝cn+1，k＝k+s进入步骤2，其中cn代表时间序列中出现变化的个数；

6b)若k+s<n，则进入步骤6；否则，表明在[k,n]时刻内，相同位置的地物没有发生改变：xk＝…＝xn，但在[1,k]时刻内，存在cn个变化，进入步骤7；

步骤7，将时序polsar影像中的像素逐个进行处理，获得研究区域内地物变化检测结果，顾及统计特性的时序polsar影像无监督变化检测方法处理完毕。

进一步的，步骤1中，所述预处理包括：1)分别对各个时相polsar影像进行几何和辐射校正；2)利用同一参考影像对各时相polsar影像进行配准；3)对各时相polsar影像进行滤波处理；

具体为，通过利用nest，envi，polsarpro分别对各个时相polsar数据进行几何和辐射校正；选取同一参考影像，并利用强度互相关法，采用从粗到精的匹配策略对多时相polsar数据进行配准，使得配准精度达到亚像素级；对配准后的polsar数据进行精致lee滤波处理，一定程度上抑制相干斑噪声对变化检测的影响。

进一步的，步骤2中，时序polsar差异影像的表达式为式1，

其中，xi表示第i个时相极化sar图像的协方差矩阵；n、p、m分别表示参与计算的时序polsar影像的个数、维度和视数。

进一步的，步骤3中，最小误差阈值选取方法如式4所示，

其中ct(·)表示的为代价函数，分别表示在阈值t和差异影像灰度直方图条件下，未变化类别u和变化类别c的后验概率密度函数；

为了表示方便，令优化后的判别函数j(·)表达为：

其中l,t,分别代表灰度级，阈值和差异影像灰度直方图；γu,θu分别代表未变化类别后验概率密度的形状函数和尺度函数；γc，θc分别代表变化类别后验概率密度的形状函数和尺度函数，г(·)代表gamma函数。

进一步的，步骤4中，最优阈值t^*如式3，

通过遍历整个灰度级l，获得最小值j对应的最优阈值t^*，最优阈值t^*的步长选择根据处理影像的大小决定。

进一步的，步骤1中初始值为k＝1，步骤6中初始值为s＝1，cn＝0。

本发明的创新之处在于：1)提出了利用时序polsar信息进行无监督变化检测，能够提供更为丰富的时间-空间信息。2)提出了利用综合检验统计方法对长时间序列polsar数据进行差异影像提取的策略，提高了时序数据差异影像获取的效率，更为真实地表达相同位置、不同时刻地表的相似程度，并为检测发生变化的时刻和变化的次数提供了条件。3)提出了利用weibull和gamma概率分布模型拟合由时序polsar数据生成的差异影像，不仅能够更为精准地反映差异影像的统计分布特性，而且能够提高改进后最小误差阈值选取方法获取最优阈值的精度。4)本方法不仅能够适用于时序polsar影像的变化检测，同时也能较好地应用于单极化和双极化sar数据的变化检测，为观测区域大范围、长时间地监测提供了基础。

附图说明

图1是本发明的流程图。

具体实施方式

以下结合附图和实施例详细说明本发明的技术方案。

步骤1，通过利用相关软件nest，envi，polsarpro分别对各个时相polsar数据进行几何和辐射校正；选取同一参考影像，并利用强度互相关法，采用从粗到精的匹配策略对时序polsar数据进行配准，使得配准精度达到亚像素级；对配准后的polsar数据进行精致lee滤波处理，一定程度上抑制相干斑噪声对变化检测的影响。假设k代表起始影像位置，且0<k<n，cn代表时间序列中出现变化的个数；令初始值为k＝1，cn＝0。

步骤2，选取能够反映polsar数据特征的协方差矩阵以xi表示，考虑sar影像的统计特性，利用综合检验统计方法获取时间序列polsar数据差异影像其中时序polsar差异影像的表达式为式1。为了简化处理默认待处理的polsar数据视数相同。

其中，xi表示第i个时相极化sar图像的协方差矩阵；n、p、m分别表示参与计算的时序polsar影像的个数、维度和视数。

步骤3，利用weibull和gamma函数概率分布模型对差异影像进行拟合，并对最小误差阈值选取方法(如式4所示)进行优化，优化后的判别函数j(·)如式2所示。

其中ct(·)表示的为代价函数，分别表示在阈值t和差异影像灰度直方图条件下，未变化类别(u)和变化类别(c)的后验概率密度函数。

为了表示方便，令优化后的判别函数j(■)表达为：

其中l,t,分别代表灰度级，阈值和差异影像灰度直方图；γu,θu分别代表未变化类别后验概率密度的形状函数和尺度函数；γc，θc分别代表变化类别后验概率密度的形状函数和尺度函数，г(■)代表gamma函数。

步骤4，通过遍历整个灰度级l，获得最小值j对应的最优阈值t^*，其中最优阈值t^*如式3。

步骤5，通过步骤2和步骤4获得的差异影像和最优阈值t^*，进行不同策略选择：若且k＝1，说明整个时间序列polsar影像相同位置的地物没有发生变化，即：x1＝x2＝…＝xn，索引为0并记录检测结果，进入步骤7；否则，说明整个时序polsar影像中至少存在一个变化，进入步骤6。

步骤6，假设s表示参与计算的影像间步长，且初始值为s＝1，利用综合检验统计(式1)计算差异影像并利用改进的最小误差阈值选取方法式2和式3计算最优阈值

6a)若则令s＝s+1进入步骤6b；否则，说明在时间区间[k,k+s-1]内，相同位置的地物没有发生改变，但在k+s时刻，相同位置地物出现变化：xk+s≠xk+s-1，则cn＝cn+1，k＝k+s进入步骤2。

6b)若k+s<n，则进入步骤6；否则，表明在[k,n]时刻内，相同位置的地物没有发生改变：xk＝…＝xn，但在[1,k]时刻内，存在cn个变化。进入步骤7。

步骤7，将时序polsar影像中的像素逐个进行处理，获得研究区域内地物变化检测结果，顾及统计特性的时序polsar影像无监督变化检测方法处理完毕。

本文中所描述的具体实例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实例做各种各样地修改或补充或采用类似的方式替代。如预处理过程中的几何和辐射校正方法可以用其它软件完成，配准可以通过其它方法代替，精致lee滤波算法可以采用其它滤波算法代替；时序polsar影像的视数可以采用其它方法估计；等等。但这并不会超出本发明所提算法的框架，不会偏离本发明的精神，或者超越所附权利要求书所定义的范围。