农田近地表风沙流风程效应变化特征研究

冯 哲， 刘瑞娟， 白宇晨， 常春平,2， 郭中领,2，李继峰,3， 王仁德， 李 庆

（1.河北师范大学地理科学学院，河北 石家庄 050024；2.河北省环境演变与生态建设实验室，河北 石家庄050024；3.河北省环境变化遥感识别技术创新中心，河北 石家庄 050024；4.河北省科学院地理科学研究所，河北 石家庄 050021）

风蚀是造成土壤退化和粉尘排放的重要原因之一［1］，准确理解风沙流的变化规律是风蚀定量研究的重要环节。风沙流从可蚀边界开始，沿着下风向运移，输沙通量随着距离的变化而变化的过程为风程效应［2］，输沙通量随着下风向运输距离的增大而增大，直到输沙率达到一个稳定值，也就是饱和输沙率（fmax）。输沙率从0 到最大平衡饱和态之间，风所运动的距离为饱和路径长度（Lsat）［3］。早在20世纪30年代，Chepil等［4］就对风程效应进行描述：风的输沙能力随着运输距离的增大而增大。人们通过野外试验、风洞试验和数值模拟等方法，对不同环境下风程效应的机理、模式和影响因素进行了大量研究，已取得了很大的进展。Chepil［5］提出“雪崩机制”解释了几十米高空的风程效应。在这个理论中，跃移的风蚀颗粒撞击地表，击溅起更多的运移颗粒。在风速稳定的情况下，输沙通量随风程距离的增加而呈指数增长，直到风沙流达到饱和状态。Owen［6］提出了一种空气动力学反馈机制，基于大量风洞试验验证，风在对地表作用过程中，将其动能转移到跃移颗粒中，增加了地表的空气动力学粗糙度。由于较高的空气动力学粗糙度，导致更大的动能传递回地表，进一步增加了风沙流的运输能力。Pähtz 等［7］引入物理学中“弛豫过程”的概念来解释这种机制。Gillette等［8］在Owen研究的基础上，通过田间试验，结合实时气象数据、输沙通量以及表层土壤性质，提出了土壤阻力机制来定量地解释更大空间尺度上的风程效应。风程效应模型应用于定量模拟研究不同机制下的风沙流。初期风程效应的探讨通常基于假定的理想状态下，其中的参数为风速、气象条件、结皮、团聚体粒径等［9-10］。许多研究通过风洞实验［6,11］、数值模拟等方法验证，这些模型在描述自然环境中风程效应的机制和影响因素时存在较大误差［12］。Bauer 等［13］基于海岸海滩构建了三角函数模型，模拟狭窄海滩风沙流的风程效应，但仅适用于海滩宽度小于Lsat的狭窄海滩环境。Stout［14］根据风沙流发育的自平衡原理提出了指数生长模型以描述某一高度上颗粒物沿沙床表面运移过程，能够较好地模拟风沙流沿程变化的“快-慢”发展模式。但其忽视了风沙流中颗粒碰撞以及颗粒撞击地表对风程效应产生的影响，因此难以描述自然环境下的沿程风沙流［15］。Andreotti［16］在Stout模型的基础上进行了修正，并通过风洞模拟和田间测量2 种方法来研究Lsat与风速的关系，结果表明，Lsat的变化和风速无关。Fryrear等［17-18］提出指数增长模型，用于描述不同下风距的风沙流，输沙率沿输沙距离的增大呈“慢-快-慢”的发展模式。

