朱欣然， 吴 波， 张 强

(空间数据挖掘与信息共享教育部重点实验室,福建 350000)

0 引言

土地利用/覆盖变化(land use/cover change，LUCC)是全球环境变化和可持续发展研究的核心问题之一[1-3]。LUCC不仅客观记录了地球表面景观的空间格局，也集中体现了自然因素和人类活动改变地球表面景观的时空变化过程[4]。因此，如何利用遥感技术快速、实时和具有周期性等特点，高精度动态获取LUCC及其状态是当前LUCC 研究的前沿和热点问题[5]。

LUCC自动更新的通用技术是首先利用变化检测方法获得变化/非变化的像元； 再选择非变化的像元对后期待更新影像进行分类器训练，并对后期影像中变化的像元进行重分类； 未变化像元则直接继承原始分类影像的类别，以此实现LUCC的更新过程[6]。Xian等[7-8]首次结合变化检测方法，使用2001年的LUCC分类数据自动更新出2006年的LUCC情况，该方法先使用变化矢量分析(change vector analysis，CVA)算法检测出变化/未变化像元，再利用决策树分类(decision tree classification，DTC)方法重新对变化像元进行分类。应用CVA-DTC方法自动更新LUCC时，需要2景获取时间为同一物候期的影像，且自动化更新程度较低，限制了该方法的应用； Jin等[9]提出多特征指数集成的综合变化检测方法(comprehensive change detection method，CCDM )，首先提取多种特征指数用于检测变化信息，然后结合现有的知识系统进一步更新变化/未变化像元，但该过程复杂并需要一定的知识或经验作为参考； Chen等[10]提出利用后验概率变化矢量分析(change vector analysis in posterior probability space，CVAPS)进行LUCC自动更新，该方法将原始影像的光谱空间转换为后验概率空间，实现了对另一期影像分类图的自动更新，并取得了较好效果。然而，CVAPS并没有考虑到遥感影像各波段之间和多时相之间的光谱相关性，因而可能会造成信息丢失[11]。多元变化检测(multivariate change detection，MAD)算法是另一种变化检测技术，能够将多波段的变化信息集中体现在少数波段中，从而最大限度地消除了多波段影像间的相关性，特别适用于多时相影像的变化检测[12]。因此，在前人研究基础上，本文综合发挥CVAPS与MAD的优势，提出改进CVAPS的LUCC自动更新技术，通过改善像元的变化检测精度，增强LUCC分类自动更新过程中训练样本的可靠性，提高LUCC分类自动更新的精度。为减少分类图中“椒盐”噪声的影响，本文进一步利用迭代马尔科夫随机场(iterative Markov random field，IR-MRF)模型进行分类后的空间邻域处理，以提高自动更新的精度。

1 技术方法

改进CVAPS的LUCC自动更新方法主要包括3大步骤： ①确定变化/未变化像元； ②采用迭代法选择样本进行重分类； ③采用IR-MRF 去除噪声。在检测变化时，首先利用CVAPS确定变化/未变化像元，然后利用MAD方法更新变化/未变化像元，提高检测未变化像元的精度，以便在迭代分类中高精度地选取未变化像元。考虑到像元的空间邻域相关特征，对每次的迭代分类结果进行IR-MRF后处理，去除类似“椒盐”的噪声现象。具体技术流程见图1。

图1 LUCC自动更新流程

在初始化训练样本时，将时相T1的所有像元作为训练样本输入。由于最大似然分类器(max-likelihood calssifer，MLC)能够输出像元属于不同类别的后验概率，本文以MLC分类器为例进行相关实验。从图1 的自动更新流程中可以看出，本文方法主要涉及多时相影像的变化检测、训练样本选择以及影像的IR-MRF后处理3个主要技术环节。

1.1 多元变化检测

MAD是一种以典型相关分析为基础、以方差最大化为判断准则的多元变化检测技术。该方法利用2个时相各波段之间的典型相关结构将2组多元变量变换成1组新的多元变量，并按原始变量之间的相关性大小重新分配； 同时构造出变化差异变量，最大限度地检测出变化信息。如果将2个时相的多波段遥感影像视为2组多元随机变量，则可引入MAD变化检测方法，对2组随机变量进行典型相关分析。经MAD变换后，可大大降低2景影像各波段间的相关性，增强变化信息并抑制噪声，有利于进一步的信息提取和地学分析[13]。

