周研來(lái),梅亞?wèn)|,張代青

(武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,湖北 武漢 430072)

水文隨機(jī)模擬[1-2],是指根據(jù)水文系統(tǒng)觀測(cè)資料的統(tǒng)計(jì)和隨機(jī)變化規(guī)律,建立能預(yù)估系統(tǒng)未來(lái)水文情勢(shì)的隨機(jī)模型,由隨機(jī)模型通過(guò)統(tǒng)計(jì)試驗(yàn)獲得大量的模擬序列,再進(jìn)行水文系統(tǒng)分析計(jì)算,解決系統(tǒng)的規(guī)劃、設(shè)計(jì)、運(yùn)行與管理問(wèn)題的方法。水文隨機(jī)模擬技術(shù)的關(guān)鍵是建立合理可靠的隨機(jī)水文模型。目前,國(guó)內(nèi)外研究的隨機(jī)水文模型主要包括回歸類模型[3]、解集類模型[4]、非參數(shù)類模型[5-6]和非線性類模型[7-9]等。Sharma等[5]用核密度估計(jì)方法估計(jì)序列的條件分布,在核函數(shù)選為高斯函數(shù)時(shí)導(dǎo)出了基于非參數(shù)方法進(jìn)行隨機(jī)模擬的NP(1)模型。王文圣等[6]又將NP(1)模型推廣到單變量和多變量的多階情形,即NP(p)模型,并研究了模型在徑流模擬中的可行性。肖義等[10]和張濤等[11]就洪水特征變量間的相關(guān)性基于Copula函數(shù)提出了一種洪水過(guò)程隨機(jī)模擬方法。閆寶偉等[9]建立了基于Copula函數(shù)的一階季節(jié)性隨機(jī)模型,嘗試將其應(yīng)用于徑流的隨機(jī)模擬,所建模型能較好地保持原序列的偏態(tài)特性、非線性相關(guān)性和概率密度特征等統(tǒng)計(jì)特性。一階自回歸模型[3]和典型解集模型[4]由于構(gòu)造簡(jiǎn)單,在水文模擬中備受青睞,為避開(kāi)一階自回歸模型和典型解集模型中對(duì)水文序列的正態(tài)轉(zhuǎn)化(該過(guò)程難免會(huì)出現(xiàn)部分信息失真),筆者嘗試用Copula函數(shù)構(gòu)造一種新的徑流過(guò)程解集隨機(jī)模型,即先用Copula函數(shù)構(gòu)造年徑流量與最大月徑流量之間的聯(lián)合分布,通過(guò)聯(lián)合分布達(dá)到隨機(jī)模擬年徑流量的目的,再優(yōu)選實(shí)測(cè)徑流過(guò)程,將年徑流量分解為各分量的徑流過(guò)程。

1 基于Copula函數(shù)的解集隨機(jī)模型

1.1 Copula函數(shù)的構(gòu)造

Copula函數(shù)是定義域?yàn)閇0,1]均勻分布的多維聯(lián)合分布函數(shù),由Sklar定理可以將多個(gè)隨機(jī)變量的分布函數(shù)連接起來(lái)構(gòu)造聯(lián)合分布,對(duì)于二維情況,可以表述為[8,12]

式中:H(x,y)為隨機(jī)變量X和Y的聯(lián)合分布函數(shù);Cθ為參數(shù)為θ的Copula函數(shù);FX(x),FY(y)分別為隨機(jī)變量X和Y的邊緣分布函數(shù)。

Copula函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展,文獻(xiàn)[12]已有詳細(xì)介紹。

徑流過(guò)程的2個(gè)特征量(最大月徑流量和年徑流量)具有一定的相關(guān)性。兩變量耿貝爾聯(lián)合(Gumbel-Hougaard Copula)函數(shù)可用于描述具有上尾相關(guān)特性多變量之間的相關(guān)關(guān)系[13]。例如它可以很好地描述當(dāng)一個(gè)變量增幅較大時(shí),另一個(gè)與之存在相關(guān)性的變量隨之明顯增大的多變量相關(guān)關(guān)系[10]?；贕umbel-Hougaard Copula函數(shù)的數(shù)學(xué)特性,采用該函數(shù)來(lái)描述最大月徑流量和年徑流量,其數(shù)學(xué)表達(dá)式為[8,10,13]

式中:FX(x)和FY(y)均為PⅢ型分布函數(shù);θ越大,說(shuō)明u,v的相關(guān)性越強(qiáng),θ的估計(jì)見(jiàn)文獻(xiàn)[9]。

