豫东北黄泛区2019年土壤侵蚀特征——以河南省兰考县为例

高睿瑜, 袁 利, 张荣华, 吴 迪, 李明宇, 李文龙

(1.山东农业大学 林学院 山东泰山森林生态系统国家定位研究站,山东 泰安 271018; 2.淮河水利委员会 淮河流域水土保持监测中心站， 安徽 蚌埠 233001)

黄泛平原风沙区(以下称“黄泛区”)属国家级水土流失重点预防区，长久以来，风蚀和水蚀造成土壤养分流失、土地生产力下降[1]，破坏区域生态环境，制约经济社会发展。黄泛区地处暖温带季风性气候区，雨热同期，风旱同季，地势整体平坦。区域在1—5月和10—12月同时存在强风和干旱天气，导致耕地、林草地和沙地上易出现风蚀。降水集中在6—9月，加之局部存在沙丘、沙岗等微起伏地貌[2]，易造成一定程度的水蚀。为摸清黄泛区风蚀规律，明确水土流失防治的的重点时间，应对区域土壤侵蚀分布格局、年内变化特征进行研究。

目前，国内关于黄泛区土壤侵蚀的研究侧重于风蚀，早期以试验分析和定位观测为主。比如赵存玉等[3]结合风洞模拟、野外调查阐明了鲁西北黄泛区农田与风蚀的作用机制；姬生勋等[4]证明了不同土地利用方式下土壤风蚀量变化规律存在差异；2017年，宋胜明等[5]通过定位观测研究了风速累积时间、农作物覆盖情况、地表粗糙度对兰考县耕地风蚀的影响。近10 a，遥感监测、模型计算被应用于黄泛区土壤侵蚀研究，内容涉及风蚀和水蚀。2012年，王友胜[2]借助遥感监测和野外调查研究了淮河流域黄泛区风水侵蚀格局及其驱动因子，为黄泛区土壤侵蚀后续研究提供了参考和有力依据；毛玉磊[6]利用遥感技术系统分析了河南黄泛区易风蚀性土地的现状及变化特征。2019年，高分辨率(≤2 m)遥感影像广泛应用于水土流失动态监测，袁利等[7]、张乐[8]以高分辨率遥感影像为数据源，利用风蚀模型、CSLE模型分别研究了淮河流域黄泛区、鲁西北黄泛区的水土流失情况。随着基础数据质量和遥感影像分辨率提高，县域尺度土壤侵蚀研究增多[9-10]，但目前黄泛区土壤侵蚀多以年为时间尺度进行现状和分布格局研究，有关土壤侵蚀年内变化特征关注较少，有待进一步丰富。

本研究以河南省兰考县为研究区域，以高分辨率遥感影像、MODIS-NDVI和SMAP等影像为基础数据，借助风力侵蚀模型和CSLE模型对兰考县土壤侵蚀进行系统研究，探索土壤侵蚀及其影响因子2019年年内半月变化特征，综合分析区域全年侵蚀特征。以期丰富黄泛区土壤侵蚀研究成果，为确定黄泛区及黄河下游地区风蚀治理的重点时间、地区提供科学依据，为全国水土流失动态监测工作提供数据支持。

1 研究区概况

河南省兰考县属于淮河流域黄泛区，地理位置位于34.44°—35.01°N和114.40°—115.16°E的范围内(具体情况见图1)，总面积为1 116 km2，包括仪封乡、三义寨乡、城关镇、东坝头乡等15个乡镇。研究区属暖温带半湿润季风性气候，四季分明，雨热同期，冬春季降水少，大风较多；冬季平均风速为1.71 m/s，春季平均风速1.89 m/s，2月风速值最高，为5.90 m/s～10.95 m/s[10]。区域地貌以平原为主，存在冲积扇、沙丘、沙岗等地貌景观；地形平坦，以平缓坡为主。豫东北黄泛区植被属暖温带落叶阔叶林带，种植农作物、农田林网、经济林、防风固沙林，常见树种为杨树、苹果等。农作物多为一年两熟制耕作，主要为小麦与花生、玉米或地瓜轮作；部分区域存在一年一熟制，是风蚀易发生区域。农田风蚀与风沙化是区域的主要水土流失问题[10]。

