不同起沙阈值判定方案在塔克拉玛干沙漠的适用性对比研究

杨兴华， 马明杰， 周成龙， 何 清

（1.山西师范大学地理科学学院，山西 太原 030032；2.新疆塔克拉玛干沙漠气象国家野外科学观测研究站，新疆 乌鲁木齐 830002）

风沙运动是全球干旱-半干旱区重要的地表过程之一。风沙运动可以引起沙漠化与土壤风蚀［1-2］，造成道路、农田与村庄等掩埋［3-4］。同时，由其引起的沙尘气溶胶排放、传输与沉降还会对区域与全球气候变化［5-8］、生态环境［9-10］及人类健康［11-12］等产生重大影响。因此，风沙运动的相关研究受到地貌、气象与气候、生态与环境等学科领域的广泛关注［13］。

风是风沙运动发生的重要动力因素。只有当风速（u）或者摩擦速度（u*）达到某一临界值时，风沙运动才会发生，该值即临界起沙风速（ut）或临界起沙摩擦速度（u*t）。临界起沙风速与临界起沙摩擦速度的表征意义相同，前者可通过风速和风沙运动的观测直接确定，后者需将风速换算成摩擦速度，然后配合风沙运动观测资料确定［1］；不同学科专业的学者在使用上存在差异，地理学科的学者多使用临界起沙风速，而气象、环境学科的学者倾向于使用临界起沙摩擦速度，本研究中将二者统称为起沙阈值。起沙阈值不仅是界定风沙运动能否发生的关键指标，也是土壤风蚀模型和沙尘暴预报模式中的核心参数，决定着风蚀强度和沙尘通量的计算精度［1,14-17］。有研究表明，同等条件下起沙阈值降低32.1%和49.5%，起沙量可增加65.8%和69.2%［18］。起沙阈值的大小与土壤（如土壤湿度、土壤粒径和土壤组分等）、植被（如植被盖度等）和大气条件（如气温和空气湿度等）等密切相关［1,16-23］。当土壤湿度大于对起沙产生影响作用的临界值0.005 m3·m-3时，起沙阈值随土壤湿度增加而增大［17］；在起沙参数化方案中，起沙阈值被表述为土壤粒径的函数，当土壤粒径为74 μm 左右时对应的起沙阈值最小，小于或者大于该粒径值时，起沙阈值也越大［16-17］；土壤中因含有盐分而形成盐壳或结皮可使起沙阈值增大［21］；植被覆盖可减少沙源供给，降低作用于地表的风剪切力，进而降低起沙发生几率，增大起沙阈值［16-18］；空气湿度的增大可提高表层土壤湿度或者增加土壤颗粒间的粘滞度，进而导致起沙阈值升高［22-23］；气温的变化可以改变空气湿度，进而影响起沙阈值［22-23］。起沙阈值的判定方法已开展了较多研究，李晓岚［24］将其分为观测法包括风洞实验［1,25-27］和外场试验［1,28-29］、统计法［30-32］及参数化［1,16,33］。各种方法对比研究的相关工作开展尚少，因此，哪种方法更有优势至今尚不明确。

我国有荒漠化土地261.16×104km2，是受风沙运动危害较为严重的国家之一［34-35］。塔克拉玛干沙漠为我国第一大沙漠，风沙活动强度在世界上首屈一指［36］。同时，其所处的塔里木盆地是新疆社会、经济、文化等的重要承载区，塔克拉玛干沙漠及周边地区分布有新疆53.5%的绿洲，居住着新疆48.07%的人口，储藏有丰富的油气等矿产资源［37］。因此，开展该区域起沙阈值研究，为风沙灾害及沙尘天气发生提供预报预警，对当地的经济社会发展和生态环境保护具有重要的意义。

