侯麗麗, 銀 山,2, 都瓦拉, 王艷琦, 彭秀清, 雅 茹

(1.內(nèi)蒙古師范大學(xué) 地理科學(xué)學(xué)院, 呼和浩特 010022; 2.內(nèi)蒙古自治區(qū)遙感與地理信息系統(tǒng)重點實驗室,呼和浩特 010022; 3.內(nèi)蒙古自治區(qū)生態(tài)與農(nóng)業(yè)氣象中心, 呼和浩特 010051)

凈初級生產(chǎn)力(NPP)指植物在單位面積單位時間內(nèi)所積累的有機(jī)物的數(shù)量，是由光合作用所產(chǎn)生的有機(jī)質(zhì)總量中扣除自養(yǎng)呼吸后的剩余部分，是維持草原生態(tài)系統(tǒng)的物質(zhì)基礎(chǔ)[1]。沙地植被NPP不僅能反映沙地生境的生產(chǎn)能力和質(zhì)量狀況，也是沙地生態(tài)系統(tǒng)結(jié)構(gòu)和功能的重要標(biāo)志[2]。沙地植被NPP的大小，對沙地碳循環(huán)的合理調(diào)控和預(yù)測生態(tài)系統(tǒng)未來變化趨勢等具有重要的理論和現(xiàn)實意義。

傳統(tǒng)的NPP實地估測方法無法滿足區(qū)域或全球尺度的測量，NPP研究手段發(fā)展成基于遙感數(shù)據(jù)的模型估測。在眾多的模型中，光能利用率CASA模型是非常有效且有潛力的研究手段[3]。1993年P(guān)otter等[4]改進(jìn)CASA模型方程進(jìn)行全球尺度NPP的估算，中國學(xué)者也利用CASA模型估算了不同尺度的植被NPP，在大尺度上對中國陸地生態(tài)系統(tǒng)進(jìn)行模擬[5]，中小尺度上包括對藏北地區(qū)、內(nèi)蒙古草地、黑河流域等進(jìn)行NPP模擬分析[6-8]，如張錦水等[9]利用CASA模型估算中國內(nèi)蒙古的植被NPP。然而利用CASA模型對干旱區(qū)沙地、荒漠生態(tài)系統(tǒng)NPP的研究很少，且模型參數(shù)模擬相對較少。

渾善達(dá)克沙地是我國四大沙地之一，也是頻繁襲擊我國華北平原乃至全國的沙塵暴災(zāi)害最主要發(fā)源區(qū)之一[10]。整個沙地生態(tài)系統(tǒng)植被NPP模擬研究特別薄弱，在生態(tài)工程實施后的大背景下，我們對沙地NPP的發(fā)展趨勢認(rèn)識不全面，而且目前對渾善達(dá)克沙地的研究也大多局限在2010年以前關(guān)于荒漠化和氣候變化的研究[11]，很難準(zhǔn)確判斷當(dāng)前沙地生態(tài)系統(tǒng)的收支情況和生態(tài)工程效益。

因此，本文以渾善達(dá)克沙地為研究對象，開展其植被NPP模擬估算，構(gòu)建沙地2000—2016年的NPP時空數(shù)據(jù)集，并結(jié)合線性回歸方法描述沙地的時空分布特征，對于全面了解生態(tài)治理工程后的植被發(fā)展趨勢，并指導(dǎo)未來沙地生態(tài)系統(tǒng)的生態(tài)保護(hù)和恢復(fù)工作具有重要意義。

1 材料與研究方法

1.1 研究區(qū)概況

渾善達(dá)克沙地位于內(nèi)蒙古中部錫林郭勒草原南側(cè)(41°46′—44°24′N，112°22′—117°57′E)，地勢西南高，東北低，平均海拔1 300 m[11]。分布在錫林郭勒盟的蘇尼特左旗、蘇尼特右旗、阿巴嘎旗、錫林浩特市、鑲黃旗、正鑲白旗、正藍(lán)旗、多倫縣、二連浩特市少部分和赤峰市的克什克騰旗10個旗縣(市)(圖1)。氣候溫和，屬中溫帶大陸性氣候，年平均氣溫為1.5℃，1月份平均氣溫-18.3℃，7月份平均氣溫18.7℃，極端最高溫度35.9℃，極端最低氣溫-36.6℃。全年降雨量為365.1 mm，而且主要集中在7月、8月、9月份，約占全年降雨量的80%～90%。植被類型以灌木、半灌木草原、草原化荒漠和草原地帶的沙地植被為主體，其生態(tài)系統(tǒng)對外界干擾非常敏感、脆弱[12]。

