南京市植被覆盖管理措施因子的时空格局动态变化

林 杰，董 波，潘 颖，杨 敏，朱 茜

（南京林业大学江苏省南方现代林业协同创新中心/江苏省水土保持与生态修复重点实验室，江苏 南京 210037）

土壤侵蚀是全球范围发生最广、危害最严重的世界十大环境问题之一。植被覆盖管理措施因子（vegetation cover and management factor，以下简称C因子）用来表征植被对土壤的防蚀作用，是评价土壤侵蚀的关键参数，通用土壤侵蚀方程（universal soil loss equation，USLE）和修正的通用土壤侵蚀方程（revised universal soil loss equation，RUSLE）模型均将其作为重要的输入参数［1］。目前，利用遥感数据大面积、快速、实时提取C因子的主要方法包括基于遥感分类的直接赋值法、基于光谱混合分析模型的估算法和基于植被覆盖度的估算法等。WIJESUNDARA等［2］基于前人在伊斯兰卡地区2001和2013年的实地测量值，根据24种土地利用类型赋予研究区不同的C值，但由于忽略了地形所导致的非均匀植被覆盖和粗略的地图分类问题，利用这种方法得出的C值精度较低。LU等［3］利用ETM影像通过混合像元分解法将地表分为裸土、绿色植被和阴影3个部分，然后建立C因子与这些组分的关系式，进而估算区域C因子值，结果表明混合像元分解法能够较为全面地反映C因子的空间变异信息，但是混合光谱分析本身有一些缺陷，如地物分类数有限、存在多重散射等。

长期以来，在研究植被防蚀机制方面，一直使用水平植被覆盖度作为 其 评价的主要指标［1，4］。VATANDAŞLAR 等［1］和 RAWAT 等［4］分别在半干旱流域和印度地区利用高分辨率影像提取归一化植被指数（normalized difference vegetation index，NDVI）计算C因子值，从而计算研究区水土流失情况，结果表明利用传统植被覆盖度提取的C因子值比利用传统赋值法更精确。但利用传统植被覆盖度时并没有综合考虑林地植被在空间分布上的复杂结构特征，无法完全地表征地表枯枝落叶层。有学者认为叶面积指数（leaf area index，LAI）是理想的土壤侵蚀评价指标［5-7］。因此，如何提高区域尺度C因子的估算精度则成为利用USLE/RULSE预测土壤侵蚀准确性的关键［8］。

林杰等［5］基于统计模型，构建了LAI与C因子的量化耦合模型，结果表明LAI作为水土保持定量估算与评价的指标是可行的。虽然统计模型法的使用简单快捷，但易受土壤背景等多种外在因素的影响，反演精度往往不高，缺少可移植性［9］。尽管物理模型方法理论基础完善，模型参数物理意义明确，并可对作用机制进行适当的数学描述，但此类模型一般是非线性的，输入参数多，方程复杂，适用性较差，且对非主要因素有过多的忽略或假定［10］。而神经网络在对复杂、非线性数据进行拟合及模式识别方面有着无可比拟的优势［11］。杨敏等［12］基于 2013年 Landsat 8 Operational Land Imager（Landsat 8 OLI）数据构建了基于后向反射（back propagation，BP）神经网络的LAI反演模型，反演精度明显高于非线性回归模型，模型有较高的空间可靠性。

在区域尺度对C因子进行准确估算仍然是个难点，而区域尺度上的高质量时间序列C因子的合理估算及空间分布特征是预测流域与区域土壤侵蚀动态变化的关键环节。南京市是长江下游地区重要的城市之一，随着经济的快速发展，各类不合理的生产建设活动造成不同程度的水土流失，土壤侵蚀状况较为严峻［13］。马力等［14］和邵方泽等［15］基于RUSLE模型计算了南京市2001—2008年和2006—2014年水土流失及空间分布情况，结果表明近年来由于南京市水土保持措施的实施和植被覆盖的增加，全市水土流失控制效果显著，生态环境明显改善，但人类活动仍增加了部分地区的土壤侵蚀强度。因此，笔者以南京市为研究区，采用BP神经网络方法，基于研究区1988—2013年遥感数据，构建LAI和C因子的反演模型，利用137Cs同位素示踪技术分析并评价C因子遥感定量反演模型的反演精度，最终获取区域尺度高质量时间序列C因子，对比分析了年际间C因子变化及其驱动因子，以期为区域尺度土壤侵蚀动态变化的遥感监测提供新方法。

