基于极点对称模态分解的三江源径流时空演变规律

李凯 任卓钰 王永强 许继军

摘要： 气候变化和人类活动加剧改变了地球的水循环系统，因此，重新认识和掌握流域水循环演变规律具有重要意义。利用极点对称模态分解（ESMD）方法同时从趋势、突变、周期3个方面对三江源径流变化进行分析，并与现有的检测方法进行对比。结果表明：① 黄河源地区由于降水增加趋势不明显，随着气温升高， 蒸散量增大，导致年径流呈现不显著的下降趋势;长江源、澜沧江源径流由于降水和冰川融雪的增加，导致年径流呈现显著增强的趋势。② 黄河源近61 a的径流过程在1990～1994年、2008年左右发生明显突变，长江源径流过程在2006年左右发生显著突变，澜沧江源径流过程在20世纪90年代和2008年左右发生突变。③ 黄河源近61 a径流过程存在12～15 a和30 a的周期性变化;长江源近60 a径流过程存在13～15 a的周期性变化;澜沧江源近49 a径流过程存在15～16 a和28 a的周期性变化，这些周期性变化主导着三江源年径流变化特性。

关 键 词： 径流变化; 周期分析; 气候变化; 极点对称模态分解; 长江; 黄河; 澜沧江

中图法分类号：  P333

文献标志码：  A

DOI： 10.16232/j.cnki.1001-4179.2021.08.013

0 引 言

近年来，随着全球气候发生显著性变化和人类活动越来越强烈，地球的水循环系统正发生着不同程度的改变，重新认识和掌握变化环境下的流域水循环演变规律具有重要意义。其中周期、趋势和突变变化是水循环过程的重要变化特性之一，利用统计分析方法理解和认识水文要素在不同时间尺度上的变化规律是变化环境下水循环演变研究的重要内容。目前，国内外学者使用的常规方法有：周期分析方法主要包括周期图法、功率谱法、Morlet小波变换、经验模态分解等;趋势分析方法主要包括滑动平均、Mann-Kendall检验、R/S分析等;突变分析方法主要包括Mann-Kendall突变检验、Pettitt突变检验、累积距平法、有序聚类法等。然而，上述方法均只能从单一的方面量化探讨气候变化对径流变化过程的影响。而极点对称模态分解（Extreme-point Symmetric Mode Decomposition，ESMD）方法能够同时从趋势、突变和周期方面对非平稳、非线性的时间序列变化进行分析[1]。现阶段，该方法已经应用于新疆温度变化分析[2]、北京的降水序列变化分析[3]以及江西干旱周期分析[4]中。将该方法应用于径流序列分析时，径流数据虽然具有较强的周期性，但在每半个周期上并不关于局部中点对称，已有的ESMD方法不仅无法客观代表数据的趋势，反而会产生较大的误差。姜铄等[5]对ESMD方法进行改进，将不具备代表性的中点替换为能反映每半个周期数据情况的局部均值后，就可以用来对径流时空演变规律进行分析。

三江源区包括长江源、黄河源以及澜沧江源，在整个中国以及亚洲东部的水资源保护、径流调节以及物种多样性的保护方面起着重要的作用。由于全球气候变化加剧和人类活动的影响，三江源地区径流演变过程可能因降水量增加、温度升高、蒸散量减少等异常现象而发生显著变化[6]。张士峰等[7]从趋势、突变分析三江源1965～2004年的径流变化特征，发现径流呈现减少的趋势，并且在1994年发生突变。商放泽等[8]分析了1956～2012年三江源的径流演变规律，发现长江源和澜沧江源径流呈现显著上升的趋势，黄河源径流呈现下降的变化趋势。刘希胜等[9]单独研究黄河源2012年以前的径流趋势变化，结果表明黄河源的径流在减少。朱海涛[10]研究长江源2012年以前的径流演变规律，发现长江源的径流呈现显著增加的趋势，在1961，1967年受到气候和人类活动的影响，径流发生突变，此外径流存在8 a和24～28 a的显著周期。王珂等[11]对澜沧江源1960～2010年的径流变化进行分析，发现澜沧江源地区年径流并未呈现显著性的变化趋势。

