青藏高原退化草地的恢復潛勢研究

王瑞涇，馮琦勝，金哲人，劉潔，趙玉婷，葛靜，梁天剛

（蘭州大學草地農業生態系統國家重點實驗室，蘭州大學農業農村部草牧業創新重點實驗室，蘭州大學草地農業教育部工程研究中心，蘭州大學草地農業科技學院，甘肅 蘭州 730020）

草地生態系統是全球范圍內分布最為廣泛的生態系統之一，約占地球陸地面積的25%［1］。同時，草地也是一種十分重要的自然資源，在氣候調節、防風固沙、生物多樣性保育、水土保持和維護生態平衡等方面具有不可替代的作用。如今世界各地的草地都面臨著退化，尤其是以高寒草地（約占比49%）為主要草地類型的青藏高原［2］，雖然近年來青藏高原高寒草地生態系統整體上呈現改善的狀態，但仍有部分草地存在不同程度的退化［3］。草地生態系統恢復的關鍵是草地植被的恢復［4］，目前青藏高原草地恢復措施主要有施肥、補播、圍欄封育、人工草地建植等。近年來，近自然恢復（close-to-nature restoration）理念受到了廣泛關注并在植被恢復中得到了有效利用，近自然恢復并不是完全摒棄傳統人工恢復，而是借助傳統人工恢復措施，主要依靠生態系統的自我調節進行可持續生態恢復，這也是一種較為切合高寒草地恢復的措施［5］。但青藏高原高寒草地近自然恢復潛勢的大小及其空間分布格局尚不明確。因此，開展青藏高原草地恢復潛勢研究對青藏高原草地生態系統保護與恢復具有重要意義。

近年來，關于恢復潛勢的研究已有一些進展，Venter 等［6］用Sen 斜率和Mann-Kendall 檢驗計算1986?2019年的增強植被指數（enhanced vegetation index，EVI）變化趨勢，分析南非的植被退化與恢復潛力；高海東等［7］與趙廣舉等［8］通過地形、土壤、植被、氣候、坡度等因素疊加分區，基于“生境越相似的區域，植被恢復潛力越接近”的相似生境法原則，以某區內植被覆蓋度的最大值評估黃土高原的恢復潛力；李海東等［9］采用層次分析法和模糊綜合評價方法，建立了高寒河谷沙地植被恢復潛力綜合評價模型來研究不同類型沙地植被的恢復潛力。潘竟虎等［10］基于潛在歸一化植被指數（normalized difference vegetation index，NDVI）數據，采用改進的CASA 模型模擬得到中國潛在植被凈初級生產力。

植被凈初級生產力（net primary productivity，NPP）可以較為準確地反映草地植被的生長狀況，且對氣候變化與人類活動的影響較為敏感［11］，適合作為評價草地恢復潛勢的指標。本研究擬在估計青藏高原草地現實與潛在NPP 的基礎上，將恢復潛勢定義為草地恢復現狀與恢復頂級之間的“距離”，依此分析青藏高原高寒草地恢復現狀和恢復潛勢。

基于此，本研究擬通過CASA 模型估計現實NPP（actual net primary productivity，NPP），通過Thornthwaite Memorial 模型估計潛在NPP（potential net primary productivity，PNPP），采用最大值合成法將19年來的年PNPP進行合成得到PNPP 最大值（maximum potential net primary productivity，PNPPm），用PNPPm與NPP 之間的差值評估草地的恢復潛勢，定量研究2001?2019年青藏高原草地恢復潛勢及其演變趨勢。本研究結果可以判定不同草地恢復現狀、恢復潛力與恢復價值，為青藏高原退化草地恢復政策的制定提供科學與理論支撐，對指導青藏高原草地保護與生態系統恢復具有重要意義。

1 材料與方法

1.1 研究區概況

青藏高原是我國的重點牧區之一，約占我國草地總面積的1/3。介于26°00′12″?39°46′50″N，73°18′52″?104°46′59″E，涉及青海、西藏、新疆、四川、甘肅和云南6 個省區共201 個縣級行政區，面積約為2.5724×106km2，約占中國陸地總面積的26.8%。青藏高原草地資源豐富，草地類型眾多，其中包括高寒草原、高寒草甸、高寒草甸草原、高寒荒漠等17 種草地類型［12］（圖1）。

圖1 青藏高原草地類型分布Fig.1 Distribution of grassland types on the Qinghai-Tibet Plateau