1 研究区概况与数据资料

1.1 研究区概况

南京市位于长江下游中部地区，江苏省西南部，介于31°14"～32°37"N、118°22"～119°14"E之间。全市行政区域总面积约为6 587.02 km2。地形以低山、丘陵为骨架，以环状山、条带山和箕状盆地为主要特色。气候属于北亚热带季风气候区，四季分明，雨量充沛，年平均气温为16℃，年平均降水时间为117 d，年平均降水量为1 106 mm，无霜期为237 d［16-18］。南京北、中部广大地区土壤为黄棕壤（地带性土壤），南部与安徽省接壤处有小面积红壤。植被类型属于常绿落叶阔叶混交林。南京地区人口密集，属于农业活动强烈区，自然植被在历史上屡遭严重破坏，几乎全部消失，现有植被多属次生性质，其中人工林面积大于自然恢复的次生林。境内现有林业用地约为840 km2，其中用材林和生态林约为570 km2，经济林和竹林约为210 km2。用材林和生态林树种主要有马尾松、黑松、杉木、国外松、麻栎、刺槐、水池杉和柏类等，经济林以茶果桑为主，竹林以毛竹为主，集中于丘陵山区。

1.2 LAI野外实测数据

LAI数据的野外测定采用美国LI-COR LAI-2200植物冠层分析仪。野外实测数据采集均在7—9月进行，与遥感影像获取时间基本同步。其中，2009年样点数为19个，主要分布在铜山林场；2010年样点数为37个，主要分布在江宁区和铜山林场；2013年样点数为111个，主要分布在聚宝山、紫金山、幕府山、将军山、吉山、方山、东善桥林场、朱门山、铜山和老山等地。所有样点测定时间均选择在6：30—9：00或16：30—19：00，尽量避免因太阳光线直射而引起的测试误差。分别在每个样地的4个角和中心位置各测量1次，取5次平均值作为测试结果。鉴于遥感影像的空间分辨率为30 m，因此设置采样间距均大于30 m，每个样点均由GARMIN手持GPS接收机定位，坐标系为WGS-84，各样点重复测定2次，取其平均值作为测试结果。

1.3 遥感影像数据及预处理

1.3.1 遥感影像数据

所选用的遥感影像数据为Landsat 5 Thematic Mapper（Landsat 5 TM）和 Landsat 8 OLI数据，共 10景影像数据。10景影像成像时间分别为1988-07-05、1994-07-22、2000-10-10、2002-07-12、2004-10-21、2006-08-08、2007-07-26、2009-10-03、2010-08-19和2013-08-11；除成像时间2013-08-11对应的传感器为LANDSAT_8_OLI外，其余成像时间对应的传感器均为LANDSAT_5_TM；10景影像的分辨率均为30 m，轨道号均为120/38。影像选取要求如下：选取时间为当年7—9月，是植被生长最茂盛且降雨最集中的月份，也是最容易发生水土流失的月份；为了保证反演精度，研究区还要尽量选取无云覆盖的高质量影像，但由于2000、2004和2009年7—9月没有符合质量要求的影像，故这3 a选择了10月的影像。

1.3.2 其他数据源

南京市1986、1996、2002和2013年4期土地利用数据中，前3期土地利用数据由中国科学院提供，2013年土地利用数据是以2013年Landsat 8 OLI数据为基础，参考Google Earth高分辨率影像和2002年土地利用数据以及国家科技基础条件平台——地球系统科学数据共享平台获取2013年土地覆被数据，进行人机交互解译获得，经实地采样验证，解译精度在88%以上。

1.3.3 遥感数据预处理

