基于EMD的岷江上中游流域流量特性分析

王俊鸿，覃光华,2，童 旭

(1.四川大学水利水电学院，成都 610041；2.水力学与山区河流开发保护国家重点实验室，成都 610041)

山区和丘陵地带河流的径流变化过程往往较平原地区更为复杂。受人为或自然等众多因素的影响，这一地区河流的径流变化越来越显著。研究径流在长时间序列中的变化规律，有利于进行水文规律的预测预报，同时指导区域河流水量调度。

国内外研究长时间序列尺度的径流变化规律常用的主要有小波形分析、多维混沌理论、BP激励函数神经网络模型和经验模态分解(EMD)等方法[1]。小波形分析[2]是运用最为广泛的一种数学-水文耦合方法，理论上是一种基于时频局部化的信号，进入新世纪以来在流量、水位、输沙率等多个水文要素时间序列研究中均取得了一定的研究成果。而具有优良自适应性的BP人工神经网络方法[3]，被广泛应用于更复杂的水文变化特征中，例如水文与气象结合的二重尺度交互变化规律[3]。多维混沌理论方法——基于杂乱无序的随时间变化的大量水文数据，分解出非线性特征并对其进行分析，从而进一步创建序列[4]。

经验模态分解方法(EMD)虽然没有以上三种方法使用广泛，但也常被用来分析非线性和非平稳时间序列[5]，它能够把复杂的时间序列组拆分成若干个本征模函数(Intrinsic Mode Function)[6]，即关于序列自身形态模式的函数。这种分解方式具有很强的自适应性，对数据本身特征进行针对性的分解，这在一定程度上提高了分解效率。岳相臣[7]在传统EMD方法的基础上利用端点效应和模态混叠提出新的改进方法，攻破了信号因接受频率而受限的问题；于本成[8]在进行熵性能研究时用到了改进了的EMD算法。另一方面，EMD方法用于水文预测时，其精度比单一神经网络更高[9,10]，因而EMD具有一定的推广价值，目前基于EMD的优化算法已经在其他各学科得到有效应用[11]，在水文领域也有涉及[12]。

岷江上中游，作为典型的山区丘陵河流，基于EMD方法的规律研究较少。本文运用EMD方法对岷江上中游流域6个代表水文站观测到的1956～2016年平均流量序列分别进行多时间尺度分析，以此获得径流序列的时间变化规律，并对其进行特征分析。

1 研究方法

EMD的基本思想是将一个频率不规则的波化为若干个单一频率的波(IMF)与残波的和的形式，本质上是对数据序列或信号的平稳化处理。在分解过程中，每个IMF均需符合以下2个性质：一是信号在时间序列上任意一个特定长度区间段的最大值或最小值的数目和过零点数目相等，或满足指定的某一标准(最大不超过1)，这个条件是为适应数据在频率区间上的密集性需求；二是由特定长度区间段的若干最大值构成的上围线和由若干最小值构成的下围线的中值为零。其目的是为使序列的波形趋于对称，从而方便对瞬时频率进行操作。

根据上述IMF的两个性质，对原始时间序列进行分解时分别采用由最大值产生的局部上包络线和最小值产生的局部下包络线确定若干个极值。当极值确定后，利用三次样条曲线将所有的局部最大值和局部最小值分别拟合成上、下两条包络线，并且需要保证所有的时间序列数值都介于这两条线之间，取m1(t)为这两条包络线的中位值。

初始目标值通过x(t)与m1(t)相减得到：

h1(t)=x(t)-m1(t)

(1)

若一系列的h1(t)不平稳，则视为不满足作为IMF分量的条件，因此要对h1(t)重复以上过程p次进行筛分，有：

hp(t)=hp-1(t)-mp(t)

(2)

直到hp(t)符合IMF分量的条件为止，此时将hp(t)作为确定的第一个IMF分量，记为imf1(t)。为简化筛选程序，可以利用限制标准偏差STDV的值来判断每次筛分结果是否为IMF分量，其定义为：

