基于Copula函數(shù)的大通站水沙特征分析

黃宇明 ，繳 健，竇希萍，郭海軍，丁 磊

(1. 南京水利科學(xué)研究院 港口航道泥沙工程交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029； 2. 武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

近年來(lái)，隨著人類活動(dòng)對(duì)自然生態(tài)系統(tǒng)影響的不斷增強(qiáng)，流域下墊面情況發(fā)生顯著變化，加之全球范圍內(nèi)的極端天氣頻發(fā)，使得水文序列表現(xiàn)出多時(shí)間尺度的隨機(jī)性和突變性[1]。流域水文要素如降雨、徑流、氣溫和來(lái)沙量等發(fā)生劇烈變化，水文一致性過(guò)程遭到破壞。水、沙要素響應(yīng)過(guò)程一直是河流研究工作者關(guān)注的重點(diǎn)，河流水沙變化的一致性關(guān)系到河床沖淤、深水航道疏浚等工程實(shí)際。因此，開(kāi)展河道水沙關(guān)系突變性分析，對(duì)水庫(kù)穩(wěn)定運(yùn)行、河道規(guī)劃管理及維持變化環(huán)境下的水沙關(guān)系動(dòng)態(tài)平衡過(guò)程具有重要意義。

目前，針對(duì)單個(gè)水文要素的突變分析，最常用的方法有Mann-Kendall秩相關(guān)檢驗(yàn)法、滑動(dòng)F檢驗(yàn)法[2]、滑動(dòng)T檢驗(yàn)法[3]及Hurst系數(shù)法[4-5]等。二變量條件下的水文變異分析則一般采用滑動(dòng)偏相關(guān)系數(shù)法[6]或Copula函數(shù)理論。其中，Copula函數(shù)作為一種能夠?qū)⒍鄠€(gè)不同變量擬合起來(lái)的簡(jiǎn)單有效的實(shí)用分析工具，已被廣泛應(yīng)用于多變量水文分析研究[7]。李藝珍等[1]采用滑動(dòng)Copula函數(shù)研究了變化環(huán)境下金溝河流域徑流與氣象因素關(guān)系的變異情況；李艷玲等[8]結(jié)合滑動(dòng)窗口算法和Copula函數(shù)理論，對(duì)流域降雨徑流關(guān)系變異進(jìn)行了研究分析；郭愛(ài)軍等[9]利用Copula函數(shù)建立了水沙聯(lián)合分布，分析了水沙組合的變化特征；何兵等[10]基于滑動(dòng)Copula函數(shù)理論分析了新疆內(nèi)陸干旱河流的水文氣象要素的變化特性。

綜上可見(jiàn)，Copula函數(shù)方法是現(xiàn)階段多變量水文變異分析最常用的。基于變量序列的相關(guān)參數(shù)突變分析，將水文序列根據(jù)突變點(diǎn)劃分為不同長(zhǎng)度的子序列并分別采用Copula函數(shù)建模，能夠有效表征各子序列的不同特性，但同時(shí)也增加了各子序列邊緣分布組合方式的多樣性，而對(duì)基于Copula理論的多變量水文分析的邊緣分布不確定性分析卻缺乏足夠的研究。因此，本文采用滑動(dòng)窗口方法計(jì)算大通水文站水沙組合1965—2019年逐月序列的相關(guān)系數(shù)，分析水沙組合變異特性，基于Copula函數(shù)理論構(gòu)建水沙聯(lián)合分布模型，計(jì)算水沙同步、異步的遭遇概率，并對(duì)不同序列水沙組合的邊緣分布函數(shù)選取不確定性問(wèn)題進(jìn)行分析，以期揭示在人類活動(dòng)影響下的大通水文站水沙組合序列演化規(guī)律及突變特征，并為流域水沙規(guī)劃、灘涂治理、深水航道泥沙沖淤平衡及水資源調(diào)控等提供參考。

1 研究方法與理論

根據(jù)Sklar理論[11]，Copula函數(shù)的一般表達(dá)式可描述如下：

