随着社会经济的快速发展，人类活动对地表植被造成了极大的破坏，水土流失的速度不断加快，严重制约了我国生态环境建设与可持续发展^[1-4]，水土流失问题成为广大学者关注的环境问题之一。水土流失的调查与评价是解决水土流失问题的基础，受到国内外学者的高度关注^[5-7]。目前对水土流失的调查与评价方面的研究主要分为定性与定量2方面，定性多以多因子综合评判和专家经验法为主，判断水土流失发生发展的可能^[8]。定量则是以经验模型和物理模型为主，分析评价区域水土流失状况。多使用通用土壤侵蚀方程(universal soil loss equation，USLE)^[9-10]、修正版的土壤侵蚀模型(revised universal soil loss equation，RUSLE)^{[5, 11]}以及中国土壤侵蚀模型(Chinese soil loss equation，CSLE)^[12]分析探讨土地侵蚀的空间分布及变化特征。但大多数研究选取定性或定量的单一方面进行水土流失评价，通过野外调查定性分析结合经验模型或物理模型的定量计算来评价区域水土流失状况的研究并不多见。

水利部《关于加强水土保持监测工作的通知》(水保〔2017〕36号)要求，“各地根据社会关注重点和实际工作需要，有针对性地开展生态脆弱地区、禁止开垦陡坡地、湖泊和水库周边植物保持带等区域的监测工作，保障生态脆弱地区不遭受破坏、陡坡地面积不增加、植物保持带生态功能不降低”。根据这一要求，笔者以九洲江流域为研究对象，采用遥感监测与野外实地调查相结合的方法，利用CSIE^[13-14]，重点围绕区域土地利用、植被覆盖度和水土流失分布特征、面积和强度等方面开展监测，以期了解和掌握九洲江流域水土流失与生态环境状况，为后期防治对策的制订提供数据支撑。

1 研究区概况

九洲江是一条跨桂粤2区(省)的入海河流，地处亚热带气候带，一年四季温暖湿润，发源于广西壮族自治区陆川县沙坡镇秦镜村的文龙径，流经广西玉林市陆川、博白2县共10个镇进入广东湛江市的鹤地水库，最后注入北部湾。其全长168 km，是典型的雨源型河流，4—9月汛期雨量占全年降雨总量的85%，年径流量平均为28亿m³，但枯水期许多河段接近断流，集雨面积3 396 km^2[15]。流域地势以低丘为主，中下游两岸为狭长的河谷平原，受洪水威胁的农田面积约5 300 hm²，下游缸瓦窑站(集水面积3 050 km²)于1906年曾出现洪峰流量5 680 m³/s。研究区位于广西境内，坐标为E 109°57′~110°23′、N 21°51′~22°22′。行政区涉及陆川县的沙湖、温泉、沙坡、大桥、横山、乌石、滩面、良田、古城9个乡(镇)，以及博白县的宁潭、英桥、凤山文地、新田、大垌6个乡(镇)(图 1)。

2 数据与方法 2.1 数据采集与处理

