不同起沙閾值判定方案在塔克拉瑪干沙漠的適用性對比研究

楊興華， 馬明杰， 周成龍， 何 清

（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方案的理論基礎較為完善，但計算過程較為繁瑣，需要獲取的參數較多，且很多系數源于風洞實驗或外場的觀測試驗，系數和參數均存在地域適用性問題。因此，至今仍缺乏一種快速、簡便、準確的起沙閾值判定方法。

與國內外不同地區沙漠的起沙閾值比較可知，塔克拉瑪干沙漠的起沙閾值最?。ū?）［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；起沙閾值符合夏季＞秋季＞春季＞冬季的變化規律。