关于塔克拉玛干沙漠起沙阈值方面的研究工作，主要表现在如下几个方面：陈渭南等［28］在沙漠北缘肖塘地区通过人工目测风沙运动发生情况，配合同步观测的风速数据，确定2 m 高度瞬时冲击临界起沙风速为5.0 m·s-1，该值被广泛应用于塔克拉玛干沙漠风沙运动研究；Ishizuka 等［38］在沙漠南缘策勒地区通过光学传感器测得3.8 m高度临界起沙风速分别为7.5 m·s-1（干沙）和9.5 m·s-1（湿沙）；Kurosaki 等［31］基于沙漠周边气象站风速和沙尘暴资料，统计得到该区域10 m高度临界起沙风速为5.2～8.2 m·s-1；Yang 等［39］利用沙漠腹地所测沙粒跃移数据，计算得到2 m 高度起沙风速为3.5～10.9 m·s-1；Zhou 等［40］基于外场观测资料对比了几种常用起沙阈值判定方法的差异。尽管如此，塔克拉玛干沙漠起沙阈值的研究仍存在以下不足之处：观测期较短，缺少对起沙阈值季节变化的认识；起沙阈值判定方法在研究区的适用性未进行评估。本研究基于塔克拉玛干沙漠腹地野外长期观测数据，评估5种常用起沙阈值判定方法的适用性，遴选出最优判定方法，进而确定新的起沙阈值，为研究区风沙运动的判定提供更精准的判据。

1 研究区概况

观测试验在塔克拉玛干沙漠腹地塔中地区开展（39°00′N，83°38′E，海拔1103 m；图1）。塔中地区下垫面均为流沙地表；地表土壤粒径集中分布在63～250 μm之间，中值粒径为147 μm［41］。塔中地区地貌为纵向沙垄与垄间地相间分布，沙垄主要为NNE-SSW 或EN-SW 走向，沙垄与垄间地的相对高度约40～50 m，垄间地宽约1～3 km，长约2～5 km，地势相对开阔平坦［42］。本研究的试验点位于2条高大沙垄间的平坦沙地上，利用塔中气象站1996—2017年气象观测资料，得到研究区年均气温为11.7 ℃，有观测记录以来最高气温为46.0 ℃，最低气温为-32.6 ℃；年均降水量为27.6 mm，3—8 月的降水约占全年降水的87.6%；年均蒸发量达3741.8 mm；年均风速为2.2 m·s-1，年均大风日数为10.7 d，主要集中在春、夏季；年均沙尘暴、扬沙和浮尘日数分别为17.0 d、68.0 d和122.0 d，春、夏季为高发期。

2 数据与方法

2.1 风沙运动数据

借助Sensit 压电式风蚀传感器和Big Spring Number Eight（BSNE）集沙仪开展风沙运动观测（图1）。当沙粒撞击风蚀传感器的晶片时，会输出脉冲信号给数据采集器（Campbell CR1000），记录下撞击颗粒的数量；无沙粒撞击时则记录为0。数据采集频率设为1 Hz，可提供1 s时间步长撞击颗粒数的连续观测数据。当某1 s的撞击颗粒数大于0时，则记录为起沙时长1 s，进而可统计每次沙尘暴天气过程的总撞击颗粒数和起沙时长。风蚀传感器安装高度为5 cm 和10 cm，本研究仅使用5 cm 高度的测量数据。Sensit 风蚀传感器经过长期野外试验验证，对研究区起沙过程的测量效果良好［43］。为了弥补风蚀传感器不能测量沙尘输送通量的缺陷，在其一侧200～300 cm 处安装BSNE 梯度集沙塔，且2 组仪器的连线与研究区主风向垂直。BSNE 集沙塔的高度为200 cm，安装有6组集沙仪，集沙盒进沙口中间距离地面分别为5 cm、10 cm、20 cm、50 cm、100 cm和200 cm；采样频率为单次沙尘暴或扬沙天气过程。2 组仪器测量数据相互配合可以获取观测点100 cm×200 cm截面上不同时间步长的沙尘水平通量，用以验证沙尘水平通量参数化方案的计算值［44］。

图1 塔克拉玛干沙漠试验点位置及观测仪器Fig.1 Location of experimental site in the Taklimakan Desert and the instruments that were used in the field experiments