(3)

式中：T表示待研究的时间序列总长度。

霍译：The words were as oil upon fire.The old lady blazed.“Who gave orders for those coffins…”

这样得到的imf1(t)应是信号的特征时间尺度最小的高频IMF分量，然后进行下一步骤：

x(t)-imf1(t)=r1

(4)

因为剩余序列r1中也许会含有比imf1(t)的特征时间尺度还大的IMF分量，所以又要把 当成新的序列进行上述筛选过程，以求得第二个IMF分量imf2(t)。

这便是EMD分解核心步骤——筛分过程，其目的是为消除畸形波并平滑细微的弯曲，保存大的振幅，简化波形；同时使序列曲线变成围绕零均值线的局部极大值和极小值基本对称的波形。

为使得筛分过程不至于持续进行下去，我们还有这样的标准[13]，在EMD过程中当满足0.2

图1 EMD模型分解步骤Fig.1 EMD model decomposition steps

各个IMF分量分别表示不同的特征时间尺度，用连续两个极值之间的时间间隔作为本征模态的时间尺度，不仅能清晰地分辨出不同尺度的震荡分量，还可用于数值全为正或负的没有穿越零点的序列。根据以上筛分过程，直到imf1(t)或rn(t)的值小于预先设计的定值，或者rn(t)不能进一步分解为IMF分量，以单调函数的形式截止。这里的rn(t)表示原时间序列的趋势，记为趋势项RES[14]。最后，原时间序列可以改写为：

(5)

2 研究区概况

岷江是长江上游的重要支流，在乐山城区以上为岷江的上中游，以大渡河为正源，其流域面积为126 478 km2，河长为1 129 km，落差达到2 795 m，这一范围大致位于龙门山、岷山等众多盆地边缘山脉与成都平原之间组成藏区高原东缘地形陡变带内，平均海拔2 640 m，河道不仅比降大，且两岸以砂壤土质为主。又因其位于我国沿海季风湿润区向高原高寒区的过渡地带，气候条件错综复杂，流域年平均降水500～950 mm，多年平均气温12～20 ℃。岷江上中游出口断面处径流较为充沛，据有关水文资料多年平均流量为1 000 m3/s，丰水期流量占全年流量的75%，而枯水期流量仅占全年流量的25%，丰枯季界限分明。并且，岷江上中游流域因其水系发达，从古至今一直以来也是成都平原的水源补给河流。

3 实例分析

青衣江和大渡河是岷江中游以上重要的一级支流，并且本次选择的这6个站点有上游的大金站、都江堰站，中游的泸定站、彭山站，下游的夹江站和沙湾站，既有山区站点也有丘陵站点，同一河道上最近的两站点间隔在100 km以上，覆盖面较广，同时这些站点都有超过60 a的长时间序列资料，因此都具有较好的代表性，通过这些代表站点的水文资料能较好地反映研究流域的流量特征。各代表站在研究流域上的分布见图2。

图2 各代表站在研究流域上的分布Fig.2 The distribution of representative gaging stations on the study basin

针对上述站点1956-2016 a实测径流资料，运用EMD方法对年平均流量进行分析[15]。经过EMD拆分过程，年均流量将得到若干个IMF分量(分别记作imfi)及一个趋势项RES。见图3，图3(a)～图3(f)中signal表示原始时间序列，imf1/2/3和res.表示原始时间序列经过EMD分解得到的3个ifm分量和1个趋势项[16]，各自的平均周期如表1所示。

表1 各站IMF平均周期Tab.1 Average period of IMF in each station

从表1可以看出，较高频率波动的 所反映的周期基本相同[17]，各站点imf1显示周期集中在2.5～3.5 a，imf2显示周期集中在5.0～7.0 a。而低频波动的imf3显示的周期差异则较大，最小的都江堰水文站周期为10.0 a，最大的泸定水文站周期为17.4 a[18]。同时，岷江上中游各站点年均流量的EMD趋势均反映出近乎直线型的下降趋势[19]，其中大金站的年均流量在2007年开始停止下降，此后近十年趋于平稳。这是由于大金站地处大渡河上游高原山区，在全球气候变化及地区地理气候、和人类环保意识加强的多重影响下，该地区的流量在经历了前几十年的削减后，在近十年得到有效缓解[20]，而这一缓解现象在未来将在岷江流域的其他站点重现。

