雪寶頂及其周邊區域NPP時空格局演變及氣候響應

邱 豪，王建華,2， 董廷旭,2，張雪茂， 嚴 霜

(1.綿陽師范學院資源環境工程學院，四川綿陽 621006；2.生態安全與保護四川省重點實驗室，四川綿陽 621006)

0 引言

圖1 雪寶頂區域地勢圖Fig.1 The area of Xuebaoding

植被凈初級生產力(NPP)是指綠色植物在一定面積和時間內通過光合作用所產生的有機物總量減去呼吸作用消耗后的積累部分，是能夠反應區域氣候變化、區域碳積蓄與陸地生態系統響應的重要參數[1-3].NPP是物質與能量轉化研究的基礎，同時能夠作為評價陸地生態系統可持續發展的一個重要生態指標[4-5].近幾十年，NPP估算的模型陸續被提出，其中，光能利用率模型認為能夠利用任何影響植被生長的要素估算NPP[6].1993年，學者Potter、Felid依次建立起以遙感數據作為估算模型數據源的CASA模型.Potter依據植被NDVI指數以及氣象因子實現了基于光能利用率模型的全球陸地生態系統NPP的估算模擬[7].Field隨后對CASA模型中光能利用率參數的取值給出了具體的界定[8]，使得CASA模型估算結果的精度和可靠度得到提高.CASA模型所需數據主要來源于遙感衛星數據，覆蓋廣泛且時間分辨率較高，且所需參數相對較少、誤差小，CASA模型已應用于區域植被NPP長期變化的監測和NPP研究中[9-13].朱玉果等通過CASA模型分析寧夏省不同草地2001-2010年NPP的時空分布特征，并探索了不同草地NPP與氣象因子的相關性[14];王耀斌等利用CASA模型結合遙感數據，分析了2000-2015年秦巴山區植被凈初級生產力，探索了NPP的時空變化及其趨動因子[15].通過遙感手段研究區域植被凈初級生產力已經成為如今一種重要手段.

雪寶頂區域位于四川省西北部，包含平武縣、松潘縣、北川縣、茂縣和黑水縣，地處四川盆地與青藏高原的過渡帶，是川西北地區重要的生態屏障，在維系川西北地區生態環境安全中有著重要作用.但是由于雪寶頂區域地形差異大，站點較少導致無法完整掌握雪寶頂區域NPP時空分布.因此，尚未見雪寶頂區域NPP相關研究，探究雪寶頂區域NPP對氣候響應的報道也相對較少.本文利用2008-2017年的遙感影像和氣象站點資料，通過CASA模型對近10a雪寶頂區域的NPP進行估算，并在此基礎上分析雪寶頂區域NPP時空分布特征，探討與氣象因子的相關性，以期填補該區域NPP研究空白，為生態修復提供科學依據.

1 數據與方法

1.1 數據來源及預處理

