利用综合变化检测方法进行土地覆盖变化制图

2018-03-06刁娇娇龚鑫烨李明诗

自然资源遥感 2018年1期

刁娇娇，龚鑫烨，李明诗

(1.南京林业大学林学院,南京 210037； 2.南京林业大学南方林业协同创新中心,南京 210037)

0 引言

土地覆盖和土地利用是影响陆地生态系统生物多样性，水、碳和养分循环，能量平衡和温室气体排放等环境问题的重要因素，是目前学术界最关注的环境问题之一[1]。土地变化类型和分布在评估景观状态、生态变化状态以及生态变化趋势等方面得到了广泛应用[2]。根据土地利用变化的类型和变化位置的调查监测结果，可以更好地了解生态变化机制，从而模拟其对地球和生态系统的影响[3-4]。因此，持续、准确及现势性强的土地覆盖数据对自然资源和生态系统管理非常重要。但在区域和全球尺度，对景观属性长时间、连续的监测研究面临巨大的时间和资金挑战。此外，在可以检测到的土地覆盖变化中，一部分是真正的变化，另一部分是伪变化(因物候引起的变化)，因此，如何区别真正的变化和伪变化是问题的关键所在。

目前，大多土地分类检测的研究仅基于生长季节的影像，很少考虑到同一季节内[5-6]和不同季节间[7-8]的变化。基于Landsat长时间序列轨道光谱变化分析的检测方法已被证实可有效检测土地利用的变化情况[9-10]，该方法可检测出森林的变化和干扰[11]，但对其他地类变化的检测效果并不明显。另外，仅基于光谱的变化检测方法通常不能实现大区域不同土地覆盖变化的检测分析。因此，结合已有的土地覆盖和土地覆盖变化数据及光谱变化值，发展更敏感、高效、精确或者半自动的方法检测区域尺度的土地覆盖变化，乃是提高土地覆盖变化检测能力的理想方法[10，12]。

本文使用Jin等[13]提出的综合变化检测方法(comprehensive change detection method, CCDM)，根据前期土地分类数据，以及2期生长季和非生长季的Landsat影像对，计算出研究区2011—2016年间土地覆盖的最大潜在变化，并利用检测样本验证算法的精度，以期为今后大区域土地覆盖变化制图提供方法基础。

1 研究区概况及数据源

1.1 研究区概况

本文的研究区位于美国中西部的科罗拉多州，中心点坐标为W104°39′58.55″，N38°53′52.47″。科罗拉多西部为多山地带，东部为平原地区；气候晴朗干燥，但海拔差异明显，高山地区夏季凉爽，冬季则寒冷多雪，部分山峰终年被深雪覆盖；年平均降水量406 mm。夏季平均最高温度32℃，平均最低温度13℃；冬季平均最高温度9℃，平均最低温度-8℃[14]。每年4—9月为植被生长季，主要土地覆盖类型为草地(55.69%)、森林(22.51%)和灌木(12.77%)，其他土地覆盖类型共占9.03%(其中，裸地3.06%、城市2.55%、农田1.49%、湿地1.26%、荒地0.37%、水体0.21%)。大多数草地和灌木丛都是森林演替后期的植被。森林生长缓慢的原因是火干扰，植被物候随温度和降水的变化而发生显著变化。图1示出2011年Landsat影像(P033/R033)基于美国国家土地覆被数据(national land cover dataset，NLCD)的土地覆盖分类结果及研究区位置，选择研究区中2个包含典型地物的区域(区域1和区域2)为例展示基于CCDM的变化检测过程及结果图。其中区域1的面积为260 km2, 在研究时间内表现出明显的森林干扰变化；区域2面积为50 km2,表现出明显的建筑用地的变化。

图12011年Landsat影像NLCD土地覆盖分类及研究区位置

Fig.1NLCDlandcoverclassificationfromLandsatimagein2011andlocationofstudyarea

1.2 数据源及预处理