綜上所述，现有关于三江源径流演变规律的研究都是在2012年以前，而欠缺对近几年三江源地区径流演变的研究，此外尚未采用一种综合的方法同时从趋势、突变和周期方面分析三江源的径流演变规律。因此，为更好认识和掌握三江源地区的水文循环规律，分析三江源径流时空演变规律，本文以三江源区直门达水文站1957～2016年（长江源）、唐乃亥水文站1956～2016年（黄河源）、昌都水文站1968～2016年（澜沧江源）实测径流序列为研究对象，利用ESMD方法同时从周期、趋势、突变3个方面研究三江源地区的径流时空演变规律，从而为径流预报和水资源管理提供切实可靠的信息。

1 研究区域概况

三江源位于青海省的南部地区，地处青藏高原腹地（31.39°N～36.12°N，89.45°E～102.23°E），是长江、黄河和澜沧江的发源地，同时也是中国淡水资源的重要补给地。三江源流域面积31.26万km2，其中长江源（直门达水文站以上）流域面积13.77万km2，占三江源区总面积的44.0%;黄河源是指黄河流域唐乃亥水文站以上区域，流域面积12.2万km2，占源区总面积的39.0%;澜沧江流域是指昌都水文站以上的流域，流域面积5.29万km2，占源区总面积的16.9%（见图1）。三江源区属于典型的高原大陆性气候，年平均气温-5.6～3.8 ℃，年降水量变化较大，介于200～800 mm之间。三江源地区是中国典型的生态脆弱区域，对气候变化和人类活动非常敏感。

2 数据与方法

2.1 研究资料

本文径流数据分别源于黄河源控制站唐乃亥水文站、长江源控制站直门达水文站以及澜沧江源控制站昌都水文站;降水数据均来自中国气象数据网，数据长度与径流保持一致，每个流域的降水数据采用泰森多边形法计算得流域面平均雨量。数据具体情况如表1所列。

2.2 研究方法

径流是一种非线性、非平稳的时间序列，其中包含多种周期成分、趋势成分和跳跃成分。极点对称模态分解方法可以有效地进行多种频率成分的提取及径流序列变化特征分析，并利用ESMD模态分解的极点对称插值法进行时频分析[12]。同时，选取线性趋势回归分析、M-K检验、累积量斜率变化比较方法[13]、Pettitt突变点检验以及Morlet小波变换、周期图法作对比，综合分析径流序列趋势变化、突变规律以及周期变化。

2.2.1 极点对称模态分解

极点对称模态分解方法是著名的希尔伯特-黄变换（包括EMD和希尔伯特谱分析）的新发展，该方法保留了EMD方法的优点的同时，将EMD方法的外包线插值改为内部极点对称插值，降低了由插值带来的不确定性;通过ESMD方法的频率分布图，可以清晰地反映各模态的频率分布情况，解决了EMD方法“模态混叠”的问题。目前，ESMD方法已成功应用在生命科学、机械工程、地震学和气候等领域[12]。

2.2.2 Mann-Kendall检验

Mann-Kendall可以诊断径流序列趋势和突变点[14]，顺序结构为

UFk=  Sk-E（Sk）   Var（Sk）   （1）

根据时间序列Xn，Xn-1，…，X1的逆序排列，得到：

UBk=-UFk′k′=m+1-k   （2）

式中：k′，k=2，3，4，…，n;Sk为用时间序列构造的统计量;E（Sk）为Sk的数学期望;Var（Sk）为Sk的方差;UBk和UFk为一条随k变化的曲线;m为xi>xj的累积数。

绘制UBk和UFk曲线，置信水平α=0.05，在置信区间中 Uα/2，当|UFk|< Uα/2时，即认为趋势不显著，反之，趋势显著。如果两条曲线出现交点，则这些交点就是径流序列突变点。

2.2.3 Morlet小波分析

在时域中对小波系数的平方值进行积分以获得小波方差[15]：

Var（a）=∫+SymboleB@-SymboleB@ Wf（a，b） 2db （3）

式中：a为时间，年;Var（a）为小波方差;Wf（a，b）为小波变换系数;b为平移参数。

信号波动能量随着尺度a的分布式变化称为小波方差图，它可以直接反映径流时间序列在不同时标上的周期性变化的强度，从而确定主要时间尺度，即主要周期。

2.2.4 误差诊断

