CN101592488A - 水文序列分析中基于信息熵理论的消噪方法 - Google Patents

Abstract

本发明公开了一种水文序列分析中基于信息熵理论的消噪方法，首先依据待分析的水文时间序列的基本特性选择小波函数和小波分解层数，然后对水文序列进行离散小波变换，得到不同时间尺度水平上的小波系数d_j，k；应用小波系数阈值优选熵准则，确定各层的小波系数阈值；对各层高频小波系数进行硬或软阈值量化处理，之后再对处理后的小波系数进行重构，得到实测水文序列中的重构主序列；实测水文序列与重构主序列之差为噪声成分，即实现噪声成分的分离。本发明基于信息熵理论和小波消噪思路，建立了小波系数阈值优选熵准则，从而利用小波分析方法可有效地分离出水文序列中的噪声成分，提高了水文数据的真实性和准确性。

Description

水文序列分析中基于信息熵理论的消噪方法

技术领域

本发明涉及一种水文时间序列分析方法，具体是一种水文序列分析中基于信息熵理论的消噪方法。

背景技术

不确定性现象是水文水资源系统中不可否认和忽视的客观存在，因此长期以来一直是随机水文学中十分重要的研究内容。信息熵理论(如POME(principle of maximumentropy)等)是研究和解决不确定性问题的有效方法。自1957年E.T.Jaynes(Jaynes ET.Information theory and statistical mechanics[J].Phys.Rev.，1957，106：620-630；108：171-190.)首次明确提出POME以来，国内外的许多学者致力于信息熵理论在水文水资源学中的应用研究，并取得令人瞩目的成果。

水文时间序列是研究和解决随机水文学问题的重要基础和依据。随机水文学理论认为，实测水文数据中含有噪声成分，使真实水文序列受到“污染”，造成其真实性和可靠性降低，并对序列成分识别和变化特性分析、水文参数识别和模拟预测等工作产生干扰和影响。因此，深入研究合理有效的水文序列消噪方法具有重要的意义。实际中，若真实水文序列或水文物理演变机制可获知，即可方便地识别并分离噪声成分。然而，目前主要的时间序列模型(随机水文学模型、混沌动力学模型等)都是对真实水文演变机制的一种近似描述，无法准确代表客观现实。因此将实测水文序列和模型模拟序列之间的差作为噪声的处理方式并不合理。

另一类重要的消噪方法主要基于波谱分析方法。传统方法中的维纳滤波、卡尔曼滤波等，只适用于线型系统，且严格依赖于状态空间函数的建立；Fourier变换仅适用于平稳时间序列。具有复杂特性的水文序列显然无法满足这些消噪方法所要求的适用条件。小波分析方法(wavelet analysis，WA)是揭示和研究时间序列多时间尺度变化特性的有效工具，但在应用WA进行水文序列消噪的过程中，也存在以下几个主要问题：(1)水文噪声成分服从的分布类型。与信号小波分析中的噪声不同，将水文噪声成分看作水文独立随机变量时，一般服从偏态分布；(2)小波系数阈值的选择。目前阈值选择方法主要有：固定阈值法(FT)、Stein无偏风险阈值法(SURE)、极大极小原理阈值法(MAXMIN)等，但实际中各种方法互有优缺点，且阈值选择结果也不尽相同；(3)应用这些方法分离出的噪声成分常具有较好的自相关性，这不符合噪声自身的特性。

水文序列中的噪声成分是随机水文过程中许多不确定性因素综合作用的体现和反映。在实际水文序列消噪过程中，需要针对水文问题的特殊性探讨相应适用的消噪方法。

发明内容

本发明所要解决的技术问题是应用信息熵理论及相关原理，并结合小波消噪思路，建立一种更为合理有效的水文时间序列消噪方法。

本发明所述的水文序列分析中基于信息熵理论的消噪方法，包括以下步骤：

1)依据待分析的水文时间序列的基本特性选择小波函数和小波分解层数，然后对水文序列进行离散小波变换，得到不同时间尺度水平上的小波系数d_j，k；

2)应用小波系数阈值优选熵准则，确定各层的小波系数阈值；

3)对各层高频小波系数进行硬或软阈值量化处理，之后再对处理后的小波系数进行重构，得到实测水文序列中的重构主序列；

4)实测水文序列与重构主序列之差为噪声成分，即实现噪声成分的分离。

上述步骤2)应用小波系数阈值优选熵准则，确定各层的小波系数阈值的过程是：

