基于Copula函數(shù)的葉爾羌河水沙豐枯遭遇研究

居金浩，彭 亮, 2，何 英, 2，娜扎凱提·托乎提,衛(wèi)仁娟

(1.新疆農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，烏魯木齊 830052； 2.新疆水利工程安全與水災(zāi)害防治重點實驗室，烏魯木齊 830052； 3.四川大學(xué) 水利水電學(xué)院 水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610065；4.四川水利職業(yè)技術(shù)學(xué)院，成都 611230)

1 研究背景

徑流與泥沙的關(guān)系涉及到水土資源規(guī)劃與管理、水利基礎(chǔ)設(shè)施建設(shè)與運行以及流域系統(tǒng)演化，一直以來都是水科學(xué)領(lǐng)域一個十分熱點的題目[1-2]。聯(lián)合國政府間氣候變化專門委員會第五次報告指出，20世紀(jì)80年代至21世紀(jì)10年代可能是北半球上千年來氣溫最高的30 a。全球氣候變化沖擊了原本平穩(wěn)的水文循環(huán)系統(tǒng)，導(dǎo)致全球范圍內(nèi)極端災(zāi)害(如干旱、洪水和高溫?zé)崂说?頻繁發(fā)生[3]，對社會經(jīng)濟、生態(tài)環(huán)境等造成重大影響，水文循環(huán)系統(tǒng)中活躍的因子徑流和輸沙也發(fā)生了明顯改變[4-5]。諸多學(xué)者針對流域內(nèi)水沙展開研究，徐金鑫等[6]分析了丹江流域徑流和泥沙單一要素在空間上的分布及在時間上的發(fā)展趨勢；也有學(xué)者分析了徑流和輸沙對氣候變暖和人類活動(水庫建設(shè)和土地利用類型改變等)的響應(yīng)[7-8]，并參照未來情景模式數(shù)據(jù)探究徑流和輸沙各自未來的發(fā)展?fàn)顩r[9-11]。以上的研究手段都將水沙割裂開來進行分析，而水沙二者之間具有一定的某種相關(guān)關(guān)系，故對年徑流量和年輸沙量兩變量分開單獨進行研究而分析所得的結(jié)論是不夠客觀的。

近年來， 在金融領(lǐng)域廣泛應(yīng)用的Copula理論廣泛應(yīng)用于水文領(lǐng)域。 Copula理論的優(yōu)勢是其不需要變量的分布嚴(yán)格相同， 不一致的邊緣分布也可以通過Copula函數(shù)構(gòu)建聯(lián)合分布模型， 而在轉(zhuǎn)換過程中不會造成變量中的信息失真[12]。 周念清等[13]采用洞庭湖流域水沙2個系列的數(shù)據(jù)， 構(gòu)建徑流量和輸沙量Copula二維聯(lián)合分布模型， 探究洞庭湖流域水沙豐枯遭遇概率， 研究結(jié)果顯示洞庭湖流域水沙運動與水沙豐枯遭遇概率有密切的聯(lián)系。 郭愛軍等[14]采用滑動相關(guān)系數(shù)法與雙累積曲線法相互對比得出水沙關(guān)系的突變點， 然后基于Copula函數(shù)建立水沙聯(lián)合分布模型， 發(fā)現(xiàn)變異后時段水沙均值均減少， 但變異前后水沙豐枯同步概率一直占主導(dǎo)地位， 變異后時段水沙豐枯異步概率增加。 諸多研究[15-17]均顯示構(gòu)建Copula水沙二維聯(lián)合分布模型能夠較好地表達水沙二者之間的相關(guān)關(guān)系。

本文基于Copula函數(shù)建立水沙二維聯(lián)合分布模型，探究葉爾羌河徑流量和輸沙量豐枯遭遇頻率，并分析水沙兩者的相互作用關(guān)系。本文研究成果可為葉爾羌河流域水沙運行調(diào)度與管理提供理論依據(jù)，并為流域的防洪減災(zāi)提供技術(shù)參考。

2 流域概況

