張子凡,熊茂秋,李福杰,劉曉煌,郝玉恒,邢莉圓,王新華,賴 明,袁鵬程
(1.中國地質大學(武漢)地理與信息工程學院, 湖北 武漢 430074;2.自然資源要素耦合過程與效應重點實驗室, 北京 100055;3.中國地質大學(武漢)地理與資源學院, 湖北 武漢 430074;4.中國地質調查局烏魯木齊自然資源綜合調查中心, 新疆 烏魯木齊 830000;5.中國地質調查局自然資源綜合調查指揮中心, 北京 100055;6.中國地質調查局煙臺海岸帶地質調查中心, 山東 煙臺 264000)
內蒙古地區(qū)是我國主要的草地資源分布區(qū),也是我國的畜牧業(yè)區(qū)和糧食的主要生產地[1]。植被凈初級生產力(net primary productivity, NPP)是指綠色植物在單位時間、單位面積內所累積的有機干物質總量[2]。NPP 是生態(tài)系統(tǒng)碳匯、碳源的主要參數(shù),是陸地生態(tài)系統(tǒng)循環(huán)中必不可少的內容,也是評估生態(tài)系統(tǒng)質量與健康狀況的重要依據(jù)[3]。目前研究植被凈初級生產力的手段有實測法和模型估算法[4],實測法主要通過抽樣測量實地單位面積和時間內的生物量進行NPP 的調查和計算。隨著遙感技術的發(fā)展,利用遙感觀測數(shù)據(jù)估算區(qū)域內,長時間序列的NPP 值時空變化成為主流的研究趨勢[5]?;谶b感數(shù)據(jù)結合模型估計NPP 的方法主要有CASA模型、BEPS 模型和BIOME-BGC 模型等[6-8]。通過長時間的遙感數(shù)據(jù),結合統(tǒng)計學上的方法,對區(qū)域內NPP 值的變化趨勢和變化可持續(xù)性研究已有較多[9-12]。研究NPP 長時間序列時空變化及其影響因子,以及影響因子的時空分異性,對研究環(huán)境的變化以及不同區(qū)域碳匯的影響因素有重要意義。
目前,利用區(qū)域的遙感觀測數(shù)據(jù),估算NPP 值,結合氣候、環(huán)境等要素,研究NPP 值變化對氣候和人類活動的響應對區(qū)域內生態(tài)環(huán)境監(jiān)測的預測有重要意義,王耠熠等[13]通過計算NPP 柵格像元點在長時間序列上的變化,利用相關性和偏相關性分析方法,研究了若爾蓋高原區(qū)NPP 值與降水量和溫度之間的關系,同時重點研究了最高溫和最低溫對NPP 的 響 應。Yin 等[14]利 用Google Earth Engine 平臺研究了森林NPP 對植被類型和氣候的響應,得出降水對NPP 的影響較強,郭睿妍等[15]利用GEE 遙感數(shù)據(jù)云平臺,采取嶺回歸分析和冗余分析的方式,研究了黃河流域森林NPP 變化及其影響因子。閆妍等[16]依據(jù)研究區(qū)內不同流域之間影響因素的不同以及CASA 模型的性質,通過從流域和省域尺度,來構建NPP 值與氣候相關要素的多元回歸模型,從時間上以月尺度來研究得出湖南省地區(qū)溫度和水汽壓對NPP 的影響較大,而降水相對較小,表明了不同區(qū)域氣候因子對NPP 值的影響不同。隨著機器學習技術的發(fā)展,在統(tǒng)計學的角度上,根據(jù)實測生物量數(shù)據(jù)和實測NPP 值結合各個因子等數(shù)據(jù)采取機器學習算法中的各類訓練的方法,對區(qū)域內NPP 和生物量進行回歸預測,可以以決策樹的方式進行非線性預測,同時可以測定各因子數(shù)據(jù)的貢獻率[17-18]。然而這些研究NPP 變化的影響因子大多是以區(qū)域每個柵格點的長時間序列的NPP 和氣候數(shù)據(jù),使用統(tǒng)計學的方法,進行相關性,偏相關性,以及整個區(qū)域的最小二乘法、嶺回歸模型和機器學習中的決策樹算法進行回歸分析。進行時間序列分析時,容易忽略對于時間變化不顯著但地理變化顯著的因子,比如地形等因子。所以只通過時間序列數(shù)據(jù)變化進行回歸分析時存在局限性,而以整個區(qū)域的柵格點為對象進行全局回歸分析,會存在空間局部不平穩(wěn)性的問題,無法反映出NPP 分布和影響因素的空間異質性。由于位于地理接近的區(qū)域,物體空間分布存在一定的自相關性和聚集性,內蒙古NPP 區(qū)域自相關性顯著[19],研究區(qū)域內NPP 時空變化的驅動因子時,進行區(qū)域全局分析多元線性回歸OLS 模型以及采用決策樹模型進行分線性回歸分析存在局部聚集等局部不平穩(wěn)性問題。而時空地理加權回歸模型,通過核函數(shù)確定時空距離的權重矩陣,可以較好地反映影響因子時空非平穩(wěn)性的問題,可以從時間和空間的角度研究驅動因素[20-21],這樣可以分析不同區(qū)域內各因子不同的回歸系數(shù)。
本文研究基于2000-2019 年內蒙古草原自然資源大區(qū)的MOD17A3-NPP 數(shù)據(jù)、landsat7,8 的NDVI數(shù)據(jù)和國家氣象網(wǎng)站的氣象站點數(shù)據(jù)。采用統(tǒng)計分析的方法,研究該區(qū)域的NPP 時空變化及其可持續(xù)性。利用時空地理加權回歸模型,研究不同時空尺度下NPP 值的驅動因子,得出不同時空尺度的回歸系數(shù),量化不同區(qū)域NPP 變化的影響因素,同時也可為NPP 的預測提供幫助。
研究區(qū)域為全國自然資源綜合區(qū)劃結果中的華北草耕自然資源大區(qū)的北部區(qū)域[22-23],命名為內蒙古草原資源大區(qū),該大區(qū)是我國草地資源分布的主要區(qū)域(37°37′12″~51°02′42″ N, 106°28′40″~121°28′03″ E),主要包括河套平原和內蒙古高原、黃土高原的部分地區(qū),位于我國內蒙古自治區(qū),面積約6.51 × 108hm2,主要的自然資源為草原,約占63.42% (圖1)。
圖1 研究區(qū)概況圖Figure 1 Overview of the study area
本研究主要使用了NPP 數(shù)據(jù)、氣象數(shù)據(jù)、NDVI數(shù)據(jù)、高程數(shù)據(jù)和土地利用數(shù)據(jù)。NPP 數(shù)據(jù)來源于2000-2019 年的MODIS 陸地4 級標準數(shù)據(jù)NPP 產品數(shù)據(jù)集MOD17A3-NPP 數(shù)據(jù),來自美國航空航天局網(wǎng)站(https://ladsweb.modaps.eosdis.nasa.gov/),其空間分辨率為500 m,坐標系為WGS_1984 坐標系,投影為UTM 投影。氣象數(shù)據(jù)來自于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn),選取了內蒙古、陜西、山西和河北省123 個氣象站點2000-2019 年的降水量、平均氣溫、相對濕度、光照時長數(shù)據(jù),采用克里金插值法進行插值,得到與NPP 數(shù)據(jù)投影坐標系相同和分辨率相同的年尺度降水、平均氣溫、相對濕度、光照時長的柵格數(shù)據(jù)。NDVI 數(shù)據(jù)來自美國陸地衛(wèi)星landsat7,8 數(shù) 據(jù) 的 反 演(https://landsat.gsfc.nasa.gov/)。高程數(shù)據(jù)和土地利用數(shù)據(jù)來自中科院資源環(huán)境數(shù)據(jù)中心(https://www.resdc.cn/),選取1 km 分辨率的土地數(shù)據(jù)和dem 數(shù)據(jù)。
2.2.1 Theil-Sen Median 方法和Mann-Kendall 顯著性檢驗
通過Theil-Sen Median 方法植被凈初級生產力的趨勢變化,利用Mann-Kendall 方法進行變化的顯著性檢驗[24]。該方法可以減少外部帶來的誤差,同時對趨勢的判斷受異常值的影響較小,有非常高的適用性,可以更準確反映長時間序列下植被凈初級生產力的變化。通過這個方法研究區(qū)域內每一個像元的時間序列變化,確定2000-2019 年NPP 的時空變化特征。Theil-Sen Median 方法公式如下:
式中:Z為NPP變化的顯著性,n為研究的時間序列,tp為每年重復的NPP值重復的個數(shù)。定義α為顯著性檢驗水平,Z(1-α)/2為標準正態(tài)方差,當|Z| >Z(1-α)/2時,表示時間序列ET 在α水平上存在顯著變化。
2.2.2 Hurst 指數(shù)
Hurst 指數(shù)是用來評估時間序列信息長時間依賴的方法[25]。在本研究中,用Hurst 指數(shù)來研究植被凈初級生產力變化的可持續(xù)性,可以判斷未來變化的趨勢。具體公式和計算步驟如下:
定義研究區(qū)間NPP 變化的時間序列:
對Hurst 指數(shù)進擬合:
Hurst 指數(shù)(H)在0.5 <H< 1 時,表明時間序列是一個可持續(xù)序列,即未來變化趨勢與現(xiàn)在變化趨勢一致,H越接近1,持續(xù)性變化越強。H= 0.5 則說明NPP時間序列為隨機序列,不具備時間相關性,0 <H< 0.5 則表明未來變化趨勢與現(xiàn)在有相反的趨勢,H越接近0,反持續(xù)性變化越強。
2.2.3 時空地理加權回歸(GTWR)模型
基于最小二乘法的多元線性回歸模型可反映全局的變量因子與因變量的關系,但對于局部聚集性比較強的數(shù)據(jù),利用全局分析無法反映因子的局部不平穩(wěn)特征。地理加權回歸模型可以反映局部空間的不平穩(wěn)性,時空地理加權回歸模型是在地理加權回歸模型基礎上,引入了時間維度[26],與地理加權回歸模型相比,能更好地揭示時空非平穩(wěn)性變化規(guī)律,不僅能夠反映空間變化趨勢,還可以提高對時間數(shù)據(jù)演變規(guī)律的認知,從而增強對時空數(shù)據(jù)的分析能力[27]。本文選擇2000-2019 年的數(shù)據(jù),研究NPP值與降水、氣溫、濕潤度、光照時長、高程、植被覆蓋度(NDVI)的相互關系,計算回歸系數(shù),同時也能根據(jù)回歸系數(shù)依據(jù)氣候NDVI 等值對未來NPP 值進行預測。時空地理加權回歸模型公式如下:
統(tǒng)計2000-2019 年各像元點NPP 均值的空間分布特征(圖2),并計算各年份區(qū)域內NPP 的均值,2000-2019 年區(qū)域NPP 平均值約為208.96 g·(m2·a)-1,NPP 值范圍為0~727.0 g·(m2·a)-1。NPP 較高的區(qū)域主要位于東部和南部,地表資源主要為森林資源和耕地資源,較低的區(qū)域主要位于西北部的荒漠地區(qū)(圖1、圖2)。從時間序列上看,2000-2019 年,內蒙古草原自然資源大區(qū)NPP 均值波動較為明顯,但整體上呈現(xiàn)較快的上升趨勢,平均每年上升約4.06 g·(m2·a)-1,NPP 均 值 在2001 年 最 低,僅 約150.55 g·(m2·a)-1,在2019 年 最 高,NPP 均 值 約 為247.79 g·(m2·a)-1。在2000-2019 年,2001-2002、2007-2008、2011-2012 和2017-2018 年,區(qū)域內NPP 均值急劇上升,在2000-2001、2006-2007 和2014-2015年下降較為明顯,其余年份變化較為平緩(圖3)。
二是松花江黑龍江洪水量級大,持續(xù)時間長。松花江發(fā)生1998年以來最大流域性較大洪水,其中嫩江上游發(fā)生超50年一遇特大洪水,松花江上游發(fā)生超20年一遇大洪水。黑龍江發(fā)生1984年以來最大流域性大洪水,下游發(fā)生超100年一遇特大洪水。嫩江、松花江干流水位超警戒歷時46天,黑龍江干流水位超警戒歷時58天。
圖3 2000-2019 年內蒙古草原自然資源大區(qū)植被凈初級生產力均值及變化趨勢圖Figure 3 Mean net primary productivity (NPP) values and their change trend for natural resources in the grassland region of Inner Mongolia from 2000 to 2019
基于2000-2019 年內蒙古草原自然資源大區(qū)的NPP 柵格數(shù)據(jù),采用Theil-sen 斜率趨勢分析方法和Mann-Kendall 顯 著 性 檢 驗 方 法,取|Z| > 1.645、1.96 和2.576,表示研究序列通過了置信度為90%、95%和99%的顯著性檢驗,分別對應為弱顯著、顯著和極顯著。得到NPP 值變化趨勢以及變化顯著性的空間分布圖(圖4a),大部分區(qū)域NPP 值呈上升趨勢,約占整個區(qū)域的94.01%,僅5.99%的區(qū)域呈現(xiàn)減少趨勢。并且大部分區(qū)域NPP 值增加趨勢表現(xiàn)為極顯著增加,約占總研究區(qū)面積的60.31%,其次為顯著增加,約占研究區(qū)域的19.57%。NPP 呈減少趨勢的區(qū)域零星分布,其中不顯著減少所占的區(qū)域最多,約5.88%,極顯著、弱顯著和顯著減少的區(qū)域僅占約0.15%,零星分布于建設用地相關的區(qū)域,僅有約1.01 × 105hm2。
圖4 2000-2019 年植被凈初級生產力值變化趨勢和變化顯著性(a)以及可持續(xù)性(b)Figure 4 Net primary productivity change trend and significance (a) and persistence of the change trend (b) from 2000 to 2019
將NPP 的變化趨勢與Hurst 指數(shù)相疊加得出NPP 值的變化趨勢和持續(xù)形式,分為持續(xù)減少(2000-2019 年呈現(xiàn)減少趨勢,H> 0.5)、潛在增加(2000-2019 年呈減少趨勢,H< 0.5)、無變化趨勢(H= 0.5)、潛在減少(2000-2019 年呈現(xiàn)增加趨勢,H< 0.5)、持續(xù)增加(2000-2019 年呈現(xiàn)增加趨勢,H> 0.5)。由圖4b 可知,在2000-2019 年,盡管內蒙古草原資源大區(qū)94%的區(qū)域呈上升趨勢,但未來潛在減少的區(qū)域最多,約占整個區(qū)域面積的74.15%,分布于該區(qū)域西部,南部的大部分地區(qū),其次為未來持續(xù)增長的區(qū)域,約占整個區(qū)域的19.45%,主要分布于該大區(qū)的東部和東北部的有林地地區(qū)。處于不確定游離狀態(tài)的區(qū)域僅有約5.83%,主要分布在荒漠和濕地地區(qū)。未來NPP 呈現(xiàn)持續(xù)減少的區(qū)域僅有約0.23%,零星分布于該大區(qū)的中部。
通過時空地理加權回歸模型,選取2000、2005、2010、2015、2019 年的降水、年平均氣溫、光照時長、相對濕度等氣候因子以及海拔高度地形因子和NDVI 植被覆蓋度因子。通過Arcgis 漁網(wǎng)工具全局選取了1 056 個采樣點,研究區(qū)域內NPP 值時空變化與各因子之間的關系,然后進行點的插值得出植被凈初級生產力值變化和各影響要素回歸系數(shù)分布(圖5、圖6)。最終通過時空地理加權回歸模型(GTWR)得出的模擬擬合度R2= 0.925,該結果表明回歸精度較好,可以進行分析研究。
圖6 2019 年植被凈初級生產力值變化、各影響要素回歸系數(shù)分布Figure 6 Changes in NPP of vegetation in 2019, distribution of regression coefficients of various influencing factors
統(tǒng)計區(qū)域2000-2019 各影響因子得到平均回歸系數(shù),2000-2019 年降水、相對濕度和NDVI 變化和NPP 值相關性較強,主要為正向相關關系,氣溫、光照時長和蒸散發(fā)主要為負向相關關系(表1)。在2000 年,降水對NPP 值變化的回歸系數(shù)最高值為1.113,平均值為0.368,呈負相關的區(qū)域約占5%,在空間分布上,主要分布在西南少部分地區(qū),其他區(qū)域呈正相關,影響程度從西南到東北存在明顯的空間分異性,從西南到東北地區(qū)呈現(xiàn)線性增加的趨勢;2000 年,氣溫對NPP 影響主要為負相關,回歸系數(shù)的平均值為-0.027,最高值為0.803,最低值為-0.631,其中正向相關影響的主要分布于該區(qū)錫林郭勒盟市的西部,占整個研究區(qū)域的20%,其他區(qū)域呈負向相關關系,并且存在明顯的空間異質性,在研究區(qū)域西南影響程度較大,東北部次之,中部最小;光照時長對NPP 值的影響相對較小,在2000 年影響平均回歸系數(shù)為-0.061,最高值為-0.051,最低值為-0.066,幾乎所有區(qū)域呈負相關性影響,影響程度從西南到東北有顯著的線性減小趨勢;濕潤度對NPP 值的影響在2000 年影響平均回歸系數(shù)為0.228,最高值為0.233,最低值為0.221,幾乎所有區(qū)域呈現(xiàn)正向相關,影響程度從西南到東北呈線性增加的趨勢;2000 年蒸散發(fā)對NPP 的影響平均回歸系數(shù)為-0.020,最高值為-0.030,影響程度由西南到東北區(qū)域逐漸遞減;NDVI 值與NPP 值主要是正向相關關系,且回歸系數(shù)較高,其中在西南地區(qū)最高,由西南向東北呈遞減趨勢;2000 年海拔高度與NPP 值回歸系數(shù)的平均值為0.080,最低值為-0.032,最高值為0.177,其中海拔高度與NPP 呈現(xiàn)負相關性的區(qū)域占37.45%,主要是在該大區(qū)的西南部,其他區(qū)域呈現(xiàn)正相關,影響程度從該區(qū)西南到東北有減小的趨勢(圖5)。
表1 2000-2019 年植被凈初級生產力影響因子回歸系數(shù)Table 1 Regression coefficients for net primary productivity influencing factors, 2000-2019
圖5 2000 年植被凈初級生產力值變化、各影響要素回歸系數(shù)分布Figure 5 Changes in NPP of vegetation in 2000, distribution of regression coefficients of various influencing factors
2019 年,降水、溫度、光照時長、濕度、蒸散發(fā)等氣候因子對NPP 值影響的回歸系數(shù)的平均值分別為0.257、-0.145、-0.201、0.081 和-0.001,海拔和NDVI 值的回歸系數(shù)分別為0.078 和0.559。與2000相比,降水、濕度和蒸散的影響程度在減少,光照、溫度和NDVI 值的影響程度有所增加,海拔影響基本不變(表1)。2019 年,降水對NPP 值的影響的回歸系數(shù)最高值為0.715,最低值為-0.158,其中呈負相關的區(qū)域,約占10%,主要在呼倫貝爾市西北和巴彥淖爾市東南,呈正相關的區(qū)域影響程度由西南向東北呈增加趨勢;氣溫對NPP 的影響,最高值為0.901,最低值為-1.109,大部分區(qū)域呈現(xiàn)負相關影響,約占59.56%,呈正相關的區(qū)域主要在該大區(qū)東北部和巴彥淖爾市西部;2019 年,光照時長對NPP 的影響回歸系數(shù)最高值為-0.197 4,最低值為-0.202 4,平均值為-0.201,所有區(qū)域光照時長與NPP 值呈現(xiàn)負相關影響,影響程度從西南到東北呈線性減弱;相對濕度對NPP 值的回歸系數(shù)最高為0.077,最低為0.084,影響程度較小,全部區(qū)域濕度與NPP 值呈正相關影響,與2000 年不同,影響程度由西南到東北逐漸減弱;2019 年地表蒸散對NPP 值影響的回歸系數(shù)最高值為-0.08,最低值為-0.01,影響程度變化不大,由西南到東北遞減;2019 年NDVI值和NPP 值的相關系數(shù)最低值為0.08,平均值為0.559,所有區(qū)域呈現(xiàn)正相關關系,其中在東北部和西南部的部分地區(qū),NDVI 值對NPP 的影響較小,影響程度存在顯著的聚集性;海拔高度對NPP 值的影響與2000 有所不同,最低值為0.048,最高值為0.118,全區(qū)域呈正向關,有顯著的空間異質性,由西南向東北遞增(圖6)。
本研究通過對2000-2019 年NPP 值的變化及其驅動因子研究得出,降水量變化是內蒙古草原自然資源大區(qū)NPP 值增加的主要氣候因素,NDVI 值與NPP 值密切相關,但NDVI 值越高的地區(qū),NPP受NDVI 變化的影響相對較小。另外內蒙古草原自然資源大區(qū)大部分區(qū)域溫度和光照時長的變化與NPP 值呈負相關關系,同時也有學者研究指出一定程度的氣溫升高可以提高光合效率和資源利用效率從而增加NPP[29-30],可能因為在內蒙古半干旱地區(qū),由于水分較少,溫度升高會引起水分的減少,從而使植被利用水分減少反而降低了光合效率,地表水分是影響NPP 變化的主要原因。越干旱的地區(qū)溫度、光照時長和蒸散發(fā)對NPP 值的負相關性就越強(圖5、圖6)。而NPP 與溫度呈正相關的區(qū)域主要在東北部和錫林郭盟市西部植被覆蓋指數(shù)較高的地方。各類要素對NPP 值的影響從該大區(qū)西南到東北有顯著的線性變化,有明顯的空間異質性,可以為區(qū)域碳儲量和自然資源的分區(qū)研究提供幫助。
目前本研究還存在許多不足,主要對NPP 值時空變化的趨勢進行分析。以及利用氣候、海拔、植被覆蓋度等因子對NPP 值進行時空地理加權回歸分析,通過回歸系數(shù)研究各因子對NPP 值的影響。主要還是以統(tǒng)計學以及統(tǒng)計學相關的模型進行分析,沒有深入研究內部變化機理,所得到的影響系數(shù)只是統(tǒng)計學上的相關性。由于2000 年和2019 年的氣溫和植被類型存在較大差異,有些區(qū)域的NPP值的影響回歸系數(shù)在2000 年與2019 年相差較大,甚至相關趨勢都不同,這可能是不同氣候對NPP 的影響存在一個閾值,這個閾值需要大量的時空數(shù)據(jù)分析,可能是以后需要研究的方向。另外NPP 的變化可能與人類的生產活動有很大相關性。綜合區(qū)域自然因子和工農業(yè)、畜牧業(yè)、人口密度、GDP 等社會經(jīng)濟因子研究NPP 的變化是今后需要解決的問題。
在2000-2019 年間內蒙古草原自然資源區(qū)域NPP 均值呈波動上升趨勢,平均每年增加約207.72 g·(m2·a)-1。在2000-2019 年94.01%的區(qū)域NPP 值呈上升趨勢,79.88%的區(qū)域NPP 值呈顯著上升(|Z| >1.96),呈顯著下降的區(qū)域極少,僅有約0.15%;利用Hurst 指數(shù)來判斷每個像元點變化的可持續(xù)性,得出盡管2000-2019 年有94.01%的區(qū)域NPP 值呈現(xiàn)上升趨勢,但在未來會持續(xù)增加的區(qū)域僅占整個研究區(qū)域的17.45%主要分布在研究區(qū)域東部和東北部地區(qū),74.15%的區(qū)域在未來有潛在下降的趨勢。
采用時空地理加權回歸模型,研究NPP 值變化的影響因子。得出在2000-2019 年,該研究區(qū)域內降水、相對濕度和NDVI 值的影響主要呈正相關,降水和NDVI 值對NPP 的影響較高,在2000 和2019年相關系數(shù)分別為0.368、0.257 和0.402、0.259。溫度、蒸散發(fā)和光照時間主要呈負相關關系。其中溫度和光照時間影響相對較大,2000 年和2019 年回歸系數(shù)平均值分別為-0.027、-0.145 和-0.061、-0.201。蒸散發(fā)影響較小,2000 年和2019 年回歸系數(shù)的平均值分別為-0.02 和-0.001。
各因子對NPP 的影響存在空間異質性。降水與NPP 值的變化主要呈正相關關系,但在2000 年和2019 年仍然有一部分區(qū)域呈現(xiàn)負相關,分別占5%和10%,主要在該大區(qū)西南和東北的部分區(qū)域。溫度主要是負相關,2000 和2019 呈正相關的區(qū)域主要在錫林郭盟西部;另外光照時長、相對濕度、蒸散量和高程對NPP 值的影響在不同區(qū)域也有所不同,從該區(qū)西南到東北有顯著的線性變化趨勢。