李莎， 倪维平， 严卫东， 吴俊政， 张晗

(西北核技术研究所，西安 710024)

0 引言

多光谱图像的变化检测是对不同时间获取的多光谱遥感图像进行定量分析，并确定图像区域内地表变化特征和过程的技术[1]。从技术处理流程来看，一般有图像预处理、变化信息发现、变化区域提取与变化类型确定等过程[2]，其中最为关键的是变化信息发现，大多数研究都是围绕这一问题开展的。近年来，国内外相继发展了许多适用于遥感图像变化检测的方法[3-6]，但没有一种方法是最优的，每种方法都只具有一定的适用范围。为了解决传统方法(图像比值法、变化向量分析法、主成分分析法[7-8]等)应用于多波段遥感图像变化检测中的不足(抑制噪声、消除相关性等)，丹麦学者Nielsen等[9-11]率先提出了多元变化检测(multivariate alteration detection，MAD)的概念和方法。廖明生等[12]和陈垒等[13]将这种方法应用于多时相多波段遥感图像变化检测，实验结果也显示出了该方法的明显优势和应用潜力。基于MAD方法，Nielsen[14]又提出了迭代加权多元变化检测(iterative weighted multivariate alteration detection，IRMAD)方法; Canty等[15]通过实验证实，该方法在干旱或沙漠地区以及场景中变化部分所占比例相对较大的多时相多波段遥感图像变化检测中，效果优于MAD方法； 但在变化部分所占比例相对较小的非干旱或沙漠地区，并不能达到很好的效果。为此，本文提出了一种基于选权迭代估计(iterative estimation with weight selection，IEWS)与非监督分类(unsupervised classification，UC)的多光谱图像变化检测方法。通过借鉴IEWS的思想[16]，并以类似于IRMAD的迭代模式进行回归估计，剔除掉大部分的非变化信息，得到初始变化信息； 对于一个非单一场景，由于地物的多样性，2个时相图像中不同地物的非变化信息之间存在不同的线性关系，因而通过基于初始变化信息的分类处理，以及对不同类别的IEWS，剔除残存的非变化信息，得到最终的变化信息。利用真实TM图像处理的结果验证了本文方法的有效性，所得到的变化信息在空间位置上同该区域相应时间段内土地覆盖类型变化情况具有很好的一致性； 同时与采用MAD，IRMAD方法的变化检测结果相比较，表明本文方法对相对较小的变化信息具有更好的变化检测能力。

1 算法的基本原理

借鉴IEWS思想，本文先用最小二乘法进行参数估计； 每次估计后，用差值D计算下一次迭代中图像每一个像元X的权值P，使非线性变化像元的权值越来越小，直到趋近于0，这样发生非线性变化的像元就被排除在估计样本之外； 最后获得的估计参数将不受变化信息的影响。IEWS的计算过程如图1所示。

图1 选权迭代估计的计算步骤

通过一次IEWS获得的变化信息将是图像中可能发生变化的所有像元值。由于地物类型的多样性，并不能一次将非变化信息完全剔除。为了从初始变化信息中找出更为精确的土地覆盖变化信息，需要对差值D做进一步处理。对于同类地物，其所在的前后2个时相的遥感图像应该具有很好的线性相关性，因而通过基于变化信息进行的分类，并对每类像元进行IEWS，则可以剔除掉每类像元中的线性变化部分。所有N种类型地物非线性变化像元合集即为2个时相图像之间真正的变化信息。变化信息获取过程如图2所示。

图2 基于选权迭代估计与非监督分类的变化检测流程

2 权函数的选取

如何确定权值是上述估计过程中需要解决的关键问题。本文利用类似于IRMAD的迭代模式[17]来确定每次迭代过程中的权值P。

2个时相图像X,Y都是多波段图像，其差值图像D包含有多个图像层。通常假设差值图像中每一层图像的像元值近似服从高斯分布，并且原始图像的每个波段之间相互独立，从而每个像元的D值的平方和归一化方差将近似服从自由度与波段数n相等的Chi-square分布(即“卡方分布”)。在曲线拟合中，可使用Chi-square作为模拟拟合优劣的一个指标，Chi-square越小越好。

对于第j个像元，cj为Chi-square分布中的一个样本，有

,

(1)

式中：n为波段数；σDi为第i个波段差值图像Di的标准差。进而可将像元不变化的概率Pr[5]表示为

Pr=1-Pχ2;n(c)

,

(2)

式中：χ2为自由度为n的Chi-square分布；c为变化的像元数，c值越小，不变化概率越大。基于此，迭代过程中将像元不变化概率Pr作为观测值的权值，进行加权迭代计算。这样，变化信息的权值将逐渐趋近于零，含有变化信息的观测值将被剔除，则最后获得观测值的权等于不包含变化信息的观测值的权。通过IEWS后获得的差值图像即为2个时相图像之间的初始变化信息。

3 分类方法的选择与实现

