基于RWEQ模型的茫崖市防風固沙功能評估及敏感地類識別

王 蕾, 趙 霞, 張琛悅, 唐文家

(1.青海師范大學 青藏高原地表過程與生態保育教育部重點實驗室, 西寧 810008;2.高原科學與可持續發展研究院高原土壤信息科學研究團隊, 西寧 810008; 3.青海省生態環境監測中心, 西寧 810007)

我國是荒漠化嚴重國家，2014年全國荒漠化土地面積為261.16×104km2，其中風蝕荒漠化的面積最大且分布最廣[1]，主要分為“三北”戈壁沙漠及沙地風蝕類型區[2]。防風固沙是我國北方干旱與半干旱地區的重要生態系統功能，精確地評估防風固沙功能可為生態保護紅線劃定、生態環境的保護與修復提供科學指導。2014年國家為保障和維護區域生態安全的底線與生命線，提出了生態保護紅線戰略，其中生態功能評價是劃定生態保護紅線的前提。2017年生態環境部頒布的《生態保護紅線劃定指南》規范了生態保護紅線中生態功能的評價方法[3-4]。《生態保護紅線劃定指南》推薦采用修正風蝕模型(RWEQ)[5-6]來評估全國尺度的防風固沙功能，但未明確規定模型中氣候因子WF等關鍵參數的獲取方法。RWEQ模型在我國西北干旱半干旱地區的研究應用多集中在農田、沙地、草地混雜的農牧交錯區[7-19]，而應用于以荒漠戈壁為主的縣級地區的研究較少。青海省茫崖市位于柴達木盆地，屬于典型的荒漠戈壁區，是我國西北典型的強烈風蝕區和生態脆弱區[20]，也是氣候變化敏感區[21-22]。本研究以青海省茫崖市為研究對象，首先通過Fq代替Wf、風場強度因子計算Wf法、風場強度因子改進計算Wf法、Fq直接代替WF法這4種方法確定適宜于荒漠戈壁區的RWEQ模型參數——氣候因子WF，然后基于RWEQ模型計算出茫崖市的實際土壤風蝕量并分析其防風固沙功能的空間變異特征，最后利用指南法、自然間斷法和文獻—實際結合法這3種不同分級方法識別茫崖市防風固沙功能的敏感區。本研究旨在通過確定RWEQ模型中的氣候因子WF參數來評估茫崖市的防風固沙功能，并利用不同分級方法識別出防風固沙功能敏感的土地利用類型(簡稱地類)，為生態保護紅線的劃定提供科學指導。

1 材料與方法

1.1 研究區概況

茫崖市地處青海省西北柴達木盆地的西北角，介于37°02′—38°59′N和90°07′—93°09′E(圖1)，面積約5×104km2，海拔在2 300 m以上，屬大陸荒漠區氣候。年日照時數3 152 h，年均氣溫4.14℃，年均降水量和蒸發量分別為50 mm和2 692 mm[23-24]，年均風速5.1 m/s，最大風速可達39 m/s，8級以上大風年均日數超過100 d。研究區內土壤類型主要有風沙土、灰棕漠土、漠境鹽土、冷鈣土和粗骨土等，植被覆蓋度不足0.12 %(表1)，浮塵、沙塵暴頻發，是西北典型的極干旱風蝕區[25]。2014年，茫崖市被納入國家級“沙化土地封禁保護區”試點范圍，重點開展保護區封禁設施和固沙壓沙等生態工程建設。

圖1 研究區概況

表1 茫崖市不同區域土壤、植被及山脈概況

1.2 RWEQ模型簡介

本研究采用《生態保護紅線劃定指南》中推薦的RWEQ模型評估茫崖市的防風固沙生態功能。該模型計算的防風固沙物質量是無植被/作物(包括殘茬)覆蓋的潛在風力侵蝕量與實際風力侵蝕量的差值[26]，主要公式和因子見公式(1—7)。

SR=SL潛-SL

(1)

(2)

(3)

Qmax潛=109.8(WF×EF×SCF×K′)

(4)

S潛=150.71(WF×EF×SCF×K′)-0.3711

(5)

Qmax=109.8(WF×EF×SCF×K′×C)

(6)

S=150.71(WF×EF×SCF×K′×C)-0.3711

(7)

