黑河中下游防风固沙功能时空变化及影响因子分析

彭婉月, 王兆云, 李海东, 柳本立

1.中国科学院西北生态环境资源研究院沙漠与沙漠化重点实验室, 甘肃 兰州 730000

2.甘肃省戈壁荒漠生态与环境野外科学观测研究站(敦煌)， 甘肃 兰州 730000

3.中国科学院大学， 北京 100049

4.生态环境部南京环境科学研究所， 江苏 南京 210042

黑河中下游是国家级防风固沙型生态功能保护区[1]，也是西北生态环境最为脆弱的地区之一. 黑河中游为人工绿洲区，是流域水资源主要消耗区[2]；下游沙漠、戈壁广布，水质恶化、沙尘暴频繁、土地荒漠化等生态环境问题频发，是我国北方主要的沙尘发源地之一[3]. 2000年以来，国家相继实施了黑河流域应急调水措施、黑河流域综合治理工程、全国第一个节水型社会试点及甘肃张掖黑河湿地国家级自然保护区等建设，以期强化防风固沙功能，改善流域的生态环境.

防风固沙是生态系统通过其结构与过程减少由于风蚀所导致的土壤侵蚀的作用[4]，这种保持土壤、抑制风蚀过程的功能即为防风固沙功能[5]. 开展黑河中下游近年来的防风固沙功能评价工作，有利于遏制沙化趋势、维护该流域的生态安全和科学开展土壤风蚀防治. 近年来， 黑河中下游防风固沙研究主要集中于小区域的生态效益[6-8]，缺少对中下游全区的防风固沙功能评估及长时间序列的空间差异性分析.

目前有多种风蚀模型可用于防风固沙功能的评价. 自20世纪60年代起，国外先后开发了风蚀方程(wind erosion equation, WEQ)[9]、修正风蚀方程(revised wind erosion equation, RWEQ)[10]、德克萨斯模型(texas erosion analysis model, TEAM)[11]、风蚀评价模型(wind erosion assessment model, WEAM)[12]、风蚀预报系统(wind erosion prediction system, WEPS)[13]等. 我国防风固沙评估的模型研究进展较为缓慢，在1998年才提出首个风蚀量统计模型，后基本以田间尺度的经验估算模型为主[14]. 目前，使用RWEQ进行评估的方法应用最为广泛，在内蒙古自治区锡林郭勒[15-16]、青海省[17]、浑善达克沙地[18]等不同区域尺度都开展了相关工作. 该方法也被作为生态保护红线划定中防风固沙功能重要性的评估方法[4]. 针对黑河流域的研究中，韩永伟等[5]利用风蚀流失模型评估了2006年黑河下游防风固沙功能，但该方法与RWEQ模型有显著不同，也未进行因子分析.

该文针对黑河中下游2000—2017年间连续多年和多年年均防风固沙功能进行评估，利用近20年最新数据分析变化趋势，并确定主要影响因子的贡献，可为区域的风沙灾害防治、生态修复及贫困地区的生态扶贫提供参考依据，提升流域内居民的生存环境和毗邻地区的环境质量，进而对国家生态文明和“美丽中国”建设具有重要意义.

1 材料与方法

1.1 研究区概况

黑河是我国西北地区的第二大内陆河，处于西北干旱区西风季风交汇地带，也是西北干旱区最具代表性的河流，发源于祁连山北麓中段，流经河西走廊，最终注入内蒙古自治区额济纳旗的东西居延海，其干流全长821 km. 黑河流域的经纬度范围介于98°E～102°E、38°N～42°N之间，流域总面积1.3×105km2，以鹰落峡、正义峡为分界线将全流域划分为上、中、下游. 黑河流域跨青海省、甘肃省和内蒙古自治区的5地(州)、11县(市、旗)，人口约1.3×106人. 上游地区以牧业为主，中游属灌溉农业经济区，所占人口比重最大，达1.21×106人；下游地区由于其地质特征，既有农田、林草灌溉区，又发展荒漠牧业. 黑河中下游地区包括甘肃省的山丹、民乐、张掖、临泽、高台、金塔县部分地区等县(市)，以及内蒙古自治区额济纳旗(见图1). 流域中下游属于蒙-甘气候区，分为温带河西走廊干旱亚区和额济纳温带极端干旱亚区. 中游年均降水量100～250 mm，年均蒸发量 2 047～2 341 mm，干燥度2.5～5.5，日照时数 2 683～3 088 h，基本特征为干旱多风、光照充足，气温日较差大；下游以西风和西北风为主，年均降水量41 mm，年均蒸发量 3 700～4 384 mm，年均气温8.3 ℃，年均风速4.2 ms，年均沙尘暴日数20 d，其气候特征为极端干旱、大风、少雨、沙尘频繁[2]，沙漠化敏感性和盐渍化敏感性高，防风固沙功能极重要[1].

