青海省草地生產力變化及其驅動因素

王亞暉，唐文家，李森*，趙鴻雁，謝家麗，馬超，顏長珍

（1.中國科學院西北生態環境資源研究院，甘肅蘭州 730000；2.中國科學院大學，北京 100049；3.青海省生態環境監測中心，青海 西寧 810007）

草地是地球上重要的生態系統，約占全球陸地面積的30%［1］。有研究表明，世界草地生態系統每年產生9.06×1011美元的生態價值，占所有生態系統服務總價值的2.72%［2］。我國草地資源豐富，總面積約4×106km2，占土地面積的41.7%［3］，因此草地生態系統質量變化對我國生產、生活和生態等影響巨大。青海省位于青藏高原東北部，是我國重要的天然草地分布區，同時也是我國重要的放牧區［4］。自1980年以來，由于過度放牧、草地開墾以及采挖蟲草等因素，青海省草地發生不同程度退化［5］。

草地凈初級生產力（net primary productivity，NPP）作為表征草地生產能力的重要指標［6］，其時空變化和驅動機制已經成為草地相關研究的重要內容［7-8］。有研究［9］表明，草地生產力變化受到氣候變化和人類活動的影響。為定量區分氣候因素和人類活動對草地生產力的影響，相繼提出了殘差趨勢分析法（RESTRENT）［10-11］、分段殘差分析法［12］、人類NPP法［13-15］等方法，用以分離人類活動對草地的影響，但是結果對驅動因素的構成和空間差異表達并不直觀。為解決這一問題，Leroux等［16］和王巨［17］先后構建多層判斷的驅動因素分類框架，分別定量分析非洲薩赫勒地區和河西地區植被變化驅動因素的空間結構差異，結果表明該方法在降水主導植被生長的干旱半干旱區［18］具有良好的適用性。在植被與溫度相關性更強的高寒地區［19-20］，該方法能否定量區分植被變化的主導因素，需要深入研究。

為闡明近年來青海省草地生產力變化趨勢，分析其驅動機制的空間異質性，本研究基于2001-2017年MODIS NPP時間序列數據產品，分析青海省草地NPP時空變化趨勢，并用優化后的驅動因素分類決策樹，量化草地NPP驅動因素的空間差異。

1 材料與方法

1.1 研究區概況

青海省地處89°35′-103°04′E和31°30′-39°19′N（圖1），東鄰四川，南接西藏，西靠新疆，北連甘肅，是青藏高原重要組成部分，總面積為72.23×104km2；地勢呈西高東低、南北高中間低，平均海拔大于3000 m；氣候為典型的高原大陸性氣候，全省年均氣溫在-5.1～9.0℃之間，絕大部地區的年均降水量低于400 mm［21］。

圖1 研究區地理位置與草地分布Fig.1 The location and grassland distribution of the study area

青海省天然草地面積41.92×104km2，其中74.70%的草地存在不同程度的退化［22］。為扭轉草地退化趨勢，國家和地方政府先后實施了退耕還草、退牧還草、三江源生態保護和建設工程等生態工程［23］，建立了三江源國家公園、祁連山國家公園等11處生態保護地。近年來，隨著高原氣候的暖濕化，加之多項生態工程的實施和自然保護區的建立，草地質量呈現好轉趨勢［24］。

1.2 主要數據和預處理

1.2.1 草地時空數據 青海省草地數據來源于中國土地覆被數據集（ChinaCover）［25］，包括2000、2005、2010和2015年等4期30 m空間分辨率的數據，其草地類型有溫帶草甸、溫帶典型草原、溫帶荒漠草原、高寒草甸、高寒草原、高寒荒漠草原等。為減少土地覆被類型變化帶來的NPP突變影響，將4期數據均為草地的區域作為研究區，并重采樣為500 m×500 m分辨率。經處理后，研究區草地面積為38.80×104km2。

