基于年径流序列的五种趋势检测方法性能对比

姜 瑶，徐宗学，王 静

（1.北京师范大学 水科学研究院，北京 100875；2.城市水循环与海绵城市技术北京市重点实验室，北京 100875；3.西藏自治区水文水资源调查局，西藏 拉萨 850000）

1 研究背景

在地球气候系统显著变化及人类活动影响下，全球水循环过程正发生着不同程度的改变[1-3]，认识和掌握变化环境下的流域水循环演变规律具有重要意义。趋势变化是水循环过程的重要变化特性之一，利用趋势识别技术理解和认识水文要素在不同时间尺度上的变化规律是变化环境下水循环演变研究的重要内容[4-6]。

趋势识别的主要内容是在所关心的时间尺度之上确定水文时间序列趋势的具体形状以及评估趋势的显著性[7]。目前，针对全球不同地区、不同流域的水文气象要素，国内外研究者已开展了大量趋势识别研究，应用或发展了多种趋势识别方法[4，8-12]。常用的趋势识别方法有线性回归方法、累积距平曲线、滑动t检验等参数统计法，以及Mann-Kendall（MK）秩次检验、Spearman 秩次检验、Sen’s斜率估计等非参数统计方法。这些方法多是时间域内的趋势分析[7]，其中线性回归方法和MK方法在水文气象序列的趋势识别中应用最为广泛。线性回归方法根据水文观测数据与其对应时间之间的回归系数判断水文时间序列的线性趋势及变化率，然而由于自然系统的非线性、非平稳性，水文时间序列在长时间尺度上往往表现出更为复杂的非单调变化趋势[13]，因此线性趋势常常不符合真实的变化过程。MK方法对数据的分布形态没有特定要求，且不受部分数据缺测的影响[14-15]，特别适用于水文气象数据的趋势检测，但在实际应用中MK分析结果会受到序列自相关性、序列长度、周期波动以及趋势幅度等因素影响，同时也无法给出趋势变化的具体形状[16-17]。

为更准确地识别水文时间序列的多时间尺度变化特性，一些时频分析方法开始广泛应用于水文气象要素的趋势分析中，其中离散小波变换（Discrete wavelet transform，DWT）和经验模态分解（Em⁃piric mode decomposition，EMD）方法最为常用。在DWT方法应用方面，Pathak 等[18]基于DWT法分析了美国中西部地区季节温度、降水和径流的变化趋势；Palizdan 等[19]利用DWT法对Langat 河流域降雨的长期趋势进行了分析；Sang 等[20]提出了基于DWT法的非单调趋势识别及显著性检验方法，并应用该方法估计了气温与潜在蒸散发的变化趋势。DWT法能更有效地识别水文序列的非单调变化趋势，但其也受到小波函数选择和分解水平选择等因素的制约[20]。EMD法依据时间序列本身的特性对序列进行分解，具有直接、后验和自适应等突出特点[21-22]，在水文要素的周期及趋势识别方面都有着广泛的应用[16，23-26]。

由于趋势识别方法所基于的原理不尽相同，且不同方法自身均具有局限性，在具体应用时往往存在分析能力和分析结果上的差异[5-6]。为了获得准确的趋势识别结果，一些研究综合应用多种方法进行对比分析，如Tabari 等[10]基于线性回归、MK检验、Pettitt检验和Sen’s斜率估计等方法分析研究了伊朗地区的气温、降水及潜在蒸发量等要素的演变特征；Nalley 等[27]联合利用MK方法和DWT方法估计了月尺度、季节尺度以及年尺度平均径流量与降雨数据的变化趋势；Kahya 等[28]采用多种非参数估计方法（MK检验、Sen’s斜率、Spearman检验和Seasonal Kendall检验）对土耳其26个流域的月径流变化趋势进行了分析。但目前关于趋势的定义不甚统一，对趋势的认识和识别仍具有局限性，实际研究中很难判断哪种方法得到的检验结果更为准确和可靠，水文时间序列的趋势识别仍然存在诸多不确定性[16，29]。

