基于CEEMD-TD模型的月徑流隨機模擬

王倩麗，馬細(xì)霞，2，趙 璐，張靜文，程 旭

（1.鄭州大學(xué)水利科學(xué)與工程學(xué)院，鄭州450001；2.鄭州大學(xué)黃河生態(tài)保護與區(qū)域協(xié)調(diào)發(fā)展研究院，鄭州450001；3.鄭州航空港興港投資集團有限公司，鄭州450001；4.商丘工學(xué)院土木工程學(xué)院，河南商丘476000）

0 引言

水文隨機模型是依據(jù)水文過程觀測資料，分析水文過程隨機變化特性的一種數(shù)學(xué)模型，其隨機模擬的大量水文序列對于水資源評價、水庫調(diào)度方案編制、各種水利參數(shù)確定等具有重要應(yīng)用價值。解集類模型形象、簡單、充分利用樣本信息，能夠進行多季節(jié)和多站點的水文序列模擬，在多個時間尺度或空間尺度上總量與分量序列的主要統(tǒng)計特性如均值、方差和協(xié)方差結(jié)構(gòu)等能夠得以保持，并且遵循了水量平衡［1］，因此被廣泛地應(yīng)用于水文模擬當(dāng)中。徐利崗［2］等運用季節(jié)性隨機模型的理論方法建立洛河狀頭站典型解集（TD）模型，經(jīng)實用性檢驗，月徑流模擬序列與實測序列各項統(tǒng)計參數(shù)吻合良好，符合精度要求。梅超［3］等以平寨水庫入庫徑流隨機模擬為例，分別采用雙層模型法和典型解集法分解模擬月徑流，其模擬結(jié)果顯示典型解集分解法較雙層模型分解法能更好地保持月徑流的統(tǒng)計特性。TD模型具有直觀、方便、能充分利用實測樣本資料信息等優(yōu)點，且解集而得的序列沒有引入任何與概率分布有關(guān)的假定，可以全面反映樣本序列的統(tǒng)計特性，但其在模擬過程中僅適用于平穩(wěn)序列的隨機模擬，具有一定的局限性［4］。

Huang 等［5］首次提出經(jīng)驗?zāi)B(tài)分解（EMD）法，該方法認(rèn)為任何一個復(fù)雜的時間序列都可由EMD 分解為有限數(shù)量的本征模態(tài)函數(shù)（IMF），分解得到的原始信號與EMD 的近似正交基是隨信號的變化而變化的，因而使得EMD 具備了更強的自適應(yīng)性，更有利于非平穩(wěn)數(shù)據(jù)處理［6］；Wu［7］等在對白噪聲進行EMD分解深入研究的基礎(chǔ)上，提出了集合經(jīng)驗?zāi)B(tài)分解（EEMD）法，以期遏制EMD 分解的模態(tài)混疊現(xiàn)象；EEMD 分解算法雖然可以通過增加白噪聲來解決頻率混疊問題并提高信號分解能力和算法穩(wěn)定性，但是不能完全降低白噪聲，因此在EEMD 算法基礎(chǔ)上，CEEMD 算法通過在信號中增加分解次數(shù)來降低對符號的相反噪聲，有效解決白噪聲對誤差精度的影響問題［8］；王燕鵬等基于CEEMD 具有非平穩(wěn)信號平穩(wěn)化的能力構(gòu)建了CEEMDBP 灌區(qū)地下水埋深預(yù)測耦合模型，結(jié)果表明CEEMD-BP 耦合模型和其他模型相比具有較好的預(yù)測效果［9］。

本文以黃河流域的3 個典型水文站為例，將對非平穩(wěn)序列有較強處理能力的互補式集合經(jīng)驗?zāi)B(tài)分解（CEEMD）法引入到徑流隨機模擬中，提出一種新的互補式集合經(jīng)驗?zāi)B(tài)分解-典型解集（CEEMD-TD）模型，通過檢驗截口的主要統(tǒng)計參數(shù)來評價隨機模型的模擬效果，以驗證該隨機模型的可靠性。

