基于多種平滑算法的三峽水庫反推入庫流量消噪

金 錚，張行南，夏達忠，舒衛民，黃鈺凱

（1.河海大學水文水資源學院，南京 210098；2.中國長江電力股份有限公司，湖北宜昌 443000）

0 引 言

三峽水庫作為典型的狹長河道型水庫，動庫容巨大，其入庫流量無法通過實測或傳統的靜庫容推算方法獲得，現階段三峽水庫入庫流量的計算采用“八段動庫容法”，基于水量方程每2 個小時由出庫流量和水庫庫容變化反推求得［1］，但是由于水庫面積較大，且水庫泄洪能力和庫容曲線等都存在一定的誤差等原因［2］，反推的入庫流量疊加了各種噪聲，其曲線表現出明顯的“鋸齒”現象，因此需要尋找合理的平滑算法以減少或去除反推流量中的誤差。

數據平滑算法是處理水文時間序列隨機誤差的常用手段，其可分為局部平滑和大范圍平滑算法兩類［3］。局部平滑算法通過連續滑動窗口并賦予擬合多項式不同權重系數來平滑窗口內數據，常見的水文序列平滑算法有移動平均法、指數平滑法及五點三次平滑法，相比于前兩種算法造成的過度修正，后者平滑的效果體現出更高的靈活性，五點三次平滑法是經過改進的Savitzky-Golay濾波器，在最小二乘法的基礎上加入了數據平滑技術，使得該算法更加注重數據的細節處理［4］，該方法是一種常用的平滑預處理方法，可在消除數據干擾成分的同時保證原有曲線特性不變，結構簡單、運算速度快，目前已有較多應用在入庫流量平滑的研究，例如武煒等［5］通過利用廣西多個電站的入庫流量進行驗證及李凱凱［6］在山西張峰水庫中的應用等，研究顯示五點三次平滑法在平滑入庫流量的問題中取得了較好的結果。大范圍平滑算法從數據整體的角度出發，利用有效信號和噪聲二者頻率的差別，通過不同的濾波器裁剪或削弱高頻噪聲信號以保留更多的低頻有效信息，包括傅里葉級平滑法、卡爾曼濾波法及離散小波消噪法等，但傅里葉級平滑僅適用于平穩和線性時間序列，而卡爾曼濾波適用于線性系統，且嚴格依賴于狀態空間的建立。相比之下，小波分析方法具有良好的時-頻特性且更適用于水文時間序列的非穩定特性［7］，近些年離散小波閾值消噪開始被用于水文時間序列的噪聲處理并初具成效，王鶴［8］利用Db2小波對長春市的月降水數據進行消噪，驗證了該算法的可行性和有效性，王秀杰等［9］選用Db4 小波對黃河頭道拐水文站徑流數據進行消噪，較好的保留了徑流序列的尖峰和突變部分。

實際工作中，三峽水庫入庫流量的平滑通常采用操作簡便且計算量小的三點平滑法進行處理，但結果顯示出明顯的滯后性，故本文通過時間序列數據平滑算法的對比選取局部平滑中的五點三次平滑法以及大范圍平滑中的離散小波閾值消噪法進行研究，希望通過選取合適的平滑算法，既達到減小“鋸齒”波動的目的，又能保證原入庫流量序列的水文特征不變，為三峽水庫入庫流量的平滑處理提供依據。

1 研究數據及方法

1.1 研究數據

本文選取三峽水庫2019年1月1日至2019年12月31日以2 h 為時間間隔反推得到的入庫流量序列為研究對象，考慮到豐水期和枯水期河道基流、河道匯流時間有所不同，結合實際水庫調度，將2019 全年入庫流量過程劃分為消落期、主汛期和蓄水期，各時段始末時間如表1所示。

表1 2019年入庫流量時段劃分Tab.1 Division of the inflow in 2019

1.2 五點三次平滑算法

設在2n+1 個等距為h的節點x-n，x-n+1，…，x0，…，xn-1，xn上，數據分別為Y-n，Y-n+1，…，Y0，…，Yn-1，Yn，作交換t=(x-x0)/h，則原節點變為t-n= -n，…，t0= 0，…，tn=n，建立m次多項式（1）擬合數據，系數由最小二乘法確定，方差和如（2）。