式中:SR為單位面積年固沙量(kg/m2)；SL潛為潛在風力侵蝕量(kg/m2)；SL為實際風力侵蝕量(kg/m2)；Qmax潛為潛在風沙最大轉移量(kg/m)；Qmax為風沙滯留量(kg/m)；S潛為潛在關鍵地塊長度(m)；S為區域關鍵地塊長度(m)；z為最大風蝕出現距離(m),通常采用下風向距離50 m[4]；WF為氣候因子(kg/m);EF為土壤可蝕因子(無量綱);SCF為土壤結皮因子(無量綱);K′為地表糙度因子；C為植被覆蓋因子。

1.2.1 氣候因子WF WF表征氣溫、降水、蒸發和風速等氣象因素對防風固沙功能的綜合影響[26]，計算見公式(8)。

(8)

式中:WF為累加12個月風力因子Wf得到的年平均氣候因子(kg/m)；ρ為空氣密度(kg/m3)；g為重力加速度(9.8 m/s2)；SW為各月的年平均土壤濕度因子；SD為雪蓋因子。

由于《生態保護紅線劃定指南》沒有給出WF的具體計算方式，本研究采用下列4種方法式來確定WF。

(1) 采用聯合國糧農組織提出的風蝕氣候因子Fq代替Wf，即Fq代替Wf法，見公式(9)。

(9)

式中:u為2 m高處的月平均風速(m/s)；u1,u2分別表示在z1,z2高度處的風速(m/s)；ETPi為月潛在蒸發量(mm)；Pi為月降水量(mm)；d為當月天數；Ti為月平均氣溫(℃)；ri為月平均相對濕度(%)。

(2) 依據風場強度因子[18,27]計算Wf，即風場強度因子計算Wf法，

Wf=u4×(u4-u3)2×Nd

(10)

式中:u3為起沙風速,本研究取5 m/s；u4為氣象站點2 m處的月均風速(m/s)；Nd為日均風速大于5 m/s的天數[4,28]。

(3) 為體現區域風速的日變化以及日風蝕量，將風場強度因子[29]公式(10)改進為公式(11)，即風場強度因子改進計算Wf法，

(11)

式中:u5為氣象站點2 m處的日均風速(m/s)；n為當月天數。

(4) 為反映氣候因子對風蝕的影響，采用聯合國糧農組織給出的風蝕氣候因子Fq直接代替氣候因子WF[30-31]，即Fq直接代替WF法。

1.2.2 土壤可蝕因子EF和土壤結皮因子SCF EF表示土壤受風力侵蝕影響的程度，由土壤質地、有機質和碳酸鈣含量等因素決定；SCF表示土壤抵抗風蝕能力的強弱[32-33]。

(12)

(13)

式中:sa為砂粒含量(%)；si為粉粒含量(%)；cl為黏粒含量(%)；OM為土壤有機質含量(%)；CaCO3為碳酸鈣含量(%)。

1.2.3 地表糙度因子K′K′指因地形影響產生的地表粗糙程度，反映了土壤風蝕的程度[34]。

(14)

式中:Kr為土壟糙度(cm),采用Smith-Carson方程計算;Crr為隨機糙度因子(cm),本次計算取0值;ΔH為距離L范圍內的海拔高程差,L為地勢起伏參數。

1.2.4 植被覆蓋因子CC表示植被覆蓋。植被覆蓋增加了地表糙度，使起沙風速增大，這在一定程度上阻礙了土壤顆粒的移動[35]，是控制土壤侵蝕的積極因素，其計算公式如下：

C=e-(ai×SC)

(15)

式中:SC為植被覆蓋度;ai為不同植被類型的系數,其中林地0.153 5,草地0.115 1,灌叢0.092 1,裸地0.076 8,沙地0.065 8,農田0.043 8。

綜上，RWEQ模型中EF，SCF，K′和C均受下墊面條件的限制，可由相應數據計算獲得。假設地面條件一致時茫崖市防風固沙功能主要受氣候因子WF的影響，而WF則需要根據實際情況確定，因此本研究采用上述Fq代替Wf、風場強度因子計算Wf法、風場強度因子改進計算Wf法以及Fq直接代替WF法這4種方式來確定WF。

1.3 數據來源與處理