圖1 渾善達(dá)克沙地地理位置

1.2 數(shù)據(jù)來源與預(yù)處理

氣象數(shù)據(jù)來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http:∥cdc.cma.gov.cn)提供的2000—2016年沙地境內(nèi)10個旗縣氣象站點的月平均溫度數(shù)據(jù)。太陽輻射數(shù)據(jù)來自于內(nèi)蒙古氣象局提供的覆蓋全內(nèi)蒙古的8個太陽總輻射站2000—2016年的觀測數(shù)據(jù)。對氣象數(shù)據(jù)采用克里金插值法進(jìn)行空間插值，獲取與NDVI數(shù)據(jù)投影、像元大小相同的250 m分辨率的柵格數(shù)據(jù)。太陽輻射數(shù)據(jù)通過反距離插值法對全內(nèi)蒙古的8個太陽輻射站的月值數(shù)據(jù)插值為1 000 m的柵格數(shù)據(jù)，然后進(jìn)行數(shù)據(jù)重采樣、掩膜獲取與NDVI數(shù)據(jù)分辨率相同的柵格數(shù)據(jù)。

遙感數(shù)據(jù)來自NASA對地觀測系統(tǒng)數(shù)據(jù)(http:∥modis.gsfc.nasa.gov/)共享平臺。本文使用涵蓋渾善達(dá)克沙地16 d合成的MOD13Q1數(shù)據(jù)產(chǎn)品以及8 d合成的MOD09A1的數(shù)據(jù)產(chǎn)品。MOD13Q1產(chǎn)品的空間分辨率為250 m，因內(nèi)蒙古大部分地區(qū)的植物在冬天不再生長或者被冰雪覆蓋等原因，選取生長季(4—10月)分析NPP的變化趨勢，運用最大值合成法(MVC)，生成月NDVI數(shù)據(jù)，降低大氣等外界因素的干擾[12]。選用分辨率為500 m的MOD09A1產(chǎn)品反演地表濕潤指數(shù)(LSWI),LSWI是近紅外波段(ρnir)和短波紅外波段(ρswir)的歸一化指數(shù)[13]，并根據(jù)地表濕潤指數(shù)來估算水分脅迫因子。將獲取的月NDVI數(shù)據(jù)和插值得到的月氣象數(shù)據(jù)、月太陽輻射數(shù)據(jù)通過CASA模型求出月NPP數(shù)據(jù)，并將每年4—10月的7個月數(shù)據(jù)進(jìn)行累加獲取2000—2016年NPP的時空數(shù)據(jù)集。

1.3 研究方法

1.3.1 CASA模型估算NPP概述 CASA模型通過植物所吸收的光合有效輻射(APAR)和實際光能利用率(ε)來估算區(qū)域的植被NPP值[14]。模型結(jié)構(gòu)一般由公式(1)—(3)表示：

NPP=APAR×ε

(1)

APAR=SOL×FPAR×0.5

(2)

ε=Tε1×Tε2×Wε×εmax

(3)

APAR由太陽總輻射(SOL)和植被對光合有效輻射的吸收比例(FPAR)來計算，常數(shù)0.5表示植被所能利用的太陽有效輻射(波長0.4～0.7 μm)占太陽總輻射的比例。ε由植被最大光能利用率(εmax)和對εmax產(chǎn)生影響的溫度脅迫因子(Tε)和水分脅迫因子(Wε)來計算[15]。

SOL由內(nèi)蒙古太陽輻射數(shù)據(jù)獲取，F(xiàn)PAR根據(jù)MOD13Q1數(shù)據(jù)產(chǎn)品的NDVI反演獲取，計算公式參考文獻(xiàn)[12]，溫度脅迫因子Tε1和Tε2的計算方法詳見文獻(xiàn)[14]。

