赵根庄,王 云,董 硕,牛志达
(1. 河北省基础地理信息中心 河北省空间地理数据工程技术研究中心,河北 石家庄 050031;2. 河北省第三测绘院,河北 石家庄050032;3. 河北师范大学,河北 石家庄 050024;4. 中建龙江建设工程集团有限公司,黑龙江 哈尔滨 150090)
植被是连接大气、水分、土壤、生物的纽带,作为生态系统存在的基础,对一个地区的生态环境状况具有一定的指示作用。植被覆盖度是指植被(包括叶、茎、枝)在地面的垂直投影面积占总面积的百分比[1]。研究植被覆盖度的时空变化,能准确地揭示地表植被的生长状况,对植被的动态变化监测及地区的生态环境评价有重要意义。
地表实测是植被覆盖度的传统测量方法,一般包括目估法、采样法、仪器法等。传统方法能够获取地区最详细、最准确的数据,但耗时耗力,只适于小范围应用。随着遥感技术的成熟,高时间分辨率和高空间分辨率的影像数据获取技术不断完善,遥感测量在大范围植被覆盖度的研究中得到广泛应用[2]。遥感测量对植被覆盖度的估算模型主要有植被指数模型、经验模型、亚像元模型以及混合光谱模型和光谱梯度差模型等[3]。近年来不少学者借用遥感技术对草原植被覆盖度做了研究:李晓松[4]等利用高光谱影像,对不同归一化差异植被指数(Normalized Difference Vegetation Index,NDVI)和最小二乘法估测荒漠化地区植被覆盖度的能力做了比较;李营[5]等基于SPOT_VGT,利用时间序列NDVI动态变化监测方法,分析了呼伦贝尔盟草原植被覆盖的时空演变特征;陈学兄[6]以SPOT4影像为数据源,利用MVC法、一元线性回归分析法和插值法,揭示了陕西省1998~2008年的植被覆盖度动态变化和空间分布规律;穆少杰[7]基于MODIS-NDVI数据反演了内蒙古2001~2010年的植被覆盖度,分析了空间格局及变化规律。
张家口坝上草原是保存完整的天然草原,是华北北部治理风沙、水土流失和建设防护林工程的重要地段。目前对坝上草原植被覆盖度的研究较少,且都是基于MODIS影像数据对NDVI的研究,把实测数据与Landsat影像数据结合起来,并研究植被覆盖度时空动态的还未见。因此本文以张家口坝上草原为研究区,通过采集2014年7月31日~8月8日的植被生长旺盛期的实测数据,结合Landsat OLI和TM遥感影像数据,采用相关分析建立遥感数据获得的NDVI值与植被覆盖度之间的相关关系模型,根据两者的关系模型和近6年Landsat遥感影像,得到坝上草原植被覆盖度空间分布图及年际变化图,并对植被覆盖度的空间分布做空间自相关分析,度量高、低植被覆盖度的聚集、离散程度,进而利用趋势分析的方法得到植被覆盖度2009~2014年变化趋势图,从而对坝上草原植被覆盖度的整体变化趋势、空间分布及其年际间的动态变化进行分析,为坝上草原草地资源的动态监测和评价提供科学依据。
坝上草原位于北京北部、内蒙古高原最南端的区域,西起张家口市的尚义,中挟康保、张北、沽源,东至承德市的丰宁、围场,位于E114°35′~116°45′,N41°00′~42°20′,平均海拔1 500 m左右,总面积18 214 km2。该地区年平均气温约1.4~4℃,无霜期80~110 d,年降水量400 mm左右[8]。由于处于华北平原和内蒙古高原交界的地方,地势骤然升高,坝上草原形成了草甸草原和干草原,干草原旱禾苗居多,草层高度为15~40 cm,草甸草原主要为多年生草本植物,高度为30~60 cm。森林植被覆盖度8%~12%,围场地带达到30%~70%。坝上草原污染极少,生态系统完整。植被的多少,直接影响到周边的生态环境,坝上地区地处农区与牧区的过渡带,东部季风气候和西北大陆性气候交汇,温带湿润气候向温带干旱气候过渡,对气候变化十分敏感。干旱寒冷、大风霜冻等自然灾害和近些年的人为破坏,使坝上地区成为中国北方生态最为脆弱的地带[9]。
野外实验工作在2014年7月31日~8月8日进行,这时植被生长最旺盛。为保证采样点选取的合理性,实验前期首先分析坝上地区土地利用现状图、植被类型图和遥感影像,实验过程中,根据地形地貌等实际情况,最终选择了52个大约250 m×250 m的样地进行实地取样,样地尽量覆盖坝上草原地区,样地分布如图1所示,并在每个样地选取3~5个直径为0.5 m的样方,样方均匀散布在样地中,共185个样方。
图1 研究区范围及样地分布Fig.1 Study area and sampling distribution
1.3.1 草原植被高光谱和草原植被覆盖度测定
本研究使用美国SOC710便携式可见/近红外高光谱成像光谱仪测定植被光谱,光谱范围400~1 000 nm,光谱分辨率4.687 5 nm。测量植被高光谱时,要求天气晴朗、无风无云,植被覆盖尽量均匀,传感器探头垂直向下,在样方中心正上方约1 m位置处,测量时正对太阳入射方向,每次测量前后都用参考板进行校正,并且用GPS在每个样方中心记录下地理坐标和海拔高度。
植被覆盖度采集的样方需与高光谱采集的样方一致。在植被高光谱和地理坐标测量完成后,用数码相机进行拍照,以计算植被覆盖度,拍照时数码相机处于样方中心上方,与地面垂直,距离植被冠层约0.5 m。与其他地表实测方法相比,数码相机测量植被覆盖度的精度最高。尤其在植被覆盖度较低的情况下,数码相机的优势更明显[10]。
1.3.2 Landsat影像的下载与处理
本文研究所用的遥感数据源于美国国家航空航天局(National Aeronautics and Space Administration,NASA)提供的Landsat-8 OLI_TIRS和Landsat4-5 TM卫星16d合成的数据,可见光波段空间分辨率为30 m,选取的数据要求云量较少,时间段在7月10日~8月20日植被的生长季。与MODIS数据相比,Landsat空间分辨率更高,对坝上草原这种小范围地区植被覆盖度的研究更加精准[11]。
下载的遥感数据产品已经过系统辐射校正和地面控制点几何校正,并且通过了高精度DEM数据的地形校正,在使用遥感数据计算NDVI时,还需在ENVI 5.0中对影像进行辐射定标、数据格式转换、大气校正等预处理。对处理好的影像进行拼接和研究区掩膜处理,得到坝上地区的地表反射率影像数据,通过波段计算求得当年的NDVI值。
在SRAnalysis预处理软件中对高光谱数据进行光谱标定、空间标定和辐亮度标定,将DN值转换为反射率,在ENVI 5.0中平滑去噪后,即可求取高光谱曲线的近红外波段和红光波段光谱反射率的平均值。由于光谱仪获得的光谱反射率数据中各个波段的范围与Landsat遥感影像不一致,因此在计算NDVI前首先在ENVI 5.0中对光谱反射率数据按照Landsat-8 OLI的波谱范围进行重采样,将近红外波段和红光波段范围分别控制在0.85~0.88 μm,0.64~0.67 μm。然后根据NDVI计算公式在波段计算工具中计算实测的NDVI,计算公式为:
式中,ρNIR为近红外波段的光谱反射率平均值,ρRED为红光波段的光谱反射率平均值。
为获得影像NDVI值,以建立实测NDVI与影像NDVI的关系模型,将采样点与Landsat影像叠加,采样点落在影像数据的某个像元中,在ArcGIS 10.0中利用点提取栅格像元值工具提取各实地测量样方对应的影像NDVI值。
在PhotoShop中切除相片的边缘部分,上、下、左、右四边分别平行切除1/5,剩余部分用于计算植被覆盖度。在ERDAS 9.2软件中,使用Modeler创建计算模型,将照片转化为灰度值,找出植被与非植被的分界点,设定阈值,将照片转化为(0,1)的二值黑白图,统计植被像元占照片总像元数的百分比,得到每个样方的实测植被覆盖度,其计算公式为[12]:
式中,PNveg为照片内植被像元数,PNsum为照片总像元数。
空间自相关是指同一变量在不同空间位置上的相关性,用于度量空间单元属性值的聚集离散程度。若某一空间单元的属性值高(低),其相邻空间单元的属性值也高(低),则为正相关;若某一空间单元与相邻空间单元的属性值一高一低,趋势相反,则为负相关[13]。空间自相关分析将研究对象的空间位置信息和属性信息综合考虑,是检验某一空间单元的属性值是否与相邻空间单元的属性值显著相关的重要指标,分为全局空间自相关分析和局部空间自相关分析,目前常用的参数主要有Moran's I指数,Geary'C指数和G系数统计量,本文中采用Moran's I指数进行分析。
2.3.1 全局空间自相关
全局空间自相关是对属性值在整个区域的空间分布特征的描述,Moran's I指数的公式如下:
式中,n表示空间单元的个数,xi和xj分别表示单元i与j的观测值,为xi的平均值,Wij为单元i与j的空间关系权重矩阵。得到Moran's I指数后需要进行显著性检验,一般采用标准化统计量z检验,z统计量公式为:
式中,E(I)和Var(I)分别为Moran's I的期望与方差。
Moran's I的值介于-1到1,大于0表示正相关,小于0为负相关,且绝对值越大空间分布的相关性越大,空间聚集性分布的现象越明显。当值趋于0时,无自相关,空间分布呈随机性[14]。
2.3.2 局部空间自相关
局部空间自相关用来分析局部空间范围内某属性值是否具有相关性,以计算空间单元与相邻空间单元的相关程度, Local Moran's I指数计算公式为:
式中对于标准状态分布S2=1,Zi和Zj是经过标准差标准化的观测值。Ii是Zi和周围点观测值平均数的积。Ii的显著性也需通过z检验加以判断。若Ii显著为正且Zi大于0,则高高集聚;若Ii显著为正且Zi小于0,则低低集聚;若Ii显著为负且Zi大于0,则高低集聚;若Ii显著为负且Zi小于0,则低高集聚[15]。
采用一元线性趋势线,分析坝上草原6年间像元尺度的植被覆盖度年际变化趋势[16]。以时间为自变量,以植被覆盖度为因变量,用最小二乘法计算出一元线性回归方程的斜率b,即每个像元的植被覆盖度年际变化值,计算公式如下:
式中,b为线性倾向值;xi为时间代表值,即2009,2010,…,2015分别为1,2,…,6;yi为像元i的植被覆盖度;-x为年份代表值的平均数;-y代表对应像元的多年平均值;n为样本数。通过时间序列和植被覆盖度序列的相关性,来分析植被覆盖度年际趋势变化的显著性。b<0说明植被覆盖度随时间变化呈减小趋势,b值越小,下降越快;b>0说明随时间变化呈增加趋势,b值越大,增长越快[17]。
草地植被覆盖度的估算模型主要有经验模型法、植被指数法、像元分解模型法和FCD模型法。经验模型法又叫回归模型法,是根据实测样方信息建立实测植被覆盖度与遥感影像数据之间的回归模型,将模型应用到整个研究区,从而计算整个研究区的植被覆盖度。植被指数法直接利用植被指数分级结果,来估算植被覆盖度。像元分解模型法中最常用的是像元二分模型,假设遥感传感器观测到的地表信息由绿色植被和无植被覆盖两部分组成,根据像元二分模型植被覆盖度估算公式计算得到植被覆盖度。本文采用两种经验模型法对坝上草原植被覆盖度进行估计,第一种方法首先对实测植被覆盖度(Fc)和实测NDVI(NDVIASD)进行回归分析,再找出NDVIASD与Landsat影像数据计算的NDVI(NDVI)之间的关系,两者结合,最终建立Fc与NDVI的回归模型;第二种方法直接采用实测植被覆盖度与Landsat影像数据计算的NDVI进行回归分析,构建模型,对比两种方法的估算精度,得到最优估算模型。采集的实验数据分为两组,一组是训练样本,样本个数为37,用于建立植被覆盖度的回归模型;另一组是测试样本,样本个数为15,进行回归模型精度检验[18]。
为验证植被覆盖度反演模型的精度,检验反演模型的模拟值与实测值之间的误差情况,本文选取相关系数(R2)、标准误差(SE)和总体预测精度(Precision)3个指标进行评价。其表达式[19-20]分别如下:
式中,n为样本个数,Xi和Yi分别为实测值和模拟值,v为残差。
3.1.1 植被覆盖度回归模型的构建
运用SPSS统计软件中的回归分析方法,选用一元线性函数、对数函数、指数函数、乘幂函数和二次多项式函数5种回归方程,来拟合实测植被覆盖度与实测NDVI、实测NDVI与Landsat影像计算得出的NDVI的相关关系。实测植被覆盖度与实测NDVI拟合结果见表1。
表1 实测植被覆盖度与实测NDVI回归方程拟合结果Tab.1 Experimental vegetation coverage and the measured NDVI regression equation fi tting results
由表1可以看出,实测植被覆盖度与实测NDVI两者具有很好的相关性,且5个模型拟合结果都能通过极显著性水平(0.01)的F检验。其中二次多项式函数拟合最好,决定系数达到0.888,一元线性函数为0.887,两者拟合结果相近,为计算简便,本文选取一次线性函数,其表达式如下:
实测NDVI与Landsat影像计算得出的NDVI拟合情况见表2。
表2 实测NDVI与Landsat影像NDVI回归方程拟合结果Tab.2 Regression equation fi tting results of NDVI measured with Landsat image
从表2中可以看出,以上5种方程对实测NDVI与Landsat影像NDVI的拟合效果都不理想,因此选用三次多项式方程拟合,相关性较好,决定系数为0.752,可以通过极显著性水平(0.01)的F检验,其回归模型的表达式如下:
将(11)式代入(10)式中,得到实测植被覆盖度与Landsat影像计算得出的NDVI的相关关系,利用两者的拟合方程可以对Landsat影像数据进行栅格计算,得到坝上草原地区各年份植被生长旺期的植被覆盖度,从而对坝上草原的植被覆盖度进行遥感动态监测。实测植被覆盖度与Landsat影像NDVI的回归模型表达式为:
通过回归模型计算出2014年Landsat影像的植被覆盖度,实测植被覆盖度与估算植被覆盖度的拟合情况如图2所示,两者相关性的决定系数R2为0.805,标准误差8.72%,可以用于植被覆盖度估算。
图2 植被覆盖度预测与实测值的回归关系(建模组)Fig.2 Regression relationship between predicted and measured values of vegetation coverage(modeling group)
3.1.2 传统植被覆盖度回归模型的构建
通过对实测植被覆盖度与Landsat影像计算得出的NDVI进行回归分析,构建传统植被覆盖度回归模型,拟合结果见表3。
表3 实测植被覆盖度与Landsat影像NDVI回归方程拟合结果Tab.3 Regression equation fi tting results of experimental vegetation coverage and the NDVI measured with Landsat image
由表3可以看出,实测植被覆盖度与Landsat影像NDVI构建的传统植被覆盖度回归模型中,二次多项式函数拟合最好,决定系数为0.642,其表达式如下:
通过回归模型计算出2014年Landsat影像的植被覆盖度,实测植被覆盖度与估算植被覆盖度的拟合情况如图3所示,两者相关性的决定系数R2为0.802,标准误差9.84%,估算精度相对较低。
图3 传统植被覆盖度预测与实测值的回归关系(建模组)Fig.3 Regression relationship between predicted and traditional measured values of vegetation coverage(modeling group)
3.2.1 植被覆盖度回归模型精度检验
通过R2、SE、Precision检验植被覆盖度回归模型的预测结果与实测植被覆盖度的误差情况及估算精度进行检验,样本检验结果表明,验证组中实测值与预测值的拟合方程决定系数R2为0.807,标准误差10.64%,估算精度为76.64%。由此可以看出,实测植被覆盖度与Landsat影像NDVI的回归模型可以较好地模拟两者的关系,对植被覆盖度估算较为准确[21]。
3.2.2 传统植被覆盖度回归模型精度检验
检验传统植被覆盖度回归模型的预测结果与实测植被覆盖度的误差情况,样本检验结果表明,验证组中实测值与预测值的拟合方程决定系数R2为0.789,标准误差11.12%,估算精度为62.69%。
植被覆盖度的实测值与两个模型预测值差异对比具体见表4。
表4 植被覆盖度实测值与模型预测值对比Tab.4 Comparison of measured values of vegetation coverage and the predicted values of the models
对比两种植被覆盖度回归模型的模拟精度,通过实测NDVI构建植被覆盖度回归模型,植被覆盖度的估算精度高于传统的植被覆盖度回归模型。
通过精度对比,选用第一种反演方法对影像NDVI进行回归模型计算,得到的坝上草原2009~2014年6年间7~8月植被生长季的平均植被覆盖度空间分布状况如图4所示。由图中可以看出,坝上草原植被覆盖度总体状况较好,从东到西呈减少趋势,坝东和坝中植被覆盖度状况明显好于坝西,坝缘优于坝中,尤以沽源、丰宁明显。围场、丰宁和沽源植被覆盖度整体处于较高水平,围场、丰宁的坝缘地区植被覆盖度最高,张北西北部较低,尚义北部和康保北部的植被覆盖度最低。
图4 2009~014年坝上草原植被覆盖度空间分布Fig.4 Spatial distribution of grassland vegetation coverage in Bashang from 2009 to 2014
坝上植被覆盖度的这种空间格局受当地的人文因素和地形地貌、水热条件等自然因素的共同影响。从气象数据来看,从西向东,坝上草原的年降水量依次递增。沽源的闪电河水库、丰宁的湿地、围场的月亮湖使得这3个地区水源丰富,草被生长茂盛。从地形来看,坝上东南高,西北低,沽源、丰宁南沿山脉,草地发育较好。康保北部是低山丘陵区,东部缓坡丘陵,植被覆盖度较低,南部为波状平原,植被覆盖度好于丘陵区。尚义、张北境内有最大的淖,尚义地处赤城深断裂带两侧,坝上为高原区,地面滩、洼、岗、丘交错,地表径流坝下比坝上多,因此尚义境内坝上地区植被覆盖度较低[22]。
通过GeoDa软件对坝上草原6年平均植被覆盖度进行全局自相关分析,计算得Moran's I指数为0.519 498,并且可以通过0.01显著水平的检验,表明坝上植被覆盖度呈显著的空间正相关,即其空间分布存在明显的聚集状态:植被覆盖度高的地区相邻,植被覆盖度低的地区相邻。
对坝上草原植被覆盖度进行局部自相关分析,LISA聚集地图和LISA显著性地图如图5a,5b所示,可以直观地反映出局部空间自相关的空间分布特征。在P≤0.05显著性水平检验的各类聚集区中,“高—高”聚集区和“低—低”聚集区占主导地位。植被覆盖度高值聚集区主要分布在沽源、丰宁的坝缘地区和围场,低值聚集区分布在坝西,主要包括康保北部、尚义北部、张北西南部。“低—高”聚集区主要在沽源、丰宁、围场的坝缘与坝中交界地带,低植被覆盖度地区被高植被覆盖度地区包围,具有空间负相关。
图5 坝上草原植被覆盖度的LISA聚集地图与LISA显著性地图Fig.5 LISA cluster map and signif i cance map of grassland vegetation coverage in Bashang
以平均植被覆盖度为指标,分析2009~2014年坝上草原植被覆盖度的多年动态特征。2009~2014年坝上草原的平均植被覆盖度年际变化曲线,如图6所示。从中可以看出,各地区和坝上全区的年际变化基本一致。由于2009年坝上地区遭遇50年一遇的大旱,年平均降水量比往年同期偏少50%~80%,植被覆盖度呈现最低值。2009~2013年,坝上植被覆盖度持续增长。2013年达到近6年最高值,为0.73,这是由于2013年气温比常年偏高0.7℃,日照接近常年,平均降水量比往年同期偏多20.4%,尤其是6月份,全市平均降水量与往年同期相比,偏多近一倍,这些都有利于植被的生长,从而提高当地的植被覆盖度。2014年入汛以来,河北省平均降水量151 mm,比往年同期减少了34%,出现了自2009年以来最严重的“卡脖旱”,坝上地区位于河北省西北部,是旱情最严重的区域,植被覆盖度下降。总体来说,坝上草原植被覆盖度在2009~2014年这6年期间,随气温、降水等条件和人为影响不断波动变化,坝上草原平均植被覆盖度2014年比2009年提高了近0.1。
图6 2009~2014年平均植被覆盖度年际变化Fig.6 Interannual variation of average of vegetation coverage from 2009 to 2014
各地区总体变化趋势一致,但各个地区之间存在着明显的空间差异,植被覆盖度由高到低依次为:围场,丰宁,沽源,张北,康保,尚义。在2009~2014年间,围场植被覆盖度最高,平均植被覆盖度为0.77;尚义最低,平均植被覆盖度为0.45,且尚义波动幅度最大,波动范围为0.20~0.64;丰宁波动幅度最小,波动范围为0.69~0.81。
2009~2014年6年间的坝上草原植被覆盖度平均值的变化趋势空间格局如图7所示。从图中可以看出,坝上全区植被覆盖度有改善的趋势。尚义、康保北部植被覆盖度改善情况最好;沽源、围场中部有轻微改善;康保西北部、张北西部植被覆盖度有所退化,尤以张北退化最严重;丰宁、围场植被覆盖度一直处于良好状态,较为稳定,但是在2009~2014年间丰宁和围场南部的植被覆盖度有轻微退化迹象。
图7 2009~2014年平均植被覆盖度变化趋势Fig.7 Trends of average of vegetation coverage from 2009 to 2014
随着气温、降水趋于正常,植被生长条件得到改善,国家推行“退耕还林还草”政策,形成了《水土保持法》《草原法》《退耕还林条例》等为主的生态建设法律体系,政府采取了一系列保护措施,极大地推动了坝上草原生态环境的建设,人们保护植被的意识也日益增强。2009~2014年平均植被覆盖度变化趋势统计结果见表4,从表4的统计结果看出,2009~2014年坝上草原的植被覆盖度得到改善。在ArcGIS 10.0中利用自然分割法对植被覆盖度的改善情况进行分级,严重退化地区面积为570.875 km2,比重占3.134%;一般退化地区面积为4 025.362 km2,比重占22.100%;无明显变化地区面积为6 907.316 km2,比重占37.923%;一般改善地区面积为4 728.686 km2,比重占25.960%;明显改善地区面积为1 981.998 km2,比重占10.882%。从整体上看,退化地区面积比重为25.23%,改善地区面积比重为36.84%,无明显变化地区面积比重为37.92%,近6年来坝上草原的植被覆盖度整体好转[23]。
表4 2009~2014年平均植被覆盖度变化趋势统计结果Tab.4 Statistical results of average vegetation coverage changes from 2009 to 2014
本文基于2014年7月31日~8月8日坝上草原植被高光谱和植被覆盖度的实测数据,结合Landsat遥感影像,建立影像NDVI和植被覆盖度的回归模型,并利用回归模型对Landsat遥感影像NDVI进行栅格计算,得到坝上草原2009~2014年6年间的植被覆盖度,对植被覆盖度的时空动态变化进行分析。Landsat遥感影像分辨率为30 m,精度较高,可以实现对坝上草原植被覆盖度的快速、精准监测。不足之处在于时间序列较短,有待进一步处理分析。通过分析,得出的结论如下:
1)利用2014年采集的实测植被覆盖度与实测NDVI构建简单的线性函数,该模型R2达到0.89;实测NDVI与Landsat影像NDVI之间三次拟合,R2为0.75。基于上述两个回归方程得到植被覆盖度的估算模型,标准误差为10.6%,预测精度达到89.80%。该模型可以用于估算坝上草原植被覆盖度。
2)坝上地区的植被覆盖度从西向东逐渐增大,丰宁、沽源的坝缘优于坝中,分布格局较为明显。植被覆盖度情况由高到低依次为围场,丰宁,沽源,张北,康保,尚义。
3)坝上草原植被覆盖度呈显著正相关,存在明显的空间聚集现象。“高-高”聚集区和“低-低”聚集区占主导地位。“高-高”聚集区主要分布在沽源、丰宁的坝缘地区和围场,“低-低”聚集区分布在坝西。
4)从各年份平均植被覆盖度可以看出,2009年由于遭遇大旱迅速降落至近6年以来的最低值,2009~2013年持续增长,2013年植被覆盖度达到近6年的最高值,2014年降水量减少,植被覆盖度再次下跌。
5)2009年以来,坝上大部分地区的植被覆盖度有改善趋势,其中改善地区面积占坝上总面积的36.84%,基本不变的占37.92%,退化的占25.23%。在空间上存在明显差异[24],沽源、丰宁植被覆盖度改善最明显,康保西北部、张北西部、尚义东部和围场植被覆盖度有所退化。