1.2 數據來源與預處理

1.2.1 氣象數據 2001?2017年的逐月降水和溫度插值數據采用的是Ge 等［13］研究中使用的數據。2018?2019年的逐月降水和溫度插值數據來自國家科技基礎條件平臺?國家地球系統科學數據中心（http：//www.geodata.cn）。

1.2.2 太陽總輻射量數據 在ArcGIS 軟件中，通過Area Solar Radiation 工具計算得到2001?2019年青藏高原每月太陽總輻射量（solar radiation，SOL）。該過程使用到的數字高程模型（digital elevation model，DEM）數據來源于NASA 的航天飛機雷達地形測繪任務數據（shuttle radar topography mission，SRTM）（http：//srtm. csi.cgiar.org/），空間分辨率為90 m［14］。

1.2.3 MODIS 數據 本研究使用了來自美國航空航天局對地觀測系統數據與信息系統（NASA’S earth observing system data and information system，EOSDIS）開發的MOD13A3 和MCD43A4 產品。MOD13A3 產品提供了2001?2019年逐月NDVI 數據，空間分辨率為1 km；MCD43A4 產品提供了2001?2019年逐日經過BRDF 校準的地表反射率數據，空間分辨率為500 m，它包括MODIS1～7 波段，本研究僅使用第2 和6 波段計算所得地表水體指數（land surface water index，LSWI）作為CASA 模型的輸入變量計算NPP。

對于以上數據本研究使用MODIS 數據重投影工具（MODIS reprojection tool，MRT）進行了格式轉換、拼接與轉投影，處理后的數據為Geo-Tiff 格式，采用WGS 1984 Albers 投影。并將所有數據重采樣至1 km，裁剪出青藏高原地區的影像數據。與此同時，本研究還采用最大值合成法（maximum value composite，MVC），合成了逐月地表水指數（land surface water index，LSWI）數據，用于后續研究［15］。

1.3 現實凈初級生產力（NPP）

本研究采用CASA 模型對2001?2019年青藏高原的NPP 進行估算，該模型是當前應用最廣泛的凈初級生產力估計模型［16］。該模型主要基于光合有效輻射（absorbed photosynthesis active radiation，APAR，MJ·m?2）和光能利用率（ε，g C·MJ?1）估計NPP，其公式為：

式中：SOL 為太陽總輻射量（MJ·m?2），是用數字高程模型（DEM）通過ArcGIS 軟件計算得出［14］；FPAR 為植被冠層對光合有效輻射的吸收比例，由各植被類型的NDVI 及其最大值、最小值，比值植被指數（SR）及其最大值、最小值以及FPAR 最大值、最小值確定，計算過程參照朱文泉等［16］的研究；ε 為光能轉化率，由月平均氣溫、植物生長最適月平均氣溫、LSWI 以及理想條件下最大光能轉化率（ε*，本研究選取全球通用值0.389 g C·MJ?1）得出，計算過程參考Xu 等［17］的研究。

1.4 潛在凈初級生產力（PNPP）

本研究采用經驗證的Thornthwaite Memorial 模型對2001?2019年青藏高原的PNPP 進行估算，該模型是在Miami 模型的基礎上對Thornthwaite 潛在蒸散發模型加以修正而建立的。是利用蒸散與碳吸收的關系，以降水和溫度為輸入，來估計PNPP，這一方法比Miami 模型更加精確［18］。其表達式為：

式中：PNPP 為年潛在凈初級生產力（g C·m?2）；v為年平均實際蒸散量（mm）；L為年平均蒸散量（mm）；t和r分別表示年平均溫度（℃）與年降水量（mm）。年平均溫度（t）和年降水量（r）使用2001?2019年研究區降水和溫度空間插值數據［13］，空間分辨率為1 km，數據為Geo-Tiff 格式，采用WGS 1984 Albers 投影。

1.5 恢復潛勢

本研究認為天然草地的恢復潛勢（R）可以定義為潛在凈初級生產力的最大值（PNPPm）與現實凈初級生產力之間的差值，PNPPm是基于2001?2019年PNPP 數據，采用最大值合成法得到的，即：

1.6 趨勢分析

