WEPS模型在乌兰布和沙漠油莎豆（Cyperus esculentus）种植区的应用

廖贵云， 吴秀芹， 谭 锦， 李 丹， 冯梦馨

（北京林业大学水土保持学院，北京 100083）

风蚀是气流（风力）作用下土壤圈或岩石圈的破损，其过程是风力作用引起的地表物质脱离地表、搬运和再堆积过程的统一［1］，是沙质荒漠化的主要成因，也是干旱、半干旱及部分半湿润地区的主要环境问题之一［2］。准确预测风蚀及风蚀引起的环境变化，对于土壤保持规划、指导与检验各种土壤风蚀防治措施、减轻风蚀引起的空气污染、维护风蚀土地的可持续利用等十分必要［3］。目前，土壤风蚀估算主要有野外直接观测［4］、风洞模拟试验［5-6］、同位素分析法［7］、粒度对比法［8］和风蚀模型模拟［9-10］等方法。通常认为，一个成功的土壤风蚀模型能够预测不同尺度和不同地表类型的风蚀情况，是估算风蚀量最科学、最合理的方法［11］。

20世纪30年代末期，现代流体力学的创立使风蚀的定量化研究成为可能；随后几十年，许多国家相继建立起适用于不同尺度的定量模型来估算风蚀强度［12］。小尺度模型适用于面积分布在100～1000 m2的区域，具有代表性的有风蚀方程（WEQ）、修正风蚀方程（RWEQ）、风蚀预报系统（WEPS）、德克萨斯风蚀分析模型（TEAM）、风蚀随机仿真模型（WESS）等［13］；中等尺度模型适用面积最大可达10000 km2，如欧洲轻质土壤风蚀模型（WEELS）、风蚀评价模型（WEAM）、综合风蚀模型系统（IWEMS）；大尺度模型适用于面积在10000 km2以上区域，目前最重要的大尺度模型为粉尘释放模型（DPM）［14］。其中，WEPS 模型由美国农业部组织开发，是目前最完整和最先进的土壤风蚀预报模型，它适用范围广，不仅针对农田区域还兼顾草原地区的风蚀预测［12］。想要大范围的运用WEPS 模型，还需要拥有适合国家级的WEPS数据库［15］。美国农业部于2016 年将WEPS 软件更新到WEPS_1.5 版本，建立起相对完备的WEPS 数据库，使WEPS 模型得到广泛运用并成为农田风蚀的主要评价工具。陈莉等［16］在2012 年首次使用WEPS 模型对天津郊区土壤风蚀起尘量及迁移量进行估算，探明了郊区起尘迁移情况对市区的影响；王燕等［17］在民勤荒漠地区运用WEPS 模型对风蚀量进行预测，发现模型在模拟无植被覆盖区域时与实测值比较接近。刘珺等［18］运用RWEQ模型和WEPS模型对北方农牧交错带潜在风蚀进行评估，得到WEPS 模型模拟效果要优于RWEQ 模型。Liu 等［19］运用WEPS 模型对中国新疆南部风蚀源PM10和PM2.5排放量进行模拟，其准确性较高。

油莎豆（Cyperus esculentus）原产地为非洲干旱沙漠区，根系发达，分蘖能力强，具有抗逆性强、耐瘠薄、适应性强等优异特性；除此之外，还具有防风固沙及地力培育等生态修复功能。通过科学种植油莎豆及合理构建油莎豆防风固沙采收模式，能有效遏制北方风沙区土壤自然侵蚀和表层沙土流动，为解决我国北方沙漠化土地风蚀沙化严重、水资源短缺、盐渍化加剧等问题提供可能［20］。有效的风蚀估算是构建油莎豆科学种植模式的前提，但传统风蚀监测手段，因工作量大、监测时段长和监测范围有限等特点，加大了风蚀监测工作的难度。研究WEPS模型在内蒙古油莎豆种植区的运用效果不仅为该区域提供一种高效的风蚀监测方法，也将为我国运用风蚀模型在小尺度区域进行风蚀预测提供理论基础。

1 研究区概况