为了减少物候和季节对土地利用变化监测的影响，分别选择2011年和2016年这2个年份生长季和非生长季的Landsat遥感影像。使用的数据包括： 2期Landsat5 TM数据和2期Landsat8 OLI数据(表1)、美国国家土地覆盖数据库的NLCD 2006和NLCD 2011。以上数据均从美国地质调查局网站(http: //glovis.usgs.gov/)免费下载，全部空间数据均统一为UTM投影。

表1 Landsat遥感数据详细信息Tab.1 Detailed information of Landsat remote sensing data

利用ENVI5.2软件对TM和OLI影像进行辐射定标和大气校正。辐射定标使用每一波段的增益和偏移系数将像元DN值转换为大气上界的辐亮度，并将影像从BSQ格式转换成BIL格式；对辐射定标后的影像使用FLAASH模块进行大气校正，并裁剪出原始影像的公共区域，以保证实验数据的一致性。以2景后期(2015/2016年)影像数据为参考影像，利用伪不变特征(pseudo invariant features，PIF)法对2景前期(2010/2011年)影像数据分别进行归一化处理。

2 研究方法

2.1 CCDM法

CCDM可以检测到研究区土地覆盖的所有变化，是发展新的土地覆盖数据库的核心技术。该方法利用前期的土地覆盖数据(NLCD 2011)并结合2011—2016年间的土地覆盖变化生成后期的土地覆盖数据(NLCD 2016)，可以检测到2011—2016年间因自然或人为因素引起的土地覆盖或土地利用的变化，通过整合多期Landsat影像、土地覆盖状态信息以及土地覆盖趋势的理论，降低由物候和年际间气候变化差异导致的非真实土地分类变化的误差。

CCDM将多指标综合变化分析模型(multi-index integrated change analysis, MIICA)和分区模型(Zone)2种基于光谱的变化检测模型整合起来，通过计算差分归一化燃烧率(differenced normalized burn ratio, dNBR)[15]、差分归一化植被指数(differenced normalized difference vegetation index，dNDVI)[16-17]、变化向量(change vector, CV)和最大相对变化向量(relative change vector maximum, RCVMAX)4个光谱指标，分别提取前期和后期影像的MIICA和Zone变化图，并将二者整合。图2为CCDM技术流程图。

图2CCDM技术流程图

Fig.2TechnicalflowchartofCCDMmodel

2.2 MIICA模型

MIICA模型可计算土地覆盖已有及潜在变化，包括变化位置和变化趋势(用生物量增加(biomass increase，BI)和生物量减少(biomass decrease，BD)表示)。多光谱指标在检测不同土地覆盖变化时可以相互补充，将dNBR，dNDVI，CV和RCVMAX这4个指标整合到MIICA模型中；同时为每个指标设定阈值，将像元分为BI和BD，从而更准确地检测2期Landsat影像之间的真正变化。dNBR和dNDVI的计算用到Landsat影像的近红外波段(NIR)、第2中红外波段(SWIR2)和红光波段(R)；而CV和RCVMAX是Landsat影像各波段间相互计算得到的。4个指标的计算方法分别为

dNBR=(B5NIR-B5SWIR2) / (B5NIR+B5SWIR2)-(B8NIR-B8SWIR2) / (B8NIR+B8SWIR2)，

(1)

dNDVI=(B5NIR-B5R) / (B5NIR+B5R)-(B8NIR-B8R) / (B8NIR+B8R)，

(2)

CV=i(B5i-B8i)2，

(3)

RCVMAX=i[(B5i-B8i) /max(B5i,B8i)2] ，

(4)

式中：B5i(i=1，…，5,7)为前期影像除第6波段外的各波段；B8i(i=2，…，7)为后期影像第2—7波段； max(B5i，B8i)为B5i和B8i的最大值。

MIICA的4个指标相互补充。CV为2期影像总波谱变化的绝对值，用来表示一般的变化模式。CV值越高，土地覆盖变化的可能性越高；但CV很容易受到土地利用和土地覆盖的影响，例如农田的反射率会受季节、物候和作物轮作周期的影响，而使CV检测出来的农田变化并非土地利用的真正变化，且CV对反射值较低的水体监测能力较低。而通过RCVMAX可以检测到2期影像中全部波段的相对变化，量化相对变化的程度，表现出一般的变化模式。RCVMAX值越高，表明变化的可能性越大，例如通过RCVMAX也可以检测出水体边界的变化。但是，RCVMAX和CV都会因地形、大气和物候的影响而产生噪点，因而都不是检测干扰类型的最佳选择，也无法提供变化趋势。归一化燃烧率(normalized burn ratio，NBR)和归一化植被指数(normalized difference vegetation index，NDVI)可弥补这个不足。NBR一般用于检测火干扰[18]，NDVI则用于检测植物生长状况和活力，二者不易受到地形和辐射的影响，对森林采伐和森林火灾等更加敏感。2期NBR和NDVI的差异(即dNBR和dNDVI)同时表明了变化幅度和变化方向，正值越大表明BD的可能性越大(即在后期影像中有较少的绿色)，负值越大表明BI的可能性越大(即在后期影像中有更多的绿色)。图3为MIICA模型流程图。用MIICA模型计算dNBR，dNDVI，CV和RCVMAX的整体平均值和标准差，以标准差作为衡量单位，根据与整体平均值的光谱距离将每一个指标排序为若干个等级。本文通过设置4个指标的取值范围，得到可以反映2期影像真正变化的阈值(表2)。

图3 MIICA模型流程图Fig.3 Flowchart of MIICA model

表2 MIICA阈值确定规则Tab.2 Rule of threshold determination of MIICA

①：mean为整体平均值； ②：std为标准差。

2.3 Zone变化检测模型

Zone变化检测模型包含对森林变化的大小和方向都很敏感的dNBR和dNDVI，因此可以弥补MIICA不能捕捉森林细微变化的不足。使用Zone模型来提高由采伐、火灾或其他原因导致森林用地发生变化的森林干扰制图精度。其流程见图4。

图4 Zone模型流程Fig.4 Flowchart of Zone model

首先根据大小(标准差)和方向(正负)将dNBR和dNDVI划分为4个区间；然后将dNDVI和dNBR得到的Zone结合起来，得到1景有16个区域的影像，本文指定高阶且dNDVI和dNBR同时为正方向的区域为BI(44区域)，指定高阶且dNDVI和dNBR同时为负方向的区域为BD(33区域)；最后用逻辑法则将2种算法的结果整合，以减少误差。

2.4 变化整合

通过计算2景影像对的MIICA和Zone，得到dNBR，dNDVI，CV和RCVMAX，以期在不会导致大的错分误差和尽量减少漏分误差的前提下，进行变化整合，得到合理的变化影像。整合过程中，以NLCD 2011为基础，将土地覆盖分为“动态类”和“稳定类”2类：前者的光谱变化与物候或季节变化有关，包括森林、灌木、草地和木本湿地，2期光谱均发生变化的像元才能代表真正的变化；后者的光谱变化可能反映了真实的土地覆盖变化，包括除“动态类”之外的其他土地覆盖类型，2期光谱中任一变化即可表示真正变化。对于某些区域(如美国西部)，草本和灌木的持续性光谱信息可能因物候和季节而变化，但并非真正的变化，因此将灌木和草本归为“稳定类”。根据不同土地覆盖类型的特点，制定不同的规则。在“动态类”使用“AND”(逻辑乘)减少误差，即要求变化影像对都发生变化；在“稳定类”使用“OR”(逻辑加)，即只要有1景影像发生变化，就记为发生变化。进行“OR”运算时，变化方向先由前期变化影像决定，再由后期变化影像决定。图5是MIICA和Zone整合流程图。

图5 整合MIICA和Zone的土地覆盖变化制图流程Fig.5 Flowchart of integrating MIICA and Zone Models to map land cover change

2.5 精度评估

评价由CCDM生成的最大潜在光谱变化图，需要进行统计抽样和对抽取样本进行解译和精度评估。用NLCD 2006建立一个持续的草本和灌木的掩模，将NLCD 2006和NLCD 2011同为灌木和草本的区域定义为不变的灌木和草本，并用此区域进行变化检测验证。Landsat P033/R033影像的NLCD 2011共有25 602 375个像元点，选择影像至少1/3的像元作为验证样本(不少于8 534 125个像元)。在ENVI 5.2中，对照2个时期的影像对，利用移动窗口法[19〗，以每个像元为中心构建293像元×293像元的移动窗口(1个移动窗口记为1个样本单元)，按分层抽样法共选取100个样本单元，用于评估变化检测结果的质量。在每个样本单元中，对2011年和2016年获取的Landsat影像进行目视解译，从而确定在2011—2016年间土地覆盖类型的变化情况。一般来说，Landsat影像的获取时间不会影响土地覆盖类型的确定，因而将对影像的解译结果用作确定该样本最终土地覆盖类型的初始依据。采样和目视解译可为精度评估提供参考数据，并计算每个研究区的精度。本文进行精度评估的依据是土地覆盖变化，而非光谱变化，因此只有当实际土地覆盖类型发生变化时，才可以将其认定为正确的变化样本。

3 结果与分析

3.1 输出变化图

以研究区2个包含典型地物的区域(区域1和区域2)为例展示基于CCDM的变化检测过程及结果图。图6和图7分别为区域1(包含森林、裸地等)和区域2(包含建筑、裸地等)CCDM的中间过程图和最终结果图，包括了2对Landsat影像、对每对影像计算得到的4个指标(CV，RCVMAX，dNBR和dNDVI)、用MIICA和Zone模型得到的变化图及NLCD 2011和整合的变化图。

(e) 前期CV (f) 前期RCVMAX(g) 后期CV(h) 后期RCVMAX

(i) 前期dNBR (j) 前期dNDVI (k) 后期dNBR (l) 后期dNDVI

图6-1区域1的变化检测过程

Fig.6-1Changedetectionprocessesdemonstratedforsubset-1area

(m) 前期MIICA(n) 前期 Zone (o) 后期MIICA (p) 后期Zone

(q) NLCD2011 (r) 整合_2 MIICA (s) 整合_2 Zone(t) 整合_2 MIIC_2 Zone

图6-2区域1的变化检测过程

Fig.6-2Changedetectionprocessesdemonstratedforsubset-1area

(e) 前期CV (f) 前期RCVMAX(g) 后期CV(h) 后期RCVMAX

(i) 前期dNBR (j) 前期dNDVI (k) 后期dNBR (l) 后期dNDVI

(m) 前期MIICA(n) 前期 Zone (o) 后期MIICA (p) 后期Zone

图7-1区域2的变化检测过程

Fig.7-1Changedetectionprocessesillustratedforsubset-2area

(q) NLCD2011 (r) 整合_2 MIICA (s) 整合_2 Zone(t) 整合_2 MIIC_2 Zone

图7-2区域2的变化检测过程

Fig.7-2Changedetectionprocessesillustratedforsubset-2area

由图6可以看出，在影像2010-149与2016-070之间，区域1内发生大面积的森林采伐。对于相同区域，BD比BI更易被捕捉，因为BD的变化指数比BI更大。对于森林的生长过程，Zone模型比MIICA模型更容易捕捉斑块内的变化，所以，本文结合MIICA模型和Zone模型来减少错分误差和漏分误差。

3.2 精度分析

用分层抽样法选取的100个样本单元中，变化类32个，未变化类68个。通过目视判别，32个变化类样本单元中，有19个样本单元在2011—2016年间发生伪变化； 68个未变化类样本单元中，有3个样本单元实际发生了变化。由于只有当实际土地覆盖类型发生变化时才为真正的变化，而CCDM对细微的变化检测不够灵敏，因此对森林演替、恢复、采伐等变化的检测尚不完整。事实上，如果考虑细微变化，可认为变化类中有15个未发生变化的样本单元实际上发生了变化。本研究区内，未变化类精度达96%，变化类精度为40%，总体精度为68%(表3)；如果将15个发生微小变化的伪变化样本单元计为变化样本，则变化类精度可提高到87.5%。

表3 Landsat P033/R033影像100个验证样本的精度评价结果Tab.3 Accuracy assessment results of 100 samples for Landsat P033/R033 image

4 结论

1)本文采用CCDM方法捕捉全方位土地覆盖/利用的变化，包括可能会导致土地覆盖变化的土地利用的干扰和变化，以及微小或渐进的变化(如森林更新)。在NLCD 2006和NLCD 2011的基础上，利用CCDM得到2011—2016年间的变化，结合实地调查结果生成NLCD 2016，可以提高工作效率和增加准确率。

2)用MIICA捕捉不同种类的土地干扰可以减少漏分误差，其优点在于可以在不引入大量错分误差的同时，捕捉发生在各种生态系统中的简单、有效且完整的干扰。但MIICA为了抵消错分误差和漏分误差设置的阈值会遗漏一些微小的变化，而且对植被、物候和某些土地类型的敏感度较低，单独的MIICA分类结果除了真实变化外可能包括了很多虚假的土地覆盖变化，变化检测结果的质量也取决于选取的陆地卫星影像的质量。