对于分类方法的选择，目前常用的遥感分类方法主要包括监督分类与非监督分类。由于监督分类法需要事先确定训练样本数据，而在很多情况下训练样本难以获取； 而非监督分类直接根据待分类数据个体间的相似性测度进行类别划分，无需训练样本数据。为了实现较高的自动化程度，本文选用非监督分类ISODATA算法进行分类。ISODATA算法是在没有先验知识的情况下对数据进行分类，能够吸取中间结果所得到的经验，具有自组织性，能够自动地进行类别的“合并”和“分裂”，其各个参数可在聚类调整中逐渐确定，并最终构建所需要的判别函数，且运算速度与精度均较好[18]。为此，本文选用ISODATA法作为分类算法。

4 实验结果与分析

本文选取香港维多利亚港周围区域(覆盖了九龙半岛及香港岛和青衣岛的一部分)进行研究。从20世纪90年代开始的香港填海造地等工程使该区域内的土地覆盖和利用均发生了很大的变化。为了尽量避免植被覆盖、气候变化等因素对真实变化的影响，用作变化检测的2景图像的成像时间处于不同年份的同一月份。2景图像均为Landsat5 TM图像，空间分辨率均为30 m，成像时间分别为1989年11月和2005年11月。选用除热红外波段(TM6)以外的6个可见光和近红外波段，对原始图像进行几何配准后截取出360像元×380像元的区域作为研究对象(图3)。

(a) 1989年11月20日(b) 2005年11月23日

对2个时相图像运用本文方法进行变化检测，通过一次IEWS得到的初始变化信息图像如图4(a)所示，颜色很暗或很亮的地方为变化区域。对初始变化信息进行非监督的ISODATA分类，并对每类像元单独进行估计，得到每类地物的变化信息，如图4(b)所示。经过分类后的再次IEWS，进一步清楚地显示出感兴趣的变化信息，使得变化区域与非变化区域的对比更为明显。

(a) 第一次估计结果(b) 最后估计结果

本文选用可分离性测度Jeffries-Matusita distance,即JM距离来评价图4中2次估计结果的有效性。根据先验知识，选取了2个典型样本区(图4)，典型样本区1在2个时相期间内由水体变为陆地，将其作为变化样本； 典型样本区2在2个时相都为林地，故将其作为非变化样本。2次估计结果在各波段上2类样本间的JM距离如表1所示，除波段6外，最后估计结果的JM距离都大于第1次估计结果，这说明最后估计结果具有更好的可分离性。

表1 2次估计结果在各波段上2类样本间的JM距离

将最后变化图像分为2部分，图5(a)为像元灰度值相对变大的情况，图5(b)为像元灰度值相对变小的区域，分别对应于第2时相遥感图像相对第1时相遥感图像变亮和变暗的区域。

(a) 灰度值相对变大区域(b) 灰度值相对变小区域

图5(a)中的变化像元显示出了1989—2005年间该区域因填海造地工程而新出现的填海土地、码头及桥梁等； 而在图5(b)中，红色圆圈标记出的地块，则显示在该期间土地覆盖情况的变化导致了遥感图像亮度变暗。这表明本文方法可以分辨出每个波段上像元灰度值是相对变大还是变小了，有助于对变化情况做进一步的分析。

对2个时相遥感图像使用MAD方法得到6个检测图像(图6)。可以看出，仅MAD5图像集中反映出区域内填海造地工程带来的变化，其他图像反映的则是少数次要变化和大量噪声成分。

(a) MAD1(b) MAD2(c) MAD3

(d) MAD4(e) MAD5(f) MAD6

采用本文方法最后得到的变化像元单波段图像与MAD5图像使用相同方法进行变化区域提取的变化图像的对比如图7所示。本文方法所得到的结果在表达影像细节方面(图7中绿色圆圈所标注的桥梁和码头等处)有更好的效果。

为进一步说明本文方法对细微目标的检测能力，选取图7(a)中间和右侧绿色标注出的2个区域作为典型区域，分别采用本文方法和MAD5图像进行比较。先对这2个典型区域进行目视解译获取其真实变化情况，再将本文方法及MAD5图像进行阈值提取后的变化结果与真实变化情况进行对比，结果如表2所示。表中MAD5的提取结果2是该方法误检率和漏检率均为最低情况下的统计结果。将其与本文方法提取结果比较可知，本文方法有更低的误检率和漏检率，能更好地区分变化与非变化信息。

(a) 变化区域提取图像 (b) MAD5变化区域提取图像

前已述及，IRMAD方法在干旱或沙漠地区以及场景中变化部分所占比例相对较大的多时相多波段遥感图像变化检测时，其性能优于MAD方法，但在变化部分所占比例相对较小的非干旱或沙漠地区，却不一定能达到很好的效果。采用IRMAD方法对

表2 典型区域提取结果对比

本文研究区检测的变量图像如图8所示。目测便可发现，其变化检测效果均逊于MAD方法和本文方法。