注： 数据来源于赵军,王建华.国家青藏高原科学数据中心(http:data.tpdc.ac.cnzh-hans)，2015.

1.2 评价方法与数据来源

1.2.1RWEQ模型

防风固沙功能主要与风速[19]、土壤[20-22]、地形[23]和植被[24-25]等因素密切相关. 根据生态环境部生态保护红线划定指南[4]和有关文献[26-28]，该研究采用修正风蚀方程(RWEQ). 计算公式：

SR=SLP-SL

(1)

(2)

S=150.71×(WF×EF×SCF×K′×C)-0.371 1

(3)

Qmax=109.8×(WF×EF×SCF×K′×C)

(4)

(5)

Qmaxp=109.8×(WF×EF×SCF×K′)

(6)

SP=150.71×(WF×EF×SCF×K′)-0.371 1

(7)

各主要因子的计算方法如下：

a) 气候因子(WF).

WF=Wf×SW×SD×ρg

(8)

(9)

式中：Wf为各月多年平均气候风蚀因子，m3s3；SW为各月多年平均土壤湿度因子；SD为雪盖因子(无积雪覆盖天数研究总天数)，定义积雪覆盖深度大于25.4 mm视为积雪覆盖；ρ为空气密度，当气温为15 ℃时为1.226 kgm3；g为重力加速度，取9.8 ms2；U2为距离地面2 m处的风速，ms；UC为临界风速，一般设为5 ms[16,30]；Nd为每月天数，d；N为每月观测风速大于临界风速的天数，d.

b) 土壤可蚀因子(EF).

(10)

式中：sa为土壤粗砂含量，%；si为土壤粉砂含量，%；cl为土壤黏粒含量，%；[OM]为土壤有机质含量，%；[CaCO3]为碳酸钙含量，%，可不予考虑.

c) 土壤结皮因子(SCF).

d) 植被覆盖度因子(C).

C=e-ai×SC

(12)

式中：SC为植被覆盖度，%；ai为不同植被类型的系数，林地、草地、灌丛、裸地、沙地、农田分别为 0.153 5、0.115 1、0.092 1、0.076 8、0.065 8、0.043 8.

e) 地表糙度因子(K′).

K′=e(1.86Kr-2.41Kr0.934-0.127Crr)

(13)

Kr=0.2×ΔH2L

(14)

式中：Kr为土垄糙度，根据Smith-Carson方程计算，cm；Crr为随机糙度因子，取0 cm；L为地势起伏参数，m；ΔH为距离L范围内的海拔高程差，m.

在充分考虑气候条件、植被覆盖状况、土壤可蚀性、土壤结皮、地表粗糙度等要素情况下，以单位面积防风固沙量(潜在土壤风蚀模数与实际土壤风蚀模数的差值)作为表征防风固沙功能的评估指标，所以下文中防风固沙功能实际指代单位面积防风固沙量. 同时，为更好地表达空间差异，依据研究区防风固沙功能值域及参考文献[5,18]，设置防风固沙功能在0～2×104tkm2之间为较低区，在2×104～2.8×104tkm2之间为一般区，在2.8×104tkm2以上为较高区.

1.2.2一元线性回归斜率分析

该研究基于最小二乘拟合直线对数组进行回归分析，模拟并预测一组数据的时间变化趋势过程[31]，计算公式：

(15)

式中：slope为变量回归方程的系数，若slope>0，表示变量趋势增加，而slope<0，则表示变量趋势减少；n为年数，该研究取值为18；Yt为第t年的变量值.