购买2019年度GF-1遥感影像(2 m分辨率)，2017、2018年度MOD13Q1(250 m分辨率)和Landsat 8 OLI(30 m分辨率)通过地理空间数据云平台(http://www.gscloud.cn/)和GloVis(https://glovis.usgs.gov/)下载。基于ENVI和ArcGIS平台，将影像投影统一转换为正轴等面积割圆锥投影(Albers投影)，坐标系统为CGCS2000，采用1985国家高程基准，1 ∶1万比例尺按3°分带。结合Google Earth、流域内土地利用类型、水土保持重点工程、生产建设项目批复等基础资料，采用面向对象+人工交互解译的方式对高分遥感影像进行初解译，根据GB/T 21010—2017《土地利用现状分类》^[16]、影像初分类结果(图斑总数量与图斑总面积)，以满足“对不小于解译结果总图斑数的抽样5%的核查对象，抽取10%作为验证样本进行实地验证”为最低目标，进行图斑抽样野外调查，利用无人机现场记录区域内的土地利用类型、植被覆盖度等信息，并与遥感影像颜色、纹理特征进行对比，建立解译标志，对九洲江流域的土地利用类型、植被覆盖度进行校正。

2.2 研究方法 2.2.1 植被覆盖度因子提取

计算1年内30 m空间分辨率OLI影像归一化植被指数(normalized difference vegetation index，NDVI，量纲为1)，与1年23个半月250 m空间分辨率MODIS影像NDVI值，通过融合、插值、重采样获取每年24个半月30 m空间分辨率的植被覆盖度^[17]。

$ {\rm{NDVI}}= \frac{{\rm{NIR}}-R} {{\rm{NIR}}+R} 。$

式中：NIR为近红外波段的反射率，nm；R为可见红光波段反射率，nm。利用以下公式将融合的24个半月30 m空间分辨率NDVI转换为相应的植被覆盖度(fractional vegetation cover，FVC，量纲为1)^[18]：

$ \mathrm{FVC}=\left(\frac{\mathrm{NDVI}-\mathrm{NDVI}_{\min }}{\mathrm{NDVI}_{\max }-\mathrm{NDVI}_{\min }}\right)^k 。$

式中：NDVI_max为当年监测区的植被覆盖度最大区域NDVI值(不能大于1)；NDVI_min为裸土NDVI值；k为系数。

2.2.2 土地利用类型因子提取

将无人机影像与高分影像进行配准，勾绘典型土地利用类型地块，建立九洲江流域土地利用遥感解译标志。基于土地利用遥感解译标志，利用归一化差分水体指数(normalized difference water index，NDWI)^[19]，量纲为1，设置阈值(NDWI＞0.1)提取九洲江流域水体，即河湖库塘地类；利用NDVI，设置阈值将剩下的非水体地类分为植被地类(NDVI＞0.1，包括耕地、园地、林地、草地)、非植被地类(NDVI≤0.1，包括建设用地、交通运输用地、其他土地)；再将植被地类与非植被地类分别设置参数，通过面向对象分类+人工交互解译的方式进行精细分类^[20]，对于土地利用和水土保持措施解译结果，总体上抽取不少于总图斑数的10%进行核查；对于核查图斑，抽取10%作为验证样本进行实地验证。验证点位置与遥感影像偏差小于1个遥感影像像元，图斑属性判对率应大于90%。研究区土地利用类型划分水田、旱地、其他园地、有林地、其他林地、其他草地、城镇建设用地、农村建设用地、采矿用地、其他建设用地、其他交通用地、河湖库塘和裸土地13类。

$ \mathrm{NDWI}=\frac{G-\mathrm{NIR}}{G+\mathrm{NIR}}。$

式中G为可见绿光波段反射率，nm。

2.2.3 水土流失因子值计算

降雨侵蚀力因子R、土壤可蚀性因子K计算采用水利部水土保持监测中心收集加密雨量站点降水资料、更新计算降雨侵蚀力因子和土壤可蚀性因子，统一下发各流域管理机构和省级水行政主管部门用于土壤侵蚀计算。

坡长因子L、坡度因子S计算基于九洲江流域10 m空间分辨率的数字高程数据(DEM)，进行坡长坡度计算，获得流域坡长因子、坡度因子栅格图。坡长因子

$ L_i=\frac{\lambda_i^{m+1}-\lambda_{i-1}^{m+1}}{\left(\lambda_i-\lambda_{i-1}\right) \cdot(22.13)^m }。$

式中：λ_i、λ_i-1分别为第i个和第i-1个坡段的坡长，m；m为坡长指数，随坡度而变，量纲为1^[21]。

$ m=\left\{\begin{array}{ll} 0.2 & \theta \leqslant 1^{\circ} \\ 0.3 & 1^{\circ}<\theta \leqslant 3^{\circ} \\ 0.4 & 3^{\circ}<\theta \leqslant 5^{\circ} \\ 0.5 & \theta>5^{\circ} \end{array}\right.。$

坡度因子

$ S=\left\{\begin{array}{ll} 10.8 \sin \theta+0.03 & \theta<5^{\circ} \\ 16.8 \sin \theta-0.5 & 5^{\circ} \leqslant \theta<10^{\circ} \\ 21.9 \sin \theta-0.96 & \theta \geqslant 10^{\circ} \end{array}\right.。$

式中θ为坡度，(°)。生成的L、S栅格数据空间分辨率为10 m。

植被覆盖与生物措施因子B计算采用融合法与单期高分影像分别计算植被覆盖与生物措施因子，根据解译的土壤侵蚀地块属性表中土地利用类型计算各地类B因子值，生成B因子栅格图层。园地、林地和草地B因子计算公式为：

$ B=\sum\limits_{i=1}^{24} \mathrm{SLR}_i \times \mathrm{WR}_{i }。$

式中：WR_i为第i个半月降雨侵蚀力占全年侵蚀力比例，取值范围为0~1；SLR_i为第i个半月园地、林地和草地的土壤流失比例，量纲为1，取值范围为0~1。茶园和灌木林地SLR_i计算公式为：

$ \mathrm{SLR}_i=\frac{1}{1.17647+0.86242 \times 1.05905^{100 \times \mathrm{FVC}} }。$

果园、其他园地、有林地和其他林地SLR_i计算公式为：

$ \begin{array}{l} \mathrm{SLR}_i=0.44468 \times \mathrm{e}(-3.20096 \times \mathrm{GD})- \\ 0.04099 \times \mathrm{e}(\mathrm{FVC}-\mathrm{FVC} \times \mathrm{GD})+0.025。\end{array} $

草地SLR_i计算公式为：

$ \mathrm{SLR}_i=\frac{1}{1.25+0.78845 \times 1.05968^{100 \times \mathrm{FVC}}}。$

式中：GD为乔木林的林下盖度，取值范围为0~1，包括除乔木林冠层以外的所有植被(灌木、草本和枯落物)构成的林下盖度，按实地调查和经验，有林地GD取值为0.6~0.8，其他园地GD取值0.4~0.6，其他林地GD取值0.2~0.4。非园地、林地、草地地类的B因子查《非园地、林地、草地的B因子赋值表》 ^①。

① 源自2020年度水土流失动态监测技术指南。

工程措施因子E计算参考2018年广西自治区水土流失动态监测工作中九洲江流域E因子成果^②，将水田地类作为工程措施，并进行赋值，生成2 m空间分辨率(对应土地利用分类结果)的E因子栅格数据。

② 源自2018年广西自治区水土流失动态监测工作中九洲江流域E因子成果

耕作措施因子T赋值采用九洲江流域所涉及区县(即玉林市陆川县、博白县)耕作措施因子0.417，对2°以下(含2°)耕地考虑等高耕作措施，因子值统一为0.431 ^③。

③ 源自《耕作轮作措施赋值表》

2.2.4 土壤侵蚀模数计算

将降雨侵蚀力R、土壤可蚀性因子K、坡长因子L、坡度因子S、植被覆盖与生物措施因子B、工程措施因子E、耕作措施因子T，代入中国土壤流失CSLE方程^[20]，进行栅格运算。方程如下：

$ A=RKLSBET。$

式中A为土壤侵蚀模数。当土地利用类型为耕地时，在植被覆盖与生物措施因子、耕作措施因子中，选取耕作措施因子与其他5个因子图层相乘；当土地利用类型为非耕地时，则选取植被覆盖与生物措施因子，与其他5个因子图层相乘，获取10 m空间分辨率土壤侵蚀模数栅格数据。采用综合评判法，依据SL190—2007《土壤侵蚀分类分级标准》^[22]，将九洲江流域土壤侵蚀模数计算结果进行分级，定性和定量相结合分析评价九洲江流域水土流失分布特征、面积和强度。

3 结果与分析 3.1 土地利用现状

2019年九洲江流域土地利用类型(图 2)以有林地和水田为主，2种土地利用类型面积分别为778.47和144.37 km²，分别占流域总面积的62.0%和11.5%，其余土地利用类型占流域总面积的26.5%(表 1)。从空间分布来看，有林地占流域面积最大，水田沿河道两侧分布。建设用地分布较为集中，形成了上、中、下游3段集聚区。

3.2 植被覆盖度

2019年九洲江流域大部分地表都有植被覆盖且植被覆盖度都较高，园地、林地和草地3种主要的植被地类总面积952.68 km²，占流域总面积的75.9%；高覆盖园地和林地占园地和林地总面积的92.3%，中覆盖度以下(包括中覆盖度)园地和林地仅占园地和林地总面积的0.5%；草地方面，高覆盖草地占草地总面积的96.5%，中覆盖度以下(包括中覆盖度)草地仅占草地总面的0.2%(图 3和表 2)。从空间分布来看，植被覆盖度沿河道两侧逐渐升高，沿河道越近，植被覆盖度越低。这与土地利用现状相一致，沿河道两侧多分布人为活动较高的建设用地，人为活动越强烈的区域植被覆盖度越低。

3.3 水土流失及变化 3.3.1 土壤侵蚀模数及土壤侵蚀量

2019年九洲江流域平均土壤侵蚀模数为532 t/(km² ·a)。从空间分布上来看，高土壤侵蚀模数集中出现在地表扰动剧烈的城镇建设、采矿、公路建设等人为活动造成的裸土地区域；中、中高土壤侵蚀模数则集中分布于流域的山地、丘陵地区；九洲江及其支流所在的中部则集中分布低、中低土壤侵蚀模数地块(图 4)。土壤侵蚀总量为64万1 834 t，其中剧烈地块土壤侵蚀量20万6 854 t，占流域总侵蚀量的32.2%；微度地块土壤侵蚀量13万3 255 t，占流域总侵蚀量的20.8%；中度地块土壤侵蚀量9万272 t，占流域总侵蚀量的14.1%(表 3)。

3.3.2 2018—2019年水土流失变化情况

与2018年研究区水土流失监测数据相比，2019年九洲江流域轻度及以上土壤侵蚀面积减少13.07 km²。在各类型土壤侵蚀强度消长变化方面，中度土壤侵蚀强度面积增加最多，为5.5%；强烈土壤侵蚀强度面积和微度土壤侵蚀面积也稍有增加，分别为2.1%和1.2%。剧烈、极强烈土壤侵蚀强度面积明显下降，分别减少27.4%和24.1%；轻度土壤侵蚀强度面积减少7.8%(表 4)。

2018年九洲江流域发生土壤侵蚀的地块(图 4)，在2019年的监测结果中都不同程度的向微度土壤侵蚀进行转化，其中轻度土壤侵蚀地块向微度转化的比例最多，占到2018年轻度土壤侵蚀总面积的78.7%；中度、强烈、极强烈、剧烈土壤侵蚀地块分别为57.2%、47.1%、48.7%、48.6%。

未发生土壤侵蚀强度变化的地块方面，有1 005.95 km²的地块在2018、2019年土壤侵蚀强度均为微度；2018年监测结果中的剧烈土壤侵蚀地块，有42.3%在2019年的土壤侵蚀强度监测结果中仍然为剧烈；极强烈、强烈、中度、轻度土壤侵蚀地块分别为41.1%、40.3%、26.4%、17.1%(表 5)。

4 结论与讨论

1) 2019年九洲江流域土地利用类型以有林地和水田为主，2种土地利用类型面积分别占流域总面积的62.0%和11.5%。大部分地表都有植被覆盖且植被覆盖度都较高。从空间分布来看，水田、建设用地主要沿河道两侧分布，沿河道距离越近，人为活动越强烈，植被覆盖度越低。

2) 2019年九洲江流域平均土壤侵蚀模数为532 t/(km² ·a)，流域土壤侵蚀总量为64万1 834 t。剧烈和极强烈的土壤侵蚀发生在人为活动强烈的采矿用地以及大规模开发建设项目造成的裸土地区域。

3) 与2018年全区水土保持动态监测结果相比，2019年九洲江流域土壤侵蚀面积减少13.07 km²，中度、强烈、微度土壤侵蚀强度面积有所增加，剧烈、极强烈、轻度土壤侵蚀强度面积下降较为明显。轻度土壤侵蚀地块向微度转化的比例最多，占到2018年轻度土壤侵蚀总面积的78.7%。有42.3%在2019年的土壤侵蚀强度监测结果中仍然为剧烈。人为活动造成的水土流失对生态环境的影响仍然较大。

研究从土地利用与植被覆盖2个角度出发，采用中国土壤侵蚀方程定性与定量分析九洲江流域水土流失现状，增进对九洲江流域水土流失的认识，为九洲江流域水土流失治理措施的制订提供了科学依据。研究的结果表明，人为活动集聚区普遍水土流失严重，而人为活动与水土流失的关系，是一个更为复杂和深入的过程，亟待更深层次的研究。另外，利用土壤侵蚀方程进行科学分析时应结合流域实际，警惕“异参同效”现象。