基于时序遥感的喀斯特山区植被覆盖研究

黎良财 ，陆灯盛 ，张晓丽 ，邓 利 ，张 莹

（1.北京林业大学 林学院，北京 100083； 2.广西生态工程职业技术学院，广西 柳州 545004；3.密歇根州立大学 全球变化与对地观测中心，密歇根 东兰辛 48823）

基于时序遥感的喀斯特山区植被覆盖研究

黎良财1，2，3，陆灯盛3，张晓丽1，邓 利2，张 莹1

（1.北京林业大学 林学院，北京 100083； 2.广西生态工程职业技术学院，广西 柳州 545004；3.密歇根州立大学 全球变化与对地观测中心，密歇根 东兰辛 48823）

为监测喀斯特山区植被覆盖度演化的过程与规律，以广西壮族自治区柳江县喀斯特山区1987年、1990年、1995年、2000年、2005年、2010年和2013年Landsat TM/ETM+/OLI时间序列遥感影像为数据源，经辐射定标、大气校正等预处理后提取归一化植被指数(Ndvi)，运用混合像元分解模型由Ndvi值反演研究区各时相植被覆盖度，再利用决策树分类法对植被覆盖度分级并计算各时期植被覆盖度的分布特征、变化趋势、转移矩阵及植被覆盖度与海拔、坡度的关系。结果表明：研究区平均植被覆盖度由1987年的68.61%下降到1990年的61.46%再上升到2013年的66.92%；1987—2013年植被覆盖度的变化服从双峰分布，主要集中在-0.25～+0.25区间；植被覆盖度的分布与海拔、坡度相关，呈现出植被覆盖度随海拔而升高，同时也随坡度而升高的现象。该方法可以准确、定量地反演植被覆盖情况，可以为喀斯特地区长期植被动态监测提供参考。

植被覆盖度；多时相遥感；像元分解模型；决策树分类；喀斯特山区

喀斯特(Karst)是指由高溶解度的岩石和充分发育的次生孔隙度相结合并产生特殊的水文和地形的地区[1-2]。在全球三大喀斯特地形集中分布区中，中国西南地区是面积最大、喀斯特发育最强烈、景观类型最多、生态环境最复杂、人地矛盾最尖锐的地区[3]。广西位于中国西南喀斯特集中分布带，全区109个县(市、区)中有75个县份有岩溶出露，根据2005年对广西喀斯特地区的遥感调查结果，广西有岩溶土地面积8.33万km2，占全区土地总面积的35.2%，是岩溶面积最多的省（区）之一。喀斯特地区具有土层薄、植被覆盖低、水土流失严重、自然生态环境脆弱等特征[4]，加上人类社会经济活动的影响，极易造成土地石漠化现象。

植被作为生态系统的重要组成部分，是连接土壤，大气和水分的自然“纽带”，它具有涵养水源、减缓径流、保持水土等重要作用；同时也是人民群众农、林、牧业生产的重要载体，成为岩溶地区可持续发展、实现小康社会的重要保障。目前，利用遥感技术对植被覆盖进行监测与评价方面，国内外已经对林区、沙区、旱区、矿区进行了许多的研究与试验[5-10]，但针对喀斯特山区的定量研究较少[11]。本文以全国岩溶石漠化综合治理重点县广西柳江县为研究区，基于多时序陆地卫星影像对近26年来地表植被覆盖进行变化监测，分析了植被覆盖的时空变化规律，为岩溶山区的植被覆盖调查评价、生态保护和综合治理等提供基本数据和理论依据。

1 研究区概况