研究区位于内蒙古自治区巴彦淖尔市磴口县（图1），地处乌兰布和沙漠东缘，土地面积4167 km2，沙漠面积2847 km2。海拔在1030～2046 m, 除北部狼山山脉外，地势整体由东南向西北降低，地貌以沙漠、平原和山地为主；研究区属温带大陆性季风气候，年平均气温7.6 ℃，年平均降水量145 mm，年平均蒸发量2398 mm，水热时空分布不均，年际变化较大。土壤质地较轻，土壤类型包括风沙土、灌淤土、灰漠土、草甸土和棕钙土等，其中以流动风沙土为主。区内植被稀疏，乔木树种有河柳（Salix chaenomeloides）、胡杨（Populus euphratica）、榆树（Ulmus pumila）和梭梭（Haloxylon ammodendron）、等；草本与灌木植物有白刺（Nitraria tangutorum）、沙打旺（Astragalus adsurgens）、多枝怪柳（Tamarix ramosissima）、柠条锦鸡儿（Caragana korshinskii）、细枝岩黄耆（Corethrodendron scoparium）和沙棘（Hippophae rhamnoides）等。

图1 研究区位置示意图Fig.1 Sketch map of study area location

2 材料与方法

2.1 试验设计

在乌兰布和沙漠油莎豆种植示范区内选择3个油莎豆纯作地块及1个油莎豆和梭梭间作地块开展研究，4 个地块相邻且独立，大小均为80 m×140 m，地块内油莎豆种植方式为带状种植，每垄2 行油莎豆，纯作种植方向为东—西走向，间作呈南—北走向。将1个纯作地块内的油莎豆进行全部采收以作对照（全采收，CK），其余地块均在中部建立了不同采收方式的试验小区，小区内的采收类型分别为：留4垄采6垄（留4采6）、留6垄采6垄（留6采6）和1垄油莎豆搭配1 行梭梭间作未采收（间作留茬）3 种采收模式，大小分别为34 m×40 m、39 m×40 m、20 m×40 m。并于2020年11月15日和12月26日，对4种采收模式（图2）地表分别进行了2次风蚀观测，各模式内呈阶梯式分别均匀布设3 个集沙仪（集沙口大小为1 cm×3 cm 矩形口、共16 个梯度）和HO-BO多剖面自计式风速仪（监测距地面高度0.5 m、1 m、1.4 m 处风速）；每个模式内选取4 个1 m×1 m 样方，测定植被直立株高、叶面积指数、带宽、带间距、植株密度，并在样方内采用五点交叉法进行土壤取样，分别采用环刀法和烘干法测定0～15 cm 土层土壤容重及土壤含水率。在下午13:30—14:00在研究区上空用无人机拍摄光学正射影像，拍摄高度10 m,焦距为9 mm，最大光圈2.97，图像空间分辨率为7 mm；利用Agisoft PhotoScan 软件进行影像拼接处理，在ENVI5.3 中采用颜色混合分析法（CMA）［21］获取植被盖度。WEPS模型中结皮面积比表示为除植被外的地表可见覆盖物使地表不起沙面积比，研究中将覆盖在各小区地表的可见滴灌带面积测算为结皮面积，通过在ENVI5.3 中对无人机影像采用支持向量机法进行监督分类计算得到。2次风蚀事件中，不同采收模式风速情况及地表特征参数分别见表1和表2。

表1 不同采收模式2次风速监测Tab.1 Twice wind speed monitoring under different harvesting modes

表2 不同采收模式地表特征Tab.2 Soil surface characteristics of different harvesting modes

图2 采收模式布局及试验布设示意图Fig.2 Schematic diagram of harvesting mode and instrument layout

2.2 WEPS模型

风蚀预报系统（Wind Erosion Prediction System，WEPS）为模块化结构，由1 个用户界面、1 个管理程序、7 个子模型和4 个数据库组成，是基于过程的每日时间步长模型，不仅可以模拟基本的风蚀过程，还可以模拟改变土壤对风蚀敏感性的田间管理和风化过程。目前，WEPS_1.5 版本无中国土壤、作物和耕作管理方式等基础数据文件，不能在软件中直接选择参数模拟风蚀过程。本文采用WEPS模型中的风蚀子模型对日风蚀量进行计算，首先确定阈值摩阻风速（模型中最小值设为0.35 m·s-1），当风速超过阈值时，计算代表该地表的一系列单个网格单元的土壤风蚀量［22］。

风蚀量预测公式：

