延河流域位于东经108°45'—110°28',北纬36°23'—37°17',面积为7 725 km²,平均坡度4.4%,河网长度4.7 km/km²。该区域属于暖温带大陆性半干旱季风气候,年均降水量为500 mm,年均气温为8.8～10.2℃。土壤类型主要以黄绵土为主,占流域面积80%以上,该土体疏松、软绵,抗侵蚀能力差,因此,延河流域是我国水土流失的重要区域之一。研究区地势西北高东南低,地貌类型主要包括,残原平梁沟壑区、黄土丘陵沟壑、黄土台地、以及黄土覆盖的山地,其中分布最为广泛的是丘陵沟壑区,约占研究区面积的90%。由于流域内水土流失严重,近年来采取了一系列的水土保持措施及植被覆盖措施,使得流域生态状况逐渐好转,图1为2003年、2014年的年NDVI空间分布特征。

图1 延河流域2003年和2014年NDVI空间分布

本研究所使用的数据:(1)研究区DEM是通过选用Sentinel-1A数据为数据源,利用INSAR技术获取的,分辨率为16.87 m,通过与水经注获取的19.11 m的DEM进行对比验证,发现该DEM精度较高。(2)土地利用数据来源于中国科学院资源科学环境数据中心,分辨率为30 m,研究中主要选用2005年、2010年土地利用数据。(3)土壤数据来源于世界土壤数据库(Harmonized World Soil Database,HWSD),通过土壤参数计算软件(Soil Parameter Calculation Software,SPCS)获取土壤参数数据。(4)气象数据来源于中国气象数据共享网,包括延河流域内及周边站点的1999—2014年的逐日降雨、温度、相对湿度和风速数据。(5)水文数据来源于《中华人民共和国水文年鉴》,包括甘谷驿水文站1999—2014年的逐日径流数据和1999—2014年逐日产沙数据。(5)梯田、水库数据来自于Sentinel-2A数据,通过解译影像得到2015年研究区梯田、水库分布情况。(6)归一化植被指数(Normalized Difference Vegetation Index,NDVI)数据来自MOD13Q1产品,该产品具有250 m的空间分辨率,每隔16 d提供一次,通过合成处理,得到2003—2014年NDVI数据。

2.1.1 SWAT模型土壤侵蚀子模块以及土壤保持量的获取 目前已有大量研究利用SWAT水文模型对黄土高原区域的产流产沙进行模拟^[11-13],因此本研究选用SWAT模型,对延河流域土壤保持量进行模拟。该模型中产沙量的计算采用改进版通用水土流失方程(Modified Universal Soil Loss Equation,MUSLE),如式(1)所示。

Sed=11.8·(Q_surf·q_peak·area_hru)^0.56·K_usle·C_usle·P_usle·LS_usle·CFRG(1)

式中:Sed为模拟产沙量; Q_surf为地表径流(mm H₂O/hm²); q_peak为径流峰值(m³/s); area_hru为HRU面积(hm²); C_usle为地表植被覆盖因子; P_usle为水土保持措施因子; K_usle为土壤可蚀性因子; LS_usle为地形因子,包括坡长因子L和坡度因子S; CFRG为直径大于2 mm的粗碎块因子。

(1)土壤可蚀性因子(K_usle)可以表示土壤对侵蚀的敏感性,SWAT中计算K_usle因子的方法,具体公式如下式所示。

K_usle=f_csand·f_cl-si·f_orgc·f_hisand(2)

式中:f_csand为沙土质地土壤侵蚀因子; f_cl-si为黏壤土土壤侵蚀因子; f_orgc为土壤有机质因子; f_hisand为高沙质土壤侵蚀因子。

式中:m_s,m_silt,m_c分别为砂粒、粉粒、黏粒的含量(%); orgC为有机碳含量(%)。

(2)植被覆盖因子(C_usle)是影响土壤侵蚀的最敏感因子,主要受到覆盖度、植被类型等的影响,SWAT模型中计算公式如下。

C_usle=exp([ln(0.8)-ln(C_usle,mn)]· exp[-0.00115·rsd_surf]+ln[C_usle,mn])(7)

式中:C_usle,mn为土地利用的最小植被覆盖度因子; rsd_surf为土壤表面残留物(kg/hm²)。

C_usle,mn因子可由已知的C_usle,aa因子用下式估计:

C_usle,mn=1.463ln[C_usle,aa]+0.1034(8)

式中:C_usle,mn为土地利用的最小植被覆盖度因子; C_usle,aa为土地利用年均植被覆盖因子。

