2. 中国科学院遥感与数字地球研究所，北京 100101；
3. 长安大学理学院，陕西 西安 710064
2. Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100101, China;
3. College of Science, Chang'an University, Xi'an 710064, China
遥感技术因具备大面积同步观测能力、获取的数据具有快速、客观等特征，而被广泛用于土地覆被变化、城市监测、灾害评估等领域。作为上述各类应用的基础，变化检测是在遥感影像配准基础上，通过一定的相似性度量确定变化区域，提取地物变化信息[1]。根据实现方法和检测结果的差异，可以分为非监督和监督两类方法[2]。
非监督方法是在构建两期影像相似性(差异)度量函数基础上，根据变化区域的光谱差异大小，通过阈值分割确定变化区域。因此，相似性度量和阈值选取是其中的关键。在相似性度量方面，差值、比值是最先使用的相似性度量；随后，克服不同时相影像间辐射特征差异的归一化植被指数[3]以及利用多波段信息的主成分分析[4]等被相继引入；变化向量分析将影像从灰度空间转化为极坐标空间，利用极径和极角分别表示变化的幅度和类型[5]；而针对仅利用光谱信息的不足，学者们发展了考虑局部区域空间特征的变化检测方法，如Hopefield神经元方法[6]；通过多证据融合，将纹理、形状等多种高级特征同时用于变化检测的相似性度量[7]。在阈值选取方面，直方图分割、双(多)阈值方法简单，且在很多应用中取得了较好的效果[8]，但变化地物与未变化类别混合严重时精度较低，对此，能够自动学习影像变化特征的水平集[9]等智能方法得到应用。
监督方法是在选取变化样本的基础上，利用神经网络[10]、蚁群优化[11]、慢特征分析[12]等进行特征训练和检测，以提取变化区域。分类后比较是其中常用的一类监督变化检测方法，它是在两期影像选取样本后进行监督分类，再通过比较分类结果来确定地表变化信息。受误差传递的影响，精确的变化检测前提是两期影像均能被正确分类[13]，为此，动态证据理论等用于减少不同时相分类结果合成的不确定性[14]。为了减少样本选择的工作量，主动学习、复合分类等方法也被用于不同时相的影像之间的亮度变化模式的学习中[15]，而潜分析[16]、频繁序列分析[17]等方法则被用于提取时间序列影像上的土地覆被变化语义信息，是一种广义的变化检测方法。
从以上分析可以看出，监督方法大部分能够获得地表覆被转移矩阵，主要面向精确的地表覆被演变分析，而非监督方法则将地表快速划分为变化和未变化区域，在缺乏研究区先验知识的情况下，是从遥感影像中快速提取变化区域的必要手段，而且可进一步挖掘和分析地表覆被的变化信息。然而，当前的非监督方法普遍是对整景影像采用相同的相似性度量和变化提取规则，成功应用的前提在于不同类型的地物变化具有相同或者相似的亮度变化幅度。而辐射条件、植被时相变化、地物变化等多种因素导致变化与未变化地物是交织在一起的，限制了传统的变化检测方法的精度。鉴于此，本文将影像转换到联合概率密度空间，以便自适应地确定变化区域，提高变化检测精度。
1 本文方法变化检测是要利用辐射特征，准确识别出地表覆被类型发生改变的区域。而根据应用的不同，一些方法探索自然地物的时间变化，而部分方法则是研究地表类型的改变，本文仅讨论后一种情形，即考虑地物类型发生变化的情况。
在不考虑地物变化的理想情况下，同一地区不同时间获取的两期影像可假定是独立同分布的(independent and identically distributed，IID)，即两期影像亮度所组成的变量序列具有相同的概率分布，且互相独立。而实际上，被动遥感成像是一个复杂的电磁辐射与地表相互作用的过程，受多种因素的综合影响。
遥感器接收的地物反射辐射度量R可表示为
式中，Es表示辐射、大气、传感器畸变等系统因素，所有像素服从相同的变化规律；Et表示植被、水体(受叶绿素、泥沙等因素的影响)等自然地物随时相的改变；Er表示人类和自然作用使得地表覆被发生改变；f表示上述因素的综合作用。
阈值法非监督变化检测是基于Er导致的亮度变化相对于Es和Et较大的假设，依据一定的优化函数来寻找最优的分割阈值。很多情况下，辐射校正后的残差(∈Es)、植被、水体等时相变化(∈Et)与地物变化(∈Er)导致的地物辐射特征变化相近，从而使得上述多种类型引起的辐射变化交织在一起，阈值法难以有效区分。
联合概率密度表达了地物在两期影像上的亮度值变化的统计规律[18]，相对于传统的相似性度量具有更丰富的信息，已成功用于变化检测[18]、大区域影像镶嵌的色彩一致性处理[19]等领域。基于此，本文提出一种联合概率空间的自适应非监督变化检测方法，试图利用联合概率空间的多阈值克服传统阈值分割方法无法区分不同因素造成像元灰度变化的缺陷，减少对自然地物时相变化的误检，提高变化检测精度。
为了叙述方便，定义如下符号：It1和It2为同一地区两期不同时间获取的影像，t1<t2，It1和It2灰度值x和y构成的空间分别为X和Y。变化检测的目标为提取It2相对于It1的变化区域(Φ)和未变化区域(Ω)，变化结果利用二值影像IR表示，0表示未变化，1表示发生变化。
1.1 联合概率密度空间及其变化定义遥感影像以扫描方式记录了地表的电磁辐射特征，受Es、Er和Et等因素的影响，可以将其看成是描述地表电磁辐射的一个随机变量。因此，两期影像It1和It2进而可看成是描述地表辐射特征的二维随机变量，表征了对应波段之间的亮度值的转换关系，其分布特征便可用联合概率密度表示。由于影像获取过程中的采样对地表辐射值进行了离散化，所设的联合概率密度事实上是一个离散型的二维随机变量，因此可表示为
式中，P(x,y)和N(x,y)表示影像It1和It2对应位置处亮度值分别为x和y的像素出现的概率和频数，概率密度P(x,y)构成了一个以X和Y为值域的二维空间，即联合概率密度空间。
图 1是某地区2006年和2008年一组ALOS影像近红外波段的联合概率密度图，受到Es、Et和Er 3种辐射特征变化因素的影响，不同区域的概率密度值随之发生变化(颜色由红到蓝表示密度逐渐降低)。为了说明问题的方便，将联合概率密度空间标记出A、B、C 3个区域。其中，在A区域，It1和It2两期影像的辐射亮度值相似；形态上以直线l为轴，基本对称分布，近似为椭圆形状；椭圆中心概率密度较高，随着远离直线l而逐渐降低。该区域主要是地表类型保持稳定的区域，仅在远离直线l的边缘有少量变化的地物。B、C区域远离直线l分布，形态不规则，例如B区域在It1上亮度值相对集中，而It2分布分散；C区域恰好与B区域相反，该区域内地物亮度变化大，概率密度较小且不规律，主要为变化地物。
|
| 图 1 联合概率密度图 Fig. 1 The joint probability space |
根据以上分析，联合概率密度分布空间上的变化地物可定义为概率值较小[18]，且分布在远离直线l的空间区域点，即离群点(outlier)。据此定义，变化检测问题转化为从联合概率密度空间进行离群点检测。为此，本文提出一种迭代方法确定离群点以提取变化像素。需要说明的是，本文为方便起见将联合概率密度空间标记为A、B、C 3个区域，并利用多阈值方法提高变化检测精度，而并不需要区分各种因素导致的地物变化。
1.2 迭代的变化区域检测如上所述，将联合概率密度空间上地物变化的像素看成是远离直线l的离群点。本文通过a倍均方差迭代识别技术实现变化像素的检测，具体实现流程包含3个关键步骤：提取对称轴l，确定l上每个亮度值处的变化边缘以及计算影像上的变化区域。详细实现流程如下。
1.2.1 提取对称轴l记对称轴所在直线l为y=kx+b，其中y为影像It2上的亮度值，x为影像It1上的亮度值。根据直线l上概率密度之和尽可能大以及直线两侧的概率尽可能均匀两个准则，采用枚举法提取直线l，即从联合概率密度空间任取两点，将满足上述两个条件的最优解作为结果，详细实现流程参见文献[19]。
1.2.2 确定l上每个亮度值处的变化边缘对离群点进行去除，过程如图 2所示，记P为直线l在联合概率密度空间上整数点集(∀pi∈P)，垂直于l且通过整数点pi的直线为l′pi。其在概率密度空间的离散点集Qpi，坐标(x′1,y′1),(x′2,y′2),…,(x′m,y′m)，对应的概率密度W=(w1,w2,…,wm),YQpi的加权均值(μy)，标准差(σy)为
|
| 图 2 变化区域提取方法 Fig. 2 The change area extraction method |
按照式(6)和(7)确定离群点
式(6)表示对于∀qi∈Qpi，计算距离μy最远的点；式(7)表示距离均值μy在a倍标准差之外的点，a是一个待设置的参数，设置方法及其对结果的影响在试验部分予以详细讨论。
将同时满足式(6)和式(7)的点作为离群点进行剔除，然后根据式(3)、(4)、(5)重新计算(μy,σy)，若没有满足式(6)和式(7)的离群点或者达到最大迭代次数，则迭代结束。记迭代中止时Qpi上、下边缘的两个点分别为qpiT、qpiB。
对于∀pi∈P，都得到一组点与之对应，将所有与pi对应的点按照以下顺序连接，构成封闭的空间A
轨迹A内部封闭区域为地物类型不变的区域，而轨迹外的区域为地物变化区域。
若采用传统的单一阈值，确定的变化区域如图 2中虚线所示，虚线外部区域为变化像素，而内部区域为未变化像素。本文方法相当于一种多阈值的方法，根据不同灰度值处的亮度分散情况(标准差)自适应确定分割阈值，从而精细地刻画不同类别地物受辐射特征变化的差异影响，以期获得更高的二值化精度。
1.2.3 计算影像上的变化区域经过步骤(3)之后，变化检测结果在联合概率密度空间上表示为
对It1和It2对应位置上亮度值为(x,y)的点进行判断，若位于A区域内，则作为不变像素，并将结果用1表示，不变像素构成的集合为Ω；反之，作为变化像素，用0表示，并构成变化像素集合Φ。
2 试验 2.1 数据与试验方法为了验证算法的有效性，选取大小分别为1000×1000像素和800×800像素的LandSat TM5影像和ALOS多光谱影像各一组进行试验，分别记为试验1和试验2。试验1的两期LandSat TM5影像获取时间分别为2002年11月4日和2009年11月23日，位于湖北省宜昌市的三峡大坝地区，土地覆被变化主要由大型水利枢纽设施修建和蓄水引起，造成地面低洼处被淹没及其相关变化；试验2的ALOS影像获取日期分别为2007年9月31日和2009年10月31日，位于浙江省诸暨市，地表覆被变化驱动因素主要为城市建设将耕地、林地、水面等开发为城市建设用地。两组试验代表了土地覆被变化和城市监测的典型应用，LandSat影像6、4、2波段的假彩色合成及ALOS影像的真彩色分别如图 3所示。两组影像均进行了几何校正，配准精度达到亚像素级。
|
| 图 3 试验使用的原始影像 Fig. 3 The original images of experiments |
由于采用直线l对称轴作为距离标准，相当于进行了一次相对辐射校正，能够消除两期影像上的线性辐射差异，因此，本试验直接采用原始灰度值而未进行相对辐射校正。使用上述3个波段逐波段变化检测，获得二值结果，然后通过逻辑乘运算得到最终检测结果，即对于某一像素所有波段均检测为变化才最终判定为变化，这种方法对于椒盐噪声具有一定的抑制作用，能获得较平滑的处理结果。
为了对变化检测结果进行评价，采用目视解译方法，手动提取地物变化信息，作为地表变化真值。试验2中地表覆被类型变化复杂，单纯的目视解译难以获取准确的变化信息，笔者选取了Google Earth上试验日期相近的2007年2月2日和2009年10月2日两期高分辨率数据来辅助判别。
2.2 试验结果及其比较试验中，设置a=2，最大迭代次数为100。其中，试验1检测结果如图 4(a)所示，红色表示变化区域而灰色表示背景(未变化区域)。可以看出：水利设施修建、蓄水造成的地表淹没、下游泥沙淤积的地表覆被变化均得到了有效检测。值得注意的是，右上侧的沙坪水库因蓄水原因造成的水面变化也得到有效检测。图 4(b)为目视方法的变化检测结果，红色和灰色分别表示正确检测的目标和背景、绿色表示漏检、蓝色表示误检，通过对比可知，孤立的误检测和漏检测的情况较少，主要出现在变化与未变化区域的过渡区域。
|
| 图 4 变化检测结果 Fig. 4 The change detection result |
试验2中地表覆被类型变化多样，如农地、林地转为建设用地(道路和建筑物修建)、河道疏通、湖泊筑坝等导致水面的变化。采用与试验1相同的渲染方法，本文方法与目视解译的检测结果如图 4(c)和图 4(d)所示。可以看出：各种地表覆被类型变化区域均得到了有效检测，误检测和漏检测情况较少。
在信息提取和目标识别中，一般用P表示目标物的数目，而N表示非目标物的数目，采用TP表示算法检测为目标物且经验证为目标物的数目，TN表示算法检测为非目标物且经验证为非目标物的数目，FP是算法检测为目标物而经验证为非目标物的数目(误检)，FN表示算法检测为非目标物但经验证为目标物的数目(漏检)。在此基础上定义了正确率(TR)、误检率(FAR)、漏检率(OAR)3个指标，即
下面将本文方法与变化向量分析法以及文献[18]中联合概率密度方法进行算法精度比对。变化向量分析[5]是通过将影像从灰度空间转为极坐标空间，并寻找阈值进行非监督变化检测，记为CVA方法；文献[18]方法也是将变化定义在联合概率密度空间，以此作为相似性度量，并采用局部混合信息指标自动确定阈值提取变化区域，记为PD方法。采用上述两种方法及本文方法，在相同的影像上进行试验，以像素为单位，统计上述3个指标，其中a值分别设置为1.5、2、2.5，结果如表 1。由表 1可以看出：本文方法在a=2时，TR、FAR、OAR 3个指标均优于对比方法；而a=1.5和a=2.5时，TR与对比方法相近。
| 方法 | a | 试验1 | 试验2 | |||||
| TR | FAR | OAR | TR | FAR | OAR | |||
| CVA | 90.6 | 9.4 | 8.8 | 85.3 | 14.7 | 8.5 | ||
| PD | 91.5 | 8.5 | 8.6 | 87.6 | 12.5 | 7.7 | ||
| 本文方法 | 1.5 | 90.0 | 10.0 | 13.3 | 88.1 | 11.9 | 12.1 | |
| 2 | 93.9 | 6.1 | 8.3 | 91.9 | 8.1 | 6.7 | ||
| 2.5 | 91.3 | 8.7 | 4.3 | 88.1 | 11.9 | 4.6 | ||
进一步选取典型区域，比较各种地物的变化检测细节，图 4(a)中L区域放大后如图 5(a)—(c)所示，可以看出本文检测结果与水体边缘结合紧密，具有较高的定位精度。特别注意白框内区域，山体阴影是水体检测中的混淆地物，CVA方法由于仅利用极径作为变化的衡量指标，造成了漏检，而本文方法充分利用不同地物变化幅度，能够有效检测；图 4中M区域如图 5(d)—(f)所示，变化为湖泊疏浚和筑坝，由于待检测影像上湖面部分区域存在植被，光谱特征变化相对于其他区域较小，PD方法利用相同阈值对整景影像检测，造成了漏检，同时左上角绿地转为裸地的变化也被漏检；图 4中N区域如图 5(g)—(i)所示，变化由建筑物修建引起，靠近边缘处辐射特征变化幅度相对较小，PD方法发生了漏检。以上说明本文方法能够区分灰度变化较小的地物，具有更好的灵敏性和稳健性。
|
|
注：(a)、(d)、(g)表示变化前的影像与本文方法检测结果的叠加；(b)、(e)、(h)表示变化后的影像与本文方法的变化检测结果的叠加；(c)、(f)、(i)为本文检测与对比方法的比较，红色区域内为对比方法结果，黄色为本文方法结果，(c)为CVA方法，(f)、(i)为PD方法) 图 5 局部区域的叠加放大显示 Fig. 5 The change detection result of local details |
参数a的值是本文方法成功的关键之一，为了分析和评价a值对结果的影响，将a设置在区间[1,3]，以步长0.1进行变化。试验1和试验2的正确率、误检率和漏检率变化情况如图 6，可以看出：①a=2.0左右(试验1，a=1.8；试验2，a=2.2)，正确率最高、误检率最低，结果具有最好的可信度；②随着a减小，漏检率逐渐降低，也造成误检率逐渐升高；③由于直线l附近联合概率密度较高，随着a值较小，漏检率、正确率和误检率变化速率减慢；④a为2附近的较宽区域内，结果的正确率均较高，与传统方法相比具有相似的精度，且具有较高的鲁棒性。
|
| 图 6 参数a对试验结果的影响 Fig. 6 The result accuracy effected by parameter a |
根据以上分析，为获得较高的检测正确率，建议设置a为2。另外，在实际应用中，也可以根据不同的应用需求，调整a的值，以获得不同的误检率和漏检率，如城市违建执法检查中需要尽可能找出靶区，避免漏检，此时可以将a值设置较大，尽可能找出发生变化的疑似区域，进而通过实地察看验证计算机检测结果。
3 结 论针对单阈值非监督变化检测的不足，本文提出了一种联合概率密度空间的自适应变化检测方法，将变化像素定义为联合概率密度空间的离群点，通过迭代确定多个最优阈值进行识别。与当前原理相近的方法比较，本文方法由于充分考虑了不同地物的辐射特征变化幅度，能够一定程度克服两期影像上的非线性辐射差异，在正确率、误检率和漏检率方面具有一定的优势，同时，所需设置的参数较少且检测结果鲁棒性较好，是一种较为理想的变化检测方法。
然而，通过上述试验探索，发现本文研究也存在一些不足，表现在：①本文方法虽然利用统计模型并取得了较好的效果，但对其数学物理意义的深入分析和讨论不足；②本文方法是逐波段进行变化检测后再通过逻辑乘获得最终结果的，更理想的方式是将n波段影像作为一个n元随机变量，然后在n×n维联合概率密度空间中去定义变化区域并发展变化提取方法；③自然地物时相变化在很多应用中是研究的对象，而本文方法不能区分该类型的变化，同时，本文方法是一种逐像素的非监督变化检测，未能充分利用对象级的变化特征及其挖掘变化所表达的地学意义；④试验虽然选择了不同分辨率、不同变化类型的数据，并表现出了较好的精度和稳健性，但遥感数据及地物变化模式多样，方法的广泛适用性仍有待进一步研究和探索分析。
| [1] | LU Dengsheng, MAUSEL P, BRONDÍZIO E, et al. Change Detection Techniques[J]. International Journal of Remote Sensing, 2004, 25(12): 2365-2401. |
| [2] | RADKE R J, ANDRA S, AL-KOFAHI O, et al. Image Change Detection Algorithms: A Systematic Survey[J]. IEEE Transactions on Image Processing, 2005, 14(3): 294-307. |
| [3] | LYON J G, YUAN Ding, LUNETTA R S, et al. A Change Detection Experiment Using Vegetation Indices[J]. Photogrammetric Engineering and Remote Sensing, 1998, 64(2): 143-150. |
| [4] | BYRNE G F, CRAPPER P F, MAYO K K. Monitoring Land-cover Change by Principal Component Analysis of Multitemporal Landsat Data[J]. Remote Sensing of Environment, 1980, 10(3): 175-184. |
| [5] | BOVOLO F, BRUZZONE L. A Theoretical Framework for Unsupervised Change Detection Based on Change Vector Analysis in the Polar Domain[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(1): 218-236. |
| [6] | GHOSH S, BRUZZONE L, PATRA S, et al. A Context-sensitive Technique for Unsupervised Change Detection Based on Hopfield-type Neural Networks[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(3): 778-789. |
| [7] | 汪闽, 张星月. 多特征证据融合的遥感图像变化检测[J]. 遥感学报, 2010, 14(3): 558-570. WANG Min, ZHANG Xingyue. Change Detection Using High Spatial Resolution Remotely Sensed Imagery by Combining Evidence Theory and Structural Similarity[J]. Journal of Remote Sensing, 2010, 14(3): 558-570. |
| [8] | 胡召玲. 广义高斯模型及KI双阈值法的SAR图像非监督变化检测[J]. 测绘学报, 2013, 42(1): 116-122. HU Zhaoling. An Unsupervised Change Detection Approach Based on KI Dual Thresholds under the Generalized Gauss Model Assumption in SAR Images[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(1): 116-122. |
| [9] | BAZI Y, MELGANI F, AL-SHARARI H D. Unsupervised Change Detection in Multispectral Remotely Sensed Imagery with Level Set Methods[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(8): 3178-3187. |
| [10] | 戴芹, 马建文, 欧阳赟, 等. 利用贝叶斯网络进行遥感变化检测[J]. 中国图象图形学报, 2005, 10(6): 705-709. DAI Qin, MA Jianwen, OUYANG Yun, et al. Remote Sensing Change Detection Using Bayesian Networks[J]. Journal of Image and Graphics, 2005, 10(6): 705-709. |
| [11] | 戴芹, 刘建波, 刘士彬. 微粒群优化方法的遥感影像变化检测研究[J]. 测绘学报, 2012, 41(6): 857-863, 885. DAI Qin, LIU Jianbing, LIU Shibing. Remote Sensing Image Change Detection Using Particle Swarm Optimization Algorithm[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(6): 857-863, 885. |
| [12] | WU Chen, DU Bo, ZHANG Liangpei, et al. Slow Feature Analysis for Change Detection in Multispectral Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(5): 2858-2874. |
| [13] | SERRA P, PONS X, SAURÍ D.Post-classification Change Detection with Data from Different Sensors: Some Accuracy Considerations[J]. International Journal of Remote Sensing, 2003, 24(16): 3311-3340. |
| [14] | LIU Zhunga, DEZERT J, MERCIER G, et al. Dynamic Evidential Reasoning for Change Detection in Remote Sensing Images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(5): 1955-1967. |
| [15] | DEMIR B, BOVOLO F, BRUZZONE L. Detection of Land-cover Transitions in Multitemporal Remote Sensing Images with Active-learning-based Compound Classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(5): 1930-1941. |
| [16] | VADUVA C, COSTACHIOIU T, PATRASCU C, et al. A Latent Analysis of Earth Surface Dynamic Evolution Using Change Map Time Series[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(4): 2105-2118. |
| [17] | JULEA A, MÉGER N, BOLON P, et al. Unsupervised Spatiotemporal Mining of Satellite Image Time Series Using Grouped Frequent Sequential Patterns[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(4): 1417-1430. |
| [18] | GUEGUEN L, SOILLE P, PESARESI M. Change Detection Based on Information Measure[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(11): 4503-4515. |
| [19] | 吴炜, 沈占锋, 李均力, 等. 联合概率密度脊提取的影像镶嵌色彩一致性处理方法[J]. 测绘学报, 2013, 42(2): 247-252. WU Wei, SHEN Zhanfeng, LI Junli, et al. Ridge of Joint Probability Density Based Color Normalization Method for Image Mosaic[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2): 247-252. |