3)CCDM中，将Zone和MIICA整合可以弥补MIICA的不足。Zone通过计算dNBR和dNDVI提高森林变化的检测精度，并可灵敏地检测到物候、季节或天气形态(如干旱)引起的变化。Zone可以捕捉植被发生的微小变化，可用于研究森林干扰和演替(如森林采伐、森林火灾、森林侵扰和森林更新)。

4)相同时段的2对Landsat影像相互补充，可以更好地检测土地覆盖真实的变化，同时减少因季节变化、云、雪和影像获取时间引起的虚假变化。2种模型的结合可以保留更多、更完整的干扰信息，进一步减少MIICA得到的虚假变化并且减少错分误差和漏分误差，方法简单、符合逻辑，并可生成精确的土地利用变化产品。将土地覆盖类型分为“稳定类”组(即物候或季节变化引起的非真实变化)和“动态组”(即真实变化)2组。使用“AND”或者“OR”结合变化信息降低错分和漏分误差。该方法的不足之处是对“稳定组”的监测可能遗漏一些真实变化，对“动态组”的监测可能会有虚假变化，因此在实际运用过程中，需要结合质量比较高的土地覆盖基础数据。

CCDM通过整合光谱的变化检测指标，可为国家土地覆盖更新和土地变化检测提供技术支撑。但使用遥感算法精确捕捉生态系统和大面积区域的土地覆盖信息十分复杂，在全国范围内使用仍需有效的统一算法，以平衡相互冲突的光谱轨道和不同土地覆盖模式的影响。