1.2.2 MODIS產品數據集 MOD17A 3HGF.006數據集由美國國家航空航天局（National Aeronautics and Space Administration，NASA）生產分發，目前提供有2001-2020年全球500 m×500 m分辨率的年尺度NPP數據。MODIS NPP產品的反演算法相對成熟［26］，其精度經過驗證，已用于區域性植被生產力的相關研究［27-28］。本研究使用的青海省MOD17A 3HGF.006 500 m NPP數據經由Google Earth Engine（https://earthengine.google.com）平臺導出下載，其時間尺度與氣象數據年限保持一致，為2001-2017年。

1.2.3 氣象數據 氣象數據集為青藏高原及周邊地區氣溫和降水格點數據（1998-2017年）［29］，來自國家青藏高原科學數據中心（https://data.tpdc.ac.cn/zh-hans/）。該數據集包括1998-2017年青藏高原逐年的平均氣溫（0.1℃）和降水量（0.1 mm）格網數據，空間分辨率為1 km。為與NPP數據集分辨率一致，將該數據集重采樣為500 m×500 m，用以表示氣溫和降水因素。

1.3 研究方法

1.3.1 趨勢分析 Theil-Sen斜率估計（又稱Sen斜率估計）和Mann-Kendall檢驗法（MK檢驗）是無須數據正態分布的非參數估計與檢驗方法，對測量誤差和小的離群點不敏感［30-31］。因此，Sen-MK趨勢分析在長時間序列數據趨勢分析中被廣泛使用［32］。本研究NPP時間序列的Sen斜率計算公式如下：

式中：NPPi和NPPj為NPP時間序列數據中第i年和第j年的NPP值；Medain（）為取中值函數；β為所有數據對斜率的中值，即時間序列的Sen斜率。當β>0時，時間序列呈增加趨勢；β<0時，時間序列呈減少趨勢。時間序列趨勢的顯著性按不考慮值重復（NPP數據在年際間很難出現重復）和時間序列長度n≥10（本研究n=17）時的MK方法檢驗［33］，其統計量S的計算公式為：

式中：x i和x j分別是時間序列的第i和第j項，sgn（θ）為符號函數，定義為：

當n≥10時，對S標準化后的檢驗統計量Z可按以下公式計算：

根據Sen-MK趨勢分析的顯著性水平|Z|>1.96（P<0.05）和|Z|>2.58（P<0.01），將NPP變化趨勢分為5類（表1）。

表1 草地NPP變化趨勢分級標準Table 1 Classification of gr assland NPP tr ends

1.3.2 相關性分析 相關性分析常用來衡量要素間關系的密切程度。Pearson相關系數是對兩個要素x與y相關程度的度量，定義為：

式中：Rxy為要素x與y的Pearson相關系數，x i和yi分別表示第i組樣本中x和y的值，和分別為要素x與y的平均值，n為樣本量。

在多元相關分析中，偏相關系數可以控制其他變量的影響，突出所研究的兩個變量間的真實關系。在有3個變量時，控制第3變量后其他兩個變量間的偏相關系數為：

式中：r12·3為變量1和變量2的偏相關系數，r12、r13、r23分別為3個變量中兩兩間的Pearson相關系數。

基于本研究數據的自由度15（n-2），在95%置信度水平上，即|Rxy|>0.482時，兩個要素之間顯著相關。而偏相關的顯著性檢驗采用標準的t檢驗（雙邊），置信度水平同樣為95%。

1.3.3 殘差趨勢分析 在未經過相關性檢驗確定NPP與年均氣溫、年降水量同時相關之前，使用一元回歸模型擬合NPP與單氣候因素間的關系；而在確定NPP與兩個氣候因素都相關之后，使用多元回歸模型擬合其關系。具體的回歸模型如下：

式中：a和b為一元回歸模型的系數，c、d和e為多元回歸模型的系數，W為某一種氣候因素，P為年降水量，T為年均氣溫。用OLS方法計算兩種回歸模型的系數，并用標準的F檢驗衡量模型的顯著性，置信度水平取95%。