(3)水土保持措施因子(P_usle)是指采取水保措施后,土壤流失量相对于未采取任何措施时的比例,取值范围为0～1^[14]。0表示不存在侵蚀现象的区域,1表示未采取水土保持措施的区域。SWAT模型中P_usle因子是采用Wischmeier和Smith给出的信息计算,为了使模拟结果更适合研究区,因此本文对P_usle因子进行了修改。SWAT模型中可对每一个HRU的P_usle因子进行修改,HRU中包含土地利用、坡度等信息,因此本文根据土地利用类型及坡度范围,并考虑前人在黄土高原的研究以及黄土高原上水保措施的实施,包括梯田、水平沟以及鱼鳞坑等的影响^[15-16],最终确定表1中P_usle因子的取值。

表1 水土保持措施因子(Pusle)的取值

(4)坡度坡长因子(LS)可以表示地形因子对土壤侵蚀的影响,SWAT模型中采用如下公式进行计算。

式中:L_hill为坡长(m); m为指数项; α_hill为坡度。

m的计算公式如下:

m=0.6(1-exp[-35.835·slp])(10)

式中:slp为HRU的坡度。

slp与α_hill的关系为:

slp=tanα_hill(11)

研究中涉及潜在土壤侵蚀量(Sed_q)、实际土壤侵蚀量(Sed_a)、土壤保持量(Sed_c)3种数据量的计算。潜在土壤侵蚀量是指在假定没有地表植被覆盖和水土保持措施下的侵蚀量,即通过将SWAT源代码和数据库中的C_usle,P_usle因子修改为1得到的产沙量。实际土壤侵蚀量是指流域实际发生的侵蚀量,即由SWAT模型直接模拟得到的产沙量。土壤保持量是二者的差值,计算公式如下式所示。

Sed_c=Sed_q-Sed_a=11.8(Q_surf·q_peak·area_hru)^0.56·

K_usle·LS_usle·CFRG·(1-C_usle·P_usle)(12)

2.1.2 数据库的构建

(1)土地利用数据库。研究中主要使用2005年、2010年的土地利用数据,将延河流域土地利用数据分为耕地、林地、园地、高覆盖草地、低覆盖草地、建设用地、水域和未利用地8类,模型在输入土地利用数据时,同时还需要输入土地利用索引表,方便将流域的土地利用类型与SWAT模型自带的植被生长模型库和农业管理数据库进行连接。因此,按照SWAT模型中土地利用的编码对各种土地利用重新进行编码,见表2,从而建立土地利用索引表。

表2 土地利用类型SWAT编码

(2)土壤数据库。本文所使用的土壤数据来源于HWSD,流域的土壤类型主要包括黄绵土、冲积土、粗骨土、红黏土、潜育始成土等,见图2。通过在HWSD中获取各种土壤的部分参数以及利用SPCS软件计算部分土壤参数,建立土壤参数数据库,并根据各类土壤在SWAT中的编码,建立土壤参数索引表。

图2 土壤类型

(3)气象数据的准备。本文选用流域内及距离流域较近的11个传统气象站点1999—2014年的逐日气象资料输入到SWAT模型中,其中降水、最高温度、最低温度、相对湿度、风速可直接在中国气象数据共享网站获取,太阳辐射数据通过天气发生器模拟得到。

2.1.3 模型的建立 本研究使用ArcSWAT 2012版本,首先是子流域的划分,基于DEM生成河网,通过设定阈值(形成子流域的最小给水面积)为10 000 hm²,设置流域出水口点,将延河流域划分为41个子流域; 其次是水文响应单元(HRU)的划分,模型在划分子流域的基础上,通过叠加土地利用数据、土壤数据和坡度数据,获得各个子流域中具有相同土地利用类型、土壤类型和坡度的HRU。最后将气象数据等导入模型,对延河流域1999—2014年产沙量进行模拟。模型以1999—2002年作为预热期,利用2003—2008年甘谷驿水文实测值进行率定,利用2009—2014年的水文实测值进行验证。由于研究时间较长,土地利用会发生一定程度地变化,因此将预热期和率定期视为整体,选用2005年土地利用进行模拟; 在验证期时选用2010年土地利用进行模拟。

2.1.4 模型的参数敏感性分析、率定、验证以及不确定性分析 本文通过选用SWAT-CUP(SWAT-Calibration and Uncertainty Programs)中的SUFI-2算法对SWAT模型进行敏感性参数分析、率定、验证以及不确定分析。在该算法中,采用判定系数(Cofficient of determination,R²)和纳什效率(Nash-Sutcliffe efficiency coefficient,NSE)对模型进行适用性评价,其中R²可评价模拟值与实测值变化趋势的一致性,NSE可反映模拟值与实测值之间的拟合程度^[17]。

