基于Copula分析生態(tài)恢復(fù)對(duì)黃土高原水沙頻率耦合關(guān)系的影響

薛 帆，易海杰，宋松柏，張曉萍,

(1.西北農(nóng)林科技大學(xué)水土保持研究所，黃土高原土壤侵蝕與旱地農(nóng)業(yè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 楊凌 712100；2.中國(guó)科學(xué)院教育部水土保持與生態(tài)環(huán)境研究中心，陜西 楊凌 712100；3.西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院，陜西 楊凌 712100)

黃土高原自然環(huán)境復(fù)雜，地形破碎，土質(zhì)疏松，植被破壞后極易發(fā)生水土流失，黃河年均輸沙量曾達(dá)到16億t，其中90%來源于黃土高原。自1970s，黃土高原開展了大規(guī)模的水土流失治理，1999年后又大力實(shí)施退耕還林(草)措施，土地利用和地表覆被發(fā)生劇烈變化，流域水沙數(shù)量、水沙關(guān)系成為研究熱點(diǎn)。

水沙關(guān)系，其內(nèi)涵是指某時(shí)空尺度上流域徑流量與輸沙或泥沙含量的關(guān)系，受自然因素，如降水、土壤、地形、植被等的影響，而大規(guī)模人類活動(dòng)導(dǎo)致的下墊面條件變化，如大規(guī)模植被恢復(fù)工程、水土流失治理、農(nóng)業(yè)灌溉等，對(duì)流域徑流和泥沙過程均產(chǎn)生深刻影響。一方面表現(xiàn)在徑流量和泥沙數(shù)量的協(xié)同增加或減少效應(yīng)及密切程度；另一方面，還表現(xiàn)在其發(fā)生頻率的協(xié)同性。目前，已有許多學(xué)者探究了黃土高原典型流域生態(tài)建設(shè)對(duì)水沙數(shù)量關(guān)系及其密切程度的調(diào)控效應(yīng)，如劉寶元、Zhang等研究表明，生態(tài)環(huán)境治理下，黃河徑流量由1919—1959年的年均426億m減少到2010年以來的近300億m，同時(shí)期輸沙量由年均15.8億t減少至不到2億t。針對(duì)黃土高原流域水沙關(guān)系變化的研究，多是基于水、沙數(shù)量的角度進(jìn)行理解與分析。Zhao等、Wang等研究表明，2000年以來黃土高原高含沙量明顯降低；Zheng等研究表明，在流域尺度上，植被和坡面水保措施僅通過減水來減沙；Zhang等、Gao等基于經(jīng)驗(yàn)回歸關(guān)系式探究水沙關(guān)系的變化，2000年后水、沙數(shù)量相關(guān)性降低。而水沙關(guān)系在發(fā)生頻率協(xié)同性上的研究，如郭愛軍等基于Copula函數(shù)建立涇河流域水沙聯(lián)合分布模型，結(jié)果表明，受生態(tài)建設(shè)影響，水沙豐枯各種遭遇組合頻率分布更為均勻。以上研究主要從數(shù)量角度分析黃土高原典型流域水沙關(guān)系對(duì)生態(tài)恢復(fù)的響應(yīng)特征，而對(duì)于其頻率耦合關(guān)系的研究開展得還不夠。

要深入剖析黃土高原流域水沙頻率耦合關(guān)系的演變特征，需全面考慮徑流、輸沙2個(gè)相關(guān)變量的聯(lián)合分布狀況。Copula函數(shù)是定義域?yàn)閇0，1]均勻分布的多維聯(lián)合分布函數(shù)，它可以根據(jù)多個(gè)隨機(jī)變量和其邊緣分布函數(shù)的相關(guān)關(guān)系構(gòu)建聯(lián)合分布模型，具有計(jì)算方便、形式簡(jiǎn)單的優(yōu)越性。目前，已經(jīng)有學(xué)者將Copula函數(shù)應(yīng)用于多變量水文分析計(jì)算的研究中，如水文隨機(jī)變量概率分布、洪水頻率分析、降雨頻率分析、水文干旱特征分析、水文隨機(jī)模擬等，均表明根據(jù)Copula理論建立的聯(lián)合分布模型模擬精度較高。