通常需要对变化趋势进行显著性检验，该研究采用F检验，通过比较两组数据的方差，以确定它们的精密度是否具有显著性差异，其计算公式：

(16)

(17)

(18)

1.2.3灰色关联度分析法

灰色关联度分析法(grey relation analysis，GRA)[32-33]是灰色系统理论中进行多因素相关程度分析的方法，通过定量描述和比较的方法分析该系统的变化发展态势. 关联度是反映事物或者因素之间相关性大小的量度值，而关联系数是子序列与母序列在各时刻的关联程度值，使用关联度可计算各因子的贡献率. 采用初值化方法[27]归一化数据序列后参照式(19)～(21)计算.

(19)

(20)

(21)

式中，Em为各影响因子贡献率，r(x0,xm)为关联度，n′为被评价对象的个数，m为影响因素的指标数，k为不同的时刻数，n″为时刻总数，x0(k)为母序列，xm(k)为要分析的影响因素子序列，r[xo(k),xm(k)]为关联系数，ξ为分辨系数且ξ∈(0，1).

1.2.4数据来源

数据来源详细说明见表1. 其中，由2000—2017年共18年的日最大风速插值得到风力因子；雪盖因子利用中国雪深长时间序列数据集进行计算；土壤特性因子利用黑河流域土壤粒径分布数据集、面向陆面模拟的中国土壤数据集及中国土壤有机质数据集得到；用MRT投影转换工具对MOD13A3进行投影和格式转换批处理后，采用最大值合成法获得NDVI数据，并基于像元二分法生产出年植被覆盖度，进而得到植被覆盖因子. 地表糙度因子由GDEM DEM 30 m分辨率数字高程数据得到. 研究区域内共有4个地级市及其所辖的10个县级行政单元，但民乐县、临泽县、肃南裕固族自治县、嘉峪关市气象站的数据资料存在缺失，仅能使用张掖市、酒泉市、祁连县、山丹县、高台县、金塔县及额济纳旗7个气象站数据(见图1). 将所有因子统一为250 m分辨率的栅格数据，保证栅格计算过程准确.

表1 数据来源介绍

2 结果与讨论

2.1 防风固沙功能空间分布特征及趋势分析

由图2可见，研究区有1.02×105km2的面积风蚀严重，占黑河中下游地区面积的一半以上，主要发生在内蒙古自治区额济纳旗地区，与区域内绿洲、荒漠的分布特征一致. 防风固沙功能较高区主要是流域中下游绿洲、山地及黑河水系沿线部分，约占研究区面积的31.54%，一般区主要分布在较高区周边及内部，约占研究区面积的20.77%，另外，分布有47.69%的较低区，主要位于荒漠区及山区.

图2 黑河中下游多年平均土壤风蚀模数及防风固沙功能空间分布

具体说来，青海省祁连县山区防风固沙功能区间为0～2.9×104tkm2，96.82%的面积为较低区，3.17%的面积为一般区，而较高区的面积趋近于0. 中游甘肃省的防风固沙功能范围为0～3.1×105tkm2，34.70%的面积为较高区，14.92%的面积为一般区，50.38%的面积为较低区. 下游内蒙古自治区额济纳旗防风固沙功能范围为0～1.1×105tkm2，其中30.02%的面积为较高区，26.30%的面积为一般区，43.68%的面积为较低区. 研究区防风固沙功能整体呈现中游较强，向下游递减的空间分布特征，主要原因是中游绿洲分布集中且是灌溉农业区，额济纳旗有大片沙地、戈壁及裸岩，大风频繁；而青海省祁连县防风固沙功能最低，主要是因为其下垫面多为林地、草地及河渠，虽然防风固沙能力高，但其风沙活动少，实际土壤风蚀量小.

防风固沙功能的变化趋势如图3所示. 在研究区内，主要是甘肃省张掖市和嘉峪关市防风固沙功能趋势增加，回归方程系数(slope)为0～26.29%，占总面积的12.51%. 趋势减弱地区主要分布于内蒙古自治区额济纳旗东北部和甘肃省高台县中部，回归方程系数为-17.17%～0，占总面积的23.30%. 趋势增强区有较多自然植被及栽培植被且土质肥沃，趋势减弱区多是缺乏植被保护的荒漠风沙土. 总体来说，土壤类型、植被覆盖和气候条件三要素形成了防风固沙功能的空间分布模式.

