基于Copula函數的日含沙量隨機模擬

張繼鵬，彭 楊，時玉龍，趙曉東，丁夢霞

(1.華北電力大學 可再生能源學院，北京 102206；2.華北水利水電大學 水利學院，鄭州 450045)

0 引 言

含沙量是一個主要的河流特征變量，綜合反映河流的徑流與輸沙特性。含沙量變化直接影響水庫泥沙淤積、河道沖淤、水質及水生環境，因此含沙量隨機模擬對流域水沙調控、水質與水環境管理等方面具有重要意義。

自回歸模型和解集模型是目前常用的水文隨機模型[1,2]，但這兩類模型的共同點是采用線性函數關系描述水文現象[3]，難以反映水文時間序列復雜的非線性特征。近年來非線性、非參數隨機模型成為水文隨機模擬研究的熱點，提出了非線性自回歸模型[4,5]、人工神經網絡模型[6,7]、基于核密度估計的非參數隨機模型[3,8]以及非參數解集模型[9,10]等，為水文非線性時間序列模擬提供了多種途徑。水文非線性時間序列模擬的核心問題是描述時間序列相鄰時段的非線性關系。Copula函數是一種描述變量相關關系的有效方法，由于其對變量的邊緣分布沒有限制，構造聯合分布較為靈活，目前已廣泛應用于水文多變量分析中，如洪水頻率[11,12]、降雨頻率[13,14]、干旱特征分析[15,16]、洪水遭遇[17-19]等。相比較而言，Copula函數在水文隨機模擬的應用成果不多，且主要集中在降雨[20]、徑流[21,22]和洪水[23]等方面。迄今為止，還未見到Copula函數應用于日含沙量時間序列的隨機模擬。因此，本文將采用Copula函數來描述相鄰兩日含沙量之間的相關結構，在此基礎上建立了日含沙量隨機模擬模型，并將該模型應用于長江屏山水文站日含沙量過程的隨機模擬中。

1 基于Copula函數的日含沙量隨機模擬模型

1.1 Copula函數

Copula函數是定義域為[0,1]的多維聯合分布函數。依據Sklar定理[24]，一個二維聯合分布函數可以由兩個獨立的邊緣分布函數和一個Copula連接函數組成，即

F(x,y)=Cθ[FX(x),FY(y)]=Cθ(u,v)

(1)

式中：Cθ為Copula函數；θ為Copula函數的參數，u=FX(x)，v=FY(y)分別為隨機變量X和Y的邊緣分布。

Copula函數分為橢圓型、二次型和Archimedean型三大類，其中Archimedean Copula函數結構相對簡單，計算方便，可以構造出多種形式、適應性強的組合變量聯合分布函數，在水文領域中應用最多。常見的二維Archimedean Copula主要有[24,25]：

(1)Gumbel-Hougaard (GH)Copula函數。

(2)

(2)Clayton Copula函數。

C(u,v)=(u-θ+v-θ-1)-1/θ,θ∈(0,∞)

(3)

(3)Frank Copula函數。

(4)

(4)Ali-Mikhail-Haq (AMH)Copula函數。

(5)

1.2 相鄰兩日含沙量聯合分布

根據式(1)，相鄰兩日含沙量的聯合分布可表示為：

FS(st,st-1)=Cθt[Fst(st),FSt-1(st-1)]=Cθt(ut,ut-1)

(6)

式中：st-1和st分別為第t-1日和第t日的含沙量；FS(st,st-1)為st-1和st的聯合分布；ut-1=FSt-1(st-1)和ut=FSt(st)分別為st-1和st的邊緣分布；θt為第t日的Copula函數的參數，根據相關性指標法[25]，可利用其與Kendall秩相關系數τ之間的關系求出；第t日的Kendall秩相關系數τt采用下式計算：

(7)

1.3 基于條件分布的隨機模擬模型

根據式(1)，第t日的含沙量st的條件分布可表示為：

(8)

式中：Gθt(ut,ut-1)為條件Copula函數。由此可得到基于Copula函數的日含沙量隨機模擬模型，即：

(9)

根據式(9)就可以實現對 的隨機模擬，具體步驟如下：

(1)采用正態分位數變換法對實測含沙量數據進行標準正態化；

2 實例應用

選取金沙江下游屏山站1955-2010年共56年的實測日含沙量資料作為建模的基本資料，采用所建模型對屏山站日含沙量進行隨機模擬。屏山站是金沙江干流出口控制站，流域控制面積458 592 km2，多年平均徑流量1 460 億m3，多年平均輸沙量2.39 億t，多年平均含沙量1.66 kg/m3。為簡化計算，模擬時每年2月份均取28 d。