黃土高原丘陵溝壑區(qū)，不僅是黃河下游多沙粗砂的主要策源地，也是黃土高原植被恢復(fù)最為顯著的區(qū)域，開展該區(qū)生態(tài)環(huán)境變化下流域水沙特征演變規(guī)律具有很好的典型性。本文基于典型丘陵溝壑區(qū)北洛河上游劉家河水文站1960—2019年的日徑流量和輸沙量數(shù)據(jù)，采用Copula函數(shù)和經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition)，建立水沙聯(lián)合分布模型，分析多時(shí)間尺度下流域水沙頻率耦合關(guān)系對(duì)水土保持與生態(tài)恢復(fù)的響應(yīng)，為開展黃土高原水土流失治理規(guī)劃、生態(tài)恢復(fù)等提供參考依據(jù)。

1 研究區(qū)概況

北洛河是渭河的一級(jí)支流，黃河的二級(jí)支流，流經(jīng)陜西省榆林、延安、銅川、渭南及甘肅省慶陽等5個(gè)地(市)的18個(gè)縣(區(qū))，于大荔縣東南注入渭河，地處34°42′—37°19′N，107°33′—110°10′E。劉家河水文站位于北洛河流域上游，其控制面積為7 325 km，約占流域總面積(26 905 km)的1/4，地貌類型為典型的黃土丘陵溝壑區(qū)，屬暖溫帶半干旱氣候，年平均降水量452 mm，多集中在7—9月；黃綿土為主要土壤類型，抗侵蝕能力差，極易發(fā)生嚴(yán)重的水土流失。自1970s以來，北洛河上游大規(guī)模實(shí)施生態(tài)建設(shè)和水土保持措施，如造林種草、修建淤地壩和梯田，2006底，上游地區(qū)各項(xiàng)水保措施面積達(dá)26.64萬hm(表1)。1999年又推行退耕還林(草)政策，林草植被蓋度由1970年的59.8%增加到2019年的81.1%。經(jīng)過持續(xù)的生態(tài)建設(shè)，與1960s相比，劉家河水文站2010s徑流量、輸沙量分別減少約1/2和9/10，人類活動(dòng)對(duì)徑流量、輸沙量變化貢獻(xiàn)分別為94.7%，91.6%；不同發(fā)生頻率的徑流、輸沙和含沙量存在不同程度的顯著減少趨勢(shì)，豐水期(1%～5%)的減少尤其顯著。

表1 北洛河上游水土保持措施累積面積[24]

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

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

本文收集了北洛河流域上游劉家河水文站1960—2019年實(shí)測(cè)日徑流量、輸沙量，并構(gòu)建月徑流量()和輸沙量()序列進(jìn)行研究分析，所用數(shù)據(jù)來源于黃河中游水文資料及年鑒。

2.2 數(shù)據(jù)處理方法

2.2.1 水沙數(shù)量時(shí)間變化趨勢(shì)診斷 Pettitt非參數(shù)檢驗(yàn)法廣泛應(yīng)用于水文序列突變檢驗(yàn)。首先建立Mann-Whitney統(tǒng)計(jì)量(,)：

(1)

式中：和=2，…，。

令-=，則sgn()值為：

(2)

檢驗(yàn)統(tǒng)計(jì)量()為：

=Max1≤≤|，|

(3)

統(tǒng)計(jì)量相關(guān)概率的顯著性檢驗(yàn)公式為：

(4)

2.2.2 二維Copula聯(lián)合分布的建立

(1)確定邊緣分布函數(shù)。通過皮爾遜III型曲線確定邊緣函數(shù)：

