喀斯特流域徑流對植被和氣候變化的多尺度響應

朱靖軒,劉 雯,*,李振煒,李笑含,張志慧,徐憲立

1 湖南師范大學資源與環境科學學院,長沙 410081 2 湖南師范大學地理空間大數據挖掘與應用湖南省重點實驗室,長沙 410081 3 中國科學院亞熱帶農業生態研究所,長沙 410125

徑流是人類生產生活用水的主要來源之一,對徑流的研究是對水資源進行有效開發利用和管理的基礎[1]。過去幾十年中,全球氣候變化促使水循環加快,改變了降水的時空變化,極端水文事件發生更加頻繁[2- 3],導致我國南方地區的洪澇災害在不斷增加,水生態環境在不斷惡化[4- 6]。喀斯特地區在全球分布十分廣泛,約占地球陸地面積的7%—12%,供應了全球約25%人口的水源[7]。我國是全球喀斯特面積最大、分布最為廣泛的國家之一,尤其是在西南地區,喀斯特連片分布面積多達54萬km2[7- 8]。由于碳酸鹽巖的廣泛分布,喀斯特地貌具有土層淺薄、土壤入滲能力強和地形復雜等特點[9- 10]。因此,喀斯特地表持水能力較弱,水資源有效利用較為困難；同時,喀斯特地區局部地下管網通暢性差,極易發生洪澇災害。喀斯特地區旱澇頻發,嚴重限制了當地社會經濟可持續發展。因此,研究喀斯特流域徑流變化特征,以及其對氣候和植被變化的響應,對合理配置喀斯特流域水資源,促進流域社會經濟和生態環境保護可持續發展,具有十分重要的現實意義。

徑流的變化規律錯綜復雜,其具有非線性、隨機性和突變性等復雜特點,主導了自然界中整個水文系統的變化[11]。其受到了氣候、植被、土壤、地形和人類活動的共同影響,其中植被和氣候變化是引起徑流變化的重要因素[12]。植被通過地上部分的林冠截留對降水進行再分配,以及通過地表枯落物和地下根系增加土壤的入滲能力等改變徑流的產生機制[13- 14]。歸一化植被指數(Normalized Difference Vegetation Index,NDVI)與植被覆蓋度、生長狀況、生物量和光合作用強度有著密切的關系,被廣泛應用于植被動態監測[15]。已有研究表明,在年內變化上,NDVI與徑流呈顯著正相關,但年際變化上,關系相對復雜,不具有明顯相關性[16]。氣候因素包括降水、潛在蒸散發以及氣溫等,降水是引起徑流變化最敏感的因子[17]。氣候變化直接導致降水的空間分布、時間變異以及地表蒸發等的變化來影響流域的徑流量[18- 19]。流域徑流與氣象要素相關性研究表明,徑流下降的主要原因是人類活動引起了流域內蒸發量增加,而且導致徑流對降水的響應程度下降[20]。很多學者仍采用線性趨勢法、Man-Kendall趨勢檢驗、滑動平均和多元回歸分析等方法研究氣候水文變化趨勢[21- 24]。然而,氣候和水文系統是一個復雜的非線性系統[25],這些因子的長期變化過程大多數是非線性、非平穩的變化過程,而且通常伴隨著多尺度或者周期性振蕩。所以,線性方法并不能準確地揭示出植被、氣候和水文系統的非線性特征。因此,一些研究者采用了小波分析和集合經驗模態分解等非線性方法。柏玲等[6]利用集合經驗模態分解方法對開都河長時間徑流序列進行了多尺度分析；董喆等[26]利用小波分析方法研究了在不同周期下,降水與植被對徑流的調節作用；吳創收等[27]通過小波分解方法,揭示了珠江流域徑流周期變化特征,同時得出徑流的變化周期與降水的變化周期具有較好的相關性。但是,集合經驗模態的參數設置沒有特定的方法和規則,針對不同的信號缺乏自適應性,且應用于多元數據分析時存在分解的本征模函數的數量不一致問題。小波分析雖然具有時頻分辨功能,但是其作用的對象是單一序列的,只能反映不同時頻域上的徑流變化規律,無法反映植被和氣候因素對徑流的影響[20]。同時,小波分解的小波基和分解層的數量需要人為的設定,自適應性較差[6],使得其在水文領域的應用具有一定的局限性。