1 CEEMD-TD模型介紹

TD模擬月徑流的思路是：先對剔除確定性成分后的各站年徑流隨機序列進行平穩(wěn)性檢驗，然后對平穩(wěn)的實測的年徑流量建立隨機模型，采用自回歸AR（p）模型進行隨機模擬，得到模擬的年徑流序列。在各站年徑流實測系列中分別找出與每個年徑流模擬量最為接近的那個值，以此為典型年，按照該年各月徑流量占全年徑流量的分配系數(shù)，對模擬的年徑流量進行縮放，得到模擬的月徑流序列。

CEEMD-TD 模型是對TD 模型的改進，其過程仍分為兩大步，即年徑流自回歸隨機模擬和年徑流月分配過程的確定。具體思路是：①建立年徑流CEEMD-AR 隨機模型，首先剔除年徑流實測序列中非周期成分與周期成分，剔除確定性成分后的隨機序列進行CEEMD 分解，得到一組固有模態(tài)函數(shù)（IMF）分量和殘差序列；對每個IMF 分量進行ADF 平穩(wěn)性檢驗，經(jīng)平穩(wěn)檢驗后的各個IMF 分量進行AR（p）模型隨機模擬，然后對殘差序列進行多項式擬合；將隨機模擬后的各IMF 分量序列與多項式擬合后的殘差序列進行重構(gòu)，并在此基礎(chǔ)上還原實測序列中的確定性成分，便得到模擬的年徑流序列；②月徑流隨機模擬，在年徑流實測系列中找出與每個年徑流模擬量最為接近的那個值，以此為典型年，按照該年各月徑流量占全年徑流量的分配系數(shù)，對模擬的年徑流量進行縮放，得到模擬的月徑流序列。

其中，CEEMD方法的具體步驟如下［11］：

（1）設(shè)原始信號為x(t)，為克服模態(tài)混疊，向其添加白噪聲，使得信號分布到合適的振蕩尺度上。白噪聲記為ωi(t)，將原始信號與噪聲疊加，得：

式中：λ0表示噪聲系數(shù)。

（2）將信號x′(t)的局部極大值點用一條曲線連接起來，定義為上包絡(luò)fmax（t）；同樣將所有的局部極小值點也用一條曲線連接定義為下包絡(luò)fmin（t）；并設(shè)上、下包絡(luò)的平均值為m（t），則有：

（3）將信號x′(t)與包絡(luò)的平均值m（t）作差運算：

（4）對信號x′(t)分解N次，即重復(fù)N次步驟（2）、（3），得到第一個IMF分量：

（5）得到第一個IMF分量IFM1（t）后，剩余信號用r1（t）表示：

（6）繼續(xù)執(zhí)行以上步驟，將信號r1(t)+λ1E1[ωj(t)]進行N次分解，第二次分解后的結(jié)果可以表示為：

（7）對分解得到的模態(tài)分量用Ei表示，則第j個剩余的分量可以表示為rj（t）。

（8）對某次分解后的信號r1(t)+λ1Ej[ωj(t)]，對其再次進行分解，可以得到j(luò)+1個分量，表示形式如下：

（9）重復(fù)執(zhí)行以上過程，直至某次IMF分量不可再分時，停止分解過程?？梢缘玫絁個分量，將最終的殘差值記為：

（10）整理以上公式可得原始信號x（t）表達如下：

2 實例應(yīng)用

2.1 研究區(qū)概況

河南省位于我國中東部、黃河中下游，界于東經(jīng)111°21′～116°39′和北緯31°23′～36°22′之間，全省總面積16.7 萬km2，占全國總面積的1.73%。黃河干流自靈寶市進入河南境內(nèi)，流經(jīng)三門峽、濟源、洛陽、鄭州等8 個省轄市和27 個縣（市、區(qū)），河道總長711 km，流域面積3.62 萬km2，占黃河流域總面積的5.1%。為深入研究河南省黃河流域的水文情勢，驗證CEEMD-TD 模型的適用性，本文分別選擇黃河流域的黑石關(guān)站、花園口站、白馬寺站等具有較長徑流實測資料且集水面積不一的典型水文站作為研究對象，各典型水文站的地理位置如圖1所示。