氣象數據來源于國家氣象信息中心(http://www.nmic.cn/)選取2008-2017年平武縣、松潘縣、北川縣及其周邊共10個氣象站的日降水量(mm)、月平均氣溫(℃)、太陽總輻射(MJ/m2)、站點編號與經緯度坐標.在ARCGIS10.2中，通過空間插值的方法得到研究區的月降水總量(mm)、月平均氣溫(℃)以及太陽總輻射(MJ/m2)分辨率為250 m*250 m的柵格圖.

遙感數據選取來源于美國宇航局(https://ladsweb.modaps.eosdis.nasa.gov/)分辨率為250 m*250 m的全球植被16日合成產品MOD13Q1，對數據進行格式轉化，將影像重投影為WGS_1984，利用ENVI軟件對重投影后的影響進行研究區裁剪，并對所需波段進行合成處理.土地覆被數據來源清華大學的土地利用/土地覆被數據，并根據NPP估算需求將其分為針葉林、闊葉林等7大類.

1.2 NPP估算方法

本文采用了基于光能利用率的CASA模型對雪寶頂自然保護區進行凈初級生產力估算，CASA模型的凈初級生產力估算公式如下：

NPP(x,t)=APAR(x，t)×ε(x,t)×0.5

式中APAR(x，t)表示像元x在t月吸收的光合有效輻射， ε表示像元x在t月的實際光能利用率.

APAR(x，t)是根據植被對紅外及近紅外波段的反射特征實現的，植被吸收的光合有效輻射取決于太陽總輻射和植被本身特征，可由以下公式獲得[7]：

APAR(x,t)=SOL(x,t)×FPAR(x,t) ×0.5

式中SOL(x,t)表示t月在像元x處的太陽總輻射量，FPAR(x,t)為植被層對入射光合有效輻射的吸收比例，常數0.5表示植被所能利用的太陽有效輻射占太陽總輻射的比例.

ε光能利用率的估算公式如下[16]：

ε=Tε1×Tε2×Wε×ε*

Tε1表示在低溫和高溫時植物內在的生化作用對光合作用的限制而降低植被凈初級生產力，Tε2表示環境溫度從最適溫度向高溫或低溫變化時植物光能利用率逐漸變小的趨勢，Wε反映了植物所能利用的有效水分影響光能轉化率程度[17].

1.3 NPP時空分異模型

1.3.1 趨勢分析 一元線性回歸分析能夠模擬每個柵格的變化趨勢，通過顯示每一時間段像元的空間變化特征，綜合反映一定時間序列的區域格局演變規律[18].該方法通過擬合不同時間段數據，消減異常因子對NPP的影響，準確地反映區域NPP的變化趨勢.Slope即雪寶頂區域的NPP近10a時間序列每個柵格像元的線性變化斜率反應了在10a的時間序列中的變化趨勢，公式如下：

式中，Slope為像元NPP回歸方程的斜率；NPPi為第i年的NPP值；n為研究時段跨度(2008-2017共10a).

當Slope<0時，表明該時間序列中，區域NPP處于下降趨勢，且值越小下降趨勢越明顯；反之則為增加.趨勢顯著性通過相關系數進行檢驗，計算公式如下：

1.3.2 Hurst指數 基于重標極差(R/S)分析方法的Hurst指數最早是由英國水文學家赫斯特在研究尼羅河水庫流量和貯存能力的關系時提出，能夠有效定量描述時間序列信息自相似性和長相關性[19].

NPP時間序列NPPi，i=1，2，3，…，n，對于任意正整數m，定義該時間序列：

(1)差分序列

ΔNPPi=NPPi-NPPi-1

(2)均值序列

(3)累計離差

(4)極差

對于比值R(m)/S(m)≌R/S，若存在如下關系R/S∝mH，則說明分析的時間序列存在Hurst現象，H稱為Hurst指數，該指數值與時間序列持續性關系見表1.

1.3.3 相關性分析 相關系數能夠定量表達兩個變量之間相關程度，本文針對NPP值與降水、太陽輻射與氣溫等氣象因子的進行相關度分析：

圖2 雪寶頂區域年均NPP分布特征Fig.2 Annual average NPP distribution characteristics in the Xuebaoding region

式中，r為X和Y之間的簡單相關系數，Xi為第i月的NPP值，Yi月的相關因子值，分別為NP平均值和對應相關因子平均值，i為時間月變量.

2 結果與分析

2.1 NPP空間分布格局

2008-2017年雪寶頂區域的年平均植被凈初級生產力如圖2所示，研究區近10a的多年NPP均值為398.34 gCm-2a-1，東部山地丘陵地區(北川、平武)多年NPP均值主要集中在400～500 gCm-2a-1，其中雪寶頂自然保護區、千佛山等區域NPP均值超過600 gCm-2a-1；西北部高原草地地區多年NPP均值主要集中在300～400 gCm-2a-1；西南部多年NPP均值分布差異較大，最低不足100 gCm-2a-1，最高超過600 gCm-2a-1；由此，雪寶頂區域的NPP分布存在較強的空間異質性，這與雪寶頂區域的地理位置和氣象特征有關.東部的丘陵山地地區處于四川盆地與青藏高原的邊緣地區，年均溫度在11～13℃左右，降水充沛.且植物群落主要以闊葉林和針葉林為主，植被覆蓋度高，凈初級生產力強.西北部處于青藏高原東緣，地形多為丘狀高原，分布以草地和沼澤為主，植被覆蓋度高，凈初級生產力較強.研究區西南部多為高山河谷，海拔差異較大，低海拔地區以闊葉林為主，植被覆蓋度高，凈初級生產力強.研究區內的雪寶頂、三奧雪山因常年覆蓋積雪造成植被覆蓋度極低，凈初級生產力弱.

2.2 雪寶頂區域NPP時間變化特征

雪寶頂區域2008-2017年NPP值變化如圖3所示.近10a來，雪寶頂區域NPP處于波動下降趨勢，下降趨勢明顯，線性增長率為-6.862 8 gCm-2a-1.2012年的NPP年均值最低，為318.65 gCm-2a-1，2008年的NPP均值最高，為403.06 gCm-2a-1.將研究區NPP均值劃分為6個等級，統計其分布面積占比(圖3)，近10a雪寶頂區域NPP值200～300 gCm-2a-1和300～400 gCm-2a-1的面積整體呈現增加的趨勢，而400～500 gCm-2a-1和≥500 gCm-2a-1整體則不斷下降.NPP值300～400 gCm-2a-1和400～500 gCm-2a-1面積占總面積的60%以上，其中NPP值300～400 gCm-2a-1面積常年超過總面積的50%.研究區NPP值年內動態變化呈現典型單峰的特性(圖4)，三月份植被生長期開始的時候，NPP值急速增加，在六月份因為太陽輻射的原因，NPP值增加在一定程度上放緩或者出現下降，到7、8月份達到峰值，9月份以后隨著生長季的結束，NPP值迅速下降.

不同地表覆被類型對應的年均NPP值的差異較大，其中，常綠闊葉林的多年NPP均值最高，達到了555.74 gCm-2a-1，主要分布在研究區東部丘陵高山地區以及西南海拔低的地區；其次是農用地、灌叢、高山草甸和高原草甸，其它地表覆被類型多年NPP均值都在300 gCm-2a-1以下，其中城市的多年NPP值最低僅為102.2 gCm-2a-1.

圖3 2008-2017多年NPP均值年際變化以及年內變化Fig.3 The interannual and annual changes of NPP from 2008 to 2017

圖4 不同地表覆被NPP年際以及年內變化Fig.4 The interannual and annual changes of NPP

2.3 雪寶頂區域NPP的空間變化特征

圖5 雪寶頂區域NPP變化趨勢Fig.5 Variation trend of NPP in Xuebaoding region

2.3.1 NPP的空間變化特征 本文通過一元線性回歸逐像元分析近10a雪寶頂區域NPP變化趨勢，從圖5中可以看出，研究區95.7%的區域，其線性斜率小于0，僅有4.3%的區域，其線性斜率大于0.從空間變化分析，研究區自東向西NPP下降趨勢逐漸減弱，其中趨勢分布在-20～-10 gCm-2a-1的區域面積最廣，占研究區的43%，而下降趨勢最明顯的區域集中研究區東部平武縣和北川縣；NPP呈現增加趨勢的地區集中分布在黑河流域和雪寶頂，特別是黑河流域，增長率最高可達55 gCm-2a-1.

圖6 NPP變化HURST指數Fig.6 HURST index of NPP

2.3.2 雪寶頂區域NPP變化持續性分析 本研究逐像元計算了雪寶頂區域近10a的Hurst指數(圖6)以分析雪寶頂區域NPP變化趨勢的持續性.雪寶頂區域的Hurst指數在0.21～0.87≥之間，均值為0.49，Hurst指數小于0.5的區域占雪寶頂區域的56%，而大于Hurst指數大于0.5的面積占雪寶頂區域的44%，表明大部分雪寶頂區域NPP具有較弱的持續性，其變化特征反向特征要高于同向特征.Hurst指數高值主要分布在研究區西北部若爾蓋大草原邊緣和西南部黑河流域地區，其值在0.75左右；低值最要分布在雪寶頂以東平武縣境內.從不同的地表覆被類型分析得出，針葉林、闊葉林、山地草地和農牧地的Hurst均值小于0.5，其中，闊葉林的Hurst均值最低，僅為0.46；灌叢草甸、草甸草原和沼澤草甸Hurst均值大于0.5，沼澤草甸Hurst均值最大，達到了0.56，其次為灌叢.通過將研究區近10aNPP變化斜率與Hurst指數疊加分析獲得近10a雪寶頂區域NPP變化持續性特征(圖6).從圖6中可以得到，NPP變化趨勢減少轉增加、持續減少、持續增加和增加轉減少的面積分別占研究區總面積的54%、40%、4%和2%，即目前NPP呈現減少的大部分地區在未來會出現逆轉，未來表現出逆轉主要分布在雪寶頂以東的丘陵山地地區，地表覆被類型以闊葉林和針葉林為主；而仍有40%的區域NPP會持續減少，持續減少的區域主要分布在研究區西北地區的丘狀高原地帶，該區域地表覆被類型以草甸草原為主.

2.4 雪寶頂區域NPP與氣象因子相關性分析

2.4.1 年際NPP與氣象因子相關性分析 本研究逐像元分析了近10aNPP與同期氣象因子(降水、氣溫和太陽總輻射)的相關關系(圖7)，并進行顯著性檢驗.結果顯示，年際NPP變化與年際氣溫與降水波動關系不大，整體呈不明顯負相關，這可能是因為該區域近10a年均降水和溫度較為穩定，變化波動不大.其中，研究區NPP年際變化與降水波動呈正相關的區域集中分布在東部低海拔區域和黑河流域沿岸，該區域多為河谷氣候，降水變化大且不穩定，地表覆被類型多為闊葉林為主；而年際NPP變化與太陽總輻射因子存在著顯著的正相關，呈現正相關的面積超過了95%，因此，對太陽輻射量的依賴性較大.通過分析發現雪寶頂區域年際變化尺度上，太陽輻射量是影響NPP年際變化的主要氣象因子.

圖7 雪寶頂區域NPP與氣象因子空間相關性 Tig.7 Spatial correlation between NPP and meteorological factors in the Xuebaoding region

2.4.2 不同地表覆被類型NPP與氣象因子相關性分析 不同的地表覆被類型NPP對氣象因子的也存在明顯不同(表2).通過表2發現，各類地表覆被類型NPP與太陽輻射相關系數均在0.65以上，且通過顯著性檢驗(P<0.05)的面積比均超過50%，說明研究區各類地表覆被NPP與太陽輻射存在顯著正相關.研究區各類地表覆被與降水、溫度的相關系數均為負值；其中，沼澤草地NPP與溫度存在顯著負相關，可能與沼澤草地中分布以富養植被為主；然而研究區大多表覆被類型與降水、氣溫存在不顯著負相關，通過顯著性檢驗(P<0.05)的面積比均在7%以內，可能與該區域常年降水穩定，氣溫波動不大，使得氣象環境穩定且適宜，使得該區域大部分植被固碳量相對穩定.

3 討論

本文采用光能利用率模型(CASA)對雪寶頂區域近10a的植被凈初級生產力進行估算，并對該區域NPP的時空分布變化以及與氣象因子的相關性進行分析，由于光能利用率模型(CASA)的機理與其它模型不同，其所估算的NPP值必然有所不同，且不可避免存在一定誤差，就雪寶頂區域而言，本文對各類地表覆被的NPP值估算與前人對平武縣、北川縣等區域的模擬值與平均值基本一致[19-21]，且空間分布與美國宇航局植被凈初級生產力產品MOD17A3基本一致.

雖然本文研究表明近10a雪寶頂區域NPP以-6.862 8 gCm-2a-1的速率下降，這可能是由于2008年地震以后，大部分區域植被遭到破壞，且期間不斷出現各種地質災害，使得植被恢復受到影響有關；但是未來雪寶頂絕大部分區域呈上升趨勢，這可能與近些年各類災后生態修復工程以及該區域地表覆被類型主要以森林、草地有關.本文受數據源限制，僅選取10a的時間段進行NPP估算，未能將災前10a的時空分布與變化進行模擬分析，且在估算NPP時未能結合該區域地災較多特性將地形因子與地災影響因子加入模型之中，這是今后需要進一步研究的問題.

4 結論

本文采用CASA模型估算雪寶頂區域2008-2017年的NPP，簡單分析了其時空變化特征與氣候因子的關系，取得以下幾點結論：(1) 雪寶頂區域NPP空間分布差異較大，以雪寶頂為中界呈東高西低的分布，東部是地表覆被多為闊葉林的丘陵山地地區區域，西部是地表覆被多為草甸草原的高原丘陵地帶；(2) 近10a雪寶頂區域年均NPP值為398.34 gCm-2a-1，呈現波動下降的趨勢，以年均-6.862 8g Cm-2a-1的速率減少；(3)雪寶頂近10a呈現整體下降的趨勢，下降NPP面積占區域95.7%，主要集中在東部低海拔地區；自東向西下降率遞減，且未來大部分區域呈現上升趨勢，占區域總面積的58%；(4)在年際尺度上，太陽輻射是影響區域NPP的主要因素；且對于各類地表覆被，太陽輻射也是主要影響因素，僅有灌叢與沼澤草地存在顯著負相關.