图1 河南省兰考县乡镇分区

2 研究方法

2.1 数据类型及来源

本研究所用数据包括遥感影像、地形地貌数据、土壤数据及气象数据和其他资料，具体情况见表1。

表1 河南省兰考县有关基础数据信息

2.2 土壤侵蚀计算

本研究侵蚀模型计算、侵蚀因子数据的获取与处理均以水利部水土保持监测中心制定的《区域水土流失动态监测技术规定(试行)》(以下称《技术规定》)为准则，并结合研究区域河南省兰考县实际情况进行。

2.2.1 侵蚀模型

(1) 风蚀模型。本研究借助第一次全国水利普查的风蚀模型计算2019年风蚀模数[7,11]，该模型源于中科院大田推广模型[11]。包括耕地、草(灌)地和沙地(漠)3种模型。模型在研究区的适用情况见表2，计算公式如公式(1)—(3)所示：

(1)

(2)

(3)

式中：Qfa,Qfg,Qfs分别为每半月耕地、草地和沙地风力侵蚀模数(t/hm2·a);W为每半月表土湿度因子,取值范围为0～1;Tj为每半月各风速等级的累积时间(min);Z0为地表粗糙度,无量纲;j为风速等级序号;Uj为第j个等级的平均风速(m/s);V为植被覆盖度(%)。

表2 研究区风蚀模型适用范围

(2) 水蚀模型。本研究运用刘宝元等[12]提出的的CSLE模型计算水蚀模数[7,12-15]，计算如公式(4)：

A=R·K·L·S·B·E·T

(4)

式中：A为土壤流失量(t/hm2·a);R为降雨侵蚀力因子(MJ·mm/hm2·h·a);K为土壤可蚀性因子(t·hm2·h/MJ·hm2·mm);L为坡长因子;S为坡度因子;B为植被覆盖与生物措施因子;E为工程措施因子;T为耕作措施因子。

2.2.2 模型因子获取及处理

(1) 土地利用与水土保持措施解译。根据遥感影像的色彩、形状、纹理等特征，结合实际进行人机交互解译，获取矢量及栅格数据。

(2) 植被覆盖度计算。借助MODIS-NDVI数据初步计算3 a平均24个半月植被覆盖度FVC数据，进行投影转换、裁剪和重采样等处理后，借助参数修订法计算，最终获取10 m分辨率栅格数据。

(3) 风蚀模型因子。依据兰考县风蚀定位观测站及周边县市2013—2019年风速数据，按1 m/s间隔统计≥临界风速的各等级风速累积时间，获取风力因子U数据。将SMAP L3数据进行投影转换、裁剪等预处理后相同半月取均值，提取表土湿度因子W数据。以土地利用数据为基础，结合兰考县风蚀观测站定位观测数据与野外调查成果，对一年两熟、一年一熟制耕地进行粗糙度赋值，生成地表粗糙度因子Z0数据。各因子数据为10 m分辨率栅格数据，以半月为1期，共24期。

(4) 水蚀模型因子。降雨侵蚀力因子R、土壤可蚀性因子K数据由淮河流域水土保持监测中心提供，经重采样、投影转换等处理后可用于模型计算。R因子数据全年共24期；以DEM数据为基础，运用北京师范大学开发工具初步提取坡长因子LS数据，对林草地LS修正后获取最终数据。依据《技术规定》初步计算园、林、草地的植被覆盖与生物措施因子B数据，其余土地利用类型直接赋值，经栅格计算、镶嵌等处理后生成B因子数据。依据土地利用数据和《技术规定》中要求，赋值生成水土保持工程措施因子E、耕作措施因子T数据。各因子为10 m分辨率栅格数据，其中R，B以半月为1期，共24期。