[1] 陈广生,田汉勤.土地利用/覆盖变化对陆地生态系统碳循环的影响[J].植物生态学报,2007,31(2):189-204.

Chen G S,Tian H Q.Land use/cover change effects on carbon cycling in terrestrial ecosystems[J].Journal of Plant Ecology,2007,31(2):189-204.

[2] Coppin P,Jonckheere I,Nackaerts K,et al.Review Article Digital change detection methods in ecosystem monitoring:A review[J].International Journal of Remote Sensing,2004,25(9):1565-1596.

[3] Homer C,Dewitz J,Fry J,et al.Completion of the 2001 national land cover database for the counterminous United States[J].Photogrammetric Engineering and Remote Sensing,2007,73(4):337-341.

[4] Xian G,Crane M.Assessments of urban growth in the Tampa Bay watershed using remote sensing data[J].Remote Sensing of Environment,2005,97(2):203-215.

[5] Pouliot D,Latifovic R,Olthof I.Trends in vegetation NDVI from 1 km AVHRR data over Canada for the period 1985—2006[J].International Journal of Remote Sensing,2009,30(1):149-168.

[6] Zhan X,Defries R,Townshend J R G,et al.The 250 m global land cover change product from the moderate resolution imaging spectroradiometer of NASA’s earth observing system[J].International Journal of Remote Sensing,2000,21(6/7):1433-1460.

