金 錚，張行南，夏達(dá)忠，舒衛(wèi)民，黃鈺凱

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

0 引 言

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

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

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

1 研究數(shù)據(jù)及方法

1.1 研究數(shù)據(jù)

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

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

1.2 五點(diǎn)三次平滑算法

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

為使方差和最小，令：

可得方程組：

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

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

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

1.3 離散小波閾值消噪算法

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

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

離散小波閾值消噪法的基本思想為通過計算閾值來劃分小波分解后的各層系數(shù)并分別對其進(jìn)行處理，最后進(jìn)行小波逆變換重構(gòu)信號，通?？蓪⒂嬎氵^程分為如下3部分。

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

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

（1）計算細(xì)節(jié)系數(shù)的噪聲強(qiáng)度σ：

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

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

（3）對元素計算風(fēng)險值Rk：

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

（5）對比小波系數(shù)wk與閾值T數(shù)值，進(jìn)行軟閾值消噪：

式中：Wk為經(jīng)過閾值處理的小波系數(shù)；T為閾值。

離散小波逆變換重構(gòu)信號。將離散小波變化分解得到的近似函數(shù)和經(jīng)過閾值量化處理的細(xì)節(jié)函數(shù)進(jìn)行重構(gòu)，得到小波消噪后的信號。小波逆變換如公式（16）所示。

1.4 平滑效果評價指標(biāo)

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

（1）相對誤差。

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

（2）平滑度。

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

2 平滑結(jié)果對比及分析

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

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

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

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

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

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

2.2 不同方法的平滑結(jié)果對比及分析

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

（1）通過選取綜合性最好的Db小波與Sym小波并基于三峽水庫入庫流量特點(diǎn)選取Stein 無偏風(fēng)險估計法對三峽水庫入庫流量進(jìn)行了離散小波閾值消噪，對比各小波平滑評價指標(biāo)結(jié)果，表明Db7小波閾值消噪的平滑效果優(yōu)于其他小波。

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

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