王 蕾, 赵 霞, 张琛悦, 唐文家
(1.青海师范大学 青藏高原地表过程与生态保育教育部重点实验室, 西宁 810008;2.高原科学与可持续发展研究院高原土壤信息科学研究团队, 西宁 810008; 3.青海省生态环境监测中心, 西宁 810007)
我国是荒漠化严重国家，2014年全国荒漠化土地面积为261.16×104km2，其中风蚀荒漠化的面积最大且分布最广[1]，主要分为“三北”戈壁沙漠及沙地风蚀类型区[2]。防风固沙是我国北方干旱与半干旱地区的重要生态系统功能，精确地评估防风固沙功能可为生态保护红线划定、生态环境的保护与修复提供科学指导。2014年国家为保障和维护区域生态安全的底线与生命线，提出了生态保护红线战略，其中生态功能评价是划定生态保护红线的前提。2017年生态环境部颁布的《生态保护红线划定指南》规范了生态保护红线中生态功能的评价方法[3-4]。《生态保护红线划定指南》推荐采用修正风蚀模型(RWEQ)[5-6]来评估全国尺度的防风固沙功能，但未明确规定模型中气候因子WF等关键参数的获取方法。RWEQ模型在我国西北干旱半干旱地区的研究应用多集中在农田、沙地、草地混杂的农牧交错区[7-19]，而应用于以荒漠戈壁为主的县级地区的研究较少。青海省茫崖市位于柴达木盆地，属于典型的荒漠戈壁区，是我国西北典型的强烈风蚀区和生态脆弱区[20]，也是气候变化敏感区[21-22]。本研究以青海省茫崖市为研究对象，首先通过Fq代替Wf、风场强度因子计算Wf法、风场强度因子改进计算Wf法、Fq直接代替WF法这4种方法确定适宜于荒漠戈壁区的RWEQ模型参数——气候因子WF，然后基于RWEQ模型计算出茫崖市的实际土壤风蚀量并分析其防风固沙功能的空间变异特征，最后利用指南法、自然间断法和文献—实际结合法这3种不同分级方法识别茫崖市防风固沙功能的敏感区。本研究旨在通过确定RWEQ模型中的气候因子WF参数来评估茫崖市的防风固沙功能，并利用不同分级方法识别出防风固沙功能敏感的土地利用类型(简称地类)，为生态保护红线的划定提供科学指导。
1.1 研究区概况
茫崖市地处青海省西北柴达木盆地的西北角，介于37°02′—38°59′N和90°07′—93°09′E(图1)，面积约5×104km2，海拔在2 300 m以上，属大陆荒漠区气候。年日照时数3 152 h，年均气温4.14℃，年均降水量和蒸发量分别为50 mm和2 692 mm[23-24]，年均风速5.1 m/s，最大风速可达39 m/s，8级以上大风年均日数超过100 d。研究区内土壤类型主要有风沙土、灰棕漠土、漠境盐土、冷钙土和粗骨土等，植被覆盖度不足0.12 %(表1)，浮尘、沙尘暴频发，是西北典型的极干旱风蚀区[25]。2014年，茫崖市被纳入国家级“沙化土地封禁保护区”试点范围，重点开展保护区封禁设施和固沙压沙等生态工程建设。
图1 研究区概况
表1 茫崖市不同区域土壤、植被及山脉概况
1.2 RWEQ模型简介
本研究采用《生态保护红线划定指南》中推荐的RWEQ模型评估茫崖市的防风固沙生态功能。该模型计算的防风固沙物质量是无植被/作物(包括残茬)覆盖的潜在风力侵蚀量与实际风力侵蚀量的差值[26]，主要公式和因子见公式(1—7)。
SR=SL潜-SL
(1)
(2)
(3)
Qmax潜=109.8(WF×EF×SCF×K′)
(4)
S潜=150.71(WF×EF×SCF×K′)-0.3711
(5)
Qmax=109.8(WF×EF×SCF×K′×C)
(6)
S=150.71(WF×EF×SCF×K′×C)-0.3711
(7)
式中:SR为单位面积年固沙量(kg/m2)；
SL潜为潜在风力侵蚀量(kg/m2)；
SL为实际风力侵蚀量(kg/m2)；
Qmax潜为潜在风沙最大转移量(kg/m)；
Qmax为风沙滞留量(kg/m)；
S潜为潜在关键地块长度(m)；
S为区域关键地块长度(m)；
z为最大风蚀出现距离(m),通常采用下风向距离50 m[4]；
WF为气候因子(kg/m);EF为土壤可蚀因子(无量纲);SCF为土壤结皮因子(无量纲);K′为地表糙度因子；
C为植被覆盖因子。
1.2.1 气候因子WF WF表征气温、降水、蒸发和风速等气象因素对防风固沙功能的综合影响[26]，计算见公式(8)。
(8)
式中:WF为累加12个月风力因子Wf得到的年平均气候因子(kg/m)；
ρ为空气密度(kg/m3)；
g为重力加速度(9.8 m/s2)；
SW为各月的年平均土壤湿度因子；
SD为雪盖因子。
由于《生态保护红线划定指南》没有给出WF的具体计算方式，本研究采用下列4种方法式来确定WF。
(1) 采用联合国粮农组织提出的风蚀气候因子Fq代替Wf，即Fq代替Wf法，见公式(9)。
(9)
式中:u为2 m高处的月平均风速(m/s)；
u1,u2分别表示在z1,z2高度处的风速(m/s)；
ETPi为月潜在蒸发量(mm)；
Pi为月降水量(mm)；
d为当月天数；
Ti为月平均气温(℃)；
ri为月平均相对湿度(%)。
(2) 依据风场强度因子[18,27]计算Wf，即风场强度因子计算Wf法，
Wf=u4×(u4-u3)2×Nd
(10)
式中:u3为起沙风速,本研究取5 m/s；
u4为气象站点2 m处的月均风速(m/s)；
Nd为日均风速大于5 m/s的天数[4,28]。
(3) 为体现区域风速的日变化以及日风蚀量，将风场强度因子[29]公式(10)改进为公式(11)，即风场强度因子改进计算Wf法，
(11)
式中:u5为气象站点2 m处的日均风速(m/s)；
n为当月天数。
(4) 为反映气候因子对风蚀的影响，采用联合国粮农组织给出的风蚀气候因子Fq直接代替气候因子WF[30-31]，即Fq直接代替WF法。
1.2.2 土壤可蚀因子EF和土壤结皮因子SCF EF表示土壤受风力侵蚀影响的程度，由土壤质地、有机质和碳酸钙含量等因素决定；
SCF表示土壤抵抗风蚀能力的强弱[32-33]。
(12)
(13)
式中:sa为砂粒含量(%)；
si为粉粒含量(%)；
cl为黏粒含量(%)；
OM为土壤有机质含量(%)；
CaCO3为碳酸钙含量(%)。
1.2.3 地表糙度因子K′K′指因地形影响产生的地表粗糙程度，反映了土壤风蚀的程度[34]。
(14)
式中:Kr为土垄糙度(cm),采用Smith-Carson方程计算;Crr为随机糙度因子(cm),本次计算取0值;ΔH为距离L范围内的海拔高程差,L为地势起伏参数。
1.2.4 植被覆盖因子CC表示植被覆盖。植被覆盖增加了地表糙度，使起沙风速增大，这在一定程度上阻碍了土壤颗粒的移动[35]，是控制土壤侵蚀的积极因素，其计算公式如下：
C=e-(ai×SC)
(15)
式中:SC为植被覆盖度;ai为不同植被类型的系数,其中林地0.153 5,草地0.115 1,灌丛0.092 1,裸地0.076 8,沙地0.065 8,农田0.043 8。
综上，RWEQ模型中EF，SCF，K′和C均受下垫面条件的限制，可由相应数据计算获得。假设地面条件一致时茫崖市防风固沙功能主要受气候因子WF的影响，而WF则需要根据实际情况确定，因此本研究采用上述Fq代替Wf、风场强度因子计算Wf法、风场强度因子改进计算Wf法以及Fq直接代替WF法这4种方式来确定WF。
1.3 数据来源与处理
气象数据来源于中国气象科学数据共享服务网(http:∥data.cma.cn/)，本研究选取2014年和2018年各气象观测站点的月均降水、气温、风速、湿度等数据，借助ArcGIS软件通过空间插值得到栅格图层；
平均土壤湿度因子采用SMAP2017年逐日数据(https:∥nsidc.org/data/smap)，空间分辨率为1 km，通过ENVI软件合成月值；
雪盖数据来源于寒区旱区科学数据中心(http:∥westdc.westgis.ac.cn)中国地区的MODIS雪盖产品数据集；
土壤质地数据(黏粒、粉粒、砂粒)来自中国1∶100万土壤数据库中的各类粒径土壤颗粒含量空间分布图，土壤有机质的图层来自FAO中国土壤数据库的有机质含量属性，本研究按土壤类型的平均有机质含量计算；
数据高程模型(DEM)来源于中国科学院资源环境科学数据中心(http:∥www.resdc.cn)，空间分辨率为30 m；
植被覆盖来源于“全国生态环境十年变化评估”项目通过遥感解译获得地表覆被数据。
所有参数因子统一采用2000国家大地坐标系统Gauss-Kruger投影，各因子重采样按空间分辨率为250 m的栅格单元进行模型计算。
2.1 模型气候因子确定及关键因素分析
图2为通过Fq代替Wf、风场强度因子计算Wf法、风场强度因子改进计算Wf法、Fq直接代替WF法这4种方法计算获得的Wf和WF空间分布情况。从图2中可知，方法1即Fq代替Wf(图2A，B)和方法3即风场强度因子改进计算Wf法(图2E，F)获取的WF空间分布格局相似，如茫崖市西南部祁漫塔格山和乌兰乌珠尔山区均属相对高值区，但这两种算法得到的WF值域相差较大，其原因是Fq经归一化后替代了Wf。方法2即风场强度因子计算Wf法(图2C，D)利用风场强度因子计算的WF由于受月均风速u4和风力大于5 m/s天数Nd的影响，抹平或削弱了该区域风速的日变化，导致日风蚀量很大程度上被低估。方法4即Fq直接代替WF法(图2G，H)中Fq直接代替WF获得的茫崖市2014年WF空间分布由东北向西南方向减小，值域范围为19.94～86.38 kg/m，最大值出现在冷湖地区，这与祁栋林获取的柴达木地区WF的空间分布与值域范围一致。因此，选择方法4进行RWEQ模型运算更为合适[36]。
图3表明，茫崖市EF值介于0.001～0.617。已有研究表明，地表25 mm范围内粒径小于0.84 mm的土壤颗粒为可蚀性颗粒[37]，而风沙土结构疏松，发育于砂性母质，砂粒(0.05～2 mm)占90%以上，其EF最高(平均为0.524，表2)，因此极易受风蚀影响。除寒冻土和寒漠土外，其他类型土壤的EF最大值均大于0.6，棕漠土EF平均值最低(为0.229)。SCF值介于0.01～1.00，棕漠土SCF平均值最小(为0.364)，所占面积较少，仅为167.1 km2，但抗风蚀能力较强，与EF的结果一致。风沙土的SCF最大(为0.788)，其抗风蚀能力最弱(表1)。K′介于0.637～1，西南部祁漫塔格山、乌兰乌珠尔和北部边缘阿尔金山脉、玉勒肯·塔赫达依以及东部赛什腾山等山脉的K′值较低，这是由于山脉纵横海拔较高、地形起伏明显，不利于发生风蚀；
而低地势、低海拔的盆地K′较高，易发生风蚀。C值介于0.005～1，其在西南角的山区较低，中、东部区较高，近90%地区C值接近1，由于茫崖市地处柴达木盆地，自然景观多为干旱荒漠，植被稀疏，对风蚀的抑制作用弱，因此更易发生风蚀。
2.2 茫崖市实际风蚀量和防风固沙量
图4为采用方法4(Fq直接代替WF法)的RWEQ模型计算出的茫崖市2014年风力侵蚀量和防风固沙量的空间分布情况。根据《土壤侵蚀分类分级标准(SL190-2007)》，RWEQ模型计算的茫崖市2014年实际风蚀量可分为6级(图4和表3)。结果显示，茫崖市2014年平均风蚀模数(单位面积实际风力侵蚀量)在57.42 kg/(m2·a)以内，平均风蚀模数总量为57.99 kg/(m2·a)，实际风蚀总量为116 814.62万t，属极剧烈等级(实际风蚀量为106 421.17万t)，占实际风蚀总量的近91%，极剧烈侵蚀区面积大于3万km2，占全市总面积的60%以上。剧烈等级的实际风蚀量次之，全市约10 %的区域平均风蚀模数为11.88 kg/(m2·a)；
除极剧烈和剧烈侵蚀外，其余等级的平均风蚀模数均较小，实际风蚀量也较少。空间分布上，2014年全市极剧烈侵蚀区主要集中在中东部，主要原因是该区降雨量小蒸发量大，植被盖度低，风场强度大，而且土壤类型主要为荒漠风沙土；
风蚀强度较轻的区域主要分布在阿尔金山脉、祁漫塔格山以及乌兰乌珠尔等西部山区。图4与图3相比，单位面积实际风蚀量的空间分布格局与EF和SCF的格局相似，表明EF和SCF对风蚀影响较大。
图2 4种算法计算的风力因子和气候因子原值
图3 防风固沙功能关键影响因素
表2 茫崖市土壤可蚀因子和土壤结皮因子
图4 茫崖市2014年实际风力侵蚀量和防风固沙量
2014年全市单位面积防风固沙量介于0～44.93 kg/m2。茫崖市90 %为未利用地，属自然景观和土地类型单一的荒漠戈壁区，仅西南角的少部分山区抵抗风力侵蚀和防风固沙的能力较强，整体上茫崖市防风固沙能力弱，属风蚀严重区，研究结果与茫崖市的区域地带性(西北干旱荒漠区)的实际情况相符。因此，进一步证实Fq直接代替WF法确定的模型参数，应用于干旱半干旱荒漠戈壁区的防风固沙功能研究中具有可行性。
表3 实际风蚀量等级及变化
2.3 防风固沙功能评估分级与敏感地类识别
为更细致地反映茫崖市防风固沙功能的空间分布特点，以及更好地识别出对防风固沙功能变化敏感的地类，本研究分别采用了《生态保护红线划定指南》推荐的三级分类方法(简称“指南法”)、自然间断法和文献—实际结合法对茫崖市防风固沙功能进行分级，分级依据见表4。指南法评估技术将防风固沙功能服务值(单位面积年固沙量SR)从高到低排序，以累加服务值占服务总值比例的50 %和80 %所对应的栅格值作为生态系统服务功能评估分级的分界点，划分为“1等：一般重要”，“2等：重要”和“3等：极重要”三级；
自然间断法是利用ArcGIS软件将优化后的RWEQ模型计算出的SR进行自然分组，使各等级之间的差异最大化，SR差异相对较大的栅格值设为分级裂点；
文献—实际结合法主要是依据茫崖市SR值域范围以及实际风蚀量等级(表3)，并结合已有研究设置对应的SR分级标准[38]。同时，通过增加分级数量来更详细地反映功能原值SR中蕴含的丰富空间分异特征。
表4 防风固沙服务功能评估分级方法
从图5(仅显示差异最明显的西南角)可知，指南法、自然间断法和文献—实际结合法3种评估分级方法所得的茫崖市防风固沙功能空间分异具有相似的规律，2014—2018年防风固沙功能等级的空间分布特征差异较小，其中文献—实际结合法(图5E，F)刻画的茫崖市西南山区的防风固沙功能识别范围(除一般重要等级)比其他两种方法要广，能够在较大程度上识别出高等级区域，避免因等级划分低而忽视相应的保护措施，并且其极重要等级分布区域面积介于其他两种方法之间。因此，在考虑实施生态保护措施的时间与经济成本时，文献—实际结合法评估分级结果具有更好的参考价值。
通过ArcGIS将土地利用类型和防风固沙功能空间分布进行叠置分析，来反映不同土地利用类型的防风固沙能力以及对功能变化敏感的地类[39]。结果表明，2014年防风固沙能力最强的土地利用类型为草地，固沙量达到1 753.9×104t，其次为未利用地和水域(主要是人工河渠、水库等)(表5)，而城乡居民工矿用地(主要是人工建设地面)、林地和耕地的防风固沙能力较弱，说明草地在防风固沙中起着主导作用。2014—2018年防风固沙能力变化最敏感的土地利用类型也是草地(敏感系数3.77)，防风固沙量增加了108.1×104t，原因可能是生态系统的固沙能力随着草地覆盖度的增加而增强(表6)。从防风固沙功能重要性等级(普适性的指南法)上看，2018年一般重要区域(1级)的主要土地利用类型是未利用地，占茫崖市面积的近90%，其次为草地、城乡居民工矿用地和水域等；
重要区域(2级)和极重要区域(3级)面积最大的土地利用类型均为草地，出现这一现象的原因是茫崖市固沙能力较强的草地面积占比相对较少，仅占7%左右，而未利用地面积虽大但保沙率却较低。2014—2018年各土地利用类型(除耕地外)的防风固沙功能等级属于一般重要(1级)，这与茫崖市气候干燥、水资源贫乏的荒漠环境有关。功能等级变化最敏感的土地利用类型是水域，其面积减小导致功能1级区域占比降低。
图5 3种防风固沙功能评估分级结果
表5 不同土地利用/覆被类型的防风固沙能力及敏感系数
RWEQ模型在我国西北干旱半干旱地区农田、沙地、草地及其交错区的应用研究报道较多[7-17]，能否用于以荒漠和戈壁为主的县级地区尚不明确，本研究表明通过选择适宜的WF获取方法，RWEQ模型也可用于以荒漠和戈壁为主的干旱半干旱地区。RWEQ模型的关键因子是气候因子WF，方法1将归一化后风蚀性气候因子Fq代替风力因子Wf，这使计算结果不能反映原始值域的变化。方法2与方法3的起沙风速统一取值为5 m/s未考虑到区域土壤质地对风力搬运土壤颗粒的影响，因此本研究最终采用方法4即Fq直接代替WF法，其空间分布与值域结果与祁栋林[40]、董玉祥等[41]获取的中国干旱与半干旱地区风蚀气候侵蚀力极强区——柴达木盆地一致，这说明茫崖市的风蚀气候因子主要受风速限制，荒漠景观是强风蚀长期作用的结果，并且人类活动也能在一定程度上通过改变植被覆盖度来影响风蚀强度[18-19]。模型中的土壤可蚀因子EF、土壤结皮因子SCF、地表糙度因子K′短期内变化不明显，其空间分异规律与前人的研究结果相似[18]。在环境相似地区的风蚀状况对比上，本文与江凌等[18]对青海省中度以上风蚀区的研究以及迟文峰等[14]对内蒙古高原荒漠生态系统的研究结果一致。
在功能分级中，《生态保护红线划定指南》规定的累计百分比分级方法能够较好地找出防风固沙功能变化明显的大致空间分布区域，但是由于三级分级方法不够精确，使得目前的分级结果不能反映功能原值中蕴含的空间分异特征，降低了评估结果的政策参考价值。因此本研究探索的文献—实际结合法不仅在分级数量上能够较详细地刻画SR的空间分布特征，还能在功能原值的基础上结合茫崖市的实际风蚀量，具有明确的功能分级依据。
茫崖市2014年的平均风蚀模数即实际风蚀量为57.99 kg/(m2·a)，与贺倩等[42]基于RWEQ模型计算出的青海三江源地区2015年实际风蚀量的最大值60.86 kg/m2相近，且最大风蚀量集中在三江源地区的西部，临近柴达木盆地的茫崖市。本研究中茫崖市的防风固沙功能空间分布特征反映出该市整体防风固沙功能较低，这与温豪[43]的模拟结果一致，证实了本研究结果的准确性。综上，RWEQ模型应用于荒漠戈壁为主的茫崖市县域研究时具有可行性，可为生态保护措施的制定与实施提供理论依据。
表6 2014-2018年土地利用/覆被类型面积变化矩阵 km2
(1) 利用Fq直接代替WF法获得的空间分布和阈值范围与实际最相符，RWEQ模型可较准确地评估茫崖市的防风固沙功能，能为县级地区防风固沙功能评估技术的应用研究提供参考。
(2) 茫崖市2014年的平均风蚀模数为57.99 kg/(m2·a)，属风力侵蚀严重区；
2014年全市单位面积防风固沙量介于0～44.93 kg/m2，防风固沙功能整体较弱。
(3) 文献—实际结合法更能详细地反映茫崖市内部防风固沙功能的空间变异特征，且高等级功能的区域范围较广，西南部的防风固沙能力略强于中东部。
(4) 在耕地、林地、草地、水域、城乡居民工矿用地、未利用土地6大类土地利用类型中，草地是茫崖市2014—2018年期间对防风固沙功能变化最敏感的土地利用类型，敏感系数最大3.77。
猜你喜欢防风固沙风蚀气候因子延怀盆地不同土地利用类型土壤风蚀物特征水土保持研究(2022年1期)2022-12-27土壤风蚀可蚀性研究进展评述中国水土保持科学(2022年1期)2022-03-25秦王川灌区种植春小麦与披碱草对耕地风蚀的影响差异绿色科技(2018年16期)2018-09-22气候因子对烤烟质量风格特色的影响江苏农业科学(2017年5期)2017-04-15基于GIS技术的山西忻州精细化酥梨气候区划科技与创新(2016年18期)2016-11-04天山北坡NDVI对气候因子响应的敏感性分析湖北农业科学(2014年21期)2014-12-22沙漠地区微波地表发射率年内变化规律与气候因子的关系分析自然资源遥感(2014年3期)2014-02-27保护性耕作对土壤风蚀的影响中国水土保持科学(2011年3期)2011-05-12