(5)

式中：為皮爾遜III型曲線的形狀參數(shù)；為尺度參數(shù)；為位置參數(shù)；Γ()為的伽馬函數(shù)。

、和的矩法估計(jì)計(jì)算公式為：

(6)

(2)邊緣分布函數(shù)檢驗(yàn)。采用檢驗(yàn)方法對(duì)邊緣分布函數(shù)進(jìn)行假設(shè)檢驗(yàn)。設(shè)樣本總體的分布函數(shù)為()且已知。將原假設(shè)為：總體的分布函數(shù)為(：，，…，)

檢驗(yàn)統(tǒng)計(jì)量可表示為：

(7)

當(dāng)為真時(shí)近似地有：

(8)

其拒絕域?yàn)椋?/p>

≥(--1)

(9)

式中：為顯著水平。

(3)Copula函數(shù)參數(shù)估計(jì)與選擇。本研究選取Frank Copula、Gumbel-Hougaard(GH)Copula、Clayton Copula 3種函數(shù)構(gòu)建二維水沙聯(lián)合分布模型：

①Frank Copula

(10)

②Gumbel-Hougaard(GH)Copula

(11)

③Clayton Copula

(12)

式中：為Copula函數(shù)的參數(shù)。

為確定3種Copula函數(shù)參數(shù)，首先需得到徑流量和輸沙量的Kendall秩相關(guān)系數(shù)():

(13)

式中：、、、為觀測(cè)數(shù)據(jù)；sgn為符號(hào)函數(shù)。

Copula函數(shù)的參數(shù)與τ存在的關(guān)系為：

①Frank Copula

(14)

②Gumbel-Hougaard(GH)Copula

(15)

③Clayton Copula

(16)

(4)Copula函數(shù)的擬合優(yōu)度評(píng)價(jià)

本研究采用3種擬合優(yōu)度評(píng)價(jià)指標(biāo)：RMSE(均方根誤差)、AIC(赤池信息量準(zhǔn)則法)、BIC(貝葉斯信息準(zhǔn)則法)，計(jì)算的評(píng)價(jià)指標(biāo)值越小，則模型模擬越好。

①RMSE：

(17)

式中：(1,2,…,)為經(jīng)驗(yàn)頻率；(1,2,…,)為理論頻率；為聯(lián)合分布函數(shù)位數(shù)；為觀測(cè)樣本長(zhǎng)度。

②AIC:

(18)

AIC=ln(MSE)+2

(19)

③BIC:

BIC=ln(MSE)+ln()

(20)

式中：為模型參數(shù)的個(gè)數(shù)。

(5)兩變量聯(lián)合概率與重現(xiàn)期。設(shè)邊緣分布函數(shù)為()=，邊緣分布函數(shù)為()=，則兩變量聯(lián)合概率分布為：

(,)=(<,<)=(,)

(21)

聯(lián)合重現(xiàn)期為：

(22)

同現(xiàn)重現(xiàn)期為：

(23)

2.2.3 經(jīng)驗(yàn)?zāi)B(tài)分解 經(jīng)驗(yàn)?zāi)B(tài)分解EMD適用于非線性、非平穩(wěn)性的信號(hào)處理，可將復(fù)雜的原始信號(hào)分解為有限個(gè)本征模函數(shù)IMF(intrinsic mode function)，分解出來的各IMF分量包含了原始信號(hào)在不同時(shí)間尺度的局部特征信號(hào)。

EMD經(jīng)驗(yàn)?zāi)B(tài)分解的步驟為：

(1)輸入需要處理的原始信號(hào)()；

(24)

(25)

(4)得到IMF1分量后，()減分量IMF1作為新的原始信號(hào)，進(jìn)行(1～3)，得到IMF2分量，以此類推，最終完成分解。

EMD分解原始序列()還得到1個(gè)趨勢(shì)項(xiàng)()：

(26)

