張驗科，張佳新，邰雨航，紀(jì)昌明，馬秋梅

（華北電力大學(xué)水利與水電工程學(xué)院，北京 102206）

1 研究背景

入庫徑流預(yù)報誤差是水文預(yù)報不確定性的一種表現(xiàn)，目前的研究一般將量化其隨機(jī)性作為研究重點，多采用參數(shù)估計［1-2］、非參數(shù)估計［3］、最大熵［4］等方法對其分布形式近似描述，取得了大量研究成果［5-9］。考慮到入庫徑流預(yù)報誤差往往包含多種統(tǒng)計特征，若采用單一分布，例如正態(tài)分布、tLocation-Scale 分布等對其進(jìn)行擬合，某些情況下雖然可呈現(xiàn)出較好的擬合效果，但由于各分布函數(shù)有時在對稱性和尾部特征上的表現(xiàn)存在差異，并不能保證總可以找到一種分布形式剛好能夠達(dá)到預(yù)期的要求；同時，已有的研究成果多集中于單一變量或單一預(yù)見時刻的徑流預(yù)報誤差的量化分析［10-11］， 而對于多個變量或多個預(yù)見時刻的徑流過程預(yù)報誤差的分析相對較少。對于預(yù)報調(diào)度工作而言，前期的入庫徑流預(yù)報作業(yè)多在某一固定時刻以某一固定間隔對預(yù)見期內(nèi)多個后續(xù)時刻的入庫流量進(jìn)行預(yù)報［12］，而根據(jù)對應(yīng)時刻的歷時實測入庫流量序列，即可分析得到相應(yīng)的徑流過程預(yù)報誤差序列。這一預(yù)報誤差序列反映了整個調(diào)度期的預(yù)報不確定性，而僅對單一預(yù)見時刻的預(yù)報誤差進(jìn)行研究，則忽略了調(diào)度期內(nèi)不同預(yù)見時刻預(yù)報誤差間的關(guān)聯(lián)性。因此，需要在尋求單一預(yù)見時刻預(yù)報誤差分布的基礎(chǔ)上，綜合考慮分析不同預(yù)見時刻預(yù)報誤差之間的相關(guān)性，才能更為全面地揭示入庫徑流過程預(yù)報誤差的變化規(guī)律。

為解決上述入庫徑流預(yù)報誤差分布研究中存在的問題，紀(jì)昌明等［13］依據(jù)高斯混合模型（Gaussian Mixture Model，GMM）在密度函數(shù)估計中的高度自適應(yīng)性，結(jié)合高維meta-student t Copula 函數(shù)在耦合多個類型邊緣分布上的優(yōu)勢，建立了多個預(yù)見時刻入庫徑流預(yù)報誤差的GMM-Copula隨機(jī)模型，可實現(xiàn)對入庫徑流預(yù)報誤差序列的隨機(jī)模擬。在該項研究中，GMM-Copula模型是以單一預(yù)見時刻徑流預(yù)報誤差分布的擬合結(jié)果為基礎(chǔ)建立的多維徑流預(yù)報誤差聯(lián)合分布，在應(yīng)用過程中仍存在兩個關(guān)鍵問題：（1）高斯分布混合個數(shù)的確定，這是一個權(quán)衡的過程，個數(shù)過少會影響擬合的精度，相反則會導(dǎo)致模型過于復(fù)雜不利于分析；（2）參數(shù)初始值的確定，由于采用的是EM（Expectation-Maximization al?gorithm，EM）法［14］求解模型擬合需要的參數(shù)，它是一類通過迭代進(jìn)行極大似然估計的優(yōu)化算法，需要對高斯混合模型的權(quán)重等參數(shù)進(jìn)行初始化，不同的初始參數(shù)值會對迭代的結(jié)果造成較大的影響，從而決定最終的擬合精度。由于原GMM-Copula 模型對于擬合過程的兩個關(guān)鍵問題都未做細(xì)致的研究，各個預(yù)見時刻的徑流預(yù)報誤差的高斯分布混合個數(shù)以及初始參數(shù)值的確定存在較大的主觀性，難以保證擬合結(jié)果的精度，從而得到的模擬誤差序列也存在進(jìn)一步優(yōu)化的空間。