將Theil-Sen Median 趨勢分析、Mann-Kendall 檢驗與Hurst 指數相結合，可以評估2001?2019年青藏高原NPP 變化趨勢分布特征。在進行Theil-Sen Median 趨勢分析時，由于幾乎不存在SNPP完全等于0 的區域，所以本研究根據實際情況，認為SNPP＜?2 的區域為退化區域；?2＜SNPP＜2 的區域為穩定不變的區域；SNPP＞2 的區域為改善區域。本研究將Mann-Kendall 檢驗在0.05 置信水平上的顯著性檢驗結果劃分為趨勢顯著變化（Z＞1.96 或Z＜?1.96）和不顯著變化（?1.96≤Z≤1.96）。將Hurst 指數小于0.5 的區域定義為不可持續性變化，Hurst 指數大于0.5 的區域為可持續性變化。具體計算方法參考袁麗華等［19］的研究。

2 結果與分析

2.1 NPP 的時空變化特征

利用2001?2019年的年NPP 數據計算得到多年平均NPP 空間分布圖（圖2）。從圖2 可以看出，青藏高原NPP 表現為東南部的甘孜藏族自治州、阿壩藏族羌族自治州以及甘南藏族自治州等地明顯較高，西北部的阿里、那曲以及海西西部等地區偏低，由此可見青藏高原草地NPP 呈現出由西北向東南逐漸增加的趨勢，即青藏高原東南部草地狀況較西北部更好。

圖2 2001-2019年青藏高原年平均NPP 空間分布Fig.2 Spatial distribution of annual mean NPP over the Qinghai-Tibet Plateau from 2001 to 2019

基于2001?2019年NPP 數據，用Theil-Sen Median 趨勢分析、Mann-Kendall 檢驗與Hurst 指數進行分析，得到了2001?2019年NPP 的變化趨勢與可持續性耦合信息（圖3 和表1）。

圖3 2001-2019年青藏高原NPP 變化特征空間分布Fig.3 Spatial distribution of NPP variation characteristics over the Qinghai-Tibet Plateau from 2001 to 2019

表1 NPP 變化趨勢統計Table 1 Statistics of NPP trend（%）

從圖3 和表1 中可以看出，耦合共分為6 種情況，即未來趨勢不確定、持續性輕微退化、持續性顯著退化、持續性穩定不變、持續性輕微改善以及持續性顯著改善。除去無法預測未來變化趨勢的區域，青藏高原草地NPP 狀況改善的區域明顯多于退化區域。其中，無法預測區域占36.21%，持續恢復區域占40.98%，持續穩定不變區域占12.72%，而持續退化區域僅占3.47%。相對而言，退化較為嚴重的草地類型為熱性灌草叢和溫性草甸草原。而青藏高原中面積較大的草地類型?高寒荒漠、高寒草甸草原、高寒草甸、溫性草原化荒漠以及溫性荒漠均呈現出較為明顯的恢復趨勢。由此可見，整個青藏高原天然草地還是以可持續的恢復狀態為主。

為了分析不同草地類型的NPP年際變化，以便更好地了解青藏高原NPP 的變化情況，本研究制作了各草地類型NPP年際變化圖（圖4）。從圖4 可以看出，在2001?2019年間大部分草地類型年NPP 都呈現波動上升趨勢，且在2018年有大幅增加，多類草地增幅達到100 g C·m?2，這可能與當地氣候變化或政策制定有關。

圖4 2001-2019年青藏高原各草地類型NPP年際變化Fig.4 Interannual variation of NPP of different grassland types in the Qinghai-Tibet Plateau from 2001 to 2019A：熱性灌草叢類Thermal shrub tussock；B：熱性草叢類Thermal tussock；C：暖性灌草叢類Warm-temperate shrub tussock；D：暖性草叢類Warmtemperate tussock；E：低地草甸類Lowland meadow；F：山地草甸類Mountain meadow；G：溫性草甸草原類Temperate meadow steppe；H：溫性草原類Temperate steppe；I：溫性荒漠草原類Temperate desert grassland；J：溫性草原化荒漠類Temperate steppe desert；K：溫性荒漠類Temperate desert；L：高寒草甸類Alpine meadow；M：高寒草甸草原類Alpine meadow steppe；N：高寒草原類Alpine grassland；O：高寒荒漠草原類Alpine desert grassland；P：高寒荒漠類Alpine desert；Q：沼澤類Marsh.下同The same below.

2.2 PNPPm的空間分布特征