式中：()為原始序列的各IMF分量。

3 結(jié)果與分析

3.1 水沙數(shù)量變化的階段性

由圖1可知，流域上游5，7月徑流量均于1980年和2000年發(fā)生明顯轉(zhuǎn)折，且均達(dá)到極顯著水平(<0.01)，而月輸沙量突變時(shí)間與月徑流量保持高度一致性。上游水沙3階段的變化特性，與1980s以來大面積開展水土流失綜合治理，以及1999年來退耕還林(草)措施實(shí)施的階段性事件相吻合。因此，在分析中，將水文序列分為P1(1960—1979年)、P2(1980—1999年)和P3(2000—2019年)3個(gè)階段。

注：Z表示不同月份P1(1960—1979年)、P2(1980—1999年)和P3(2000—2019年)的Mann-Kendall趨勢(shì)檢驗(yàn)統(tǒng)計(jì)量值；NS表示不顯著。

3.2 水沙聯(lián)合分布模型

假設(shè)P1、P2和P3的徑流量()和輸沙量()均服從皮爾遜III型分布，并采用檢驗(yàn)對(duì)其進(jìn)行假設(shè)檢驗(yàn)，當(dāng)顯著水平為0.05時(shí)，自由度為2的臨界值為5.99，而3個(gè)階段、的邊緣分布檢驗(yàn)結(jié)果均小于該臨界值，所以通過原檢驗(yàn)。表明各階段的、均服從皮爾遜III型分布(表2)。

表2 皮爾遜III型分布參數(shù)及檢驗(yàn)結(jié)果

本文選用Frank、Gumbel-Houggard(GH)和Clayton 3種Copula函數(shù)分別構(gòu)建劉家河站P1、P2和P3的水沙聯(lián)合分布模型。從表3可以看出，P1、P2、P3階段GH Copula函數(shù)的AIC、BIC和RMSE均小于Frank和Clayton Copula函數(shù)的計(jì)算結(jié)果。表明P1、P2和P3的水沙聯(lián)合分布模型均采用GHCopula函數(shù)擬合精度最優(yōu)。

表3 二維Copula函數(shù)參數(shù)估計(jì)

1(1960—1979年)：

2(1980—1999年)：

3(2000—2019年)：

由圖2可知，P1、P2和P3徑流量()和輸沙量()的經(jīng)驗(yàn)聯(lián)合頻率與GH Copula得到的理論聯(lián)合頻率擬合系數(shù)分別為0.999 0，0.998 1和0.993 5，表明各階段所選的Copula函數(shù)與經(jīng)驗(yàn)值擬合程度較好，選擇合理。

圖2 P1、P2和P3經(jīng)驗(yàn)累積頻率與理論累積頻率的一致性比較

3.3 水沙豐枯遭遇頻率分析

由表4可知，1960—2019年重現(xiàn)期為5～100年對(duì)應(yīng)的徑流量()和輸沙量()均呈減少趨勢(shì)。與P1相比，P2和P3各重現(xiàn)期的減小幅度分別為2.5%～26.9%和21.5%～55.7%，的減小幅度分別為31.4%～59.5%和72.7%～88.7%，的減小程度更為劇烈。造成上游和大幅降低的主要原因可以歸結(jié)為2個(gè)方面：一是P2階段的水土流失綜合治理；二是P3退耕還林(草)工程實(shí)施以來林草植被覆蓋率大幅提升(表1)。P2重現(xiàn)期為100年的徑流量相較于P1表現(xiàn)為增加趨勢(shì)，造成這一現(xiàn)象的原因可能是上游1994年遭遇了千年一遇的極端降水。