5)首先给定一个较小的第一层小波系数阈值，其余各层的阈值依次按2^-1/2倍的速率递减而确定，然后对各层高频小波系数d_j，k进行硬或软阈值量化处理；

6)利用阈值量化处理之后的小波系数d’_j，k，求解对应重构主序列在各尺度水平j下的能量E_j：

$E_j=\underset{k}{\mathrm{\Σ}}{d^{\′}}_{j,k}^2,j=\mathrm{1,2},...,M---\left(1\right)$

M为小波变换的最大尺度水平，

重构主序列在各尺度水平j下的能量概率分布为P_j：

$P_j=E_j/\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{E}_j,j=\mathrm{1,2},...,M---\left(2\right)$

根据信息熵的定义，将重构主序列各尺度水平下能量概率分布的熵值称为信息量系数ICF，

$\mathrm{ICF}=-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{P}_j\lg P_j---\left(3\right)$

7)对阈值量化处理之后的小波系数进行重构得到主序列，其与原序列之差作为噪声成分，再根据水文物理成因机制选择水文分布线型计算噪声成分的不确定度H值；

8)依次不断增大第一层小波系数阈值的取值，其增大幅度最好取1％，最大值取小波系数中的绝对值最大值。然后确定其余各层阈值，并应用其分别对各层下的高频小波系数进行硬或软阈值量化处理和噪声成分分离，从而得到对应的一系列噪声成分不确定度H值和主序列复杂度ICF值；

9)当小波系数阈值逐渐增大时，噪声成分的H值逐渐增大并趋于稳定，且随着噪声成分的分离，主序列的复杂度会逐渐降低，ICF值也会逐渐减小并趋于稳定；

10)将最终噪声成分H值和主序列ICF值稳定时各层的小波系数阈值作为最优的小波系数阈值。

本发明的有益效果是：水文资料的真实性和可靠性是准确研究水文问题的前提和基础，本发明基于信息熵理论和小波消噪思路，建立了小波系数阈值优选熵准则，从而利用小波分析方法可有效地分离出水文序列中的噪声成分，提高了水文数据的真实性和准确性。算例分析结果显示了该方法有较好的适用性。

附图说明

图1含正态随机成分模拟序列消噪结果，

图2含正态随机成分模拟序列对应的H和ICF变化曲线，

图3含偏态随机成分模拟序列消噪结果，

图4含偏态随机成分模拟序列对应的H和ICF变化曲线，

图5Gauss分布时大山口站月径流序列对应的H和ICF变化曲线，

图6P-III分布时大山口站月径流序列对应的H和ICF变化曲线，

图7Gauss分布时大山口站月径流序列噪声成分分离结果，

图8P-III型分布时大山口站月径流序列噪声成分分离结果，

图9Gauss分布时利津站年径流序列对应的H和ICF变化曲线，

图10P-III型分布时利津站年径流序列对应的H和ICF变化曲线，

图11利津站年径流序列噪声成分分离结果。

具体实施方式

1.1熵函数求解方法

求解概率密度函数f(x)熵函数值的方法如下：(1)根据信息熵的定义(式1)，推导f(x)的熵函数表达式H(x)；(2)应用POME，根据实测序列估计f(x)中的相关参数值；(3)将参数估计值的求解结果代入熵函数表达式H(x)，可确定熵H的具体数值。

H(x)＝∫f(x)lnf(x)dx    (1)

将水文序列中噪声成分看作水文独立随机变量时，一般以P-III分布或Gauss分布描述。两种概率密度函数的熵函数表达式如表1所列(Singh V P.Entropy-based parameterestimation in Hydrology[M]，Kluwer Academic Publishers(Boston/London)，1998.)。

表1 Gauss分布和P-III型分布熵函数表达式

1.2信息量系数ICF求解方法

对水文序列进行小波变换，可得到不同时间尺度上的小波系数，它们是水文序列在不同时间尺度和时间位置上的投影，可用来刻画和描述水文序列的时频结构和多时间尺度变化特性。

应用小波变换系数d_j，k，可求解水文序列在各尺度水平j下的能量E_j：

$E_j=\underset{k}{\mathrm{\Σ}}{d}_{j,k}^2,j=\mathrm{1,2},...,M---\left(2\right)$

M为小波变换的最大尺度水平。水文序列在各尺度水平j下的能量概率分布为P_j：

$P_j=E_j/\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{E}_j,j=\mathrm{1,2},...,M---\left(3\right)$