广西壮族自治区柳江县地处桂中岩溶盆地地区，109°58′～ 110°45′E，24°21′～ 24°26′N（见图1A）。地势由东、西南两面向中部倾斜；地貌类型以岩溶丘陵为主，平原、丘陵、洼地、谷地、孤峰、残丘交错分布，岩溶漏斗及落水洞极其发育（见图1C）。气候属南亚热带季风气候区，常年温暖潮湿，雨热充沛，年均降水量1 482.9 mm，特别是夏季的持续降雨容易造成水土流失。岩溶地区居住有瑶、壮等少数民族，各民族有多彩的文化背景，部分保留有刀耕火种、烧山放牧的习惯；随着区域人口增长和经济发展，人地矛盾不断加剧，特别是开垦采矿等人为活动的影响，致使有限的石山植被遭受破坏，加速了石漠化的形成和自然灾害的发生，严重制约了柳江县经济社会的可持续发展，阻碍了全面建设小康社会目标的实现。近年来，随着人们的环保意识增强和政府保护力度加大，许多石漠化治理的政策和措施在该县执行。

图1 广西柳江县位置、卫星影像和数字高程模型Fig.1 Location map, satellite image and DEM of Liujiang county

2 材料与方法

2.1 数据源选择

本研究采用的遥感影像为Landsat5/7/8卫星的TM/ETM+/OLI数据Level 1T（L1T）级别产品，L1T数据经过系统辐射校正和地面控制点几何校正，并且通过数字高程模型（DEM）进行了标准地形校正[9]，采用统一横轴墨卡托（UTM）地图投影、世界大地系统（WGS）84坐标，空间分辨率30 m。遥感影像获取时间分别为1987年10月26日（TM）、1990年 9月 16日（TM）、1995年11月17日（TM）、2000年11月6日（ETM+）、2005年10月11日（TM）、2010年11月10日（TM）和2013年12月4日（OLI）；辅助数据有柳江县行政区划矢量图层（Shape fi le格式）和ASTER GDEM 30米分辨率高程模型（IMG格式）。

2.2 数据处理与信息提取

2.2.1 影像裁剪

利用柳江县行政区划矢量图作为掩膜，经过与遥感图像投影坐标统一后，对1987-2013年7个时期影像及数字高程模型进行裁剪，得到研究区范围的影像数据，面积为253 723 hm2，见图1B、1C。

2.2.2 辐射定标

常用遥感影像的像元值是经过量化的、无量纲的DN值，利用DN值只能进行影像内部的相对比较，而对多时相、多区域、多种传感器的遥感数据进行定量分析时，需要将DN值转换成对应的辐射亮度值，这个计算过程即为辐射定标，是保持数据一致性和数据精度的一项重要的基础工作[9]。由于L1级产品的DN值是由辐射亮度经线性变换得到的，因此从L1级产品计算辐射亮度只需利用相关参数进行线性反变换即可，见公式（1）[9,12-13]：

式（1）中，Lλ：辐射亮度；DN：像元灰度值；Gain：增益；Bias：偏置。

2.2.3 大气校正

卫星在接受讯号的过程中受到大气的影响，致使接受的讯号失真，为了获得真实的地表反射率，必须去除大气对讯号的影响[9,14]。常见的大气校正方法有基于辐射传输模型的MODTRAN模型、LOWTRAN模型、ATCOR模型和6S模型等；基于简化辐射传输模型的黑暗像元法；基于统计学模型的反射率反演等[9]。本研究使用了MODTRAN4+辐射传输模型，可对多光谱、高光谱数据及航空影像进行快速大气校正，能有效消除大气和光照等因素对地物反射的影响，获得地物较为准确的反射率和辐射率、地表温度等[9,15]。

2.2.4 Ndvi提取

植被指数是利用叶冠的光学参数提取的独特的光谱信号，是遥感领域中用来表征地表植被覆盖、生长状况的一个简单又有效的度量参数[9,16]。绿色植物在近红外波段具有高反射而在红外波段具有强吸收，在利用多光谱图像提取植被信息时多采用红外/红波段图像作比值，比值处理以后可以消除部分辐射误差[9]。本研究采用归一化植被指数（Normalized Difference Vegetation Index，Ndvi）模型，计算公式如下：

式（2）中，NIR为近红外波段反射率，R为红色波段反射率。Ndvi取值介于（-1，1）之间，正值表示有植被覆盖，且随覆盖度增大而增大[9]。

2.2.5 植被覆盖度反演

