基于耦合協(xié)調(diào)度的大理河流域徑流和輸沙關(guān)系分析

賈 路，任宗萍※，李占斌，李 鵬，徐國(guó)策，張鐵鋼，楊媛媛

（1. 西安理工大學(xué)省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710048；2. 旱區(qū)生態(tài)水文與災(zāi)害防治國(guó)家林業(yè)和草原局重點(diǎn)實(shí)驗(yàn)室，西安 710048；3. 水利部牧區(qū)水利科學(xué)研究所，呼和浩特010020）

0 引 言

黃河是中華民族的母親河，隨著氣候變化與人來(lái)活動(dòng)的影響，黃河流域生態(tài)環(huán)境狀況及自然條件發(fā)生了巨大變化，黃河流域的健康和高質(zhì)量成為一項(xiàng)重大國(guó)家戰(zhàn)略，黃河流域的發(fā)展關(guān)系著流域內(nèi)人民的生產(chǎn)和生活的長(zhǎng)遠(yuǎn)發(fā)展[1-2]。黃河中游地區(qū)占據(jù)了黃河流域的大部分面積，該區(qū)域位于世界水土流失最為嚴(yán)重的地區(qū)之一-黃土高原，因此治理黃河流域的水沙變化問(wèn)題一直是治黃工作的核心問(wèn)題[3-4]。自20 世紀(jì)50 年代開始，中國(guó)在黃土高原地區(qū)實(shí)施了大規(guī)模的水土保持工程建設(shè)，特別是淤地壩工程的建設(shè)，有效地?cái)r截了黃土高原地區(qū)產(chǎn)生的泥沙，調(diào)節(jié)了徑流，阻止泥沙進(jìn)入黃河，黃河徑流和泥沙顯著減少[5]。同時(shí)，近幾十年全球氣候發(fā)生顯著變化，全球變暖的趨勢(shì)日益加劇，極端事件頻發(fā)，已經(jīng)受到國(guó)際社會(huì)的普遍關(guān)注，氣候變化進(jìn)一步加劇黃河流域水資源的供需矛盾[6-7]。氣候變化有可能會(huì)加速全球水文循環(huán)，這對(duì)水沙關(guān)系的變化將產(chǎn)生非常重大的影響。研究氣候變化條件和人類活動(dòng)作用下水文過(guò)程和水沙變化的狀況是目前水文學(xué)和水土保持的熱點(diǎn)問(wèn)題之一。

徑流和泥沙作為流域極具重要價(jià)值的資源，在流域社會(huì)、經(jīng)濟(jì)和生態(tài)環(huán)節(jié)中扮演著重要的角色。許多研究調(diào)查了氣候變化和人類活動(dòng)影響下的徑流和輸沙的變化特征。曹麗娟等[8]對(duì)未來(lái)氣候變化條件下黃河和長(zhǎng)江流域極端徑流影響預(yù)估進(jìn)行了研究。歐陽(yáng)潮波等[9]對(duì)60 a 來(lái)黃河河龍區(qū)間水沙變化特征及人類活動(dòng)的影響進(jìn)行了評(píng)價(jià)。高海東等分析了2000—2017 年河龍區(qū)間輸沙量銳減的原因[10]。盡管關(guān)于徑流和輸沙對(duì)氣候變化和人類活動(dòng)時(shí)空變化響應(yīng)及歸因分析的研究已經(jīng)被進(jìn)行了研究，并取得了很多的研究成果。但是有一個(gè)事實(shí)是至關(guān)重要的，由于過(guò)去的研究大多數(shù)是基于單變量角度對(duì)流域水文變量間的關(guān)系包括水沙關(guān)系進(jìn)行的研究，所以這成為造成識(shí)別的變量之間的變化存在巨大的差異[11-12]。在自然界中，地表水、大氣、植被、土壤和地下水等環(huán)境要素之間存在復(fù)雜的相互作用[13]，使用單變量的分析方法可以讓事情變得簡(jiǎn)單、便捷。然而，單變量的分析方法只關(guān)注了單個(gè)變量的變化（例如徑流或輸沙），這與現(xiàn)實(shí)差別較大，不能完全揭示對(duì)氣候變化和人類活動(dòng)的水沙響應(yīng)。因此，有必要基于多變量的視角研究徑流和輸沙的變化特征。