圖1 典型水文站示意圖Fig.1 Schematic diagram of a typical hydrological station

2.2 數(shù)據(jù)來源

收集各典型水文站的徑流實測資料作為基礎(chǔ)數(shù)據(jù)資料，以期為徑流序列的模擬與分析提供理論支持，詳細(xì)信息見表1。資料來源于《河南省水情手冊》，數(shù)據(jù)年限較長，資料經(jīng)過審查，且數(shù)據(jù)全面、合理。

表1 典型水文站信息表Tab.1 Typical hydrographic station information table

2.3 年徑流序列平穩(wěn)性檢驗

CEEMD-TD 模型隨機模擬的基本思路是先對年徑流隨機序列（不含非周期成分與周期成分）進行隨機模擬，再將隨機模擬得到的大量年徑流序列分解為月徑流序列。因此，開展CEEMD-TD 模型隨機模擬前，需剔除年徑流實測序列中非周期成分與周期成分，得到年徑流隨機序列，對剔除確定性成分后的序列隨機項進行CEEMD分解和ADF平穩(wěn)性檢驗。

2.3.1 非周期性成分識別

（1）趨勢診斷項。采用Mann-Kendall 檢驗法［12］、Spearman秩相關(guān)檢驗法［13］分別對黑石關(guān)站、花園口站和白馬寺站的年徑流實測序列進行趨勢項檢驗，結(jié)果見表2。

表2 各站趨勢檢驗結(jié)果Tab.2 Trend test results of each station

結(jié)果表明：Mann-Kendall 檢驗法下黑石關(guān)站等3 個水文站年徑流實測序列統(tǒng)計變量Z的絕對值均大于臨界值Z1—α/2，且變化趨勢β均小于0，說明各站年徑流實測序列均具有顯著的減小趨勢；此外，Spearman秩相關(guān)檢驗法計算得各站統(tǒng)計量T的絕對值均大于自由度為n-2 的t分布在α=0.05 顯著水平下的臨界值t1-α/2，n-2，說明該趨勢顯著，結(jié)果表明Spearman 秩相關(guān)檢驗法的檢驗結(jié)果與Mann-Kendall檢驗法一致。

剔除趨勢項成分后，黑石關(guān)站等3 個水文站年徑流分解序列和實測序列的線性趨勢對比見圖2。

圖2 各站年徑流序列分解前后過程圖Fig.2 Process diagram before and after decomposition of annual runoff sequence of each station

根據(jù)圖2 趨勢項診斷結(jié)果，可看出黑石關(guān)站等3 個水文站分解后的年徑流過程線更加平緩，線性趨勢更加平穩(wěn)，趨近水平線，并在0.05顯著性水平下對3個水文站分解后的各年徑流序列進行一致性診斷，結(jié)果表明分解后的各序列不存在顯著性趨勢。

（2）跳躍診斷項。采用Mann-Kendall 突變檢驗法對黑石關(guān)站等3 個站分解后的年徑流序列進行跳躍項檢驗，檢驗結(jié)果見圖3。選取統(tǒng)計曲線UB和UF在臨界區(qū)間內(nèi)的交點作為最可能變異點。

由圖3 可知，黑石關(guān)站、花園口站、白馬寺站年徑流序列的最可能變異點分別為2002年、2009年及2007年，引起變異的主要原因是在20世紀(jì)90年代后分別發(fā)生了3次大水災(zāi)和旱災(zāi)，以及大規(guī)模水利水電工程、農(nóng)田設(shè)施建設(shè)、城市化進程和水土保

圖3 各站年徑流序列變異檢驗結(jié)果Fig.3 Test results of annual runoff sequence variation at each station

持工程等人類活動的影響，改變了下墊面水文過程的形成條件以及水循環(huán)的途徑與速度。根據(jù)跳躍項診斷結(jié)果，對存在跳躍項的序列進行預(yù)處理，即將變異點之后的序列加上變異前后兩序列均值之差，從而剔除確定性成分中的跳躍項部分。