葉爾羌河是塔里木河的重要源流之一，起源于喀喇昆侖山北坡，河源至河口全長1 078 km，流域地理位置介于74°28′E—80°54′E、34°50′N—40°31′N之間。葉爾羌河灌區(qū)處于西北干旱和半干旱氣候區(qū)，排名全國第四、新疆第一大灌區(qū)。葉爾羌河卡群站實測多年平均年徑流量67.10億m3，年徑流量變差系數(shù)0.18;多年平均懸移質(zhì)輸沙量3 028萬 t，最大日含沙量698 kg/m3(1998年6月3日)，多年平均最大月含沙量6.41 kg/m3(8月份)。

葉爾羌河上游地區(qū)氣溫較低，冰凍期較長[18]，降水稀少，造成地面植被稀疏和土壤水土保持能力不足，故河道泥沙含量較高，水沙問題突出，給灌區(qū)內(nèi)水生態(tài)修復(fù)和流域防洪減災(zāi)等工作造成很大影響。

3 資料來源與研究方法

3.1 資料來源

本文研究數(shù)據(jù)采用葉爾羌河流域卡群水文站1954—2015年、庫魯克欄干站1979—2015年徑流量和輸沙量數(shù)據(jù)，均來自喀什水文水資源勘測局。水文站具體位置見圖1。

圖1 葉爾羌河流域水系Fig.1 Water system in Yarkant River Basin and locations of hydrological stations

3.2 模比系數(shù)差積曲線法

(1)

(2)

(3)

3.3 滑動t檢驗法

(4)

根據(jù)顯著性水平α和t分布表查得臨界值tα，若|ti|>tα，則判定出現(xiàn)了變異。

3.4 滑動相關(guān)系數(shù)法

設(shè)有徑流X=(x1,x2,…,xn)和輸沙Y=(y1,y2,…,yn)2個序列，確定窗口大小為m年，時間長度為1 a的滑動，首先計算第1個窗口X1=(x1,x2,…,xm)和Y1=(y1,y2,…,ym)2個子序列的Kendall相關(guān)系數(shù)r1,然后計算第2個窗口X2=(x2,x3,…,xm+1)和Y2=(y2,y3,…,ym+1)2個子序列的Kendall相關(guān)系數(shù)r2，依次得出其后子序列的Kendall值，直至得到Xn-m+1=(xn-m+1,xn-m+2,…,xn)和Yn-m+1=(Yn-m+1,yn-m+2,…,yn)2個子序列的Kendall值rn，繪制相關(guān)系數(shù)-年份曲線。Kendall相關(guān)系數(shù)是反映X和Y兩個序列的相關(guān)程度，故繪制的相關(guān)系數(shù)-年份曲線能很好地體現(xiàn)X和Y兩個序列在不同時段的相關(guān)系數(shù)?；诖?，我們可以認(rèn)為曲線上明顯的轉(zhuǎn)折點即為X和Y的變異點。

3.5 二維Copula函數(shù)

Copula理論在1959年Sklar[19]提出之后一直局限在理論創(chuàng)新而沒有參與到各個行業(yè)的實踐應(yīng)用中。直到Embrechts等[20]首次將其應(yīng)用于金融領(lǐng)域，并取得了較好結(jié)果，隨后也在水文領(lǐng)域廣泛應(yīng)用。Copula是定義域在[0,1]上的多維聯(lián)合分布函數(shù)，通過Copula可聯(lián)接分布不一致的變量建立多元聯(lián)合分布模型。以二元Copula舉例：假設(shè)H為二維聯(lián)合分布函數(shù)，F(xiàn)和G為其邊緣分布，那么有Copula函數(shù)C使得

H(x,y)=Cθ(F(x),G(y)) 。

(5)

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

3.6 Copula的參數(shù)估算與優(yōu)度檢驗

3.6.1 Copula參數(shù)估算

本文采取非參數(shù)估計法計算Copula函數(shù)中的參數(shù)，通過θ與 Kendall 系數(shù)τ之間的關(guān)系求解。Copula參數(shù)θ和Kendall系數(shù)τ具體的對應(yīng)關(guān)系(θ對應(yīng)T,α對應(yīng)τ)見表1。

表1 函數(shù)類型與參數(shù)Table 1 Types and parameters of functions

3.6.2 Copula優(yōu)度檢驗

各類Copula都擁有各自的特點，因此需要選擇擬合度較高的Copula函數(shù)，以期提高水沙二維聯(lián)合分布模型的精確性。本文選用AIC、OLS和圖形分析等方法進行Copula模型評價。