隨著信號檢測手段的不斷進步,有關于信號處理的方法在不斷出現。多元經驗模態分解(multivariate empirical mode decomposition,MEMD)是一種與小波分析相對應的另外一種多尺度分析方法,是經驗模態分解(Empirical Mode Decomposition,EMD)在多元空間上的擴展[28]。它可以依據數據自身的尺度特征,經驗性地將多元數據分解成多個表征尺度,有效解決了線性和平穩性假設的問題[28]。EMD和EEMD在處理多元數據時,存在分解的本征模函數的數量不一致問題,而MEMD可以將多元數據分解成具有相同個數的多元本征模函數分量。雖然 MEMD方法具有很多明顯的優點,但是在地質結構背景復雜的喀斯特地區徑流研究中尚未得到應用。因此,本文利用MEMD方法,基于1982—2015年植被、氣象和水文數據,分析喀斯特流域徑流及其影響因子的非線性變化特征,研究徑流在不同時間尺度對這些因子的響應,進而揭示徑流對植被、氣象因子變化的響應過程,并對徑流深進行預測,為喀斯特流域復雜地質背景下徑流相關研究提供幫助。

1 研究區概況

西江是華南地區最長河流,發源于云南省曲靖市烏蒙山余脈馬雄山的東麓,先后流經云南、貴州、廣西和廣東等省份(圖1)。西江全長2214 km,流域(21°36′17″—27°00′21″N, 102°16′52″—113°23′51″E)總面積約為3.53×105km2,地形以高原、山地和平原為主,最低和最高海拔為2 m和2888 m,西部海拔高于東部。氣候為典型的亞熱帶和熱帶季風濕潤氣候,多年平均氣溫約為18 ℃,多年平均降水量約為1470 mm,主要集中在每年的5月至10月份,流域的降水由東南向西北逐漸減少。土地利用以林地為主,約占流域面積的62%。流域內喀斯特地貌分布較為廣泛,約占流域面積的45%。研究區內選擇郁江、紅水河、潯江和梧州4個典型喀斯特流域(圖1),其喀斯特地貌面積占比分別為42%、64%、51%和47%。

圖1 研究區域位置圖Fig.1 Location map of research area

2 數據與方法

2.1 數據來源與處理

(1)徑流數據：1982—2015年西江流域4個水文站(南寧、遷江、大湟江口、梧州)的年徑流數據,來自水利部編制的水文年鑒。

(2)NDVI數據：采用1982—2015年GIMMS NDVI數據,該數據是由美國國家航天航空局(NASA)推出的全球植被變化數據。時間分辨率為15d,空間分辨率為8 km×8 km[29]。利用最大值合成法(Maximum Value Composite,MVC)生成年NDVI最大值,可以有效地減少因為云、陰影、氣溶膠、太陽高度角及視角等帶來的影響[30]。

(3)1982—2015年西江流域36個氣象站的降水、氣溫等數據。來自中國氣象科學數據共享服務網。通過空間插值方法獲得流域面氣象數據。此外,1982—2015年的潛在蒸散發數據是根據研究區內氣象站點的觀測資料數據,利用聯合國糧農組織修正的Penman-Monteith公式[31]計算得到。

2.2 研究方法

2.2.1Mann-Kendall趨勢檢驗

Mann-Kendall 趨勢檢驗是研究水文變化趨勢研究的良好方法[32- 33],其主要的原理如下：

(1)

其中：

(2)

(3)

(4)

式中,Xi代表第i年的數值,Xj代表第j年的數值。當iZ(1-α/2),則拒絕無趨勢的原假設,即在α置信水平上,時間序列的數據有明顯的上升或者下降趨勢。如果Zs為正值,則表示為上升趨勢,反之,則表示為下降趨勢。若|Zs|>1.96,則表示通過了置信度為95%的顯著性趨勢檢驗。

2.2.2MEMD