根据信息熵的定义，将水文序列各尺度下能量概率分布的熵值称为信息量系数(information cost function，ICF)。ICF反映了系统的能量分布，表征了系统的复杂程度。当水文时间序列越复杂越无序，ICF值越大；反之，当序列的能量集中于某一频带，时间域相依性强时，复杂性就越弱，其ICF值也相对较小。

$\mathrm{ICF}=-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{P}_j\lg P_j---\left(4\right)$

1.3水文时间序列消噪新方法

应用小波分析方法可以实现信号中不同成分的分离。水文序列中不同成分在小波变换之后表现出不同的特性：(1)确定性成分(趋势、周期等)的能量在时域上分布不均匀，其主要集中在少数的几个频率值上(对应着周期值)，且小波变换模极大值随时间尺度水平的增大而增大；(2)噪声成分对应的能量均匀分散，随着时间尺度增大小波变换模极大值减小，且大致按2^-1/2倍的速率衰减，一般情况下经两次分解后的高频成分可包含信号中所有的噪声成分。基于上述两点，确定合理的阈值对小波系数进行处理之后再重构，可实现噪声成分的分离，此即小波消噪的基本思路。

依据上述小波消噪思路进行水文时间序列中噪声成分分离时，需要解决如下两个关键问题：小波系数阈值的合理确定和消噪结果优劣的判别，下面分述之。

1.3.1小波系数阈值优选熵准则

为便于叙述，本发明将低频确定性成分构成的序列称为主序列，即水文序列主要分为主序列和噪声两部分。一方面，将水文序列中的噪声成分看作水文随机变量，用某种分布线型(一般以P-III型分布或Gauss分布描述，选用时应具体分析确定)加以描述，并用熵函数值H描述噪声成分的不确定程度；另一方面，用信息量系数ICF描述主序列的复杂程度。

对于待分析水文序列，其含有的噪声成分的不确定度和主序列的复杂度是确定的。因此，根据上述主序列和噪声成分对应小波系数的不同特性，可按下述步骤优选合理的小波系数阈值：(1)选择合理的小波函数(桑燕芳，王栋.水文序列小波分析小波函数选择方法[J].水利学报，2008，39(3)：296-300，306.)和小波分解层数，对水文序列进行小波分解；(2)给定一个较小的第一层小波系数阈值，其余各层的阈值依次按2^-1/2倍的速率递减而确定；(3)依次不断增大第一层小波系数阈值的取值，然后确定其余各层阈值，并应用其分别对各层下的高频小波系数进行阈值量化处理和噪声成分分离，从而得到对应的一系列噪声成分不确定度H值和主序列复杂度ICF值；(4)随着噪声成分逐渐被完全分离，H值逐渐增大并趋于稳定。且随着噪声成分的分离，主序列的复杂度会逐渐降低，ICF值也会逐渐减小并趋于稳定；(5)将最终噪声成分H值和主序列ICF值稳定时各层的小波系数阈值作为最优的小波系数阈值。此即本文提出的水文时间序列消噪时的小波系数阈值优选熵准则。

具体的小波系数阈值量化过程中，可采用效果较优的软阈值处理方法。且实际中分析水文时间序列时，应同时应用P-III型分布和Gauss分布描述分离出的噪声成分，然后综合分析并最终确定合理的小波系数阈值。

1.3.2序列消噪结果优劣的判别依据

若水文序列消噪结果合理，原序列、主序列和噪声成分的特征值(x、σ、C_s、r₁)应符合以下规律：(1)主序列和原序列的x相差不应过大，噪声成分的x为原序列和主序列x之差；(2)相比原序列，消除噪声成分干扰后主序列的σ应有所减小；(3)主序列和原序列的C_s应相近；(4)主序列的r₁应较原序列增大；噪声成分无相关性，所以其r₁应接近于0。本文中，应用该四条规律对不同方法对应序列消噪结果的合理性和优劣性进行分析和判断。

2、算例分析

为验证水文序列消噪时，本文建立的小波系数阈值优选熵准则及消噪方法的合理性和适用性，分别对模拟序列和实测水文序列进行分析。

2.1模拟序列分析

为讨论该阈值优选熵准则对不同分布类型噪声成分的适用性，生成含有正态(Gauss分布)和偏态(P-III型分布)噪声成分的模拟序列，分别进行分析。

2.1.1模拟序列1

应用Monte-Carlo方法生成一500a的模拟序列(图1)，其中含有一线性趋势项和20a、50a两个周期项，另叠加有x为0，σ为10的Gauss随机成分。应用“dmey”小波取水平4进行离散小波分解，然后进行序列消噪。其中，应用Gauss分布描述分离出的噪声成分。