植被覆盖度是乔木、灌木和草本植物地上部分（枝、叶、茎）垂直投影的面积占地面的比率。遥感测定植被覆盖度的方法大致有3种：经验模型法、植被指数法和像元分解模型法[17]。其中，像元分解模型法不依赖实测数据，可以广泛使用；像元分解模型中又以像元二分模型最常见，它假定一个像元的信息可以表达为由绿色植被和无植被覆盖（裸土）两部分组成，像元二分模型估算植被覆盖度时多采用Ndvi数据，此时植被覆盖度Fc的计算公式为[18]：

式（3）中， [Ndvi]为影像中各像元的Ndvi值；[Ndvi]_soil为全裸土覆盖像元的Ndvi值；[Ndvi]_veg为纯植被覆盖像元的Ndvi值。

对于大多数类型的裸地表面，理论[Ndvi]_soil值应该接近0，然而由于大气、地表温湿度、粗糙度、土壤类型、土壤颜色等条件的影响， [Ndvi]_soil会随着时间和空间而变化[9,19]。[Ndvi]_veg代表全植被覆盖像元的最大值，由于植被类型的不同，季节的变化，叶冠背景的污染等因素，[Ndvi]_veg也会随时间和空间而改变[9,20]。因此，对于各个时期的遥感影像不宜采用固定的[Ndvi]_soil和[Ndvi]_veg值。本文对每期Ndvi图像进行统计，以1%和99%为置信度，计算得各年度[Ndvi]_soil和[Ndvi]_veg阈值见表1。将[Ndvi]_soil和[Ndvi]_veg值代入（3）式，即可反演出1987—2013年广西柳江县植被覆盖度。

表1 纯裸土和纯植被像元的Ndvi阈值Table 1 Threshold value in Ndvi image for pure soil and pure vegetation pixel

2.2.6 植被覆盖度分级

参照岩溶地区石漠化分级标准中关于植被覆盖度的划分，将研究区植被覆盖度分为如下级别，即低覆盖度（LVC，VC＜30%，对应重度石漠化），中低覆盖度（LMVC，30%≤VC＜50%，对应中度石漠化），中高覆盖度（MHVC，50%≤VC＜70%，对应轻度石漠化），高覆盖度（HVC，70%≤VC＜90%，对应潜在石漠化），全覆盖度（FVC， VC≥90%，对应无石漠化）。在建立分级的实际工作中，由于全覆盖度（FVC）阈值区间较窄可合并入高覆盖度（HVC），因此可将研究区植被覆盖度分为LVC、LMVC、MHVC和HVC（70%≤VC≤100%）4个等级。基于以上分级阈值，构造简单的决策树算法对研究区植被覆盖度进行分级，结果见图2。

图2 植被覆盖度分级决策树算法及分级结果Fig.2 Decision tree algorithm and result for vegetation coverage images classi fi cation

3 结果与分析

3.1 植被覆盖度特征分析

在ENVI软件支持下对研究区1987—2013年植被覆盖度、植被覆盖度等级图进行统计，可以得出各时期植被覆盖度分级分布与平均分布情况，结果见表2和图3。研究区植被覆盖度大致经历了三个演变阶段：（1）1987年到1990年大幅下降时期，植被覆盖度从68.61%下降到61.46%，减少7.15%；（2）1990年到2005年缓慢上升时期，植被覆盖度从61.46%持续上升到67.80%，增加6.34%；（3）2005年到2013年的平稳波动期，植被覆盖度先降低2.05%再上升1.17%，总体减少0.88%。

表2 1987—2013年广西柳江县植被覆盖度分级统计Table 2 Statistics table of different vegetation coverage in Liujiang county from 1987 to 2013 hm2

对研究区植被覆盖度等级的面积和比例分析可知，从1987到2013年各个时期均以高植被覆盖度（HVC）占优势，达到总面积的40%～53%左右；且植被覆盖度等级表现出相同的趋势，即高植被覆盖度（HVC）＞中高植被覆盖（MHVC）＞中低植被覆盖度（LMVC）＞低植被覆盖度（LVC）；各植被覆盖度等级之间的动态变化则引起研究区平均植被覆盖度的变化。