式中：Xi為不同的隨機(jī)變量；F(xi)=ui為不同隨機(jī)變量的邊緣分布函數(shù)；θ為Copula 函數(shù)的相關(guān)參數(shù)。表1列舉了4個(gè)不同的Copula函數(shù)類型，其中Clayton、Frank和G-H是3種常見(jiàn)的阿基米德Copula函數(shù)；Gaussian是將標(biāo)準(zhǔn)正態(tài)分布與二重積分結(jié)合的Copula函數(shù)，Ф?1是標(biāo)準(zhǔn)正態(tài)分布函數(shù)的逆函數(shù)。

表1 4種不同的Copula 函數(shù)Tab. 1 Four different Copula functions

采用歐氏距離d評(píng)價(jià)Copula函數(shù)擬合優(yōu)度，其表達(dá)式如下：

歐氏距離d越小，則Copula函數(shù)模型的擬合度越高。分別是水沙序列的理論頻率和經(jīng)驗(yàn)頻率，經(jīng)驗(yàn)頻率的計(jì)算如下：

對(duì)于常見(jiàn)的二變量水文分析如洪水峰、量組合，降雨、徑流組合或者干旱強(qiáng)度、烈度等，需要考慮的有聯(lián)合(OR)、同現(xiàn)(AND)、Kendall及生存Kendall[12]等4種不同的設(shè)計(jì)重現(xiàn)期。其中，聯(lián)合(OR)重現(xiàn)期是指二變量中至少1個(gè)變量大于其對(duì)應(yīng)設(shè)計(jì)值的概率，表達(dá)式如下：

與傳統(tǒng)的單變量水文分析不同，在多變量水文分析中，在給定設(shè)計(jì)重現(xiàn)期T或設(shè)計(jì)頻率P的情況下，同一個(gè)設(shè)計(jì)頻率對(duì)應(yīng)著無(wú)數(shù)個(gè)變量組合，最大可能組合是其中最受關(guān)注的一個(gè)。

式中：c(u,v)、f(x)和f(y)分別是Copula函數(shù)和邊緣分布函數(shù)的概率密度函數(shù)。

以大通水文站1965—2019年的55年徑流量和輸沙量資料為基礎(chǔ)樣本，進(jìn)行滑動(dòng)窗口計(jì)算，具體如下：（1）確定滑動(dòng)步長(zhǎng)L，滑動(dòng)窗口長(zhǎng)度W；（2）保持滑動(dòng)窗口W不變，從水、沙序列的第1個(gè)數(shù)據(jù)開(kāi)始，以基本步長(zhǎng)L滑動(dòng)窗口，計(jì)算每個(gè)窗口內(nèi)兩個(gè)序列的相關(guān)參數(shù)；（3）根據(jù)水沙序列的相關(guān)參數(shù)變化情況進(jìn)行水沙變異判斷；（4）采用Copula函數(shù)對(duì)每一段序列進(jìn)行建模分析。

2 研究區(qū)域及數(shù)據(jù)

長(zhǎng)江口是連接長(zhǎng)江干流與外海的重要部分，廣義的長(zhǎng)江口是指從安徽大通水文站到外海50 m等深線長(zhǎng)660 km的河段（見(jiàn)圖1）。 長(zhǎng)江口是典型的喇叭形河口，由于地球自轉(zhuǎn)偏向力的影響，長(zhǎng)江河口北岸沖刷而南岸泥沙淤積，從而形成了“三級(jí)分汊，四口入海”的河口地貌特征[13-14]。大通水文站是長(zhǎng)江口區(qū)域的重要水文站點(diǎn)，多年數(shù)據(jù)統(tǒng)計(jì)分析表明，大通站以下的長(zhǎng)江干流入流量?jī)H為大通水文站徑流量的3%，因此可以認(rèn)為大通水文站的水沙情況基本上反映了長(zhǎng)江干流的來(lái)水來(lái)沙特征[15]。本文收集了1965—2019年大通水文站逐月徑流及相應(yīng)的來(lái)沙量資料，作為實(shí)測(cè)樣本數(shù)據(jù)進(jìn)行分析計(jì)算。

圖1 長(zhǎng)江口示意圖Fig. 1 Schematic diagram of the Yangtze River estuary

3 水沙變化特征分析

3.1 水沙關(guān)系突變判別

