基于小波多孔算法的黑河径流变异规律分析

摘 要：

关键词:

径流; 变异; 二次样条小波; 多孔算法; 黑河；

朱新丽(1989—),女，博士研究生，主要从事河川径流演变。E-mail:zhuxl890@126.com;*李彦彬(1973—),男，教授，博士，主要从事农业水土工程、水资源高效利用及水旱灾害管理、河川径流演变。E-mail:Liyb101@sina.com;

基金：

国家自然科学基金(51779093);

引用：

朱新丽, 李彦彬, 李红星, 等. 基于小波多孔算法的黑河径流变异规律分析[ J] . 水利水电技术(中英文) , 2021, 52( 10) : 46- 58.

ZHU Xinli, LI Yanbin, LI Hongxing, et al. Analysis of Hei River runoff variation law based on A′tuous algorithm[J]. Water Resources and Hydropower Engineering, 2021, 52(10): 46- 58

---

0 引 言

小波分析，具有时频多分辨功能，被誉为数学“显微镜”,是进行信号处理的理想工具。自小波理论发展以来，在各领域得到广泛应用。由于水文时间序列的非平稳性、复杂性，以及时域和频域随时间变化特性，小波分析为分析水文时间序列变化特性奠定了坚实的基础。小波分析应用于水文学领域只有30a左右时间，且主要在水文时间序列多时间尺度分析、变化特性分析、预测预 等方面。而在其他领域，应用于信号奇异性检测或故障检测等方面也较广泛。但水文学领域研究较少：2001年薛小杰等应用小波分析将水文时间序列分解成不同尺度下低频和高频成分，然后对低频成分进行重构，判断水文时间序列趋势变化，并应用于黄河上中下四个水文站，证明多分辨分析用于水文时间序列趋势分析的合理可行性；2003年王文圣等通过观察小波多尺度分析的年径流量的丰枯转折点，得出水文序列的突变特征。但是都没有形成系统。目前，随着小波理论的完善，小波分析在突变检测方面存在的突变点漂移等问题也得到改善。在多分辨率分析中,多孔算法的平移不变性，以及能够保证采样信号的完整性,奇异点检测的准确性,而逐渐被广泛用于突变点检测。其分解得到的近似系数，可以反映信号的趋势特征，细节系数保留了信号的局部特性和噪声特性，有效地克服了Mallat算法需要大量的采样数据信息和不能直接对分解信号单支重构处理的问题。B样条小波，具有紧支集、全正性和线性相位性，能够满足奇异点检测所需小波的光滑条件,尤其是二次B样条小波，具有良好的局部化特性，是奇异点检测的最佳选择。

黑河作为我国西北地区第二大内陆河，其流域水资源短缺已成为困扰其生存和发展的突出问题。因此，各位学者对黑河流域径流变化过程研究仍旧比较关注。王钧等对黑河莺落峡和正义峡径流量变化周期、突变和渐变趋势进行分析；刘浏等采用多种方法，对黑河流域53年水文气象要素变化特征进行分析，并指出单独1种或2种变异点检测方法不能完全检测出序列存在的突变点；聂晶和吴彦昭研究了黑河近60年来的干流径流变化特征，并指出：整体上干流径流处于稳定增长趋势，且突变年份在2004年左右；CHENG等在研究气候变化和人类活动在不同潜在蒸发条件下引起的径流变化特征时，提出未来的研究中应充分考虑气候变化和人类活动对径流变化的影响；李芳等人对黑河径流变化成因的分析表明，人类活动对流域径流变化影响大于气候。

1 主要研究方法

1.1 STL季节分解

STL是一种可有效处理非线性、局部趋势的非参数统计方法,是基于Loess将时间序列分解为加性变化分量的滤波季节分解方法。该方法由CLEVELAND等提出，可以稳健处理任何类型的季节性数据，且在数据出现少数异常观测值时，不会影响模型对趋势周期和季节性因素的估计。STL将时间序列X分解为季节项、趋势项和随机项，分解表达式为

式中，X是月径流量(亿m3);XS是分解出的高频季节项；XT是分解出的低频趋势项；XR是随机分量。在实际分析中，序列的长期趋势和无固定周期的循环波动，以及季节性变化和具有固定周期的循环波动难以被准确地区分。趋势项包括长期趋势和无固定周期的循环波动，季节项即为所有具有稳定周期的循环波动。除季节项和趋势项之外的其他因素的综合影响为随机项。随机波动项属于不确定性成分，常表现出异方差性。