大理河流域作為黃河的二級(jí)支流[14]，位于中國(guó)典型的干旱和半干旱區(qū)域。同時(shí)，由于大理河流域地處中國(guó)黃土高原地區(qū)，地形破碎，溝壑區(qū)密度大，水土流失比較嚴(yán)重，耕地資源和水資源非常有限。這一情況對(duì)該區(qū)域工農(nóng)業(yè)發(fā)展來(lái)說(shuō)，是十分不利的。眾所周知，大理河流域的徑流呈現(xiàn)顯著下降[15]，這使得當(dāng)?shù)厮Y源的供需矛盾進(jìn)一步加劇，同時(shí)伴隨著嚴(yán)重的泥沙問(wèn)題。在過(guò)去的許多時(shí)間，關(guān)于大理河流域的氣候變化和人類活動(dòng)影響下的徑流變化和輸沙變化的研究已經(jīng)取得了顯著的成果。但是，大部分研究重視探索造成徑流變化和輸沙變化等單個(gè)變量變化的氣候變化和人類活動(dòng)因素，關(guān)于徑流和輸沙之間關(guān)系的研究還需要不斷探索和深入研究。徑流和輸沙之間的水沙關(guān)系作為流域物質(zhì)循環(huán)中至關(guān)重要的環(huán)節(jié)，二者之間的關(guān)系具有緊密的聯(lián)系，這是水土保持研究關(guān)注的重點(diǎn)問(wèn)題之一[16]。一直以來(lái)研究徑流輸沙之間的關(guān)系對(duì)于優(yōu)化流域水土流失治理工作，評(píng)估流域水土保持效益都是必要的。一旦流域徑流和輸沙和之間的關(guān)系存在突變，這意味著流域的水文機(jī)制和水沙關(guān)系已經(jīng)發(fā)生了改變，對(duì)于工程設(shè)計(jì)工作來(lái)說(shuō)也必須進(jìn)行改變，才能有效地興利除害。因此，基于多變量視角對(duì)大理河流域徑流和輸沙的變化關(guān)系進(jìn)行探索至關(guān)重要，對(duì)在氣候變化和人類活動(dòng)的巨大影響條件下指導(dǎo)當(dāng)?shù)剡M(jìn)行水資源規(guī)劃和水土保持具有重要意義。

本文通過(guò)Mann-Kendall 檢驗(yàn)研究徑流輸沙的變化趨勢(shì)，基于耦合協(xié)調(diào)度理論和Pettitt 檢驗(yàn)方法構(gòu)建流域徑流和輸沙之間關(guān)系的突變點(diǎn)的識(shí)別方法并進(jìn)行診斷，并通過(guò)Copula 方法進(jìn)行了驗(yàn)證，探討了造成大理河流域徑流和輸沙關(guān)系變化的可能驅(qū)動(dòng)因素。研究有助于進(jìn)一步深化對(duì)大理河流域水沙關(guān)系變化的認(rèn)識(shí)，為大理河流域水土保持工作開展提供一定的科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

大理河流域位于中國(guó)黃河流域中游（109°14′～110°13′E，37°30′～37°56′N），綏德水文站為流域出口水文站[17]，如圖1 所示。大理河為無(wú)定河的第二大支流，起源于靖邊縣白于山東測(cè)的五臺(tái)山，流經(jīng)橫山縣、子洲縣、綏德縣，在綏德縣東北部注入無(wú)定河。大理河全長(zhǎng)159.9 km，流域面積3 904.24 km2，比降3.16‰，流域地面坡度較大，大部分在15°以上，主要支流有小理河、岔巴溝、駝耳巷溝等。流域西高東低，地形起伏較大，海拔796～1 744 m。流域主要土壤類型為黃綿土和新積土，屬于黃土峁?fàn)钋鹆隃羡謪^(qū)。流域內(nèi)坡長(zhǎng)多大于200 m，溝谷切割深度多大于50 m，使得流域形成了許多坡陡溝深的小流域。流域土壤侵蝕為嚴(yán)重，經(jīng)過(guò)幾十年的水土保持綜合治理，水土流失量明顯下降。

1.2 數(shù)據(jù)來(lái)源

