基于TM图像的农业区域植被覆盖变化检测

王晓东，何 浩，侯 东，孙冠楠，朱文泉

(1.北京师范大学资源学院，北京 100875;2.北京师范大学地表过程与资源生态国家重点实验室，北京 100875)

基于TM图像的农业区域植被覆盖变化检测

王晓东1，2，何 浩1，2，侯 东1，2，孙冠楠1，2，朱文泉1，2

(1.北京师范大学资源学院，北京 100875;2.北京师范大学地表过程与资源生态国家重点实验室，北京 100875)

以交叉相关光谱匹配(cross correlogram spectral matching，CCSM)为基础构建土地覆盖变化强度指标，利用华北农业植被覆盖区2期不同时相的TM图像计算该地区土地覆盖变化强度图像。认为变化强度图像任意二阶邻域中像素的变化强度服从隐马尔可夫模型，用马尔可夫随机场－最大后验估计(maxium a posteriori estimation of markov random field，MRF－MAP)的方法从变化强度图像中提取植被变化区域。实验证明:该方法能够有效识别各种外源噪声造成的农业植被覆盖区域同物异谱的现象，可准确提取植被变化区域;但对于水体区域存在误判现象。

植被覆盖;变化检测;交叉相关光谱匹配(CCSM);交叉相关系数;隐马尔可夫随机场模型

0 引言

农业植被覆盖区域变化检测对于及时掌握农作物种植状况、了解区域土地利用变化具有重要意义。卫星传感器获取的遥感图像包含有大量的地物特征信息。利用同一地区、不同时期的遥感图像来确定地物状态变化的过程称为变化检测。目前已经提出的遥感影像变化检测算法主要包括分类后比较法和直接比较法［1－4］。直接比较法是通过定量评价不同时相、同一像素的光谱曲线获得不同时期地表土地覆盖变化的强度，进而分割变化强度图像提取变化区域。由于能够综合利用变化像元的强度信息和空间信息，直接比较法成为变化检测的一个重要研究方向。

直接比较法包括检测指标和图像分割两个关键技术［5］。当前的检测指标主要从不同时期地物光谱在数值、曲线形状、相似性及归属概率等特征上的差异出发进行设计，包括欧式距离［6］、夹角余弦［7］、相关系数［8］、交叉相关系数［9］、变化向量分析［6］、变化概率［10］等。目前用于图像分割［11］和描述遥感图像数据相关关系的统计模型包括有限高斯混合模型［12］、马尔可夫模型［13］、隐马尔可夫模型［14］、模糊马尔可夫模型［15］等。隐马尔可夫模型将有限高斯混合模型和马尔可夫随机场模型两者的优点相结合，是目前比较理想的差异图像先验模型。

华北农作物种植区田块面积较大、种植规模成片、种类相对单一，而TM遥感数据具有分辨率适中、单景覆盖面积大的特点，基本能够满足这种农作物种植区变化区域提取的需要。农作物种植区土地覆盖的光谱受到较大的人为干扰，农作物轮作造成的种类差异、灌溉造成的土壤背景影响以及农作物本身的物候特征，都导致了植被覆盖区较多的同物异谱现象，而TM遥感数据较高的光谱分辨率能够较好地反映地表植被覆盖的连续的光谱曲线，为准确提取变化区域提供了光谱信息保障。

本研究以2期TM遥感图像为数据源，由于华北农作物种植区的植被光谱受到各种干扰，导致各波段光谱曲线形状相同而光谱值统一偏高或偏低。针对这一特点，本文选用可比较光谱形状差异的交叉相关系数为地表变化强度的指标;利用隐马尔可夫模型来描述八邻域像素中各像素之间的影响;按照贝叶斯最小规则，根据区域变化强度灰度图像来估计类别场最大后验概率，然后自动分割图像，提取植被变化区域。

1 研究方法

1.1 变化检测指标设计

检测指标的性能是变化检测的关键。就遥感图像来说，检测指标要能够准确反映地表覆盖的微小变化，同时要抑制遥感图像获取过程中大气和传感器的噪声以及农田灌溉作业等地表农事活动对土地覆盖光谱造成的干扰。