近地表风沙流运输机制受多方面因素的影响。Lynch 等［19］在北爱尔兰海滩进行研究，对风的负反馈机制进行描述。实验表明输沙通量在短距离内增长迅速，随着地块长度逐渐增加，输沙通量稳定衰减。Zhang等［20］在腾格里沙漠东南部对Lsat进行研究，发现输沙通量随着地块长度的增加，先增大、后减小、再增大。Zobeck 等［21］和Stout 等［22］结合地表因子和气象数据进行实地研究，发现可蚀性土壤干团聚体颗粒和地表空气动力学粗糙度、以及测量高度是Lsat的主要影响因素。在风洞试验中，有充足沙源供给，干松土壤特性［23］、持续吹风时间和风速大小［24］是Lsat的主导因素。在盐湖和戈壁，地表湿度、土壤物理（生物）结皮也会影响Lsat［2,8］。在海岸沙滩地带，海滩几何形状和风速对Lsat起着主导作用［13］。这些机制揭示了地表可蚀性因子状况与风因子状况决定着风程效应的发展变化过程。

近地表风沙流的精准定量观测仍是一个迫切需要解决的问题。以往的研究中，无论是经过多次改进的Bagnold 集沙仪［9］，还是国际广泛应用的MWAC集沙仪［25］和BSNE集沙仪［26］都属于机械型的集沙仪，必须在一场风蚀事件或一段时间后，人工收集风蚀物样品，带回实验室称量才能得出风沙流观测数据［27］。受仪器水平限制，风程效应的研究通常是基于较大的时间尺度，一些风程效应主要研究的相关时间尺度如表1所示：

表1 风程效应主要研究的相关时间尺度Tab.1 Main studied related temporal scale for monitoring the fetch effect of aeolian sand transport

研究发现，采用较大的时间尺度研究风程效应，一些关键变化特征难以避免地因平均化处理而被忽略［37］，野外环境中风沙流的平衡饱和态是不断变化的，风程效应研究受时间尺度的影响比较显著，选取较小的时间尺度进行观测能够更加详细地捕捉风程效应参数变化特征。基于此，本研究使用高分辨率自动连续称重式集沙仪［25］，分析近地表5 cm 高度的风沙流5 min 时间尺度的时空变化特征，并对风程效应的影响因素进行讨论。

1 材料与方法

1.1 研究区概况

河北坝上地区属于北方农牧交错带，地势相对较高，全境平均海拔1450 m，位于中温带亚干旱区，属典型温带大陆性气候。年均温1.2 ℃，年均降水量350 mm，年蒸发量1772 mm，无霜期114 d，年均风速2.99 m·s-1，是风蚀的易发区，主要发生在春季［38］。该区域土地利用方式以耕地、林地、草地为主，是传统农业区与畜牧业区耦合作用区［39］，农田风蚀灾害严重，是京津冀地区重要沙源地之一。农田地表土壤层薄，质地松散且较为均一，以沙质土和沙壤土为主，风蚀强烈［37］。本研究于2017、2018 年和2021年观测了4 次风蚀事件，观测场地位于河北坝上地区康保县境内典型旱作农田，场地1 位于六段村（41°55′37″N、114°27′28″E，海拔1480 m），场地2 位于马场（42°07′54″N、114°48′22″E，海拔1303 m）（图1a），具体观测时间及地点如表2 所示。在后文中，将2017 年4 月19 日试验称为风蚀事件1，2018 年4月30日试验称为风蚀事件2，2021年4月15日试验称为风蚀事件3，2021年5月4日试验称为风蚀事件4。

表2 风蚀事件基本信息Tab.2 Basic information of wind erosion events

图1 试验样地位置及仪器布设Fig.1 Location of the test sites and layout of the test instruments

按照当地耕作习惯，试验场地为旱地，于每年3月下旬进行翻耕耙平，实验时段内地表均无植被覆被，是典型翻耕农田，4次试验前地表的可蚀性颗粒（＜0.84 mm）分别占80.26%、78.28%、83.71%、70.19%，说明4 次风蚀事件的地表均具有高度可蚀性。试验场地大致为西北—东南走向。6段试验场地的上风向为天然草地，北面为小麦田，下风向为耕地。马场试验场地的西、北、南侧为灌丛沙堆，下风向为耕地。为了避免灌丛沙堆影响迎风气流，在试验中采用防尘网制作的防尘屏障建造了一个20 m 宽的人为不可蚀边界，布置在集沙仪上风向10 m处。