由表4還可看出，單變量的重現(xiàn)期大于兩者的聯(lián)合重現(xiàn)期(即發(fā)生大水或大沙事件)，小于兩者的同現(xiàn)重現(xiàn)期(即大水大沙同時(shí)發(fā)生的事件)。如P1時(shí)段5年一遇的和的聯(lián)合重現(xiàn)期為4.6年，同現(xiàn)重現(xiàn)期為8.1年，即較大的徑流量可能伴隨著較大的輸沙量，也可能攜帶著較小的輸沙量；徑流量和輸沙量同時(shí)超過實(shí)測(cè)最大值的事件發(fā)生概率減小。同時(shí)Coupla函數(shù)對(duì)3時(shí)期相應(yīng)徑流量和輸沙量重現(xiàn)期的估算結(jié)果表明，隨著時(shí)期進(jìn)展，重現(xiàn)期相同的情況下，5～20年重現(xiàn)期的水沙事件，其聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期大多表現(xiàn)出小幅增加趨勢(shì)，而重現(xiàn)期為50～100年的極端事件，在P1～P3階段聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期均大幅增加。表明水土流失治理和生態(tài)環(huán)境建設(shè)等下墊面條件變化，對(duì)20年重現(xiàn)期以上的水沙極端事件發(fā)生概率的影響程度，要遠(yuǎn)大于20年重現(xiàn)期以內(nèi)事件發(fā)生概率的影響。

表4 徑流量(Q)與輸沙量(QS)不同組合下的重現(xiàn)期

采用頻率法將水沙頻率分為豐、平、枯3種狀態(tài)，本文以=25%，=75%作為水沙豐枯劃分的頻率，豐枯遭遇組合：

水豐沙豐：=(≥，≥)

水豐沙平：=(≥，<<)

水豐沙枯：=(≥，≤)

水平沙豐：=(<<，≥)

水平沙平：=(<<，<<)

水平沙枯：=(<<，≤)

水枯沙豐：=(≤，≥)

水枯沙平：=(≤，<<)

水枯沙枯：=(≤，≤)

從表5可以看出：(1)P1、P2、P3的水沙豐枯遭遇頻率變化范圍分別為0.6%～41.3%，1.1%～38.6%和0.6%～28.4%。(2)P1、P2、P3的水沙豐枯遭遇頻率均表現(xiàn)為同步遭遇頻率均大于異步遭遇頻率；同步遭遇頻率中，同平表現(xiàn)最大，同豐和同枯相近；異步遭遇頻率中，水豐沙平與水平沙豐、水豐沙枯與水枯沙豐、水平沙枯與水枯沙平分別相等，水豐沙枯與水枯沙豐最小。表明北洛河上游水沙序列之間具有較強(qiáng)的相關(guān)性，出現(xiàn)極端情況的概率較小。(3)與P1相比，P2水沙豐枯各異步組合頻率均增加；同步組合中同枯頻率增加，而同豐和同平頻率均減小。(4)與P1相比，P3水豐沙枯和水枯沙豐頻率降低，其余異步組合頻率均增大；各同步組合頻率均減小。(5)總體來看，與P1相比，P2、P3水沙豐枯異步頻率均增加，同步頻率均降低，反映水土流失綜合治理和生態(tài)建設(shè)等人類活動(dòng)對(duì)流域水沙頻率一致性的影響在逐漸擴(kuò)大。

表5 不同時(shí)段水沙豐枯遭遇頻率

3.4 水沙聯(lián)合概率的多時(shí)間尺度特征

從表6和圖3可以看出，北洛河上游水沙聯(lián)合概率可分解為5或6個(gè)具有不同波動(dòng)周期的本征模態(tài)分量和1個(gè)趨勢(shì)項(xiàng)，表明北洛河上游存在復(fù)雜的水沙過程。P1具有平均周期為3.9，9.0，15.9，27.2，47.1，60.0月6個(gè)不同尺度的本征模態(tài)分量，P2具有平均周期為3.4，8.7，15.3，25.6，43.5，57.0月6個(gè)不同尺度的分量，P3則有5個(gè)本征模態(tài)分量，平均周期分別為3.3，8.1，13.8，24.8，39.7月。相同本征模態(tài)分量的平均波動(dòng)周期，隨時(shí)間推移均表現(xiàn)出減小的趨勢(shì)。