识别水文要素演变规律的方法众多，且各自方法原理不同，不同方法自身具有局限性，在具体应用时往往存在分析能力和分析结果上的差异。为了获得准确的识别结果，设置误差范围，计算不同方法间诊断误差的允许范围。具体公式如下：

δ= k X（T）  （4）

β=  x（i）-x（j）  X（T）  （5）

式中：k为误差指数;X（T）为水文时间系列总长度;δ为诊断误差范围;x（i），x（j）为不同方法的识别结果;β为诊断误差，当β<δ时，则认为在该点处诊断结果是有效的。

3 结果与讨论

3.1 趋势分析

在ESMD方法中，通过判断自适应全局均线（AGM）的变化趋势来反映出径流总体的变化趋势。由图2可知：黄河源61 a径流呈现不显著下降趋势;长江源60 a径流变化呈现显著上升的趋势;澜沧江源49 a径流呈现显著上升趋势;这些变化趋势在20世纪90年代以后表现尤为明显。ESMD的结果与线性回归趋势分析、M-K检验结果保持一致（见表2）。从线性回归趋势分析结果可以看出：黄河源区流量以每年1.05 m3/s的速率在减小，而长江源径流、澜沧江源流量分别以每年1.26 m3/s和3.55 m3/s的速率显著升高（见图3）。

分析导致三江源地区径流变化的原因，黄河源地处人口稀少的高海拔地区，人类活动对于该地区的径流的影响相对较小，导致径流减少的原因更多取决于气候变化。黄河源的径流补给来源主要是降水和冰雪冻土融水，而黄河源区的降水与径流的相关系数较大，降水的丰枯年与径流的丰枯年相对应，如图3（a）所示，这也说明黄河源降水是径流变化的主导因子。分析黄河源地区降水径流变化情况可知，黄河源的年降水量呈不显著上升趋势，增幅0.38 mm/a，径流却呈现出减少的趋势。根据已有的研究表明：三江源地区气温增幅在0.35  ℃/a[16]，气温升高导致黄河源地区蒸散发量增大，在20世纪90年代以后黄河源蒸散发以5.3 mm/a的趋势上升，其上升趋势远大于降雨的上升趋势，而20世纪90年代以后降雨类型主要以小雨为主，并且降水历时短，蒸散量增大导致流域的产流能力下降。利用累积量斜率变化比较法分析出20世纪90年代以后降水对于径流减少的贡献率仅占11.17%，其他因素导致径流减少的贡献占88.83%。

长江源区径流呈现显著上升趋势，一方面是长江源区降水与径流有较好的相关性，降水有明显的上升趋势，线性趋势增幅达0.84 mm/a，降水增多直接导致径流也相应增加，如图3（b）所示。而另一方面是因為长江源地区气温增幅显著，气温增加的同时导致冰川融水过程增强。根据已有研究表明：长江源地区冰川区产流量以每10 a 1.26×106 m3的速度增加，导致径流呈现显著上升趋势[17]。通过累积量斜率变化比较法定量分析20世纪90年代以后导致径流显著增加的因素，发现降水占主导地位，其贡献率为98.31%，其他因素占1.69%。

澜沧江源受到人类活动的影响较小，导致径流显著上升的主要影响因素是降水和气温。澜沧江源地区降水有明显的上升趋势，增幅达2.11 mm/a，并且降水与径流的分配规律保持一致，在降水丰沛的年份，流域径流也大，降水少的月份同期径流也较少，如图3（c）所示。在该地区气温影响径流主要表现在：气温升高引起冰川和积雪消融，导致河道径流增加。利用累积量斜率变化比较法定量的分析20世纪90年代以后导致径流显著增加的因素，发现降水占主导作用，其贡献率为95.34%，其他因素占4.66%。

3.2 突变分析

在ESMD方法中，利用直线插值法得到各模态分量的频率与振幅时变图。该图直观地体现了径流序列分解得到的各模态分量振幅与频率的时变性，根据图中低频、大振幅或高频、小振幅振荡的时刻，表明分解得到的模态分量中存在异常时段与频段，从而判断此刻径流序列发生突变。利用ESMD方法诊断黄河源、长江源和澜沧江源多年径流的突变位置（见图4～6），并将结果与M-K突变检验、Pettitt法检验结果相比较（见图7、表3）