水分脅迫因子Wε以往的計算方法通過對實際蒸散量和潛在蒸散量模型或者土壤濕度函數(shù)進(jìn)行估算，這些計算方法考慮的降水?dāng)?shù)據(jù)和土壤模型的經(jīng)驗系數(shù)較多，模擬計算的Wε精度難以保證，本文采用遙感反演地表濕潤指數(shù)(LSWI)的方法進(jìn)行Wε的計算(公式4)，可以較大程度地降低計算過程的不確定性，更好地反映地表的濕潤情況，從而提高Wε的精確性[16]。

(4)

式中：LSWImax為每個像元中全年最大的地表濕潤度指數(shù)。LSWI的取值范圍為-1～1，Wε的取值范圍為0～1。

植被最大光能利用率εmax采用朱文泉等[17]中國典型植被最大光能利用率模擬的結(jié)果作為參數(shù)，渾善達(dá)克沙地的植被類型主要是草原和荒漠植被，因此將草原和荒漠植被的模擬結(jié)果0.542 g/MJ作為研究區(qū)的最大光能利用率。

1.3.2 一元線性回歸分析法 采用一元線性回歸分析法計算2000—2016年研究區(qū)NPP的整體變化趨勢，擬合NPP與時間的線性變化趨勢(公式5)，用斜率來反映多年NPP值的變化情況[1]。

(5)

式中：slope為趨勢斜率；n為年序列時間跨度，本文為17；i為時間變量，i=1～17；NPPi為第i年的NPP值。

1.3.3 皮爾遜相關(guān)系數(shù)法 采用皮爾遜相關(guān)系數(shù)定量描述渾善達(dá)克沙地NPP的變化規(guī)律(公式6)，該系數(shù)可以表示NPP的長期趨勢變化的程度，是m個時刻的NPP與自然數(shù)列1～m的相關(guān)系數(shù)[16]。

(6)

2 結(jié)果與分析

2.1 NPP模擬結(jié)果

以2016年8月為例，依據(jù)CASA模型建模思路，通過對SOL，F(xiàn)PAR，Tε，Wε眾因子的估算，得出2016年8月的渾善達(dá)克沙地NPP取值范圍為0.057～184.098 g/m2，平均值為73.434 g/m2(圖2)。

圖2 渾善達(dá)克沙地2016年8月NPP估算結(jié)果

2.2 模型估算結(jié)果的驗證

由于本文缺乏實測數(shù)據(jù)，故采用間接驗證的方式對估算結(jié)果進(jìn)行驗證。將本文CASA模型估算的NPP值與其他人估算的NPP的均值或部分?jǐn)?shù)據(jù)進(jìn)行對比驗證[12,15,18-22]，來檢測估算數(shù)據(jù)的準(zhǔn)確性與可行性(表1)。

表1 渾善達(dá)克沙地NPP模擬值與其他模擬結(jié)果比較

注：“—”表示參考文獻(xiàn)中無月值實測數(shù)據(jù)。

本文估算的2000—2016年渾善達(dá)克沙地植被NPP多年均值為282.42 g/(m2·a)，與表1對比分析，估算結(jié)果介于內(nèi)蒙古地區(qū)植被NPP的模擬結(jié)果之間。且7月、8月份NPP均值為78.31，67.31 g/m2，與李剛等[18]、穆少杰等[20]、元志輝[12]7月、8月計算的實測數(shù)據(jù)相符合。此外，由于生態(tài)治理工程的推進(jìn)，2010年后研究區(qū)內(nèi)植被NPP值不斷增大，本研究結(jié)果高于元志輝[12]估算的渾善達(dá)克沙地NPP值，該估算結(jié)果在一個合理的區(qū)間內(nèi)，因此可以用于后續(xù)的結(jié)果與分析。

2.3 NPP年際變化特征

利用CASA模型估算的渾善達(dá)克沙地17 a生長季NPP的變化趨勢，其結(jié)果表明沙地生長季NPP呈波動上升趨勢，年均增長率為2.136 2 g/(m2·a)，從2000年的276.82 g/(m2·a)，增長到2016年的292.7 g/(m2·a)。在2000—2010年內(nèi)，NPP值波動變化明顯，2003年、2004年、2006年出現(xiàn)高值，2002年、2005年、2007年出現(xiàn)低值。而在2010年以后，每年NPP值不斷增大，在2013年達(dá)到峰值，NPP值為339.96 g/(m2·a)(圖3A)。

