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

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

1 湖南師范大學(xué)資源與環(huán)境科學(xué)學(xué)院,長沙 410081 2 湖南師范大學(xué)地理空間大數(shù)據(jù)挖掘與應(yīng)用湖南省重點(diǎn)實(shí)驗(yàn)室,長沙 410081 3 中國科學(xué)院亞熱帶農(nóng)業(yè)生態(tài)研究所,長沙 410125

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

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

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

1 研究區(qū)概況

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

圖1 研究區(qū)域位置圖Fig.1 Location map of research area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源與處理

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

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

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

2.2 研究方法

2.2.1Mann-Kendall趨勢檢驗(yàn)

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

(1)

其中：

(2)

(3)

(4)

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

2.2.2MEMD

MEMD是對標(biāo)準(zhǔn)經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical mode decomposition,EMD)的多元擴(kuò)展,是通過多維空間超球面的方向向量求均值包絡(luò),從而實(shí)現(xiàn)不同頻率尺度下不同變量的共同模式分解[28,34]。自然界中的各種事物的整體變化受到不同過程影響,而且在不同尺度下以不同的強(qiáng)度發(fā)生[35]。在同一個(gè)尺度下發(fā)生的過程可以分解成為相同的本征模函數(shù)(Intrinsic mode function,IMF)。每個(gè)變量分解出的本征模函數(shù)(IMF)數(shù)量相同,每層的IMF頻率相同,IMF頻率高低由分解的次序所決定,先分解的頻率高,后分解的頻率低。不同變量對應(yīng)的IMF按尺度對齊,組成多元IMF。

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

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

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

(5)

(6)通過D(s)=V(s)-M(s)得到本征模函數(shù)D(s)。假如D(s)滿足多元IMF的判斷標(biāo)準(zhǔn),則將D(s)=V(s)-M(s)的結(jié)果作為步驟(2)的輸入變量。如果不滿足,則將D(s)作為步驟(2)的輸入變量,繼續(xù)進(jìn)行迭代運(yùn)算。

3 結(jié)果與分析

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

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

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

NDVI: Normalized Difference Vegetation Index

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

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

3.2 徑流深及其影響因子的多尺度關(guān)系

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

圖3 基于多元經(jīng)驗(yàn)?zāi)B(tài)分解的徑流深及其影響因子的本征模函數(shù)和殘差Fig.3 Intrinsic mode function (IMF) and residues of runoff and its influencing factors based on multivariate empirical mode decomposition

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

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

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

各流域徑流深及其影響因子的IMF所表征的尺度方差貢獻(xiàn)率不同(表3)。郁江流域徑流深方差貢獻(xiàn)率主要分布在IMF1和IMF2處(32%和36%),所對應(yīng)尺度為3年和5年。紅水河流域徑流深方差貢獻(xiàn)率主要分布在IMF3和IMF4處(26%和28%),對應(yīng)的尺度約為10年和22年。潯江和梧州流域徑流深方差貢獻(xiàn)率主要分布在IMF1和IMF4處(35%和31%、36%和34%),對應(yīng)尺度約為3年和22年。從其余影響因子的IMF方差貢獻(xiàn)率可以看出,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%)。郁江和紅水河流域的潛在蒸散發(fā)的表征尺度主要分布在IMF1和IMF2處,對應(yīng)的特征尺度約為3年和5年；潯江和梧州流域主要分布在IMF1和IMF4處,對應(yīng)的特征尺度約為3年和22年。4個(gè)流域的氣溫表征尺度主要在IMF1和IMF2處,即主要分布在年際尺度上。

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

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

圖4 基于多元經(jīng)驗(yàn)?zāi)B(tài)分解的各本征模函數(shù)(IMF)和殘差的徑流深與其影響因子的相關(guān)性分析圖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：歸一化植被指數(shù)Normalized difference vegetation index；PRE：降水 Precipitation；PT：潛在蒸散發(fā) Potential evaporation；T：氣溫 Temperature；**表示顯著性水平小于0.01；*表示顯著性水平小于0.05