本文所使用的徑流和輸沙數(shù)據(jù)來(lái)源于黃河水文年鑒，時(shí)間序列為1960—2015 年。降水?dāng)?shù)據(jù)是來(lái)源于國(guó)家地球系統(tǒng)科學(xué)數(shù)據(jù)中心的柵格降水量數(shù)據(jù)（http：//loess.geodata.cn），時(shí)間序列為1960—2017 年，空間分辨率為1 km×1 km，時(shí)間分辨率為月，該數(shù)據(jù)使用全國(guó)496 個(gè)獨(dú)立氣象觀測(cè)點(diǎn)數(shù)據(jù)進(jìn)行了驗(yàn)證，驗(yàn)證結(jié)果可信，該數(shù)據(jù)可以準(zhǔn)確的反映全流域降水的變化分布情況。淤地壩數(shù)據(jù)來(lái)源于2011 年第一次全國(guó)水利普查。中國(guó)年度植被指數(shù)（NDVI）空間分布數(shù)據(jù)集是基于連續(xù)時(shí)間序列的 SPOT/VEGETATION NDVI 衛(wèi)星遙感數(shù)據(jù)（http：//www.resdc.cn/data.aspx?DATAID=257），采用最大值合成法生成的1998 年以來(lái)的年度植被指數(shù)數(shù)據(jù)集。該數(shù)據(jù)集有效反映了全國(guó)各地區(qū)在空間和時(shí)間尺度上的植被覆蓋分布和變化狀況，對(duì)植被變化狀況監(jiān)測(cè)、植被資源合理利用和其他生態(tài)環(huán)境相關(guān)領(lǐng)域的研究有十分重要的參考意義。

圖1 陜西省大理河流域地理位置 Fig.1 Location of the Dali River Basin, Shaanxi province

1.3 研究方法

1.3.1 Mann-Kendall 檢驗(yàn)

Mann-Kendall（M-K）檢驗(yàn)[18]是一種非參數(shù)的判斷數(shù)據(jù)序列變化趨勢(shì)的分析方法，常被廣泛使用在氣象和水文序列的趨勢(shì)檢驗(yàn)中，本文采用Mann-Kendall 檢驗(yàn)捕捉徑流量、輸沙量、NDVI 以及其他變量的變化趨勢(shì)。

1.3.2 耦合協(xié)調(diào)理論

耦合最早是作為一個(gè)物理學(xué)的概念，指2 個(gè)或2 個(gè)以上系統(tǒng)或運(yùn)動(dòng)形式通過(guò)各種相互作用而彼此影響的現(xiàn)象[19]。隨著研究的發(fā)展，耦合概念被逐步引入到生態(tài)學(xué)領(lǐng)域中，并被廣泛使用[20]。多個(gè)系統(tǒng)相互作用耦合度模型可以用一下模型表示[21]

式中Cn為n 元系統(tǒng)的耦合度；u1…un分別為第一個(gè)子系統(tǒng)到第n 個(gè)子系統(tǒng)對(duì)總系統(tǒng)有序度的貢獻(xiàn)，計(jì)算方法如下

式中ui為第i 個(gè)子系統(tǒng)對(duì)總系統(tǒng)有序度的貢獻(xiàn)；uij為第i個(gè)子系統(tǒng)中第j 個(gè)指標(biāo)的歸一化值；wij為第i 個(gè)子系統(tǒng)中第j 個(gè)指標(biāo)的權(quán)重，每個(gè)子系統(tǒng)中指標(biāo)的權(quán)重計(jì)算使用熵權(quán)法進(jìn)行計(jì)算。

在計(jì)算每個(gè)子系統(tǒng)的熵權(quán)時(shí)，必須先進(jìn)行歸一化處理，這里采用最大最小值法進(jìn)行歸一化處理：

當(dāng)uij為數(shù)值越大對(duì)系統(tǒng)越好時(shí)（正向歸一化）

當(dāng)uij為數(shù)值越小對(duì)系統(tǒng)越好時(shí)（負(fù)向歸一化）

式中xij為子系統(tǒng)第i 個(gè)指標(biāo)的第j 個(gè)時(shí)序的值。

根據(jù)多元系統(tǒng)的耦合度模型，容易得到徑流和輸沙2個(gè)子系統(tǒng)之間的耦合度模型，可以表示成以下形式

式中C 為徑流輸沙二元系統(tǒng)的耦合度。在本研究中徑流和輸沙子系統(tǒng)的指標(biāo)為12 個(gè)月的月徑流量和月輸沙量，每個(gè)指標(biāo)的權(quán)重計(jì)算采用熵權(quán)法，假定徑流越大對(duì)系統(tǒng)越好，輸沙越少對(duì)系統(tǒng)越好。

由于耦合度指標(biāo)在有些情況下卻很難反映出子系統(tǒng)整體“功效”與“協(xié)同”效應(yīng)，耦合度各子系統(tǒng)指標(biāo)的上下限取自各指標(biāo)的極值，而極值存在動(dòng)態(tài)、不平衡的特性，單純依靠耦合度判別有可能產(chǎn)生誤導(dǎo)。為此，進(jìn)一步采用耦合協(xié)調(diào)度模型來(lái)評(píng)判徑流子系統(tǒng)和輸沙系統(tǒng)交互耦合的協(xié)調(diào)程度，計(jì)算公式如下