模型通過檢驗后，基于模型自變量和模型參數獲得NPP預測值，用以表示單一氣候因素或降水+氣溫對NPP的影響。而NPP殘差是NPP實際值與NPP預測值之間的差值，用來表示自變量以外的其他因素的影響。當NPP殘差序列在時間上有顯著的趨勢時，認為NPP受回歸模型中自變量以外的因素影響；否則認為NPP的變化僅受回歸模型內部自變量的影響。本研究用扣除氣候因素影響后的NPP殘差表示人類活動對NPP的影響部分。為提高非線性趨勢的檢測，NPP殘差時間序列的趨勢顯著性同樣使用Sen-MK方法檢驗。

1.3.4 基于決策樹的驅動因素分析 在自然條件下，除氣候因素外，地形、土壤等自然因素，在短期內難以發生較大變化，因此假定引起草地NPP變化的自然因素主要為氣溫和降水兩個氣候因素，而其他變化是由人類活動影響造成的。基于上述假定，綜合利用Sen-MK趨勢分析、相關性分析、殘差趨勢分析，構建影響草地NPP變化的主導因素的分類決策樹框架（圖2）。通過該框架，可結合草地NPP趨勢分析結果，將人類活動、氣溫和降水3個驅動因素組合為人類活動、氣溫、降水、人類活動+氣溫、人類活動+降水、氣溫+降水、人類活動+氣溫+降水等7種主導因素類型，并檢驗各區域草地NPP顯著變化是何種主導因素類型造成的。

圖2 草地NPP顯著變化的主導因素類型分類決策樹Fig.2 The decision tree to distinguish combination types of dominant factors driving significant grassland NPP changes

同時，基于草地NPP變化趨勢及其主導因素類型空間分類結果，可以從決策樹分類過程中獲取人類活動、氣溫和降水3個因素分別對草地NPP變化的作用類型：在氣候因素主導或參與主導草地NPP顯著變化的區域，疊加相關性分析結果和草地NPP變化趨勢，獲得氣候因素對草地NPP變化的4種作用類型，包括正相關促進、負相關促進、正相關抑制、負相關抑制。其中，正相關促進表示氣候因素的升高（增加）導致草地NPP增加；負相關促進表示氣候因素的降低（減少）導致草地NPP增加；正相關抑制表示氣候因素的升高（增加）導致草地NPP減少；負相關抑制表示氣候因素的降低（減少）導致草地NPP減少。在人類活動主導或參與主導草地NPP顯著變化的區域，根據Sen-MK趨勢分析和殘差趨勢分析結果，得到人類活動對草地NPP變化的2種作用類型，包括促進和抑制。其中，促進表示無氣候因素參與主導時草地NPP的Sen-MK趨勢為增加或扣除氣候因素后的NPP殘差趨勢為增加；抑制則反之。

以上研究方法及其結果中各指標作用區域面積和面積占比等的統計分析，基于Python 3.6，使用gdal和numpy等庫實現。

2 結果與分析

2.1 2001-2017年青海省草地NPP時空變化

2001-2017年，青海省草地NPP呈自東向西逐漸減低的趨勢，具有明顯的地域分異性（圖3），除柴達木盆地外，其空間分布與年均氣溫相對一致。其中，NPP高值區主要分布在青海湖周邊和龍羊峽上游黃河干流以東地區；NPP低值區主要分布在可可西里、昆侖山以及柴達木盆地邊緣。研究時段內，全省平均草地NPP緩慢波動上升（圖4），總體上與平均氣溫波動趨勢近于一致，比平均降水滯后一個階段。同時期內，青海省草地NPP與氣溫的關系更密切。

圖3 2001-2017年青海省草地多年平均NPP空間分布Fig.3 Spatial distribution of multi-year average grassland NPP in Qinghai Province during 2001-2017

圖4 2001-2017年青海省平均草地NPP、平均氣溫和平均降水的年際變化Fig.4 Inter-annual var iations of aver age grassland NPP，aver age temperatur e and average pr ecipitation in Qinghai Pr ovince during 2001-2017