氣象數據來源于中國氣象科學數據共享服務網(http:∥data.cma.cn/)，本研究選取2014年和2018年各氣象觀測站點的月均降水、氣溫、風速、濕度等數據，借助ArcGIS軟件通過空間插值得到柵格圖層；平均土壤濕度因子采用SMAP2017年逐日數據(https:∥nsidc.org/data/smap)，空間分辨率為1 km，通過ENVI軟件合成月值；雪蓋數據來源于寒區旱區科學數據中心(http:∥westdc.westgis.ac.cn)中國地區的MODIS雪蓋產品數據集；土壤質地數據(黏粒、粉粒、砂粒)來自中國1∶100萬土壤數據庫中的各類粒徑土壤顆粒含量空間分布圖，土壤有機質的圖層來自FAO中國土壤數據庫的有機質含量屬性，本研究按土壤類型的平均有機質含量計算；數據高程模型(DEM)來源于中國科學院資源環境科學數據中心(http:∥www.resdc.cn)，空間分辨率為30 m；植被覆蓋來源于“全國生態環境十年變化評估”項目通過遙感解譯獲得地表覆被數據。

所有參數因子統一采用2000國家大地坐標系統Gauss-Kruger投影，各因子重采樣按空間分辨率為250 m的柵格單元進行模型計算。

2 結果與分析

2.1 模型氣候因子確定及關鍵因素分析

圖2為通過Fq代替Wf、風場強度因子計算Wf法、風場強度因子改進計算Wf法、Fq直接代替WF法這4種方法計算獲得的Wf和WF空間分布情況。從圖2中可知，方法1即Fq代替Wf(圖2A，B)和方法3即風場強度因子改進計算Wf法(圖2E，F)獲取的WF空間分布格局相似，如茫崖市西南部祁漫塔格山和烏蘭烏珠爾山區均屬相對高值區，但這兩種算法得到的WF值域相差較大，其原因是Fq經歸一化后替代了Wf。方法2即風場強度因子計算Wf法(圖2C，D)利用風場強度因子計算的WF由于受月均風速u4和風力大于5 m/s天數Nd的影響，抹平或削弱了該區域風速的日變化，導致日風蝕量很大程度上被低估。方法4即Fq直接代替WF法(圖2G，H)中Fq直接代替WF獲得的茫崖市2014年WF空間分布由東北向西南方向減小，值域范圍為19.94～86.38 kg/m，最大值出現在冷湖地區，這與祁棟林獲取的柴達木地區WF的空間分布與值域范圍一致。因此，選擇方法4進行RWEQ模型運算更為合適[36]。

圖3表明，茫崖市EF值介于0.001～0.617。已有研究表明，地表25 mm范圍內粒徑小于0.84 mm的土壤顆粒為可蝕性顆粒[37]，而風沙土結構疏松，發育于砂性母質，砂粒(0.05～2 mm)占90%以上，其EF最高(平均為0.524，表2)，因此極易受風蝕影響。除寒凍土和寒漠土外，其他類型土壤的EF最大值均大于0.6，棕漠土EF平均值最低(為0.229)。SCF值介于0.01～1.00，棕漠土SCF平均值最小(為0.364)，所占面積較少，僅為167.1 km2，但抗風蝕能力較強，與EF的結果一致。風沙土的SCF最大(為0.788)，其抗風蝕能力最弱(表1)。K′介于0.637～1，西南部祁漫塔格山、烏蘭烏珠爾和北部邊緣阿爾金山脈、玉勒肯·塔赫達依以及東部賽什騰山等山脈的K′值較低，這是由于山脈縱橫海拔較高、地形起伏明顯，不利于發生風蝕；而低地勢、低海拔的盆地K′較高，易發生風蝕。C值介于0.005～1，其在西南角的山區較低，中、東部區較高，近90%地區C值接近1，由于茫崖市地處柴達木盆地，自然景觀多為干旱荒漠，植被稀疏，對風蝕的抑制作用弱，因此更易發生風蝕。

2.2 茫崖市實際風蝕量和防風固沙量

圖4為采用方法4(Fq直接代替WF法)的RWEQ模型計算出的茫崖市2014年風力侵蝕量和防風固沙量的空間分布情況。根據《土壤侵蝕分類分級標準(SL190-2007)》，RWEQ模型計算的茫崖市2014年實際風蝕量可分為6級(圖4和表3)。結果顯示，茫崖市2014年平均風蝕模數(單位面積實際風力侵蝕量)在57.42 kg/(m2·a)以內，平均風蝕模數總量為57.99 kg/(m2·a)，實際風蝕總量為116 814.62萬t，屬極劇烈等級(實際風蝕量為106 421.17萬t)，占實際風蝕總量的近91%，極劇烈侵蝕區面積大于3萬km2，占全市總面積的60%以上。劇烈等級的實際風蝕量次之，全市約10 %的區域平均風蝕模數為11.88 kg/(m2·a)；除極劇烈和劇烈侵蝕外，其余等級的平均風蝕模數均較小，實際風蝕量也較少。空間分布上，2014年全市極劇烈侵蝕區主要集中在中東部，主要原因是該區降雨量小蒸發量大，植被蓋度低，風場強度大，而且土壤類型主要為荒漠風沙土；風蝕強度較輕的區域主要分布在阿爾金山脈、祁漫塔格山以及烏蘭烏珠爾等西部山區。圖4與圖3相比，單位面積實際風蝕量的空間分布格局與EF和SCF的格局相似，表明EF和SCF對風蝕影響較大。