2.2.3 侵蚀计算及强度评价 借助ArcGIS的Raster Calcularor功能，叠加各因子数据进行模型计算(风蚀模型需将3种模型计算结果进行镶嵌)，最终获取兰考县2019年24个半月及全年的风蚀模数、水蚀模数。将24个半月及全年累计的风蚀、水蚀模数按《土壤侵蚀分类分级标准》(SL190-2007)[16]判断侵蚀强度(当区域同时存在风蚀和水蚀时，比较各栅格水蚀、风蚀强度，保留强度高的侵蚀类型；水蚀、风蚀强度相同时，保留水蚀)，得24个半月及全年土壤侵蚀数据和空间分布图，并在Origin，ArcGIS中进行数据分析。

2.3 各半月土壤侵蚀比例

为分析兰考县2019年24个半月土壤侵蚀变化特征，本研究参考植被覆盖与生物措施因子B相关研究中的WRi的计算方法，计算WAi和WQi，以明确区域风蚀、水蚀较强时间。数据结果在Excel和Origin中进行分析，计算公式如下：

(5)

(6)

式中：WAi为某半月侵蚀(风蚀、水蚀)模数占全年累计模数比例(%);Ai为某半月侵蚀(风蚀、水蚀)模数(t/hm2·a); WQi为某半月水土流失面积占全年累计面积比例(%);Qi为某半月水土流失面积(km2)。

3 结果与分析

3.1 侵蚀影响因子时间变化特征

气候、地形和植被等因素对土壤侵蚀发生、发展有重要影响，地形和土壤等因素短时间内不会变化，而气候、植被一年内动态变化显著。在遥感监测中，风力因子、降雨侵蚀力因子等数据是侵蚀影响因素的直观反映，侵蚀影响因子值的变化影响土壤侵蚀模数和强度。若研究区域土壤侵蚀特征，应先分析侵蚀影响因子时间变化特征。

风蚀模型和CSLE模型中，U，W，R及林园草的B因子数据以半月为单元，共24期，且存在各半月的动态变化。K，L，S因子短时间内变化的可能性较小，其余因子则依据土地利用解译结果赋值生成。因此，本研究重点关注U，W，R及林园草的B。在ArcGIS中获取其各半月的数值情况(图2—3)，分析其变化特征。

(1) 风蚀模型因子。分析图2可知，第13—18个半月(7—9月)表土湿度值较高；第1—4，23—24个半月(1—2月、12月)值较低。第15个半月表土湿度值最高，变化范围为0.19～0.28，均值为0.23；第3个半月值最低。同样由图2知，第3—9，21—22个半月(2—5月、11月)风速≥5 m/s累积时间较高，第14—18个半月(7—9月)时间较低。第6个半月风速累积时间处于1 442.96～1 549.52 min的范围内，均值为1 485.83 min，高于其他时间，第16个半月值最低。

图2 兰考县2019年24个半月风蚀模型因子

风速≥5 m/s累积时间全年呈“M”型变化；W值呈“波浪型”变化，风蚀模型因子变化与黄泛区季风性气候相关。2—5月及11—12月，区域干旱多风，此时风时间较长、土壤含水量较低，造成U值较高，W值较低。7—9月区域降水增加，且风力弱、风速小，还存在农田灌溉，此时起风时间短，土壤含水量升高，导致W值较高，U值较小。

(2) 水蚀模型因子。由图3可知，第9—16个半月(5—8月)林园草的B平均值较高；第1—3，23—24个半月(1—2月、12月)值较低。第13个林园草B的均值最高，为0.005；第24个半月值最低。同样由图3a知，第12—18个半月(6—9月)R值较高，第1—7，10—24个半月(1—4月、10—12月)值较低。第14个半月R值最高，变化范围为522.51～562.81 MJ·mm·(hm2·h·a)，均值为538.42 MJ·mm·(hm2·h·a)；第24个半月R值最低。

注：SLRi为第i个半月园地、林地和草地的土壤流失比例； WRi为第i个半月降雨侵蚀力占全年侵蚀比例。

R值和林园草B平均值呈“先升高再降低”的趋势，其变化与黄泛区气候及植被特点有关。1—4月及10—12月，区域降水较少，植被处于生长阶段，此时降雨造成的侵蚀较弱，植被覆盖度较低，造成R值和林园草B的均值较低。6—9月降水增加，植物进入茂盛期，此时植被覆盖度增加，且降雨造成侵蚀较强，导致R值和林园草B的均值较高。