MEMD是對標準經驗模態分解(Empirical mode decomposition,EMD)的多元擴展,是通過多維空間超球面的方向向量求均值包絡,從而實現不同頻率尺度下不同變量的共同模式分解[28,34]。自然界中的各種事物的整體變化受到不同過程影響,而且在不同尺度下以不同的強度發生[35]。在同一個尺度下發生的過程可以分解成為相同的本征模函數(Intrinsic mode function,IMF)。每個變量分解出的本征模函數(IMF)數量相同,每層的IMF頻率相同,IMF頻率高低由分解的次序所決定,先分解的頻率高,后分解的頻率低。不同變量對應的IMF按尺度對齊,組成多元IMF。

(1)獲取n維空間向量集X；

(2)計算時間序列V(s)沿著給定方向Xθk的投影pθk(s)；

(5)利用以下公式計算包絡線的均值M(s)；

(5)

(6)通過D(s)=V(s)-M(s)得到本征模函數D(s)。假如D(s)滿足多元IMF的判斷標準,則將D(s)=V(s)-M(s)的結果作為步驟(2)的輸入變量。如果不滿足,則將D(s)作為步驟(2)的輸入變量,繼續進行迭代運算。

3 結果與分析

3.1 徑流深及其影響因子的變化趨勢

4個流域中,年均徑流深最小的為紅水河流域,最大的為梧州流域(表1)。4個流域間年均NDVI、降水量、潛在蒸散發和氣溫之間差異均較小。流域的徑流深和降水量的變異系數均呈中等程度變異(10% < Cv < 100%),而NDVI、潛在蒸散發和氣溫均呈弱變異(Cv < 10%)。在4個流域中,徑流深與降雨量呈顯著正相關(P<0.01),與潛在蒸散發呈顯著負相關(P<0.01),而與NDVI和氣溫無顯著相關性(P>0.05)(圖2)。

表1 徑流深及其影響因子平均值及變化趨勢

NDVI: Normalized Difference Vegetation Index

圖2 徑流深及其影響因子相關性矩陣圖Fig.2 Correlative Matrix plot of runoff and its influencing factorsRF：徑流深 Runoff；NDVI：歸一化植被指數Normalized difference vegetation index；PRE：降水 Precipitation；PT：潛在蒸散發 Potential evaporation； T：氣溫 Temperature； R：相關系數Correlation coefficient；P：顯著性水平 Significance level

通過Mann-Kendall趨勢檢驗,發現流域的NDVI和氣溫均表現為顯著增加的趨勢(P<0.05；表1)。流域的徑流深和降水均表現為下降的趨勢,但只有紅水河流域的徑流深通過了0.05信度檢驗。流域的潛在蒸散發均表現為增加趨勢,但僅郁江流域通過了0.05信度檢驗,呈顯著增加趨勢,其余流域增加趨勢不顯著。

3.2 徑流深及其影響因子的多尺度關系

利用MEMD方法將4個流域的徑流深及其影響因子的多元數據分解成4個本征模函數IMF和1個殘差(圖3)。與EMD方法相似,各變量的震蕩周期隨著IMF增大而增大。對于一組經過分解的多元數據,在相同IMF中,不同的變量通常具有相同的震蕩數量,所以具有相似的震蕩頻率[36]。本文使用Hilbert 變換來計算不同IMF之間的瞬時頻率,得到不同的尺度。

圖3 基于多元經驗模態分解的徑流深及其影響因子的本征模函數和殘差Fig.3 Intrinsic mode function (IMF) and residues of runoff and its influencing factors based on multivariate empirical mode decomposition

每個流域不同因子的IMF對應的時間尺度的平均值代表該IMF對應的特征尺度(表2)。郁江流域IMF1、IMF2、IMF3和IMF4對應的特征尺度為2.9、5.0、8.7、19.5年；紅水河流域IMF1、IMF2、IMF3和IMF4對應的特征尺度為3.1、5.0、9.6、22.4年；潯江流域IMF1、IMF2、IMF3和IMF4對應的特征尺度為3.0、5.6、10.0、21.5年；梧州流域IMF1、IMF2、IMF3和IMF4對應的特征尺度為3.0、5.4、10.0、22.4年。除郁江流域外,其他3個流域各因子的IMF4所對應尺度的變異系數較大,表明該尺度下無可以代表徑流深與其影響因子的共同尺度。4個流域中,4個本征模函數IMF對應的平均時間尺度比較接近,而且徑流深的尺度也比較相似,該結果表明控制徑流的過程發生在相似的尺度。