图3 2000—2017年黑河中下游年防风固沙功能变化趋势

2.2 区域土壤风蚀和防风固沙功能年际变化特征

由图4可见，2000—2017年研究区土壤风蚀强度波动变化，平均土壤风蚀量约为4.6×109t，平均土壤风蚀模数为3.48×104tkm2. 其中2000—2001年风蚀最严重，尤以2001年土壤风蚀量达到评估期顶峰(9×109t)，土壤风蚀模数为6.84×104tkm2；而在2011年土壤风蚀量达到最小值(1.3×109t)，土壤风蚀模数为1×104tkm2. 土壤风蚀模数降低，代表整体土壤风蚀活动减弱，地表物质区域稳定，有利于开展生态工程的建设. 其中下游内蒙古自治区额济纳旗多年平均土壤风蚀量和多年平均土壤风蚀模数分别为3.1×109t和4.4×104tkm2，中游所含甘肃地区分别为1.5×109t和2.5×104tkm2，青海省祁连县山区分别为9×105t和368 tkm2(见图5). 土壤风蚀模数趋势减小，年土壤风蚀量平均减少1.36×107t，年均变化率为1.60%.

图5 代表性年份的黑河中下游土壤风蚀模数变化

图4 2000—2017年黑河中下游土壤风蚀模数变化

研究区年均防风固沙功能于2002年和2017年达到最低值和最高值，分别为8.2×103～3.7×104tkm2(见图6)，多年平均值为2.44×104tkm2. 其中，中游祁连县多年平均防风固沙功能为8×103tkm2，甘肃省为2.7×104tkm2，下游内蒙古自治区额济纳旗为2.3×104tkm2. 防风固沙总量在1.1×109t～4.9×109t之间变动，多年年均值为3.2×109t，其中，祁连县为1.96×107t，甘肃省为1.57×109t，额济纳旗为1.63×109t(见图7). 防风固沙功能趋势增加，年防风固沙量平均增加6.67×107t，年均变化率为1.85%.

图7 代表性年份的黑河中下游防风固沙功能变化

图6 2000—2017年黑河中下游防风固沙功能变化

2000—2015年，额济纳旗的土壤风蚀模数从8.89×104tkm2减至4.34×104tkm2(见图6)，防风固沙功能从3.2×104tkm2减至2×104tkm2(见图7)，导致整个下游区域的结果变化较大. 这表明这段时间内，流域上游天然植被封育、中游湿地保护及下游生态移民等生态保护工程取得巨大成效.

与韩永伟等[5]运用董治宝建立的风蚀流失模型[39]得到的2006年黑河下游防风固沙功能结果相比，笔者计算的黑河中下游年均防风固沙功能(2.44×104tkm2)稍高于文献[5]中黑河下游低覆盖草地、中覆盖草地、灌木林、有林地的防风固沙功能(分别为1.23×104、1.91×104、2.26×104、2.27×104tkm2)，但与荒漠区相当. 这是由于近20年来研究区气候因子及植被覆盖发生显著变化，导致流域防风固沙功能趋势增加.

2.3 防风固沙功能评估主要影响因子分析

由图8可见，各因子空间分布与防风固沙功能具有差异化的相似性. 气候因子值超过282 kgm的较强区集中在荒漠区，占研究区面积的48.40%；在212～282 kgm之间的一般区主要分布于山区、绿洲区及下游黑河带，占研究区面积的43.53%；小于212 kgm的较低区主要集中在张掖市绿洲平原区，占研究区面积的8.07%. 在空间上，额济纳旗东北部气候因子趋势强度减小，张掖市、酒泉市小部分区域趋势增加，但由于研究区气象站点较少(见图1)，空间分析误差可能较大，在此不做气候因子趋势变化的空间面积统计. 多年平均植被覆盖度超过65%的区域仅占研究区面积的5.71%，分布在中游山地、绿洲区及黑河水系边缘；73.51%的研究区面积多年平均植被覆盖度较低，在0～11%之间，主要为荒漠地带. 54.14%的研究区面积土壤可蚀性指数与土壤结皮指数均超过0.50，主要分布在荒漠区，表明易受到风沙侵蚀. 土壤可蚀性指数与结皮指数皆低于0.30的区域占研究区面积的2.40%，集中于山区，土壤抗风蚀能力较好. 研究区北部及中部较平坦，地表糙度因子多在0.96～1之间；占研究区面积的54.96%，祁连山脉地表起伏较大，地表糙度因子在0.09～0.34之间，占研究区面积的5%.