本文基于AIC（Akaike Information Criterion，AIC）與BIC（Bayesian Information Criterion，BIC）準(zhǔn)則［15］選取最優(yōu)高斯分布混合個數(shù)，同時引入數(shù)據(jù)挖掘中的K-means++算法確定高斯混合模型的初始參數(shù)值，對GMM-Copula 模型中的GMM 部分進(jìn)行改進(jìn)，以此建立基于IGMM-Copula（Improved Gauss?ian Mixture Model Copula）的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型。以雅礱江流域錦屏一級水庫短期入庫徑流過程預(yù)報誤差分析為例，基于IGMM-Copula 模型，首先對預(yù)見時刻為6 h、12 h、18 h 和24 h的入庫徑流預(yù)報誤差分布進(jìn)行擬合分析；其次，將擬合結(jié)果作為邊緣分布從而建立四維入庫徑流預(yù)報誤差的聯(lián)合分布，據(jù)此對誤差序列進(jìn)行隨機(jī)模擬與統(tǒng)計分析，并與GMM-Copula模型得到的結(jié)果進(jìn)行了對比，驗證模型的可行性與合理性。

2 GMM-Copula中高斯混合模型描述及其存在的問題

現(xiàn)假設(shè)預(yù)報誤差為ex(i)t(j)，表示當(dāng)預(yù)報起始時刻為x(i)時，對未來t(j)時刻進(jìn)行徑流預(yù)報所產(chǎn)生的預(yù)報誤差，其中i表示歷史預(yù)報的次數(shù)，i∈(1,2,…，n)；j表示預(yù)報的時刻數(shù)，j∈(1,2,…，m)?？啥x為：

式中：Sx(i)t(j)、Hx(i)t(j)分別為預(yù)報流量值和實測流量值。

當(dāng)進(jìn)行n次預(yù)報后，可得到由n個預(yù)報誤差序列組成的誤差數(shù)據(jù)矩陣E：

針對時刻t(j)的預(yù)報誤差ex(i)t(j)，可利用高斯混合模型求其分布，設(shè)模型中待估計的參數(shù)為θ，則高斯混合模型的概率密度函數(shù)可表示為：

式中：p(ex(i)t(j);θ)為高斯混合分布的概率密度函數(shù)；αk、uk、σ2k分別為第k個高斯分布的參數(shù)，分別代表權(quán)重、均值和方差，且K表示高斯混合模型中的高斯混合個數(shù)。

從式（3）中可知，若要求解單一預(yù)見時刻誤差分布擬合的表達(dá)式，則需要對參數(shù)αk、uk、σ2k進(jìn)行估計，而直接對原式進(jìn)行最大似然估計則會產(chǎn)生對數(shù)的和，求解十分困難，因此采用EM法［14］對參數(shù)進(jìn)行推求。具體的求解步驟如下。

（1）進(jìn)行EM 算法中的E 步（求期望值），提出一個隱變量γik，γik∈{0,1}，表示第i個樣本是否來自于第k個高斯分布，k=1，…，K。將γik引入式（3），得到單個樣本完全數(shù)據(jù)的對數(shù)似然函數(shù)，如式（4）。對參數(shù)θ進(jìn)行初始化，其中ak=1/K，uk設(shè)為隨機(jī)數(shù)，σ2k取1，根據(jù)利用貝葉斯定理求解γik的后驗概率P(γik=1|ex(i)t(j),θ)，見式（5）。

（2）進(jìn)行EM 算法中的M 步（求最大值），由式（4）得到全部樣本完全數(shù)據(jù)的對數(shù)似然函數(shù)，見式（6）。依據(jù)式（6）分別對αk、uk、σ2k求偏導(dǎo)，令似然函數(shù)值為0，求出更新后的參數(shù)分別如式（7）所示。

（3）不斷迭代步驟（1）和（2），重復(fù)更新αk、uk、σ2k，直到前后兩次迭代結(jié)果的變化幅度小于一個設(shè)定值ε，則終止迭代，即|θ-θnew|≤ε，θnew代表更新后的參數(shù)值，ε通常取10-4。