將各站年徑流序列從變異點處開始重新繪制過程線和趨勢線，并與還原變異前的序列進行對比，結(jié)果如圖4所示。

圖4 各站年徑流序列變異前后過程圖Fig.4 Process diagram before and after variation of annual runoff sequence of each station

對黑石關(guān)站等3個水文站剔除趨勢項和跳躍項后的年徑流序列進行跳躍項診斷，結(jié)果表明經(jīng)上述預(yù)處理后的各序列不存在顯著性變異。

2.3.2 周期性成分識別

對經(jīng)過預(yù)處理后的各水文站年徑流序列xt采用方差線譜法進行周期分析，根據(jù)Fourier 級數(shù)的理論，分別計算出各站ωj對應(yīng)的。

取顯著性水平α=0.05，對各站查F分布表分別得到相應(yīng)的F0.05（2，n-2-1），以黑石關(guān)站為例，黑石關(guān)站僅F7（4.23）大于F0.05（2，60），即其第7 個諧波顯著，則黑石關(guān)站年徑流序列有一個9年的顯著周期。各站年徑流序列周期項的相應(yīng)統(tǒng)計量見表3，其中，aj、bj、Aj為第j個諧波的Fourier 系數(shù)，ωj為對應(yīng)諧波的角頻率，Tj為頻率ωj對應(yīng)的周期。

表3 各站年徑流序列周期項統(tǒng)計量Tab.3 Periodic statistics of annual runoff sequence of each station

分別將各水文站相關(guān)統(tǒng)計量代入周期項公式，識別出各水文站年徑流序列的周期項。以黑石關(guān)站為例，該站年徑流序列的周期項為：

在剔除非周期成分的年徑流序列基礎(chǔ)上，剔除周期項，最終可得序列隨機成分。根據(jù)各站得到的年徑流隨機序列重新繪制過程線和趨勢線，并與還原周期成分前的年徑流序列進行對比，結(jié)果如圖5所示。

圖5 各站年徑流序列還原周期前后過程圖Fig.5 Process diagram of the annual runoff sequence of each station before and after the restoration cycle

對各站剔除確定項成分后的年徑流序列進行周期成分識別，結(jié)果表明各站經(jīng)上述處理后的序列均不存在周期成分。

2.3.3 單位根平穩(wěn)性檢驗

對剔除確定性成分后的各站年徑流隨機序列進行ADF 單位根法檢驗其平穩(wěn)性，結(jié)果發(fā)現(xiàn)：黑石關(guān)站年徑流隨機序列為非平穩(wěn)序列，不能直接采用AR 模型進行隨機模擬；花園口站、白馬寺站的年徑流隨機序列均為平穩(wěn)序列，均可采用AR 模型進行隨機模擬。

為解決黑石關(guān)站年徑流隨機模擬問題，同時也為了提出一種新的、精度高于AR 模型的年徑流隨機模型，本文根據(jù)確定性成分識別成果，將剔除趨勢項、跳躍項和周期項后的黑石關(guān)、花園口和白馬寺3個水文站年徑流隨機序列進行CEEMD 分解，總體平均次數(shù)N取為200 次，噪聲標(biāo)準(zhǔn)差為原始序列標(biāo)準(zhǔn)差的20%。經(jīng)CEEMD 分解后分別得到3 個水文站的各階固有模態(tài)函數(shù)IMF序列及殘差序列，如圖6所示。

圖6 各站年徑流隨機序列CEEMD分解結(jié)果Fig.6 CEEMD decomposition results of annual runoff random sequence at each station

經(jīng)CEEMD 分解后黑石關(guān)站等3 個站的IMF 序列的|t|值均大于顯著性水平α=0.05 時對應(yīng)的臨界值|t0.05|，且p值均小于0.05，即各站IMF 序列均為平穩(wěn)序列，可采用AR 模型進行隨機模擬。

2.4 年徑流隨機模擬

