基于多种平滑算法的三峡水库反推入库流量消噪

2022-03-23张行南夏达忠舒卫民黄钰凯

中国农村水利水电 2022年3期

金铮，张行南，夏达忠，舒卫民，黄钰凯

（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小波阈值消噪法对原入库流量锯齿问题的处理结果更好，该方法减小了原序列波动的振幅，使平滑结果基本贯穿了入库流量的中心轴线，可以认为该方法的平滑结果较理想。