本研究以Van der Meer［9］提出的交叉相关系数作为检测指标的设计基础。交叉相关系数反映两条光谱曲线在不同匹配位置上的相似度，其定义为

式中:a={ai|i∈K}，b={bi|i∈K}，a 和 b 为需要比较的两条光谱曲线数据，K为波段集合数;Rm为b光谱曲线移动m位置后和a光谱曲线的交叉相关系数，对b光谱曲线进行移动变换，其中m为b光谱曲线的整体移动程度，m=0表示没有移动，m=1表示b曲线整体向右移动一个波段，并将最后一个波段补在第一波段的前面;i为a光谱和移动变换后的b光谱中的波段位置;n为移动后两曲线的重叠波段数。

所有匹配位置的交叉相关系数可以组成一条差异曲线，即R={Rm|m∈K}。当两条比较光谱采用同一图像的光谱时，差异曲线反映土地覆盖完全不变情况下的标准曲线;当两条比较光谱采用不同时相的光谱时，则差异曲线反映两个时相之间的土地覆盖现状变化的实际曲线。通过比较标准曲线和实际曲线的差别，可以刻画土地覆盖变化的强度。越接近标准曲线，土地覆盖变化情况越小，反之，土地覆盖的变化程度越大。土地覆盖变化的检测指标可以通过实际曲线与标准曲线之间的均方根差(RMS)进行表达［16］，即

式中Rm，R'm分别代表标准曲线和实际曲线。显然，土地覆盖的变化必然导致地物光谱曲线的变化，进一步引起不同时相实际曲线和标准曲线的差异，并最终反映到RMS的数值上。因此，RMS综合代表了不同时相土地覆盖变化的强度。

1.2 隐马尔可夫随机场模型

差异灰度图像可以表示为一个变化强度的集合Y={y1，…，yi，…，yN|yi∈D，i∈S}，是可以观测到的随机场。其中，S={1，2，…，N}，表示所有像素的位置;D表示所有检测指标值。令L代表类别集，则L={0，1}，其中0代表没有变化，1代表有变化。图像分割就是对图像集合Y中每个像素i进行判断并标记为L中的一个类别(l)的过程，所有像素的标记结果就是类别图像。可以表示为类别的集合X={x1，…，xi，…，xN|xi∈L，i∈S}，是无法直接观测到的隐随机场，也是变化检测的结果［17］。

根据有限高斯模型，图像的概率可以用几个分量组成的混合模型来表示，即

式中:xi和yi分别代表类别图像集合X和强度图像集合Y中的元素，i∈S;θ为类别参数;p(xi=l)为混合比例参数，由类别数量决定，由灰度图像统计特征决定，没有考虑不同像素的空间属性。

马尔可夫随机场(markov random field，MRF)很好地描述了当前像素与其邻域中像素之间的相互关系［18］。当前点的邻域定义为 N={Ni，i∈S}，其中邻域不包括当前点本身i∉Ni，并且邻域关系相互对称，即i∈Nj⇔j∈Ni。MRF表达了当前像素点的灰度仅与其邻域像素的灰度有关，与其他位置像素的灰度无关，即

有限高斯模型通过马尔可夫随机场加入空间信息就可以得到隐马尔可夫随机场模型，即

式中P(xi=l|xNi)为混合比例参数，由像元在类别场X中的邻域xNi来决定。隐马尔可夫模型不仅能够描述灰度图像统计特征，还可以描述图像的空间结构，产生连续的空间聚类，因而可以去除单独的噪声，得到更准确的分类结果。

1.3 MRF－MAP分类器

基于灰度场的条件独立性，可以得到

式中θ为类别参数，且 θ=(μxi，σxi)在递归过程中获得并更新。根据Hammersley－Clifford定理，MRF分布可以等效地刻画为Gibbs随机场分布［19］，即

式中:U(X)是能量函数;Z为配分函数，即

将式(8)代入式(7)后取对数并推导，可以得到

关系式。其中，令

则类别场X的后验概率最大式(6)等价于后验势能的最小，即

1.4 检测流程