3.2 植被覆盖度与地形关系分析

图3 1987—2013年广西柳江县植被覆盖度等级分布及变化趋势Fig.3 Variable trend and grade distribution of vegetation coverage from 1987 to 2013

为了揭示植被覆盖度与地形的关系，在ArcGIS软件空间分析工具支持下，对研究区数字高程模型提取海拔、坡度图层，再将海拔、坡度图层进行分级；海拔分级根据研究区高程区间，按照自然裂点法分为5级；坡度分级参考《水土保持综合治理规划通则GB_T15772—1995》之标准，对研究区划分为4个坡度级：微坡I（小于5°），缓坡 II（5°～ 15°），较陡坡 III（15°～25°），陡坡IV（大于25°）；然后将海拔、坡度级图层和植被覆盖度图层进行叠置分析与统计，统计结果见表3、表4。两个统计表揭示了一个共同的规律：从1987年到2013年，广西柳江县植被覆盖度随着海拔、坡度的上升也随之升高；这与贵州赤水河流域[21]、黄土高原文安驿流域[22]的同类研究结果不同。主要是因为广西柳江县地处亚热带南缘，水热丰富，整体海拔较低，非常适宜植被生长，而人类的生产经营活动更多在低海拔、缓坡度地区进行，说明人类活动对研究区植被覆盖度起决定性作用。

表3 不同海拔高的植被覆盖度统计Table 3 Average vegetation coverage in different elevation classes %

3.3 植被覆盖度变化监测与分析

为了分析广西柳江县植被覆盖度的动态变化，探索变化规律，通过变化检测工具对1987—2013年以来7个时期的植被覆盖度图做初期末期差值运算，将变化比例按照表5的标准进行分级，并统计各级变化比例，如表5所示。

表4 不同坡度级的植被覆盖度统计Table 4 Average vegetation coverage in different slope classes %

分析表5可以得出图4的植被覆盖度变化趋势，即从1987—2013年除了河流等地物外其余地表的植被覆盖度都处于动态变化中，变化趋势服从双峰分布，变化主要集中在Change（+1）和Change（-1）区间，占柳江县总面积的71.66%～88.92%；而变化幅度越大，变化的面积比例越小。

3.4 植被覆盖度转移矩阵分析

为了明确广西柳江县1987—2013年间各级植被覆盖度的相互转移面积，利用变化检测工具将1987—2013年植被覆盖等级图导入统计，可以得到植被覆盖度等级相互转化的转移面积矩阵数量关系，见表6。

从1987年到2013年，呈现出高覆盖度（HVC）向其他低覆盖度（LVC，LMVC，HMVC）转移的趋势，净转移面积5 929 hm2，其中：低覆盖度（LVC）增加3 569 hm2，中低覆盖度（LMVC）增加1 338 hm2，中高覆盖度（HMVC）增加1 022 hm2，故研究区整体植被覆盖度降低。

表5 1987—2013年植被覆盖度变化统计Table 5 Statistics table for Vegetation coverage change from 1987 to 2013 %

图4 1987—2013年植被覆盖度前后期变化等级分布Fig.4 Change level distribution of vegetation coverage from 1987 to 2013

表6 1987—1990年研究区植被覆盖度等级转移面积矩阵Table 6 Transfer area matrix for vegetation coverage from 1987 to 1990 单位：hm2

4 结论与讨论

通过对中国南方广西柳江县喀斯特山区1987—2013年TM/ETM+/OLI遥感图像进行辐射定标、大气校正、植被指数提取等预处理，然后以像元分解模型反演地表植被覆盖度，并采用决策树分类算法对植被覆盖度进行分级分类，最后在ENVI和ArcGIS软件的变化检测及空间分析功能支持下，探讨了研究区植被覆盖度的分布特征、变化趋势、转移矩阵及植被覆盖度与海拔、坡度的关系。研究结论表明：