STL方法由内循环和外循环两部分组成。在内循环中，长期趋势项和季节项通过连续的局部加权曲线提取。内循环完成之后，在外循环中使用 Robust 估计，即偏离正常值越多的异常值在估计中所占的权重越低，在最大程度上，保证在后续内循环中长期趋势项和季节项免受异常值的影响。当所有成分提取完成且验证收敛后，迭代程序停止。

STL季节分解过程中的季节项平滑参数(ns)和趋势项平滑参数(nt)是最重要的两个参数，分别控制提取趋势项和季节项时取样窗口的宽度。一般情况下，季节项平滑参数与数据的时间跨度相关，一般取略大于时间跨度的奇数。趋势项平滑参数的取值满足

式中，np表示一个周期内样本数量。

1.2 小波奇异点识别原理

1.2.1 突变点检测

小波变换卷积形式可表示为

式中，WTax(t)可看成是信号x(t)经过冲激响应ψa(t)系统后的输出。

信号的突变性检测是先对原信号在不同尺度上进行“磨光”,使信号平滑化，再检测“磨光”后信号的一阶(二阶)导数极值点(过零点)。常用的磨光函数(平滑函数)θ(t)满足

若平滑函数二次可导，其一阶、二阶导数表达式ψ(1)(t)和ψ(2)(t)分别为

且ψ(1)(t)和ψ(2)(t)满足小波的容许性条件

运用平滑函数导数分别与信号做卷积可以表现不同尺度a下突变信号的突变特征，表达式为

信号与光滑函数的卷积将信号在突变点处的部分高频信息滤去，使信号在突变点处变得光滑，取平滑函数的一阶导数作为小波函数时，小波变换系数WT(1)ax(t)的模局部极值点即为信号的突变点，其对应二阶导小波变换系数WT(2)ax(t)的模过零点，这就是小波变换奇异点检测原理。由于|WT(2)ax(t)|的过零点易受到噪声的干扰，难以区分突变点变化程度，一般不用作信号的突变点检测。而|WT(1)ax(t)|的极大值点对应信号快速突变点，极小值点对应信号缓慢突变点，用于突变点检测更具有优越性。

1.2.2 趋势识别

水文时间序列随着时间的增长，会在相当长时期里向上或向下缓慢地发生变化，这种有一定规则的变化称之为趋势。当趋势只出现在序列中一段时期时，叫局部趋势；出现在序列整个过程中的，叫整体趋势。小波分析具有非常强大的多尺度分辨功能，能够识别水文时间序列中各种高低不同的频率成分。不同尺度下的低频(近似)成分表示水文时间序列在该尺度下的变化趋势。因此，通过小波变换，得到不同尺度下水文时间序列的近似系数，由近似系数的变化过程可分析出水文时间序列的发展趋势。随着尺度的增加，由于时间分辨率降低，噪声的影响减小，信号的发展趋势会越来越清晰，从而使信号的发展趋势显现出来。

利用小波分析识别趋势成分的步骤为：(1)选择合适的小波函数，利用离散小波变换对水文时间序列进行小波分解，得到各种尺度下小波分解近似系数；(2)对各尺度下近似系数序列进行单支重构，得到不同尺度下的低频序列(针对于“二抽取”离散小波变换);(3)对不同尺度下的近似系数序列进行分析，判断其趋势性质，进一步运用传统方法定量确定趋势成分。

1.3 二次样条小波多孔算法

多孔算法也称之为稳定小波变换(Stationary Wavelet Transform),它不需要对采样点数据进行“二抽取”,只需对滤波器各相邻点间进行补零插值,再与定频信号做卷积。

小波分析中，尺度函数?(t)、小波函数ψ(t)以及滤波器系数h(低通滤波器)和g(高通滤波器)满足关系式

式(10)是双尺度方程，尺度函数满足∫∞-∞?(t)dt≠0, ‖?(t)‖=1;式(11)是小波方程，满足∫∞-∞ψ(t)dt=0。对信号f(t)进行二进小波变换

原始信号x(t)通过小波变换可以分解为细节系数{dj}1≤j≤J和近似系数aJ,初始条件下取a0(n)=x(t),离散二进小波变换多孔算法细节系数和近似系数可表达为

2 研究区域及数据

图1 区域水系

3 结果与讨论

3.1 季节分解