為了直觀地反映青藏高原2001?2019年PNPP 的空間分布以及比較各草地類型的具體情況，基于2001?2019年PNPP 數據，用最大值合成法制作了PNPPm的空間分布圖（圖5），從圖中可以看出，PNPPm在青藏高原的分布呈現出明顯的南高北低的態勢，其中日喀則地區、阿壩藏族羌族自治州、甘南藏族自治州和甘孜藏族自治州等地是青藏高原PNPPm較高的幾個地區，這與NPP 的空間分布相近。統計不同草地類型的PNPPm（圖6），可以看出熱性灌草叢、熱性草叢、暖性灌草叢和暖性草叢是青藏高原中PNPPm較大的4 種草地類型，其PNPPm約為1300～1500 g C·m?2，而低地草甸、溫性荒漠草原、溫性草原化荒漠、溫性荒漠、高寒荒漠草原以及高寒荒漠PNPPm較小，約為250～750 g C·m?2，其他草地類型PNPPm在1000 g C·m?2左右。這表明在僅考慮氣象因素的情況下，青藏高原東南與西南部草地會具有更好的生長狀況，熱性灌草叢、熱性草叢、暖性灌草叢和暖性草叢應為長勢較好的4 種草地類型。

圖5 2001-2019年青藏高原PNPPm空間分布Fig.5 Spatial distribution of PNPPm over the Qinghai-Tibet Plateau from 2001 to 2019

圖6 2001-2019年青藏高原各草地類型PNPPm統計Fig.6 PNPPm statistics of grassland types in the Qinghai-Tibet Plateau from 2001 to 2019

2.3 恢復潛勢的空間分布特征

基于2001?2019年的PNPPm與NPP 數據得到恢復潛勢分布圖，并分析了19年平均恢復潛勢的空間分布（圖7）。從圖7 可以看出，恢復潛勢的分布與PNPPm的分布情況相似，東南與西南兩部分較高，北部普遍偏低。其中日喀則地區與甘南地區恢復潛勢最大，而阿里北部以及那曲北部等地區恢復潛勢較低。因此在日喀則地區與甘南地區開展草地恢復與保護的價值最大，前景最優。

圖7 2001-2019年青藏高原恢復潛勢空間分布Fig.7 Spatial distribution of recovery potential over the Qinghai-Tibet Plateau from 2001 to 2019

基于2001?2019年青藏高原恢復潛勢數據，分草地類型制作了19年平均恢復潛勢的柱狀圖（圖8）。從圖中可以看出，熱性灌草叢、熱性草叢、暖性灌草叢、暖性草叢、山地草甸以及沼澤是青藏高原恢復潛勢較大的6 種草地類型，而低地草甸、溫性草原化荒漠、溫性荒漠、高寒荒漠草原、高寒荒漠、溫性荒漠草原則恢復潛勢較小，這一結果與PNPPm的情況相近，在開展青藏高原草地恢復時應優先考慮熱性灌草叢、熱性草叢、暖性灌草叢、暖性草叢、山地草甸以及沼澤這6 種草地類型。

圖8 2001-2019年青藏高原各草地類型恢復潛勢統計Fig.8 Recovery potential statistics of grassland types in Qinghai-Tibet Plateau from 2001 to 2019

2.4 氣溫、降水與恢復潛勢的關系

本研究基于2001?2019年的年均降水與年均溫數據得到青藏高原降水與氣溫空間分布圖，并分析了它們與恢復潛勢的關系（圖9 和圖10）。從圖中可以看出，青藏高原氣溫較高的區域為東部、南部以及西南部，降水則具有明顯的由東南向西北逐漸降低的態勢，這與青藏高原西南與東南兩部分恢復潛勢較高的情況不謀而合，這說明青藏高原恢復潛勢或許受到降水和溫度的影響，且溫度的影響可能較降水更大一些，這與德吉央宗等［20］研究得出的結論一致。

圖9 2001-2019年青藏高原氣溫空間分布Fig.9 Spatial distribution of temperature over the Qinghai-Tibet Plateau from 2001 to 2019

圖10 2001-2019年青藏高原降水空間分布Fig.10 Spatial distribution of precipitation over the Qinghai-Tibet Plateau from 2001 to 2019

2.5 恢復潛勢的變化趨勢