将ESMD方法应用于黄河源61 a径流过程中得到3个高频、小振幅的时刻，分别是1961，1983，2008年，1个低频高振幅点1992年，说明黄河源径流序列在这些位置发生了突变（见图4）。利用公式（4）、（5）将ESMD与M-K、Pettitt法检验的结果相结合，确定黄河源径流序列在1990～1994年以及2008年左右发生突变。径流第一次突变的原因是20世纪90年代黄河源地区蒸散发量显著增强，同时该地区降水减少，导致径流量显著减少。第二次突变是由于中国在2005年开始实施三江源生态保护和建设工程，三江源地区植被情况得到明显好转，并且与工程修复前相比，唐乃亥站年径流量增加了36.9亿m3。

利用同样的方法判断长江源径流序列的突变位置，ESMD方法得到2个高频、小振幅的位置，分别是1983，2006年（见图5）。通过M-K检验判断长江源径流序列并无明显的突变位置，而Pettitt法同样检验出长江源径流序列在2006前后发生突变。3种方法综合确定长江源径流序列在2006年左右发生突变，分析发现：1957～2016年平均流量增加了39.53 m3/s。进一步分析长江源地区径流在2006年以后显著性增强的原因：一方面是实施三江源生态保护和建设工程[18]，另一方面是该地区降水的增加。

同理，判断澜沧江源49 a径流序列发生突变的位置，ESMD方法得到3个高频、低振幅的位置，1个大振幅、低频点，分别是1995，1997，2008年和1991年（见图6）。通过对比其他两种方法，发现在20世纪90年代以及2008年左右澜沧江源径流序列可能发生突变。从图3（c）中可以看出，在20世纪90年代后，澜沧江源径流减少趋势显著，1980～1989年多年平均径流量为959.23 m3/s，而1990～1999年多年平均径流量为889.59 m3/s，降低7.25%左右。已有的研究表明：20世纪90年代降水量减少、温度升高是直接导致澜沧江源地区径流发生突变的原因[11]。澜沧江源同黄河源一样在2008年左右径流序列也发生突变，导致这次突变的原因同样是2005年开始实施的三江源生态保护和建设工程[19]。

3.3 周期分析

ESMD分解得到不同频率模态分量，再利用快速傅里叶变换（FFT）计算径流序列各个模态对应的功率谱，根据幅值的大小获得径流信號的频率，从而计算出各模态分量的平均周期。不同周期性变化的识别方法都存在一定的误差，为避免单一方法造成的识别误差，将EMSD结果与Morlet小波分析、周期图结果对比。利用 ESMD方法将黄河源、长江源、澜沧江源多年径流序列进行多时间尺度分解，将非平稳的径流序列转化为平稳的各模态分量，再利用FFT周期图法分别计算年径流序列分解得到的各模态分量的平均周期。

利用ESMD分析得到黄河源区61 a径流序列具有7，15 a和30 a左右的周期性变化，通过与Morlet小波分析（见图8（a）、图9（a））、周期图结果对比（见表4），确定黄河源径流序列存在12～15，30 a左右的周期振荡（见表4），其中12～15 a左右的周期振荡最为明显，在这个周期段内，黄河源径流发生16次丰枯交替，通过丰枯交替可以判断出2016以后的一段时间内黄河源的径流是偏枯的。而在30 a的时间尺度下，黄河源径流序列经历了“丰-枯-丰-枯-丰”5次丰枯交替变化，同样可以判断在未来一段时间内黄河源径流是偏枯的。

同理，利用ESMD方法分析得到长江源60 a径流序列具有8，15，30 a周期变化，其中30 a的周期变化包含在15 a的周期变化里面;周期图法说明该径流序列存在8，15 a的周期性变化。通过Morlet小波分析（见图8（b）、图9（b））长江源存在13 a的周期性变化。通过公式（4）、（5）结合3种方法分析结果，确定长江源径流序列存在13～15 a的显著周期性变化，在这个周期尺度内径流经历了14次的丰枯交替变化，并且可以判断在2016年以后的一段时间内，长江源地区径流量会有增加的趋势，这与长江源地区的径流趋势分析结果保持一致。