為使方差和最小，令：

可得方程組：

當n取2（即5 個節點）、m取3 時，將計算所得的a0，a1，a2帶入公式（1），同時令t= -2，- 1，0，1，2，得到五點三次平滑公式。設三峽水庫的入庫流量序列為Q1，Q2，…，Qt，…，Qn-1，Qn（n＞5），將Qn代入五點三次平滑公式中的Yn，則得到五點三次平滑計算公式如下［5］。當節點個數超過5時，除在序列兩端4個節點分別使用式（5）、式（6）、式（8）及式（9）進行平滑，其余節點均用式（7）進行平滑。

為保證水量平衡，需要對五點三次平滑結果進行二次修正。首先計算各時段累計誤差，由于一定范圍內的誤差是正常的，因此本文根據允許誤差設置累計誤差閾值，當累計誤差絕對值大于閾值時，認為該時段平滑值誤差較大，需要進行修正；當累計誤差絕對值小于等于閾值時，則認為該時段平滑結果誤差合理，不需要進行修正。公式如下。

式中：Et為t時段平滑后序列的累計誤差，m3/s；et為第t時段的誤差值，m3/s；Eh為累計誤差閾值，m3/s；b為允許誤差，%；b的選取由各水庫根據經驗以及對修正流量的平滑程度要求選取；Wa為原始序列總洪量（2×3 600 m3）。

1.3 離散小波閾值消噪算法

離散小波變換在水文時間序列上的研究主要為識別和分離序列的不同組成成分，同時為小波閾值消噪和小波模擬預報等過程提供方便［7］。離散小波變換定義為：

式中：f(t)為原始信號，本文中為入庫流量，m3/s為小波基函數ψt的復共軛函數；Wf(j，k)為離散小波變換系數，可分為近似函數和細節函數；j為分解水平，也稱時間尺度水平；k為時間位置因子；a0和b0均為常數。

離散小波閾值消噪法的基本思想為通過計算閾值來劃分小波分解后的各層系數并分別對其進行處理，最后進行小波逆變換重構信號，通常可將計算過程分為如下3部分。

離散小波變換分解信號。通過選取合適的小波基函數和分解水平，依據公式（11）對原始數據進行離散小波變換，從而得到各分解水平上的離散小波變換系數。水文序列可由確定性成分（趨勢、周期等）和噪聲成分組成，經過小波分解，噪聲成分被包含在各分解水平的高頻成分中，對應著較小的小波系數，而確定性成分則在各分解水平的低頻成分中，對應著較大的小波系數［7，9］。本文分別選擇綜合適用性較好及被較廣泛使用的Symlets（SymN，取N=5，6，7，8）小波及Daubechies（DbN，取N=5，6，…，10）［12］小波對入庫流量進行消噪處理，一般對于波動性強的序列分解水平一般不超過3 層［13］，結合入庫流量的特點本文將分解水平設置為3。

閾值確定及閾值量化處理方法。此計算的目的是利用噪聲成分和確定性成分小波系數的差異將噪聲削弱或去掉，同時保留有用信號。估計小波系數閾值的方法包括無偏風險估計閾值法、固定閾值法和極大極小閾值法等。由于本文的研究序列時間間隔為2 小時，小于日尺度，屬于高頻序列，故選取更適合來處理含有噪聲的高頻序列方法——Stein 無偏風險估計法來計算閾值［12］，過程如下步驟（1）～（4）所示。閾值量化處理方法可理解為利用閾值將小波系數分類，認為絕對值小于閾值的小波系數為噪聲并用零替代它，而超過閾值的小波系數則經過縮減對其重新取值［10］，通常閾值處理方法可分為軟閾值法和硬閾值法。硬閾值法處理結果較粗糙，其函數的不連續性會造成信號的震蕩，相比之下軟閾值處理的結果更平滑，更符合本文的研究目的，其計算過程可見步驟（5）。以某一層細節系數的計算過程如下。

（1）計算細節系數的噪聲強度σ：

式中：wk為分解得到的小波系數中第k個系數；N為序列長度。

（2）將細節系數中各元素取絕對值后升序排列，再取各元素平方值，得到新的系數1 ≤k≤N)。