1.2 模型算法

許多水文要素是累加而成的,如年水量是由月水量累加而得,15 d洪水總量是由15 d的日洪水量累計(jì)而得。一般而言,將各分量累加起來(lái)的總量叫做聚集;反之,從總量中可以分解出分量。解集隨機(jī)模型的實(shí)質(zhì)在于將總量隨機(jī)解集成各分量,其顯著的特點(diǎn)在于保持水量平衡,即各分量的水量相加嚴(yán)格等于總水量。Copula函數(shù)對(duì)最大月徑流量 X和年徑流量Y進(jìn)行聯(lián)合描述,通過(guò)聯(lián)合分布的隨機(jī)抽樣方法對(duì)存在相關(guān)性的最大月徑流量和年徑流量成對(duì)地取樣,從而達(dá)到隨機(jī)模擬最大月徑流量和年徑流量的目的。參照文獻(xiàn)[12]的峰量隨機(jī)抽樣方法以及典型解集模型[14]得出基于Copula函數(shù)的解集隨機(jī)模型模擬徑流過(guò)程的步驟如下:

步驟1 根據(jù)最大月徑流量和年徑流量的聯(lián)合分布可以得到當(dāng)最大月徑流量為某指定值時(shí)年徑流量的條件分布:

步驟2 產(chǎn)生服從[0,1]均勻分布的2個(gè)獨(dú)立隨機(jī)數(shù)r1和r2。

步驟3 令r1為最大月徑流量發(fā)生的不超過(guò)x值的概率FX(x),即 u=FX(x)=r1,根據(jù) x=(u)得到x值。

步驟4 令r2為最大月徑流量X=x時(shí)年徑流量y的條件概率分布值,即r2=CX(y X =x),從而可根據(jù)y=(r2X=x)計(jì)算得到 y。

步驟5 由步驟3和步驟4即可抽取1對(duì)最大月徑流量和年徑流量的組合(x,y)。

步驟6 以模擬的最大月徑流量和年徑流量組合(x,y)為控制條件,選擇最大月徑流量和年徑流量的比值與x/y值接近的實(shí)測(cè)徑流過(guò)程,并采用解集[11]方法進(jìn)行分解,得到模擬的徑流過(guò)程:

式中:Xi(t),xi(t)分別為第i年第t個(gè)截口的模擬徑流值和實(shí)測(cè)徑流值;Yi為第i年的模擬年徑流量;A為將年徑流量(億m3)轉(zhuǎn)換為總流量(m3/s)的單位換算系數(shù);Ki(t)為第i年第t個(gè)截口的分解系數(shù)。

步驟7 重復(fù)步驟2至步驟6共M次,即可模擬出M年的徑流過(guò)程。

2 模型應(yīng)用研究

利用黃河中上游河口鎮(zhèn)水文站1919—1988年共70年的實(shí)測(cè)月徑流量資料,檢驗(yàn)基于Copula函數(shù)的解集隨機(jī)模型的可行性和有效性。由于所選資料序列較長(zhǎng),某種意義上可以代表整體的統(tǒng)計(jì)特性,以此來(lái)驗(yàn)證模型,可信度較高。同時(shí)選取季節(jié)性一階自回歸模型和典型解集模型作為對(duì)照,以說(shuō)明本文所建模型的實(shí)用性。

在實(shí)測(cè)徑流資料中,每年選擇1對(duì)最大月徑流量和年徑流量數(shù)據(jù),構(gòu)成樣本容量為70的聯(lián)合觀測(cè)值系列,并據(jù)此建立最大月徑流量X和年徑流量Y的兩變量聯(lián)合分布。X和Y的邊緣分布采用PⅢ型分布,分別為FX(x)和 FY(y),采用線性矩法估計(jì)FX(x)和FY(y)的參數(shù),得到實(shí)測(cè)年徑流量的統(tǒng)計(jì)特征值,見(jiàn)表1。將估計(jì)的參數(shù)假定為河口鎮(zhèn)站最大月徑流量和年徑流量總體分布的參數(shù)。依據(jù)Kendall秩相關(guān)系數(shù)和Copula函數(shù)的參數(shù) θ的解析關(guān)系[8,12],得到最大月徑流量和年徑流量的參數(shù)θ=4.25。圖1為聯(lián)合觀測(cè)值的經(jīng)驗(yàn)分布和理論分布的相關(guān)關(guān)系。據(jù)圖1可知,θ=4.25時(shí),最大月徑流量和年徑流量聯(lián)合觀測(cè)值的經(jīng)驗(yàn)分布和理論分布擬合較好,θ=4.25可作為聯(lián)合總體分布的參數(shù)估計(jì)量,所建立的X和Y的聯(lián)合分布是合理的。由Sklar定理及式(2)得到X和Y的聯(lián)合分布F(x,y):

