PLS与EM算法支持下的遥感影像变化检测

郑志峰，曹建农，张雯佼

(1. 长安大学地球科学与资源学院，陕西 西安710054； 2. 陕西省测绘地理信息局，陕西 西安 710054)

早在1995年，土地利用、覆盖变化的研究就成为全世界学者研究的焦点和趋势[1]。主要目的是增强人们认识全球土地利用和土地覆盖变化的力度，而且重点提高估计这种变化的能力。随着卫星系统应用范围的不断扩大，遥感影像变化信息提取和变化检测已成为地理国情监测领域的重要内容[2]。其中，遥感影像变化检测技术作为测绘技术的重要分支，提升其技术水平已成为地理国情监测的迫切需要[3]。

遥感影像变化信息提取中构造差异影像与变化信息的阈值提取是两个主要的研究方向。对于普通的多光谱影像来说，影像各个波段间存在数据冗余与数据相关等问题，同时影像变化信息提取的效率与试验结果对差异影像的生成也带来一定的挑战。经典的遥感影像变化检测方法有：代数计算如影像比值、影像差值法，影像变化向量分析等[4-5]。Singh利用影像比值、影像差分、影像回归、归一化植被指数差分等方法进行热带雨林植被的变化信息检测，并采用不同时相的MSS影像对不同方法的结果进行对比分析，但直接对所有的波段影像进行运算没有消除多波段影像间的相关性。影像变换法主要包括缨帽变换(K-T变换)、主成分分析法(PCA)[6]、独立成分分析法(ICA)[7]、多元变化探测(MAD)[8]等。该类方法主要通过变换分析的思想，将影像中的差异信息或变化信息集中到几组不相关的变量中，进而提取影像间的变化信息，但变化检测结果对多时相影像间的辐射差异敏感、对变量的检测尺度依赖较大，但当多时相影像数据间有多重共线的情况时提取变化信息效果较差。对于差异影像变化检测阈值的提取，传统的方法是人工选择样本提取变化检测阈值，但样本选择的合理性难以确定，结果人为干扰因素比较大。文献[9—10]采用SVM算法确定变化区域与非变化区域的最优分类超平面，完成遥感影像变化信息的提取，但提取结果对样本数据的选择比较敏感。文献[11—12]以EM算法提取，采用PCA构造的差异影像变化检测阈值，提取变化信息,但国外学者BAZI[13]通过试验证明若提取区域非变化比例较大，变化检测阈值难以准确获取。考虑遥感成像复杂环境的影响因素及现实中地表覆盖变化率较小的情况，变化检测的具体阈值经常难以确定，如果仅仅考虑差异影像灰度值的整体统计问题，而不考虑影像整体变化率的影响，算法的效率、精度等较差，无法满足地理国情监测的需求。

针对上述问题，本文提出PLS和EM算法结合的局部遥感影像确定变化检测阈值的变化检测算法。首先，采用偏最小二乘算法去除多波段影像间的冗余信息，提取波段间的独立成分，有效地构造差异影像和集中变化信息，并依据PLS变换分析获得差异信息计算获得差异影像；然后，对整幅影像分块选取变化率最大的区域，采用最大期望的方法将多维数据进行分类，避免传统方法选择训练样本效率低的问题，提高变换检测的自动程度与试验结果的精度；最后，采用形态学算子对提取的变化信息进行后处理，减弱“椒盐现象”的影响，得到的变化检测结果更加符合地表的真实情况。

1 原理与方法

1.1 偏最小二乘原理

近年来，快速发展的一种新的多元数据统计分析PLS模型是通过最小化误差平方和找到一组数据的最佳函数匹配，在理论、方法和应用方面都有了快速的发展。PLS通过对多时相影像数据多元线性回归分析，降低影像数据本身的相关性，同时解决前后两组影像数据间的相关问题，有效解决多源变量间的多重共线问题。

假设有m个自变量{x1,x2,…,xm}和n个因变量{y1,y2,…,yn}。统计K个样本组成的数据集自变量X=[x1,x2,…,xm]K×m和因变量Y=[y1,y2,…,yn]K×n间的关系。为减少数据间差异，首先分别对前后影像数据进行预处理，获得处理后的前期影像数据矩阵E0和后期影像数据矩阵F0。在E0和后期影像数据矩阵F0中提取第一对主成分t1和μ1，且t1和μ1满足：

(1) 确保t1和μ1分量尽可能多地提取原始影像的信息。

(1)

回归方程的残差矩阵采用E1和F1表示；回归向量系数p1、r1应满足

(2)

(3)

1.2 EM算法

期望最大算法是无需外来数据和先验知识对有数据丢失的数据集进行最大似然估计的常用算法[14-15]。针对小比例变化量区域进行变化检测，通过样本选择的方法选择整幅影像中具有大比例变化量的局部区域作为训练样本，利用EM算法确定变化检测阈值。假设选择的局部影像变换后的数据集存在变化与非变化两个相互独立且整体上服从混合高斯分布的数据集。因此，运用EM算法求解极大似然进行参数估计的过程如下：