运用STL方法对径流数据进行时间序列分析，可以剔除干扰因素，从而更为准确地反映各水文站径流量变化趋势。图2是黑河四个水文站的月径流量STL季节分解图，由于时间序列为月数据，因此np取12,各水文站季节项平滑参数取值如表2所列。

图2 各水文站月径流量季节分解示意

3.2 基于小波多孔算法的径流变异规律分析

图3(a)—(d)分别为札马什克、祁连、莺落峡和正义峡站，小波分解近似系数和细节系数图。在进行渐进趋势分析时，在小尺度上，高频信息占主要成分，不易分析整体趋势变化特征。因此要结合大尺度(即j=J)下的近似系数，分析时间序列整体趋势变化过程，且尺度越大细节成分越少，整体趋势越明显。相反，在进行突变趋势分析时，大尺度上低频信息占主要成分，只能大致确定各水文站突变点发生时间，要结合小尺度上细节系数模极大值，才能进一步确定突变点发生年份，以减小低频信息对检测结果的影响。

图3 各水文站年径流量小波分解示意

根据图3(a)—(d)年径流量小波分解近似系数图(左侧),可得各水文站年径流量趋势特征，札马什克水文站年径流量整体上呈上升状态；局部趋势有先下降后波动上升的特点，并在20世纪70年代末达到最小。祁连水文站年径流量整体呈现波动上升趋势，随着分解尺度的增加，整体趋势越来越明显；局部上，前期呈上升趋势，在80年代中后期有一个突然下降趋势，随后又呈稳定上升趋势，且20世纪90年代前上升趋势比90年代后明显。莺落峡水文站径流量，整体也呈上升趋势；局部上，20世纪50年代到70年代呈下降趋势，随后波动上升，且在20世纪70年代到90年代上升趋势比21世纪初上升的快，最小径流量在1970年左右。正义峡水文站年径流量整体趋势不明显，但是在20世纪90年代初期，径流量在逐渐下降，且80年代下降速率相比60—70年代快，20世纪90年代中期径流量又突然上升，在1992年左右达到最小。

不同尺度下，小波分析细节系数图如图3右侧所示。札马什克站年径流量，大尺度下显示极值点主要发生在1965年左右、1990年和1995年左右、以及2001—2005年期间，结合小尺度，可确定可能突变点主要发生1958年、1967年、1989年、1996年和2005年。祁连水文站，大尺度下小波细节系数模极大值主要在1975年、1980年、1990年左右和2000年左右，结合小尺度分析结果，确定突变点发生时间在1977年、1982年、1989年和2002年左右。莺落峡水文站，综合分析不同尺度下小波分解细节系数模极大值，确定其突变点发生时间大概在1958年、1967年、1975年和1989年及2003年左右。正义峡水文站小波分解细节系数模极大值，主要集中在1974年左右、1983年左右、1989年和2001年左右。

3.3 基于传统方法变异规律分析

3.3.1 渐进趋势识别

季节调整后的黑河上中游各水文站实测年径流量线性回归和五年滑动平均图，如图4中虚线和点划线所示。为方便后期对比五年滑动平均与小波分析方法，将重构后小波分析近似系数曲线图在图4以直线形式显示出来，分解尺度为j=J。

图4 各水文站年径流量趋势

根据图4中虚线(线性回归)可以看出：札马什克水文站年径流量，整体呈现上升趋势，但是没有通过0.05显著性F检验，上升趋势不明显(|r|=1.591<4.027);祁连水文站年径流量，整体呈现显著的上升趋势(|r|=8.854>4.079);莺落峡水文站年径流量，整体呈上升的趋势特征，线性回归分析中相关系数|r|=7.090>4.016超过了0.05显著性F检验标准，上升趋势显著；正义峡水文站年径流量，整体有下降趋势，F检验相关系数|r|=3.387<4.052显著性水平未达到0.05,总体趋势不明显。

图4中五年滑动平均结果：札马什克水文站，径流量在20世纪70年代波动下降后上升，并在20世纪90年代出现最低点，但实际径流量在20世纪70年代初达到最小；祁连水文站，年径流量呈波动上升，在20世纪80年代末达到最大值；莺落峡水文站，年径流量在20世纪80年代前持续下降，并在70年代初达到最小，随后呈波动上升；正义峡水文站，年径流量呈持续缓慢下降变化特点，并在20世纪70年代到达最低值，随后波动上升，上升趋势较明显。

3.3.2 突变点检测

图5 各水文站年径流量M-K突变检验

