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

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

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

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

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

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

1 材料與方法

1.1 研究區概況

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

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

1.2 數據來源

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

1.3.1Mann-Kendall趨勢檢驗與Hurst指數

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

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

1.3.2RDR指數

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

(1)

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

1.3.3重心轉移模型

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

(2)

(3)

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

1.3.4耦合協調度理論

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

(4)

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

(5)

(6)

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

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

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

(7)

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

(8)

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

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

D=(C×T)1/2

(9)

T=a×u1+b×u2

(10)

式中,D為耦合協調度；C為耦合度；T為NDVI子系統與降水子系統的綜合調和指數,它反映NDVI子系統與降水子系統的整體協同效應或貢獻；a、b為待定系數,實際中常常認為兩個子系統的重要性相同,所以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時認為數據中存在均值變異點。

2 結果與分析

2.1 NDVI變化特征

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

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

2.2 NDVI時空變化差異

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

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

圖3 1998—2017年黃土高原地區NDVI的RDR指數和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子系統與的降水子系統權重分布變化

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

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

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

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

圖5 NDVI與降水的耦合協調度變化及其突變點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年黃土高原地區NDVI顯著增加區域(P<0.05)和不顯著變化區域(P>0.05)逐年平均年降水量的時間變化過程。NDVI顯著增加區域的降水在1998—2017年間均高于NDVI變化不顯著的區域,這說明降水確實在一定程度上影響著黃土高原NDVI的變化。但根據1998—2017年黃土高原地區NDVI和降水量的相關性分布圖可知,黃土高原大部分地區NDVI和降水之間不存在顯著的相關關系,部分面積NDVI和降水存在顯著正相關關系(P<0.05),主要集中在丘陵溝壑區、高原溝壑區北部和東部河谷及土石山區北部。這說明,降水確實在黃土高原地區的局部地區發揮了促進植被增加的作用。

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

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

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

圖7 黃土高原地區多年平均年降水量空間分布；年降水量變化趨勢；年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年,黃土高原地區73.49%面積NDVI存在顯著變化趨勢,大面積區域植被在未來依舊呈現增加的趨勢,東部河谷及土石山區的植被覆蓋狀況最好,北部風沙區和農灌區的植被覆蓋狀況最差；

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

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

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