本流程的实现分为4步:①分别对两期TM图像进行辐射校正和几何纠正等预处理，获取位置准确的地表覆盖表观反射率;②逐像素计算交叉相关系数的均方根差，得到图像范围内的土地覆盖变化强度图像;③假设土地覆盖变化强度和空间关系符合隐马尔可夫模型，采用MRF－MAP分类器逐像素判断土地类型是否变化，其中变化强度图像的初始分类采用K均值算法，后验势能最小化采用条件迭代模式(iterated conditional mode，ICM)［20］，隐马尔可夫模型的隐含参数模拟采用期望最大值算法(expection－maxium，EM)［17］;④计算整体后验势能，与 MRF－MAP分类前的整体后验势能进行比较，如果2次分类前后整体后验势能未达到稳定，重新进行逐像素MRF－MAP分类，如果2次分类前后整体后验势能稳定，则停止计算，输出变化区域提取结果。整个变化检测流程如图1所示。

图1 植被覆盖区域变化检测流程Fig.1 Flow chart of vegetation cover change detection

2 应用实例

2.1 研究区概况与数据源

研究区范围在 E 115.95°～116.06°，N 39.67°～39.59°之间，位于北京市大兴区内。土地覆盖类型包括耕地、林地、草地、裸地和城镇用地。为验证本文方法的实际效果，采用2006年4月7日和2006年12月3日TM图像进行试验。图像时相选择在冬小麦覆盖区域光谱与土壤背景对比较大的抽穗和分蘖阶段，以消除冬小麦季节变化对检测结果的影响。

2.2 数据预处理和变化区域提取

对遥感图像进行辐射校正和几何纠正，去除与本研究无关的热红外波段，利用FLAASH大气校正模块，去除气溶胶和视场角的影响。将其他6个波段的DN值转化为地表反射率［21］，然后选择二次多项式进行几何纠正，几何误差控制在一个像素之内。两期TM地表反射率图像如图2所示。

图2 研究区TM地表反射率图像Fig.2 TM surface reflection images in study area

采用1.1节所述的变化检测指标逐像素进行计算，获得研究区土地覆盖变化强度图，如图3(a)所示。以此为基础，用1.3节所述的MRF－MAP分类算法进行变化区域提取，结果如图3(b)所示。

图3 变化强度图像及MRF－MAP提取结果Fig.3 Change intensity images and MRF －MAP extract results

2.3 提取结果定性对比

为了更有效地对比分析提取结果，在实验区选择了有代表性的4个典型区域，对比了典型区域的原始地表反射率、变化强度图像和变化区域MRFMAP提取结果，如表1所示。

表1 典型区域对比Tab.1 Comparison of typical region

2.3.1 检测指标性能分析

区域A前后两期土地覆盖为裸地和植被，在图像光谱上体现出较大的差异，在变化强度图上灰度值较高，说明针对土地类型不同、光谱曲线差异的土地覆盖变化，本文采用的检测指标能够准确地指示出该变化的强度。

区域C前后两期土地覆盖主要为植被，实际并未发生变化，但后时相由于进行了灌溉作业，土壤背景发生变化，从而导致了地物在图像上的色调也出现差异。在变化强度图上区域C的灰度值较区域A灰度低，与城镇区域B的灰度值相似，说明本文检测指标对于同类地物中由于噪声导致的光谱曲线差异并不敏感，能有效反映地表覆盖的变化情况。

区域D前后两期土地覆盖为水体，实际未发生变化。在变化强度图上区域D的灰度值较高，说明检测指标对于水体的覆盖变化存在较大误差。原因是本文检测指标主要对各波段按指定顺序形成的光谱曲线形状差异程度进行评价，而水体区域在各波段上反射率都很低，大气辐射的影响体现得尤为明显，况且微量的大气辐射就能导致水体光谱曲线形状发生变化，因此区域D水体形成了较高的变化强度。

2.3.2 MRF－MAP提取结果分析

区域A在后时相冬小麦种植区域内存在线状的田间道路，与前时相土地覆盖类型一致，变化强度图中线状道路灰度值低，准确体现了地表覆盖未变化的情况。MRF－MAP提取方法在判断像素变化与否的时候考虑了相邻像素的情况，由于田间道路周围都是变化的像素，所以MRF－MAP将田间道路误认为是变化区域中的噪声，并判断为变化区域。说明MRF－MAP方法在去除点状噪声的同时，没有很好地保留呈现线状特征的不变区域。