将ESMD方法结合Morlet小波分析、周期图法判断澜沧江源49 a径流序列的周期性变化。ESMD方法检测出径流序列具有15，8 a的周期性变化，Morlet小波分析得到澜沧江源径流序列具有16，28 a的周期性变化（见图8（c）、图9（c）），其中以28 a左右的周期性变化最为强烈;周期图法分析得到澜沧江源径流序列具有16，24 a的周期性变化。通过公式（4）、（5）结合3种方法的分析结果可知，澜沧江源径流存在15～16 a左右的不显著周期变化和28 a左右的显著性周期性变化。从图8（c）可以看出，15～16 a这个周期尺度变化并不显著，而在28 a左右显著性周期尺度内径流经历了6次丰枯交替，最近的一次丰枯交替是发生在2015年前后，径流由枯转丰，并且在未来的一段时间内，澜沧江源径流将仍然保持偏丰状态，这与趋势分析结果保持一致。

4 讨 论

（1） ESMD方法综合实现了Morlet小波变换对时间序列的周期分析、M-K趋势检验对时间序列趋势分析和M-K突变检验对时间序列突变分析的功能。同时，ESMD与Morlet小波变换相比，具有更高的周期分辨率，与M-K趋势检验验相比，可以精细地刻画径流的具体变化过程，与M-K突变检验相比，避免了得到较多虚假的突变年份，为径流时空演化规律的分析开辟了新思路。

（2） 黄河源区多年径流下降趋势不显著，而在20世纪90年代以后下降趋势变得尤为显著，通过对黄河源区降水、气温、蒸散发分析发现，导致径流减少原因主要是因为降水的增加趋势不明显，而气候变暖引起黄河源地区温度升高，蒸散发量显著上升，导致降雨产流能力下降，使得黄河径流出现下降的趋势。同样地，长江源和澜沧江源地区温度升高，蒸散发量有明显的增强，但由于长江源、澜沧江源地区的降水上升趋势比较明显，所以20世纪90年代以后的径流不仅没有下降趋势，反而有明显的上升趋势。定量分析导致三江源地区径流变化因素可知，降水是导致长江源和澜沧江源径流增加的主要原因，贡献率高达90%以上，在黄河源区，降水不是导致径流减少的主要因素，降水对径流减少的贡献率仅占11.17%，其他因素占88.83%。

（3） 三江源径流过程在20世纪90年代以后和2005年以后发生两次明显的突变，第一次突变可能是20世纪90年代以后降水显著减少，导致三江源地区径流普遍减少。第二次发生明显的突变可能是2005年开始实施的三江源生态保护和建设工程，三江源地区植被覆盖率明显提高、湖泊水域和湿地明显增加，河道的径流有显著提高。

（4） 本文在分析三江源径流演变规律时，缺少对2016年以后的径流资料的分析。在探讨三江源径流变化原因时，没有考虑到冻土退化对于径流的影响，这些问题还需要进一步开展研究。

5 结 论

利用ESMD方法从趋势、突变、周期3个方面，全方位地研究径流序列的时空演化规律。并采用M-K检验和Pettitt检验、Morlet小波变化和周期图法等方法对比分析，得到如下结论：

（1） 黄河源年径流呈现逐年不显著的下降趋势，长江源、澜沧江源年径流呈现显著上升趋势。

（2） 黄河源年径流在1990～1994年以及2008年左右发生突变;长江源年径流过程2006年左右发生突变;澜沧江源年径流过程在20世纪90年代以及2008年左右可能发生突变。

（3） 黄河源径流序列存在12～15 a和30 a的周期性变化;长江源逐年径流序列存在13～15 a的周期性变化;澜沧江源多年径流序列存在15～16 a与28 a的周期性变化。

参考文献：

[1]  王金良，李宗军.可用于气候数据分析的ESMD方法[J].气候变化研究快报，2014（3）：1-5.

[2] 赵直，徐晗.极点对称模态分解下中国新疆温度变化趋势的区域特征[J].地理研究，2014，33（12）：2358-2366.

[3] 段志鹏，李继清.基于极点对称模态分解的北京市降水特征分析[J].中国农村水利水电，2018（3）：17-21.

[4] LIN Q，WU Z，SINGH V P，et al.Correlation between hydrological drought，climatic factors，reservoir operation，and vegetation cover in the Xijiang Basin，South China[J].Journal of Hydrology，2017，549：512-524.

[5] 姜鑠，王金良.关于非对称数据的ESMD改进算法[J].应用数学进展，2018，7（12）：1500-1505.