假设2景不同时相的影像为X和Y，波段数分别为p和q，E(X)和E(Y)为期望值，E(X)=0，E(Y)=0，p≤q。使用2组线性组合系数向量a和b，分别构造2个线性组合，即

aTX=a1X1+…+apXp，

(1)

bTY=b1Y1+…+bqYq，

(2)

其中X=[X1,…,Xp]T，Y=[Y1,…,Yq]T，这2景影像的线性离差为

d=aTX-bTY。

(3)

MAD变化检测法就是以离差d作为衡量变化信息的测度，将d的方差最大化作为判断准则，即满足max(Var{aTX-bTY})，在单位方差约束条件下Var{aTX}=Var{bTY}=1。

Var{aTX-bTY}=Var{aTX}+Var{bTY}-

2Cov{aTX,bTY}=2(1-Corr{aTX,bTY})。

(4)

MAD变换的核心就是利用式(4)找出使aTX和bTY相关性最小的系数向量a和b，而这可以借助多元统计中的典型相关分析[14-15]来求解，在此不再赘述。此时a和b即为典型相关系数，aTX和bTY为2个时相的典型变量。经典型相关分析获取的相关系数依次递减，第一对典型变量之间的相关系数最大，而最后一对的最小。只需把相关系数按从小到大重新排序，即可得到MAD差异变量，即

(5)

式中P为自由度。

MAD近似符合高斯正态分布，对MAD变量的平方标准化后的Tj则服从自由度为P的χ2分布[13]，即

(6)

(7)

ρi=Corr{aTX,bTY}，

(8)

式中σMADi为第i个MAD变量的方差。

1.2 样本自动选择

LUCC自动更新的关键是对影像中的变化像元进行重分类，前提是准确地确定影像在需更新时期内的变化与非变化像元，以便从非变化的像元中选取合适的样本用于后期影像的重分类。本文采取的策略是利用CVAPS检测出初始的变化/未变化像元，再结合MAD变化检测结果来更新未变化像元，从而确定未变化区域。由于需要从未变化的像元中选择样本，而单次变化检测的精度会直接影响样本的选择，从而影响再分类的精度。为此，本文引入迭代条件模型(iterated conditional model，ICM)自动选择未变化像元的样本，并根据2次迭代间变化/未变化像元的一致性比率确定迭代是否终止。具体的步骤如下：

1)利用T1时相的参考分类图初始化训练样本Ω=C1,C1为已知的参考分类图。

2)采用MLC分类器以初始样本Ω分别对T1和T2时相影像进行分类，同时获取2个时相的后验概率矢量P1和P2。

3)采用自动化阈值选取方法对后验概率变化矢量自动选取阈值，确定2个时相间的变化像元Θc和未变化像元Θu，这里应用直方图熵的方法自动选取阈值[16]。

4)使用IR-MRF模型更新T2时相分类图，消除变化区域中的“椒盐”现象，并更新变化像元Θcm和未变化像元Θum。

5)应用MAD算法检测T1与T2时相间的变化/未变化像元，结合步骤4中的检测结果对未变化像元进行逻辑并运算，更新2个时相间未变化像元，用于下次迭代选取样本。

6)从未变化像元中选择样本，执行步骤1—5的操作，直至2次迭代过程中变化与未变化像元趋于稳定，一致性比率大于99.9%； 最后输出T2时相的分类影像。一致性比率的定义为

(9)

式中：nunchange和nchange分别为2次迭代过程中都不变化和都变化的像元数；N为像元总数；pconsistence为2次迭代过程中的一致性比率。

1.3 IR-MRF模型后处理

MRF模型考虑了相邻像元间的空间相关特征，因而可以提高影像自动更新的精度。本文采用二阶8邻域系统N(i,j)来定义MRF模型，像元(i,j)属于类别Cl的概率由该像元的光谱信息和其邻域像元所决定[17]，即

(10)

式中：Z为归一化常数；Ucontext和Uspect分别为像元(i,j)处的邻域能量函数和光谱能量函数，即