采用滑動(dòng)窗口算法計(jì)算水沙組合的Spearman秩和Kendall秩相關(guān)系數(shù)，基本步長(zhǎng)L取12個(gè)月，6個(gè)不同的滑動(dòng)窗口W取為24、36、48、60、120 和240個(gè)月。根據(jù)相關(guān)參數(shù)的波動(dòng)情況來(lái)判斷水沙組合的變異特性，計(jì)算結(jié)果如圖2所示。從圖2可以看出，當(dāng)滑動(dòng)窗口W較?。╓=24）時(shí)，相關(guān)參數(shù)上下波動(dòng)較為頻繁，且波動(dòng)范圍也相對(duì)較大；相反，當(dāng)滑動(dòng)窗口W較大（W=240）時(shí)，相關(guān)參數(shù)波動(dòng)范圍較小。因此，在滑動(dòng)窗口取值偏小或偏大的情況下，并不能準(zhǔn)確地判斷水沙組合的突變年份。而當(dāng)滑動(dòng)窗口W取為36、48、60及120的情況下，可以發(fā)現(xiàn)水沙組合的相關(guān)參數(shù)在1979年和2000年波動(dòng)較大且具有一致性。

圖2 變化滑動(dòng)窗口下大通站徑流量與來(lái)沙量相關(guān)系數(shù)Fig. 2 Correlation coefficient of runoff and sediment discharge under the changing sliding window of Datong Station

3.2 水沙關(guān)系突變前后分析

根據(jù)水沙關(guān)系變異分析結(jié)果，以1979和2000年作為劃分點(diǎn)，將大通站水沙序列分為1965—1979年、1980—2000年和2001—2019年3個(gè)階段。采用PE3分布對(duì)分段后的徑流量序列和輸沙量序列進(jìn)行擬合，并計(jì)算3個(gè)階段重現(xiàn)期為1~100 a對(duì)應(yīng)的水、沙特征值，計(jì)算結(jié)果見(jiàn)表2。由表2可見(jiàn)，從1965到2019年，徑流量與輸沙量均呈減少趨勢(shì)；在2000年之前，大通站水沙序列表現(xiàn)出相對(duì)穩(wěn)定的狀態(tài)；2000年后，水、沙序列均表現(xiàn)出減少的趨勢(shì)，其中，來(lái)沙量減少幅度達(dá)到200%以上。造成長(zhǎng)江干流輸沙量在2000年以后大幅減少的主要原因可以歸結(jié)兩方面：一是以三峽工程為主的大型水利工程相繼建成運(yùn)行；二是流域森林植被覆蓋率大幅提升。研究表明，相較于1985—2000年，在2001—2013年間由于水庫(kù)建設(shè)使得輸沙量下降達(dá)到83.1%，而長(zhǎng)江流域植被覆蓋率上升導(dǎo)致的輸沙量下降占實(shí)際輸沙量下降的18.0%[16]。同時(shí)也可看出，1965—1979年和2001—2019年的水沙序列均有較好的變化一致性，而1980—2000年的水沙序列相較于前一個(gè)階段卻表現(xiàn)出了水增沙減的趨勢(shì)，造成這一現(xiàn)象的原因可能是由于在1980年到2000年，長(zhǎng)江流域處于水能資源開(kāi)發(fā)的高潮時(shí)期，同時(shí)流域環(huán)境破壞使得植被覆蓋率降低，且在20世紀(jì)90年代后期遭遇了極端洪水年份，破壞了這一階段水沙關(guān)系變化的協(xié)調(diào)一致性。

表2 水沙關(guān)系突變前后特征值分析Tab. 2 Analysis of eigenvalues before and after the mutation of runoff-sediment relationship

采用指數(shù)方程對(duì)3個(gè)階段的水沙序列進(jìn)行函數(shù)擬合，計(jì)算結(jié)果如圖3所示。3個(gè)階段的R2分別是0.785 6、0.745 0和0.708 5，但1980—2000年階段的擬合函數(shù)表達(dá)式過(guò)于復(fù)雜，且擬合曲線趨于直線，可能由于數(shù)據(jù)離散化太嚴(yán)重導(dǎo)致了過(guò)擬合，與根據(jù)表2計(jì)算結(jié)果分析的水沙數(shù)據(jù)不協(xié)調(diào)相一致。