大气校正是通过遥感影像获取地面真实反射率必不可少的步骤，对植被定量遥感尤为重要。采用FLAASH（fast line-of-sight atmospheric analysis of spectral hypercubes）模型进行大气校正。同时，以校正好的南京市2007年TM遥感影像（高斯投影，椭球体为krosovsy）为标准底图，采用二次多项式拟合法，对影像进行几何精纠正，误差控制在0.5个像元内。最后，用南京市行政区划界线矢量数据进行裁切，得到研究区10景影像数据。

2 研究方法

2.1 LAI反演算法

Landsat 8 OLI多光谱遥感影像数据包括ETM+（enhanced thematic mapper）传感器所有波段，只对部分波段进行调整，增加蓝波段（Band 1，0.433～0.453 μm）和短波红外（Band 9，1.360～1.390 μm）2个波段。其中，调整比较大的波段是OLI Band 5（0.845～0.885 μm），排除0.825 μm处水汽吸收特征。因此，基于2013年Landsat 8 OLI数据所建的模型适用于前25 a的TM-LAI估算。应用隐含层为2层的BP神经网络对1988—2013年南京市LAI进行反演［12］，利用2009和2010年LAI实测数据检验反演精度。为了排除不同时期影像间传感器的差异和大气等因素的影响，将BP神经网络训练和模拟的归一化部分做了相应调整，即对每景影像分别做归一化处理，以减小反射率差异对反演精度造成的影响。归一化公式如下：

式（1）中，Xi，max和 Xi，min分别为第 i个神经元各输入分量的最大值和最小值；Xi和Xi"分别为第i个神经元预处理前、后的输入分量。

2.2 植被覆盖度

2.2.1 传统植被覆盖度

传统植被覆盖度是基于NDVI根据像元二分模型法［19］计算得到的，计算公式为

式（2）中，fc为传统植被覆盖度；INDV，min为裸土或无植被覆盖像元的NDVI值；INDV，max为全植被覆盖像元的NDVI值，即纯植被像元的 NDVI值。INDV，max和 INDV，min通常根据图像大小、图像清晰度等情况，以置信度截取NDVI的上、下限阈值分别近似代表INDV，max和INDV，min，在一定程度上消除遥感影像噪声所带来的误差［20］。根据实际数据情况最终采用5%和95%的累计百分比为置信区间来计算 INDV，max和 INDV，min。

2.2.2 植被方向性覆盖度

植被方向性覆盖度与LAI之间存在一个经典物理意义上的转换公式。由CHEN等［21］和KUUSK等［22］提出，计算公式为

式（3）中，Fcover（θ）为植被方向性覆盖度；P（θ）为冠层空隙率；G（θ）为叶片在太阳入射方向上的投影，表征叶子截光能力的大小；ILA为叶面积指数；θ为太阳入射天顶角，当植被叶片呈球形随机分布时，G（θ）＝0.5，式（3）可简化为

因此，只要从遥感影像上反演出高精度的长时间序列LAI，即可准确地估算出植被乔、灌、草植被垂直覆盖度的动态变化。

2.3 植被覆盖管理措施因子C

蔡崇法等［23］通过径流小区的多场人工降雨和天然降雨资料与地表覆盖度之间的相关关系，采用回归分析方法建立C因子与植被覆盖度fc（以百分数表示）之间的关系式：

通过LAI估算植被垂直覆盖度，进而估算植被覆盖管理措施因子C。

2.4 植被覆盖管理措施因子C值的实测值

2.4.1 基于137Cs的土壤侵蚀计算模型

土壤样品采集于2014年5月，采样点选在南京市铜山林场的山坡和附近耕地农田，共32个采样点，其中，背景值点8个，水稻田、菜地、林地、竹林和茶园的采样个数分别为4、4、10、3和3，共160个土壤样品。样品137Cs含量测定均采用标准源对比法，采用中国原子能研究院标准源质量和放射性活度。根据式（6）计算样品中137Cs的放射性活度，进而分耕作和非耕作土壤测定土壤流失速率［24］，然后根据耕作或非耕作土壤侵蚀模型［25-26］计算年土壤流失厚度。

式（6）中，A样为样品放射性活度，Bq·g-1；A标为标准源放射性活度，Bq·g-1；T标为标准源测量时间，s；W标为标准源质量，g；S标为标准源计数；S样为样品测量计数；T样为样品测量时间，s；W样为样品质量，g。