Uspect[Cl(i,j)]=-ln{pspect[Cl(i,j)]} ，

(11)

Ucontext计算式为

(12)

(13)

式中β为邻域像元信息的权重系数，一般取值为1.6[10]。

在使用MRF模型对迭代分类结果进行后处理时，利用IR-ICM可进一步优化分类效果，大大减少分类结果中的影像噪声，提高区域连续性。其基本计算步骤如下：

1)对分类结果及每个像元的后验概率进行初始化，按照式(12)计算每个像元属于类别Cl的最大概率，即最小邻域和光谱能量。

2)按照像元所在最大概率时的类别编码更新输入的分类结果。

3)重复步骤1和2，直到2次迭代间类别一致性比率大于99.9%或者迭代次数超过10为止。

2 实验研究

2.1 研究区与实验数据

为了测试和对比所提方法的有效性，选取覆盖福建省长汀县为研究区。研究区位于E115°59′～116°39′，N25°18′～26°02′之间，东西宽66 km，南北长80 km，面积共3 090 km2。长汀县的植被覆盖率为60%，森林结构单一，林下植被稀疏； 常有大雨暴雨，降水强度大，水土流失较为严重； 主要地表覆盖类型为植被、居民地、农田、水体和裸地等。

实验用遥感数据为2013年10月26日获取的Landsat8 OLI影像(图2(a))和2003年10月4日的Landsat5 TM影像(图2(b))，成像季节一致，影像大小为2 682像元×2 210像元，空间分辨率为30 m。结合2013年野外实地考察数据、高空间分辨率遥感影像等辅助数据，对2013年Landsat8 OLI 影像进行目视解译，最终获取了该地区2013年LUCC分类图(图2(c))。

(a) 2013年Landsat8 OLI (b) 2003年Landsat5 TM (c) 2013年LUCC分类图

图2研究区原始遥感影像及参考分类图

Fig.2Originalremotesensingimagesandreferencesclassificationmapofstudyarea

2.2 分类结果自动更新

以2013年的原始遥感影像数据和对应的参考分类图为基础，可自动更新出2003年的LUCC分类结果(本文没有从2003年自动更新到2013年，主要是因为2013年有比较详细的野外调查数据，而2003年的野外调查数据很难获取)。本文选取CVAPS方法为对照，分别采用CVAPS和本文方法对2003年的LUCC进行自动更新，结果如图3所示，其中图3(a)为CVAPS的自动更新结果，图3(b)为本文方法的结果。

(a) CVAPS方法(b) 本文方法

图3LUCC自动更新的地表覆盖分类结果

Fig.3AutomaticupdatingofLUCCclassificationresults

从图3可以看出，由于长汀县植被覆盖面积较大，这2种方法对2003年的LUCC更新分类结果整体上比较接近，在图3(a)和(b)中植被和居民地的分类精度较一致； 但对于某些居民地、裸地和水体，本文方法要比CVAPS的分类精度高。

为定量评价图3中的分类结果，本文随机选择了5 043个样本用于精度验证。各类样本的选取情况如表1所示，其中包括2 309个变化像元和2 734个未变化像元。

表1 变化检测及分类精度评价样本Tab.1 Samples for precision evaluation ofchange detection and classification (个)

基于上述样本，分别计算出CVAPS(表2)与本文方法(表3)分类结果的混淆矩阵。

表2 CVAPS分类混淆矩阵Tab.2 Confusion matrix of CVAPS classification

表3 本文方法分类混淆矩阵Tab.3 Comfusion matrix of classification using method proposed in this paper

从表2和表3可以看出，CVAPS方法自动更新分类的总体精度和Kappa系数分别为76.76%和0.691，而本文方法的分类总体精度达到80.11%，Kappa系数为0.733，比CVAPS方法分别提高了近3.4%和0.04。从表3中可以看出，水体因其光谱稳定，分类精度最高，用户精度达到了98.50%； 植被和农田因二者具有相似的光谱特征而造成了一定的误分，用户精度分别为88.10%和85.67%； 由于2003年长汀县的水土流失现象较为严重，部分植被因水土流失变成极易与裸地混淆的低密度稀疏植被，因此裸地的自动更新分类精度仅仅在50%左右。对照表2和表3可以看出，本文方法对裸地的分类精度较CVAPS方法有较大改善，其分类精度由50.80%提高到56.70%。