圖3 2000-2016年生長季NPP年際變化和偏差分析

偏差分析法直觀展示每年生長季NPP平均值與17年間NPP多年均值的偏離情況[22]，結(jié)果表明：2000—2016年，生長季NPP的偏差值呈上升趨勢。2010年以前，除2003年、2004年、2006年以外，其余年份的生長季NPP值均低于多年NPP均值。2010年后，每年的生長季NPP值均大于多年均值，2014年NPP值低于2013年，但總體處于上升狀態(tài)(圖3B)，說明渾善達(dá)克沙地植被狀況得到改善和好轉(zhuǎn),與杭玉玲等[23]關(guān)于錫林郭勒草原植被覆蓋狀況轉(zhuǎn)好的研究結(jié)果趨于一致。

2.4 NPP季節(jié)與月變化特征

從季節(jié)尺度分析，渾善達(dá)克沙地2000—2016年春季NPP值呈下降趨勢，夏季和秋季均呈上升趨勢(圖4)。春季氣溫升高,植物開始返青，NPP值逐漸增大，渾善達(dá)克沙地春季NPP值波動較大，2003年、2009年、2013年NPP值達(dá)到高值期，其中2003年達(dá)到峰值，值為57.5 g/(m2·a)，研究時間內(nèi)多數(shù)年份NPP值低于平均水平，NPP值以0.158 8 g/(m2·a)的速率緩慢下降(圖4A)，比較符合近年來內(nèi)蒙古植被返青期推遲的趨勢[24]。夏季，氣溫、降水、太陽輻射等氣候因子都為植物的生長提供良好的環(huán)境，NPP值達(dá)到最高，以2.364 6 g/(m2·a)的速率增長，2010年以后NPP上升趨勢呈指數(shù)型增長，到2013年達(dá)到峰值，值為248.83 g/(m2·a),2014年以后又呈現(xiàn)下降趨勢，但總體水平還是呈增長狀態(tài)，這也與全年NPP變化趨勢相一致(圖4B)。秋季，氣溫降低降水減少，植被開始枯黃凋零，NPP值減小。2000—2016年，秋季NPP值呈現(xiàn)波動上升趨勢，每年以0.150 3 g/m2的速率緩慢增長，秋季NPP值在2000年、2004年、2012年、2013年出現(xiàn)明顯的高值期(圖4C)。

從月變化來看，渾善達(dá)克沙地17 a內(nèi)月平均NPP值變化明顯(圖4D)。由于不同月份的水熱組合不同，各月的NPP值也各不相同，4—7月NPP值呈上升趨勢，到7月達(dá)到峰值，值為78.31 g/m2，8—10月份，NPP值呈下降趨勢，10月份，NPP值降到生長季最低值，夏季6—8月NPP值整體處于高值。其中，7月、8月份積累的NPP值占總量的51.57%，4月和10月的NPP值僅占總量的4.09%，3.05%。

圖4 NPP季節(jié)與月變化特征

2.5 NPP的空間分布特征

將2000—2016年的NPP值取平均值來反映渾善達(dá)克沙地17 a來NPP的整體空間分布特征，渾善達(dá)克沙地NPP值整體呈由西向東逐漸遞增的趨勢(附圖2)，東部的克什克騰旗、多倫縣和正藍(lán)旗處于高值區(qū)，且大部分地區(qū)NPP值高于350 g/(m2·a)。西部的二連浩特市、蘇尼特右旗、蘇尼特左旗和鑲黃旗處于低值區(qū)，占地面積為研究區(qū)的46.4%，NPP值均低于250 g/(m2·a)。從不同等級NPP的空間分布來看，NPP值為0～50 g/(m2·a)的面積占研究區(qū)的0.28%，主要分布在蘇尼特右旗、蘇尼特左旗、阿巴嘎旗、正鑲白旗和正藍(lán)旗的湖泊水域部分。NPP值大于650 g/(m2·a)的面積占沙地面積的0.25%，主要分布在克什克騰旗的東部以及東南部。NPP值介于250～350 g/(m2·a)的區(qū)域占總面積的26.06%，屬于研究區(qū)的均值區(qū)，分布在阿巴嘎旗的中西部地區(qū)、正鑲白旗的南部區(qū)域、正藍(lán)旗的西部地區(qū)和錫林浩特市的南部區(qū)域。NPP值介于350～450 g/(m2·a)的面積占沙地總面積的17.65%，集中分布在正藍(lán)旗、錫林浩特市北部、克什克騰旗和多倫縣。NPP值為450～650 g/(m2·a)的高值區(qū)占總面積的9.35%，集中分布在沙地的克什克騰旗和多倫縣。