(1)AIC又稱赤池信息量準(zhǔn)則。AIC基于熵理論而建立，可以評價所估計模型的復(fù)雜度和此模型對數(shù)據(jù)擬合的情況，AIC值較大說明模型的擬合度較差。AIC可以表示為

AIC=2K-lnL。

(6)

式中：K是參數(shù)的數(shù)量；L是似然函數(shù)。

(2)OLS又稱離差平方和最小準(zhǔn)則。通過離差平方和最小準(zhǔn)則評估Copula函數(shù)的擬合情況，選取OLS最小值作為最優(yōu)Copula。計算公式為

(7)

式中Ci和Pi分別為理論累積概率和經(jīng)驗累積概率。

(3)Genest-Rivest圖形分析法是Genest和Rivest提出的一種評價Copula函數(shù)擬合效果的方法，分別計算理論估計值Kc(t)和Ke(t)，然后點繪Kc(t)和Ke(t)關(guān)系圖，如果圖上的點都落在45°對角線上，那么表明Kc(t)和Ke(t)完全相等，即Copula函數(shù)擬合得很好。

4 結(jié)果與分析

4.1 水沙年際變化

在氣候變暖和人類活動加劇兩方面影響下，流域水沙均發(fā)生顯著變化[21-22]。葉爾羌河卡群站、庫魯克欄干站水沙多年豐枯交替情況見圖2。對比卡群站和庫魯克欄干站水沙模比系數(shù)差積曲線，可以發(fā)現(xiàn)兩站年徑流量和年輸沙量的豐枯遭遇總體上同步，均以1993年為節(jié)點，從趨于枯水(枯沙)狀態(tài)轉(zhuǎn)變?yōu)橐粋€大的豐水期(豐沙期)。

圖2 模比系數(shù)差積曲線Fig.2 Curves of product of difference of modulus coefficient

卡群站1954—1992年為年徑流量、年輸沙量下降的枯水時段，其模比差積曲線整體呈下降趨勢， 1993—2015年曲線陡升，是一個顯著的豐水(豐沙)期。庫魯克欄干站1979—1992年為年徑流量、年輸沙量下降的枯水(枯沙)時段，其模比差積曲線整體呈下降趨勢， 1993—2015年曲線陡升，是一個顯著的豐水(豐沙)期。在兩站的枯水(枯沙)期和豐水(豐沙)期，輸沙量的變化幅度遠顯著于徑流量的變化幅度。

為了驗證模比系數(shù)差積曲線判別水沙序列突變點的精準(zhǔn)度，運用滑動t檢驗法再次對兩站水沙序列的突變點進行分析。采取子序列長度為8，在α=0.05的顯著性水平下的分析結(jié)果見圖3。由圖3可見，卡群站和庫魯克欄干站水沙均在1993年發(fā)生突變，對比模比系數(shù)差積曲線分析的結(jié)果，確定1993年為兩站水沙序列的變異點。

給予大劑量LPS刺激后，Kupffer細(xì)胞胞質(zhì)組織蛋白酶B活性檢測結(jié)果顯示：分別與各自對照組（NS組）比較，WT LPS組與TLR4-/-LPS組Kupffer細(xì)胞胞質(zhì)中組織蛋白酶B活性顯著增加（P＜0.05，圖4）；WT LPS組與TLR4-/-LPS組相比，Kupffer細(xì)胞胞質(zhì)中組織蛋白酶B活性增加程度無明顯差異（圖4）。表明TLR4缺失并不影響大劑量LPS刺激下胞質(zhì)中組織蛋白酶B活性的增加。

圖3 滑動t檢驗結(jié)果Fig.3 Result of moving t-test

4.2 水沙關(guān)系變異點診斷

為確定水沙關(guān)系的突變點，采用滑動相關(guān)系數(shù)法對葉爾羌河流域卡群站的年徑流量和年輸沙量進行窗口大小為1 a、時間長度為4 a的滑動，削減異常年份對水沙產(chǎn)生的變化。圖4為基于滑動相關(guān)系數(shù)法檢驗的葉爾羌河卡群站水沙關(guān)系突變點。由圖4可知1979年、1993年是卡群站水沙相關(guān)關(guān)系較為明顯的轉(zhuǎn)變點?？ㄈ赫?979年之前平均滑動相關(guān)系數(shù)為0.70，1979—1993年平均滑動相關(guān)系數(shù)為0.53，1993年之后平均滑動相關(guān)系數(shù)為0.69。