?

圖1 聯(lián)合觀測(cè)值經(jīng)驗(yàn)分布和理論分布的比較

繪制最大月徑流量和年徑流量的兩變量聯(lián)合分布函數(shù)及概率密度函數(shù)三維圖,如圖2、圖3所示。

圖2 X和Y的聯(lián)合分布函數(shù)

圖3 X和Y的聯(lián)合概率函數(shù)

此外,從聯(lián)合分布隨機(jī)抽取7萬(wàn)對(duì)最大月徑流量和年徑流量數(shù)據(jù),點(diǎn)繪得到它們的頻率散點(diǎn)圖,如圖4所示。由圖4可見(jiàn),當(dāng)模擬的最大月徑流量為某一頻率時(shí),可能發(fā)生多種頻率的年徑流量,最大月徑流量和年徑流量之間不存在嚴(yán)格的一一對(duì)應(yīng)關(guān)系;從回歸關(guān)系和相關(guān)關(guān)系來(lái)看,最大月徑流量和年徑流量頻率的回歸系數(shù)為0.92,相關(guān)系數(shù)為0.85,說(shuō)明最大月徑流量和年徑流量組合既存在較好的同頻率特性又存在一定的隨機(jī)性;從頻率分布狀況來(lái)看,最大月徑流量和年徑流量偏小值的頻率居多,中等值的頻率次之,而較大值頻率分布較少,這與水文要素的枯、平、豐水年的頻率分布狀況十分吻合。

圖4 隨機(jī)模擬X和Y頻率的散點(diǎn)圖

根據(jù)圖4中任意一點(diǎn)對(duì)應(yīng)的最大月徑流量和年徑流量,可求得最大月徑流量與年徑流量的比值,依據(jù)該比值選擇具有相似最大月徑流量與年徑流量之比的實(shí)測(cè)徑流過(guò)程,采用典型解集方法進(jìn)行分解,模擬多年月徑流過(guò)程(如7萬(wàn)年),其年徑流量的統(tǒng)計(jì)參數(shù)見(jiàn)表1。

統(tǒng)計(jì)分析河口鎮(zhèn)水文站1919—1988年實(shí)測(cè)的月徑流過(guò)程資料,可知1967年為豐水年(最大月徑流量為3401m3/s,年徑流量為420億m3),見(jiàn)圖5。圖5還給出了以1967年實(shí)測(cè)月徑流過(guò)程為基礎(chǔ)進(jìn)行縮放后得到的2條月徑流過(guò)程線。由圖5可見(jiàn),模擬的徑流過(guò)程在形狀上與實(shí)測(cè)過(guò)程十分相似。因此,隨機(jī)模擬的徑流過(guò)程系列能夠反映實(shí)測(cè)徑流的形狀。

圖5 以1967年實(shí)測(cè)徑流過(guò)程模擬的徑流過(guò)程線

對(duì)基于Copula函數(shù)的解集隨機(jī)模型是否能很好地保持實(shí)測(cè)系列的年徑流量統(tǒng)計(jì)特征進(jìn)行檢驗(yàn),并與水文領(lǐng)域常用的季節(jié)性一階自回歸模型和典型解集模型[14]進(jìn)行比較。對(duì)于每個(gè)模型,在其模擬的徑流過(guò)程系列對(duì)應(yīng)取出7萬(wàn)個(gè)年徑流量,分別組成1000個(gè)樣本容量為70的年徑流量樣本。對(duì)每個(gè)年徑流量樣本分別估計(jì)統(tǒng)計(jì)特征參數(shù)(均值,Cv和Cs),并分別計(jì)算統(tǒng)計(jì)特征參數(shù)的均值。采用綜合相對(duì)均方誤差RMSE[9-10]作為模型評(píng)價(jià)指標(biāo):