式中：Qm为单宽输沙量（kg·m-1）；U*为摩阻风速（m·s-1）;U*ts为临界摩阻风速（m·s-1）。

摩阻风速U*的计算：

（1）计算气象站摩阻风速U*f

式中：U为高度为Z时的实测风速，为10 m高度风速（m·s-1）；Z0f为气象站空气动力学粗糙度，在WEPS中取25 mm。

（2）计算观测地点摩阻风速

①计算无植被覆盖时，摩阻风速

式中：Z0为当地空气动力学粗糙度，在本文分别对应4种采收模式近地表的空气动力学粗糙度。

②计算有植被覆盖时，摩阻风速

当地表有直立植被时，需考虑直立植被对摩阻风速的削弱作用，通过判断有效植被拖曳系数的大小选择合理的公式进行计算。

（a）有效植被拖曳系数

式中：BRcd为有效生物量拖曳力系数；BRlai为直立植被叶面积指数（m2·m-2）；BRsai为茎面积指数（m2·m-2），表示茎面积与水平地表面积的比值，油莎豆是草本植物［23］，在此忽略不计。

（b）植被冠层空气动力学粗糙度

根据有效植被拖曳系数是否大于0.1，判断直立植被冠层空气动力学粗糙度计算公式如下。

式中：Z0v为直立植被冠层空气动力学粗糙度（m）；BZ为直立植被高度(m)。

（c）植被冠层上部摩阻风速

（d）植被冠层下部摩阻风速

临界摩阻风速U*ts的计算：

（1）光滑平坦地表摩阻风速

式中：SFcv为地表覆盖土块/结皮或石块的面积百分比。

（2）地表有倒放植物引起的临界起动摩阻风速增加量

式中：SCcv为倒放植被引起的临界起动摩阻风速的增量系数；BFcv为倒放植被覆盖率；UC*ts为倒放植被引起的临界起动摩阻风速增量（m·s-1）。

（3）含水率引起的临界起动摩阻风速增加量

式中：H0wc为地表含水率（kg·kg-1）；HR15wc为1.5 MPa地表含水率（kg·kg-1）。

（4）总临界起动摩阻风速

WEPS 中临界起动摩阻风速的最小值设为0.35 m·s-1。因此，可根据计算所得摩阻风速与总临界起动摩阻风速的关系，采用式（1）计算单宽输沙量。

3 结果与分析

3.1 实测风蚀结果

3.1.1 输沙量随高度分布规律 从2次风蚀事件的风沙流结构特征（图3）来看，随着油莎豆留茬数量的增加，风沙流结构特征由曲折状逐渐变为“1”字形。第1次风蚀事件中，在0～0.21 m输沙高度内，全采收模式地表不同高度的输沙量均高于其他3种留茬模式；但在0.21～0.39 m内，4种模式地表不同高度的输沙量差异较小，且随高度增加，各层输沙量趋近于0（图3a）。与全采收模式相比，3种留茬模式均能有效降低0.21 m以下的地表输沙量，表现出较好的固沙效果。第2次风蚀事件的风沙流结构特征与第1 次相似，但3 种留茬模式的输沙高度由0.21 m降低至0.15 m 以下（图3b）。结合2 次风蚀事件来看，随着留茬数量的增加，不同高度处的输沙量呈现：全采收＞留4采6＞留6采6＞间作留茬的趋势（图3）。结果表明，间作留茬模式的固沙能力最强，增加留茬数量可以提高固沙能力。

图3 不同采收模式输沙量随高度分布特征Fig.3 Sediment discharge distribution rules with height under different harvesting modes

从不同高度与输沙量的最优拟合关系来看（表3），在2 次风蚀事件中，全采收、留4 采6 和留6 采6模式地表的输沙量随高度分布均呈指数关系，且全采收模式的拟合度最高，其次为留4采6模式，最后为留6 采6 模式；间作留茬模式呈对数规律，第2 次风蚀事件的拟合度比第1 次低。表明在留茬区域，输沙量随高度呈指数分布，但随着留茬数量的增加，拟合度逐渐降低，并向对数函数方向变化。

表3 不同采收模式输沙量与高度的拟合关系Tab.3 Fitting relationship between sediment transport and height in different harvesting modes