GMM-Copula 模型在對單一預(yù)見時刻誤差分布進(jìn)行估計的過程中，由于涉及到高斯混合個數(shù)K的選取，而通過人為確定的常規(guī)方法主觀性太強(qiáng)，尤其在于相鄰時段的預(yù)報誤差分布可能具有不同的特征，若高斯混合個數(shù)選取不當(dāng)，勢必會對最終擬合的精度造成影響；同時，模型采用隨機(jī)抽取確定初始參數(shù)值的方式同樣不嚴(yán)謹(jǐn)，初始值的不同會影響最終迭代得到的參數(shù)，從而影響最終分布擬合的結(jié)果。本文基于AIC與BIC準(zhǔn)則選取最優(yōu)的高斯混合個數(shù)，同時采用K-means++算法確定高斯混合模型參數(shù)的初始取值，建立基于IGMM-Copula 的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型。模型在優(yōu)化單一預(yù)見時刻徑流預(yù)報誤差分布擬合方法的基礎(chǔ)上，考慮了各個預(yù)見時刻之間的相關(guān)關(guān)系，并結(jié)合隨機(jī)抽樣理論，可實現(xiàn)對多個預(yù)見時刻入庫徑流過程誤差的隨機(jī)模擬。

3 基于IGMM-Copula的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型

模型主要從兩個方面敘述，提出GMM-Copula 模型改進(jìn)的兩個方面；2.3 節(jié)描述了如何利用IGMM-Copula模型建立多個預(yù)見時刻的徑流預(yù)報誤差聯(lián)合分布以及誤差序列隨機(jī)模擬的步驟。

3.1 基于AIC與BIC準(zhǔn)則的最優(yōu)高斯分布混合數(shù)選取AIC與BIC都作為衡量統(tǒng)計模型擬合優(yōu)良性的標(biāo)準(zhǔn)，與AIC不同的是，BIC引入了與模型參數(shù)個數(shù)相關(guān)的懲罰項，考慮了樣本數(shù)量，更適合樣本數(shù)量較多時的情況。對應(yīng)的計算公式如下：

式中：K為模型的參數(shù)個數(shù)，在本文中表示為高斯混合模型中的高斯混合個數(shù)；n為預(yù)報誤差樣本個數(shù)；L為模型參數(shù)求解中所對應(yīng)的似然函數(shù)。

當(dāng)模型的參數(shù)變動時，AIC 及BIC 的值越小則表示模型的擬合效果越好。由于高斯混合模型中K≥2，因此選擇以K=2作為起點，通過逐一列舉的方式計算模型擬合預(yù)報誤差分布的AIC 及BIC 的值，并選取對應(yīng)AIC及BIC值最小的K值作為當(dāng)前預(yù)見時刻誤差分布最優(yōu)的高斯分布混合個數(shù)。

3.2 基于K-means++的參數(shù)初始值確定K-means++算法是一種確定聚類迭代初始起點的算法［16］，一般用于初始化K-means++聚類算法的聚類中心。高斯混合模型本質(zhì)上作為一種聚類模型，其原理是通過計算每個數(shù)據(jù)所屬的高斯類從而實現(xiàn)聚類的效果，因此同樣可以利用K-means++算法進(jìn)行預(yù)分類，并根據(jù)預(yù)分類的結(jié)果計算高斯混合模型的初始參數(shù)值，其具體步驟如下：

（1）從單一預(yù)見時刻誤差數(shù)據(jù)集中隨機(jī)選定與最優(yōu)高斯分布混合數(shù)K相同數(shù)目的初始聚類中心；

（2）計算每個誤差數(shù)據(jù)ex(i)t(j)與之最近一個聚類中心的馬氏距離D(ex(i)t(j))；

（4）重復(fù)（2）（3），直至每一個聚類中心不再變化為止；

（5）依據(jù)得到的聚類中心將誤差數(shù)據(jù)聚類，并計算每個類的均值u′k與方差(σ2k)′作為迭代的初始值，其中權(quán)值的初始值