評(píng)價(jià)指標(biāo)RMSE的值越小,說(shuō)明隨機(jī)模型越能夠保持實(shí)測(cè)系列的統(tǒng)計(jì)特征。表1給出了實(shí)測(cè)系列和模擬系列的年徑流量統(tǒng)計(jì)特征值及評(píng)價(jià)指標(biāo)值。由表1可知,季節(jié)性一階自回歸模型模擬徑流過(guò)程的年徑流量統(tǒng)計(jì)特征效果相對(duì)較差;基于Copula函數(shù)的解集隨機(jī)模型的RMSE均低于季節(jié)性一階自回歸模型和典型解集模型的RMSE。究其原因,基于Copula函數(shù)的解集隨機(jī)模型避開(kāi)了季節(jié)性一階自回歸模型和典型解集模型中對(duì)水文序列的正態(tài)轉(zhuǎn)化(難免會(huì)出現(xiàn)部分信息失真),而是直接隨機(jī)模擬年徑流總量,然后選用實(shí)測(cè)徑流過(guò)程分解為模擬的徑流過(guò)程,因此能夠很好地保持實(shí)測(cè)系列的徑流過(guò)程形狀及年徑流量的統(tǒng)計(jì)特征。

3 結(jié) 語(yǔ)

構(gòu)造了最大月徑流量和年徑流量的聯(lián)合分布,采用基于聯(lián)合分布的隨機(jī)抽樣方法模擬了大樣本的徑流過(guò)程,從而提出了一種新的徑流過(guò)程隨機(jī)模擬方法,建立了基于Copula函數(shù)的解集隨機(jī)模型。該模型避開(kāi)了季節(jié)性一階自回歸模型和典型解集模型中對(duì)水文序列的正態(tài)轉(zhuǎn)化和相依形式的主觀假定,很好地考慮了最大月徑流量和年徑流量之間的相關(guān)性及兩者頻率組合的隨機(jī)性,選用實(shí)測(cè)徑流過(guò)程作為隨機(jī)模擬徑流過(guò)程的時(shí)程分配形式,能很好地保持實(shí)測(cè)系列的徑流過(guò)程形狀及年徑流量的統(tǒng)計(jì)特征。實(shí)例計(jì)算結(jié)果表明,基于Copula函數(shù)的解集隨機(jī)模型的模擬效果至少可以與常用的典型解集模型相當(dāng),適用于水文時(shí)間序列的隨機(jī)模擬,從而為水文水資源的隨機(jī)模擬提供了一種新途徑。

[1]芮孝芳,劉方貴,邢貞相.水文學(xué)的發(fā)展及其所面臨的若干前沿科學(xué)問(wèn)題[J].水利水電科技進(jìn)展,2007,27(1):75-79.

[2]胡慶芳,王銀堂,劉克琳,等.基于改進(jìn)的兩參數(shù)月水量平衡模型的月徑流模擬[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2007,35(6):638-642.

[3]葛朝霞,王會(huì)容.黃河上游月徑流量與海溫的關(guān)系研究[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2003,31(6):605-608.

[4]KOUTSOYIANNIS D,MANETAS A.Simple disaggregation accurate adjusting procedures[J].Water Resources Research,1996,32(7):2105-2117.

[5]SHARMA A,TARBOTON D G,LALL U.Streamflow simulation:a nonparametric approach[J].Water Resources Research,1997,33(2):192-308.

[6]王文圣,丁晶.基于核密度估計(jì)的多變量非參數(shù)隨機(jī)模型初步研究[J].水利學(xué)報(bào),2003,34(2):9-14.

[7]PRAIRIE J R,RAJAGOPALAN B,FULP T J,et al.Modified K-NN modelfor stochastic streamflow simulation[J].Journal of Hydrologic Engineering,2006,11(4):371-378.

[8]SALVADORI G,de MICHLELE C.Frequency analysis via copulas theoretical aspects and applications to hydrological events[J].Water Resources Research,2004,40(12):1-17.

[9]閆寶偉,郭生練,劉攀,等.基于Copula函數(shù)的徑流隨機(jī)模擬[J].四川大學(xué)學(xué)報(bào):工程科學(xué)版,2010,42(1):5-9.

[10]肖義,郭生練,熊立華,等.一種新的洪水過(guò)程隨機(jī)模擬方法研究[J].四川大學(xué)學(xué)報(bào):工程科學(xué)版,2007,39(2):55-60.

[11]張濤,趙春偉,雒文生.基于Copula函數(shù)的洪水過(guò)程隨機(jī)模擬[J].武漢大學(xué)學(xué)報(bào):工學(xué)版,2008,41(4):1-4.

[12]郭生練,閆寶偉,肖義,等.Copula函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展[J].水文,2008,28(3):1-7.

[13]NELSON R B.An introduction to Copula[M].New York:Springer,1999.

[14]中華人民共和國(guó)能源部,中華人民共和國(guó)水利部.水利水電工程設(shè)計(jì)洪水計(jì)算手冊(cè)[M].北京:中國(guó)水利水電出版社,2001.