式中D 為耦合協(xié)調(diào)度；C 為耦合度；T 為徑流子系統(tǒng)與輸沙子系統(tǒng)的綜合調(diào)和指數(shù)，它反映徑流子系統(tǒng)與輸沙子系統(tǒng)的整體協(xié)同效應(yīng)或貢獻(xiàn)；a、b 為待定系數(shù)，實(shí)際中常常認(rèn)為兩個(gè)子系統(tǒng)的重要性相同[21]，所以a=b=0.5。

1.3.3 邊際分布函數(shù)與Copula 函數(shù)

邊際分布函數(shù)[22]是統(tǒng)計(jì)學(xué)中一種描述單變量與其概率分布的函數(shù)，不同的分布函數(shù)可以刻畫不同的變量的特征。在本研究中，指數(shù)分布函數(shù)、耿貝爾分布函數(shù)、皮爾遜三型分布函數(shù)和廣義極值分布函數(shù)被使用來(lái)擬合徑流和輸沙數(shù)據(jù)，同時(shí)通過(guò)計(jì)算4 種理論分布函數(shù)累積概率和經(jīng)驗(yàn)累積概率的R2來(lái)優(yōu)選出最優(yōu)的理論分布函數(shù)。

Copula 函數(shù)是一種可以用來(lái)描述多個(gè)變量聯(lián)合分布的一種連接函數(shù)[22]，由于它可以連接任何2 個(gè)邊際分布函數(shù)，因此使用起來(lái)十分方便。Copula 函數(shù)具有很多家族，其中阿基米德家族Copula 函數(shù)因?yàn)榫哂休^少的參數(shù)，所以被廣泛使用。本研究中通過(guò)AIC 準(zhǔn)則，從阿基米德家族Copula 函數(shù)的Frank Copula、Gumbel Copula 和Clayton Copula 3 種Copula 這種選出最優(yōu)的Copula 來(lái)連接大理河流域徑流和輸沙的最優(yōu)邊際分布函數(shù)，構(gòu)建出最優(yōu)的徑流輸沙最優(yōu)Copula 分布函數(shù)。

1.3.4 Pettitt 突變點(diǎn)檢驗(yàn)

Pettitt 檢驗(yàn)法[23]采用Mann-Whitney 中Ut,n值檢驗(yàn)同一總體中2 個(gè)樣本X1,…,Xt和Xt+1,…,XN，Pettitt 檢驗(yàn)的零假設(shè)為沒(méi)有變化點(diǎn)，當(dāng)|Ut,n|取最大值時(shí)對(duì)應(yīng)的Xt被認(rèn)為是可能的突變點(diǎn)。當(dāng)P≤0.05 時(shí)認(rèn)為數(shù)據(jù)中存在均值變異點(diǎn)。其顯著性水平可由下式計(jì)算

當(dāng)P≤0.05 時(shí)認(rèn)為數(shù)據(jù)中存在均值變異點(diǎn)。

2 結(jié)果與分析

2.1 徑流和輸沙變化趨勢(shì)

大理河流域1960—2015 年各月及年徑流量和輸沙量的變化趨勢(shì)計(jì)算結(jié)果如表1 中所示。徑流在3、4、5、8、10、11 月以及年尺度均呈現(xiàn)顯著的下降趨勢(shì)，其中3、4、10 月以及年尺度的M-K 統(tǒng)計(jì)量絕對(duì)值較大，表明徑流下降幅度較大，1 月的徑流呈現(xiàn)不顯著的增加趨勢(shì)，其他月份的徑流呈現(xiàn)不顯著的下降趨勢(shì)。在1960—2015年間，各月輸沙以及年尺度輸沙均呈現(xiàn)顯著的下降趨勢(shì)，相比同一月份的徑流變化趨勢(shì)而言，輸沙的M-K 統(tǒng)計(jì)量絕對(duì)值均高于徑流，說(shuō)明輸沙的下降幅度更大。在1960—2015 年間大理河流域的徑流與輸沙變化均存在明顯的下降趨勢(shì)，但徑流與輸沙的變化幅度呈現(xiàn)出不協(xié)調(diào)的特點(diǎn)。

表1 1960—2015 年大理河流域徑流量和輸沙量變化趨勢(shì) Table 1 Trends of runoff and sediment load in the Dali River Basin from 1960 to 2015

2.2 徑流輸沙關(guān)系變異診斷