1.2 实验仪器与方法

4次风蚀事件所采用的集沙仪均为河北师范大学自主设计研发的自动连续称重式集沙仪［25］（Cy⁃clone type continuous-weigh sand trap, CCST），如图1e、g 所示，集沙仪由连续称重模块、数据采集模块和电池模块构成。称重模块最大量程为400 g，分辨率0.01 g，CR6 数据采集器（Campbell Scientific Inc.）采集频率可调节，最高达50 Hz，野外试验采用1 Hz的采集频率，能够实现对风沙流的高频连续观测。CCST集沙口朝向可以根据风向自行调整，试验中6个集沙口均沿主风向设置，高度均设置为0.05 m。风蚀事件1 和2 中的CCST 集沙仪均布设在距不可蚀边界下风向20 m、40 m、70 m、100 m、130 m 和160 m 处，风蚀事件1 的集沙口朝向设为270°，风蚀事件2 的集沙口朝向设为310°。风蚀事件3 和4 的CCST 集沙仪布设在距不可蚀边界下风向10～110 m处，每个集沙单元间距20 m，以4×6 矩阵布设，集沙口朝向设为315°。由于CCST集沙仪在使用中可能出现集沙口堵塞、被风吹歪或数采失灵情况，因此，在使用中将选取数据质量较好的一列进行分析。

野外观测风速采用2种风速仪，风蚀事件1和2的观测采用北京师范大学研制的便携式梯度风速仪［40］，如图1f 所示，9 个风杯高度设置为0.05 m、0.1 m、0.3 m、0.5 m、0.8 m、1.0 m、2.0 m、3.0 m 和4.0 m，风速分辨率为0.1 m·s-1，并采用CR1000数据采集器记录1 min 平均风速数据。梯度风速仪计算方法如下：

将实测风速数据通过最小二乘法进行拟合，得到计算风速的公式：

式中：UZ是高度为Z时的风速（m·s-1）；A、B为回归系数。

风蚀事件3 和4 的观测采用三维超声风速仪（Met one）［41］，如图1h 所示，测量高度为1.2 m，测量频率为50 Hz，数据采集频率为1 s。仪器均布设于距不可蚀边界60 m 的位置。Met one 风速仪计算风速的过程如下：

由于电子仪器在野外应用中存在不稳定现象，在数据分析前，首先计算出1 min的平均值，并将无效值剔除，然后根据水平风速计算出风向偏角：

将Met one 风速仪坐标根据偏角调整至来流风向与主风向一致，进而得到新的坐标系，计算过程如下：

式中：u1、v1、w1为x、y坐标轴调整偏转角度后的三维风速（m·s-1）；ux、vy、wz为实测坐标系下的三维风速（m·s-1）；u、v为x、z坐标轴调整偏转角度后的三维风速（m·s-1）。

以1 min 为时间间隔，计算主风向风速值（U，m·s-1），计算公式如下：

式中：uˉ、vˉ为坐标轴调整偏转角度后水平和垂直方向风速的平均值（m·s-1）。

Fryrear 和Saleh 提出的风程效应模型适用于模拟实际环境下的风沙流沿程变化，模型表达了输沙通量的沿程变化遵循慢-快-慢的变化规律，呈S 型曲线增长，并且其关键路径长度不同［17］。风程效应模拟方程如下：

式中：f（x）为距离x时的输沙率（g·cm-2·min-1）；fmax为饱和输沙率（g·cm-2·min-1），fmax≈0.63f（x）；b为关键地块长度（m）。

应用R 语言软件的最小二乘非线性回归（Gauss-Newton 迭代算法）计算5 min 时间尺度下4次风蚀事件的输沙通量，得到fmax和b，并采用决定系数（R2）和均方根误差（RMSE）检验模型的拟合效果。饱和路径长度（Lsat）是风程效应模拟中的一个重要参数，指当输沙通量接近饱和时的地块长度，即当fmax-f（x）＜fmax/100时相应的地块长度［1］。

2 结果与分析