（3）對元素計算風險值Rk：

（4）求風險值為最小時的，得到閾值T如下式：

（5）對比小波系數wk與閾值T數值，進行軟閾值消噪：

式中：Wk為經過閾值處理的小波系數；T為閾值。

離散小波逆變換重構信號。將離散小波變化分解得到的近似函數和經過閾值量化處理的細節函數進行重構，得到小波消噪后的信號。小波逆變換如公式（16）所示。

1.4 平滑效果評價指標

本文將分別從兩方面對平滑結果進行判別，一是從對水文的洪水過程影響重大的三個要素出發，即洪峰、洪量和峰現時間，通過計算其相對誤差評價平滑效果；二是從曲線平滑程度的角度出發，因標準差可以用于描述平滑曲線的光滑度，且其值越小曲線越光滑［11］，故引入平滑度指標評價入庫流量的平滑程度，該指標由平滑后流量與原始流量的一階差分的標準差比值構成，平滑度范圍在0～1 之間，且越接近1 表明平滑后流量序列相比原始流量的平滑程度更高。

（1）相對誤差。

式中：EBIAS為相對誤差；Hori，i和Hde，i分別為平滑前后的洪峰，m3/s、洪量，m3和峰現時間，h。

（2）平滑度。

式中：Esm為平滑度，其值越大表示曲線越平滑；Qori，i、Qori，i+1分別為第i、i+1時段平滑前的入庫流量，m3/s；Qde，i、Qde，i+1分別為第i、i+1時段平滑后的入庫流量，m3/s。

2 平滑結果對比及分析

2.1 離散小波閾值去噪算法研究

將入庫流量序列為Q1，Q2，…，Qt-1，Qt代入公式（11）中的f(t)，經過離散小波變換分解、Stein 無偏風險估計法計算閾值、硬閾值法量化處理小波系數以及逆變換重構信號，得到取兩種小波中各小波函數時的去噪入庫流量序列，計算結果見表2。

表2 Sym小波和Db小波對全年過程的消噪結果Tab.2 The denoising results of the whole year process by Sym wavelet and Db wavelet

由表中計算結果可得，各小波消噪的峰現時間誤差均為0，可以很好的保證峰現時間的準確性。同時，原始序列的洪峰流量為46 448 m3/s，由洪峰的誤差計算可以看出，相對誤差為-2.074%（Sym7）和-0.010%（Db10）的實際洪水絕對誤差分別為-963 和-5 m3/s，由計算結果可以看出，除Sym7、Db5、Db9 外其他小波的洪峰誤差均在±1.000%區間內，較好的保留了洪峰洪量的量級。在剩余各小波中，除Sym8和Db8小波外的其余小波的洪量誤差均保持在±0.001%內，當洪量的相對誤差在允許范圍內，可以認為該過程維持了水量平衡。從平滑度指標的結果顯示可得，Sym5 小波及Db7 小波的平滑度為分別為0.56 和0.75，表明其平滑效果更好。綜合各指標結果可認為Sym5小波及Db7小波的平滑結果更理想。

為更進一步對比兩種小波的平滑結果，分別評價其在消落期、主汛期及蓄水期時段中的消噪情況，截取局部表現如圖1所示。對比各時期兩個小波表現所得，Db7小波較Sym5平滑程度更大，即消除的噪聲部分更多，且平滑結果的波動更小，基本達到了平滑鋸齒的目的，同時在流量過程中轉折點的平滑上，其在時間點的表現上更準確。故本文選取Db7 小波作為離散小波閾值消噪法的小波函數。

圖1 各劃分時期Sym5小波和Db7小波的局部平滑對比Fig.1 Partial smoothing comparison between Sym5 wavelet and Db7 wavelet

2.2 不同方法的平滑結果對比及分析

本文在三點平滑計算的過程中，將每個窗口平滑得到的流量替換其原始流量，并放入下一個窗口的平滑計算中依次計算。根據水量平衡，應用公式（10）對五點三次法平滑結果進行修正，選取閾值為0.001%，修正后的累積誤差如圖2所示。

圖2 五點三次平滑法經水量平衡修正后的流量累積誤差Fig.2 Flow accumulation error corrected by water balance with five-point cubic smoothing method