當(dāng)初始參數(shù)值確定后，利用EM算法迭代計算則可得到最終的αk、uk、σ2k，帶入式（3）即可得到單一預(yù)見時刻的誤差分布表達(dá)式。

3.3 多個預(yù)見時刻徑流預(yù)報誤差聯(lián)合分布建立及應(yīng)用meta-elliptic Copula 函數(shù)族中的高維meta-stu?dent t Copula［17］、高維meta-Gaussian Copula［18］是水文中常用的兩類分布，由于在預(yù)報模型在預(yù)報過程中受到的干擾因素難以量化，可能會在某一時間點產(chǎn)生與均值偏離較大的預(yù)報誤差，而高維meta-Gaussian Copula 不具有尾部相關(guān)性，因此本文在推求各個預(yù)見時刻誤差分布的前提下，選用高維me?ta-student t Copula建立多個預(yù)見時刻徑流預(yù)報誤差間的聯(lián)合分布，可表示為：

將其展開可得：

式中：C(φ1,φ2,…,φm) 為m維隨機(jī)變量聯(lián)合分布，φ1，φ2，…，φm分別表示m個預(yù)見時刻的入庫徑流預(yù)報誤差；tΣν為ν個自由度的t分布函數(shù)，其協(xié)方差矩陣為Σ；t-1ν(·)為自由度為ν的t分布的逆函數(shù)；Γ(·)為伽馬分布；X為不同預(yù)見時刻的預(yù)報誤差變量矩陣，X=(φx(i)t(1),φx(i)t(2),…，φx(i)t(m))；φ為被積函數(shù)變量矩陣，φ=[φ1,φ2,…，φm]T。

當(dāng)對整個入庫徑流預(yù)報過程進(jìn)行不確定性分析時，需要依據(jù)入庫徑流過程預(yù)報誤差聯(lián)合分布對誤差序列進(jìn)行隨機(jī)抽樣，主要分為如下兩步。

（1）依據(jù)聯(lián)合分布C(φ1,φ2,…，φm)生成隨機(jī)序列矩陣ω=[ω1,ω2,…，ωm]，ωj表示對應(yīng)預(yù)見時刻誤差分布的累積分布概率。

（2）推求對應(yīng)預(yù)見時刻誤差分布的逆累積分布函數(shù)Fj-1(φi)，將隨機(jī)序列矩陣ω帶入該函數(shù)，可得到模擬的入庫徑流過程誤差序列，表示為e′=(F1-1(ω1) ,F2-1(ω2),…，F(xiàn)m-1(ωm))。

基于IGMM-Copula的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型的研究流程如圖1所示。

圖1 基于IGMM-Copula的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型的流程

4 基于IGMM-Copula的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型

錦屏一級水電站水庫在雅礱江下游“錦官電源組（錦屏一級、錦屏二級、官地）”梯級水電站水庫群中起控制性作用，其入庫徑流預(yù)報采用的是三水源新安江模型［19］，具體是以每天早上的8 點作為預(yù)報起點，以連續(xù)預(yù)報的形式對當(dāng)天14點、20點、次日2點及次日的8點進(jìn)行徑流預(yù)報（依次對應(yīng)預(yù)見時刻為6 h、12 h、18 h和24 h），一次預(yù)報可得到未來一天的預(yù)報流量序列，以此制訂水電站水庫的調(diào)度計劃［20-21］，預(yù)報的實施過程見圖2。自投產(chǎn)以來，存在因來水預(yù)測誤差造成水庫運行水位超出控制范圍而需要調(diào)整發(fā)電計劃的情況［22］，本文選用錦屏一級水電站水庫2014年1月1日至2017年9月23日每日的實測及預(yù)報徑流資料，根據(jù)式（1）計算得到未來24 h 內(nèi)的預(yù)報誤差序列（每個序列包含4 個預(yù)見時刻的預(yù)報誤差），共得到1362 條預(yù)報誤差序列樣本，而后利用本文所提方法建立基于IGMM-Copula 的入庫徑流誤差隨機(jī)模擬模型，一方面通過求解各預(yù)見時刻的徑流預(yù)報誤差分布，分析了確定性預(yù)報模型產(chǎn)生預(yù)報誤差的統(tǒng)計特征和分布規(guī)律，同時也考慮了不同預(yù)見時刻下預(yù)報誤差的相關(guān)性，實現(xiàn)了對預(yù)報誤差的多元聯(lián)合模擬，為根據(jù)預(yù)報資料制訂水庫調(diào)度計劃，提供了更多參考信息。