（1）殘差序列的多項式模擬。采用多項式隨機模擬法對黑石關(guān)站等3 個水文站經(jīng)CEEMD 分解后的趨勢項序列分別進行模擬，模擬結(jié)果如下：

式中：Rti為各站趨勢項模擬值（i=1 表示黑石關(guān)站；i=2 表示花園口站；i=3表示白馬寺站）；t為年份序號。

（2）各站年徑流序列AR 模型建立。建立IMF1～IMF5 序列對應(yīng)的AR 模型，采用RIC 準(zhǔn)則進行AR 模型模式識別，確定模型階數(shù)。其中AR（p）模型的數(shù)學(xué)表達式見下式，各站年徑流序列AR模型參數(shù)見表4。

表4 IMF1～IMF5序列AR模型參數(shù)Tab.4 AR model parameters of IMF1～IMF5 sequence

式中：μ為序列均值；σ為序列標(biāo)準(zhǔn)差；φ1，φ2，…，φp為自回歸系數(shù)；εt為純隨機序列。

（3）年徑流序列生成。將黑石關(guān)站等3 個水文站的各IMF序列的AR模型分別進行隨機模擬，均模擬出100組與實測年徑流序列樣本容量相同的序列；然后將各IMF 序列模擬結(jié)果與趨勢項模擬結(jié)果求和，并在此基礎(chǔ)上還原實測序列中的確定性成分，即可分別得到基于CEEMD-AR 模型的各水文站100組年徑流模擬序列。

為檢驗隨機模擬生成的年徑流序列是否具有實測序列相近的統(tǒng)計特征，本文選取年徑流主要統(tǒng)計參數(shù)：均值、均方差σ、離勢系數(shù)Cv、偏態(tài)系數(shù)Cs、一階自相關(guān)系數(shù)r1和二階自相關(guān)系數(shù)r2為檢驗指標(biāo)，通過隨機生成序列與實測徑流序列統(tǒng)計參數(shù)的相對誤差評價模型精度。分別計算黑石關(guān)、花園口和白馬寺站年徑流模擬序列各統(tǒng)計參數(shù)的平均值，并與各站年徑流實測序列相應(yīng)統(tǒng)計參數(shù)進行對比分析，其相對誤差的絕對值A(chǔ)RE（以%表示）如表5所示。由表5可知，黑石關(guān)站等3個站年徑流模擬序列的各項統(tǒng)計參數(shù)與年徑流實測序列相應(yīng)統(tǒng)計參數(shù)的相對誤差均小于10%，說明各站的模擬序列均能較好地反映年徑流實測序列的統(tǒng)計特性。

表5 各站年徑流實測、模擬序列統(tǒng)計特性比較Tab.5 Comparison of statistical characteristics of annual runoff measured and simulated series at each station

為進一步分析CEEMD-AR 模型年徑流模擬的精度，本文對剔除確定性成分后的花園口站和白馬寺站的年徑流隨機序列建立AR 隨機模型，并將隨機模擬結(jié)果與實測序列中的確定性成分疊加，得到基于AR模型的花園口站和白馬寺站100組年徑流模擬序列。從兩水文站年徑流模擬序列統(tǒng)計參數(shù)相對誤差的絕對值A(chǔ)RE（見表5），可以發(fā)現(xiàn)：花園口站CEEMD-AR 模型年徑流模擬序列統(tǒng)計參數(shù)的相對誤差均小于AR 模型，且其最大的相對誤差僅為4.4%，白馬寺站中除CEEMD-AR 模型的均值的相對誤差略大于AR 模型之外，其余統(tǒng)計參數(shù)的相對誤差均小于AR 模型，說明CEEMD-AR 模型相較于AR 模型可以更好地反映年徑流實測序列的統(tǒng)計特性。

2.5 月徑流模擬序列生成