三點平滑法、五點三次平滑法及Db7 小波閾值消噪法對全年的平滑結果如表3所示。

表3 3種平滑方法對全年序列的處理結果Tab.3 The processing results of three smoothing methods for the whole year sequence

由表計算結果可得，五點三次平滑的峰現時間出現了一個時段的推遲，通過比較洪峰和洪量的平滑誤差可以看出Db7 小波閾值消噪法更能保證原入庫流量的水文特征值，在平滑度的計算上，該方法也表現出了明顯的優勢，可以認為其消除了更多劇烈的波動。

為進一步比較平滑效果，分別選取全年三次峰值較大的洪水過程，即8月4日-8月15日、7月17日-7月29日及6月21日-6月28日的洪水過程（記為洪水1、2、3）進行平滑，3種方法平滑結果局部如圖3所示，其中圖3（d）、（e）和（f）分別為三次洪水的洪峰平滑過程，圖3（a）為洪水2 的起漲階段平滑過程，圖3（b）為洪水1 和2 之間平滑過程，圖3（c）為洪水3 的落水階段平滑過程。

由圖3 可以看出，五點三次的平滑效果與三點平滑的結果基本重合，但仍存在較大的波動，尤其是在原入庫流量波動劇烈的部分，例如洪水的漲水及落水階段，而Db7 小波消噪對此的平滑結果相對更好，可以看出其基本消除了幅度較大的鋸齒，還原了自然狀態下的洪水過程，雖有波動但是幅度微小，在流量變化較小的時段里認為合理。由于入庫流量存在波動使得計算峰值偏大，故經過平滑其峰值將有所下降，同時考慮三點平滑法有過度平滑的劣勢，而五點三次平滑結果與其基本重合，所以在峰值上的平滑本研究認為Db7小波消噪效果更好。

同時，根據本文的時段劃分對比3種方法在消落期、主汛期及蓄水期的平滑結果分別如表4 及圖4 所示，其中圖4（a）、（b）為消落期局部平滑過程，圖4（c）、（d）為蓄水期局部平滑過程。

由于入庫流量時間序列在3 個時期呈現不同的特征，即消落期和蓄水期的變化頻率較快，振幅較小，而主汛期的變化則為頻率較小，振幅較大，故而3 種方法在3 個時期的平滑效果也略有不同。由表4 計算結果可得，Db7 小波閾值消噪法在消落期及蓄水期的平滑效果更好，既維持了原序列的洪量量級又提高了平滑結果的平滑度，而在主汛期中，五點三次平滑算法的平滑度較小波閾值消噪法稍大。同時由于庫區較大，實際過程中其對來水有一定的調蓄作用，故實際入庫流量不會出現較大的波動，由圖4 中各方法的平滑結果可以明顯的看出，Db7小波閾值消噪法對原入庫流量鋸齒問題的處理結果更好，該方法減小了原序列波動的振幅，使平滑結果基本貫穿了入庫流量的中心軸線，可以認為該方法的平滑結果較理想。

圖4 消落期及蓄水期局部平滑結果Fig.4 Partial smoothing results of fluctuation period and water storage period

3 結 論

本文基于三點平滑算法、五點三次平滑算法及離散小波閾值消噪算法，對2019年三峽水庫的入庫流量中的“鋸齒”問題進行平滑，同時采用相對誤差及平滑度指標對其平滑效果進行評價。

（1）通過選取綜合性最好的Db小波與Sym小波并基于三峽水庫入庫流量特點選取Stein 無偏風險估計法對三峽水庫入庫流量進行了離散小波閾值消噪，對比各小波平滑評價指標結果，表明Db7小波閾值消噪的平滑效果優于其他小波。

（2）基于水庫入庫流量在消落期、蓄水期及主汛期的特點，分別對比三點平滑算法、五點三次平滑算法及Db7 小波閾值消噪算法平滑效果，結果表明Db7 小波閾值消噪算法呈現出更理想的平滑效果，保證了原入庫流量水文特征的同時，消除了更多的鋸齒。

（3）離散小波閾值消噪算法在應用過程中，平滑效果隨小波基函數選取的不同而有差別，如何高效地選取小波基函數對三峽入庫流量進行平滑有待進一步的研究。 □