注：圖(a)為上游Copula水沙聯(lián)合概率分布，圖(b～e)為基于上游水沙聯(lián)合分布概率，使用經(jīng)驗(yàn)?zāi)B(tài)分解方法分解概率IMF1、IMF2、IMF3、IMF6/7；圖(f)為經(jīng)驗(yàn)?zāi)B(tài)分解概率以后的殘差圖(趨勢(shì)項(xiàng))。

表6 不同時(shí)段經(jīng)驗(yàn)?zāi)B(tài)分解后各IMF分量的平均周期 單位：月

IMF1表明，無論是年、月尺度，P1水沙聯(lián)合概率波動(dòng)總體較穩(wěn)定。與P1相比，P2水沙聯(lián)合概率的年內(nèi)波動(dòng)頻率略微增加，振幅呈減小趨勢(shì)；然而1994—1996年波動(dòng)頻率明顯減小，振幅顯著增加，可能是受到1994年劉家河千年一遇暴雨的影響。P3水沙聯(lián)合概率年內(nèi)呈現(xiàn)2～4次波動(dòng)；較大的波動(dòng)出現(xiàn)在2000—2002年，之后波動(dòng)逐漸放緩，趨于平穩(wěn)。

IMF2表明，P1水沙聯(lián)合概率年內(nèi)基本為1次波動(dòng)，年際間波動(dòng)總體較穩(wěn)定。P2水沙聯(lián)合概率年內(nèi)主要呈現(xiàn)1次波動(dòng)，伴隨著少數(shù)的2次年內(nèi)波動(dòng)，且振幅不明顯。P3水沙聯(lián)合概率年內(nèi)為1～2次波動(dòng)，主要表現(xiàn)在年內(nèi)豐枯月份的交替，年際間波動(dòng)較穩(wěn)定。

IMF3表明，P1水沙聯(lián)合概率主要存在1～3年的周期，到P1后期，周期明顯變短，振幅顯著增加。P2水沙聯(lián)合概率主要表現(xiàn)為1～3年的周期，振幅總體較穩(wěn)定。P3水沙聯(lián)合概率初期波動(dòng)較大，后期波動(dòng)比較穩(wěn)定。值得注意的是，如2012—2013年，振幅明顯減小，年際間變化極不明顯，可能是因?yàn)樵摃r(shí)期徑流量明顯增大，但是輸沙量卻無顯著變化，表現(xiàn)為水豐沙枯，且出現(xiàn)該現(xiàn)象的概率較小(表5)。

IMF5表明，在更大時(shí)間尺度上，水沙聯(lián)合概率呈現(xiàn)3～6年的波動(dòng)周期，且隨著時(shí)間推移，周期逐漸變短。圖3(f-1)表明上游P1水沙聯(lián)合概率總體呈上升趨勢(shì)，圖3(f-2)、(f-3)表明P2、P3水沙聯(lián)合概率總體呈下降趨勢(shì)。

4 討 論

濫砍濫伐和陡坡開墾等不合理的土地利用導(dǎo)致黃土高原生態(tài)環(huán)境惡化，侵蝕加劇。為有效治理水土流失，1970s以來，國(guó)家相繼在黃河中游開展了大規(guī)模的水土流失治理和生態(tài)建設(shè)工程。一系列的生態(tài)恢復(fù)等人類活動(dòng)使流域土壤侵蝕明顯減弱，河流輸沙量顯著降低。與P1(1960—1979年)相比，P2(1980—1999年)和P3(2000—2019年)相同重現(xiàn)期的徑流量和輸沙量均明顯減小，且輸沙量的減少程度更加顯著(表4)。表明生態(tài)恢復(fù)措施在相同重現(xiàn)期表現(xiàn)出減水更減沙的水土保持功能。