2.1 近地表风沙流风程效应特征

分析发现，输沙通量随着集沙时间的推移，有逐渐增大的趋势（图2）。风蚀事件1、2、4在初始1～2 h 左右，输沙通量增长速度达到最大，然后逐渐放慢，能够较好地反映风沙流输沙率沿程的变化规律。输沙通量增幅大的时刻与风速增大的时刻基本对应，说明风速是影响风沙流搬运能力的重要因素。输沙通量随风程距离的增加而增加，但风蚀事件270 m处累积输沙通量小于100 m处，这可能是由于100 m 处集沙口朝向在大风条件下发生偏移，影响集沙效率。风蚀事件3试验场地11:30开始下雪，地表含水量增大，根据风速仪和集沙仪数据反映，在12:00—14:00 虽然风速较大，但输沙通量几乎没有增长；14:00后，随着风速的增大，累积输沙通量波动增长。

图2 风速和不同距离处累积输沙通量Fig.2 Wind velocity and cumulative sediment flux at different distances

分析2017年4月19日10:40—14:10，2018年4月30日12:10—15:45，2021年4月15日13:25—14:55，2021年5月4日14:45—16:15的沿程输沙率数据（图3～5）发现，5 min 时间尺度内，近地表输沙率大体上随着地块长度的增加而增加。并且，输沙率增长呈波动上升的趋势，总体上沿程变化特征为“慢-快-慢”曲线，大致呈现出S型曲线的增长规律，这与Fryrear等［17］提出的风程效应模型相吻合。

图3 风蚀事件1输沙率的沿程变化Fig.3 Variations of aeolian transport rate along the path of wind erosion event 1

2.2 饱和输沙率与风速的关系

通过风程效应模型的预测能够较精确地反映该农田地块的风沙流沿程变化特征。利用R语言软件中的非线性回归分析对4 次风蚀事件5 min 时间尺度的近地表fmax进行预测，采用Fryrear模型［17］进行拟合，并且使用拟合回归R2和RMSE 检验模型的拟合情况。经过相关性检验后，将R2＜0.40 的数据剔除。其中，事件3 和事件4 由于仪器原因，只取第2列数据进行计算。计算结果如表3所示：

表3 5 min时间尺度的模型拟合结果及统计学参数特征Tab.3 Model fitting results and statistical parameter characteristics at 5 min time scale

R2越接近1，拟合效果越好，RMSE 越接近0，误差越小。剔除R2＜0.40的数据后，Fryrear模型预测结果的平均R2为0.86，平均RMSE 为0.25。模型对于风蚀事件2 的风程效应有着较好的模拟效果，对于其他3 次的模拟效果一般。研究发现，在风蚀事件发生过程中，fmax与Lsat变化波动较大，风速对fmax的沿程变化有显著的影响。

3 讨论

3.1 近地表风沙流风程效应的时间尺度

通过分析4次风蚀事件中风沙流的沿程变化发现，近地表累积输沙通量随着风速的增大而增大，风蚀物运移主要分为2 个阶段，初始阶段为增长阶段，以流体侵蚀为主，随着地块长度的增加，输沙率迅速增大；第二阶段为弛豫阶段，流体侵蚀率随着下风距离的增大而减小，跃移沙粒撞击侵蚀率逐渐增大到极大值，输沙率达到饱和状态，地表风蚀物运输逐渐达到平衡，风蚀作用显著，风沙流在这一阶段达到fmax，Lsat就是风沙流弛豫过程的最后阶段［17］。并且，根据集沙仪数据显示，随着地块长度的增加，输沙通量有逐渐增大的趋势，这与Taylor 等［43］的研究结果相一致。