分析结果显示(图2)：第一层小波系数阈值增大到12之后H值和ICF值均趋于稳定。因此，选择12作为第一层的最优小波系数阈值，其余各层阈值依次按2^-1/2倍进行递减确定，然后进行噪声成分分离，结果如图1所示。

表2 不同方法对模拟序列1消噪结果对比分析

图1中，上部曲线为含有Gauss随机成分的模拟序列；中部为不含随机成分的真实序列；下部为消噪之后的主序列。真实序列和主序列对比可以看出，消噪结果良好。

然后统计并分析原模拟序列、主序列和噪声成分的特征值，进一步判别消噪结果的优劣，并与常用小波消噪方法的结果进行对比(表2)。

表2中数据显示：(1)应用熵准则去除噪声成分后，三个序列(原序列、主序列和噪声成分)特征值之间的大小关系符合前述的四条规律，噪声成分的r₁、C_s接近于零，服从Gauss分布，表明了消噪结果的合理性；(2)FT、SURE、MAXMIN三种方法分离出的噪声成分r₁值均较大，即自相关性较好，这不符合独立随机变量的特征，表明噪声成分分离结果中含有确定性成分，结果存在一定误差；(3)对比分析显示，熵准则的噪声成分分离结果好于其它三种方法。

2.1.2模拟序列2

同样应用Monte-Carlo方法生成一500a的模拟序列(图3)，其与模拟序列1中的真实序列相同，但此处叠加均值为0，σ为10，C_s为0.8的P-III型分布的随机成分。应用“dmey”小波函数，取水平4进行离散小波分解。然后进行序列消噪。其中，应用P-III型分布描述分离出的噪声成分。

由图4可以看出，H值和ICF值在24之后分别趋于稳定。因此，取24作为第一层最优小波系数阈值，其余各层阈值依次按2^-1/2倍递减确定，最后进行噪声成分分离(图3)。

图3中，上部曲线为叠加P-III型随机成分的模拟序列。中部曲线为真实序列，下部曲线为消噪之后的主序列。对比分析结果同样显示了消噪结果良好。

分别统计原模拟序列、主序列和噪声成分的特征值，并与其他消噪方法的结果对比(表3)。结果显示：(1)熵准则的噪声成分分离结果更为合理，优于其它三种方法；(2)其它三种方法分离出的噪声成分具有一定的自相关性。SURE的结果相对优于FT和MAXMIN方法；(3)四种方法分离出的噪声成分都表现出偏态特性，但熵准则对应噪声成分的C_s值更接近于设计值0.8。

表3 不同方法对模拟序列2消噪结果对比分析

两模拟序列分析结果表明：(1)时间序列噪声成分分离过程中，本文建立的小波系数阈值优选熵准则及消噪方法对含有正态、偏态随机成分序列的分析结果均良好，因此该熵准则有较好的适用范围；(2)该消噪方法的分析结果优于常用的小波消噪方法。

2.2实测水文序列分析

2.2.1实例1

选用大山口水文站实测的20年(1978-1997年)月径流序列为例(数据取自左其亭，高峰.水文时间序列周期叠加预测模型及3种改进模型[J].郑州大学学报，2004，25(4)：67-73.)，选取“demy”小波函数，取水平3进行离散小波变换。分别用Gauss分布和P-III型分布描述分离出的噪声成分。由于该实测月径流序列受人类活动影响较小，噪声成分对真实序列的干扰程度也较小，因此实际消噪过程中仅对第一层高频小波系数进行阈值量化处理。

应用Gauss分布和P-III型分布描述噪声成分时的第一层小波系数阈值优选结果分别为78和83(图5和图6)。然后对该序列进行消噪。图7和图8中，上部曲线均为实测的月径流序列，中部曲线为消噪之后的重构主序列，下部为分离出的噪声成分。两阈值对应的消噪结果相近，且与原序列相比，重构主序列均很好地保持了原序列的变化规律，而噪声成分表现出不规则的高频震荡。对原实测序列、重构主序列和噪声成分的特征值进行分析，同样验证了消噪结果的合理性(表4)。

表4 不同方法对大山口站月径流序列消噪结果对比分析

2.2.2实例2

选用黄河利津站实测的54年年径流序列为例，选取“bior3.5”小波函数，取水平3进行离散小波变换。分别用Gauss分布和P-III型分布描述分离出的噪声成分。