3 结果分析与讨论

由图1示出的技术流程可知，本文方法主要涉及以下3个方面问题： ①所采用的变化检测方法不同； ②迭代样本选择不一致； ③后处理的方式不同。对上述3方面问题对LUCC自动更新的影响进一步分析。

3.1 变化检测

图4分别示出MAD方法、CVAPS法和本文方法对2个时相变化像元的检测结果(红色部分代表变化像元)。

(a) MAD方法 (b) CVAPS方法 (c) 本文方法

图43种变化检测结果对比

Fig.4Comparisonamongthreekindsofchangedetectionresults

从图4(a)可以看出，MAD方法能够检测出新建的居民地、道路和部分火烧迹地等变化。图4(b)中CVAPS方法检测的长汀县城和沿汀江两岸的变化像元(蓝色框内)明显是错误的。在长汀县城内，地物并没有发生明显的变化，而CVAPS方法却将居民地内部区域检测为变化像元； 对于沿汀江两岸，2003年与2013年的土地覆盖类型都是植被，并没有土地覆盖类型的变化。图4(c)中的检测结果则明显避免了4(b)中的错误，并且对道路和火烧迹地的变化检测也很准确，这表明在CVAPS自动更新的模型中引入MAD的变化检测技术对提高模型的自动更新精度具有重要意义。为了定量评价上述3种方法的变化检测精度，同样利用表1中的样本数据，计算各方法变化/未变化像元检测结果的精度(表4—6)。从表中可以看出，本文方法综合了CVAPS方法和MAD方法的优势，将总体精度和Kappa系数分别提高到83.90%和0.674，总体精度比CVAPS提高了近2.6%，Kappa系数提高了近0.06。表5表明，CVAPS方法的未变化像元检测精度比较高，而变化像元的检测精度较低，这可能与该方法仅仅利用后验概率迭代选择未变化像元有关。由于影像中的变化像元数量通常远小于未变化像元的检测数量，较低的变化像元的检测精度直接导致了迭代过程中的未变化像元在自动选择样本过程中逐渐累积，造成分类过程中变化/未变化像元数量的不均衡和检测精度的降低。表4表明，MAD方法中的变化/未变化像元检测精度比较均衡。因此，本文利用MAD方法在每次算法迭代之前通过逻辑运算更新了未变化像元，可以保持变化/未变化像元数量在迭代过程中的动态均衡性，因而最终检测出的变化/未变化像元的错分和漏分误差都相对较低，从而有效提取出了地表变化信息。

表4 MAD变化/未变化像元检测结果Tab.4 Results of changed/ unchangedpixels detected by MAD

表5 CVAPS变化/未变化像元检测结果Tab.5 Results of changed/ unchanged pixelsdetected by CVAPS

表6 本文方法变化/未变化像元检测结果Tab.6 Results of changed/ unchanged pixels detected bymethod proposed in this paper

3.2 迭代样本选择

在影像自动迭代分类中，每次迭代都需要更新影像中变化/未变化像元的分布以及分类结果，本文设置的迭代终止条件是邻近2次迭代过程中变化/未变化像元的一致性比率大于99.9%。由于迭代次数对变化/未变化像元以及分类结果都会造成显著影响，因而需要分析本文方法和CVPAS的迭代次数与一致性比率以及变化检测和自动分类精度评价指标之间的关系(图5)。

(a) 变化检测 (b) 分类精度 (c) 二者关系

图5迭代次数与一致性比率曲线

Fig.5Curvesofiterativetimesandconsistencyratio