2.1 日含沙量聯合分布建立

圖1 各日與的秩相關系數Fig.1 Daily rank correlation coefficients of

(10)

式中：Pei，Pi分別為經驗聯合頻率和理論聯合頻率；i和n分別為樣本序號和樣本容量。

各日的OLS值與對應的最優Copula函數分別如圖2和圖3所示。由圖2可見，在一年內，OLS值的最大值為0.025 1，最小值為0.008 1，平均值為0.013 7，說明所選定的二維Copula聯合分布計算出的理論聯合分布概率與經驗頻率擬合較好。圖4給出了屏山站日OLS值接近平均值時聯合分布的經驗概率與理論概率的擬合情況，從圖4中可以看出，二者擬合情況較為理想。

2.2 日含沙量的模擬及檢驗

采用上述模型對屏山站日含沙量過程進行了長序列模擬，模擬了100組與實測樣本容量(56 a)相等的日含沙量序列，并對模擬結果進行檢驗。

圖2 各日OLS值Fig.2 Daily OLS values

圖3 各日最優Copula函數Fig.3 Daily best-fit copula functions

圖4 經驗與理論聯合概率對比Fig.4 Comparison of empirical and theoretical joint probabilities

圖5 屏山站模擬與實測序列各日統計參數Fig.5 Simulated and measured statistics of daily sediment concentrations in the Pingshan station

圖5給出了屏山站實測和模擬含沙量序列各日統計參數，即均值、標準差Sd、變差系數Cv、偏態系數Cs以及滯時為1的秩相關系數τ。表1統計了屏山站日含沙量模擬序列各日統計參數絕對誤差和相對誤差的平均值以及相對誤差≤10%的通過率。由圖5和表1可見，模擬序列各日的統計參數與實測序列吻合較好，各統計參數的通過率均在98%以上，最大平均相對誤差4.22%，其中均值和秩相關系數τ的通過率均達到100%，最大相對誤差分別為0.83%和0.49%。因此，基于二維Copula函數的屏山站日含沙量過程隨機模擬結果能夠反映實測日含沙量序列的統計特征，模擬精度也較高。

表1 各統計參數通過率及平均相對誤差Tab.1 Passing rate and average relative error of statistical parameters

2.3 與其他模型比較

為進一步說明本文所建模型的模擬效果，本文還將模擬結果與分期平穩AR模型的模擬結果進行比較。表2給出這兩種模型的日含沙量模擬序列各統計參數與實測序列的相對均方根誤差RMSE[27,28]，其中RMSE采用下式計算：

(11)

式中：δo,j為實測序列第j日的統計參數；δs,j為模擬序列第j日的統計參數；m為總日數。式(11)表明，RMSE值越小，模擬序列與實測序列之間的擬合程度越好。

由表2可見，基于Copula函數的日含沙量模擬的均值和Cv值對應的RMSE值與分期平穩AR模型的基本一致，但Sd、Cs和τ值對應的RMSE值明顯小于分期平穩AR模型，表示基于Copula函數的日含沙量模擬序列在偏態性和非線性相關關系方面明顯優于分期平穩AR模型。因此，與分期平穩AR模型相比，基于Copula函數的日含沙量模擬序列的擬合結果基本一致，但在偏態性和非線性相關關系方面具有明顯的優勢。

表2 兩種模型統計參數RMSE值比較Tab.2 Comparison of RMSE values of statistical parameters in two models

3 結 論

本文考慮相鄰兩日含沙量之間的相關性，建立了基于Copula函數的日含沙量隨機模擬模型。針對日含沙量模擬時段較多，日含沙量邊緣分布難以統一確定的特點，采用正態分位數變換對實測含沙量資料進行了標準正態化處理。將該模型應用到金沙江下游屏山站日含沙量隨機模擬，并將計算結果與分期平穩AR模型計算結果進行比較，結果表明，本文所建模型能較好地保持實測日含沙量資料的統計特性，模擬精度較高，模擬序列的均值，標準差，變差系數，偏態系數以及秩相關系數τ的通過率均在98%以上，平均相對誤差最大不超過4.5%，且在偏態性和非線性相關關系方面優于分期平穩AR模型，說明本文所建模型可用于隨機模擬長序列日含沙量過程，并能為其他水文資料的長序列模擬方法提供參考。