黃宇明 ,繳 健,竇希萍,郭海軍,丁 磊
(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)控等提供參考。
根據(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)行建模分析。
長(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
采用滑動(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
根據(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
根據(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
多變量水文分析中,導(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
本文基于大通水文站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ù)選取的不確定性。