圖4 基于滑動相關(guān)系數(shù)的水沙關(guān)系變異點判斷Fig.4 Points of abrupt changes in the relation between water and sediment by moving correlation analysis

流域水沙序列突變是由氣候和地表過程等多種因素綜合作用造成的。自20世紀(jì)90年代起，研究區(qū)年降水量呈現(xiàn)平緩增加趨勢[23]，并未有突躍式的增長，且研究區(qū)降水稀少，雨雪對徑流補給僅占多年平均徑流量的13.4%，故降水變化對徑流的突變影響較小。研究區(qū)冰川融雪對徑流的補給較大，故溫度變化對其水沙序列影響顯著[24]。有研究[25]表明葉爾羌河流域的氣溫在50多年來有明顯的上升趨勢，氣溫上升均為突變產(chǎn)生，且研究區(qū)內(nèi)多地氣溫突變均發(fā)生在20世紀(jì)90年代。也有學(xué)者指出1993年前后新疆地區(qū)的夏季零度層高度顯著上升，雪線上移坡面產(chǎn)流量增加的同時凍土退化。一方面徑流侵蝕能力增強，另一方面下墊面抗侵蝕能力減弱，最終造成水土流失加劇，河道泥沙量增加。

4.3 變化環(huán)境前后水沙關(guān)系特征對比分析

綜合考慮滑動相關(guān)系數(shù)法、雙累積曲線法的分析結(jié)論，以及人類活動加劇和氣候變暖對葉爾羌河河流域水沙關(guān)系的影響，本文確定水沙關(guān)系于1993年變異，水沙序列劃分為變異前(1954—1992年)與變異后(1993—2015年)2個階段，采用線性矩法統(tǒng)計卡群站水沙頻率分布函數(shù)的特征值，見表2。

由表2可知，年徑流量在水沙關(guān)系變異前后有10%左右的增長，年輸沙量變異后(1993—2015年)較變異前(1954—1992年)有更為明顯的增長，變幅大約為40%；年徑流量和年輸沙量的離勢系數(shù)Cv在變異前后均未發(fā)生明顯變化，表明徑流序列和輸沙序列的波動在水沙關(guān)系變異前后無明顯變化，輸沙序列值的變化程度強于徑流序列值；變異前(1954—1992年)年徑流量和年輸沙量的偏態(tài)系數(shù)Cs值較變異后(1993—2015年)變化幅度分別為-67.00%和-35.77%，表明變異后的年徑流量和年輸沙序列較變異前穩(wěn)定。

表2 水沙關(guān)系變異前后年徑流量和年輸沙量特征值與頻率分析Table 2 Characteristic values and frequencies of runoff and sediment before and after abrupt change

深入分析年徑流量和年輸沙量在不同頻率下的特征值于變異點前后的變化。水沙關(guān)系變異后(1993—2015年)較水沙關(guān)系變異前(1954—1992年)年徑流量呈現(xiàn)增長趨勢，增長幅度大約為10%；年輸沙量在1%≤P≤90%時均呈現(xiàn)增加趨勢，且增幅較大，除P=90%外，其余特征值下增幅都達到了40%以上。

4.4 變化環(huán)境前后水沙豐枯遭遇分析

4.4.1 聯(lián)合分布模型的建立

根據(jù)上文水沙關(guān)系變異前后特征值與頻率分析結(jié)論，應(yīng)用頻率法定義9種水沙組合類型,見表3。

表3 水沙組合類型Table 3 Types of runoff and sediment combinations

圖5 變異前后Kc-Ke關(guān)系Fig.5 Coherence comparisons between empirical cumulative frequency and theoretical cumulative frequency before and after variation

表4 Copula函數(shù)參數(shù)估計結(jié)果及擬合優(yōu)度檢驗Table 4 Parameters of Copula functions and evaluation result of goodness of fit