2001-2017年，青海省70.59%的草地區域NPP變化不顯著，主要分布在玉樹州東部以及果洛州南部（圖5）。NPP呈極顯著增加的草地面積為4.61×104km2，占全省草地面積的11.88%，在空間上分布較為集中，主要在可可西里三江源區、扎陵湖和鄂陵湖的北岸及東岸、龍羊峽上下游河段、青海湖北部匯水區以及祁連山山區；這些區域多位于國家公園或者生態保護工程區內，草地生態得到較好的保護和恢復。NPP呈顯著增加的草地面積為6.69×104km2，占全省草地面積的17.25%，其空間分布與極顯著增加的區域一致，但主要在極顯著增加草地的周邊，而且更為破碎。NPP極顯著減少和顯著減少的草地面積分別為0.03×104和0.08×104km2，占全省草地面積的0.08%和0.20%；二者在空間上分布少而集中，主要在治多縣與曲麻萊縣交界的通天河兩側、黃河源的西北和東南部分區域、瑪多縣北部邊界、久治縣境內以及青海湖東南部日月山口附近，而且NPP極顯著減少的草地多位于NPP減少草地的內部，草地生產力退化區域的中心比邊緣嚴重。

圖5 2001-2017年青海省草地NPP變化趨勢空間分布Fig.5 Spatial distr ibution of NPP trends in Qinghai Pr ovince gr asslands during 2001-2017

2.2 草地NPP變化主導因素類型空間分布

青海省11.41×104km2草地NPP顯著變化區域中，不同因素類型主導的草地NPP顯著變化的面積和空間分布存在明顯差異（圖6）：草地NPP顯著變化的主導因素類型為氣溫的區域面積最大，其占比為60.66%，主要分布在唐古拉山、祖爾肯烏拉山和可可西里等高寒地區，以及鄂陵湖南岸和青海湖西北匯水區；其中高寒地區，人類活動稀少，溫度較低，氣溫是草地生長的限制性因子，而在鄂陵湖南岸和青海湖西北匯水區，由于水分充足，相對降水來說，草地生長對氣溫更加敏感。其次是人類活動主導的區域，面積占比為23.45%，主要分布在扎陵湖與鄂陵湖北岸、柴達木盆地西側、共和盆地南部、湟水谷地和祁連山山區；這些區域主要為人口相對密集，并且人類活動相對強烈的地區。第三是人類活動+氣溫主導的區域，面積占比為9.49%，主要分布在沱沱河以北到可可西里、哈拉湖東南部、黑河上中游等地區。降水、人類活動+降水和人類活動+氣溫+降水3種結構類型分別主導的區域面積較小，其占比均不足3.00%，主要在共和盆地、青海湖周邊沿岸和柴達木東側邊緣一帶。氣溫+降水主導草地NPP顯著變化的面積最小，僅有3.50 km2。

圖6 2001-2017年青海省草地NPP顯著變化的不同主導因素類型空間分布Fig.6 Spatial distribution of the combination types of dominant factors driving significant grassland NPP changes in Qinghai Province during 2001-2017

2.3 不同驅動因素對草地NPP變化的影響

2001-2017年，氣溫、降水和人類活動3個因素對青海省草地NPP顯著變化的主導區域和作用類型有著明顯差異（圖7）。

圖7 2001-2017年青海省氣溫（a）、降水（b）和人類活動（c）分別對草地NPP變化影響的空間分布Fig.7 Spatial distribution of effect of temperature（a），precipitation（b）and human activities（c）on grassland NPP in Qinghai Province dur ing 2001-2017

氣溫主導或參與主導草地NPP顯著變化的面積為8.11×104km2。其中，氣溫升高導致草地NPP增加的面積占99.81%，主要分布在高寒地區；氣溫對青海省草地NPP的正相關促進作用十分顯著。溫度降低導致草地NPP增加的面積占0.16%，主要分布在柴達木盆地的東南邊緣；這些地區氣候干旱、人類活動相對頻繁，相對于溫度，降水和人類活動對草地生長的限制性更強。而氣溫變化導致草地NPP減少的區域極少，僅有0.03%。