（1）非耕作土壤侵蚀模型

在自然条件下，137Cs在未扰动稳定的土壤剖面中的深度分布是随着深度逐渐减小的，选择张信宝［25］的非农耕地土壤137Cs深度分布指数衰减形态模型，土壤流失速率（RSL，i， t·km-1·a-1）计算公式为

式（7）～（8）中，Δh为年土壤侵蚀厚度，cm·a-1；ρb为土壤容重，kg·cm-3；T为采样年份，a；λ为137Cs分布的剖形指数，cm-1，取值0.27；A和Aref分别为单位面积土壤中137Cs含量和土壤137Cs背景值，Bq·m-2。

（2）耕作土壤侵蚀模型

选择杨浩等［26］基于质量平衡的农耕地土壤侵蚀模型。由于人为翻耕，137Cs在耕层中呈均匀分布，但分布深度大多比耕层深，则无侵蚀耕作剖面中Cref和侵蚀耕作剖面中Ct均分为两个部分：一部分是均匀分布于耕层的Cref1（或Ct1），另一部分是入渗于耕层以下的Cin，即：

式（9）～（10）中，ρb为土壤容重，kg·cm-3；Q为质量活度浓度，Bq·kg-1；Hc为耕层厚度，cm；h为土壤流失总厚度，cm。

而137Cs流失量△C=Cref-Ct＝ρb×Q×h，整理得：

t时间以来土壤年均流失厚度Δh"（mm）和相应的年均侵蚀速率即土壤流失速率（RSL，i"）为

式（12）～（13）中，T为采样年份，a；T＞1 963 a时，1 963 a为137Cs输入最大年份；Δh"为年均土壤侵蚀厚度，cm·a-1；ρb为土壤容重，kg·cm-3。

2.4.2 C因子的实测值

根据唐克丽［27］对C因子的定义，其计算公式为

根据式（6）～（14），计算得到铜山林场1963—2014年土壤流失速率和C因子值（表1）。如表1所示，铜山林场林地、茶园、竹林、水稻田和菜地C因子均值分别为 0.074 1、0.061 2、0.044 8、0.308 6和0.756 5。各土地利用类型C因子值由大到小依次为菜地、水稻田、茶园、竹林和林地。

2.5 精度评价方法

平均相对误差（MAPE，EMAP）、均方根误差（RMSE，ERMS）和相关系数（R）被用来衡量和刻画模型精度，其计算公式分别为

式（15）～（17）中，xi(c)为模型反演值C；xi为C因子实测值；-x为C因子实测值的平均值；n为样本个数。

表1 铜山林场不同土地利用类型C因子值计算结果Table 1 Calculation results of C factor values of different land use types in Tong Mountain Forest Farm

3 结果与分析

3.1 植被覆盖管理措施因子C的精度评价

传统植被覆盖度（fc）和植被方向性覆盖度（Fcover）计算结果见图1～2，基于NDVI和基于LAI的C因子估算结果见图3～4。利用2013年铜山林场基于137Cs同位素示踪获取的C因子实测值检验2种不同方法反演的精度（表2）。

由图1～4可知，基于NDVI计算的fc在整体上偏大，得到的C值在整体上偏小，主要是由于当植被覆盖度较高时，NDVI对植被的敏感性降低，而且基于NDVI计算的fc反映的是水平面上的植被覆盖度，不包含植被的垂直结构信息。而南京市多为次生植被，以人工林为主，如马尾松纯林、杉木纯林等，林下植被覆盖较少，植被结构简单，虽然在水平面上有较好的植被覆盖度，但也会遭受严重的土壤侵蚀。而LAI不仅包含植被的水平覆盖信息，还包括植被的垂直结构信息。

图1 基于NDVI的传统植被覆盖度分布Fig.1 The traditional vegetation coverage fcdistribution map based on NDVI

图2 基于LAI的植被方向性覆盖度分布Fig.2 The directional vegetation coverage Fcover distribution map based on LAI

图3 基于NDVI的C因子分布Fig.3 The C factor distribution map based on NDVI