（1）研究区植被覆盖度大致经历了从1987年到1990年大幅下降时期、1990年到2005年缓慢上升时期和2005年到2013年的平稳波动期三个演变阶段。（2）通过分析海拔、坡度与植被覆盖度的关系，得出从1987年到2013年研究区植被覆盖度随海拔而升高，同时也随坡度而升高的现象。（3）从1987年到2013年各个时间段植被覆盖度处于动态变化中，变化趋势服从双峰分布，变化主要集中在Change（+1）和Change（-1）区间。（4）从植被覆盖度等级互相转移的方向来看，从1987年到1990年，植被覆盖度由高向低转移的面积较多；从1990年到1995年、2000年、2005年，植被覆盖度由低向高转移的面积较多；从2005年到2013年，植被覆盖度先由高向低转移而后由低向高转移。

分析1987—2013年26年间研究区政治、经济和政策的发展变化阶段，可以发现对应于植被覆盖度变化的一些驱动力：（1）1987—1990年处于改革开放前期，该时期对资源能源的需求大增以及环保意识薄弱使植被遭受滥砍乱伐、毁林开垦、采石挖矿等人类活动破坏，植被覆盖呈加速退化趋势，特别是在海拔低、坡度平缓的地区受影响更严重。（2）1990年广西全区开始“造林灭荒、绿化达标”工程，1998年国家启动退耕还林、公益林保护与补偿等工程，2010年启动石漠化治理和农村沼气建设等工程，一方面保护增加了地表的植被覆盖，另一方面减少了山区群众对薪柴能源的需求，使植被覆盖度稳步增加并维持在高位波动。喀斯特山区的生态环境非常脆弱，需要更加重视生态文明建设，使经济社会能够可持续发展。

致谢：文中部分数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站（http://www.gscloud.cn）。

[1]D C Ford, P W Williams. Karst Geomorphology and Hydrology[M]. London: Unwin Hyman, 1989: 547-582

[2]D.C.福特, P.威廉姆斯, 宋林华.喀斯特的定义及其发展[J].地理译报, 1990(4): 6-11.

[3]王德炉, 朱守谦, 黄宝龙. 石漠化的概念及其内涵[J]. 南京林业大学学报(自然科学版), 2004, 28(6): 87-90.

[4]王敬贵, 杨德生, 余顺超, 等. 珠江上游喀斯特地区土地石漠化现状遥感分析[J]. 中国水土保持科学, 2007, 5(3): 1-6.

[5]J. Duncan, D. Stow, J. Franklin, et al. Assessing the relationship between spectral vegetation indices and shrub cover in the Jornada Basin, New Mexico[J]. International Journal of Remote Sensing, 1993,14(18): 3395-3416.

[6]Larsson H. Linear regressions for canopy cover estimation in Acacia woodlands using Landsat-TM,-MSS and SPOT HRV XS data[J]. International Journal of Remote Sensing, 1993, 14(11):2129-2136.

[7]Luciano Telesca, Rosa Lasaponara. Quantifying intra-annual persistent behaviour in SPOT-VEGETATION NDVI data for Mediterranean ecosystems of Southern Italy[J]. Remote Sensing of Environment, 2006, 101 (1): 95-103.

[8]Lasaponara R. On the use of principal component analysis(PCA)for evaluating inter-annual vegetation anomalies from SPOT/NDVI temporal series[J]. Ecological Modelling, 2006, 194(4):429-434.

[9]黎良财, 邓 利, 曹 颖, 等. 基于NDVI像元二分模型的矿区植被覆盖动态监测[J]. 中南林业科技大学学报, 2012,32(6): 18-23.

[10]游浩妍, 骆成凤, 刘正军, 等. 基于MODIS植被指数估算青海湖流域植被覆盖度研究[J]. 遥感信息, 2012, 27(5): 55-60.

[11]易 胜. 基于RS和GIS岩溶地区植被覆盖度分析——以湖南洛塔为例[D]. 桂林: 广西师范大学, 2005.