区域B是土地覆盖类型基本未发生变化的城镇区域，但该区域中存在植被绿化等现象，因此零星覆盖类型发生的变化导致了变化强度图中该区域存在较多的点状变化区。由于区域B中的变化区域周围为未变化的房屋或道路，综合邻域像素的信息后，MRF－MAP法将点状变化区域判断为未变化，与实际情况一致。

2.4 提取结果精度评价

为了进一步定量评价本文方法的有效性，在全图范围内随机选择600个点作为评价总体。通过目视解译获取评价总体的真实变化情况，并根据本方法的提取结果分别计算总体精度、总体Kappa系数、用户精度和生产者精度4个评价指标，结果如表2所示。

表2 提取结果精度评价Tab.2 Accuracy evaluation of the extract results

本文方法总体Kappa系数为0.40，说明本文方法的提取结果和真实情况接近中等吻合的程度。总体精度达到0.81，说明在以像素为基本单位的精度评价体系下，本方法的精度还存在较大的提升空间，主要错分出现在线状区域和水体区域。

3 结论

本文从比较不同时期土地覆盖变化导致的地物光谱曲线形状差异出发，以交叉相关系数为基础构建检测指标，对土地覆盖的变化强度进行表达，生成土地覆盖变化强度图像;再以隐马尔可夫模型为基础，利用变化强度和变化区域空间分布信息对变化强度图像进行图像分割，提取变化区域。本方法初步实现了变化区域的自动提取，能够避免传统方法的阈值选取和伪像素人工去除两大问题，提高了变化区域提取的效率，为快速、准确获取农作物种植区植被覆盖变化区域提供了借鉴。通过定性分析和精度评价，可以得出如下结论:

1)检测指标能够避免地面灌溉和不同时期农区植被长势差异造成的干扰，可准确反映农区土地覆盖的变化强度;

2)能够有效利用土地覆盖变化的空间信息，抑制城镇和农田中的点状噪声对提取结果的干扰;

3)整个过程无需先验知识和人工交互，实现了自动化;

4)由于本方法是以地物光谱曲线形状的差异来衡量土地覆盖变化强度，而水体本身的反射率接近于零，很小的干扰就能引起光谱曲线形状的变化，因此本文检测指标对于水体区域的检测结果误差较大。今后需在提高检测指标在水体区域的适应性或者对未变化水体区域进行掩模等方面开展进一步研究。

［1］Lu D，Mausel P，Brondizio E，et al.Change Detection Techniques［J］.International Journal of Remote Sensing，2004，25(12):2365 －2401.

［2］李淑坤，李培军，程 涛.加入多时相纹理的遥感变化检测［J］.国土资源遥感，2009(3):35－40.

［3］谢仁伟，牛 铮，孙 睿，等.基于多波段统计检验的土地利用变化检测［J］.国土资源遥感，2009(2):66－70.

［4］Foody G M.Monitoring the Magnitude of Land － cover Change Around the Southern Limits of the Sahara［J］.Photogrammetric Engineering and Remote Sensing，2001，67(7):841 －847.

［5］李月臣，陈 晋，宫 鹏，等.基于NDVI时间序列数据的土地覆盖变化检测指标设计［J］.应用基础与工程科学学报，2005，13(3):261－275.

［6］陈 晋，何春阳，史培军，等.基于变化向量分析的土地利用/覆盖变化动态监测(I)——变化阈值的确定方法［J］.遥感学报，2001，5(4):259 －266.

［7］许卫东，尹 球，匡定波.地物光谱匹配模型比较研究［J］.红外与毫米波学报，2005，24(4):296 －300.

［8］申邵洪，万幼川，龚 浩，等.遥感影像变化检测自适应阈值分割的Kriging方法［J］.武汉大学学报:信息科学版，2009，34(8):902－905.

［9］Van der Meer.Spectral Curve Shape Matching with a Continuum Removed CCSM Algorithm［J］.International Journal of Remote Sensing，2000，21(16):3179 －3185.

［10］Chen X H，Jin C，Shen M G，et al.Land － use/Land － cover Change Detection Using Change－vector Analysis in Posterior Probability Space［C］//Proc.SPIE 7144，714405，2008:10.1117/12.812671.