由表4可知，變異前(1954—1992年)AIC檢驗最優(yōu)的是Frank Copula，OLS檢驗最優(yōu)的是Gumbel Copula，圖形分析法檢驗最優(yōu)的是Gumbel Copula；變異后(1993—2015年)AIC檢驗最優(yōu)的是Gumbel Copula，OLS檢驗最優(yōu)的是Clayton Copula，圖形分析法檢驗最優(yōu)的是Gumbel Copula。故葉爾羌河流域變異前和變異后的水沙關(guān)系耦合模型均選擇Gumbel Copula函數(shù)，即

H(u,v)=exp{-[(-lnu)×2.581 2+(-lnv)×2.581 2]1/2.581 2} 。

(8)

式中u、v分別代表2個邊緣函數(shù)。

4.4.2 水沙豐枯遭遇結(jié)果

葉爾羌河流域的水沙豐枯遭遇結(jié)果見表5。 由表5可以看出：①水沙關(guān)系突變前后水沙豐枯同步概率均大于豐枯異步概率，變異點后水沙的豐枯同步概率更高，且同平概率最大；②流域水沙遭遇的豐枯異步概率中的水豐沙枯與水枯沙豐、水平沙枯與水枯沙平、水豐沙平與水平沙豐遭遇概率相近。水豐沙枯和水枯沙豐遭遇概率為0，表示這2種類型遭遇的概率極低;③環(huán)境變化后(1993—2015年)水沙豐枯同步概率為82.61%，相對環(huán)境變化前增加23.64%，豐枯異步概率減小，有20.62%。變化環(huán)境前后水沙同豐組合類型所占的比例變化幅度較大，增加幅度為22.74%。變異前年不同遭遇組合發(fā)生概率變化范圍為0～38.46%，變異后年不同遭遇組合發(fā)生概率變化范圍為0～43.48%。

表5 變異前后水沙豐枯遭遇概率Table 5 Probabilities of the combinations of runoff and sediment before and after abrupt change

圖6為流域變異前概率密度函數(shù)圖和概率分布函數(shù)圖，圖7為流域水沙關(guān)系變異后(1993—2015年)概率密度函數(shù)圖和概率分布函數(shù)圖。從圖中可以查出流域水沙各種遭遇組合的概率，為流域水沙治理提供決策支持。

圖6 1954—1992年水沙聯(lián)合概率分布函數(shù)和等值線Fig.6 Joint probability distribution of runoff and sediment and the isolines during 1954-1992

圖7 1993—2015年水沙聯(lián)合概率分布函數(shù)和等值線Fig.7 Joint probability distribution of runoff and sediment and the isolines during 1993-2015

5 結(jié) 論

基于葉爾羌河卡群站1954—2015年、欄桿站1979—2015年徑流量和年輸沙量，分析葉爾羌河水沙相關(guān)程度的變化情況，進一步構(gòu)建不同時段下不同類型的Copula水沙模型，分析葉爾羌河流域水沙豐枯遭遇狀況，得出以下結(jié)論：

(1)葉爾羌河年徑流量與年輸沙量均呈現(xiàn)增長趨勢，采用模比系數(shù)差積曲線圖法和滑動t檢驗法均診斷出徑流和輸沙序列分別在1993年發(fā)生變異，盡管水沙均呈現(xiàn)增加趨勢，但在突變點后輸沙量的增加趨勢遠顯著于徑流量的增加趨勢。

(2)通過滑動相關(guān)系數(shù)法判斷出葉爾羌河流域水沙關(guān)系同樣是在 1993年發(fā)生變異，氣候變暖引起融雪量的增加可能是導(dǎo)致水沙序列突變的重要因素。

(3)水沙關(guān)系變異前后的水沙聯(lián)合分布相差較大，設(shè)計概率下特征值變化較大，P≤90%時徑流量、輸沙量均增加，輸沙增加趨勢更為顯著;P>90%時，徑流量增加，輸沙量減少。

(4)兩時段中豐枯同步的概率都大于豐枯異步概率，“水豐沙枯”和“水枯沙豐”2種組合概率為0。變異前(1954—1992年)豐枯遭遇概率(38.46%)最高的類型是同平，水枯沙平(15.38%)次之；變異后(1993—2015年)遭遇概率(43.48%)最高的類型是同平，同豐(30.43%)次之，變異后的水沙豐枯同步概率比水沙關(guān)系變異前的豐枯同步概率高23.64%。