[12]韦玉春, 黄家柱. Landsat5图像的增益、偏置取值及其对行星反射率计算分析[J]. 地球信息科学, 2006, 8(1): 110-113.

[13]Zhang Zhaoming, He Guojin. Radiometric calibration of Landsat 5 TM Data [J]. Review Science & Technology, 2008,26(7):54-58.

[14]茅荣正, 倪绍祥, 蒋建军. LANDSAT-7 ETM+影像大气校正算法IDL的实现[J]. 测绘通报, 2004(1): 8-10.

[15]邓书斌. ENVI遥感图像处理方法[M]. 北京: 科学出版社,2010.

[16]Guo Ni. Vegetation Index and its advances[J]. Arid meteorology,2003,21(4):71-75.

[17]Purevdorj Ts, Tateishi R, Ishiyama T, et al. Relationship between percent vegetation cover and vegetation indices[J]. International Journal of Rote Sensing, 1998,19(18): 3519-3535.

[18]程红芳, 章文波, 陈 锋. 植被覆盖度遥感估算方法研究进展[J].国土资源遥感, 2008(1): 13-18.

[19]Toby N, Carlson, David A, et al. On the relation between NDVI,fractional vegetation cover and leaf area Index[J]. Remote Sensing of Environment, 1997, 62(3): 241-252.

[20]Kaufman Y J, Tanre D. Atmospherically resistant vegetation index (ARVI) for EOS-MODIS[J]. IEEE Transactions on Geoscience and Remote Sensing, 1992, 30(2): 261-270.

[21]蔡 宏, 何政伟, 安艳玲, 等. 基于RS和GIS的赤水河流域植被覆盖度与各地形因子的相关强度研究[J]. 地球与环境,2014, 42(4): 518-524.

[22]程圣东, 李占斌, 鲁克新, 等. 文安驿流域植被覆盖度时空分异及其与地貌因子关系研究[J]. 西安理工大学学报, 2011,27(2): 145-150.

Study on the dynamics of vegetation coverage in the karst mountainous area based on multi-temporal remote sensing images

LI Liangcai1,2,3, LU Dengsheng3, ZHANG Xiaoli1, DENG Li2, ZHANG Ying1
(1.College of Forestry, Beijing Forestry University, Beijing 100083, China; 2.Guangxi Ecological Engineering Vocational and Technical College, Liuzhou 545004, Guangxi, China; 3.Center for Globe Change and Earth Observations, Michigan State University, East Lansing,MI 48823, USA)

This article explores the distribution characteristics and change tendency of vegetation coverage in the karst mountains area of Southwestern China, taking the Liujiang County of the Guangxi Zhuang Autonomous Region as an example. Landsat TM&ETM+&OLI time sequence images from 1987 to 2013 were used as the data source, and radiation calibration and atmospheric correction carried out.Then normalized difference vegetation index (NDVI) was extracted, and vegetation coverage inversed using mixed pixels decomposition model. And then vegetation coverage was divided into different grades by using the decision tree classi fi cation algorithm. Results showed that the vegetation coverage in overall study area dropped from 68.61% in 1987 to 61.46% in 1990, and then increased to 66.92% in 2013. Vegetation coverage changed between -0.25～+0.25 and obeyed a bimodal distribution from 1987 to 2013. It increased with altitude and slope. This method could accurately and quantitatively inverse vegetation coverage dynamics and can be used as reference for long-term dynamic monitoring of vegetation in karst areas.

vegetation coverage; multi-temporal remote sensing; mixed pixel decomposition model; decision tree classi fi cation; karst mountainous area

S771.8

A

1673-923X(2017)07-0011-07

10.14067/j.cnki.1673-923x.2017.07.002

2016-03-15

国家“863”项目（2012AA102001）；广西林业厅项目（201438）

黎良财，副教授，博士研究生

张晓丽，教授，博士生导师；E-mail：Zhang-xl@263.net

黎良财，陆灯盛，张晓丽，等. 基于时序遥感的喀斯特山区植被覆盖研究[J].中南林业科技大学学报，2017, 37(7): 11-17.

[本文编校：吴 毅]