为此，本文选取5种原理不同的趋势识别方法（线性回归方法、累积距平方法、MK方法、EMD方法和DWT方法），基于已知成分的人工生成序列和多个流域的实测年径流序列，综合对比5种趋势识别方法的性能差异，探讨不同方法对年径流序列趋势识别的适用性，以期提高对常用水文要素趋势识别方法的认识，为水文要素的趋势分析及其方法选择提供参考。

2 方法介绍

2.1 线性回归法线性回归法是构建水文变量x与其时间t之间的一元线性回归方程：

式中：x（t）为样本量为n的某一水文变量；t为x（t）所对应的时间；a为截距；b为斜率。

a和b可通过最小二乘法求出，其中，回归系数b 即可反映水文变量的变化速率，其正负表示水文变量的增减变化趋势。

回归系数b的显著性可通过t检验进行判断，构造T 统计量：

2.2 累积距平法累积距平法是一种由曲线直观判断变化趋势的方法。对于样本长度为n的水文变量x（t），其某一时刻t的累积距平值cm表示为：该方法的核心是当数据持续大于平均值时，累积距平值增加，曲线呈上升趋势，反之则呈下降趋势[6]。根据曲线的起伏变化，可以判断时间序列的长期演变趋势及发生突变的大致时间。

2.3 MK法Mann-Kendall（MK）检验[14-15]是基于原始时间序列秩次的一种非参数统计方法，其无需数据服从一定的分布，同时不受少数异常值的影响，对于非正态分布的水文气象数据分析具有突出的适用性，在水文要素的趋势识别中得到广泛应用。对于一个长度为n的水文时间序列x（t），原假设H0认为原数据序列是n个独立同分布的随机变量，没有趋势存在。MK检验的统计量按下式计算：

其中，

统计量S的正（负）值表示序列的上升（下降）趋势。当数据长度n＞8时，S近似服从标准正态分布[14-15]，其期望E（S）和Var（S）可表示为：

式中：m为数据序列内含有的相等数据的组数；kl为第l组内相等数据的个数。

标准化的统计检验量Z 按下式计算：

2.4 EMD法经验模态分解（Empirical model decomposition，EMD）方法的本质是设法识别并分离原始序列固有的多尺度振荡特征，从而形成一组本征模态函数（Intrinsic mode function，IMF）分量和一个残余项，各IMF表征了原序列不同时间尺度的振荡变化，残余项一定程度上表现原序列的总趋势[21]。EMD法的分解是基于序列时间尺度的局部特性，具有自适应性，同时物理意义明确，特别适用于非线性非平稳序列的分析处理[23]。EMD 通过自适应均值筛选过程实现，详细描述见文献[21，23]。经过分解后的原序列可以表示为：

式中：N为从原序列x（t）分解得到的IMF总个数；Cj为第j个IMF分量；rn为残余项，是一个单调函数或无法继续分离IMF分量的唯一极值函数。

各IMF分量的显著性检验可利用Wu 等[23]所建立的基于均匀分布白噪声统计特征的显著性检验方法。

2.5 DWT法令L（2R）是定义在实轴上、可测的平方可积函数空间，对于时间序列 f (t)∈L2(R)，其离散小波变换表示为：

其中，

式中：Wf(j,k)为离散小波系数；ψ*(t)为母小波ψ(t)的复共轭函数；a、b均为常数；j为分解水平；k为时间位移因子。

在实际分析中，常使用二进制离散小波变换，即设a=2和b=1[30-31]，将ψ(j,k)

t表示为：

通常，分解水平j的理论最大值M可根据时间序列f（t）的长度L求得[32]：

当选定合理的小波分解水平之后，可应用离散小波分解和重构将水文序列分解为不同分解水平上的子序列：

各分解水平下子序列的显著性检验参考Sang 等[34]提出的基于噪声能量分布的显著性检验方法进行。

3 案例分析