[7] 黄维,黄进良,王立辉,等.基于PCA的变化向量分析法遥感影像变化检测[J].国土资源遥感,2016,28(1):22-27.doi:10.6046/gtzyyg.2016.01.04.

Huang W,Huang J L,Wang L H,et al.Remote sensing image change detection based on change vector analysis of PCA component[J].Remote Sensing for Land and Resources,2016,28(1):22-27.doi:10.6046/gtzyyg.2016.01.04.

[8] 常潇,肖鹏峰,冯学智,等.近30年长江中下游平原典型区耕地覆盖变化[J].国土资源遥感,2014,26(2):170-176.doi:10.6046/gtzyyg.2014.02.27.

Chang X,Xiao P F,Feng X Z,et al.Change of cropland in typical area of Middle-Lower Yangtze Plain over the past 30 years[J].Remote Sensing for Land and Resources,2014,26(2):170-176.doi:10.6046/gtzyyg.2014.02.27.

[9] Huang C Q,Song K,Kim S,et al.Use of a dark object concept and support vector machines to automate forest cover change analysis[J].Remote Sensing of Environment,2008,112(3):970-985.

[10] Kennedy R E,Yang Z Q,Cohen W B.Detecting trends in forest disturbance and recovery using yearly Landsat time series:1.LandTrendr-Temporal segmentation algorithms[J].Remote Sensing of Environment,2010,114(12):2897-2910.