3.2 兰考县2019年风蚀模数特征

侵蚀模数是模型计算的原始数据，其时间和空间分布特征可直观反映区域风蚀强弱。

(1) 时间特征。分析2019年兰考县各半月侵蚀模数情况(图4)知，第4—6，9，21—22个半月(2—3月、5月和11月)模数较高；第13—18个半月(7—9月)较低；模数均值、最大值呈“M型”变化。风蚀侵蚀模数全年累计值的变化范围为0～5 186.31 t/(km2·a)，均值为153.26 t/(km2·a)。其中第9个半月侵蚀模数均值最高，为18.36 t/(km2·a)；第6个半月模数处于0～659.61 t/(km2·a)范围内，存在全年最大值；第17个半月模数均值最低，为0.91 t/(km2·a)。

风蚀模数的变化与风速累计时间、表土湿度、地表粗糙度相关。第4—9，21—22个半月(2—5月和10—12月)，黄泛区风速大、风力强，且降水较少，造成风速累积时间长、表土湿度低；此外，此时耕地存在休耕、翻耕，部分区域粗糙度较小。因此，此时间段内侵蚀模数模数越高，风蚀较强。

图4 兰考县2019年24个半月风蚀模数

分析图5可知，风蚀全年WAi值处于0%～59.12%范围内，第3—6，9—12，19—22个半月(2—3月、5—6月及10—11月)的时间段内WAi平均值相对较高；第13—18个半月(7—9月)较低。其中第9个半月WAi均值、最大值均高于其他时间；第13个半月WAi均值最低；第14，18个半月WAi最大值低于其他时间。

图5 兰考县2019年风蚀24个半月WAi

WAi最大值、均值变化趋势与风蚀模数相似，呈“M型”，证明了2—3月、5月和11月兰考县风蚀较强，7—9月风蚀较弱。

(2) 空间特征。由图6可知，东坝头乡、许河乡、红庙镇、南彰镇存在侵蚀模数极高区域，仪封乡、张君墓、三义寨乡镇整体模数较高；谷营镇、城关镇、孟寨乡整体模数较低。综合土地利用分析[10]可知，风蚀多发生在水浇地、采矿用地和沙地上；兰考县水浇地约占全县面积的60.78%[10]，且部分区域冬季休耕，使水浇地易出现风蚀；而采矿用地、沙地周围堆积的沙土易被强风吹蚀搬运，导致范围内风蚀较强。城镇/农村建设用地、道路及水域不参与模型计算，风蚀模数为0；园地、林地和草地植被覆盖相对较高，模数不为0但数值偏小，不易出现风蚀。

图6 兰考县2019年风蚀侵蚀模数空间分布

3.3 兰考县2019年水蚀模数特征

(1) 时间特征。由图7可知，第12—18个半月(6—9月)模数较高，水蚀相对较强；第1—7，19—24个半月(1—4月、10—12月)较低，水蚀较弱；模数均值、最大值呈“升高—降低”的变化。水蚀模数时间特征与风蚀相反。模数全年累计值处于0～8 028.86 t/(km2·a)范围内，均值为9.54 t/(km2·a)。第14个半月模数最高，最高值为1 518.16 t/(km2·a)，均值为1.82 t/(km2·a)；第3个半月模数最低，均值为0.01 t/(km2·a)。

分析图8可知，水蚀全年WAi值处于0%～77.46%范围内，第9—18个半月(5—9月)WAi值较高；其余时间较低。第13个半月WAi最值高于其他时间，最大值为77.46%；第14个半月WAi均值最高，为20.34%。WAi最大值、均值变化趋势与水蚀模数相似，证明6—9月水蚀相对较强，其余时间水蚀较弱。

(2) 空间特征。由图9可知，东坝头乡、城关镇、仪封乡存在侵蚀模数较高区域；爪营乡、张君墓镇、孟寨乡整体模数较低。黄泛区非典型水蚀区，水蚀与自然条件相关。兰考县地势平坦，地形起伏较小，土层深厚，夏季降水多，生产建设活动频繁，故模数较高的区域多为采矿用地。

