张 帅，原伟杰，刘 卉，王海霞，官昊慧，王鹿振

（1.中国林业科学研究院华北林业实验中心/北京九龙山暖温带国家森林长期科研基地，北京 102300；2.沃德兰特（北京）生态环境技术研究院有限公司，北京 100080；3.青海省水利水电勘测规划设计研究院有限公司，青海 西宁 810000）

风力侵蚀是指在气流冲击作用下土粒、沙粒脱离地表，被搬运和堆积的过程［1］。风蚀造成了严重的自然灾害和环境问题，尤其是在干旱半干旱地区，往往成为制约农业生产和导致环境恶化的主导因素，影响了这些地区的资源开发和社会经济的持续发展［2-3］。为有效降低风沙灾害，研究人员开展了多种风蚀防护措施的研究，如防护林、生物结皮、固沙剂及沙障等［4-7］。其中，沙障因其布设简便、防护效果好及防护时间长等优点，是我国防沙治沙工作中的重要工程措施［8-9］。在风蚀防护措施的防护效益评价中，地表粗糙度因其对下垫面改变的敏感度最高，得到了广泛的应用［10-11］。

由于地表起伏以及粗糙不光滑的性质，气流流经地表时受摩擦力的影响，导致近地表处风速降低，在垂直高度上呈现出不同分布规律，这种风速随高度分布的曲线称为风速廓线［12］。在风速廓线上风速为零的位置并不出现在高度为零的地表，而是高出一定距离，这一高度被称为地表粗糙度，又称空气动力学地表粗糙度［13］。近代流体力学奠基人普朗特指出粗糙度概念的本质表明了气体具有粘性，并给出了粗糙地面的关系式，后经尼古拉斯的实验测定，得出普适速度分布公式［14］：

式中：u为速度（m·s-1）；V*为摩阻流速（m·s-1）；h为高度（m）；v为运动粘性系数（m2·s-1）。

随后风沙物理学奠基人拜格诺将上述理论引入到风沙运动中，并提出了固定地表的风速廓线公式［15］：

式中：uh为高度h处的风速（m·s-1）；Z0为粗糙度（m）。

进一步的，可以推导出粗糙度的计算公式：

式中：Z1、Z2为2个不同的高度（m）；v1、v2分别为高度Z1和Z2处的风速（m·s-1）。

地表粗糙度反映了粗糙地表对气流的阻滞作用，表明地表不同下垫面与近地表气流的相互作用和物质能量交换，在风蚀防护评价中得到了广泛的应用［16-18］。但在近几十年的国内外研究中，科研人员大多只是采纳这一参数，并没有进行更深入的研究，甚至有学者提出现有的计算方式只是对粗糙度的估算，计算结果并不准确［19］。通过分析粗糙度的计算公式，本研究发现，粗糙度由2个不同高度的风速计算获得，而不同的高度组合会得到不同的数值，因此粗糙度数值具有不唯一性，进而只适用于当次实验，无法用于多实验间的对比分析。地表粗糙度在特定下垫面条件下并不是一个常数，而是随着风速的变化产生明显的变化，野外气流变化较为复杂，仅用2 个高度的数据进行计算所得到的误差也较大。

通过对风速廓线的研究，本研究提出利用附面层位移厚度来进行评价防护效益。附面层是指由于气流具有粘性而不可忽略的地表流动层［20］，附面层位移厚度是指相比于光滑通道，相同体积的粘性流体在相同时间内流过的通道所需加宽的厚度［21-22］，代表了粗糙地表对气流能量的消耗。将风速廓线进行积分，计算函数曲线与坐标轴间的面积，与其面积相等的矩形的长度即为附面层位移厚度（图1），进一步可推导出附面层位移厚度的计算公式［11］：

图1 附面层位移厚度示意图Fig.1 Diagram of displacement thickness of boundary layer

式中：Δd为附面层位移厚度；为风速廓线函数；a0为理想状态下的风速值（m·s-1）。