图5(a)表明，本文方法通过6次迭代，而CVAPS经过4次迭代，它们的一致性比率均超过了99.9%。从图5(b)可以看出，本文方法第一次迭代时的变化检测总体精度只有74.62%，经6次迭代之后提升到83.90%； 相应地，Kappa系数也从0.466提高到0.674，相比CVAPS变化检测的总体精度和Kappa系数分别提高了近1.5%和0.05。图5(c)为迭代次数与自动更新分类结果的总体精度和Kappa系数的关系，从该图可以看出，本文方法随着迭代次数的增加，总体精度和Kappa系数也逐渐提高，总体精度从72.92%提高到80.11%，Kappa系数从0.617提高到0.733； 与CVAPS方法相比，总体精度提升了近4%，而Kappa系数也提高了约0.04。因此，本文的自动迭代选择样本技术可以较显著地提升影像自动更新过程中变化检测和分类的精度。

3.3 IR-MRF分类后处理

为了提高影像分类的空间一致性效果，减少影像分类结果中产生的“斑点”或“椒盐”噪声，本文进一步利用IR-MRF模型进行分类后的空间邻域处理，以提高自动更新的精度。图6 是采取和未采取IR-MRF模型自动更新结果的局部对比。

(a) 2003年原始影像 (b) 未采用IR-MRF分类结果 (c)采用IR-MRF分类结果

图6IR-MRF模型改进自动更新分类的结果

Fig.6ResultofautomaticupdatingandclassificationimprovedbyIR-MRFmodel

从图6(c)中可以看出，使用IR-MRF模型分类结果的局部连续性比较好，并且消除了局部区域内类别不一致的孤立点，因而整体视觉结果较好，大大降低了图6(b)中未采用后处理分类结果的“椒盐”现象。图7中的直方图定量统计了IR-MRF模型对变化检测和分类精度的影响， 从图7中可以看出，采用了IR-MRF模型方法的整体分类精度比未采用IR-MRF方法提高了近5%，Kappa系数大约提高了0.1。这表明采用IR-MRF分类后处理技术还可以有效改善影像自动更新的精度。

(a) 总体精度 (b) Kappa系数

图7采用/未采用IR-MRF模型的精度对比直方图

Fig.7HistogramforprecisioncomparisonbetweenIR-MRFmodelusedandunused

4 结论

1)本文发展了一种基于CVAPS和IR-MRF的LUCC自动更新方法，主要包括多时相影像的变化检测、训练样本选择及影像的IR-MRF后处理3个主要技术环节。

2)针对LUCC自动更新的关键(即对变化像元重分类)，本文方法综合了CVAPS与MAD变化检测的优势，同时引入ICM模型对未变化像元进行分类更新，能够有效提高变化检测结果精度。

3)基于2013年长汀县的参考分类影像对2003年的TM影像进行土地利用类型自动更新实验，实现了长汀县LUCC的自动更新。结果表明，本文方法的变化检测结果的总体精度为83.90%，Kappa系数为0.674，分别比CVAPS提高了近2.6%和0.06； 自动更新的分类结果显示，本文方法的总体精度达到80.11%，Kappa系数为0.733，比传统CVAPS的分类结果分别提高了近3.4%和0.04。

参考文献(References)：

[1] Turner II B L,Skole D L,Sanderson S,et al.Land-Use and Land-Cover Change Science/Research Plan[R].Stockholm,Geneva:[s.n.],1995.

[2] 李秀彬. 全球环境变化研究的核心领域——土地利用/土地覆被变化的国际研究动向[J].地理学报,1996,51(6):553-558.

Li X B.A review of the international researches on land use/land cover change[J].Acta Geographica Sinica,1996,51(6):553-558.

[3] 王秀兰,包玉海.土地利用动态变化研究方法探讨[J].地理科学进展,1999,18(1):81-87.

Wang X L,Bao Y H.Study on the methods of land use dynamic change research[J].Progress in Geography,1999,18(1):81-87.

[4] Turner II B L,Meyer W B,Skole D L.Global land-use/land-cover change:Towards an integrated study[J].Ambio,1994,23(1):91-95.

[5] 孙晓霞,张继贤,燕 琴,等.遥感影像变化检测方法综述及展望[J].遥感信息,2011,26(1):119-123.

Sun X X,Zhang J X,Yan Q,et al.A summary on current techniques and prospects of remote sensing change detection[J].Remote Sensing Information,2011,26(1):119-123.

[6] Chen X H,Chen J,Shi Y S,et al.An automated approach for updating land cover maps based on integrated change detection and classification methods[J].ISPRS Journal of Photogrammetry and Remote Sensing,2012,71:86-95.