图4 基于LAI的C因子分布Fig.4 The C factor distribution map based on LAI

基于LAI计算的Fcover是水平面上植被覆盖和植被垂直结构信息的综合反映。表2显示，当基于fc反演的C值为0时，C实测值并不为0，说明基于fc来定量评价土壤侵蚀并不准确。有些地区通过模型计算的土壤流失量很小，土壤侵蚀强度很轻，但事实上由于地表缺乏灌木和草本覆盖，土壤裸露程度较高，仍然会发生中度甚至强烈以上等级的土壤侵蚀，这与赵其国［28］提出的“远看青山在，近看水土流”的情况相吻合。而基于LAI计算的C值与实测值接近，平均RMSE为30.017%，C值反演精度更高，说明LAI能综合反映林地植被在空间分布上的复杂结构特征，可以完全表征地表枯枝落叶层，研究结果与杨勤科等［6］提出的观点一致。

表2 基于BP神经网络的C因子反演精度Table 2 The accuracy ofCfactor inversion based on BP neural network

通过叠加2013年南京市土地利用图（图5）并进行统计分析发现：C值大于0.3的区域主要分布于建筑物较为密集、植被稀疏且植被结构简单甚至无植被覆盖的市区；C值小于0.05的区域主要分布在植被密集且结构复杂的丘陵山区；有较高C值的区域主要集中在土地利用类型为居民点和城镇建设用地的地区，尤其是新开发的生产建设用地，其次是山地丘陵区的工矿用地，未利用地也是土壤侵蚀较严重的地区；C值基本为0的区域主要是植被覆盖较好且植被结构复杂的山区。由此可见，南京市C值的分布与植被覆盖和土地利用类型有着十分密切的关系。

3.2 长时间序列植被覆盖管理措施因子C的反演

利用BP神经网络和Fcover反演C值方法得到1988—2013年南京市C因子图。2013年南京市C因子遥感估算结果见图6，其他年份C因子图从略。

图5 2013年南京市土地利用分布Fig.5 The land use distribution map of Nanjing in 2013

图6 2013年南京市C值遥感估算结果Fig.6 Vegetation cover and management factor C by remote sensing inversion in 2013

为了定量反映1988—2013年不同数值范围C值的变化趋势，统计分析了全市C＜0.05和C＞0.3的区域面积变化趋势（图7）。如图7所示，从全市来看，1988—2013年C＜0.05的抵抗土壤侵蚀能力较强的区域面积先由南京市总面积的15.66%（1988年）减小到9.43%（2006年以前），后逐渐增大到12.07%；而C＞0.3的抵抗土壤侵蚀能力较弱的区域面积先由南京市总面积的7.29%缓慢增大到9.22%（2002年）后迅速增大到12.31%（2002—2006年），然后缓慢减小至11.77%，这主要是受城市发展扩张的影响。20世纪80—90年代，南京市城区不断扩张，面积扩大，土地开发等生产建设活动使植被遭受严重破坏，植被覆盖度降低，植被垂直结构遭到破坏，林分类型较好的区域逐渐减少，导致C＜0.05的区域面积不断减少，而居民地等开发建设用地区域面积不断增加，C＞0.3的区域面积不断增加；进入21世纪后，城市化速度加快，房地产开发加剧，城市面积迅速扩张，导致C＞0.3的区域面积迅速增加；而到2006年前后随着城市化的发展，国家政策以及人们意识的改变，“和谐社会”和“生态文明建设”等社会发展战略目标的提出，在加速城市建设的同时，注重植被保护和恢复，植被覆盖度逐年增加，林下植被也得到一定恢复，使得C＜0.05的区域面积不断增加。由于已建成居民地等区域面积不会减少，随着增加绿化面积等措施的实施，各类植被类型的水土保持功能有所增强，C＞0.3的区域面积缓慢减少。

图7 1988—2013年C值像元数百分比的年际变化趋势Fig.7 The interannual change trend of the pixels number percentage of C value from 1988 to 2013