此外,利用MEMD方法進(jìn)行數(shù)據(jù)分解有兩點(diǎn)值得注意。第一,一些變量的殘差部分的方差貢獻(xiàn)率較大,比如郁江流域的NDVI,殘差部分的方差貢獻(xiàn)率為32.1%,這可能是因?yàn)檠芯繒r(shí)段較短,更大的表征尺度未被發(fā)現(xiàn)；第二,每個(gè)變量的方差貢獻(xiàn)率相加之和不等于100%,主要是因?yàn)槊總€(gè)特征尺度下的生態(tài)過程不是獨(dú)立的,變量的不同IMF缺少完全正交性[36]。

3.3 基于MEMD的徑流深預(yù)測

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

表4 基于多元模態(tài)分解的每個(gè)本征模函數(shù)和殘差的徑流深多元逐步回歸預(yù)測模型及回歸統(tǒng)計(jì)特征(可調(diào)整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原始數(shù)據(jù)-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原始數(shù)據(jù)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原始數(shù)據(jù)-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原始數(shù)據(jù)-287.998+0.711(0.887)PRE0.78117.66

NDVI：歸一化植被指數(shù)Normalized difference vegetation index；PRE：降水 Precipitation；PT：潛在蒸散發(fā)Potential evaporation；T：氣溫 Temperature；括號中的數(shù)字代表標(biāo)準(zhǔn)偏回歸系數(shù)

表5 多元模態(tài)分解與多元逐步回歸預(yù)測徑流深效果

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

圖5 多元經(jīng)驗(yàn)?zāi)B(tài)分解前后多元逐步回歸預(yù)測徑流與實(shí)際徑流對比Fig.5 Observed vs. estimated runoff using stepwise multiple linear regression models before and after multivariate empirical mode decomposition

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

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

4 討論

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

徑流深與降水和潛在蒸散發(fā)在不同的時(shí)間尺度始終呈顯著相關(guān)(P<0.05),而與NDVI和氣溫在某些尺度上并不顯著。說明與降水和潛在蒸散相比,氣溫和植被對于徑流的影響機(jī)制相對復(fù)雜。在短時(shí)間尺度內(nèi),氣溫主要通過影響蒸散發(fā)等因素間接影響徑流的變化。對于降水補(bǔ)給型河流來說,徑流主要受到降水量的影響,氣溫的影響較小。但在較長的時(shí)間尺度上,氣溫變化會對氣候系統(tǒng)產(chǎn)生長期影響,進(jìn)而通過影響降水的時(shí)空分配格局引起流域徑流變化。植被對于徑流影響機(jī)制存在復(fù)雜性,主要原因是植被一方面可以減少表面徑流,提高蒸散發(fā)；另一方面,氣溫和降水可以通過影響植被覆蓋間接影響徑流[39]?？λ固亓饔蛴捎趶?fù)雜的雙層水文地質(zhì)結(jié)構(gòu),可能更加劇了上述影響機(jī)制的復(fù)雜性。

從預(yù)測的效果來看,經(jīng)過多元模態(tài)分解后的數(shù)據(jù)進(jìn)行徑流預(yù)測精度高于原始數(shù)據(jù)直接預(yù)測徑流深,主要原因在于多元經(jīng)驗(yàn)?zāi)B(tài)方法在某些特征尺度上可以獲得相關(guān)因子對徑流的影響程度。該方法在土壤水分含量[36,40]預(yù)測上取得較好的效果。由于多元經(jīng)驗(yàn)?zāi)B(tài)分解的理論基礎(chǔ)建立時(shí)間相對較短,所以在徑流領(lǐng)域的應(yīng)用并不多見。

5 結(jié)論

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

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

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

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

(4)利用MEMD方法對徑流深進(jìn)行預(yù)測,發(fā)現(xiàn)4個(gè)流域徑流深的預(yù)測精度均優(yōu)于直接利用多元逐步回歸擬合的結(jié)果,因此該方法可用于西南喀斯特區(qū)徑流量的預(yù)測。