(1) 初始参数的确定：由于初始参数对EM算法结果影响较大，本文选取图上少量样本获得数据集的期望步的参数，即均值μk、协方差Σk、混合比例πk。

(2) 计算期望(E-step):假设模型参数已知的情况下求隐含变量Z，分别取Z1、Z2、…、Zk的概率，即在GMM中求数据点由各个成分生成的概率为

(4)

(3) 最大化期望(M-step):最大化E步得到的最大似然值重新估计分布参数为

(5)

(4) 极大似然函数估计

(6)

检查参数与似然函数是否收敛，迭代运算步骤(2)—(3)直至收敛，实现变化与非变化的聚类。

2 结果分析

2.1 试验区数据

本文采用杭州市2013年和2016年的高分一号卫星的多光谱影像数据进行试验。选取空间分辨率为2 m的全色波段与空间分辨率为8 m的多光谱数据进行融合。为了满足变化检测影像一一对应关系，首先，对前后时相影像进行空间配准，配准误差控制在0.5像素以内；然后，以2016年的遥感影像数据为基准，对2013年影像采用直方图匹配算法进行相对辐射校正；最后，参考文献[16]将配准后的多时相的影像进行叠加，并对叠加后的影像进行多尺度分割，将分割结果拆分成叠加前的遥感影像，由此分别得到2个时相影像的分割结果，综合考虑了2个时相影像的光谱和空间信息，保证了不同时相相同影像位置所对应的像斑包含的像元相同，并在各自的时相中具有光谱、空间同质性，利于进行像斑级的比较和分析。图1为杭州市2013年及2016年OLI传感器多光谱与全色影像采用NNDiffuse Pan Sharping法融合后321波段融合生成的真彩色影像。

图1 合成真彩色影像

2.2 PLS变换结果分析

图2 成分对应变量系数

分量均值的标准差标准差的均值典型相关系数信噪比PLS_13.3546.41×10-1035.128.57PLS_22.1298.06×10-1031.827.32PLS_32.1564.24×10-1029.396.78PLS_40.5121.47×10-109.160.96

2.3 不同方法的对比分析

由PLS_1、PLS_2、PLS_3提取的影像采用差值法构造差异影像，影像可以看作是由符合变化和非变化两类混合分布的像元组成,采用这个变化强度影像来进行阈值确定的试验，利用文中公式计算出局部EM算法的变化阈值为92，小于该阈值的影像标记为非变化。图3(e)为采用本文方法得到的变化检测结果图，其中白色表示变化区域，黑色表示非变化区域。试验结果表明，本文提出的变化检测算法能较很好地提取出耕地/荒地-道路、耕地/荒地-建筑物、水体-耕地植被覆盖等变化。

为了验证本文提出变化检测算法的有效性，将PCA-EM、MAD-EM及PLS-全局EM变化检测算法同本文试验算法进行对比。将试验结果同人工目视解译结果进行对比分析，采用漏检率、误检率评价变化检测精度[17]。图4为不同变化检测算法的试验结果图。

图3 不同算法的变化检测结果

表2给出了不同算法的变化检测精度。由表2可知，文中提出的基于PLS和局部EM算法变化检测结果的正确率最高，漏检率和误检率最低。这是由于PCA方法仅考虑影像内部的相关性，缺少对多时相影像间的相关性的考虑，而且PCA-EM算法对多时相影像相对辐射校正的结果有较高要求;对多时相影像的多波段的重组影像进行MAD变换，虽然消除了多时相影像间的相关性，解决了构造差异影像问题，但是难以解决变化信息有效集中的问题。PLS变换考虑多时相影像间的线性关系，提取影像间相关性较高的成分，有效地消除了影像间的多重共线性，从而解决了变化信息有效集中的问题，为整体上提高变化检测正确率提供必需信息。另外，将局部EM(最大期望)确定阈值方法与全局EM确定阈值方法进行对比，试验结果表明全局EM阈值确定方法对差异影像直方图拟合效果差，而局部EM阈值方法对于小比例变化区域，能准确确定变化检测阈值，减少了自动阈值算法的不确定性。总之，本文提出的PLS-局部EM算法变化信息基本都能检测出来，更加符合真实情况,提取信息结果比较理想。

表2

3 结 语

本文针对多时相影像间多重相关性提出了一种PLS与局部EM结合的高分辨率遥感影像变化检测算法。该方法首先采用多时相影像PLS变换分析提取主成分影像构造差异影像，然后利用局部变化比例较大的影像确定变化检测阈值，最后根据构造的差异影像与局部EM算法提取的阈值提取变化信息。试验表明，文中提取的算法能够提取绝大多数地物类别变化信息，与传统的算法相比检测精度最高，同时解决了小比例变化区域阈值确定问题。但文中未对地物类别信息进行变化判断及变化趋势的预测，这些内容有待于进一步研究。