圖3 不同階段大通站徑流量與來(lái)沙量關(guān)系Fig. 3 Relationship between runoff and sediment discharge in different periods of Datong Station

4 水沙聯(lián)合特征及不確定性分析

4.1 水沙聯(lián)合分布特征分析

根據(jù)水沙突變分析結(jié)果，將水沙序列分為1965—1979年、1980—2000年和2001—2019年等3段。采用PE3、 廣義帕累托(GPD)和廣義Logistic (GLD)3種分布函數(shù)分別對(duì)水沙序列進(jìn)行邊緣分布函數(shù)擬合，并采用AIC信息準(zhǔn)則和RMSE對(duì)邊緣分布函數(shù)擬合優(yōu)度進(jìn)行評(píng)價(jià)，計(jì)算結(jié)果見(jiàn)表3。根據(jù)AIC值和RMSE值越小擬合越優(yōu)的評(píng)價(jià)準(zhǔn)則，1965—1979年階段的水沙序列最優(yōu)邊緣分布是PE3-PE3，而1980—2000年和2001—2019年階段的水沙序列則是PE3-GPD擬合最優(yōu)。

表3 水沙邊緣分布選取Tab. 3 Selection of marginal distribution of runoff and sediment

表4是水沙聯(lián)合分布模型優(yōu)度評(píng)價(jià)計(jì)算結(jié)果，1965—1979 年階段的Clayton 、1980—2000年階段的Clayton及2001—2019年階段的Frank對(duì)應(yīng)的歐式距離d分別是0.045 7、0.045 6和0.052 6，為各階段最小。因此水沙序列3個(gè)階段的聯(lián)合分布模型分別由Clayton、Clayton和Frank Copula函數(shù)來(lái)擬合。

表4 水沙聯(lián)合模型選取Tab. 4 Selection of runoff-sediment joint model

根據(jù)Copula函數(shù)優(yōu)度評(píng)價(jià)結(jié)果，3個(gè)階段水沙組合的聯(lián)合模型為式（6）~（8），由式（4）計(jì)算水沙組合不同序列的風(fēng)險(xiǎn)模型及聯(lián)合(OR)重現(xiàn)期等值線見(jiàn)圖4。圖4為3個(gè)階段在重現(xiàn)期水平100、50、20、10、5和2 a聯(lián)合重現(xiàn)期等值線，同一個(gè)重現(xiàn)期對(duì)應(yīng)著無(wú)數(shù)個(gè)水沙組合。不同重現(xiàn)期等值線分布較為規(guī)律，說(shuō)明各階段的水沙組合并未出現(xiàn)較大的波動(dòng)。

圖4 3個(gè)階段水沙聯(lián)合分布及聯(lián)合重現(xiàn)期等值線Fig. 4 Joint runoff-sediment distribution and joint recurrence contours in three periods

對(duì)于1965—1979年，有：

對(duì)于1980—2000年，有：

對(duì)于2001—2019年，有：

二變量水文分析中，各變量之間的變化是否具有一致性，反映了變量之間的相關(guān)關(guān)系。本文以pf=25%、pk=75%作為劃分頻率，將水沙概率分為豐、平、枯3種狀態(tài)，9種不同組合方式[17]，對(duì)水沙組合的變化特征進(jìn)行分析，計(jì)算結(jié)果如表5所示。從整體上看，3個(gè)不同序列的水沙同步頻率分別為84.52%、84.40%和83.20%，要遠(yuǎn)大于水沙異步的頻率，說(shuō)明大通水文站水沙序列之間具有較強(qiáng)的相關(guān)性，但水沙同步概率在3個(gè)階段依次減小也反映了人類活動(dòng)對(duì)流域水沙一致性的影響在逐步擴(kuò)大。在水沙同步頻率中，水平沙平的占比較高（分別為42.26%、42.20%和41.60%），而水豐沙豐、水枯沙枯的比例則較為接近。其次，在水沙異步頻率中，水豐沙枯（或水枯沙豐）的概率均為0，說(shuō)明來(lái)水來(lái)沙條件維持在相對(duì)穩(wěn)定的平衡狀態(tài)，極端情況較少。