3.1 数据描述综合使用人工生成序列和实测径流序列对比5种趋势检测方法的性能差异。3组人工生成序列分别考虑无趋势、上升趋势和下降趋势，并在此基础上叠加周期项和随机波动后生成（图1）。各序列长度均为200，其基本组成见表1。实测序列采用3个不同水文站点的年径流量序列进行对比分析，分别为雅鲁藏布江干流奴下水文站年径流量实测序列（RN序列）、黑河流域干流莺落峡水文站年径流量实测序列（RY序列）和黄河干流上游石嘴山水文站年径流量实测序列（RH序列），序列长度均为60年（1956—2015年）（表1）。从直观上看，3组实测数据（RN序列、RY序列和RH序列）分别呈不明显变化、明显下降和明显上升趋势（表1）。为便于对比分析，人工生成序列和实测年径流序列均采用平均值及标准差进行标准化处理。

图1 3组人工生成序列实际变化及其真实趋势

表1 3组人工生成序列与3组实测径流序列的主要信息

表2 人工生成序列的线性回归法趋势检测结果

图2 3组人工生成序列的MK分析结果

3.2 人工生成序列趋势识别结果表2为人工生成序列的线性回归法趋势检测结果，线性回归结果表明，3组人工序列（S1、S2和S3）的平均变化率分别为-0.0046、0.0077和-0.0099，总体上分别呈减小、增加和减小的显著变化趋势。

MK法得到3组序列的趋势检测结果见图2，统计量UFi=200分别为-3.632（下降）、6.439（上升）和-8.327（下降），平均变化率分别为-0.0035、0.0083和-0.0104（图2），与线性回归法结果一致。然而，两种方法下S1序列的检验结果与其真实趋势不一致，而S2和S3序列的趋势检测结果虽然正确，但仅是序列的总体变化趋势，未反映其真实的非线性变化过程。

图3为3组人工生成序列的累积距平曲线，其呈现了3组人工序列的趋势变化过程及其拐点。图3与图2中MK法所得结果相似，但由于受序列周期变化的影响，其描述与序列的真实变化均有差异，尤其S1序列，其累积距平值受异常点和周期变化的影响而与实际变化差异较大。

图3 3组人工生成序列的累积距平曲线

图4 3组人工生成序列的EMD分析结果及其显著性

图5 3组人工生成序列的DWT分析结果及其显著性

应用EMD法分别对3组人工序列进行分解，得到其不同时间尺度下的变化特征曲线（即IMF1-5分量和残余项），其中残余项的变化代表序列在最大时间尺度下的趋势。图4描述了3组序列分解后残余项与IMF5分量的变化特征及其显著性检验结果。EMD法识别出S1序列总体表现为减小但趋势和周期变化均不显著（图4（a）（d）），但由于序列随机波动及个别异常点的影响，所得S1序列变化仍与其真实情况不完全一致。对于S2和S3序列，EMD法所得趋势项与其真实趋势基本一致，但识别出的序列周期变化仍存在不平滑的问题，与序列的实际变化存在一定偏差（图4（b）（c）（d））。

图5给出了基于DWT方法得到的不同分解水平下的子序列变化及其显著性检验结果，可分别代表原序列不同尺度下的波动特征。对比EMD法，DWT法检测出S1序列的趋势不显著但周期变化显著，与S1序列的实际变化基本一致（图5（a）（b））。同时，DWT方法识别出了S2序列和S3序列的线性变化趋势及其真实的非线性变化趋势，且趋势均显著，得到的序列变化曲线也与其实际变化更为接近（图5（b）（c））。

总体上，累积距平法只能描述序列的变化过程及拐点，无法判断序列的总体变化趋势，且结果受周期变化和异常点的影响较大，很难得到准确的变化拐点。线性回归和MK法能给出序列的总体变化趋势、平均变化率及其显著性，对于变化趋势显著的序列具有较准确的检验结果，但对周期变化明显而趋势不显著的序列，其识别结果很容易受周期变化的影响而出现偏差。EMD和DWT方法能对原始序列的不同成分进行有效分解，从而较准确地识别出原序列的周期与趋势变化，因此更适用于具有复杂变化特征的水文时间序列趋势分析，其中DWT法的趋势识别效果较EMD法更优。