2.2 近地层微气象和沙尘浓度数据

近地层微气象借助10 m梯度探测平台开展观测（图1）。该平台架设有0.5 m、1 m、2 m、4 m和10 m 5层梯度风速（Metone 010C）、风向（Metone 020C）和空气温湿度（Vaisala HMP45C）传感器，并在3 m 高度架设了沙尘浓度测量仪（Grimm1.108）。同时该平台还设有1.5 cm、5 cm和10 cm 3层土壤湿度探头（Campbell CS616）。微气象测量数据统一存储在CR1000 型数据采集器中（Campbell），可提供1 s、1 min、30 min 和1 h 时间步长的测量数据。借助Grimm1.108 开展沙尘浓度的测量，该仪器可测量20 μm以下沙尘粒子的浓度，测量频率为6 s，然后平均为1 min和5 min的值。沙尘天气数据源于塔中气象站2008年2月—2018年3月天气现象观测记录。

2.3 起沙阈值判定方法

基于试验内容，选取Stout［29］方案代表野外观测法；Kurosaki and Mikami（KM）［31］方案、李晓岚和张宏升（LZ）［32］方案代表统计法；Marticorena and Bergametti（MB）［16］方案、Shao and Lu（SL）［33］方案代表参数化。具体如下：

（1）观测法

Stout［29］基于Sensit风蚀传感器测量的风沙运动数据，发展了起沙风速的判定方法：

式中：ut为分钟临界起沙风速（m·s-1）；uˉ为分钟平均风速（m·s-1）；σ为分钟平均风速的标准差；γ为风沙运动强度系数，用单位时间内起沙发生时间所占比表示，取值范围在0～1 之间，当γ=1 时，表示风沙运动持续发生；γ=0时，表示无风沙运动发生；Ф(γ)为γ的正态分布函数。该方法已被用于多个地区起沙阈值的判定［21,39］。

（2）统计法

Kurosaki and Mikami［31］基于风速和沙尘暴观测数据，发展了起沙阈值的判定方法：

式中：Pi为百分比；ni为风速等级为i时风沙运动发生的频次；Ni为风速等级为i的总频次；风速等级i的分组间隔为0.2 m·s-1；Pi为50%时的风速等级定义为临界起沙风速ut。

李晓岚和张宏升［32］基于内蒙古科尔沁沙地起沙试验数据，将一定观测高度沙尘浓度开始迅速且连续增大时对应的风速判定为临界起沙风速。

（3）参数化

Marticorena and Bergametti［16］，Shao and Lu［33］发展的起沙阈值参数化方案是当前沙尘暴模式中使用最为广泛的方案。方案均可用下式表达：

式中：u*ts(d)为光滑地表临界起沙摩擦速度；fλ(λ)为地表粗糙元修正方程；fw(w)为土壤湿度修正方程。B=adx+b，a、b和x分别为1331 cm-x、0.38 和1.56；K=［(σgd/ρ)0.5(1+0.006/σg(d)2.5)］0.5，g为重力加速度，取值9.81 m·s-2，σ为沙粒密度，取值2650 kg·m-3，ρ为空气密度，取值1.02 kg·m-3，d为沙粒粒径，取值147 μm；z0为地表粗糙度，取值0.2478 cm；z0s为光滑地表粗糙度，取值0.01372 cm；w为土壤重量含水量；w′为土壤重量含水量是否影响起沙的临界值［w′=0.0014(%clay)2+0.17(%clay)］［18］。

在SL方案中，

式中：AN、ε、mr、σr和βr为系数，取值分别为0.0123、0.165 g·s-2、0.5、1 和90；λ=-0.35ln(1-a)，a为植被覆盖指数，取值0；w为土壤体积含水率。

由于2 组方案计算的均为临界起沙摩擦速度u*t，可利用下式将其换算为ut，以方便与其他方法确定的起沙阈值进行比较：

式中：k为常数，取值0.4；z为测量高度，取值2.0 m。

2.4 沙尘水平通量确定方法