在研究中发现，5 min 时间尺度下，近地表的输沙率沿程变化呈现出明显的波动性，风程效应显著。4次风蚀事件的输沙率均表现出随着地块长度增加而增长，在达到峰值后出现降低趋势的特征。风蚀事件1的输沙率沿程变化呈现出“慢-快-慢”的增长趋势（图3），近似凸型。30 m距离内，输沙率增长较慢，30～130 m 处迅速增长达到峰值，并且风沙流达到饱和状态，130～160 m 距离内的输沙率明显降低，风沙流重新趋于稳定；10:40—12:10输沙率表现为直线增长，12:10—14:10 输沙率增长过程表现出轻微的波动。风蚀事件2 的输沙率表现出“快-慢-快-慢”的阶梯型增长趋势（图4），波动性更加显著；输沙率在70 m处达到第一次小峰值，随后降低，在100～130 m 处输沙率迅速增长，直到第2 次更高的峰值出现，130～160 m处输沙率降低。风蚀事件3在13:25—14:25 的变化过程与第一次相似（图5a～c），14:25—14:55 输沙率在测量的120 m 范围内呈现出迅速增长的趋势，并且输沙率始终保持较高水平，根据模型的计算结果表明，该段时间内fmax的变化范围在5.21～18.26 g·cm-2·min-1之间，Lsat的变化范围在96～280 m 之间，风沙流搬运能力非常强烈，风蚀作用显著。风蚀事件4 的输沙率呈现出“快-慢-快”的增长趋势（图5d～f），符合S型增长曲线，在70 m范围内，输沙率保持稳定增长，70～90 m的输沙率有所降低，90～110 m 的输沙率增长速度较快，输沙率始终保持较低水平。研究发现，5 min时间尺度下，fmax明显波动。风速是影响输沙率的重要因素，拟合回归分析4 次风蚀事件中的fmax及其影响因素风速之间的关系发现，fmax随着风速的增加而增加（图6）。根据统计学中R2的计算结果显示，fmax与风速之间呈现出显著的幂函数关系，并且时间尺度越大，其相关性越强。

图4 风蚀事件2输沙率的沿程变化Fig.4 Variations of aeolian transport rate along the path of wind erosion event 2

图5 风蚀事件3和风蚀事件4输沙率沿程变化Fig.5 Variations of aeolian transport rate along the paths of wind erosion event 3 and event 4

图6 饱和输沙率与风速的幂函数拟合关系Fig.6 Power function relationship between the saturated aeolian transport rate and wind velocity

根据模型计算结果表明，Fryrear 模型［17］对阶梯型增长趋势的风沙流有着较好地模拟效果。但在风蚀事件4 中，输沙率在80 m 内基本达到饱和状态，实测输沙率仍呈现持续增长的趋势。并且根据表1，风蚀事件4的R2为0.62，RMSE值为0.05。统计学上讲，R2越接近1，RMSE 越接近0，预测值与实测值的拟合效果越好，这表明Fryrear模型［17］可能无法很好地模拟5 min 时间尺度下S 型增长趋势的风沙流。风沙流沿程变化过程受风速、地表粗糙度和土壤性质影响。5 min时间尺度观测研究发现，风沙流的风程效应明显，Lsat波动性较强，除了受风速的影响，地表土壤干团聚体构成、土壤含水率、地表粗糙度等因素的变化对风程效应的影响也不可忽视。然而，目前对风蚀事件过程中地表风蚀影响因子的观测方法尚不成熟，尤其是土壤干团聚体构成和地表粗糙度（微地貌）的观测尚处于探索阶段［1］。因此，未来农田地类风程效应研究中，应着重关注地表风蚀影响因子的实时变化对风程效应的影响。

3.2 Lsat的影响因子

对Lsat影响因素的探究，多数是基于有定常风速、地表平坦，颗粒物松散的可控风洞环境下得出的，Lsat的变化机理在自然环境下进行探究一直是一个亟待解决的问题。自然环境中土壤表面特征存在多样性，以及风的湍流状况，导致风程效应的发展变化通常表现得非常复杂。出于对土壤表面状况实时监测难度较大的限制［11］，目前进行实时定量的Lsat的研究存在不足。本研究基于农田地类的4次风蚀事件中，5 min 时间尺度下近地表风沙流的Lsat进行探究。数据显示，Lsat有着显著的变化趋势，波动起伏较大，如表1 所示，变化范围在11～280 m之间。将其与可能的影响因素风速进行探讨，发现其发展变化不依赖于风速。