式中:Q_mi为模拟值; Q_oi为实测值; Q_m为模拟均值; Q_o为实测均值。

本文选用甘谷驿水文站1999—2014年的逐月径流数据和泥沙数据作为实测值,由于未能获取枯水期时泥沙数据,因此选用0值代替。在SWAT-CUP中进行敏感性参数分析后,完成对模拟结果的率定和验证。图3为研究区径流和泥沙实测值与模拟值的比较结果。一般研究认为,R²>0.6,NSE>0.5时,模拟结果较可靠,校准期和验证期模拟结果都较好时,参数的最佳校准值适用于研究区域^[18]。

从图3(A1—A2)可以看出,径流量在率定期(2003—2008年)和验证期(2009—2014年)模拟结果的R²和NSE均达到模型要求; 从图3(B1—B2)可以看出,泥沙量在率定期(2003—2008年)和验证期(2009—2014年)模拟结果的R²和NSE也均达到模型要求。因此,SWAT模型可以对延河流域的水文状况进行较好地模拟。从而可得到延河流域基于HRU的2003—2014年逐月产沙量。

图3 延河流域甘谷驿水文站点月径流量和月产沙量在率定期、验证期的模拟结果

3.2.1 时间尺度变化 从图4可以看出,在时间尺度上,潜在土壤侵蚀量、实际土壤侵蚀量和土壤保持量变化趋势相同,均会产生不同幅度的波动,主要原因是受降水因子的影响,各个月份的降水量、降水强度均不同,其中在2013年7月份3种数据量均产生了大幅度增加且均达到多年最大值,主要原因是该月降雨总量为435 mm,远远超过其他月份的降水量。因此,可以发现延河流域土壤保持在时间尺度上的变化主要受降水因子的影响。

延河流域土壤保持量在2003—2014年,整体上呈现出波动的趋势,但土壤保持量仍有所增加,其中2003年土壤保持总量为677.1万t,2014年土壤保持总量为1 345.02万t。二者的空间分布格局也产生了很大的变化。由图5可以看出,2003年土壤保持量主要集中在流域下游,主要原因是该下游区域,地势平坦,侵蚀轻微,因此产生的土壤保持量较高,而在2014年土壤保持量在整个流域均有分布,主要原因是该区域近年来实行了一系列的水保措施及植被覆盖措施产生了良好的生态效应。

图4 潜在土壤侵蚀量、实际土壤侵蚀量和土壤保持量随时间变化趋势

为了更加明确土壤保持服务在不同年份空间分布差异性的原因,本文尝试分析研究区水土保持措施及植被覆盖措施的变化及其对水土保持服务的影响。因此,本文通过解译高分辨率影像得到研究区梯田、水库分布情况,并用此来表示研究区水土保持措施实施的状况。由于梯田数据的提取需要高分辨率影像的支持,因此选用10 m分辨率的Sentinel-2A数据,但该数据是最早发行时间为2015年,因此通过解译影像得到研究区2015年梯田、水库数据,通过观察可知,经过多年来对坡耕地的改造,到2015年为止流域上中游梯田及水库分布较广泛,下游分布相对较少。本文通过处理MOD13Q1数据得到的研究区2003年、2014的年NDVI数据,并用此来表示研究区植被覆盖状况的变化。可以看出,经过多年来实施的退耕还林(草)工程,整个流域的植被覆盖状况发生了明显地好转,尤其是流域的上中游区域。由此可知,梯田、水库和NDVI数据在空间上的变化可以直接证明2003年和2014年土壤保持量在空间上的差异,也进一步说明水土保持措施及植被覆盖措施在流域治理过程中发挥了至关重要的作用。

图5 延河流域2003年和2014年土壤保持量

3.2.2 空间尺度变化 从图6可以看出,在空间尺度上,土壤保持量主要受地形因子的影响。2003—2014年延河流域年均土壤保持量为1 002.02万t,在流域上中下游区域表现不同,其中在下游年平均土壤保持量最大,上游次之,中游最低,造成这种现象的主要原因是,流域下游的区域主要为残原平梁沟壑区,该区域黄土侵蚀方式以雨滴击溅侵蚀为主,侵蚀轻微,且地面平坦、土层肥厚^[19],因此有着最高的土壤保持量; 流域上游主要为黄土覆盖的山地区域,该区域25°以上的陡坡占总面积的70%以上^[19],虽易受土壤侵蚀的影响,导致该区域的潜在土壤侵蚀量大于其他区域,但由于近年来这个区域实施退耕还林力度最大,使得25°以上的区域全部退耕,进而植被不断恢复,植被覆盖度不断增大,该区域实际土壤侵蚀量随之不断减少,所以该地区的土壤保持量相对中游较高; 流域中游主要为梁峁丘陵沟壑区,该地区人类活动频繁,垦殖严重,因此该区域的土壤保持量最低。