［11］刘永学，李满春，毛 亮.基于边缘的多光谱遥感图像分割方法［J］.遥感学报，2006，10(3):350 －356.

［12］Bilmes J A.A Gentle Tutorial of the EM Algorithm and Its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models［R］.Technical Report:ICSI.TR － 97 －02，International Computer Science Institute，Berkeley CA，USA，1998.

［13］刘伟强，陈 鸿，夏德深.基于马尔可夫随机场的遥感图像分割和描述［J］.东南大学学报:自然科学版，1999，29(11):11－15.

［14］钟家强，王润生.基于自适应参数估计的多时相遥感图像变化检测［J］.测绘学报，2005，34(4):331 －336.

［15］郑 玮，康戈文，陈武凡，等.基于模糊马尔可夫随机场的无监督遥感图像分割算法［J］.遥感学报，2008，12(2):246 －252.

［16］Wang L，Chen J，Gong P，et al.Land Cover Change Detection with a Cross－ correlogram Spectral Matching Algorithm［J］.International Journal of Remote Sensing，2009，30(12):3259 －3273.

［17］Zhang Y，Brady M，Smith S.Segmentation of Brain MR Images Through a Hidden Markov Random Field Model and the Expectation－ maximization Algorithm［J］.IEEE Transactions on Medical Image，2001，20(1):45 －57.

［18］Liu D S，Kelly M，Gong P.A Spatial－temporal Approach to Monitoring Forest Disease Spread Using Multi－temporal High Spatial Resolution Imagery［J］.Remote Sensing of Environment，2006，101(2):167－180.

［19］冯衍秋，梁 斌，陈 明，等.基于gibbs随机场的有限混合模型改进与脑部MR图像的稳健分割［J］.中国生物医学工程学报，2003，22(3):193 －198.

［20］Besag J.On the Statistical Analysis of Dirty Pictures［J］.Journal of the Royal Statistical Society:Series B(Methodological)，1986，48(3):259－302.

［21］Matthew M W，Adler Gdden S M，Berk A，et al.Atmospheric Correction ofSpectralImagery:Evaluation oftheFLAASH Algorithm with AVIRIS Data［C］//Presented at SPIE Proceeding，Algorithm and Technologies for Multispectral，Hyperspectral，and Ultraspectral Imagery IX，2003.

Vegetation Cover Change Detection in the Cropping Area Based on TM Image

WANG Xiao － dong1，2，HE Hao1，2，HOU Dong1，2，SUN Guan － nan1，2，ZHU Wen － quan1，2
(1.College of Resource Sciences＆ Technology，Beijing Normal University，Beijing 100875，China;2.State Key Laboratory of Earth Surface Processes and Resource Ecology，Beijing Normal University，Beijing 100875，China)

In this paper，a change intensity indicator of land cover based on cross correlogram spectral matching(CCSM)technique was employed to generate the change intensity image of the cropping vegetation cover area in North China between two TM images in different periods.It was first considered that the change intensity of image pixel of the two－order neighbor in the change intensity image obeyed the hidden markov random field model，and then the vegetation cover change area was extracted from the change intensity image using maximum a posteriori estimation of markov random field(MRF － MAP)model.The experiment has proved that the proposed method could precisely extract vegetation cover change and inhibit effectively the same object with different spectra due to exogenous noises in the cropping vegetation cover area.However，this method seems to perform unsatisfactorily over the water area.

vegetation cover;change detection;cross correlogram spectral matching(CCSM);cross－correlation coefficient;hidden markov random field model

TP 79

A

1001－070X(2012)02－0092－06

王晓东(1983－)，男，博士，主要从事统计遥感及GIS软件工程方面研究。E－mail:wangxd22@gmail.com。

何 浩(1980－)，男，博士研究生，主要从事资源环境遥感及3S应用技术研究。E－mail:hehaowhu@163.com。

(责任编辑:邢 宇)

10.6046/gtzyyg.2012.02.17

2011－07－28;

2011－09－09

国家高技术研究发展计划项目(编号:2006AA120101)和国家自然科学基金项目(编号:40871191)共同资助。