2.6 NPP的空間變化趨勢

為進(jìn)一步分析渾善達(dá)克沙地植被NPP的空間變化規(guī)律，本文采用一元線性回歸法，基于像元尺度分析2000—2016年生長季NPP的空間變化趨勢(附圖3A)，結(jié)果表明：沙地內(nèi)78.2%的面積呈現(xiàn)增長趨勢，21.8%的面積呈現(xiàn)降低趨勢。NPP上升最明顯的區(qū)域主要分布在多倫縣、正藍(lán)旗東南部、阿巴嘎旗北部以及錫林浩特市，變化速率超過15 g/(m2·a)，研究區(qū)大多數(shù)區(qū)域NPP值的增長趨勢為0～5 g/(m2·a)。蘇尼特右旗東部、蘇尼特左旗中部、錫林浩特市以及克什克騰旗中部區(qū)域NPP值呈現(xiàn)下降趨勢，其中克什克騰旗中部地區(qū)植被NPP下降十分明顯，值介于-10～-5 g/(m2·a)。

通過對皮爾遜相關(guān)系數(shù)進(jìn)行顯著性檢驗來分析渾善達(dá)克沙地植被NPP的變化趨勢，rxt>0表示NPP呈上升趨勢，rxt<0表示NPP呈下降趨勢，若rxt通過0.05的顯著性水平檢驗(p<0.05)，則表示NPP呈顯著增加或顯著減少，若rxt通過0.01的顯著性水平檢驗(p<0.01)，則表示NPP呈增長極顯著或降低極顯著。未通過檢驗的為不顯著[1]，基于此6個等級進(jìn)一步分析NPP的變化情況(附圖3B)，研究區(qū)內(nèi)12.17%的區(qū)域增長顯著或增長極顯著，主要分布在多倫縣、正鑲白旗北部、正藍(lán)旗東部以及克什克騰旗東北區(qū)域，表明這些區(qū)域植被生長環(huán)境得到改善。研究區(qū)內(nèi)2.01%的區(qū)域降低顯著或極顯著，主要零星分布在克什克騰旗的中部以及蘇尼特左旗的東部區(qū)域。渾善達(dá)克沙地56.07%的面積屬于增長不顯著，19.79%的面積屬于下降不顯著，主要分布在沙地的中西部地區(qū)。

3 討論與結(jié)論

3.1 討 論

本文研究結(jié)果中渾善達(dá)克沙地植被整體呈現(xiàn)上升趨勢，該研究結(jié)果與國內(nèi)許多學(xué)者相一致，高曉霞等[25]基于logistic回歸模型分析渾善達(dá)克沙地的動態(tài)變化趨勢，整體呈現(xiàn)上升趨勢，孫建蕓[26]基于遙感數(shù)據(jù)分析渾善達(dá)克沙地的荒漠化情況，表明近年來渾善達(dá)克沙地植被覆蓋度明顯提高，荒漠化狀況得到改善，元志輝[12]分析2000—2013年渾善達(dá)克沙地植被凈初級生產(chǎn)力也表明渾善達(dá)克沙地NPP呈上升趨勢。從NPP的年內(nèi)變化來看，夏季NPP值占全年比例最大，且7月、8月份積累的NPP值最高，表明NPP的變化主要與植被的生長旺季有關(guān)。研究區(qū)植被東部NPP值大于中西部，原因可能是沙地呈東西走向，具有一定的地帶性，東部降水多，植被類型豐富，另外最初2000年實施生態(tài)工程后，多倫縣、正藍(lán)旗和錫林浩特市的白音錫勒牧場是最早的防沙治沙試驗區(qū)[11]，因而植被NPP值較高，中西部地區(qū)降水稀少，且處于沙地腹部，生態(tài)工程治理難度大，NPP值較低。