应用Gauss分布和P-III型分布描述噪声成分时的第一层小波系数阈值优选结果均为9.2(图9和图10)，然后对第一、二层的高频小波系数进行阈值量化处理和噪声成分分离。图11中，上部曲线为利津站实测年径流序列，中部曲线为主序列，下部为分离出的噪声成分。噪声成分的分离结果和各序列特征值分析结果(表5)同时验证了所选阈值的合理性和有效性。

表5 不同方法对利津站年径流序列消噪结果对比分析

本发明应用信息熵理论并结合小波消噪的基本思路，建立了水文时间序列消噪过程中小波系数阈值优选熵准则和序列消噪新方法。分别对不同特性的模拟序列和不同的实测水文序列进行分析，验证了该方法的有效性和适用性。综合分析可得到以下几点认识：

(1)该阈值优选熵准则基于信息熵理论和小波消噪思路而建立，其通过分析噪声成分H值和主序列ICF值的变化规律，进而确定合理的小波系数阈值。根据原序列、主序列和噪声成分特征值之间的变化规律，并与常用小波消噪方法的分析结果进行对比，验证了该方法的合理性。对不同模拟序列(含正态和偏态随机成分)和不同实测序列的分析结果显示了该方法的适用性。

(2)实测序列的分析结果显示，水文序列中的噪声成分具有偏态特性。因此本文中应用偏态分布线型(P-III型分布)对分离出的噪声成分进行描述相对更为合理。

(3)由于利津站54年年径流序列中噪声成分的偏态特性较弱，因此应用Gauss分布和P-III型分布描述噪声成分时优选的小波系数阈值结果相同。

(4)对本发明中各序列进行分析时，常用小波消噪方法分离出的噪声成分都有较好的自相关性，这不符合噪声自身的特性，因此结果具有一定的不合理性。

(5)由于本文建立的小波系数阈值优选熵准则基于信息熵理论，因此最终优选的小波系数阈值是整体最优的。算例分析结果也显示其优于常用的小波消噪方法。

Claims (3)

1、一种水文序列分析中基于信息熵理论的消噪方法，其特征在于包括以下步骤：

1)依据待分析的水文时间序列的基本特性选择小波函数和小波分解层数，然后对水文序列进行离散小波变换，得到不同时间尺度水平上的小波系数d_j，k；

2)应用小波系数阈值优选熵准则，确定各层的小波系数阈值；

3)对各层高频小波系数进行硬或软阈值量化处理，之后再对处理后的小波系数进行重构，得到实测水文序列中的重构主序列；

4)实测水文序列与重构主序列之差为噪声成分，即实现噪声成分的分离。

2、根据权利要求1所述的水文序列分析中基于信息熵理论的消噪方法，其特征在于上述步骤2)的确定过程是：

5)首先给定一个较小的第一层小波系数阈值，其余各层的阈值依次按2-1/2倍的速率递减而确定，然后对各层高频小波系数d_j，k进行硬或软阈值量化处理；

6)利用阈值量化处理之后的小波系数d’_j，k，求解对应重构主序列在各尺度水平j下的能量E_j：

$E_j=\underset{k}{\mathrm{\Σ}}{d^{\′}}_{j,k}^2$ j＝1，2，...，M    (1)

M为小波变换的最大尺度水平，

重构主序列在各尺度水平j下的能量概率分布为P_j：

$P_j=E_j/\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{E}_j$ j＝1，2，...，M    (2)

根据信息熵的定义，将重构主序列各尺度水平下能量概率分布的熵值称为信息量系数ICF，

$\mathrm{ICF}=-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{P}_j\lg P_j---\left(3\right)$

7)对阈值量化处理之后的小波系数进行重构得到主序列，其与原序列之差作为噪声成分，再根据水文物理成因机制选择水文分布线型计算噪声成分的不确定度H值；

8)依次不断增大第一层小波系数阈值的取值，然后确定其余各层阈值，并应用其分别对各层下的高频小波系数进行硬或软阈值量化处理和噪声成分分离，从而得到对应的一系列噪声成分不确定度H值和主序列复杂度ICF值；

9)当小波系数阈值逐渐增大时，噪声成分的H值逐渐增大并趋于稳定，且随着噪声成分的分离，主序列的复杂度会逐渐降低，ICF值也会逐渐减小并趋于稳定；

10)将最终噪声成分H值和主序列ICF值稳定时各层的小波系数阈值作为最优的小波系数阈值。

3、根据权利要求2所述的水文序列分析中基于信息熵理论的消噪方法，其特征在于上述步骤8)增大幅度取1％，最大值取小波系数中的绝对值最大值。

Images