2.2.1 徑流子系統(tǒng)與輸沙子系統(tǒng)的指標(biāo)變化

本研究中認(rèn)為徑流越大對(duì)系統(tǒng)越好，輸沙越少對(duì)系統(tǒng)越好，選取每年12 個(gè)月的徑流量和輸沙量分別作為兩個(gè)子系統(tǒng)的構(gòu)成指標(biāo)，對(duì)每個(gè)指標(biāo)即每個(gè)月的徑流量和輸沙量進(jìn)行了歸一化處理，使用熵權(quán)法計(jì)算每個(gè)指標(biāo)的權(quán)重，如表2。

表2 的結(jié)果顯示，1960—2015 年間6 月的徑流量在一年中權(quán)重最大，11 月的徑流量權(quán)重最小，4、5、7、8、9、10 月的徑流量權(quán)重較大。在1960—2015 年間，輸沙量在8 月的權(quán)重最大，2 月和7 月的權(quán)重較大。徑流的權(quán)重值變異系數(shù)為0.58，輸沙的權(quán)重值變異系數(shù)為0.51，徑流量年分配權(quán)重比輸沙更為不穩(wěn)定。

表2 1960—2015 年大理河流域徑流量和輸沙量權(quán)重 Table 2 Weights of runoff and sediment load in the Dali River Basin from 1960 to 2015

2.2.2 徑流和輸沙耦合協(xié)調(diào)度的變化趨勢(shì)與突變點(diǎn)

根據(jù)耦合協(xié)調(diào)度理論計(jì)算了大理河流域徑流子系統(tǒng)和輸沙子系統(tǒng)的耦合協(xié)調(diào)度，并通過(guò)Pettitt 檢驗(yàn)方法識(shí)別了徑流輸沙耦合協(xié)調(diào)度的突變點(diǎn)年份，它反映了大理河流域徑流和輸沙之間關(guān)系的變異特征（圖2）。

圖2a所示的是1960—2015年間大理河徑流輸沙耦合協(xié)調(diào)度的時(shí)間序列圖，最大值為0.75，最小值為0.42，根據(jù)M-K 趨勢(shì)檢驗(yàn)結(jié)果表明徑流輸沙耦合協(xié)調(diào)度呈現(xiàn)顯著下降趨勢(shì)，說(shuō)明大理河流域徑流輸沙關(guān)系協(xié)調(diào)程度不斷下降。圖2b 表明大理河流域徑流輸沙耦合協(xié)調(diào)度的突變點(diǎn)為1996 年，在1996 年徑流和輸沙之間存在顯著突變點(diǎn)，通過(guò)對(duì)突變點(diǎn)前后的2 個(gè)耦合協(xié)調(diào)度的子序列進(jìn)一步進(jìn)行Pettitt 檢驗(yàn)，發(fā)現(xiàn)不存在顯著突變點(diǎn)，這表明大理河流域徑流和輸沙關(guān)系只在1996 年存在明顯的不同。

2.2.3 基于Copula 函數(shù)的徑流輸沙突變點(diǎn)驗(yàn)證分析

為了對(duì)基于耦合協(xié)調(diào)度理論識(shí)別的徑流和輸沙關(guān)系變異點(diǎn)診斷結(jié)果進(jìn)行驗(yàn)證，本研究采樣Copula 函數(shù)構(gòu)建了徑流和輸沙之間的最優(yōu)聯(lián)合分布累積概率，并識(shí)別了突變點(diǎn)（圖3）。

圖2 1960—2015 年大理河流域徑流量和輸沙量耦合協(xié)調(diào)度變化及其突變點(diǎn) Fig.2 Changes of coupling coordination of runoff and sediment load and its change-points in the Dali River Basin from 1960 to 2015

圖3 1960—2015 年大理河流域徑流和輸沙之間的最優(yōu)聯(lián)合分布累積概率 Fig.3 Cumulative probability of optimal Copula joint distribution of runoff and sediment load in Dali River Basin from 1960 to 2015