表2 徑流深及其影響因子本征模函數(IMF)的尺度

IMF：本征模函數 Intrinsic mode function；括號中的數字代表徑流深及其影響因子本征模函數所對應尺度的變異系數(%)

各流域徑流深及其影響因子的IMF所表征的尺度方差貢獻率不同(表3)。郁江流域徑流深方差貢獻率主要分布在IMF1和IMF2處(32%和36%),所對應尺度為3年和5年。紅水河流域徑流深方差貢獻率主要分布在IMF3和IMF4處(26%和28%),對應的尺度約為10年和22年。潯江和梧州流域徑流深方差貢獻率主要分布在IMF1和IMF4處(35%和31%、36%和34%),對應尺度約為3年和22年。從其余影響因子的IMF方差貢獻率可以看出,NDVI除了紅水河流域主要分布在IMF1和IMF4處,且其殘差值較小(2.4%),其余流域主要分布在IMF1和IMF4及殘差部分。郁江和紅水河流域降水的主要表征尺度分別分布在IMF1和 IMF2、IMF1和IMF3處(43%和21%、51%和16%),潯江、梧州流域主要分布在IMF1和IMF4處(50.1%和18.6%、50.6%和20.4%)。郁江和紅水河流域的潛在蒸散發的表征尺度主要分布在IMF1和IMF2處,對應的特征尺度約為3年和5年；潯江和梧州流域主要分布在IMF1和IMF4處,對應的特征尺度約為3年和22年。4個流域的氣溫表征尺度主要在IMF1和IMF2處,即主要分布在年際尺度上。

徑流深及其影響因子多尺度相關性如圖4所示。每個流域降水量和潛在蒸散發與徑流深的多尺度關系最顯著,在所有表征尺度均呈顯著相關。徑流深與NDVI和氣溫在一個表征尺度上呈顯著正或負相關,而在另外一個尺度上呈顯著負或正相關。這主要是因為在一個表征尺度上呈現正或負相關,會降低在另外一個表征尺度上的負或正相關程度,同時導致了徑流深與NDVI和氣溫在年尺度上相關性不顯著。在郁江流域,NDVI與徑流的相關性主要表現在IMF1、IMF2、IMF4和殘差處；紅水河流域為IMF1、IMF3、IMF4和殘差；潯江流域為IMF2、IMF3、IMF4和殘差；梧州流域為IMF2、IMF4和殘差。郁江流域的氣溫與徑流深的相關性主要體現在IMF3、IMF4和殘差處,其余流域主要體現在IMF3和殘差處。一般情況下,隨著表征尺度的增大,徑流與其影響因子之間的相關性有增大的趨勢,即這些因子對徑流深的影響程度有增大趨勢。除此之外,在殘差部分,徑流深與其影響因子均存在顯著相關性(P<0.05),這在一定程度上表明這些影響因子與徑流深之間存在一定關聯。徑流深與NDVI和氣溫在年尺度上不存在顯著相關性,但是徑流深與兩者通過MEMD分解,在某些表征尺度上呈顯著相關,這說明了研究徑流深與其影響因子的時間序列多尺度關系可以獲得更多的信息。

表3 徑流深及其影響因子的本征模函數(IMF)和殘差占原始數據方差的百分比

圖4 基于多元經驗模態分解的各本征模函數(IMF)和殘差的徑流深與其影響因子的相關性分析圖Table 4 Figure of correlation coefficients between runoff and its influencing factors for each Intrinsic mode function (IMF) and residue based on multivariate empirical mode decompositionNDVI：歸一化植被指數Normalized difference vegetation index；PRE：降水 Precipitation；PT：潛在蒸散發 Potential evaporation；T：氣溫 Temperature；**表示顯著性水平小于0.01；*表示顯著性水平小于0.05

此外,利用MEMD方法進行數據分解有兩點值得注意。第一,一些變量的殘差部分的方差貢獻率較大,比如郁江流域的NDVI,殘差部分的方差貢獻率為32.1%,這可能是因為研究時段較短,更大的表征尺度未被發現；第二,每個變量的方差貢獻率相加之和不等于100%,主要是因為每個特征尺度下的生態過程不是獨立的,變量的不同IMF缺少完全正交性[36]。

3.3 基于MEMD的徑流深預測