图7 兰考县2019年24个半月水蚀模数

图8 兰考县2019年24个半月水蚀WAi

图9 兰考县2019年水蚀侵蚀模数空间分布特征

3.4 兰考县2019年土壤侵蚀特征

兰考县同时存在风蚀和水蚀，研究区域土壤侵蚀特征时应依据相关要求和技术规范，叠加风蚀、水蚀后进行综合分析，最终得兰考县2019年半月及全年土壤侵蚀特征。

(1) 半月特征。分析兰考县2019年内各半月水土流失情况(图10)可知，除第2，10—11，19—20，23—24个半月外、其余时间均产生水土流失，出现水土流失的半月数占全年的70.83%。第4—9个半月(2—5月)水土流失面积较多，面积维持8.61 km2不变，且均为风蚀；第12—18个半月(6—9月)水土流失较少。产生水土流失的时间内，WQi值变化范围为0.01%～3.18%。上述数据表明第4—9，21—22个半月(2—5月、11月)是兰考县土壤侵蚀较强时间。

图10 兰考县2019年24个半月水土流失面积

(2) 全年特征。由表3可知，兰考县2019年水土流失面积共271.66 km2，占全县面积的24.34%，主要由风蚀造成。其中，风蚀造成水土流失270.04 km2，占总面积的99.40%；水蚀面积造成水土流失仅1.62 km2。风蚀以轻度侵蚀为主，面积261.47 km2，占风蚀水土流失面积的96.83%；其次为中度侵蚀和强烈侵蚀。水蚀以轻度侵蚀为主，面积1.58 km2，占水蚀水土流失面积的97.53%；中度侵蚀最少，仅0.04 km2。

表3 兰考县2019年水土流失面积及比例

分析兰考县2019年土壤侵蚀图(图11)可知，轻度风蚀在全县广泛分布，仪封乡、张君墓镇、谷营镇和东坝头乡等面积相对较多；孟寨乡、南彰镇分布相对较少；红庙镇、许河乡、东坝头乡和城关镇存在一定面积中度风蚀。水蚀面积较少，分布较为分散，在仪封乡相对较多。

注：中度水蚀和强烈风蚀的斑块分散、微小，正常比例尺下较难分辨，故图中没有表达。

4 讨论与结论

(1) 兰考县风力因子半月变化呈“M”型，第3—9，21—22个半月值较高；表土湿度因子变化呈“波浪型”，第13—18个半月值较高。降雨侵蚀力因子、林园草植被覆盖与生物措施因子B均值的变化为“先升高再降低”，前者第9—16个半月值较高，后者第12—18个半月值较高。

(2) 兰考县半月风蚀模数变化呈“波浪型”；全年累计风蚀模数在0～5 186.31 t/(km2·a)范围内，均值为153.26 t/(km2·a)；第4—9，21—22个半月风蚀较强；仪封乡、张君墓镇等风蚀模数较高；风蚀多在水浇地、采矿用地上。与前人相似研究比较[7,10]，风蚀空间分布特征相近。

(3) 兰考县半月水蚀模数变化趋势为“先升高再降低”；全年累计水蚀模数在0～8 028.86 t/(km2·a)范围内，均值为9.54 t/(km2·a)；第12—18个半月水蚀相对较强；东坝头乡、城关镇存在水蚀模数较高区域；水蚀多在采矿用地上。

(4) 2019年兰考县共17个半月出现水土流失，占全年的70.83%；第4—9个半月水土流失面积最多；均为风蚀。全年水土流失面积为271.66 km2，占全县面积的24.34%；主要为轻度风蚀，占水土流失面积的96.83%。风蚀主要在仪封乡、张君墓镇、谷营镇和东坝头乡，水蚀主要在仪封乡。水土流失监测与防治应重点关注此时间内水浇地、采矿用地上的侵蚀。本研究所得兰考县水土流失面积与前人相似的研究[7,10]比较有所差异，由研究区域差异导致；但风蚀空间分布特征相似，证明了本研究结果的准确性、合理性。