沙尘水平通量利用Owen［45］方案计算：

式中：Qi(d)为粒径为i的粒径组沙尘水平通量（kg·m-1·min-1），与粒径组i的起沙阈值对应，由于笔者试验获取的起沙阈值代表了整体粒径组，因此Qi(d)即为Q；E为可发生起沙的地表比重，研究区取值1；c为系数，取值0.8。

3 结果与分析

3.1 起沙阈值判定方法适用性评估

根据试验数据的完整性，本研究选择2009 年7月17日的沙尘暴过程为例，验证各起沙阈值判定方法在研究区的适用性，图2 给出了本次沙尘暴过程风速与风沙运动变化情况，由图2可知，本次沙尘暴过程中风速从4:00—5:00 左右开始增大，在6:00 左右达到了5.0 m·s-1，此时已可监测到较弱的起沙活动；风速在9:00左右增加到了6.0 m·s-1以上，风沙运动强度迅速增大，持续至19:00左右，随着风速的降低，风沙运动减弱至结束。本次沙尘暴过程共发生风沙运动735 min，对应的最小分钟平均风速为2.5 m·s-1，最大为11.2 m·s-1，均值为7.5 m·s-1。

图2 2009年7月17日沙尘暴过程风沙运动与风速变化Fig.2 Changes of wind-blown sand movement and wind speeds during sandstorm process on July 17,2009

图3分别给出了Stout、KM、LZ、MB和SL方案确定的临界起沙风速，Stout方案确定的分钟临界起沙风速的变化范围为3.1～9.1 m·s-1，均值为6.4 m·s-1；分钟临界起沙摩擦速度的变化范围为0.19～0.54 m·s-1，均值为0.38 m·s-1；以分钟临界起沙风速为风沙运动判据，本次沙尘暴过程风沙运动发生时长203 min，为实际发生时长的27.6%；以均值为判据，风沙运动发生时长517 min，为实际发生时长的70.3%。根据KM 的判定方法，本次沙尘暴过程中临界起沙风速约为4.8 m·s-1，临界起沙摩擦速度约为0.29 m·s-1；以该值为风沙运动判据，本次沙尘暴过程风沙运动发生时长756 min，为实际发生时长的102.9%。LZ方案判定的临界起沙风速约为6.8 m·s-1，临界起沙摩擦速度约为0.41 m·s-1；以该值为判据，本次沙尘暴过程风沙运动发生时长457 min，为实际发生时长的62.2%。由于研究区下垫面和土壤湿度变化极小，MB方案计算的分钟临界起沙风速几乎无变化，值约为6.6 m·s-1，临界起沙摩擦速度约为0.39 m·s-1；以该值为风沙运动判据，本次沙尘暴过程风沙运动发生时长495 min，为实际发生时长的67.3%。SL方案计算的分钟临界起沙风速变化范围为5.4～6.0 m·s-1，均值为5.6 m·s-1；分钟临界起沙摩擦速度变化范围为0.32～0.36 m·s-1，均值为0.33 m·s-1；以分钟临界起沙风速为判据，本次沙尘暴过程风沙运动发生时长614 min；以均值为判据，风沙运动发生时长594 min，分别为实际发生起沙时长的83.5%和80.8%。上述结果表明，除了KM 方案，其余4种方案均不同程度地低估了风沙运动时长。

图3 2009年7月17日沙尘暴过程Stout、KM、LZ、MB和SL方案确定的临界起沙风速Fig.3 Threshold velocities for Stout,KM,LZ,MB and SL schemes during sandstorm process on July 17,2009