3.3 实测序列分析分别应用上述5种方法对3个流域的实测年径流序列（RN、RY和RH）进行分析，进一步对比5种趋势识别方法在径流变化趋势分析中的性能差异。3组径流序列分别具有不同的趋势和周期变化特征（表1），表3统计了基于5种方法得到的3组实测径流序列的趋势检测结果。线性回归法和MK法的结果均表明，RN、RY和RH序列在60年尺度上的趋势分别为上升但不显著、显著上升和显著下降（图6和图7）。两种方法计算得到的平均变化速率也相差不大，线性回归结果为0.0015/a、0.0268/a和-0.0234/a，MK法为0.0015/a、0.0281/a和-0.0191/a（表3）。同时，由图7可知，RN序列在1964年和2000年可能存在突变点，这与实测RN序列在1964年和1998—2000年左右的峰值点相对应，而RY和RN序列分别在1975年和1978之后变为显著增加趋势（图7）。相较于线性回归法，MK方法定量化地显示出更多序列变化的细节。累积距平曲线显示，RN序列总体有2个明显的趋势变化拐点，分别在1964年和1998年，而对于RY和RH序列，其累积距平曲线变化的拐点分别在2000年和1985年左右，该结果进一步验证了MK方法所得突变点（图8）。

表3 3组实测年径流序列趋势检测结果

图6 3组实测径流序列的线性回归法趋势分析结果

图7 3组实测径流序列的MK分析结果

图8 3组实测径流序列的累积距平曲线

EMD和DWT方法的趋势识别结果与前述3种方法（线性回归、MK法和累积距平法）所得结果差异较大，尤其对于RN序列。EMD法识别结果显示，在60年尺度上，RN和RH序列呈显著下降趋势，RY序列为显著增加趋势（表3），同时RN序列在第5分量上（IMF5）表现为显著非单调变化（图9）。与EMD法相似，DWT法将原序列分解为6个水平下的子序列，其中各实测年径流序列在分解水平6下的子序列变化趋势与EMD结果总体一致（图9和图10）。然而显著性分析表明，RN序列的单调下降（D6）不显著，在第5分解水平上的子序列为确定成分，表现出非单调变化趋势；RY序列的单调增加趋势（D6）显著，为确定的趋势成分；而RH序列的单调减小趋势（D6）和弱非单调变化趋势（D5）均显著（图10和表3）。

图9 3组实测径流序列的EMD分析结果

图10 3组实测径流序列的DWT分析结果

综上所述，累积距平法直观地描述了年径流量的变化特征，识别出的变化拐点与MK法结果相似。对于RY和RH序列，由于其趋势变化非常明显，因此5种方法中除累积距平法外，均识别出了相似的变化趋势且趋势显著（表3）。然而，对于RN序列，DWT法和EMD法结果均表明该序列存在显著的非单调变化趋势，受其影响，基于线性回归法和MK法得到的序列整体表现为增加趋势但不显著的结果，这表明仅线性回归或MK方法很难准确地识别出RN序列的真实变化趋势。

3.4 讨论基于上述分析，5种趋势识别方法的特点及适用性如表4所示。从特点上看，累积距平法、线性回归法和MK法均为时域内分析方法，能描述出序列整体随时间的变化特征，但也易受序列内在复杂变化成分（如不同尺度周期变化）的影响。其中，累积距平法仅能反映序列某一区间的变化特征（上升或者下降），线性回归法可描述出时间序列总体的线性趋势，MK法可定量给出序列的总体变化趋势、拐点及显著性，但无法确定趋势的具体形状。EMD法和DWT法均为时频域内分析方法，可对序列不同时间尺度的变化特征进行有效分解，同时从时域和频域揭示时间序列的局部特性，从而较准确地识别出序列的真实趋势成分及其变化的具体形状（线性和非线性）。