圖2 錦屏一級水電站水庫徑流預(yù)報實施過程

4.1 各預(yù)見時刻徑流預(yù)報誤差分布求解基于IGMM-Copula模型中的分布擬合部分（此處稱為IGMM）對各預(yù)見時刻徑流預(yù)報誤差經(jīng)驗分布進(jìn)行擬合，首先需要確定各時刻預(yù)報誤差的最優(yōu)高斯混合數(shù)，在原GMM 模型中，由于沒有相應(yīng)的指標(biāo)作為參考，作者偏向于參數(shù)較少，更容易求解的模型結(jié)構(gòu)，因此最終選取的高斯混合數(shù)為K=2?，F(xiàn)以K=2 為起點，根據(jù)逐一列舉法計算不同K值下的AIC 及BIC的值。其計算結(jié)果如表1。

從表1可以看出，以6 h、12 h、18 h、24 h 徑流預(yù)報誤差為基礎(chǔ)建立的IGMM 在AIC 及BIC 的值取得最小時，對應(yīng)的最優(yōu)高斯混合數(shù)為K=3，據(jù)此可建立對應(yīng)權(quán)重的IGMM 對各預(yù)見時刻徑流預(yù)報誤差進(jìn)行擬合，表2為基于IGMM 的各預(yù)見時刻徑流誤差分布參數(shù)估計值。為體現(xiàn)IGMM 在誤差分布估計方面的優(yōu)勢，將原GMM 與單高斯分布（SGM）作為比較模型。其中，GMM 高斯混合個數(shù)的選取及參數(shù)初始值的確定均按照原文獻(xiàn)里介紹的方法進(jìn)行計算。

表1 不同高斯混合數(shù)下各預(yù)見時刻徑流預(yù)報誤差I(lǐng)GMM的AIC及BIC的值

表2 基于IGMM的各預(yù)見時刻徑流誤差分布參數(shù)估計值

圖3是3 種模型對不同預(yù)見時刻徑流預(yù)報誤差的擬合效果圖，其中實測預(yù)報誤差樣本個數(shù)n≥1000，據(jù)此做出的頻率直方圖總體呈現(xiàn)出對稱的形狀，誤差數(shù)據(jù)大部分集中在0值附近，SGM 擬合得到的概率密度函數(shù)較為矮胖，在峰部的集中度與頻率直方圖相差較遠(yuǎn)，而IGMM 與GMM 由于自適應(yīng)性較強(qiáng)，兩者的擬合結(jié)果在數(shù)據(jù)集中的位置更契合頻率直方圖，在圖中體現(xiàn)為概率密度函數(shù)的形狀更為尖瘦，因此能更大程度的反映實測預(yù)報誤差樣本所展現(xiàn)的規(guī)律，這也印證了高斯混合模型在概率密度估計中的優(yōu)越性。同時，由于不同預(yù)見時刻的徑流預(yù)報誤差具有不同的特征，而人為確定GMM 中高斯混合數(shù)與參數(shù)初始值難以達(dá)到最佳的擬合效果，尤其對于預(yù)報誤差數(shù)據(jù)峰部的擬合，基于IGMM 得到的擬合曲線的擬合效果明顯優(yōu)于GMM，其原因主要在于IGMM 對高斯混合數(shù)進(jìn)行了優(yōu)選，并對參數(shù)初始值進(jìn)行了預(yù)計算，因此擬合曲線更貼近于頻率直方圖。

圖3 三種模型對不同預(yù)見時刻徑流預(yù)報誤差的擬合效果