[6] 易湘生，尹衍雨，李国胜，等.青海三江源地区近50年来的气温变化[J].地理学报，2011，66（11）：1451-1465.

[7] 张士锋，华东，孟秀敬，等.三江源气候变化及其对径流的驱动分析[J].地理学报，2011，66（1）：13-24.

[8] 商放泽，王可昳，黄跃飞，等.基于BUDYKO假设的三江源径流变化特性及量化分离[J].同济大学学报（自然科学版），2020，48（2）：305-316.

[9] 刘希胜，李其江，段水强，等.黄河源径流演变特征及其对降水的响应[J].中国沙漠，2016，36（6）：1721-1730.

[10]  朱海涛.长江源区长序列径流变化规律及其与气象要素的关系分析[J].中国农学通报，2019，35（22）：123-129.

[11] 王珂，蒲焘，史晓宜，等.澜沧江源区气温与降水对径流变化的影响[J].气候变化研究进展，2020，16（3）：306-315.

[12] 王金良，李宗军.极点对称模态分解方法：数据分析与科学探索的新途径[M].北京：高等教育出版社，2015.

[13] 郭爱军，畅建霞，王义民，等.近50年泾河流域降雨-径流关系变化及驱动因素定量分析[J].农业工程学报，2015，31（14）：165-171.

[14] 姜瑶，徐宗学，王静.基于年径流序列的五种趋势检测方法性能对比[J].水利学报，2020，51（7）：845-857.

[15] 张少文，丁晶，廖杰，等.基于小波的黄河上游天然年径流变化特性分析[J].四川大学学报（工程科学版），2004（3）：32-37.

[16] 强安丰，魏加华，解宏伟.青海三江源地区气温与降水变化趋势分析[J].水电能源科学，2018，36（2）：10-14.

[17] 何秋乐，匡星星，梁四海，等.1966～2015年长江源冰川融水变化及其对径流的影响：以冬克玛底河流域为例[J].人民长江，2020，51（2）：77-85，130.

[18] 邵全琴，曹巍，樊江文，等.三江源生态保护和建设一期工程生态成效评估[J].地理学报，2016，71（1）：320.

（编辑：江 文）

引用本文：

李凱，任卓钰，王永强，等.

基于极点对称模态分解的三江源径流时空演变规律

[J].人民长江，2021，52（8）：84-91.

Spatial-temporal evolution of runoff based on extreme-point symmetric mode decomposition

in Three-River Headwaters Region

LI Kai1，REN Zhuoyu2，WANG Yongqiang1，XU Jijun1

（ 1.Water Resources Department，Changjiang River Scientific Research Institute，Wuhan 430010，China; 2.Department of Hydraulics，Changjiang River Scientific Research Institute，Wuhan 430010，China，Wuhan，Hubei 430010 ）

Abstract：

Climate change and the strengthening of human activities have influenced the earths water cycle system.Therefore，it is of great significance to re-understand and master the evolution of water cycle in the basin.We use the extreme-point symmetric mode decomposition （ESMD） method to analyze the runoff variation in the Three-River Headwaters Region from aspects of trend，sudden change，and cycle，and compare the ESMD with existing detection methods.The results showed that：①  in the source region of the Yellow River，the increasing trend of precipitation was not obvious and the evapotranspiration increased as the temperature increased，which resulted in an insignificant downward trend of annual runoff.The annual runoff in the source region of the Changjiang River and the Lancang River increased significantly due to the increase of precipitation and glacial snow melting.② The 61 years runoff process of the Yellow River source region experienced significantly abrupt changes between 1990～1994 and around 2008.The runoff process of the Changjiang River source region experienced significantly abrupt change around 2006.The runoff process of the Lancang River source region had an abrupt changes in the 1990s and around 2008.The 61-years runoff process of the Yellow River source region had cyclic variation of 12～15a and 30a.The 60-years runoff process of the Changjiang River source region had 13～15 a cyclic variation.The 49-years runoff process of the Lancang River source region had cyclical changes of 15～16 a and 28a.And these cyclical changes dominate the characteristics of the annual runoff variation of the Three-River Headwaters Region.

Key words：

runoff variation;cycle analysis;climate change;ESMD;Changjiang River;Yellow River;Lancang River