图8 黑河中下游防风固沙功能主要影响因子空间分布

2000—2017年植被覆盖度以微度增加为主，26.55%的区域趋势增强，主要位于研究区南部及黑河沿线地区. 趋势减弱区零星分布在下游额济纳旗，仅占研究区面积的1%(见图9).

图9 黑河中下游2000—2017年植被覆盖度变化趋势

张掖市防风固沙功能变化趋势增加而额济纳旗变化趋势减小(见图3)，主要原因是潜在土壤风蚀量与实际土壤风蚀量随气候因子与植被覆盖因子共同作用而变化. 张掖市植被覆盖度略增加，实际土壤风蚀量减小，而气候因子增加强度更大，导致潜在土壤风蚀量增加，因此防风固沙功能增加. 植被覆盖度增加的区域主要集中在额济纳旗黑河沿线及湖泊，而其气候因子减弱，土壤风蚀模数降低，潜在土壤风蚀量明显减小，导致防风固沙功能减小.

由于土壤特性及地表糙度变化微小，该研究视为长期不变. 由灰色关联分析结果得出，防风固沙功能的影响因子关联系数分别为0.80(风力因子)、0.65(积雪覆盖因子)、0.64(土壤湿度因子)、0.56(植被覆盖因子)(见图10)；贡献率排序表现为风力因子(30.04%)>积雪覆盖因子(24.57%)>土壤湿度因子(24.26%)>植被覆盖因子(21.13%). 可见，风力是气候因子的最主要部分，其对防风固沙功能的影响最大，其余各因子占比相当.

图10 各主要影响因子与防风固沙功能关联系数

3 结论

a) 2000—2017年黑河中下游平均土壤风蚀量约为4.6×109t，平均土壤风蚀模数为3.47×104tkm2，年均防风固沙量为3.2×109t，年均防风固沙功能达2.44×104tkm2；整体风蚀减轻，年均减小1.36×107t，年均变化率为1.60%，其中下游额济纳旗减小显著，表明生态保护工程取得成效. 土壤风蚀状况排序表现为下游额济纳旗>中游甘肃省>青海省祁连县山区.

b) 黑河中下游防风固沙功能整体呈现中游较强，向下游递减的空间分布特征，各区域差异明显. 整体防风固沙功能提升、生态环境向好. 防风固沙功能较高区、一般区和较低区分别约占研究区面积的31.54%、20.77%、47.69%. 张掖市和嘉峪关市防风固沙功能趋势增加，回归方程系数在0～26.29%之间，占研究区面积的12.51%；额济纳旗东北部和高台县中部趋势减弱，回归方程系数在-17.17%～0之间，占研究区面积的23.30%.

c) 风力因子是黑河中下游防风固沙功能变化的最主要影响因子，贡献率为30.04%，积雪覆盖因子、土壤湿度因子、植被覆盖因子的贡献率分别为24.57%、24.26%和21.13%.

d) 该风沙区的土壤风蚀防治工程应综合考虑气候变化、植被覆盖、土壤特性及人类活动的复合影响，协调生态环境保护与经济社会发展的关系，实行具有区域适宜性的方案，制定有针对性的生态保护和移民安置等工程，因地制宜调整产业结构及建设美丽村镇.

e) 研究区个别气象站点数据缺失，且研究方法未能考虑不同土壤类型临界起沙风速的差异性，导致研究仍存在限制因素. 在今后研究中，可结合实地观测数据，提高气候和土壤参数的分辨率，从而提升单位面积防风固沙量的模拟精度，并在下一步工作中结合气候预测数据进行该区域未来防风固沙功能的分析，以增进对防风固沙功能时空变化的科学认识.