为了更清楚延河流域土壤保持量受地形因子的影响,将多年平均土壤保持量进行标准化,参照蒋春丽等人的研究^[20],将标准化后的土壤保持量分为低(0～0.2)、较低(0.2～0.4)、中(0.4～0.6)、较高(0.6～0.8)、高(0.8～1.0)5个等级,见表3。延河流域主要以低保持为主,面积为5 432.94 km²,占总面积的73%。

图6 延河流域多年(2003-2014)年平均土壤保持量

表3 土壤保持量分级统计

3.3.1 坡度与土壤保持分布特征 经统计得到各土壤保持量在不同坡度等级所占面积比例(表4)。由表可知,各级土壤保持量随坡度的增大均呈先增加后减小的趋势,其中在<25°区域,低保持占比较多。主要是因为耕地在这一区域分布范围较广^[21],人类活动比较频繁,易发生土壤侵蚀。>25°区域,低保持占比较小,较低、中和高保持量占比较多。这主要是因为延河流域自退耕还林以来,坡度是执行退耕的首要标准,>25°以上的坡耕地必须退耕^[21],这一措施导致>25°的区域耕地面积逐渐下降,草地、林地面积相应增加,中度和强度以下侵蚀的面积占比显著增加,剧烈侵蚀面积显著下降。而>35°的区域各级保持量均低于25°～35°区域,主要是因为大于35°属于急陡坡,其耕作类型虽发生改变,但该区域地形陡峭,易受土壤侵蚀的影响。

总之,流域土壤保持随坡度变化规律明显,25°的坡面是土壤保持量的关键带,也是今后水土保持工作应重点关注的区域。

3.3.2 坡向与土壤保持分布特征 通过将土壤保持量分级图与坡向进行叠加分析,得到土壤保持量空间分布的坡向特征图(图7)。由此可知,在平地上,基本只存在低、较低和中保持量,其他各级保持量几乎为0。经统计可得,土壤保持量在不同坡向上的规律表现为半阳坡>阴坡>阳坡>半阴坡>平地。这主要是因为对于黄土丘陵沟壑区,平地上人类活动频繁,因此基本只存在低、较低和中保持量,且占比都较小; 阳坡所获得的太阳辐射总量高于阴坡,但阳坡土壤湿度低于阴坡,且阳坡上人类耕作现象较多^[22-24],并参照汤巧英等^[24]的研究成果,可知延河流域2000年、2010年植被覆盖度总体表现为阴坡>阳坡>平地,因此阴坡上的土壤保持量高于阳坡。又根据李勉等人^[25]的研究成果可知,黄土高原不同坡向土壤侵蚀速率大小表现为半阴坡>阳坡>半阳坡,从这一规律可以看出半阳坡的土壤保持量应相对较高。由此可知,土壤保持量在不同坡向的分布与坡向的植被覆盖、土壤湿度、侵蚀速率等有着紧密的联系。

表4 不同坡度土壤保持面积分布情况

总之,流域土壤保持随坡向变化规律明显,主要是与坡向上的植被覆盖、湿度、光照、侵蚀速率等因子关系密切。半阳坡为土壤保持量最大的区域,半阴坡为土壤保持量较小的区域,二者均为日后开展水保工作的重点对象。

图7 土壤保持量坡向分布特征

3.3.3 高程与土壤保持分布特征 将流域高程分级图与土壤保持量分级图进行叠加分析,得到土壤保持量空间分布的高程特征图(图8)。由此可知,随着高程的增加,土壤保持量大致呈先增大后减小的现象,土壤保持量在低、较低、中和较高保持的高程分布特征大致相同,其峰值均在1 250～1 300 m高程范围内,土壤高保持峰值主要分布在950～1 000 m的高程范围内,主要原因是,在1 250～1 300 m高程范围内,主要为黄土丘陵沟壑区,是延河流域典型的地貌类型,该种地形上虽易发生土壤侵蚀^[26],但近年来在丘陵沟壑区开展了淤地坝工程,给该地区带来了良好的生态和经济效益,使得该地区的土壤保持量不断增加; 在950～1 000 m的高程范围内,一般为黄土台地,土壤肥厚,地面平坦,不易发生土壤侵蚀现象,因此存在高土壤保持量。

图8 土壤保持量高程分布特征

总之,流域土壤保持随高程变化分异特征明显,950～1 000 m和1 250～1 300 m是土壤保持量最多的地带,也应是今后水土保持的关注重点。