根据土地利用分类标准，利用南京市江宁区2005年土地利用图，采用ArcGIS 10.2软件将其与C因子图进行叠加，统计林地、疏林地和灌木林3种用地类型的平均C值。由于2000、2004和2009年遥感影像成像时间为10月，为了避免季节不同引起的C值差异，统计分析了 1988、1994、2002、2006、2007、2010和2013年7 a的变化趋势，结果见图8。

图8 南京市江宁区不同林地类型C平均值的年际变化趋势Fig.8 The interannual change trend chart of C mean value under different woodland types in Jiangning District，Nanjing City

如图8所示，不同林地利用类型C因子由小到大依次为林地、灌木林和疏林地，基本符合C因子的实际分布特征。一般来讲，森林的水土保持功效最大，灌木林次之，疏林地最差。从整体来看，1988—2013年江宁区林地、灌木林和疏林地C平均值大体上呈现先增大后减小趋势，表明3种林地水土保持功效先降低后增加。根据C值统计结果，结合江宁区2005年土地利用分类图可知，在2005年为疏林地的区域，1988、1994和2002年C因子的变化范围最小值为0，说明该区域在2005年之前可能是林分类型较好的林地，植被结构复杂，植被覆盖度较高，导致这3 a平均C值小于2005年的计算结果。同样地，2005年之后的C值变化范围最大值可达1，表明该区域在2005年之后可能已经不是林地，土地处于裸露状态，导致平均C值增大。影响C值变化的主要原因是城市快速发展，尤其是2000年12月，南京市江宁县经国务院批准，撤县立区，城市快速发展，房地产开发项目大量增加，城市空间不断扩张，导致不同类型植被遭受不同程度的破坏，C值逐渐增大。而在2006年之后，随着“绿色南京”工程、生态旅游建设、美丽乡村建设等政策的实施以及人们环保意识的增强，植被有所恢复，植被种类开始增多，各植被类型保持水土的功能有所增加，C值逐渐减小。因此，C值的变化在不同程度上反映了江宁区土地利用的变化。

4 讨论与结论

基于NDVI和LAI计算了fc和Fcover，分别对C进行反演，并利用基于137Cs同位素示踪技术计算的C因子实测值检验反演精度，最后，基于LAI对南京市1988—2013年10期遥感影像进行长时间序列的植被覆盖与C的反演，得出以下结论：

（1）基于NDVI计算的fc整体上偏大，C值整体上偏小，而基于LAI计算的C值与实测值接近，平均RMSE为30.017%，能更好地反映实际土壤侵蚀状态。

（2）C值大于0.3的区域主要分布于建筑物较为密集、植被稀疏且植被结构简单甚至无植被覆盖的市区；C值小于0.05的区域主要分布于植被密集且植被结构复杂的山区，南京市C值的分布与植被覆盖和土地利用类型关系密切。

（3）从全市来看，1988—2013年C＜0.05的抵抗土壤侵蚀能力较强的区域面积先减小（2006年之前）后增大，C＞0.3的抵抗土壤侵蚀能力较弱的区域面积先缓慢增大（2002年之前）后迅速增大（2002—2006年），然后又缓慢减小。

（4）南京市江宁区不同土地利用类型C平均值由小到大依次为林地、灌木林和疏林地，符合C因子的实际分布特征；江宁区林地、灌木林和疏林地C因子平均值大体上呈现先增大后减小趋势，表明3种林地的水土保持功效先降低后增加，符合南京市植被覆盖和土地利用的发展变化规律。

由于TM和OLI遥感数据只能提供一个方向的数据，缺乏足够的信息表征地表植被结构特征，多角度的反射光谱对地物（特别是植被）结构特征的估算及类型鉴别比垂直光谱具有明显优势。多角度遥感可以获取地面不同方向的反射数据，同时可以模拟和反演植被二向反射率函数，丰富地物的立体结构特征，有助于分析和提取植被垂直分层覆盖度信息，因而能够有效地提高LAI等植被结构参数的反演精度。因此，利用多角度遥感数据及物理模型，建立多角度LAI反演方法，构建不同植被结构LAI与C因子耦合模型，可为区域土壤侵蚀定量遥感监测提供新途径，也将是未来研究趋势。