为了进一步验证各起沙阈值判定方案在研究区的适用性，本研究基于各方案确定的临界起沙阈值（Stout 和SL 方案为均值），结合Owen［45］沙尘水平通量计算方案，计算了2009年7月17日沙尘暴过程分钟沙尘水平通量，并与实测结果进行了对比（图4）。本次沙尘暴过程实测沙尘水平通量为271.0 kg·m-1，基于Stout、KM、LZ、MB 和SL 方案确定的起沙阈值计算的沙尘水平通量分别为178.3 kg·m-1、314.1 kg·m-1、122.9 kg·m-1、179.2 kg·m-1和244.7 kg·m-1，分别为实测值的65.8%、115.8%、45.3%、66.1%和90.2%。与风沙运动时长相似，基于KM方案确定的起沙阈值计算的沙尘水平通量在一定程度上高估了本次沙尘暴过程中风沙运动强度，其余4 种方案均不同程度地低估了本次沙尘暴过程中风沙运动强度。相关性分析表明，KM方案效果最优，相关系数r为0.90，依次为SL方案（r=0.87）、Stout方案（r=0.86）、MB 方案（r=0.82）、LZ 方案（r=0.80）。研究表明，脉动风速是导致沙粒运动发生的重要动力机制，且脉动风速与起沙量具有较好的一致性［46］。本研究中使用的风速数据时间步长为1 min，这在一定程度上会屏蔽掉部分风速脉动信息，进而导致估算的起沙时长与沙尘水平通量与实际情况存在一定差异。利用上述5 种起沙阈值判定方案，进一步统计计算了2009 年研究区其他8 次沙尘暴天气过程起沙时间与沙尘水平通量（表1 和表2），结果与2009 年7 月17 日沙尘暴过程一致。综上所述，KM方案最为适合研究区起沙阈值的判定。由于土壤和植被等下垫面条件存在区域差异，KM 方案在其他下垫面（如荒漠、戈壁、干涸河/湖床等）是否仍然为起沙阈值的最优判定方案有待验证，因此，未来有必要开展不同下垫面、不同起沙阈值判定方案的对比研究工作。

图4 2009年7月17日沙尘暴天气过程观测与计算的沙尘水平通量比较Fig.4 Comparisons of observed and calculated horizontal dust fluxes with Stout,KM,LZ,MB and SL schemes during sandstorm process on July 17,2009

表1 塔克拉玛干沙漠不同沙尘暴天气过程观测与计算的起沙时间比较Tab.1 Occurrence times of sand movement from observed and calculated with emission thresholds determining methods during different sandstorms in the Taklimakan Desert

表2 塔克拉玛干沙漠不同沙尘暴天气过程观测与计算的沙尘水平通量比较Tab.2 Horizontal dust fluxes from observed and calculated with emission thresholds determining methods during different sandstorms in the Taklimakan Desert

3.2 塔克拉玛干沙漠起沙阈值

利用KM方法确定了研究区2008年3月—2018年2月101次沙尘天气过程的起沙阈值（图5），结果表明：2 m高度临界起沙风速变化范围为4.0～6.0 m·s-1，均值为5.0 m·s-1；临界起沙摩擦速度变化范围为0.24～0.36 m·s-1，均值为0.30 m·s-1；临界起沙风速分布较为集中，73.2%的值位于4.6～5.2 m·s-1之间。起沙阈值具有显著的季节变化，春、夏、秋、冬季的临界起沙风速分别为4.9 m·s-1、5.4 m·s-1、5.1 m·s-1和4.6 m·s-1，临界起沙摩擦速度分别为0.29 m·s-1、0.32 m·s-1、0.30 m·s-1和0.27 m·s-1，符合夏季＞秋季＞春季＞冬季的变化规律（表3）。前人的研究表明，土壤湿度、空气湿度和气温的升高均能使起沙阈值增大［17,22-23］。研究区土壤湿度、水汽压和气温的季节变化与起沙阈值一致，使得起沙阈值呈现上述季节变化。由于塔克拉玛干沙漠降水稀少，土壤湿度与空气湿度的季节差异较小，导致该区域起沙阈值季节变化较微弱。Kurosaki and Mikami［31］利用塔克拉玛干沙漠周边气象站1988—2005 年风速和沙尘暴观测数据，统计得到该区域10 m高度临界起沙风速的变化范围为5.2～8.2 m·s-1，均值为6.7 m·s-1；本研究的临界起沙风速换算成10 m 高度为5.0～7.6 m·s-1，均值为6.3 m·s-1，略低于Kurosaki and Mikami 的统计结果。分析其原因，Kurosaki and Mikami 使用的数据为沙尘暴时时风速，而沙尘暴是由强烈的风沙运动引起的，其对应的风速必然高于临界起沙风速；本研究的结果是基于5 cm高度风沙运动数据与同步的风速数据统计得出，可更加真实地反映研究区起沙阈值的变化情况。

