基于耦合協(xié)調度的黃土高原地區(qū)NDVI與降水關系的變異診斷

賈 路,于坤霞，*,徐國策,2,任宗萍,高海東,李占斌,李鵬,2

1 西安理工大學 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室,西安 710048 2 西安理工大學 旱區(qū)生態(tài)水文與災害防治國家林業(yè)局重點實驗室,西安 710048

當前,黃土高原地區(qū)水土流失依然嚴重[1],生物多樣性遭到破壞,由于氣候變化加劇[2-3],該地區(qū)面臨極端水文事件頻發(fā)帶來的一系列環(huán)境問題[4-7],生態(tài)環(huán)境治理遭受嚴重挑戰(zhàn)。隨著我國社會經濟的發(fā)展,人民對美好生活更加向往,生態(tài)環(huán)境保護刻不容緩[8]。黃河流域是中華文明的重要誕生地,在氣候變化和人類活動的共同作用下,流域生態(tài)環(huán)境與自然條件顯著變化,對當?shù)厝嗣竦纳詈蜕a帶來巨大的影響,黃河流域的健康和高質量成為一項重大國家戰(zhàn)略[9-10]。黃河中游地區(qū)主要位于中國黃土高原地區(qū),該區(qū)域水土流失極為嚴重,生態(tài)脆弱,環(huán)境惡劣,是入黃泥沙的主要來源地[11-14],因此黃土高原地區(qū)的水土流失治理工作關系到黃河流域的健康發(fā)展。

黃河流域治理的關鍵問題是水沙變化[15]。從20世紀50年代開始,大規(guī)模的水土保持工程在黃土高原地區(qū)進行了實施,包括林草措施等[16]。特別是在1999年,退耕還林政策在黃土高原地區(qū)進行了大規(guī)模實施[17-20],區(qū)域植被顯著增加,有效的減輕流域土壤侵蝕,黃河流域泥沙顯著減少[16,21]。植被在陸-氣系統(tǒng)之間發(fā)揮著能量傳輸、維持和優(yōu)化生態(tài)系統(tǒng)服務的重要作用,是全球生態(tài)環(huán)境變化的重要“指示器”[22-23]。隨著全球氣候變化與人類活動影響對區(qū)域發(fā)展影響的課題受到學術界的高度關注,陸地生態(tài)系統(tǒng)與氣候變化的關系研究成為當前國際研究的熱點科學問題之一[24-25]。由于歸一化植被指數(shù)(NDVI)能夠指示大尺度的區(qū)域植被覆蓋狀況,表征植被活動的強弱,因此常被廣泛應用于研究生態(tài)變化的研究中[26-27]。許多研究成果表明植被變化對氣候因子變化具有明顯的響應,在區(qū)域尺度上氣候條件和背景環(huán)境的空間異質性對植被變化具有決定作用[28-29],干旱區(qū)的植被凈初級生產力(NPP)更容易在降水的影響下增加[30]。雖然過去的研究已經取得了豐富的成果,但是以往氣候因子與植被變化的研究主要關注的是氣候變化的趨勢或者線性變化對植被的影響,缺乏關于氣候因子和植被之間二元關系變化的理論和模型等研究。有一個至關重要事實是,由于過去對于植被變化的研究大多數(shù)是基于單變量統(tǒng)計學手段對其進行分析,例如變異系數(shù)和Mann-Kendall趨勢檢驗等,對于植被與氣候因子的響應關系多基于線性方法進行簡單的研究,例如線性相關檢驗、線性回歸模型等等,所以這容易造成變量之間的變化存在巨大的差異。自然界中的地表水和地下水,植被,大氣,土壤等環(huán)境要素之間互相作用,非線性關系極為復雜[31],使用單變量的分析方法有利于簡化模型,減輕工作量。然而,單變量的分析方法只關注了單個變量的變化(例如降水或植被),這不符合現(xiàn)實的具體情況,很難完全揭示對氣候變化和植被變化之間的響應關系。因此,有必要構建植被與降水之間的雙變量耦合關系理論,診斷植被與降水之間變異特征。