降水主導或參與主導草地NPP顯著變化的面積相對較少，為0.73×104km2。其中，由降水增加導致青海省草地NPP增加的面積占95.96%，主要分布在青海湖北岸、共和盆地至柴達木盆地東部一線以及龍羊峽上下游黃河兩側等海拔相對較低的地區。由降水增加導致青海省草地NPP減少的面積占3.74%，主要分布在久治縣境內。由降水減少導致青海省草地NPP增加的面積占0.30%，大多零星分布在唐古拉山附近，并且伴有人類活動因素的主導。而降水對NPP變化的正相關抑制作用并沒有被檢測到。

人類活動主導或參與主導草地NPP顯著變化的草地面積為4.15×104km2。其中，人類活動對草地NPP變化起促進作用的面積占95.40%，在空間上占了人類活動和人類活動+氣溫兩種因素類型主導區域的絕大部分。人類活動對草地NPP變化起抑制作用的面積占4.60%，大于草地NPP減少趨勢區域的面積，而在空間上除了與草地NPP減少趨勢區域的重疊部分，主要分布在共和盆地南側與西側、青海湖北岸、祁連縣東部的黑河北側等地。

3 討論

本研究基于MODIS NPP數據產品，使用Sen-MK趨勢分析，發現2001—2017年青海省29.41%的草地NPP存在顯著變化趨勢，其空間分布與劉洋洋等［34］的研究結果相一致，但范圍更大，這是本研究時間范圍更長導致的；與劉旻霞等［35］對青海省未來植被NPP變化的預測相符。然而，2000年以前，青海省草地已經發生嚴重的退化，草地NPP的顯著增加并不表示草地退化的消失，而且70.59%的草地NPP變化不顯著表明青海省的草地質量沒有完全恢復，仍存在不同程度的草地退化。

將驅動因素分類決策樹框架應用到青海省草地NPP變化研究中，為確定氣溫優先作為第二層考慮因素的可行性，本研究同時以降水優先作為第二層考慮因素進行了計算。結果發現兩種方式獲得的草地NPP主導因素類型圖中有95.21%的相同區域，具有高度的相似性：其中主導因素類型為氣溫的區域面積分別以60.66%和57.11%的占比均位列第一；其次都是人類活動單一主導的區域，且其面積大小與空間分布完全一致，表明人類活動的主導作用突出時，主導因素類型為人類活動的區域結果對氣溫和降水的順序并不敏感；而人類活動+氣溫主導的區域面積占比僅相差0.32%；人類活動和氣溫在兩種方式的結果中，均表現為引起草地NPP顯著變化的主要因素。而兩種方式結果有差異的區域，主要分布在青海湖周邊和共和盆地部分地區。使用Pearson相關分析和偏相關分析方法分析有差異地區的草地NPP與氣溫或降水間的關系，結果表明65.00%區域的NPP與氣溫的相關系數更高；但如果將降水作為決策樹框架的優先考慮因素，其中有42.69%區域的氣溫主導作用被完全忽略。另外，已有學者研究表明氣溫對青藏高原植被生長的驅動作用強于降水［36］。因此，本研究首先分離氣溫因素的結果是可靠的，而且優于降水優先的結果。而偏相關分析與多元回歸分析的引入，則可以充分利用上一層的判斷結果，使各層的分析方法更有針對性。同時，對比金凱等［37］和李傳華等［38］使用貢獻率區分氣候因素和人類活動的方法，兩者有所不同：用貢獻率可以量化同一空間上3個驅動因素的貢獻并進行比較，但不易突出驅動植被變化的主導因素；而驅動因素分類決策樹框架建立在整體Sen-MK趨勢顯著之上，主要探究主導草地生產力發生顯著變化區域的驅動因素及其空間上的差異，并在空間上得到一幅較為可靠和更加直觀的主導因素類型圖，對理解人類活動、氣溫和降水在草地長期趨勢中的驅動力構成有很大幫助。而且，相對于金凱等［37］研究得到的青海省人類活動對植被恢復基本無影響的結果，本研究方法更能體現出氣候不可控的情況下，國家和地方政府一系列的生態保護與修復工程對草地恢復的人類活動作用，肯定了生態工程對草地恢復的正面意義。但是本研究未對草地變化不顯著區域和短期內的波動變化做深入研究，也未分解多種驅動因素共同主導草地變化區域的單因素分量。對此，可以基于本研究的結果，在不同類型的驅動因素結構中使用適應的貢獻率計算方法，以進一步細化各因素間的差異。