(a) IRMAD1(b) IRMAD2(c) IRMAD3

(d) IRMAD4(e) IRMAD5(f) IRMAD6

5 结论

1)本文方法利用第一次估计获得的初步变化结果进行分类，并对分类后的图像再次进行估计，能更为精确地剔除变化检测结果中的非变化信息，保留变化信息。通过计算2次估计结果的JM距离，获得的最终估计结果对变化信息和非变化信息具有更好的可分离性。

2)本文通过与MAD及IRMAD方法的实验结果相比较，证实了本文提出方法的有效性和优越性。

3)今后还有待于针对不同的数据源和不同的地形地貌情况做进一步的研究，以验证该方法在不同情况下的适用性。

参考文献(References)：

[1] 赵英时.遥感应用分析原理与方法[M].北京:科学出版社,2003:241-242.

Zhao Y S.The Theory and Method of Remote Sensing Application and Analysis[M].Beijing:Science Press,2003:241-242.

[2] 张继贤,杨贵军.单一时相遥感数据土地利用与覆盖变化自动检测方法[J].遥感学报,2005(3):294-299.

Zhang J X,Yang G J.Automatic land use and land cover change detection with one temporary remote sensing image[J].Journal of Remote Sensing,2005(3):294-299.

[3] Sing A.Digital change detection techniques using remotely-sensed data[J].International Journal of Remote Sensing,1989,10:989-1003.

[4] Mouat D A,Mahin G C,Lancaster J.Remote sensing techniques in the analysis of change detection[J].Geocarto International,1993,8(2):39-50.

[5] Arzandeh S,Wang J.Monitoring the change of phragmites and distribution using satellite data[J].Canadian Journal of Remote Sensing,2003,29(1):24-35.

[6] Li D R.Remotely sensed images and GIS data fusion for automatic change detection[J].International Journal of Image and Data Fusion,2010,1(1):99-108.

[7] Kwarteng A Y,Chavez P S.Change detection study of Kuwait City and environs using multi-temporal Landsat Thematic Mapper data[J].International Journal of Remote Sensing,1998,19(3):1651-1661.

[8] 张辉,王建国.一种基于主分量分析的SAR图像变化检测算法[J].电子与信息学报,2008,30(7):1728-1730.

Zhang H,Wang J G.A SAR image change detection algorithm based on principal component analysis[J].Journal of Electronics and Information Technology,2008,30(7):1728-1730.

[9] Nielsen A A,Conradsen K,Simpson J J.Multivariate alteration detection(MAD)and MAF postprocessing in multispectral,bitemporal image data:New approaches to change detection studies[J].Remote Sensing of Environment,1998,64(1):1-19.

[10]Canty M J,Nielsenb A A,Schmidtc M.Automatic radiometric normalization of multi-temporal satellite imagery[J].Remote Sensing of Environment,2004,91(3/4):441-451.

[11]Canty M J,Nielsen A A.Visualization and unsupervised classification of changes in multispectral satellite imagery[J].International Journal of Remote Sensing,2006,27(18):3961-3975.

[12]廖明生,朱攀,龚健雅.基于典型相关分析的多元变化检测[J].遥感学报,2000,4(3):197-201.

Liao M S,Zhu P,Gong J Y.Multivariate change detection based on canonical transformation[J].Journal of Remote Sensing,2000,4(3):197-201.

[13]陈垒,马润赓,申维.基于典型相关分析的遥感影像变化检测[J].地质通报,2007,26(7):916-920.

Chen L,Ma R G,Shen W.Detection of remote sensing image alteration based on canonical correlation analysis[J].Geological Bulletin of China,2007,26(7):916-920.

[14]Nielsen A A.The regularized iteratively reweighted MAD method for change detection in multi- and hyperspectral data[J].IEEE Transactions on Image Processing,2007,16(2):463-478.

[15]Canty M J,Nielsen A A.Automatic radiometric normalization of multi-temporal satellite imagery with the iteratively re-weighted MAD transformation[J].Remote Sensing of Environment,2008,112(3):1025-1036.

[16]邱卫宁,陶本藻,姚宜斌,等.测量数据处理理论与方法[M].武汉:武汉大学出版社,2008:123-145.

Qiu W N,Tao B Z,Yao Y B,et al.The Theory and Method of Surveying Data Processing[M].Wuhan:Wuhan University Press,2008:123-145.

[17]Marpu P R,Gamba P,Canty M J.Change detection using iteratively reweighted regression with neural networks[C]//2010 IEEE International Geoscience and Remote Sensing Symposium.Honolulu,HI:IEEE,2010:2563-2566.

[18]顾有林,韩帮春.快鸟数据在光学遥感器成像仿真中的应用[J].系统仿真学报,2008,20(22):6265-6267.

Gu Y L,Han B C.Application of high resolution quickbird satellite data in imaging simulation of optical remote sensor[J].Journal of System Simulation,2008,20(22):6265-6267.