本研究的主要目的是,通過Mann-Kendall趨勢檢驗和Hurst指數(shù)分析黃土高原地區(qū)1998—2017年NDVI的變化趨勢與持續(xù)性,利用RDR指數(shù)與重心轉移模型分析黃土高原地區(qū)NDVI的時空變化差異,基于耦合協(xié)調度理論和Pettitt檢驗方法構建區(qū)域NDVI和降水關系的突變點的識別方法并進行診斷,探討了造成黃土高原地區(qū)NDVI和降水關系變化的可能驅動因素。研究成果有助于進一步深化對黃土高原地區(qū)植被變化和降水關系變化的認識,有利于為黃土高原水土保持工作開展提供一定的科學依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

黃土高原地區(qū)作為世界上黃土覆蓋面積最大的高原(34°—42°N,101°—114°E),黃土平均厚度在50—80 m,部分地區(qū)在100 m以上,總面積大約為64萬km2。黃土高原地區(qū)位于屬于我國地勢上的第二階梯,海拔大約81—4918 m,平均海拔為1404 m。區(qū)域內主要包括甘肅省、寧夏回族自治區(qū)、青海省、內蒙古自治區(qū)、山西省及河南省等省區(qū)的46個地區(qū)(盟、州、市),282個縣。黃土高原地區(qū)河水的主要來補給源是降水,該區(qū)域干旱缺水。黃土高原地區(qū)特有的氣候特征和濕陷性黃土以及歷史上的植被破壞,使得該地區(qū)水土流失異常嚴重,溝壑交錯,流域支離破碎,生態(tài)環(huán)境脆弱。黃土高原地區(qū)的地貌分區(qū)大致可分為4種,依次為北部風沙區(qū)和農灌區(qū)、高原溝壑區(qū)、丘陵溝壑區(qū)和東部河谷及土石山區(qū),如圖1所示。

圖1 黃土高原地區(qū)地理位置與地貌分區(qū) Fig.1 Geographical location and geomorphic zoning of the Loess Plateau

1.2 數(shù)據(jù)來源