聚乳酸沙袋沙障是近年来应用较多的新材料沙障，具有绿色可降解、防护时间长、成本低、布设简单等优点，且实现了以沙治沙的理念，在风沙治理中发挥着重要的作用［23-24］。因此，本研究以聚乳酸沙袋沙障（图2）为研究对象来分析地表粗糙度和附面层位移厚度在评价防护措施防护效果时的准确性，通过风洞模拟实验，分别计算地表粗糙度和附面层位移厚度，并以风速频数、风速流场及防风效能等指标为依据来衡量这2个参数的准确性和适用性。

图2 聚乳酸沙袋沙障Fig.2 Polylactic acid sandbag barrier

1 材料与方法

1.1 风洞

本研究中应用的风洞模拟实验系统属于直流式风洞（图3a）。风洞全长约为24.3 m，其中实验段长度为12 m，实验段的横截面积为0.6 m×0.6 m；动力段长度为1.48 m，为实验提供气流；稳定段长度为1.5 m，内置蜂窝网，将大团空气打散，转变为细流；收缩段长度为1.5 m，通过横截面积的急剧收缩加速气流速度，三者可提供风速在3～40 m·s-1范围内连续可调的稳定气流。在实验中，我们利用粗糙元（4 cm×5 cm×6 cm，图3b）来拟合风速廓线，布设方式如图3c所示，并通过以下公式计算雷诺数（Re）［20］：

图3 风洞示意图Fig.3 Wind tunnel schematic diagram

式中：u为流速（m·s-1）；d为特征长度（m）；v为空气粘度（m2·s-1），取值。参考野外气象资料，本研究设置实验风速为8 m·s-（1高度1 m处），此时雷诺数为3.2×105。

为准确测量不同位置的风速，我们采用了三维移测系统，可在实验段自由移动，精度为1 mm。测量设备我们采用了热膜风速仪（IFA300，TSI，Shoreview，MN，USA）和热线风速仪（VT-200，KIMO，Montpon Ménestérol，France）。

1.2 模型制作

聚乳酸沙袋沙障常用障体高度为10 cm，障格边长为1.0 m、1.5 m和2.0 m。本研究采用聚乳酸材料灌装沙土的方式按照1:5 的比例制作模型，模型高度为2 cm，相应障格边长分别为20 cm、30 cm 和40 cm，铺成10 行3 列网格，此外为接近野外地表实际状况，在实验段底板上粘附一层沙粒。

风洞实验测点布设如图4所示。沿风向在洞体中间网格的中线布设测点，初始位置位于障前20 cm，水平方向测点间隔为障格的十等分距，垂直高度测点间隔为1 cm，垂直高度测量范围为1～20 cm。沿风向测定不同障格内风速，直至网格内风速值趋于稳定，并在下一个网格内布设水平截面测点。水平截面测点每行、列分别布设9个测点，布设高度分别为1 cm（低于障高）、2 cm（与障高平齐）和4 cm（高于障高）。

图4 实验测点分布Fig.4 Distribution of experimental measurement points

1.3 风速频数统计

统计风速稳定后障格内不同高度的风速频数，分析最小值、最大值、平均值、标准差、峰度、偏度及变异系数等指标。峰度是表征与正态分布曲线相比，概率密度分布曲线分布的陡峭程度，正值表示平均值处分布更多的数值，负值表示数值分布较为集中，计算公式为［25］：

式中：K为峰度；xi为第i个统计数值；为样本平均值；n为样本数量。

偏度是指相对于正态分布曲线，概率密度分布曲线分布的对称情况。偏度为正值表明数据更多的分布在左侧，样本中位数小于均值；偏度为负值表明数据更多的分布在右侧，样本中位数大于均值。计算公式为［26］：

式中：S为偏度；xi为第i个统计数值；为样本平均值；n为样本数量。

样本峰度和偏度可以分析样本频数分布特征，变异系数则可以分析样本的离散程度。一般认为，数值大于0.15时样本数据变异程度较大，数据可靠性较差。计算公式为［27］：

式中：Cv为变异系数；sd 为标准差；mean 为样本均值。

1.4 空间相关性

