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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

圖1 研究區(qū)概況

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

1.2 RWEQ模型簡介

本研究采用《生態(tài)保護紅線劃定指南》中推薦的RWEQ模型評估茫崖市的防風固沙生態(tài)功能。該模型計算的防風固沙物質(zhì)量是無植被/作物(包括殘茬)覆蓋的潛在風力侵蝕量與實際風力侵蝕量的差值[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潛為潛在風沙最大轉(zhuǎn)移量(kg/m)；Qmax為風沙滯留量(kg/m)；S潛為潛在關(guān)鍵地塊長度(m)；S為區(qū)域關(guān)鍵地塊長度(m)；z為最大風蝕出現(xiàn)距離(m),通常采用下風向距離50 m[4]；WF為氣候因子(kg/m);EF為土壤可蝕因子(無量綱);SCF為土壤結(jié)皮因子(無量綱);K′為地表糙度因子；C為植被覆蓋因子。

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

(8)

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

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

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

(9)

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

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

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

(10)

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

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

(11)

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

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

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

(12)

(13)

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

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

(14)

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

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

C=e-(ai×SC)

(15)

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

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

1.3 數(shù)據(jù)來源與處理

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

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

2 結(jié)果與分析

2.1 模型氣候因子確定及關(guān)鍵因素分析

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

圖3表明，茫崖市EF值介于0.001～0.617。已有研究表明，地表25 mm范圍內(nèi)粒徑小于0.84 mm的土壤顆粒為可蝕性顆粒[37]，而風沙土結(jié)構(gòu)疏松，發(fā)育于砂性母質(zhì)，砂粒(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的結(jié)果一致。風沙土的SCF最大(為0.788)，其抗風蝕能力最弱(表1)。K′介于0.637～1，西南部祁漫塔格山、烏蘭烏珠爾和北部邊緣阿爾金山脈、玉勒肯·塔赫達依以及東部賽什騰山等山脈的K′值較低，這是由于山脈縱橫海拔較高、地形起伏明顯，不利于發(fā)生風蝕；而低地勢、低海拔的盆地K′較高，易發(fā)生風蝕。C值介于0.005～1，其在西南角的山區(qū)較低，中、東部區(qū)較高，近90%地區(qū)C值接近1，由于茫崖市地處柴達木盆地，自然景觀多為干旱荒漠，植被稀疏，對風蝕的抑制作用弱，因此更易發(fā)生風蝕。

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

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

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

圖3 防風固沙功能關(guān)鍵影響因素

表2 茫崖市土壤可蝕因子和土壤結(jié)皮因子

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

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

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

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

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

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

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

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

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

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

3 討 論

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

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

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

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

4 結(jié) 論

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

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

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

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