本文所使用的降水數(shù)據(jù)是來源于國家地球系統(tǒng)科學數(shù)據(jù)中心的柵格降水量數(shù)據(jù)(http://loess.geodata.cn),時間序列為1901—2017年,空間分辨率為1 km×1 km,時間分辨率為月,該數(shù)據(jù)使用全國496個獨立氣象觀測點數(shù)據(jù)進行了驗證,驗證結果可信。本研究通過ArcGIS軟件裁剪了黃土高原地區(qū)1998—2017年逐月的降水量數(shù)據(jù),并通過對逐月數(shù)據(jù)進行合成,獲得了黃土高原地區(qū)1998—2017年逐年年降水量柵格數(shù)據(jù)。植被歸一化指數(shù)(NDVI)是一個長期的SPOT/VEGETATION NDVI衛(wèi)星遙感數(shù)據(jù)集,數(shù)據(jù)集由中國科學院地理科學與資源研究所網站提供(http://www.resdc.cn),數(shù)據(jù)時間為1998—2017年,空間分辨率為1 km×1 km,時間分辨率為年。通過使用ArcGIS軟件將NDVI數(shù)據(jù)與柵格降水量數(shù)據(jù)進行投影轉換,獲得逐年相互匹配的NDVI數(shù)據(jù)和柵格降水量數(shù)據(jù)。

1.3.1Mann-Kendall趨勢檢驗與Hurst指數(shù)

Mann-Kendall(M-K)趨勢檢驗[32]作為一種非參數(shù)的時間序列趨勢檢驗方法,由于簡單實用常被廣泛使用在氣象和水文序列的趨勢檢驗中,本文中使用Mann-Kendall趨勢檢驗識別降水量、NDVI以及其他變量的變化趨勢。

Hurst指數(shù)(H)常被用在定量表征時間序列的持續(xù)性或長期相關性[33],本研究中使用Hurst指數(shù)來描述NDVI的未來變化趨勢。其中H代表Hurst指數(shù)。H的值在[0,1]的范圍內。當H=0.5時,表明NDVI的時間序列是隨機序列,是不可持續(xù)的。當H>0.5時,NDVI的變化與目前的趨勢基本一致,表明植被的可持續(xù)性是正向的。H<0.5表示負可持續(xù)性,未來植被變化將與當前趨勢相反。

1.3.2RDR指數(shù)

相對發(fā)展率(RDR)被用來分析研究區(qū)內[34],不同柵格NDVI在研究時段初期與末期的變化量與研究區(qū)整體NDVI在研究時段初期與末期變化量的變化差異,計算公式如下:

(1)

式中:Y2i為2017年黃土高原地區(qū)第i個柵格NDVI；Y1i為1998年黃土高原地區(qū)第i個柵格NDVI；Y2為2017年黃土高原地區(qū)NDVI平均值；Y1為1998年黃土高原地區(qū)NDVI平均值。

1.3.3重心轉移模型

重心轉移模型[34]常被用來描述地理變量分布重心在空間上的變化特征,可以被用來研究NDVI空間重心隨時間的變化,重心轉移模型在本研究中被用來研究黃土高原地區(qū)1998—2017年NDVI的重心變化,計算公式如下:

(2)

(3)

式中:Xm和Ym為NDVI空間分布重心的經度和緯度；cmi為第i個柵格的NDVI；xi為第i個柵格的經度；yi為第i個柵格的緯度。

1.3.4耦合協(xié)調度理論

兩個或兩個以上的系統(tǒng)或運動形式通過各種相互作用而彼此影響的現(xiàn)象被稱為耦合[35-37],隨著科學理論的深入研究與應用,這一理論逐漸被廣泛使用在生態(tài)學領域。NDVI和降水兩個子系統(tǒng)之間的耦合度模型,表示成以下公式[38]:

(4)

式中:C為NDVI和降水二元系統(tǒng)的耦合度,u1,u2分別為NDVI子系統(tǒng)和降水系統(tǒng)對總系統(tǒng)有序度的貢獻。在本研究中NDVI和降水子系統(tǒng)的指標為黃土高原地區(qū)所有柵格1998—2017的NDVI和降水量,假定NDVI和降水量越大對系統(tǒng)越好。

(5)

(6)

式中ui為第i個子系統(tǒng)對總系統(tǒng)有序度的貢獻；uij為第i個子系統(tǒng)中第j個指標的歸一化值；wij為第i個子系統(tǒng)中第j個指標的權重,每個子系統(tǒng)中指標的權重計算使用熵權法進行計算。

在計算每個子系統(tǒng)的熵權時,必須通過進行歸一化先對數(shù)據(jù)處理,常用的方法是極值法:

當uij越大對系統(tǒng)越有利(正向歸一化):

(7)

當uij越小對系統(tǒng)越有害(負向歸一化):

(8)

式中uij為子系統(tǒng)第i個指標的第j個時序的歸一化值；xij為子系統(tǒng)第i個指標的第j個時序的值。

為了更好的反映出子系統(tǒng)整體“功效”與“協(xié)同”效應,耦合協(xié)調度模型用來評判NDVI和降水兩系統(tǒng)的交互耦合的協(xié)調程度,計算公式如下:

D=(C×T)1/2

(9)

T=a×u1+b×u2

(10)

式中,D為耦合協(xié)調度；C為耦合度；T為NDVI子系統(tǒng)與降水子系統(tǒng)的綜合調和指數(shù),它反映NDVI子系統(tǒng)與降水子系統(tǒng)的整體協(xié)同效應或貢獻；a、b為待定系數(shù),實際中常常認為兩個子系統(tǒng)的重要性相同,所以a=b=0.5。

1.3.5Pettittt檢驗

Pettitt檢驗法[39]采用Mann-Whitney中Ut,n值檢驗同一總體中兩個樣本X1,…,Xt和Xt+1,…,XN,Pettitt檢驗的零假設為沒有變化點,當|Ut,n|取最大值時對應的Xt被認為是可能的突變點。其顯著性水平可由下式計算:

(11)

當P≤0.05時認為數(shù)據(jù)中存在均值變異點。

2 結果與分析

2.1 NDVI變化特征

黃土高原地區(qū)NDVI在1998—2017年73.49%面積呈現(xiàn)顯著增加趨勢(P<0.05),變化趨勢時空分布如圖2所示。北部風沙區(qū)和農灌區(qū)、高原溝壑區(qū)、丘陵溝壑區(qū)和東部河谷及土石山區(qū)的NDVI均有大面積呈現(xiàn)顯著增加趨勢(P<0.05),4個分區(qū)NDVI顯著增加面積占各分區(qū)面積的比例為54.32%、72.91%、93.97%和74.98%(P<0.05)。北部風沙區(qū)和農灌區(qū)、高原溝壑區(qū)、丘陵溝壑區(qū)和東部河谷及土石山區(qū)的多年平均NDVI分別為0.34、0.53、0.55和0.58,東部河谷及土石山區(qū)的植被覆蓋狀況最好,北部風沙區(qū)和農灌區(qū)的植被覆蓋狀況最差(圖2)。1998—2017年黃土高原NDVI的Hurst指數(shù)空間分布最大值為0.75,最小值為0.25,黃土高原只有5.84%面積NDVI的Hurst指數(shù)低于0.5,主要集中在北部風沙區(qū)和農灌區(qū)。黃土高原大部分區(qū)域NDVI在未來依舊呈現(xiàn)增加趨勢。

圖2 黃土高原地區(qū)1998—2017年NDVI變化趨勢Fig.2 Trend of annual NDVI in the Loess Plateau from 1998 to 2017

2.2 NDVI時空變化差異

黃土高原地區(qū)NDVI在1998—2017年間變化存在時空差異,如圖3所示。RDR指數(shù)最大值為8.14,最小值為-6.42,RDR指數(shù)大于1的區(qū)域主要集中在丘陵溝壑區(qū),RDR指數(shù)小于1的區(qū)域主要集中在北部風沙區(qū)和農灌區(qū)以及黃土高原邊界區(qū)域。在高原溝壑區(qū)也有較大面積NDVI的RDR指數(shù)大于1,這表明在1998—2017年間黃土高原地區(qū)NDVI在丘陵溝壑區(qū)與高原溝壑區(qū)的增加幅度遠遠大于黃土高原整體的增加幅度,北部風沙區(qū)和農灌區(qū)以及黃土高原邊界區(qū)域的NDVI增加相對總體來說更為滯后。

圖3的結果顯示,1998—2017年間黃土高原NDVI重心在不斷發(fā)生轉移。從1998—1999年、1999—2000年、2000—2001年、2001—2002以及2015—2016年轉移距離比較長,這表明這期間植被在黃土高原地區(qū)大范圍增加。從總體看,NDVI重心由1998年的109.36°E 35.79°N變化為2017年的109.36°E 35.81°N,呈現(xiàn)從波動變化的趨勢,表明黃土高原地區(qū)植被在年際間存在明顯的變化差異。

圖3 1998—2017年黃土高原地區(qū)NDVI的RDR指數(shù)和NDVI重心轉移Fig.3 RDR index of annual NDVI in the Loess Plateau and shifting of the center of gravity of annual NDVI

2.3 NDVI子系統(tǒng)與的降水子系統(tǒng)權重分布變化

根據(jù)耦合協(xié)調度理論,構建了黃土高原地區(qū)1998—2017年間NDVI子系統(tǒng)與降水子系統(tǒng)之間的二元耦合系統(tǒng),NDVI子系統(tǒng)與降水子系統(tǒng)的構建指標分別為黃土高原地區(qū)1998—2017年每個柵格的NDVI時間序列和降水量時間序列,使用熵權法計算了每個指標的權重如圖4所示。NDVI的權重分布最大值為90.6×10-7,最小值為4.7×10-7,北部風沙區(qū)和農灌區(qū)、高原溝壑區(qū)、丘陵溝壑區(qū)和東部河谷及土石山區(qū)的NDVI子系統(tǒng)權重均值分別為17.9×10-7、16.66×10-7、14.29×10-7和15.05×10-7,高原溝壑區(qū)和北部風沙區(qū)和農灌區(qū)的NDVI變化對黃土高原地區(qū)整體的重要性較大(圖4)。圖4表明北部風沙區(qū)和農灌區(qū)、高原溝壑區(qū)、丘陵溝壑區(qū)和東部河谷及土石山區(qū)的降水子系統(tǒng)權重均值分別為13.88×10-7、19.50×10-7、13.95×10-7和15.21×10-7,高原溝壑區(qū)和東部河谷及土石山區(qū)的降水變化對黃土高原地區(qū)整體的重要性較大。

圖4 NDVI子系統(tǒng)權重分布和降水子系統(tǒng)權重分布Fig.4 Distribution of NDVI subsystem weights and distribution of precipitation subsystem weights

2.4 NDVI與降水的耦合協(xié)調度變化及其突變點

根據(jù)耦合協(xié)調理論計算,1998—2017年間黃土高原地區(qū)NDVI與降水的耦合協(xié)調度變化如圖5所示。由M-K趨勢檢驗分析可知,NDVI與降水之間的耦合協(xié)調度呈現(xiàn)顯著增加趨勢,1998年耦合協(xié)調度為0.61,2017年耦合協(xié)調度為0.81,增加幅度為32.79%。由Hurst指數(shù)分析可知,1998—2017年間黃土高原地區(qū)NDVI與降水的耦合協(xié)調度Hurst指數(shù)為0.65,未來依然呈現(xiàn)增加趨勢,這表明黃土高原地區(qū)NDVI與降水之間的耦合協(xié)調程度在持續(xù)增強。由于耦合協(xié)調度可以反映黃土高原地區(qū)NDVI與降水之間的相互作用關系,因此通過Pettitt方法識別了NDVI與降水的耦合協(xié)調度變化的突變點。1998—2017年間黃土高原地區(qū)NDVI與降水的耦合協(xié)調度在2006年發(fā)生了顯著突變(P<0.05),這表明黃土高原地區(qū)NDVI與降水之間的關系可能在2006年存在變異。

圖5 NDVI與降水的耦合協(xié)調度變化及其突變點Fig.5 Change of coupling coordination degree of NDVI and precipitation and its change-point

3 討論

3.1 降雨對NDVI的影響

降水作為重要的氣候因子,對植被的生長和變化具有不可忽視的作用[34,40-41]。許多研究成果表明,降水可以促進植被更好的生長和存活[28]。圖6所示為1998—2017年黃土高原地區(qū)NDVI顯著增加區(qū)域(P<0.05)和不顯著變化區(qū)域(P>0.05)逐年平均年降水量的時間變化過程。NDVI顯著增加區(qū)域的降水在1998—2017年間均高于NDVI變化不顯著的區(qū)域,這說明降水確實在一定程度上影響著黃土高原NDVI的變化。但根據(jù)1998—2017年黃土高原地區(qū)NDVI和降水量的相關性分布圖可知,黃土高原大部分地區(qū)NDVI和降水之間不存在顯著的相關關系,部分面積NDVI和降水存在顯著正相關關系(P<0.05),主要集中在丘陵溝壑區(qū)、高原溝壑區(qū)北部和東部河谷及土石山區(qū)北部。這說明,降水確實在黃土高原地區(qū)的局部地區(qū)發(fā)揮了促進植被增加的作用。

圖6 黃土高原地區(qū)NDVI與降水相關性以及NDVI不同變化趨勢區(qū)域的逐年降水量Fig.6 Distribution of correlation between NDVI and precipitation in the Loess Plateau and annual precipitation in different NDVI trends zones

3.2 造成NDVI與降水之間突變點的原因

植被的生長確實需要充足的水分供應,但是許多研究表明黃土高原地區(qū)大部分區(qū)域年降水并未發(fā)生明顯的增加趨勢[42],造成黃土高原地區(qū)植被在1998—2017年間增加的主要驅動力應該是退耕還林、植樹造林等水土保持工程[21]。黃土高原地區(qū)1998—2017年多年平均年降水量為445 mm,空間變化范圍為109—878 mm,北部風沙區(qū)和農灌區(qū)、高原溝壑區(qū)、丘陵溝壑區(qū)和東部河谷及土石山區(qū)四大區(qū)域的多年平均年降水量分別為271、483 、459 、546mm,北部風沙區(qū)和農灌區(qū)的多年平均年降水量最小,東部河谷及土石山區(qū)的多年平均年降水量最大(圖7)。根據(jù)M-K趨勢檢驗分析,1998—2017年黃土高原地區(qū)年降水的變化趨勢如圖7所示。黃土高原地區(qū)年降水顯著增加的區(qū)域主要集中在丘陵溝壑區(qū)和東部河谷及土石山區(qū)中部(P<0.05),這可能是該區(qū)域NDVI顯著增加的一個重要原因,但是黃土高原地區(qū)大部分區(qū)域的降水并不存在顯著變化趨勢。圖7顯示了 1998—2017年間黃土高原地區(qū)NDVI均值的變化過程。根據(jù)M-K趨勢檢驗,黃土高原地區(qū)NDVI均值在1998—2017年呈現(xiàn)顯著增加趨勢。由Pettitt檢驗可知,黃土高原NDVI均值在2006年出現(xiàn)顯著突變點,這可能是造成NDVI和降水關系在2006年發(fā)生顯著突變點的主要原因(7),而1999年中國政府開始在黃土高原地區(qū)實施了退耕還林政策,極大的改善了黃土高原地區(qū)的生態(tài)環(huán)境,這可能是造成NDVI在2006年發(fā)生突變的可能原因,這同時表明黃土高原地區(qū)整體NDVI發(fā)生顯著突變相對于黃土高原地區(qū)實施退耕還林政策開始實施的時間滯后了7年,植被恢復需要一個較長的時間周期。

圖7 黃土高原地區(qū)多年平均年降水量空間分布；年降水量變化趨勢；年NDVI平均值變化；年NDVI平均值得突變點Fig.7 Spatial distribution of multi-year average annual precipitation in the Loess Plateau;trend of annual precipitation;change of average value of annual NDVI;change-point of average value of annual NDVI

4 結論

(1)在1998—2017年,黃土高原地區(qū)73.49%面積NDVI存在顯著變化趨勢,大面積區(qū)域植被在未來依舊呈現(xiàn)增加的趨勢,東部河谷及土石山區(qū)的植被覆蓋狀況最好,北部風沙區(qū)和農灌區(qū)的植被覆蓋狀況最差；

(2)1998—2017年間黃土高原地區(qū)丘陵溝壑區(qū)與高原溝壑區(qū)的NDVI增加幅度大于黃土高原地區(qū)整體的增加幅度,北部風沙區(qū)和農灌區(qū)以及黃土高原地區(qū)邊界區(qū)域的NDVI增加滯后于整體,NDVI逐漸從東南向西北增加；

(3)本研究構建了基于耦合協(xié)調度的NDVI降水關系突變點診斷方法,識別出黃土高原地區(qū)NDVI降水關系在2006年存在顯著突變點,高原溝壑區(qū)和北部風沙區(qū)和農灌區(qū)的NDVI變化對黃土高原地區(qū)整體的重要性較大,同時高原溝壑區(qū)和東部河谷及土石山區(qū)的降水變化在黃土高原地區(qū)占有重要地位。

(4)降水對黃土高原地區(qū)NDVI的增加具有積極的促進作用,在丘陵溝壑區(qū)、高原溝壑區(qū)北部和東部河谷及土石山區(qū)北部NDVI和降水存在顯著正相關關系,但是黃土高原地區(qū)大部分區(qū)域的降水并不存在顯著變化趨勢,因此造成黃土高原地區(qū)NDVI與降水關系在2006年發(fā)生顯著突變的主要原因應該是退耕還林等水土保持人類工程措施,黃土高原地區(qū)整體NDVI發(fā)生顯著突變相對于黃土高原地區(qū)實施退耕還林政策開始實施的時間滯后了7年,植被恢復需要一個較長的時間周期。