根据表3显示，黑河各水文站的突变点发生时间：札马什克水文站主要分布在20世纪60年代和21世纪初期；祁连水文站突变点主要发生在20世纪70年代末80年代初和21世纪初；莺落峡水文站突变点主要发生在20世纪60年代末、70世纪中期和21世纪初；正义峡水文站年径流量突变点主要发生在20世纪90年代初和21世纪初。

为了确保结果的可靠性，运用两种方法，确定各水文站突变点年份。图5(a)为札马什克M-K突变检验图，图中UF与UB曲线在α=0.05显著性水平线内交点为1958年和2006年，即可能存在的突变点为1958年和2006年，结合表3滑动T检验结果，确定实测年径流量在2006年发生突变，这项结果与丁霞[46]相同。图5(b)为祁连水文站M-K突变检验图，UF与UB曲线在α=0.05置信水平线内交点发生在1982—1984年和1992年、1994年，结合滑动T检验结果，确定祁连水文站年径流量在1982年发生突变。图5(c)为莺落峡水文站M-K突变检验图，图中UF与UB曲线在α=0.05置信水平线内交点在2003年、2004年、和2005年，而滑动T检验结果显示，莺落峡水文站年径流量突变点主要发生在20世纪60年代末、70世纪初和21世纪初，则确定可能突变时间在2004年左右。图5(d)为正义峡水文站M-K突变检验图，图中UF与UB曲线在α=0.05置信水平线内交点在1969年、1972年和1973年左右，滑动T检验结果显示(见表3)可能突变点主要发生在20世纪末和21世纪初，则M-K突变检验结果与滑动T检验结果存在差异，不能确定突变时间。

3.4 对比分析

本节将对小波分析结果与传统方法变异规律分析结果进行对比，以验证小波分析对于时间序列变异规律分析的合理性和可靠性。由于在进行时间序列突变分析时，小波分析只能大致确定突变点，若要进一步确定突变时间，需要结合其他突变点检测方法。因此，本节将小波分析模极大值检测结果与M-K突变检验结果相结合，确定各水文站突变点发生时间，并与3.3节中突变分析结果进行对比。

3.4.1 渐进趋势识别

根据图3中小波分解近似系数图可以看出，随着尺度增加，小波分解高频信息逐渐被过滤，时间序列渐变趋势特征越来越显著。在大尺度j=J时，高频信息基本已完全过滤。对比图4中小波曲线与五年滑动平均曲线，近似系数曲线变化过程与五年滑动平均过程基本一致，且相比五年滑动平均，大尺度下的小波近似系数更能从整体上反映序列趋势特征。五年滑动平均虽然削弱了细节变化，但小波分析近似系数可以在不同尺度下对细节信息进行过滤，不同尺度下小波近似系数，可以在不同程度上显示时间序列的渐进趋势特征。同时在图4中还可看出，小波分解近似系数整体趋势与线性回归趋势相同，但相比线性回归，小波分析不能分析显著性特征。因此，通过分析小波分解近似系数变化过程，可以同时判断时间序列的整体和局部趋势特征。而且根据3.2和3.3趋势识别结果，可以看出，小波分析对于局部趋势识别更具有可靠性。

3.4.2 突变点检测

根据图5中M-K突变检验结果和图3小波分解细节系数结果，可进一步确定各水文站突变点发生时间：札马什克突变点主要发生在2006年；祁连突变点主要发生在1982年；莺落峡水文站年径流量突变在2003年；正义峡年径流量在1973年发生突变。而3.3节中并未检测出正义峡年径流量发生突变。

4 结 论

(1)在渐变趋势分析中，不同尺度下的小波近似系数，可以在不同程度上显示时间序列的渐进趋势特征，小尺度下小波分析结果与五年滑动平均相同，大尺度小波近似系数可反映时间序列整体趋势特征，这与线性回归趋势一致。

(2)对于渐变趋势分析，小波分析法可以有效地对各水文站的整体和局部趋势做出判断，但是不能体现出趋势变化的显著性特征。

(3)在突变分析中，小波分析不像M-K突变检验法一样会产生漏点现象，也不像滑动T检验法一样，要求时间序列满足一定的分布特征，同时不会引入过多伪突变点。

(4)小波分析可以作为突变点检测初步诊断工具，若需准确确定各水文站突变点发生时间，需将小波分析与其他检测方法相结合。

---

水利水电技术（中英文）

水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也 道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。