3.1.2 不同采收模式单宽输沙量 单宽输沙量为某一时段内通过单位宽度的总输沙量，可以表征风蚀强度，该值越大，地表所受风蚀作用越强。通过计算1 d 内4 种不同采收模式0～0.39 m 高度范围内的输沙总量得到，计算公式为：

式中：QS为单宽输沙量（g·cm-1·d-1）；Qm为第m个梯度输沙量（g·cm-1·d-1）。

经式（12）计算得到4种不同采收模式下地表单宽输沙量（表4）。在2次风蚀事件中，不同采收模式间的地表单宽输沙量有显著差异，全采收模式单宽输沙量最大，间作留茬模式单宽输沙量最小；在第1次风蚀事件中，与全采收模式相比，留4采6、留6采6 和间作留茬模式分别减少了74.61%、84.50%和98.76%的单宽输沙量；第2次风蚀事件中，留4采6、留6 采6 和间作留茬分别减少了69.64%、84.86%和96.57%，与第1次风蚀事件相比，间作留茬和留4采6的固沙能力轻微减弱。

表4 不同采收模式下地表单宽输沙量Tab.4 Unit-width sediment discharge of different harvesting modes /(g·cm-1·d-1)

3.2 WEPS模型预测结果

3.2.1 模型输入预测参数 2 次风蚀监测均在油莎豆枯萎期，无灌溉，天气晴朗，气温较高，研究在使用WEPS 模型时，忽略土壤水分的影响。根据3 个风速梯度数据，用对数廓线法计算不同采收模式地表的空气动力学粗糙度。10 m 高度处的风速通过全采收模式地表观测到的风速值拟合风速廓线得到，因4种模式处于同一研究区，10 m高度处的风速几乎不受地表植被的影响，故取值一致。模型中输入不同采收模式地表的风蚀预测参数见表5。

表5 模型中输入参数Tab.5 Parameters about soil surface used in WEPS

3.2.2 模型风蚀预测 通过将各采收模式地表所测的空气动力学粗糙度、叶面积指数、植被高度等指标带入式（1）～式（11），计算得到2 次风蚀事件中，4种采收模式的摩阻风速、有效植被拖拽系数和植被冠层空气动力学粗糙度等预测参数（表6）。

表6 模型参数输出结果Tab.6 Parameters calculated with WEPS

模型中的摩阻风速是指上层风速经过植被及其障碍物削弱后作用于地表的风速，能直观反映植被及其障碍物在垂直方向上对风速的削弱能力，其值越小表示植被及障碍物对风速的削弱能力越强。模型计算各采收模式下摩阻风速大小表现为：间作留茬＞留4采6＞留6采6＞全采收，其中间作留茬模式摩阻风速最大，主要是因为油莎豆-梭梭间作留茬模式的植被结构较均匀，疏透性好，穿过冠层内部的气流较多所导致。在3 种植被覆盖模式中，间作留茬有效植被拖拽系数和植被冠层空气动力学粗糙度最大，留6采6最小。

WEPS模型预测的临界摩阻风速是地表沙粒开始移动的风速阈值。在2 次风蚀预测中，临界摩阻风速大小均为：间作留茬＞留6 采6＞留4 采6＞全采收。综合分析模型各项风蚀预测参数得到，3 种留茬模式均有一定的防风作用，其中，间作留茬模式防风效果最优，留4采6防风效果最差。

3.3 实测单宽输沙量与模型预测值

从4 种采收模式地表单宽输沙量的预测情况（图4）可以看出，不同采收模式的单宽输沙量实测值与模型值差异较大；其中间作留茬模式地表差异最大，模型预测结果最大为实测值的10.16倍，最小为实测值的0.58 倍。就模型预测效果而言，第1 次风蚀事件中，模型对全采收模式地表预测最优；在第2 次风蚀事件中，模型低估了全采收模式地表的单宽输沙量，对留4 采6 模式地表预测最优；但2 次风蚀事件中，模型对间作留茬模式地表预测均为最差，说明WEPS 模型不适用于植被覆盖度较高的地表风蚀预测。虽然，预测效果在2 次风蚀事件中有轻微差异，但模型预测地表单宽输沙量结果与实测基本同步变化，均呈现全采收＞留4采6＞留6采6＞间作留茬的趋势。分析实测值与模型值的拟合关系（图5），采用决定系数（R2）和均方根误差（RMSE）对模型精度进行综合评价。从图5 中可以看出，模型值随实测值的增加呈幂函数增加，2 次风蚀事件R2在0.91 以上，模型对第1 次风蚀事件预测的均方根误差高于第2次风蚀事件。表明WEPS模型预测结果与实测结果具有较高的相关性，模型对第2 次风蚀事件的预测精度高于第1次。