[11] Li M S,Huang C Q,Shen W J,et al.Characterizing long-term forest disturbance history and its drivers in the Ning-Zhen Mountains,Jiangsu Province of eastern China using yearly Landsat observations(1987-2011)[J].Journal of Forestry Research,2016,27(6)1329-1341.

[12] Latifovic R,Pouliot D.Multitemporal land cover mapping for Canada: methodology and products[J].Canadian Journal of Remote Sensing,2005,31(5):347-363.

[13] Jin S M,Yang L M,Danielson P,et al.A comprehensive change detection method for updating the National Land Cover Database to circa 2011[J].Remote Sensing of Environment,2013,132:159-175.

[14] Weatherbase.Denver,Colorado Travel Weather Averages[EB/OL].[2013-07-10].http://www.weatherbase.com/weather/weather.php3?s=096427&refer=&cityname=Denver-Colorado-United-States-of-America.

[15] Key C H,Benson N C.Measuring and remote sensing of burn severity[C]//Neuenschwander L F,Ryan K C,Gollberg G E.Proceedings of Joint Fire Science Conference and Workshop:Crossing the Millennium:Integrating Spatial Technologies and Ecological Principles for A New Age in Fire Management.Boise,Idaho:University of Idaho and the International Association of Wildland Fire,1999:284.

[16] 魏新彩,王新生,刘海,等.HJ卫星图像水稻种植面积的识别分析[J].地球信息科学学报,2012,14(3):382-388.

Wei X C,Wang X S,Liu H,et al.Extraction of paddy rice coverage based on the HJ satellite data[J].Journal of Geo-Information Science,2012,14(3):382-388.

[17] 谭柳霞,曾永年,郑忠.林火烈度遥感评估指数适应性分析[J].国土资源遥感,2016,28(2):84-90.doi:10.6046/gtzyyg.2016.02.14.

Tan L X,Zeng Y N,Zheng Z.An adaptability analysis of remote sensing indices in evaluating fire severity[J].Remote Sensing for Land and Resources,2016,28(2):84-90.doi:10.6046/gtzyyg.2016.02.14.

[18] López García M J,Caselles V.Mapping burns and natural reforestation using Thematic Mapper data[J].Geocarto International,1991,6(1):31-37.

[19] 吴俊,孟庆岩,占玉林,等.一种基于移动窗口的城市绿度遥感度量方法[J].地球信息科学学报,2016,18(4):544-552.

Wu J,Meng Q Y,Zhan Y L,et al.A measure of urban green index in urban areas based on moving window method[J].Journal of Geo-information Science,2016,18(4):544-552.