空间自相关统计量是用于度量空间数据的一个基本性质，即某位置上的数据与其他位置上的数据间的相互依赖程度，模型包括块金值、基台值、变程等参数。其中，变程指相邻2个位置具有空间相关性的最大距离，块金值与基台值的比值即为空间相关度，一般认为区域中的测点间空间相关度在0～0.25范围内具有极强的相关性，在0.25～0.75范围内具有中等程度的相关性，大于0.75时空间相关性较弱［28］。

1.5 防风效能

防风效能是指气流经过防护措施后风速下降的程度。计算公式为［29］：

式中：a为防风效能；v0为原风速（m·s-1）；v1为气流经过防护措施后的风速（m·s-1）。

2 结果与分析

2.1 地表粗糙度与附面层位移厚度

通过测量风速达到稳定后的障格内不同高度处的风速值，拟合出不同边长沙障的风速廓线（图5），各风速廓线均能较好地拟合为以自然指数为底的指数函数，且决定系数（R2）均大于0.97。

图5 风速稳定后沙障网格内风速廓线Fig.5 Wind speed profile within the sand barrier grid after stable wind speed

根据风速廓线可以计算出不同边长的障格内风速稳定后的粗糙度和附面层位移厚度（表1），其中在地表粗糙度的计算中选取了实际高度为15 cm和100 cm 处的风速值。从表1 可见，地表粗糙度随着障格边长的增加呈现出先减小后增大的趋势，边长为1.0 m 沙障网格的地表粗糙度数值最高，分别为1.5 m和2.0 m的2.98和1.80倍。与地表粗糙度不同的是，附面层位移厚度呈现出随着障格边长的增大而减小的规律。

表1 不同障格边长的沙障地表粗糙度和附面层位移厚度Tab.1 Surface roughness and boundary layer displacement thickness of sand barriers with different side lengths

2.2 风速频数统计

统计分析不同边长沙障风速稳定后的障格内风速频数（表2）可见，不同高度的风速频数均值随着障格边长的增大均呈现出上升的趋势。边长为1.5 m和2.0 m的沙障障格内风速频数均值在低于障格高度处分别是1.0 m 障格的1.28 和1.69 倍，在障格高度处分别是1.0 m 障格的1.14 和1.20 倍，在高于障格高度处分别是1.0 m 障格的1.06 和1.11 倍，整体平均分别是1.0 m 障格的1.13 和1.24 倍。从增长的幅度可以看出，随着高度的增加，风速趋向于稳定。但是数据标准差较大，表明均值受极值的影响较大。且从数据分布的峰度、偏度可以看出，风速频数离散程度较大，表明障格内的风速变化较为剧烈且复杂。此外，近半数样本的变异系数达到0.15以上，这表明风速均值受极值影响较大，需要更进一步的处理和分析。

表2 不同边长沙障障格内风速频数分布Tab.2 Frequency distribution of wind speed within different side length barriers

2.3 空间自相关性

为了进一步研究沙障防护效果和障格边长之间的关系，需要拟合障格内的风速流场。在拟合之前，为分析障格内测点之间是否具有空间相关性，利用地统计学的方法对不同边长的障格风速频数进行空间自相关分析（表3）。结果表明，不同边长的沙障网格内的风速频数均能拟合为高斯模型，且空间相关度均小于0.25，表明障内测点间具有极强的空间相关性，风洞模拟实验数据具有较高的可靠性和准确性。

表3 不同边长的沙障障格空间自相关分析Tab.3 Autocorrelation analysis of barrier grid space with different side lengths

2.4 风速流场

数值模拟不同障格内风速稳定后流场如图6所示。不同边长障格内风速均小于原风速，表明沙障起到了良好的防护作用。测量高度为5 cm 处的风速流场可以看到，边长为1.0 m 的障格在障后左右形成2个蓝色减速区，中间形成黄色加速区，核心有部分橙色区域；边长为1.5 m 障格后形成较大范围的蓝色低速区，随后风速逐渐恢复，并在障格后半部分形成2 个橙色涡流加速区域，；边长为2.0 m 障格在障后1/3范围内形成2个蓝色低速区，中后部则是大范围的橙色加速区。测量高度为10 cm的风速流场可以看到，边长为1.0 m障格后风速逐渐降低，两侧速度降低更快；边长为1.5 m 障格内形成中部低速，向两侧增大的结构；2.0 m 障格内流场结构与测量高度5 cm 时相似，但颜色变浅。测量高度为20 cm处不同边长障格均呈现出障后风速下降的趋势，且障格边长越大，风速越大。

