董显聪,李晓洁
(1.中国科学院东北地理与农业生态研究所,吉林 长春 130012;2.中国科学院大学,北京 100049)
植被覆盖度(fractional vegetation cover,FVC)是衡量地表植被状况的重要指标,也是区域化生态系统稳定的重要指示,在水文、生态及区域变化等方面都具有重要意义[1]。随着遥感技术的出现和发展,植被覆盖度的估算已经从地面测量方法演变为可以实现大面积和快速采集地表植被覆盖度的遥感估算方法[2]。
植被覆盖度的地面测量可以提供基于地面的测量数据,主要测量方法有目估法、采样法、仪器法和模型法[3]。植被覆盖度的遥感估算方法有回归模型法和混合像元模型法。回归模型法一般不涉及复杂的遥感机理分析和遥感物理模型,操作简单,但是这种方法有局限性,一方面是区域限制,另一方面是植被类型的限制。根据测量的植被覆盖度与植被指数的拟合关系,该方法仅适用于特定区域和植被类型,不能扩展到更大区域的植被覆盖度估计[2]。混合像元分解模型方法是在假定情况下,将遥感图像中的实际像元分解成由多个分量组成的遥感数据信息,并利用这些遥感信息构建像元分解模型来估算植被覆盖度[4]。常用的混合像元模型包括像元二分模型法、Carlson模型和Baret模型。
在利用混合像元模型对植被覆盖度进行估算研究方面,文献[5]采用像元二分模型对吉林省进行了植被覆盖度年变化估算。文献[6]运用像元二分模型对三峡库区植被覆盖度年变化进行了研究分析。文献[7]运用像元二分模型对密云水库上游植被进行了植被覆盖度遥感估算研究,模型估算精度为85%。文献[8]运用像元二分模型对整个中国北方地区植被覆盖情况进行了遥感估算和变化分析。上述研究均表明像元分解模型在地表植被覆盖度估算方面具有良好的适用性和精确性,可以适用于地表植被丰富、山地丘陵较多,以及林地和灌木丛的地表区域。在草原植被覆盖度的研究方面,文献[9]运用像元二分模型对呼伦贝尔草原地区植被覆盖度年变化进行了分析。文献[10]运用回归分析方法对河北坝上草原进行了估算分析,精度为76.64%。目前研究在针对草原地表的植被覆盖度遥感估算模型选择方面精度不一,这将对估算和监测草原的植被覆盖度精确性产生直接影响。本文旨在探寻更适于草原地区的植被覆盖度遥感估算方法,对像元二分模型、Baret模型和Carlson模型在草原地区的植被覆盖度估算精度和适用性进行比较,最后对Baret模型的参数进行了最优化处理,以进一步提高其在草原地区估算精度。
本文选择内蒙古自治区东北部大兴安岭西部的呼伦贝尔大草原为研究区(43.3°N—49.5°N,116.9°E—121.3°E),该区域海拔550~1000 m,年平均温度0℃左右,无霜期85~155 d,温带大陆性气候,属于半干旱区,年降水量250~350 mm。地表植被为典型温带草原植被,以针茅属植物为主。试验时间选择草原植被覆盖最稳定的七八月。
选择2017年7月24、26和28日与地面试验采样时间相对应的Landsat 8 OLI卫星数据作为遥感数据源。在遥感数据共享中心(http:∥ids.ceode.ac.cn/query.html)下载了10幅研究区的Landsat 8 OLI卫星影像,数据格式为TIFF格式。对影像进行如下处理:首先,对获取的影像进行辐射校正,以消除由于外界、数据获取和传输系统等因素的影响而产生的系统的、随机的辐射失真或畸变;其次,利用ENVI软件的FLAASH 大气校正模型,根据Landsat 8 OLI卫星影像的图像采集时间,并结合Landsat 8 OLI光谱响应函数对卫星影像进行大气校正,以获得每个波段的地表反射率。
2.2.1 像元二分模型
像元二分模型是线性像元分解模型方法的简化模型,也是线性像元模型中最简单和最常用的模型。它将一个像元包含的信息分为两部分:植被与非植被。即假设任一像元的信息值就是植被部分反射率Sv和非植被部分反射率Ss的线性加权和
S=Sv+Ss
(1)
在仅包含植被和非植被的像元二分模型中,假设与像元对应的地表上的植被覆盖比例为fc,与之相对应的非植被比例即为1-fc。设纯植被像元的信息贡献值为Sveg,非植被像元信息贡献值为Ssoil,则整个混合像元的信息值即为
S=fc×Sveg+(1-fc)×SSoil
(2)
通过改变公式可以得到植被覆盖度公式
fc=(S-Ssoil)/(Sveg-Ssoil)
(3)
文献[11]将每个像元的NDVI值看成是植被部分的NDVI值和非植被部分NDVI值的线性加权平均,其中植被部分的NDVI的权重就是该像元的植被覆盖度
NDVI=fNDVI∞+(1-f)NDVIs
(4)
式中,NDVIs为非植被像元的NDVI值;NDVI∞为纯植被像元的NDVI值。对公式进行变换得到植被覆盖度计算公式
fc=(NDVI-NDVI∞)/(NDVI∞-NDVIs)
(5)
从理论上讲,对于纯裸地,NDVI的值应该接近于零,但由于土壤类型、土壤粗糙度、土壤颜色和土壤含水量等因素使得NDVI值受时间和空间的影响,其变化范围一般为-0.1 ~ 0.2[12]。同理,纯植被的NDVI值也会随着植被覆盖地季节等发生变化。因此,采用一个固定的NDVI数值是不准确的。
本文利用文献[7]对密云水库流域进行植被覆盖度提取中使用的确定NDVI∞和NDVIs值的方法,取观测NDVI最大值与最小值为图像给定置信度区间内的最大值和最小值,且可以在一定程度上消除遥感图像噪声引起的误差[7]。NDVI∞和NDVIs值分别取影像中纯植被覆盖的NDVImax值0.82和纯裸土样地的NDVImin值0.15,运用上述公式估算卫星影像的植被覆盖度。
2.2.2 Carlson 植被覆盖度估算模型
文献[12]从植被、土壤和大气之间的辐射传输模型方面进行研究,获得归一化植被指数NDVI和植被覆盖度FVC之间平方关系
(6)
2.2.3Baret植被覆盖度估算模型
Baret植被覆盖度估算模型[13]的原理是建立植被覆盖度与植被垂直间隙率之间的关系。植被冠层的垂直孔隙率可以表示为叶面积指数的指数函数,文献[2]由叶面积指数函数和垂直间隙率推算出来的植被覆盖度表达式为
(7)
式中,Kp为消光系数,取决于植被结构;KVI取决于植被冠层结构、太阳天顶角和观测角和植被叶片的光学特性;VIs和VI∞分别为裸土和叶面积指数无限大时对应的植被指数值。文献[13]提出了归一化植被指数NDVI所对应的Kp/KVI的实验值为0.617 5。则植被覆盖度计算公式为
(8)
本次试验采用照相法获取地表植被覆盖度值[14],地面试验数据采集时间为2017年7月23—29日。在采集地表植被图像时,将数码相机固定在梯子上,距离地面2.5 m,保持拍摄期间数码相机镜头垂直向下,设置相机自动曝光。选择80 cm×80 cm的正方形采样框架置于采样点,拍摄得到地物照片。照片格式为JPG格式,记录可见光的R、G、B 3个波段,同时记录每个采样框架中心点的地理坐标。将采集获得的图像进行特征分析,利用改进过绿指数算法[15]计算植被和土壤的过绿特征指数(excess green index,ExG)[16-18]。文献[19]通过计算一些植物和植被阴影的对比指数值来区分植被和非植被,得出2G-R-B对比指数对植被和非植被的区分效果最好。因此,利用下式[17]计算地面植被的平均ExG
(9)
式中,Rx,y、Gx,y、Bx,y为R,G,B 3个波段在像元点(x,y)的像元值;nx和ny分别为所摄照片的行数和列数;N为照片像元的个数。利用下式计算每个像元点(x,y)的改进过绿指数值
MExGx,y=2Gx,y-Rx,y-Bx,y
(10)
式中,x、y分别为图像的行坐标和纵坐标。
图1为不同草地类型(深绿色植株和浅绿色植株)的MExG灰度直方图。
从图1可以看到,两幅图像的波峰和波谷存在明显差异,造成这种现象的原因是不同地块上土壤和草地类型存在差异,其土壤颜色和植被的绿色指数也就不同。图1(a)的植被颜色较深,呈深绿色;图1(b)的植被颜色较浅,呈浅绿色。笔者针对不同的草地类型绘制出的灰度直方图来选择不同的分类阈值。从图1可以看出,第1种草地类型MExG图的两个波峰分别位于5和40左右,波谷在20左右;第2种草地类型的两个波峰分别在10和80左右,波谷在40左右。因此,将上述影像区分植被和土壤的阈值分别定为20和40。对于一些单子叶草地类型及植株较高的草地,由于叶子结构和自身高度遮挡阳光易导致一些土壤也被误分为植被,使得MExG灰度直方图产生漏分和错分现象。因此,在照片处理过程中采用删除小面积对象方法和膨胀腐蚀操作,如图2所示。
从图2可以看到,处理之后的图像消除了细小的噪声,填充了叶片中的小空洞,并平滑了边界。相比而言,双子叶草地及植株较矮的草地类型则基本不存在这种现象。本文地面植被覆盖度获取程序的流程如图3所示。
为了比较3种植被覆盖度估算模型在草原地区的估算精度和适用性,将这3种模型估算的影像植被覆盖度估算结果与地面实测值进行比较。
首先,从总体上分析这3种植被覆盖度估算模型的估算精度。本文取50%为高低植被覆盖度分界线,覆盖度在50%以下属于低植被覆盖度,50%以上为高被覆盖度。图4为研究区植被覆盖度地面测量数据和3种模型估算值的散点图。从图中可以看出,像元二分模型估算值多位于y=x线上方,有高估植被覆盖度的现象;Carlson模型在植被覆盖度较低时,点位多位于y=x线下方,有低估植被覆盖度的现象;Baret 模型的估算值比较接近y=x线。总体上看,随着植被覆盖度的提高,点位开始分散。Baret模型估算值更加接近于地面实测植被覆盖度。
为了定量化衡量估算误差,笔者计算了3种像元分解模型估算植被覆盖度的相关系数和均方根误差,结果见表1。从表中可以看到,3种模型的相关系数相近,与地面实测的植被覆盖度都具有良好的相关性;从均方根误差RMSE来看,Baret模型的精度最高,其RMSE明显小于其他两个模型。
表1 植被覆盖度估算精度
其次,对比3种植被覆盖度估算模型在低植被覆盖度下的估算精度。图5为研究区低植被覆盖度地面测量数据和3种模型估算值的散点图。从图中可以看到,在植被覆盖度较低的情况下,像元二分模型估算值位于y=x线上方,高估了植被覆盖度;与之相反,Carlson模型则低估了植被覆盖度;Baret 模型的估算值和地面测量值最接近y=x线,对低植被覆盖度区域的估算最准确。
表2为3种像元分解模型在低植被覆盖度的情况下估算植被覆盖度的相关系数和均方根误差结果。从表中可以看到,3种模型的相关系数相近;从均方根误差RMSE来看,相比较于总体的植被覆盖度的RMSE,3种模型的RMSE都明显下降;Baret模型的精度依然最高,其RMSE明显小于另外两个模型。
表2 低植被覆盖度估算精度
最后,对比上述3种植被覆盖度估算模型在高植被覆盖度下的估算精度。图6为研究区高植被覆盖度地面测量数据和3种模型估算值的散点图。从图中可以看到,像元二分模型仍然高估植被覆盖度,Carlson模型略高估覆盖度值,相比来说,Baret模型效果较好。
表3为3种像元分解模型在高植被覆盖度的情况下估算植被覆盖度的均方根误差结果。从表中可以看到,3种模型的RMSE数值均偏高,证实了上述的结论,在植被覆盖度较低时,3种模型的准确度较高,但随着覆盖度的提高,估算精度开始下降。相比之下,Baret植被覆盖度模型估算更准确。
表3 高植被覆盖度估算精度
文献[13]在研究中选取的归一化植被指数NDVI值所对应的KP/KVI的实验值和模拟值是针对甜菜植被确定的,分别为0.617 5和0.463 1。但是由于不同的植被类型,其植被结构也会不同,叶面积指数、叶片倾角及冠层结构等植被参数也都会随着植被类型的不同和植被的生长发生改变。因此,针对草原地区不同的地表植被类型,需要确定新的KP/KVI值。
本文通过调节模型参数使Braet模型估算值的RMSE最小,来确定KP/KVI值。设置KP/KVI的取值范围为[0.5,5],取间隔0.001。图7为参数优化前后的Baret模型估算值与地面真值的拟合情况。在低植被覆盖度下,Baret模型的参数值改进为0.586 0;在高植被覆盖度下,Baret模型的参数值改进为0.756 5;将高低植被覆盖度结合计算得出的模型参数值为0.654 5。将Baret模型与参数优化后的Baret模型进行对比可见,改进后的模型更接近于y=x线,更加符合于不同的植被覆盖情况,分别提升了Baret模型在高、低植被覆盖度的精度。
表4为Baret模型参数优化前后估算的误差结果。从表中可以看到,优化后的RMSE变小,不同的植被覆盖度精度都得到了提升。
表4 Baret模型参数优化前后精度情况
本文对像元二分模型、Carlson模型、Baret模型在草原地区的植被覆盖度估算精度和适用性进行了验证和比较。以呼伦贝尔大草原为研究对象,结果表明:
(1)在草原植被覆盖度的估算中,像元二分模型方法有高估植被覆盖度的现象;Carlson模型在低植被覆盖区域低估植被覆盖度,高植被覆盖区域高估植被覆盖度;Baret模型更加接近于地面实测值,估算效果最好,整体精度为92%。
(2)对Baret模型在草原地区植被覆盖度估算的参数进行优化后,模型在不同的植被覆盖度下精度都得到进一步的提升,低覆盖度情况下的均方根误差值由0.58降低为0.57,高覆盖度情况下由0.14降低为0.13,更适用于不同的植被覆盖情况。
本文研究验证了混合像元分解模型在草原区域可以精确获取地表植被覆盖度情况,可以被用来大范围获取及监测草原植被覆盖度和年度演变,但估算精度随着植被覆盖度的提高会略有下降。