利用多元回歸模型可以得到每個IMF徑流深的預測模型(表4)。從表中可以看出,這些預測模型的可調整R2的范圍在0.67到1之間,而且一般隨著IMF增大而增大。可以確定,表征尺度越大,對徑流深的預測精度越高,并且模型中所選因子對徑流深的影響程度更大。基于所有IMF和殘差的徑流深預測結果,利用多元逐步回歸方法進行徑流深預測。4個流域徑流深預測值和實測值的R2,RMSE,NSE和1∶1線均優于直接利用多元逐步回歸擬合的結果(表5；圖5)。由此可見,利用MEMD方法對徑流深的預測精度要高于利用原始數據直接進行多元逐步回歸預測的結果。這進一步說明了在原始時間序列上的簡單分析難以解釋NDVI、溫度與徑流之間關系的復雜性。

表4 基于多元模態分解的每個本征模函數和殘差的徑流深多元逐步回歸預測模型及回歸統計特征(可調整R2和F值)

Table 4 Predictive modeling and regression statistics (F-value and adjustedR2) for runoff for each intrinsic mode function and residue using stepwise multiple linear regression based on multivariate empirical mode decomposition

流域WatershedIMF公式FunctionR2F郁江IMF1-3.587+0.581(0.801)PRE-74.669(-0.235)T-13.443(-0.201)NDVI0.7941.79IMF2-6.785+0.63(0.575)PRE-2.888(-0.816)PT+159.869(0.581)T0.8773.34IMF33.104+0.601(0.57)PRE-302.788(-0.471)T0.88122.04IMF4-2.091+0.386(0.871)PRE+247.387(0.661)T+13.458(0.513)ND-VI0.92126.74殘差951.909-1.309(-1.345)PT+67.625(0.65)T-4.882(-0.305)ND-VI1.0756358原始數據-363.473+0.702(0.836)PRE0.6974.23紅水河IMF1-0.13+0.289(0.614)PRE-1.063(-0.398)PT+40.846(0.219)T0.8251.58IMF20.142+0.604(0.815)PRE-22.39(-0.209)T0.7448.34IMF31.4+0.669(0.766)PRE-64.934(-0.242)T0.6735.04IMF4-0.398+0.593(0.568)PRE-1.589(-0.491)PT+85.17(0.187)T1.02393.07殘差2228.644-1.914(-1.35)PT+3.889(0.35)NDVI1.02927483原始數據575.359+0.509(0.726)PRE-0.653(-0.232)PT0.7345.00潯江IMF10.573+0.543(0.847)PRE0.7180.93IMF2-0.219+0.676(0.783)PRE+30.041(0.249)NDVI0.8593.73IMF32.009+0.632(0.55)PRE-1.711(-0.376)PT+17.616(0.245)NDVI0.90103.72IMF4-2.934+0.993(0.992)PRE0.981880.85殘差-2140.752+1.479(6.295)PRE+0.813(5.441)PT1.0106793原始數據-273.478+0.688(0.885)PRE0.78115.66梧州IMF1-0.717+0.574(0.845)PRE0.7179.79IMF2-1.885+0.624(0.719)PRE+42.341(0.3)NDVI+40.873(0.192)T0.8044.10IMF30.871+1.05(0.874)PRE-0.804(-0.146)PT0.94257.17IMF4-0.92+1.09(1.051)PRE+0.19(0.061)PT1.06412殘差-1320.21+32.976(0.92)T+1.052(0.112)P1.0613063296原始數據-287.998+0.711(0.887)PRE0.78117.66

NDVI：歸一化植被指數Normalized difference vegetation index；PRE：降水 Precipitation；PT：潛在蒸散發Potential evaporation；T：氣溫 Temperature；括號中的數字代表標準偏回歸系數

表5 多元模態分解與多元逐步回歸預測徑流深效果

R2：決定系數 Coefficient of determination；RMSE：均方根誤差 Root mean squared error；NSE：納什效率系數 Nash-Sutcliffe efficiency coefficient

圖5 多元經驗模態分解前后多元逐步回歸預測徑流與實際徑流對比Fig.5 Observed vs. estimated runoff using stepwise multiple linear regression models before and after multivariate empirical mode decomposition