指數(shù)分布函數(shù)、耿貝爾分布函數(shù)、皮爾遜三型分布函數(shù)和廣義極值分布函數(shù)等4 種理論分布函數(shù)被使用來(lái)擬合大理河流域1960—2015 年的年徑流和年輸沙數(shù)據(jù)，通過(guò)K-S 檢驗(yàn)分別檢驗(yàn)了年徑流和年輸沙數(shù)據(jù)，發(fā)現(xiàn)年徑流和年輸沙都通過(guò)了K-S 檢驗(yàn)，這表明4 中理論分布函數(shù)都可以用來(lái)描述年徑流和年輸沙，同時(shí)使用線性矩法估計(jì)了4 種理論分布函數(shù)的參數(shù)，分別通過(guò)計(jì)算年徑流和年輸沙的經(jīng)驗(yàn)累計(jì)概率與4 種理論分布函數(shù)的累積概率之間的R2，根據(jù)R2最大原則選取該分布函數(shù)為最優(yōu)的分布函數(shù)。結(jié)果表明，大理河流域1960—2015 年年徑流的最優(yōu)分布函數(shù)為皮爾遜三型分布，年輸沙的最優(yōu)分布函數(shù)為指數(shù)分布。圖3a、3b 分別為大理河流域1960—2015 年年徑流量、輸沙量的最優(yōu)分布函數(shù)累計(jì)概率與經(jīng)驗(yàn)累積概率的時(shí)間序列圖?；贏IC 準(zhǔn)則，從阿基米德家族Copula 函數(shù)的Frank Copula、Gumbel Copula 和Clayton Copula 三種Copula 中選出Frank Copula 作為最優(yōu)的Copula 來(lái)連接大理河流域徑流和輸沙的最優(yōu)邊際分布函數(shù)，計(jì)算出來(lái)徑流和輸沙之間的聯(lián)合分布函數(shù)累積概率（圖3c）。根據(jù)M-K 趨勢(shì)檢驗(yàn)對(duì)徑流和輸沙之間的聯(lián)合分布函數(shù)累積概率的分析可知，大理河流域徑流和輸沙之間的聯(lián)合分布函數(shù)累積概率呈現(xiàn)顯著下降趨勢(shì)，Pettitt 檢驗(yàn)對(duì)徑流和輸沙之間的聯(lián)合分布函數(shù)累積概率分析結(jié)果表明大理河流域徑流和輸沙之間關(guān)系在1996 年存在顯著突變點(diǎn)（圖3d）。2 種徑流輸沙關(guān)系突變點(diǎn)的診斷方法表明，基于耦合協(xié)調(diào)度的徑流輸沙關(guān)系變異診斷方法可以準(zhǔn)確的識(shí)別出徑流輸沙的關(guān)系突變點(diǎn)。

2.3 變異前后徑流輸沙變化特征分析

2.3.1 徑流輸沙變化特征分析

對(duì)比年徑流量-年輸沙量關(guān)系突變點(diǎn)前后，徑流與輸沙的年內(nèi)和年際的變化特征，可以發(fā)現(xiàn)大理河流域1960—1996 年年徑流量為1.50 億m3和1996—2015 年徑流量為1.08 億m3，1996—2015 年期間相比1960—1996 年期間徑流量減少27.92%，減少幅度較小。1960—1996 年年輸沙量為0.43 億t 和1996—2015 年輸沙量為0.17 億t，1996—2015 年期間相比1960—1996 年期間輸沙量減少57.11%，輸沙的減少幅度比徑流的減少幅度高29.19 個(gè)百分點(diǎn)。1996—2015 年各月徑流量、輸沙量與1960 —1996 年各月徑流量、輸沙量相比較發(fā)現(xiàn)，除了12月徑流量相比突變前有所增加，其余各月徑流量均在減少；輸沙量在12 個(gè)月中均發(fā)生了明顯的減少現(xiàn)象，特別是在8 月份，徑流和輸沙的減少量最大（圖4）。

圖4 突變點(diǎn)前后徑流量和輸沙量年內(nèi)分配 Fig.4 Annual distribution of runoff and sediment load before and after the change-point

2.3.2 徑流輸沙關(guān)系變異

徑流和輸沙關(guān)系突變點(diǎn)前后徑流和輸沙的散點(diǎn)圖，如圖5。由徑流和輸沙關(guān)系突變點(diǎn)前的徑流和輸沙相關(guān)圖（圖5a）可以看出，徑流輸沙關(guān)系可以被線性回歸方程描述，徑流和輸沙的決定系數(shù)R2為0.867 3（P<0.05），在統(tǒng)計(jì)學(xué)意義上，決定系數(shù)R2表示自變量可以解釋因變量的程度，因此，在1996 年之前大理河徑流對(duì)輸沙貢獻(xiàn)程度高達(dá)86.73%，單位的徑流輸沙能力較強(qiáng)。圖5b 突變點(diǎn)后大理河的徑流和輸沙相關(guān)圖表明，在1996—2015 年徑流輸沙關(guān)系也可以被線性回歸方程很好的描述，徑流和輸沙的決定系數(shù)R2為0.823 1（P<0.05），在1996 年后大理河徑流對(duì)輸沙貢獻(xiàn)程度降低到82.31%，下降幅度為5.10%，單位的徑流輸沙能力降低。