人類活動對草地NPP變化起抑制作用的面積為0.19×104km2。但由于氣溫和降水的影響，這些區域內草地NPP變化并非完全為減少趨勢，其中，有41.56%區域的主導因素類型為單一的人類活動，并貢獻了草地NPP減少趨勢區域的72.36%，主要分布在居民點周邊，過度放牧等人類活動對草地的破壞作用掩蓋了氣候因素的影響；2.94%的區域與氣候因素起抑制作用的區域重合，共同貢獻了草地NPP減少趨勢區域的5.13%；55.42%的區域與氣候因素起促進作用的區域重合，在這些區域內人類活動對草地NPP變化的抑制作用小于氣候因素對草地NPP變化的促進作用，最終使草地NPP變化仍表現為增加趨勢。人類活動對草地NPP變化的促進作用，主要源于對人類破壞草地生態行為的限制和人類對草地生態的主動修復。其中，可可西里地區在建立自然保護區和國家公園后，已退化草地一方面在人類活動減少的影響下緩慢恢復，另一方面受氣溫升高的影響，草地生產能力得到提高，在主導因素類型分類結果中表現為人類活動和氣溫的共同主導。扎陵湖和鄂陵湖北部、共和盆地、龍羊峽上游、柴達木盆地東部邊緣等地區，是青海省沙漠化十分嚴重的地區［39-41］，自20世紀70年代以來，前人對這些沙漠化地區做出了眾多研究，并實施了一系列有效的沙漠化治理措施［42］，人類活動對草地的恢復作用得到突出體現；這些地區在本研究結果中主要表現為以單-人類活動為主導的草地NPP顯著增加。河湟谷地農業發達，人類活動密集，生態政策對人類活動的約束效果明顯，該地區退化草地在良好的水熱條件下快速恢復，表現為人類活動主導的草地NPP顯著增加。

4 結論

1）2001-2017年，青海省有11.41×104km2的草地NPP存在顯著變化趨勢，占全省草地面積的29.41%。其中NPP極顯著增加和顯著增加的草地面積占比分別為11.88%和17.25%，主要分布在可可西里、扎陵湖-鄂陵湖流域、祁連山山區、柴達木盆地東側以及青海湖周邊各流域；NPP極顯著減少和顯著減少的草地分布較少，面積占比分別為0.08%和0.20%，主要在曲麻萊縣和久治縣境內。17年來，全省草地生產力明顯提高。

2）青海省草地NPP顯著變化地區的7種主導因素類型中，氣溫主導的區域最大，占比60.66%；其次是人類活動和人類活動+氣溫，分別占比23.45%和9.49%。氣溫和人類活動是推動青海省草地NPP顯著變化的最主要的兩個驅動因素，而降水的作用只在柴達木盆地東側、共和盆地、龍羊峽下游黃河兩側等局部區域較為明顯；并且三者對青海省草地NPP變化的作用均以促進為主。

3）人類活動對青海省草地NPP變化起促進作用的區域占草地NPP增加趨勢區域面積的35.06%。在三江源國家公園與祁連山國家公園等自然保護區和沙漠化治理生態工程區等地區的人類活動正向效果顯著，人類活動對生態保護的作用不容忽視。但同時，草地NPP減少趨勢區域的77.49%與人類活動密切相關，人類活動仍是青海省草地退化的主要因素。

本研究能夠直觀可靠地表現出草地發生顯著變化區域的驅動因素空間異質性，對植被變化因素分析研究、生態工程效果評價、生態環境精準保護等具有重要意義。但在植被變化趨勢不顯著和驅動因素構成復雜的地區有待進一步的因素分解研究，以便更好地理解各因素間以及與植被變化間的作用機制。