對基于AR 模型與CEEMD-AR 模型的100年徑流模擬序列，在各站年徑流實測系列中分別找出與每個年徑流模擬量最為接近的那個值，以這一年為典型年。按照該年各月徑流量占全年徑流量的分配系數(shù)，將與之對應(yīng)的年徑流模擬量進行解集，最終黃河流域的三個典型水站相應(yīng)的隨機模型均模擬出100組與各站月徑流實測序列樣本容量相同的月徑流隨機模擬序列。

為驗證CEEMD-TD 模型的可靠性，選擇月徑流序列各截口主要統(tǒng)計參數(shù)：均值、均方差σ、離勢系數(shù)Cv、偏態(tài)系數(shù)Cs、一階自相關(guān)系數(shù)r1和二階自相關(guān)系數(shù)r2為檢驗指標(biāo)，模型模擬序列與實測序列各截口統(tǒng)計參數(shù)如圖7、8 和圖9所示。同時，為了便于分析，對于平穩(wěn)序列的花園口站和白馬寺站也采用TD模型進行了月徑流模擬，其模擬序列的統(tǒng)計參數(shù)如圖8 和圖9所示。由圖7、8 和圖9 可知，3 個典型水文站中CEEMD-TD 模型模擬序列和實測序列各截口統(tǒng)計參數(shù)的變化趨勢非常相近，由此可知，該模擬序列總體上均能較好地保持實測序列的主要統(tǒng)計特性。

圖7 黑石關(guān)站月徑流實測及各模擬序列統(tǒng)計特性對比圖Fig.7 Comparison diagram of monthly runoff measured at Heishiguan Station and statistical characteristics of each simulated series

圖8 花園口站月徑流實測及各模擬序列統(tǒng)計特性對比圖Fig.8 Comparison diagram of monthly runoff measured at huayuankou Station and statistical characteristics of each simulated series

統(tǒng)計各截口統(tǒng)計參數(shù)相對誤差發(fā)現(xiàn)：3 個典型水文站中CEEMD-TD 模型模擬序列和實測序列的截口統(tǒng)計參數(shù)相對誤差的絕對值均較小，其中黑石關(guān)站的均值、均方差σ、離勢系數(shù)Cv和偏態(tài)系數(shù)Cs與月徑流實測序列的截口統(tǒng)計參數(shù)的ARE均不超過10%，且均值的ARE 僅3.09%；花園口站和白馬寺的均值的ARE分別為3.05%和3.9%。

與TD 模型的ARE對比分析發(fā)現(xiàn)，除了花園口站的一階自相關(guān)系數(shù)r1高于TD模型1.89%和白馬寺站的一階自相關(guān)系數(shù)r1高于TD 模型4.24%外，其余CEEMD-TD 模型模擬序列的ARE均低于TD模型。因此CEEMD-TD 模型模擬序列總體上均能保持實測序列的主要統(tǒng)計特性，且優(yōu)于TD模型。

3 結(jié)論

（1）本文結(jié)合典型解集（TD）模型，以河南省黃河流域3 個典型水文站為例，將對非平穩(wěn)序列有較強的處理能力的互補式集合經(jīng)驗?zāi)B(tài)分解（CEEMD）法引入到徑流隨機模擬中，提出一種新的互補式集合經(jīng)驗?zāi)B(tài)分解-典型解集（CEEMD-TD）模型。通過模擬表明CEEMD-TD 模型能較好地模擬黃河流域3個典型水文站的月徑流序列，解決典型解集（TD）模型僅適用于平穩(wěn)序列的局限性。

（2）通過對各月徑流序列截口統(tǒng)計參數(shù)來對模擬序列進行評價，結(jié)果表明：CEEMD-TD 模型模擬序列和實測序列的相對誤差的絕對值均較小，且通過與TD 模型模型的模擬序列對比分析，表明在各典型水文站的模擬序列總體上均能較好地保持實測序列的主要統(tǒng)計特性。

（3）由于CEEMD-TD 模型是先模擬總量，再分解成分量，從而導(dǎo)致前一年最后分量與后一年第一分量的自相關(guān)結(jié)構(gòu)不一致的問題，因此對于CEEMD-TD 模型在月徑流的模擬有待后續(xù)進一步研究。□