表3 塔克拉玛干沙漠起沙阈值季节变化Tab.3 Seasonal variations of threshold velocities and friction velocities in the Taklimakan Desert

图5 研究区临界起沙风速频率累积分布Fig.5 Cumulative distribution of threshold velocities in study area

4 讨论

本研究基于野外观测数据评估了5种常用起沙阈值判定方法，虽然上述5 种方法均在世界不同地区的风沙运动研究中得以应用，然而将各方法用于同一区域起沙阈值的判定工作开展较少。通过比较，发现各方法之间结果的差异性较为显著，为风沙运动的判定带来较大不确定性，因此开展各方法的适用性评估具有必要性。本研究中，Stout方案利用野外风沙运动资料可以快速获取起沙阈值，但该方法是基于风速脉动的正态分布建立的，在真实野外状态下风速脉动的分布是否符合正态分布有待商榷；KM 方案选择将某一风速发生风沙运动的概率达到50%时定为临界起沙风速，从统计学的角度具有一定的合理性，但是该方案需要大量样本数据；LZ方案基于沙尘浓度变化（时间）来判定起沙阈值，具有一定的主观性；MB和LS方案的理论基础较为完善，但计算过程较为繁琐，需要获取的参数较多，且很多系数源于风洞实验或外场的观测试验，系数和参数均存在地域适用性问题。因此，至今仍缺乏一种快速、简便、准确的起沙阈值判定方法。

与国内外不同地区沙漠的起沙阈值比较可知，塔克拉玛干沙漠的起沙阈值最小（表4）［20,24,31,47-49］，且变化范围较窄。虽然降水、植被覆盖条件与敦煌沙地类似，但仍略低于敦煌沙地7.0 m·s-1和0.5 m·s-1的起沙阈值。与其他地区的沙漠相比，塔克拉玛干沙漠降水量较小，气候干燥，且无植被覆盖，为风沙运动的发生提供了有利的气候与下垫面条件；塔克拉玛干沙漠地表土壤粒径以细砂和极细砂为主，平均粒径在100 μm 左右，正处于起沙阈值的低值区［16,33,50］，为起沙提供了有利的土壤条件，上述因素是塔克拉玛干沙漠起沙阈值较小的形成原因。

表4 不同沙漠/沙地起沙阈值比较Tab.4 Comparisons of threshold velocities and friction velocities of different deserts and sands

5 结论

通过野外观测获取的风沙运动时长和沙尘水平通量数据评估了5种常用起沙阈值判定方案在塔克拉玛干沙漠腹地的适用性，并基于最优判定方案确定了塔克拉玛干沙漠腹地的起沙阈值。具体结论如下：

（1）Stout、LZ、MB 和SL 4 组起沙阈值判定方案均不同程度高估了研究区起沙阈值，使用上述起沙阈值确定研究区风沙运动时长和沙尘水平通量时将存在不同程度的低估现象；基于KM 方案判定的起沙阈值确定的风沙运动时长与实测值最为接近，沙尘水平通量与实测值接近且相关性最好。由此可以认定KM方案为判定塔克拉玛干沙漠腹地起沙阈值的最优方案。

（2）利用KM 方案统计得到塔克拉玛干沙漠腹地2 m高度临界起沙风速变化范围为4.0～6.0 m·s-1，均值为5.0 m·s-1；临界起沙摩擦速度变化范围为0.24～0.36 m·s-1，均值为0.30 m·s-1；起沙阈值符合夏季＞秋季＞春季＞冬季的变化规律。