水土保持工程措施可以通過減小含沙量，進(jìn)而實(shí)現(xiàn)水土流失的治理。北洛河上游1980s開展的水土流失綜合治理，以修建淤地壩、梯田為主(表1)。其中，淤地壩的攔截效應(yīng)使得排出淤地壩的洪水含沙量顯著變小。與P1相比，P2大水伴隨大沙的頻率降低，水沙豐枯遭遇的異步頻率略有增加；徑流量和輸沙量同豐、同平的頻率明顯降低(表4和表5)。而且同時(shí)期水沙聯(lián)合概率的波動(dòng)周期降低，振幅呈減小趨勢(shì)(圖3)。表明水土保持工程措施能夠有效減小水沙聯(lián)合概率和同現(xiàn)概率，降低高徑流含沙量，調(diào)節(jié)水沙頻率耦合關(guān)系。

研究表明，植被可以在一定程度上削減洪水的洪峰流量，有效降低降水及地表徑流對(duì)土壤的侵蝕，改變水沙極端事件的發(fā)生概率。1999以來，黃土高原大力實(shí)施退耕還林(草)工程，北洛河上游植被蓋度由2001年的30%增加到2019年的60%左右。與P1相比，2000年以來重現(xiàn)期為50～100年的水沙極端事件，聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期均大幅增加，水沙豐枯遭遇的同步組合頻率均明顯降低。然而，重現(xiàn)期為5年的水沙聯(lián)合重現(xiàn)期減小，水豐沙平和水平沙枯的頻率增加，這是由于植被能顯著增加河流枯季徑流量。林草植被蓋度在小于40%～50%時(shí)，減水減沙、調(diào)節(jié)水沙極端事件發(fā)生概率的效果非常明顯，但植被蓋度達(dá)到60%左右時(shí)，產(chǎn)流、產(chǎn)沙系數(shù)均趨于穩(wěn)定。P3水沙聯(lián)合概率初期波動(dòng)較大，后期無論是波動(dòng)周期、振幅，均比較穩(wěn)定(圖3)。表明P3退耕還林(草)措施顯著影響水沙聯(lián)合概率，尤其是在2000年初期，而后期水沙聯(lián)合概率逐漸穩(wěn)定，這可能與植被蓋度持續(xù)增加有關(guān)。

流域水沙頻率耦合關(guān)系的變化，不僅受劇烈的人類活動(dòng)影響，還與極端天氣的發(fā)生息息相關(guān)。1994年劉家河發(fā)生了千年一遇的暴雨，使得P2重現(xiàn)期為100年時(shí)的徑流量相比于P1表現(xiàn)為增加趨勢(shì)(表4)，而且導(dǎo)致1994—1996年水沙聯(lián)合概率的波動(dòng)周期和振幅明顯增加(圖3)。說明極端暴雨對(duì)流域水沙聯(lián)合概率有顯著的影響，其內(nèi)在影響機(jī)制還需進(jìn)一步探討。

5 結(jié) 論

(1)Gumbel-Houggard Copula能很好地模擬上游水沙聯(lián)合概率分布。

(2)隨時(shí)段推移，5～100年重現(xiàn)期的徑流量、輸沙量均明顯減小，各設(shè)計(jì)頻率的水沙數(shù)量在發(fā)生概率上的異步性增強(qiáng)。

(3)水土流失治理和生態(tài)環(huán)境建設(shè)等下墊面條件變化，大大降低20年以上重現(xiàn)期水沙極端事件發(fā)生的概率。

(4)水沙聯(lián)合概率在不同時(shí)間尺度上具有非常復(fù)雜的變化特征。

通過Copula函數(shù)建立的水沙聯(lián)合分布模型，可以進(jìn)一步了解生態(tài)恢復(fù)對(duì)黃土高原水沙頻率耦合關(guān)系的影響，為黃土高原的水沙調(diào)控提供參考依據(jù)。