图4 单宽输沙量模型值与实测值对比Fig.4 Comparison of model and measured results of unit-width sediment discharge

图5 单宽输沙量模型值与实测值函数关系Fig.5 Functional relationship between model and measured value of unit-width sediment discharge

4 讨论

根据乌兰布和沙漠油莎豆种植区4种不同采收模式地表实际发生的风蚀过程，对WEPS 风蚀模型预测效果进行验证。单宽输沙量的模型值与实测值存在一致性，这种一致性是非线性的，与学者们在德国和美国运用风蚀预测系统对农田风蚀量模拟研究结论基本一致［24-25］。模型不仅能够有效预测土壤质地较为复杂的农田地表风蚀量变化趋势，而且对以沙粒为主的油莎豆种植区地表也同样具有较好的模拟效果。同时，模型预测的单宽风蚀量随地表植被盖度的增加而降低，与学者们利用WEPS模型模拟不同情景地表下风蚀量的研究相吻合［26］，模型能够较为真实地反映油莎豆种植区风力对不同植被盖度下地表物质的侵蚀机制，也能够较准确地预测乌兰布和沙漠油莎豆种植区不同采收方式地表的固沙能力。

从模型对地表风蚀量预测精度来看，模型预测质量具有不确定性。首先，在同1次风蚀事件中，模型预测的单宽输沙量与实测结果具有明显差异；其次，在2次风蚀事件中，模型预测不同采收模式地表风蚀量的精度差异较大。一方面可能是对风蚀事件的观测次数不足，以及在同1 次风蚀事件中对不同采收模式地表的风蚀观测重复较少，从而引起实测数据不稳定导致；另一方面，由于阵风的风向、风速和风频的不稳定性影响实际获取的风蚀量与真实情况不符［27］，所产生的监测误差导致；且WEPS建模实施每小时步长，并未考虑阵风的侵蚀效应，可能造成模型预测结果偏低［28］；此外，相对于模型建立在植被空间结构均匀的背景而言，研究区地表呈密集条状的植被结构能有效防止“狭管效应”的发生［29］，从而减少气流穿过植被冠层孔隙时，流速增大，加速对地表侵蚀的现象出现，这也是导致模型预测单宽输沙量偏高的重要原因之一。

综合来看，WEPS 模型预测的单宽风蚀量对风速、地表和植被特征响应较为敏感。模型在小尺度预测不同区域的风蚀量变化趋势上是可行的，在指导农田风蚀管理中具有重要作用。而在风蚀定量化研究上，WEPS模型预测精度不佳，可通过缩小风蚀预测时段、增加观测频次和选择持续大风天气等方法来提高模型预测精度，但这些方法都是基于改变外部环境来更好地响应模型，对模型的运用和发展具有一定局限性；未来要真正提高WEPS 的普适性，还需要加强对模型在多气候、地表植被覆盖多类型，预测多时序以及其他与风蚀相关因素中的研究，并建立数据库，从而根据实际风蚀环境对模型参数及公式进行修正。

5 结论

（1）与全采收相比，3种留茬模式均能有效降低风力对种植区地表的侵蚀作用，其中，间作留茬防风固沙能力最强，留6采6最弱，地表防风固沙能力随留茬数量的增加而提高。

（2）全采收、留4 采6 和留6 采6 输沙量随高度增加呈指数降低的变化趋势，间作留茬输沙量随高度增加呈对数降低的分布规律；随留茬数量的增加指数拟合度逐渐降低，输沙量随高度的分布规律向对数方向演变。

（3）WEPS预测与实测结果具有一定差异，模型值最大为实测值的10.16 倍，最小为实测值的0.58倍；模型预测质量具有不确定性，在植被覆盖度较高的地表预测效果较差，但存在合理的一致性，模型预测的单宽输沙量随实测值的增加呈幂函数增加。模型能够较为准确地预测不同地表特征的风蚀量变化趋势，对于风蚀定量化估算而言，WEPS模型还需建立数据库，根据实际风蚀环境对公式及参数进行修正。