侯麗麗, 銀 山,2, 都瓦拉, 王艷琦, 彭秀清, 雅 茹
(1.內(nèi)蒙古師范大學(xué) 地理科學(xué)學(xué)院, 呼和浩特 010022; 2.內(nèi)蒙古自治區(qū)遙感與地理信息系統(tǒng)重點(diǎn)實(shí)驗(yàn)室,呼和浩特 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)實(shí)意義。
傳統(tǒng)的NPP實(shí)地估測方法無法滿足區(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)工程實(shí)施后的大背景下,我們對沙地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ù)工作具有重要意義。
渾善達(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á)克沙地地理位置
氣象數(shù)據(jù)來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http:∥cdc.cma.gov.cn)提供的2000—2016年沙地境內(nèi)10個旗縣氣象站點(diǎn)的月平均溫度數(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的變化趨勢,運(yùn)用最大值合成法(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 CASA模型估算NPP概述 CASA模型通過植物所吸收的光合有效輻射(APAR)和實(shí)際光能利用率(ε)來估算區(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í)際蒸散量和潛在蒸散量模型或者土壤濕度函數(shù)進(jìn)行估算,這些計算方法考慮的降水?dāng)?shù)據(jù)和土壤模型的經(jīng)驗(yàn)系數(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)
以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é)果
由于本文缺乏實(shí)測數(shù)據(jù),故采用間接驗(yàn)證的方式對估算結(jié)果進(jìn)行驗(yàn)證。將本文CASA模型估算的NPP值與其他人估算的NPP的均值或部分?jǐn)?shù)據(jù)進(jìn)行對比驗(yàn)證[12,15,18-22],來檢測估算數(shù)據(jù)的準(zhǔn)確性與可行性(表1)。
表1 渾善達(dá)克沙地NPP模擬值與其他模擬結(jié)果比較
注:“—”表示參考文獻(xiàn)中無月值實(shí)測數(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í)測數(shù)據(jù)相符合。此外,由于生態(tài)治理工程的推進(jìn),2010年后研究區(qū)內(nèi)植被NPP值不斷增大,本研究結(jié)果高于元志輝[12]估算的渾善達(dá)克沙地NPP值,該估算結(jié)果在一個合理的區(qū)間內(nèi),因此可以用于后續(xù)的結(jié)果與分析。
利用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é)果趨于一致。
從季節(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é)與月變化特征
將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%,集中分布在沙地的克什克騰旗和多倫縣。
為進(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)行顯著性檢驗(yàn)來分析渾善達(dá)克沙地植被NPP的變化趨勢,rxt>0表示NPP呈上升趨勢,rxt<0表示NPP呈下降趨勢,若rxt通過0.05的顯著性水平檢驗(yàn)(p<0.05),則表示NPP呈顯著增加或顯著減少,若rxt通過0.01的顯著性水平檢驗(yàn)(p<0.01),則表示NPP呈增長極顯著或降低極顯著。未通過檢驗(yàn)的為不顯著[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ū)。
本文研究結(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年實(shí)施生態(tài)工程后,多倫縣、正藍(lán)旗和錫林浩特市的白音錫勒牧場是最早的防沙治沙試驗(yà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個觀測點(diǎn)NPP實(shí)測數(shù)據(jù),模擬各植被類型的εmax,該研究結(jié)果得到許多學(xué)者的應(yīng)用[8,19]。本文對εmax的取值也參照這一研究成果,取值為0.542 g C/MJ。但是由于εmax的取值隨著植被類型的不同而變化[15],沙地實(shí)施生態(tài)治理工程后,植被類型和群落組成均發(fā)生了重大的變化[11],在后續(xù)的研究中有必要實(shí)地考察采樣,獲取最新的沙地植被類型分布,結(jié)合地面實(shí)測NPP數(shù)據(jù),計算研究區(qū)主要植被建群種類型的最大光能利用率,提高NPP模擬的精度。
另外,本文方法中使用的氣象數(shù)據(jù)來自氣象站的離散點(diǎn)數(shù)據(jù),該方法的優(yōu)點(diǎn)在于其簡單且易于實(shí)施,缺點(diǎn)是隨著距觀測點(diǎn)的距離增加精度降低。而且插值方法本身可能會影響結(jié)果。因此,在后續(xù)的研究中可以考慮結(jié)合遙感數(shù)據(jù)反演光合有效輻射值[29],提高研究結(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è)的力度。