圖5 大理河流域突變點(diǎn)前后徑流輸沙關(guān)系 Fig.5 Relationship between runoff and sediment load before and after the change-point in Dali River Basin

2.4 氣候變化和人類活動(dòng)對(duì)徑流輸沙變化的貢獻(xiàn)率

對(duì)徑流和輸沙變化影響的主要驅(qū)動(dòng)因素包括氣候因素和下墊面因素[24]。本研究使用雙累積曲線分析兩者對(duì)大理河流域1960—2015 年徑流和輸沙變化的影響。

如圖6 所示為1960—2015 年大理河流域年降水量和年徑流量雙累積曲線（圖6a）、年降水量和年輸沙量（圖 6b）。從圖6a 中可以看出，大理河流域年降水量和年徑流量之間存在較弱的非一致性。圖6b 表明，大理河流域年降水量和年輸沙量之間存在明顯的非一致性關(guān)系。結(jié)合圖6，將降水和徑流、輸沙關(guān)系劃分為3 個(gè)時(shí)段，分別為1960—1971、1972—1996、1997—2015 年。將1960—1971 作為基準(zhǔn)期，認(rèn)為該時(shí)期降水徑流關(guān)系和降水輸沙關(guān)系未發(fā)生變化，具有一致性，定量分析氣候變化和人類會(huì)動(dòng)對(duì)徑流和輸沙變化的貢獻(xiàn)率，如表3 所示。在1972—1996、1997—2015 年2 個(gè)時(shí)期人類活動(dòng)對(duì)徑流和輸沙減少的貢獻(xiàn)率均遠(yuǎn)遠(yuǎn)高于氣候變化的影響，特別是在1997—2015 年時(shí)期，氣候變化對(duì)徑流和輸沙變化表現(xiàn)出增加的作用，這表明下墊面的人類活動(dòng)是導(dǎo)致大理河流域徑流和泥沙減少的主要驅(qū)動(dòng)力，這與劉曉燕等[25-26]人的研究結(jié)果一致。

表3 氣候變化和人類活動(dòng)對(duì)大理河流域徑流輸沙貢獻(xiàn)率計(jì)算結(jié)果 Table 3 Calculation results of climate change and human activities' contribution to runoff and sediment load in the Dali River Basin

圖6 1960—2015 年大理河流域年降水量-徑流量和降水量-輸沙量雙累積曲線 Fig.6 Double accumulation curve of annual precipitation-runoff and precipitation-sediment load in the Dali River Basin from 1960 to 2015

3 討 論

3.1 淤地壩建設(shè)對(duì)徑流輸沙的影響

研究表明，黃河流域的年降水量并未發(fā)生顯著增加[27]，根據(jù)1960—2015 年大理河流域年降雨的M-K 趨勢(shì)檢驗(yàn)結(jié)果表明，M-K 趨勢(shì)檢驗(yàn)統(tǒng)計(jì)量的絕對(duì)值均小于1.96，這意味著全流域年降水量均呈不顯著的變化趨勢(shì)（圖7b）。對(duì)于黃土高原地區(qū)來(lái)說(shuō)，人類活動(dòng)變化中影響最為廣泛的就是20 世紀(jì)50 年代開始在黃土高原地區(qū)實(shí)施的水土保持工程，該工程主要包括林草工程、梯田建設(shè)和淤地壩工程的修建[28-29]。其中，因?yàn)橛俚貕尉哂袛r沙淤地的直接經(jīng)濟(jì)效益，因而被黃土高原地區(qū)廣泛推廣使用，淤地壩工程的修建歷時(shí)最長(zhǎng)，修建規(guī)模最大[30]。圖8 中所示為大理河流域1960—2011 年徑流量和輸沙量變化，其中實(shí)測(cè)輸沙量為綏德水文站的觀測(cè)值，實(shí)際輸沙量為淤地壩逐年攔沙量和實(shí)測(cè)輸沙量的加和。實(shí)測(cè)輸沙量的變化基本與實(shí)際輸沙的變化具有一致性。因此，實(shí)際輸沙量與徑流量的關(guān)系和實(shí)測(cè)輸沙量與徑流量的關(guān)系基本保持一致。根據(jù)李宗善等[31]的研究，黃土高原地區(qū)淤地壩已經(jīng)超過(guò)10 萬(wàn)座以上，效益明顯，直接發(fā)揮了攔截泥沙淤泥造地的功能，目前已經(jīng)截留了280 億t 水土流失總量，2002 年淤地壩農(nóng)田規(guī)模達(dá)到3 200 km2，糧食產(chǎn)量是梯田的2～3 倍。在大理河“7·26”特大暴雨洪水事件中，建設(shè)有完備淤地壩系的韭園溝流域相比對(duì)照流域洪峰流量消減達(dá)到8 倍，洪水含沙量只有對(duì)照流域的三分之一。根據(jù)2011 年全國(guó)第一次水利普查，大理河流域共有骨干淤地壩279 座（圖7a）。淤地壩的修建高峰期在20 世紀(jì)70 年代，同時(shí)2000 年后水利部實(shí)施淤地壩亮點(diǎn)工程，淤地壩的修建又出現(xiàn)了第二次高峰。淤地壩的修建直接攔截了大理河流域大量的泥沙，對(duì)流域徑流產(chǎn)生了再調(diào)節(jié)和再分配的作用，徑流的變化造成產(chǎn)沙輸沙能力的減弱，因而間接減少了輸沙量。