此外，在CASA模型選取的參數(shù)因子中，光能利用率ε對NPP的估算結(jié)果影響較大[13]。然而準(zhǔn)確建立水分脅迫因子對ε的影響函數(shù)是艱難的，傳統(tǒng)的Wε計算都是基于蒸散量和土壤系數(shù)的函數(shù)，其局限性在于空間降水和土壤質(zhì)地數(shù)據(jù)的空間異質(zhì)性表達(dá)[16]。本文考慮到利用MOD09A1數(shù)據(jù)提取近紅外波段和短紅外波段計算地表濕潤指數(shù)可以減少較為復(fù)雜且涉及到較多氣象數(shù)據(jù)資料的過程[27]。此外，基于MOD09A1計算的地表濕潤指數(shù)能夠很好地表現(xiàn)空間異質(zhì)性,比較準(zhǔn)確地反映地表的濕潤狀況,有效減少了蒸散量和土壤系數(shù)估算的人為誤差[13]。

最大光能利用率εmax的取值對最終的NPP估算結(jié)果也有較大影響[17]，目前，由于研究區(qū)的植被類型不同和氣候環(huán)境有所差異，以往研究在區(qū)域尺度和全球尺度上求取的εmax值各不相同，主要波動在0.389～1.044，如王鶯等[16]將甘南地區(qū)草地的εmax取值為0.608 g C/MJ；衛(wèi)亞星等[28]計算青海省εmax為0.403～0.908 g C/MJ；朱文泉等[17]根據(jù)誤差最小的原則，利用中國690個觀測點NPP實測數(shù)據(jù)，模擬各植被類型的εmax，該研究結(jié)果得到許多學(xué)者的應(yīng)用[8,19]。本文對εmax的取值也參照這一研究成果，取值為0.542 g C/MJ。但是由于εmax的取值隨著植被類型的不同而變化[15]，沙地實施生態(tài)治理工程后，植被類型和群落組成均發(fā)生了重大的變化[11]，在后續(xù)的研究中有必要實地考察采樣，獲取最新的沙地植被類型分布，結(jié)合地面實測NPP數(shù)據(jù)，計算研究區(qū)主要植被建群種類型的最大光能利用率，提高NPP模擬的精度。

另外，本文方法中使用的氣象數(shù)據(jù)來自氣象站的離散點數(shù)據(jù)，該方法的優(yōu)點在于其簡單且易于實施，缺點是隨著距觀測點的距離增加精度降低。而且插值方法本身可能會影響結(jié)果。因此，在后續(xù)的研究中可以考慮結(jié)合遙感數(shù)據(jù)反演光合有效輻射值[29]，提高研究結(jié)果的精確性。

3.2 結(jié) 論

(1) 渾善達(dá)克沙地2000—2016年NPP呈現(xiàn)波動上升趨勢，NPP多年均值為282.42 g/(m2·a)，年均增長量為2.136 2 g/(m2·a)，從2000年的276.82 g/(m2·a)增長到2016年的292.7 g/(m2·a)，表明近17 a渾善達(dá)克沙地植被固碳能力增加，植被長勢趨于良好。

(2) NPP季節(jié)變化明顯，春季NPP值呈下降趨勢，夏季和秋季呈上升趨勢，夏季NPP變化趨勢與全年NPP變化趨勢基本一致。從月變化來看，NPP值在7月達(dá)到最高值，7月、8月份積累的NPP值占總量的51.57%。

(3) 空間上，NPP整體呈由西向東遞增趨勢，沙地內(nèi)78.2%的面積呈增長趨勢，21.8%的面積呈現(xiàn)降低趨勢，2.17%的區(qū)域增長顯著或增長極顯著，2.01%的區(qū)域降低顯著或極顯著，表明國家在對渾善達(dá)克沙地進(jìn)行生態(tài)治理工程后取得了一定的效果，但是渾善達(dá)克沙地局部地區(qū)植被惡化狀況仍然存在，下一步有必要加強(qiáng)對沙地中西部地區(qū)生態(tài)工程建設(shè)的力度。