基于CASA模型的渾善達克沙地植被NPP模擬及時空分析

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

(1.內蒙古師范大學 地理科學學院, 呼和浩特 010022; 2.內蒙古自治區遙感與地理信息系統重點實驗室,呼和浩特 010022; 3.內蒙古自治區生態與農業氣象中心, 呼和浩特 010051)

凈初級生產力(NPP)指植物在單位面積單位時間內所積累的有機物的數量，是由光合作用所產生的有機質總量中扣除自養呼吸后的剩余部分，是維持草原生態系統的物質基礎[1]。沙地植被NPP不僅能反映沙地生境的生產能力和質量狀況，也是沙地生態系統結構和功能的重要標志[2]。沙地植被NPP的大小，對沙地碳循環的合理調控和預測生態系統未來變化趨勢等具有重要的理論和現實意義。

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

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

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

1 材料與研究方法

1.1 研究區概況

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

圖1 渾善達克沙地地理位置

1.2 數據來源與預處理

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

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

1.3 研究方法

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

NPP=APAR×ε

(1)

APAR=SOL×FPAR×0.5

(2)

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

(3)

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

SOL由內蒙古太陽輻射數據獲取，FPAR根據MOD13Q1數據產品的NDVI反演獲取，計算公式參考文獻[12]，溫度脅迫因子Tε1和Tε2的計算方法詳見文獻[14]。

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

(4)

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

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

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

(5)

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

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

(6)

2 結果與分析

2.1 NPP模擬結果

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

圖2 渾善達克沙地2016年8月NPP估算結果

2.2 模型估算結果的驗證

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

表1 渾善達克沙地NPP模擬值與其他模擬結果比較

注：“—”表示參考文獻中無月值實測數據。

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

2.3 NPP年際變化特征

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

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

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

2.4 NPP季節與月變化特征

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

從月變化來看，渾善達克沙地17 a內月平均NPP值變化明顯(圖4D)。由于不同月份的水熱組合不同，各月的NPP值也各不相同，4—7月NPP值呈上升趨勢，到7月達到峰值，值為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季節與月變化特征

2.5 NPP的空間分布特征

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

2.6 NPP的空間變化趨勢

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

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

3 討論與結論

3.1 討 論

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

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

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

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

3.2 結 論

(1) 渾善達克沙地2000—2016年NPP呈現波動上升趨勢，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渾善達克沙地植被固碳能力增加，植被長勢趨于良好。

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

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