表5 大通水文站水沙豐枯遭遇計(jì)算Tab. 5 Calculation of flood and drought of runoff and sediment at Datong Hydrological Station

4.2 邊緣函數(shù)選取不確定性分析

多變量水文分析中，導(dǎo)致結(jié)果不確定性的主要原因有兩方面：一是由于樣本實(shí)測(cè)數(shù)據(jù)長(zhǎng)度不夠，導(dǎo)致樣本容量代表總體時(shí)出現(xiàn)偏差；二是模型函數(shù)的選取及函數(shù)參數(shù)的確定[18]。本文開(kāi)展了針對(duì)邊緣分布函數(shù)選取的不確定性分析，分別計(jì)算了聯(lián)合重現(xiàn)期為20 a條件下每個(gè)序列不同邊緣分布對(duì)應(yīng)的1 000組最可能組合，進(jìn)行邊緣函數(shù)選取不確定性分析，結(jié)果如圖5和表6所示。

圖5 3個(gè)階段不同邊緣分布組合對(duì)應(yīng)的50%、75%和95%置信區(qū)間面積Fig. 5 Areas of 50%, 75%, and 95% confidence intervals for different combinations of marginal distributions in three periods

由圖5可以看出1965—1979年和1980—2000年階段的PE3-PE3組合方式對(duì)應(yīng)的最可能組合相對(duì)集中，因此置信區(qū)域呈細(xì)長(zhǎng)型，其他4個(gè)則呈橢圓形，6種情況下的95%置信區(qū)間均跨越了各自對(duì)應(yīng)的10 a和20 a聯(lián)合重現(xiàn)期等值線。從表6計(jì)算結(jié)果可以發(fā)現(xiàn)，不同的邊緣分布組合條件下，Copula函數(shù)參數(shù)θ的變化幅度不大；對(duì)比3個(gè)序列的95%置信區(qū)間面積可以發(fā)現(xiàn)，采用PE3-PE3組合方式對(duì)應(yīng)的置信區(qū)間面積明顯小于PE3-GPD組合的置信區(qū)間面積，說(shuō)明在本文中PE3- PE3邊緣函數(shù)組合方式可以減少函數(shù)選取的不確定性。

表6 邊緣分布選取不確定性評(píng)價(jià)指標(biāo)的95%置信區(qū)間取值Tab. 6 Values of 95% confidence interval of uncertainty evaluation index selected by marginal distribution

5 結(jié) 語(yǔ)

本文基于大通水文站1965—2019年的逐月徑流量與來(lái)沙量資料，采用滑動(dòng)Copula函數(shù)對(duì)水沙聯(lián)合序列進(jìn)行變異分析，通過(guò)對(duì)分段序列進(jìn)行Copula函數(shù)建模，分析水沙變異的一致性，并對(duì)各段序列的邊緣函數(shù)選取進(jìn)行不確定性分析研究，得出主要結(jié)論如下：

（1）大通水文站多年水沙序列關(guān)系變異發(fā)生在1979年和2000年。1965—2019年間徑流量和輸沙量均呈減小趨勢(shì)；相比1965—1979年階段，1980—2000年間水沙變化波動(dòng)出現(xiàn)不一致現(xiàn)象；2000年后徑流量和輸沙量減少趨勢(shì)較為明顯，其中輸沙量減少達(dá)200%以上。

（2）1965—1979年、1980—2000年和2001—2019年3段序列的最優(yōu)Copula函數(shù)模型分別是Clayton、Clayton和Frank Copula。水沙豐枯一致性計(jì)算結(jié)果表明，大通水文站多年水沙序列同步頻率分別是84.52%、84.40%和83.20%，遠(yuǎn)大于異步頻率，反映了大通站水沙序列具有較強(qiáng)的一致性。

（3）邊緣分布函數(shù)選取不確定性分析表明，采用PE3-PE3邊緣函數(shù)的組合方式的95%置信區(qū)間遠(yuǎn)小于采用PE3-GPD組合的置信區(qū)間面積，說(shuō)明在本文中PE3-PE3邊緣函數(shù)組合方式可以減少函數(shù)選取的不確定性。