4.2 模型適用性檢驗適用性檢驗包括擬合優(yōu)度檢驗與后驗檢驗，擬合優(yōu)度檢驗的目的是用歷史數(shù)據(jù)來檢驗?zāi)Ｐ褪欠衲軌蚍从硰搅髡`差的隨機(jī)特性［23］，具體包括χ2檢驗、t-檢驗、K-S 檢驗等，其中K-S 檢驗是一種檢驗單個總體是否服從某一理論分布的非參數(shù)檢驗方法，相比于χ2檢驗?zāi)軌蚋浞值乩脴颖拘畔ⅲm用于連續(xù)和定量數(shù)據(jù)［24］。K-S檢驗統(tǒng)計量Dn定義見下式：

式中：Dn越小說明擬合優(yōu)度越高；Fn(·)為經(jīng)驗分布函數(shù)；Ft(·)為假設(shè)的理論分布函數(shù)。分別對改進(jìn)高斯混合模型、核密度估計模型以及單高斯模型進(jìn)行K-S 檢驗，顯著性水平設(shè)置為0.01，其計算結(jié)果見表3。

從表3可看出，IGMM 及GMM 都通過K-S 檢驗，即兩者與經(jīng)驗分布在顯著性水平為0.01 的前提下沒有顯著性差異，驗證了IGMM 與GMM 模型在徑流預(yù)報誤差密度估計中的可行性，同時IGMM 在各預(yù)見時刻下誤差擬合曲線的Dn值都為最小，說明其對于單一預(yù)見時刻徑流預(yù)報誤差經(jīng)驗分布的擬合優(yōu)度最高。

表3 模型K-S檢驗指標(biāo)值

后驗檢驗的目的是定量評估概率模型與數(shù)據(jù)觀測分布（即統(tǒng)計直方圖）之間的差異［25］，采用式（13）的均方根誤差εRMSE與式（14）的平均誤差百分?jǐn)?shù)εMAPE作為指標(biāo)對求得的各預(yù)見時刻徑流誤差分布進(jìn)行檢驗，得到的指標(biāo)越小，說明擬合模型與數(shù)據(jù)觀測分布之間的差異越小，其計算結(jié)果見表4。

表4 模型后驗檢驗指標(biāo)值

由表4可見，對于所有預(yù)見時刻而言，基于IGMM 得到的誤差分布的εRMSE都保持在0.06 及以下，εMAPE都保持在1.8 以下，相較于SGM 與GMM，在兩指標(biāo)上的表現(xiàn)最好，說明基于IGMM 得到的各預(yù)見時刻誤差分布與數(shù)據(jù)觀測分布之間的差異相差最小。

為分析不同預(yù)見時刻間徑流預(yù)報誤差分布的差異，圖4將IGMM 擬合得到的各徑流預(yù)報誤差分布曲線匯總，從圖中可觀察到，隨著預(yù)見時刻的增加，誤差分布形狀逐漸從“尖瘦型”變化為“矮胖型”，表面預(yù)報不確定在隨時間的增加而增大，符合確定性預(yù)報模型的預(yù)報規(guī)律。

圖4 各預(yù)見時刻徑流預(yù)報誤差分布曲線

4.3 徑流預(yù)報誤差序列隨機(jī)模擬首先分析預(yù)見期內(nèi)各時刻徑流預(yù)報誤差兩兩間的相關(guān)性，et(6)、et(12)、et(18)、et(24)分別表示6 h、12 h、18 h、24 h 的徑流預(yù)報誤差，采用Kendall 相關(guān)系數(shù)作為相關(guān)性的度量并進(jìn)行雙側(cè)顯著性檢驗，設(shè)置零假設(shè)為：變量et(6)、et(12)、et(18)、et(24)沒有相關(guān)性，顯著性水平為0.05，其計算及檢驗結(jié)果如表5所示。

從表5可知，雙側(cè)假設(shè)檢驗的P 值均小于顯著性水平，因此拒絕零假設(shè)，認(rèn)為各預(yù)見時刻的預(yù)報誤差均存在正相關(guān)性，考慮利用高維t-Copula建立徑流過程預(yù)報誤差的聯(lián)合分布并據(jù)此進(jìn)行誤差序列抽樣，實現(xiàn)對預(yù)報誤差序列的隨機(jī)模擬。

表5 各預(yù)見時刻徑流預(yù)報誤差間的Kendall相關(guān)系數(shù)

為檢驗高維t-Copula 對多維誤差的適用情況，根據(jù)圖形分析法原理，以IGMM 求得的6 h、12 h、18 h、24 h徑流預(yù)報誤差分布為例，分別按式（11）與式（15）計算誤差的理論聯(lián)合分布概率和經(jīng)驗聯(lián)合分布概率，并點繪兩者的關(guān)系圖，見圖5。

圖5 理論聯(lián)合分布概率與經(jīng)驗聯(lián)合分布概率關(guān)系

式中：Nq為同時滿足X1≤ex(1)t(j)，X2≤ex(2)t(j)，…，Xn≤ex(n)t(j)的樣本個數(shù)；N為樣本總量。

其中，由于實測預(yù)報誤差樣本個數(shù)較大，因此繪出的散點圖比較密集（圖5中藍(lán)色數(shù)據(jù)點），形狀趨向于一條線，由趨勢線的斜率可看出，圖中的點都大致落在45°對角線兩側(cè)，理論聯(lián)合分布概率和經(jīng)驗聯(lián)合分布概率的RMSE為0.00907，表明利用高維t-Copula建立多個預(yù)見時刻徑流預(yù)報誤差的聯(lián)合分布是合理可行的。

現(xiàn)根據(jù)IGMM 求出的各預(yù)見時刻徑流預(yù)報誤差分布為邊緣分布，以高維t-Copula 作為連接函數(shù)，建立基于IGMM-Copula 的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型，并與GMM-Copula 模型進(jìn)行對比，按照2.3 節(jié)的步驟分別利用兩模型對預(yù)報誤差序列進(jìn)行500 000 次隨機(jī)模擬，計算IGMM-Copula 與GMM-Copula模擬預(yù)報誤差的統(tǒng)計參數(shù)，計算結(jié)果見表6和表7。

表6 GMM-Copula模擬預(yù)報誤差與實測預(yù)報誤差的統(tǒng)計參數(shù)值

表7 IGMM-Copula模擬預(yù)報誤差與實測預(yù)報誤差的統(tǒng)計參數(shù)值

由表中可知，相較于GMM-Copula模型，IGMM-Copula模型各預(yù)見時刻模擬預(yù)報誤差的統(tǒng)計參數(shù)更貼近于實測預(yù)報誤差，而在18 h 時，由于誤差均值較小，模型模擬預(yù)報誤差與實測預(yù)報誤差的變差系數(shù)產(chǎn)生了較大的差異，而IGMM-Copula 模型明顯更優(yōu)于GMM-Copula 模型，這也反應(yīng)了不同的預(yù)報誤差邊緣分布會對抽樣結(jié)果造成一定程度的影響。從統(tǒng)計參數(shù)的整體變化趨勢分析，IGMM-Copula 模型的模擬預(yù)報誤差能與實測預(yù)報誤差保持一致，說明模型能夠很好的反映徑流預(yù)報誤差序列的統(tǒng)計特性，可作為對預(yù)報結(jié)果進(jìn)行修正的參考依據(jù)。

5 結(jié)論

本文在GMM-Copula 模型的基礎(chǔ)上，對高斯混合數(shù)的選取以及初始參數(shù)值的確定進(jìn)行了改進(jìn)，建立了基于IGMM-Copula 的入庫徑流過程預(yù)報誤差隨機(jī)模擬模型，模型可運用于單個預(yù)見時刻及多個預(yù)見時刻的徑流預(yù)報誤差分析與模擬。以錦屏一級水庫為例的結(jié)果表明：對于單一預(yù)見時刻的徑流預(yù)報誤差的描述，IGMM-Copula 模型在擬合優(yōu)度檢驗及后驗檢驗中皆優(yōu)于單高斯與GMMCopula 模型；通過IGMM-Copula 模型模擬得到的誤差序列的均值、方差和變差系數(shù)相較于GMMCopula 模型更貼近于實測誤差序列，統(tǒng)計參數(shù)的變化規(guī)律基本一致，驗證了模型的可行性，不僅為研究徑流預(yù)報誤差提供一種新的方法，也為水庫短期優(yōu)化調(diào)度計劃的制訂提供了更為豐富的參考信息。