基于2001?2019年的青藏高原恢復潛勢數據得到青藏高原恢復潛勢變化趨勢圖，并分析了青藏高原恢復潛勢的變化趨勢以及成因（圖11）。從圖中可以看出，青藏高原草地恢復潛勢提升的區域主要集中在西南和東北部，尤其是西南部，降低的區域主要為青藏高原中部，而青藏高原恢復潛勢整體為提升、降低與不變三者均有的局面。青藏高原恢復潛勢變化趨勢的這種空間分布格局也可以作為研究青藏高原草地恢復價值的一個標準，即西南部與東北部草地恢復價值更大。

圖11 2001-2019年青藏高原恢復潛勢變化趨勢Fig.11 Change trend of recovery potential of Qinghai-Tibet Plateau from 2001 to 2019

3 討論

2001?2019年青藏高原NPP 呈現東南高、西北低的空間分布特征，且各草地類型NPP 均呈現波動上升的趨勢，這與陳舒婷等［21］的結論一致。本研究基于CASA 模型計算NPP 時，光能轉化率（ε）選用了全球通用值0.389 g C·MJ?1，然而，目前對于光能轉化率的看法不一，如王保林等［22］認為我國草原最大光能轉化率的取值區間為0.608～1.000 g C·MJ?1；Raymood 等［23］認為其最大值為3.5 g C·MJ?1。本研究的取值可能會造成對NPP 的低估，因此還需在后續研究中進一步探討更適合估算青藏高原NPP 的最大光能轉化率。

本研究使用Thornthwaite Memorial 模型估算青藏高原2001?2019年PNPPm為822.70 g C·m?2。潘竟虎等［10］利用改進的CASA 模型對全國植被PNPP 進行了模擬，其均值為468.94 g C·m?2；呂振剛等［24］采用改進的CASA 模型估算2001?2010年中國華北落葉松的PNPP，其均值為342.7 g C·m?2；任正超等［25］采用改進的綜合順序分類系統以及改進的CASA 模型模擬1982?2012年中國自然植被PNPP，其均值為278.7 g C·m?2。本試驗的研究結果較上述結果偏大，最主要的原因是本研究采用最大值合成法合成了19年年最大PNPP 得到了本研究所使用的PNPPm；最后本研究采用Thornthwaite Memorial 模型，而前人的研究多采用改進的CASA 模型，模型的選擇差異也是導致研究結果有差距的重要原因。因此，在之后的研究中還可以探討其他計算PNPP 的模型對于青藏高原草地恢復潛勢的計算精度，如改進的CASA 模型、綜合自然植被凈第一生產力模型、NPP 分類指數模型等。

本研究認為恢復潛勢應為某地有可能達到的最佳恢復程度，即可能達到的PNPP 與NPP 之間的差值，因此計算恢復潛勢時使用了2001?2019年最大值合成的PNPPm，而前人研究大多認為恢復潛勢為年PNPP 與年NPP之間的差值，這會導致本試驗的PNPPm較其他相關試驗偏大，恢復潛勢也會因此偏大，但這并不影響研究區恢復潛勢的空間分布以及得出的青藏高原各地恢復價值的大小以及恢復措施。

本研究發現青藏高原恢復潛勢呈現明顯的西南部最高、東南部次之而北部偏低的分布情況，分析認為形成這種分布格局的原因可能與青藏高原的氣溫、降水、海拔以及人類活動有關，在今后的研究中應著重探討青藏高原恢復潛勢分布的成因并采用多種方法估算青藏高原恢復潛勢，減少系統誤差。

4 結論

本研究采用PNPPm與NPP 之間的差值作為恢復潛勢，著重探討青藏高原天然草地NPP、PNPP 與恢復潛勢的空間分布與時空變化，以期為指導青藏高原草地恢復提供理論支持，主要結論如下：

1）2001?2019年青藏高原草地NPP 由西北向東南逐漸增加，青藏高原草地整體以可持續的恢復狀態為主，其中熱性草叢、低地草甸、溫性荒漠以及高寒荒漠草原等草地類型恢復態勢明顯。

2）2001?2016年間青藏高原草地NPP 波動上升，上升趨勢并不明顯，而在2016?2018年期間青藏高原各草地類型NPP 大幅增加。

3）2001?2019年青藏高原草地恢復潛勢呈現西南與東南部較高，北部偏低的分布情況。由此可見日喀則地區、阿里南部地區、阿壩藏族羌族自治州以及甘南藏族自治州等地草地具有較高的恢復價值，在這些地區開展草地恢復的前景更好。