在风因子状况不稳定的野外环境下，Stout等［22］认为，地表特性和风力因子的变化影响着1 h 时间尺度下风程效应的变化规律。Lynch 等［19］以15 min时间尺度对海岸海滩进行研究，研究表明海岸海滩的Lsat主要受沙滩表面湿度、风速和海滩几何形状的影响。笔者认为，在农田条件下，近地表的Lsat受多方面因素的影响，包括风因子、跃移颗粒位移距离和地表状况（土壤干团聚体粒径组成、地表湿度、田间垄高和地表植被覆被等）等，但最终起决定性因素的是其中最慢的限制性因子［15］。在风蚀发生过程中，土壤表面性状也在实时发生变化，直接影响着大气湍流，Lsat也因此而发生改变。在本研究中，由于仪器及观测技术的限制，对地表性状的监测尚存在不足，因此，在后续的研究中，要对地表可蚀性因子进行更多的研究，在风蚀事件发生过程中，实时地测量农田土壤表层的变化。此外，在对风程效应的研究中发现，仪器精度对实验有着很大的影响。本次实验采用的CCST集沙仪属于电子类集沙仪，优点是数据采集频度高，能够实现实时连续地监测输沙通量的变化，但在大风条件下会发生抖动，产生震动误差［25］，在数据计算前需进行平均化处理，Lsat是输沙通量达到饱和时的路径长度，一些较大的输沙通量数据可能会在数据处理中被过滤，从而影响fmax和Lsat的计算结果。总结以上，在研究中要重视时间变异性，对近地表的农田尺度上对各类影响因子展开讨论，建立风沙流的三维模型，对多种的土壤类型进行更丰富的研究，以揭示风蚀农田地表Lsat的影响因素与沿程变化规律。

4 结论

本研究基于河北坝上农田地类2017—2021 年的4 次典型风蚀事件，对近地表5 cm 高度风沙流沿程变化进行分析，并应用Fryrear模型模拟近地表的fmax和Lsat。主要结论如下：

（1）近地表的风沙流输沙通量受下风向距离的影响显著。初始阶段输沙通量随着地块长度的增加迅速增长，当风程距离达到一定长度时，风沙流达到稳定状态，输沙通量逐渐达到饱和，增长速率缓慢，风蚀作用显著。

（2）4次事件风程效应有着明显差异，沿程输沙率呈现出凸型、阶梯型及S 型3 种主要变化趋势。根据统计学相关系数R2和RMSE 值表明，Fryrear 模型能够较好地拟合输沙率呈现阶梯型和凸型变化趋势的风程效应，对于S 型曲线的风程效应的拟合效果较差。

（3）风蚀事件1、2、3、4中fmax的变化范围分别为0.57～5.97 g·cm-2·min-1、0.12～1.18 g·cm-2·min-1、0.07～18.26 g·cm-2·min-1、0.08～2.02 g·cm-2·min-1，分析4次风蚀事件的观测结果，揭示了5 min时间尺度下的fmax主要影响因素，fmax与风速呈幂函数关系，并且时间尺度越大，两者间的幂函数相关关系越强。

（4）风蚀事件1、2、3、4中Lsat的变化范围分别为67～104 m、41～258 m、58.51～280 m、11～127 m。4 次风蚀事件中Lsat的发展变化具有明显波动性，并且均不依赖于风速，说明在野外环境中，Lsat的实时定量变化是多种因子共同影响的结果，其变化与地表可蚀性因子及地表微地貌的的实时变化有着紧密联系。

（5）在野外非定常风场环境中，Lsat的发生机制更为复杂，其影响因素有待进一步研究。综上所述，需要在丰富的环境中对不同土壤类型和质地进行更深入的研究，进一步探究Lsat的影响因素及其发展变化机理。