表4 5种趋势识别方法特点及适用性

从适用性上看，对于变化趋势明显的年径流序列（如RH序列和RY序列），5种方法中除累积距平法外，均得到了较一致的趋势识别结果，适用于演变特征相对简单且趋势明显的年径流序列的趋势检测。其中，线性回归法和累积距平法简单直观，可用于年径流序列变化趋势的初步分析，MK法与线性回归法的趋势识别结果差异不大，但MK法受序列缺测值、异常值影响较小，且能定量地给出序列趋势变化的显著性及突变点，性能优于线性回归法。但MK法与线性回归法均易受序列周期或非线性特征的影响，不适用于具有明显周期波动成分的水文时间序列的趋势识别（如月尺度径流序列），而EMD法和DWT法均能对原始序列的不同成分进行有效分解，从而较准确地识别出原序列的周期与趋势变化及其形状（线性和非线性），其性能明显优于其它3种方法，可用于具有多时间尺度变化、非线性、非平稳特性的水文时间序列分析，其中DWT法表现更优。例如在案例分析中，由于人工生成序列S1的周期变化影响，其线性回归和MK法识别结果均与真实趋势不符，而EMD和DWT法不仅识别出了S2和S3序列真实的非线性变化趋势，也较好地识别出了S1序列的趋势和周期，其中DWT法的结果与真实变化相符更好，而EMD法仍会受个别异常值的影响而与真实变化产生一定偏差。因此，对于演变机理及趋势变化复杂（不明显或非单调）的年径流序列（如RN序列），仅用线性回归或MK方法很难得到准确的趋势判断，有必要在序列分解基础上进一步分析，此时EMD法和DWT法相对更适用。

5种趋势检验方法是水文要素趋势分析的常用方法，在国内外不同流域的降水、径流、蒸散发等的变化规律分析中均有应用。研究基于具有不同变化趋势的人工序列和实测序列，评价了各趋势检验方法对不同趋势特征的识别能力，所得结论也可以从以往相关研究中得到验证[6，16-17]，可为不同流域水文要素的趋势分析及其方法选择提供一定的参考。但同时，由于本文仅以3组年径流序列为实例，从应用角度对比分析了不同方法在水文时间序列趋势分析中的性能差异，对不同方法的性能评价存在局限性，有待进一步深入研究。

4 结论

本文选取5种基本原理不同的趋势识别方法（线性回归、累积距平、MK、EMD和DWT方法），基于3组已知成分的人工生成序列对比了5种方法在趋势检测方面的性能差异，并以3组具有不同趋势变化特征的实测年径流序列为例，进一步探讨了不同方法在水文时间序列趋势识别中的适用性。

实例分析结果表明，累积距平法、线性回归法和MK法能在时域内描述出序列整体的变化特征，但也易受序列内在复杂变化成分（如不同尺度周期变化）的影响，适用于趋势单调且无明显周期波动的水文时间序列，对于演变特征相对简单且趋势明显的年径流序列具有较准确的趋势检测结果。其中，累积距平法仅能表明序列某段变化趋势，且结果对序列的异常值较敏感；MK法与线性回归法的趋势识别结果差异不大，但MK法受序列异常值影响较小，且能直观、定量地给出序列趋势变化的显著性及突变点，性能优于线性回归法。EMD法和DWT法可基于时频分析准确地识别并分离出序列自身含有的不同尺度波动成分，从而避免序列周期、非线性等变化特征对趋势识别的影响，对于具有多时间尺度变化特征及非平稳特性的水文要素序列的趋势检测性能更优，其中DWT法所得结果与真实变化更相符。因此，对于演变机理及变化特征复杂的年径流序列，仅用线性回归或MK方法很难得到准确的趋势判断，有必要在序列成分分解基础上进行进一步分析，以避免序列趋势程度、周期、非线性等特征的影响，而其中DWT方法应为首选方法，EMD方法则是另一种较为可行和值得推荐的方法。