圖2 4種算法計算的風力因子和氣候因子原值

圖3 防風固沙功能關鍵影響因素

表2 茫崖市土壤可蝕因子和土壤結皮因子

圖4 茫崖市2014年實際風力侵蝕量和防風固沙量

2014年全市單位面積防風固沙量介于0～44.93 kg/m2。茫崖市90 %為未利用地，屬自然景觀和土地類型單一的荒漠戈壁區，僅西南角的少部分山區抵抗風力侵蝕和防風固沙的能力較強，整體上茫崖市防風固沙能力弱，屬風蝕嚴重區，研究結果與茫崖市的區域地帶性(西北干旱荒漠區)的實際情況相符。因此，進一步證實Fq直接代替WF法確定的模型參數，應用于干旱半干旱荒漠戈壁區的防風固沙功能研究中具有可行性。

表3 實際風蝕量等級及變化

2.3 防風固沙功能評估分級與敏感地類識別

為更細致地反映茫崖市防風固沙功能的空間分布特點，以及更好地識別出對防風固沙功能變化敏感的地類，本研究分別采用了《生態保護紅線劃定指南》推薦的三級分類方法(簡稱“指南法”)、自然間斷法和文獻—實際結合法對茫崖市防風固沙功能進行分級，分級依據見表4。指南法評估技術將防風固沙功能服務值(單位面積年固沙量SR)從高到低排序，以累加服務值占服務總值比例的50 %和80 %所對應的柵格值作為生態系統服務功能評估分級的分界點，劃分為“1等：一般重要”，“2等：重要”和“3等：極重要”三級；自然間斷法是利用ArcGIS軟件將優化后的RWEQ模型計算出的SR進行自然分組，使各等級之間的差異最大化，SR差異相對較大的柵格值設為分級裂點；文獻—實際結合法主要是依據茫崖市SR值域范圍以及實際風蝕量等級(表3)，并結合已有研究設置對應的SR分級標準[38]。同時，通過增加分級數量來更詳細地反映功能原值SR中蘊含的豐富空間分異特征。

表4 防風固沙服務功能評估分級方法

從圖5(僅顯示差異最明顯的西南角)可知，指南法、自然間斷法和文獻—實際結合法3種評估分級方法所得的茫崖市防風固沙功能空間分異具有相似的規律，2014—2018年防風固沙功能等級的空間分布特征差異較小，其中文獻—實際結合法(圖5E，F)刻畫的茫崖市西南山區的防風固沙功能識別范圍(除一般重要等級)比其他兩種方法要廣，能夠在較大程度上識別出高等級區域，避免因等級劃分低而忽視相應的保護措施，并且其極重要等級分布區域面積介于其他兩種方法之間。因此，在考慮實施生態保護措施的時間與經濟成本時，文獻—實際結合法評估分級結果具有更好的參考價值。

通過ArcGIS將土地利用類型和防風固沙功能空間分布進行疊置分析，來反映不同土地利用類型的防風固沙能力以及對功能變化敏感的地類[39]。結果表明，2014年防風固沙能力最強的土地利用類型為草地，固沙量達到1 753.9×104t，其次為未利用地和水域(主要是人工河渠、水庫等)(表5)，而城鄉居民工礦用地(主要是人工建設地面)、林地和耕地的防風固沙能力較弱，說明草地在防風固沙中起著主導作用。2014—2018年防風固沙能力變化最敏感的土地利用類型也是草地(敏感系數3.77)，防風固沙量增加了108.1×104t，原因可能是生態系統的固沙能力隨著草地覆蓋度的增加而增強(表6)。從防風固沙功能重要性等級(普適性的指南法)上看，2018年一般重要區域(1級)的主要土地利用類型是未利用地，占茫崖市面積的近90%，其次為草地、城鄉居民工礦用地和水域等；重要區域(2級)和極重要區域(3級)面積最大的土地利用類型均為草地，出現這一現象的原因是茫崖市固沙能力較強的草地面積占比相對較少，僅占7%左右，而未利用地面積雖大但保沙率卻較低。2014—2018年各土地利用類型(除耕地外)的防風固沙功能等級屬于一般重要(1級)，這與茫崖市氣候干燥、水資源貧乏的荒漠環境有關。功能等級變化最敏感的土地利用類型是水域，其面積減小導致功能1級區域占比降低。