2.5 防风效能

通过分析风速稳定后不同边长沙障障格内的防风效能（图7）可见，不同边长障格内的防风效能均大于0.4，且随着边长的增大，防风效能呈现出下降的趋势。边长为1.0 m障格内的防风效能均大于0.6，边长为1.5 m 障格内的防风效能均大于0.5，边长为2.0 m障格内防风效能均大于0.4。边长为1.0 m和1.5 m 的障格内防风效能下降较为均匀，边长为2.0 m 的障格内防风效能在0.4～0.6 范围内急剧下降，表明边长较大的网格防风蚀能力较差，网格内风速较高。

图7 不同边长沙障防风效能Fig.7 Wind prevention efficiency of sand barriers with different side lengths

2.6 相关性分析

地表粗糙度和附面层位移厚度与各指标间的相关性分析（图8）可见，地表粗糙度与障格边长和风速呈负相关关系，与防风效能呈正相关关系，但相关性不显著。附面层位移厚度与障格边长和风速呈负相关关系，与防风效能呈正相关关系，且相关性极显著，相关系数均大于0.96。

图8 各指标间相关性分析Fig.8 Correlation analysis between various indicators

3 讨论

3.1 地表粗糙度的准确性

由表2、图6 和图7 可见，聚乳酸沙袋沙障的防护效果随着障格边长的增大而减弱。从表1可以看到，附面层位移厚度这一指标准确地反映出沙障防护效果的规律性，地表粗糙度的结果并不完全准确，且图8 相关性分析也表明地表粗糙度与各参数间的相关性低于附面层位移厚度。为了验证地表粗糙度的适用性，本研究选取多个高度的风速值与实际1 m 高度处的风速组合计算地表粗糙度（表4）。从表4 可见，实际计算高度在15～20 cm 时，地表粗糙度并没有反应出真实情况，而其他高度计算出的结果则有较好的指示性。

表4 不同高度组合计算所得的地表粗糙度Tab.4 Surface roughness calculated by combining different heights

本研究进一步拟合出沙障障格沿风向的纵截面流场图（图9）。由图9可见，气流在经过沙障时出现明显的抬升加速作用，沙障上方区域内抬升气流与原空气叠加作用明显，风速变化较为剧烈，而在抬升区域之上，气流又趋于稳定。表4 中地表粗糙度异常值的计算高度正好对应于模型中气流急剧变化的3～4 cm的高度，因此可以判定气流的剧烈变化是导致地表粗糙度结果不准确的主要原因。

图9 沿风向纵截面流场Fig.9 Flow fields of longitudinal section along the wind direction

3.2 附面层位移厚度的适用性

在实际中，地表附面层的高度要远超出防护措施的有效防护高度，附面层位移厚度也超过本研究的测量范围，因此在应用于风蚀防护评价时需要进行一定的改进。在式（4）中，理想状态下的风速值（a0）无法测量，而在应用中选取超过防护措施有效防护高度的风速值即可。由图9 可见，不同边长的沙障障格内实际高度超过0.5 m 时，风速趋近于原风速，因此可以选取1.0 m 作为评价沙障防护效果时附面层位移厚度的计算高度。此外，在进行防护林、生物结皮等其他防护措施的效果评价时也可以进行相应调整。

4 结论

（1）聚乳酸沙袋沙障具有较好的防风作用，且防护效果随着沙障规格的增大而减弱。

（2）地表粗糙度是由2 个不同高度的风速计算得到，当高度选取在风速急剧变化的区间范围内时，计算得到的地表粗糙度准确性较差，因此在应用时要注意气流结构特征。

（3）附面层位移厚度是对有效防护高度内的风速廓线进行计算，弥补了区间气流急剧变化带来的影响，在评价防护措施的效果时准确性较高，可以作为地表粗糙度的辅助参考指标。