岷江上中游流域径流量的下降趋势通过乐山市区的五通桥水文站1956-2016年观测得到的数据也能够发现。五通桥水文站地处岷江中游干流上，是岷江上中游流域的出口断面，因而该水文站近半个世纪的年均流量观测数据能较好地反映岷江上中游流域径流量的发展趋势。将所研究的6个水文站和五通桥水文站1956-2016年观测计算得到的年均流量用折线图表示，见图3，图3中虚线为趋势线，趋势线斜率均为负值。

图3 岷江上中游典型水文站及五通桥水文站1956-2016年均流量Fig.3 The average flow of hydrological station of the upper and middle reaches of the Minjiang River and Wutongqiao Hydrological Station from 1956 to 2016

imfi(i=1,2,3)所对应的各站点标准偏差见表2。

表2 imfi所对应各站点的标准偏差Tab.2 Standard deviation of each station corresponding to imf1/2/3

通过计算imfi(i=1,2,3)各站点的STDV，分别为0.38、0.69、2.45，这反映了6个站点在imf1，imf2，imf33种情况下均存在离散情况，并且随着频率波动的减小，其离散程度相应增大。这是因为大渡河、青衣江、岷江三条河流的流域下垫面存在差异，大渡河流域主要以山地为主，青衣江和岷江主要以丘陵为主，这决定了河道的汇流方式和流量的年际分布，所以在高频率波动的imf1，imf2中地势差异起伏大的站点STDV普遍高于地势平缓的站点；而低频率波动的imf3中由于受气象和下垫面等因素的耦合作用[21]，这种差异被放大，同时受人为干预或高原复杂地下水演变的影响，都江堰站和大金站显著偏小，导致imf3的离散程度最大。为试图通过imfi之间的STDV规律来预测岷江上中游流域流量的趋势(图4)，可以发现STDV(imf2)约为STDV(imf1)的2倍，STDV(imf3)约为STDV(imf1)与STDV(imf2)和的2倍。

图4 岷江上中游代表水文站年均流量EMD结果Fig.4 EMD results of annual average flow of representative hydrological stations on the upper and middle reaches of Minjiang River

4 结 语

采用经验模态分解(EMD)对岷江上中游流域6个代表水文站的1956-2016年均流量进行分析，并通过乐山市区五通桥水文站的同时期资料加以复核，得到相应的周期及趋势和统计规律，结论如下。

(1)岷江上中游流域年均流量存在2.5～3.5 a、5.0～7.0 a的短周期和10.0～18.0 a的长周期，且近半个世纪来岷江上中游径流量持续减少。

(2)除水文规律外，还发现随着频率波动的减小(周期增大)，多个站点的年均流量标准偏差相应增大，且低频周期的离散程度在高频周期的基础上并不是逐级增加，而是呈倍数放大，且倍比与高频周期的离散程度密切有关，STDV(imf2)约为STDV(imf1)的2倍，STDV(imf3)约为STDV(imf1)与STDV(imf2)和的2倍。这一规律过去很少有研究涉及，这表征了水文规律的预测不仅可以从单一序列自身规律着手，还可以从一个流域的多个序列为对象，通过他们之间潜在的关系来进行预测，这符合水文实际情况，同时也使预测结果更客观全面。

但由于本研究仅就岷江上中游流域作为研究区，这一地区的地形、气候、植被[22]等变化呈规律性阶梯分布，并且受人为影响明显，因此这一结果在未对更多流域进行验证时仅适用于岷江流域的诸多河流。望今后更多学者在基于标准偏差的多站点水文规律研究中有所涉足。