圖5 3種防風固沙功能評估分級結果

表5 不同土地利用/覆被類型的防風固沙能力及敏感系數

3 討 論

RWEQ模型在我國西北干旱半干旱地區農田、沙地、草地及其交錯區的應用研究報道較多[7-17]，能否用于以荒漠和戈壁為主的縣級地區尚不明確，本研究表明通過選擇適宜的WF獲取方法，RWEQ模型也可用于以荒漠和戈壁為主的干旱半干旱地區。RWEQ模型的關鍵因子是氣候因子WF，方法1將歸一化后風蝕性氣候因子Fq代替風力因子Wf，這使計算結果不能反映原始值域的變化。方法2與方法3的起沙風速統一取值為5 m/s未考慮到區域土壤質地對風力搬運土壤顆粒的影響，因此本研究最終采用方法4即Fq直接代替WF法，其空間分布與值域結果與祁棟林[40]、董玉祥等[41]獲取的中國干旱與半干旱地區風蝕氣候侵蝕力極強區——柴達木盆地一致，這說明茫崖市的風蝕氣候因子主要受風速限制，荒漠景觀是強風蝕長期作用的結果，并且人類活動也能在一定程度上通過改變植被覆蓋度來影響風蝕強度[18-19]。模型中的土壤可蝕因子EF、土壤結皮因子SCF、地表糙度因子K′短期內變化不明顯，其空間分異規律與前人的研究結果相似[18]。在環境相似地區的風蝕狀況對比上，本文與江凌等[18]對青海省中度以上風蝕區的研究以及遲文峰等[14]對內蒙古高原荒漠生態系統的研究結果一致。

在功能分級中，《生態保護紅線劃定指南》規定的累計百分比分級方法能夠較好地找出防風固沙功能變化明顯的大致空間分布區域，但是由于三級分級方法不夠精確，使得目前的分級結果不能反映功能原值中蘊含的空間分異特征，降低了評估結果的政策參考價值。因此本研究探索的文獻—實際結合法不僅在分級數量上能夠較詳細地刻畫SR的空間分布特征，還能在功能原值的基礎上結合茫崖市的實際風蝕量，具有明確的功能分級依據。

茫崖市2014年的平均風蝕模數即實際風蝕量為57.99 kg/(m2·a)，與賀倩等[42]基于RWEQ模型計算出的青海三江源地區2015年實際風蝕量的最大值60.86 kg/m2相近，且最大風蝕量集中在三江源地區的西部，臨近柴達木盆地的茫崖市。本研究中茫崖市的防風固沙功能空間分布特征反映出該市整體防風固沙功能較低，這與溫豪[43]的模擬結果一致，證實了本研究結果的準確性。綜上，RWEQ模型應用于荒漠戈壁為主的茫崖市縣域研究時具有可行性，可為生態保護措施的制定與實施提供理論依據。

表6 2014-2018年土地利用/覆被類型面積變化矩陣 km2

4 結 論

(1) 利用Fq直接代替WF法獲得的空間分布和閾值范圍與實際最相符，RWEQ模型可較準確地評估茫崖市的防風固沙功能，能為縣級地區防風固沙功能評估技術的應用研究提供參考。

(2) 茫崖市2014年的平均風蝕模數為57.99 kg/(m2·a)，屬風力侵蝕嚴重區；2014年全市單位面積防風固沙量介于0～44.93 kg/m2，防風固沙功能整體較弱。

(3) 文獻—實際結合法更能詳細地反映茫崖市內部防風固沙功能的空間變異特征，且高等級功能的區域范圍較廣，西南部的防風固沙能力略強于中東部。

(4) 在耕地、林地、草地、水域、城鄉居民工礦用地、未利用土地6大類土地利用類型中，草地是茫崖市2014—2018年期間對防風固沙功能變化最敏感的土地利用類型，敏感系數最大3.77。