圖6 徑流深預測值與單個本征模函數或殘差處預測值的決定系數Fig.6 Coefficient of determination between runoff predicted value and predicted IMFs (residue)

徑流深的預測值與每個IMF或者殘差處的預測值的決定系數如圖6所示。從該圖中可以看出每個IMF和殘差對于徑流深預測值的重要性程度。郁江流域主要貢獻者為IMF2；紅水河、潯江和梧州流域主要貢獻者為IMF4。即郁江流域的年際尺度上對徑流預測起主要作用,而紅水河、潯江和梧州流域的年代際尺度上對徑流預測起主要作用。

4 討論

多元經驗模態分解量化了流域徑流深與植被和氣候因子在不同尺度上的相關性,可以從復雜的角度初步理解喀斯特流域徑流與其影響因子的內在作用機制,加深對此區域徑流變化以及驅動因子的研究。運用傳統的線性分析方法,無法揭示徑流與其影響因子在不同尺度上的相關性,所以,以往的研究中,有些學者運用了小波分析和集合經驗模態分解(EEMD)等方法。李寶富等[37]運用了EEMD方法研究了沂河流域徑流對氣候變化的響應,得出了在不同尺度上氣候因子對于徑流的影響機制存在明顯不同,在較大尺度上往往與徑流的關系更加密切。凌紅波等[38]通過小波分析方法,揭示了天山地區徑流和氣象因子之間存在明顯的多時間尺度相關關系。這些研究與本文研究結果較為一致,為區域水資源合理利用提供了重要的理論支撐。

徑流深與降水和潛在蒸散發在不同的時間尺度始終呈顯著相關(P<0.05),而與NDVI和氣溫在某些尺度上并不顯著。說明與降水和潛在蒸散相比,氣溫和植被對于徑流的影響機制相對復雜。在短時間尺度內,氣溫主要通過影響蒸散發等因素間接影響徑流的變化。對于降水補給型河流來說,徑流主要受到降水量的影響,氣溫的影響較小。但在較長的時間尺度上,氣溫變化會對氣候系統產生長期影響,進而通過影響降水的時空分配格局引起流域徑流變化。植被對于徑流影響機制存在復雜性,主要原因是植被一方面可以減少表面徑流,提高蒸散發；另一方面,氣溫和降水可以通過影響植被覆蓋間接影響徑流[39]。喀斯特流域由于復雜的雙層水文地質結構,可能更加劇了上述影響機制的復雜性。

從預測的效果來看,經過多元模態分解后的數據進行徑流預測精度高于原始數據直接預測徑流深,主要原因在于多元經驗模態方法在某些特征尺度上可以獲得相關因子對徑流的影響程度。該方法在土壤水分含量[36,40]預測上取得較好的效果。由于多元經驗模態分解的理論基礎建立時間相對較短,所以在徑流領域的應用并不多見。

5 結論

本文主要分析了4個喀斯特流域徑流及其影響因子1982—2015年的變化趨勢,以及徑流對這些因子變化的多尺度響應。主要結論如下：

(1)通過Mann-Kendall趨勢檢驗方法,表明4個流域的氣溫和NDVI均呈顯著增加的趨勢；除紅水河流域的徑流深呈顯著下降趨勢外,其余流域的下降趨勢并不顯著；同時4個流域降水的下降趨勢也不顯著。

(2)利用MEMD將4個流域的徑流深及其影響因子的多元數據分解成4個本征模函數IMF和1個殘差,發現郁江流域徑流深方差貢獻率主要分布在3年和5年尺度上；而紅水河流域徑流深方差貢獻率主要分布在10年和22年尺度上；潯江和梧州流域徑流深方差貢獻率主要分布在3年和22年尺度上。

(3)徑流深與其影響因子在不同時間尺度上的相關性與原始數據相關性不同,隨著尺度的不同而存在差異性。在不同尺度下,4個流域降水量和潛在蒸散發均與徑流深之間存在顯著相關性,而與NDVI和氣溫在原始時間序列上不存在顯著相關性,但是其與兩者的時間多尺度關系在某些表征尺度上存在顯著相關性。

(4)利用MEMD方法對徑流深進行預測,發現4個流域徑流深的預測精度均優于直接利用多元逐步回歸擬合的結果,因此該方法可用于西南喀斯特區徑流量的預測。