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

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

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

1 研究背景

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

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

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

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

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

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

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

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

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

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

（1）進(jìn)行EM 算法中的E 步（求期望值），提出一個(gè)隱變量γik，γik∈{0,1}，表示第i個(gè)樣本是否來自于第k個(gè)高斯分布，k=1，…，K。將γik引入式（3），得到單個(gè)樣本完全數(shù)據(jù)的對數(shù)似然函數(shù)，如式（4）。對參數(shù)θ進(jìn)行初始化，其中ak=1/K，uk設(shè)為隨機(jī)數(shù)，σ2k取1，根據(jù)利用貝葉斯定理求解γik的后驗(yàn)概率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é)果的變化幅度小于一個(gè)設(shè)定值ε，則終止迭代，即|θ-θnew|≤ε，θnew代表更新后的參數(shù)值，ε通常取10-4。

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

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

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

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

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

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

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

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

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

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

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

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

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

將其展開可得：

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3是3 種模型對不同預(yù)見時(shí)刻徑流預(yù)報(bào)誤差的擬合效果圖，其中實(shí)測預(yù)報(bào)誤差樣本個(gè)數(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ù)的形狀更為尖瘦，因此能更大程度的反映實(shí)測預(yù)報(bào)誤差樣本所展現(xiàn)的規(guī)律，這也印證了高斯混合模型在概率密度估計(jì)中的優(yōu)越性。同時(shí)，由于不同預(yù)見時(shí)刻的徑流預(yù)報(bào)誤差具有不同的特征，而人為確定GMM 中高斯混合數(shù)與參數(shù)初始值難以達(dá)到最佳的擬合效果，尤其對于預(yù)報(bào)誤差數(shù)據(jù)峰部的擬合，基于IGMM 得到的擬合曲線的擬合效果明顯優(yōu)于GMM，其原因主要在于IGMM 對高斯混合數(shù)進(jìn)行了優(yōu)選，并對參數(shù)初始值進(jìn)行了預(yù)計(jì)算，因此擬合曲線更貼近于頻率直方圖。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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