[7] Xian G,Homer C,Fry J.Updating the 2001 national land cover database land cover classification to 2006 by using Landsat imagery change detection methods[J].Remote Sensing of Environment,2009,113(6):1133-1147.

[8] Xian G,Homer C.Updating the 2001 national land cover database impervious surface products to 2006 using Landsat imagery change detection methods[J].Remote Sensing of Environment,2010,114(8):1676-1686.

[9] 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.

[10] Chen J,Chen X H,Cui X H,et al.Change vector analysis in posterior probability space:A new method for land cover change detection[J].IEEE Geoscience and Remote Sensing Letters,2011,8(2):317-321.

[11] 张 路.基于多元统计分析的遥感影像变化检测方法研究[D].武汉:武汉大学,2004.

Zhang L.Change Detection in Remotely Sensed Imagery Using Multivariate Statistical Analysis[D].Wuhan: Wuhan University,2004.

[12] 盛 辉,廖明生,张 路.基于典型相关分析的变化检测中变化阈值的确定[J].遥感学报,2004,8(5):451-457.

Sheng H,Liao M S,Zhang L.Determination of threshold in change detection based on canonical correlation analysis[J].Journal of Remote Sensing,2004,8(5):451-457.

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

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,246.

[14] 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.

[15] 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.

[16] Kapur J N,Sahoo P K,Wong A K C.A new method for gray-level picture thresholding using the entropy of the histogram[J].Computer Vision,Graphics,and Image Processing,1985,29(3):273-285.

[17] Bruzzone L,Prieto D F.Automatic analysis of the difference image for unsupervised change detection[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1171-1182.

[18] 赵 鹏.新街台格庙矿区开发初期土地利用动态遥感监测[J].国土资源遥感,2015,27(4):144-149.doi:10.6046/gtzyyg.2015.04.22.

Zhao P.Land use dynamic remote sensing monitoring at the initial stage of exploitation of the Xinjie Taigemiao mining area[J].Remote Sensing for Land and Resources,2015,27(4):144-149.doi:10.6046/gtzyyg.2015.04.22.

[19] 国巧真,宁晓平,王志恒,等.地形地貌对半山区土地利用动态变化影响分析——以天津市蓟县为例[J].国土资源遥感,2015,27(1):153-159.doi:10.6046/gtzyyg.2015.01.24.

Guo Q Z,Ning X P,Wang Z H,et al.Impact analysis of landform for land use dynamic change of the partly mountainous area:A case study of Jixian County in Tianjin City[J].Remote Sensing for Land and Resources,2015,27(1):153-159.doi:10.6046/gtzyyg.2015.01.24.

[20] 李 莎,倪维平,严卫东,等.基于选权迭代估计与非监督分类的多光谱图像变化检测[J].国土资源遥感,2014,26(4):34-40.doi:10.6046/gtzyyg.2014.04.06.

Li S,Ni W P,Yan W D,et al.Change detection of multi-spectral images based on iterative estimation with weight selection and unsupervised classification[J].Remote Sensing for Land and Resources,2014,26(4):34-40.doi:10.6046/gtzyyg.2014.04.06.

[21] 杜培军,陈 宇,谭 琨.江苏滨海湿地土地利用/覆盖变化与地表温度响应遥感监测[J].国土资源遥感,2014,26(2):112-120.doi:10.6046/gtzyyg.2014.02.19.

Du P J,Chen Y,Tan K.The remote sensing monitoring of land use/cover change and land surface temperature responses over the coastal wetland in Jiangsu[J].Remote Sensing for Land and Resources,2014,26(2):112-120.doi:10.6046/gtzyyg.2014.02.19.

[22] 孙雷刚,刘剑锋,徐全洪.河北坝上地区植被覆盖变化遥感时空分析[J].国土资源遥感,2014,26(1):167-172.doi:10.6046/gtzyyg.2014.01.28.

Sun L G,Liu J F,Xu Q H.Remote Sensing based temporal and spatial analysis of vegetation cover changes in Bashang Area of Hebei Province[J].Remote Sensing for Land and Resources,2014,26(1):167-172.doi:10.6046/gtzyyg.2014.01.28.