3.2 植被恢復(fù)對(duì)徑流輸沙變化的影響

自從1999 年中國(guó)政府在黃土高原地區(qū)實(shí)施了退耕還林政策，大規(guī)模的植樹造林工程顯著的增加了黃土高原的植被覆蓋度[32-33]。特別是在大理河流域，NDVI 從1999—2015 年明顯增加（圖9）。根據(jù)M-K 檢驗(yàn)分析，1999—2015 年大理河流域NDVI 均值呈現(xiàn)顯著增加趨勢(shì)，與徑流和輸沙的變化趨勢(shì)相反。根據(jù)ULSE 土壤侵蝕的模型的原理，隨著植被的增加土壤侵蝕程度可以降低，植被在一定程度上阻止了泥沙的產(chǎn)生，減少了輸沙量[34]。根據(jù)梁越等[35]的研究，退耕還林后黃河河龍區(qū)間土壤侵蝕強(qiáng)度減弱，退耕還林工程對(duì)河龍區(qū)間減少作用從大至小呈現(xiàn)：中部-南部-北部。同時(shí)，研究表明由于近60年河龍區(qū)間降雨和降雨侵蝕力并未發(fā)生明顯變化趨勢(shì)，退耕還林后植被增加可能是小流域泥沙減少的主要推動(dòng)力。

圖7 陜西省大理河流域骨干淤地壩空間分布與降雨量變化分布 Fig.7 Spatial distribution of check dam and precipitation in the Dali River Basin, Shaanxi province

圖8 大理河流域徑流量與輸沙量變化 Fig.8 Changes of runoff and sediment load in Dali River Basin

圖9 1999—2015 年大理河流域NDVI 均值變化 Fig.9 NDVI change in Dali River Basin from 1999 to 2015

4 結(jié) 論

1）在1960—2015 年，大理河流域的徑流與輸沙在各月及年尺度均存在顯著下降趨勢(shì)，而且徑流與輸沙的變化趨勢(shì)不協(xié)調(diào)不同步。

2）本研究構(gòu)建了基于耦合協(xié)調(diào)度的徑流輸沙關(guān)系突變點(diǎn)診斷方法，識(shí)別出大理河流域徑流輸沙關(guān)系在1996年存在顯著突變點(diǎn)，同時(shí)基于徑流輸沙Copula 聯(lián)合分布累計(jì)概率驗(yàn)證了該突變點(diǎn)確實(shí)存在，表明1996 年是大理河流域徑流輸沙關(guān)系發(fā)生變化的重要節(jié)點(diǎn)，本文方法具有很好的識(shí)別準(zhǔn)確性。

3）大理河流域徑流輸沙關(guān)系突變后，徑流和輸沙均有所減少，徑流量減少27.92%，輸沙量減少57.11%，輸沙減少幅度大于徑流。徑流輸沙關(guān)系突變后，大理河徑流對(duì)輸沙貢獻(xiàn)程度相比徑流輸沙關(guān)系突變前，下降幅度為5.10%，徑流輸沙能力降低。

4）人類活動(dòng)是影響大理河流域徑流輸沙變化的主要驅(qū)動(dòng)力，淤地壩的修建影響著徑流和輸沙的變化，自退耕還林政策實(shí)施后，從1999—2015 年大理河流域NDVI 均值呈現(xiàn)顯著增加趨勢(shì)，對(duì)徑流和輸沙變化產(chǎn)生重要作用。