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

王倩麗，馬細霞，2，趙 璐，張靜文，程 旭

（1.鄭州大學水利科學與工程學院，鄭州450001；2.鄭州大學黃河生態保護與區域協調發展研究院，鄭州450001；3.鄭州航空港興港投資集團有限公司，鄭州450001；4.商丘工學院土木工程學院，河南商丘476000）

0 引言

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

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

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

1 CEEMD-TD模型介紹

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

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

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

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

式中：λ0表示噪聲系數。

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

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

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

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

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

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

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

（9）重復執行以上過程，直至某次IMF分量不可再分時，停止分解過程。可以得到J個分量，將最終的殘差值記為：

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

2 實例應用

2.1 研究區概況

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

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

2.2 數據來源

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

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

2.3 年徑流序列平穩性檢驗

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

2.3.1 非周期性成分識別

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3.2 周期性成分識別

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

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

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

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

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

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

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

2.3.3 單位根平穩性檢驗

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

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

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

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

2.4 年徑流隨機模擬

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

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

（2）各站年徑流序列AR 模型建立。建立IMF1～IMF5 序列對應的AR 模型，采用RIC 準則進行AR 模型模式識別，確定模型階數。其中AR（p）模型的數學表達式見下式，各站年徑流序列AR模型參數見表4。

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

式中：μ為序列均值；σ為序列標準差；φ1，φ2，…，φp為自回歸系數；εt為純隨機序列。

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

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

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

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

2.5 月徑流模擬序列生成

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

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

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

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

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

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

3 結論

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

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

（3）由于CEEMD-TD 模型是先模擬總量，再分解成分量，從而導致前一年最后分量與后一年第一分量的自相關結構不一致的問題，因此對于CEEMD-TD 模型在月徑流的模擬有待后續進一步研究。□