三峡调蓄前后径流变化的多尺度分析

任逍迪，李继清，纪昌明

(华北电力大学可再生能源学院，北京 102206)

0 引 言

水文时间序列规律性的探讨一直以来都是水文研究的重点，径流时间序列复杂多变，既存在着趋势、突变等非周期性变化特征，又存在着以年、月或日不同尺度的周期性变化特征。多时间尺度分析(multiple time scales)指水文系统的变化表现出多层次的时间尺度周期结构及局域化特征，即水文系统并不表现为一种恒定的周期，其一段时间以一种周期进行变化，下一时刻又表现为另一种不同的周期，即便是在同一时段内，也会表现出不同时间尺度的周期[1]；对径流进行多时间尺度的研究，可以为我们提供水文径流序列更合理的时间变化特征结果，更准确的揭示径流实际变化的多层次及局域化特征，为水文分析及水文预测等提供重要依据[2]，水文径流的多尺度分析也因此成为了21世纪水文学研究的前沿课题，在李眉眉成功将混沌分析方法引入径流的多尺度分析上后[3]，杨明金、李舒及孙青雪等利用Mann-kendall法相继对黑河、窟野河及青山库区流域径流进行研究并取得了相应多时间尺度研究成果[4-6]，王文圣、刘俊萍、董前进等则利用小波分析方法分别进行了宜昌站、金沙江及黄河流域的径流的多尺度研究[7-11]；径流变化过程具有极高的复杂性和多层次性，传统径流方法并不能够做到基于数据的自适应性分析，且受到较多如傅里叶变换的束缚，近些年提出的极点对称模态分解(ESMD)方法是一种具有极高灵活性和自适应性的数据分析方法，目前被广泛的应用于信息科学、海洋和大气科学、生态学、经济学、地震学和医学等领域数据的分析研究上[12-15]，王飞、段志鹏等分别将该方法应用于黄河SPEI序列及北京市降水序列的多尺度分析中[16,17]。但目前该方法在水文径流分析上的应用仅处于起步阶段，鉴于此本文采用ESMD对宜昌站1950-2016年实测径流进行多尺度分析，并结合传统径流分析方法中的M-K法、小波分析法及希尔伯特黄变换法进行对比，得到三峡调蓄前后径流的变化特征，进而为附近水库的水文预报、洪水调度及相应流域水资源的调配管理提供依据。

1 研究方法

利用极点对称模态分解(ESMD)方法对径流趋势、突变及周期进行多尺度分析；同时利用传统径流分析方法中的Mann-Kendall(M-K)法、小波变换法及希尔伯特-黄变换法分别对趋势、突变及周期规律也进行了分析并与ESMD法进行对比验证。以上方法具体介绍如下。

1.1 极点对称模态分解(ESMD)

极点对称模态分解(ESMD)方法是一种将美国工程院黄锷院士提出的希尔伯特-黄变换方法进行改进后提出的新方法。其采用内部极点对称方式插值，得到序列中内在的不同尺度或周期性振荡和趋势分量，利用“最小二乘法”思想确定最佳筛选次数从而解决了“模态混叠问题”[18]，现被广泛的应用于信息科学、海洋和大气科学、生态学、经济学和医学等领域。该方法具有更高的灵活性及更强的表现力，十分适合非线性非平稳数据的多尺度分析。该方法包括模态分解和时频分析两部分：对一系列原始序列利用极点对称原理进行模态分解，得到多个模态函数及一条代表原始序列趋势项的最佳自适应全局均线(AGM)，所得不同模态独立，代表不同尺度(周期)，基于数据的离散型特征，极点对称模态分解(ESMD)方法的时频分析则采用了“直接插值法”(Direc Interpolating Method)或“旋转生成法”(Rotary Generation Method)，利用此方法不但可以得到各模态的振荡强度与变化快慢程度，还可明确得知总能量的变化。其中模态分解过程为：

(1) 所研究原始序列记为X，找出X的全部极值点(包括极大值点及极小值点)依次记为Ai(i=1,2,3,…,n)。

(2)对相邻极值点用线段连接，线段中点依次记为Fi(i=1,2,3,…,n-1)。

(3)通过一定方式补充左、右边界中点并记为F0，Fn(本文利用线性插值法)。

(4)利用上述步骤所获取的n+1个中点构造p条插值线L1,L2,…,Lp(p≥1)，并计算它们的均值曲线为L*=(L1+L2+…+Lp)/p。

(5)对X-L*序列重复上述步骤，直到|L*|≤ε(ε是预设的允许误差) 或筛选次数达到预先设定的最大值K，由此完成分解得到第一个经验模M1。

(6)对X-M1重复上述步骤，获得M2,M3,… ，直到剩余序列R只剩余一定数量(预设)的极点，便可分别得到经验模M2,M3,… 。

(7)让最大筛选次数K值在限定区间[Kmin，Kmax]内改变，重复上述6个步骤。分别计算不同K值对应方差比率σ/σ0，找出σ/σ0最小时所对应的K0，以K0作为限制条件据此重复前六步得到最终分解结果。

1.2 Mann-Kendall(M-K)法

Mann-Kendall 法(简称M-K法)是一种普遍应用在水文研究中的非参数统计检验方法，利用该方法可分析得出待研究序列趋势项。现已知x1,x2,…,xn为所研究的一系列时间序列，n为序列长度，构造秩序列Sk如下：

(1)

其中：

(2)

定义统计变量：

(3)

式中：

(4)

(5)

当其UFk值大于0时，研究序列变化趋势表现为递增，UFk值小于0时，变化趋势表现为递减；其中所定义的UFk满足标准正态分布，故对序列趋势项的显著性需通过统计量UFk的假设检验，即在给定显著性水平α下，若|UFk|>Uα/2，则表明序列存在较为显著的变化趋势，反之则不显著。趋势的显著性通过其与临界直线的关系来表征，小于临界直线表示变化趋势不显著，反之则明显。

1.3 小波变换法

小波变换法是20世纪后期发展起来的一种数据时频分析方法，目前被广泛的应用在信号处理、模式识别、故障诊断、气象分析及水文研究中。其频率具有极佳的时变性，故其十分有利于多时间尺度周期的分析。利用小波变换方法进行多时间尺度周期分析的具体步骤：

(1)对待研究序列Xi(i=1,2,…,n)进行预处理。

(2)取定恰当的小波基函数，本文选用Morlet复小波函数；其中，Morlet复小波函数是一种高斯包络下的函数，其表现出单频率性及复正弦性，该函数形式表示为：

ψ(t)=eiω0 te-t2/2

(6)

式中：ω0为常数，i表示虚数；该函数在傅里叶变换下可表示为：

由此Morlet小波的周期T与时间尺度α则表现出如下关系：

(7)

(3)进行小波变换；将实数区域内的平方可积函数空间定义为L2(R) ，所研究序列Xi(i=1,2,…,n) 为所定义函数空间L2(R)内的序列，则连续性小波变换为：

式中：ψ(t) 为复共轭函数；Wf(a,b)为小波系数，其隐含待研究序列的时频等特征，可用来分析序列所含不同尺度时频特征。

(4)绘制小波方差图，得到研究序列多尺度周期变化规律。

1.4 希尔伯特黄变换

希尔伯特-黄变换(EMD)法是由美国工程院黄锷院士等提出的一种自适应性数据处理分析的新方法，现被广泛应用于科研、工程等数据分析相关领域[19]。利用希尔伯特黄变换法可分析序列各时间尺度的频率变化，判断各时间尺度下突变发生的时间。希尔伯特-黄变换由经验模态分解和希尔伯特谱分析组成；其中经验模态分解过程如下：

(1)找出原始序列x(t) 的所有局部极值点，利用三次样条函数计算该序列的上包络线u1(t)和下包络线u2(t)的局部均值为：m1(t)=1/2[u1(t)+u2(t)]。

(2)计算h1(t)=x(t)-m1(t)；

(3)将h1(t)作为新的序列重复上述操作过程，直到h1(t)的数学期望变为0。

(4)上述步骤后得到的零均值过程h1(t) 为第一个模态分量c1(t)，其代表原始序列的最高频分量。

(5)将原始序列x(t) 减去第一个模态c1(t)，得到一个新的信号过程序列r1(t)=x(t)-c1(t)。

(6)将r1(t) 作为新的原始数据，重复(1)～(4)步骤，依次可得到c2(t)、c3(t)、c4(t)等，满足迭代条件循环停止，原始信号重构为：

式中：rn(t)为剩余趋势项可代表序列平均趋势。

模态分解完成后，利用希尔伯特谱分析(HHT时频分析)方法对各本征模态进行时频分析得到序列的瞬时频率变化过程。

2 应用实例

2.1 研究区概况及数据来源

宜昌水文站位于三峡工程下游距三峡水利工程约44 km 处，地处长江上游与中游的交界处，被视为干流上中游的咽喉。已知宜昌水文站控制流域面积占到整个流域面积的55%，该站所测径流多年平均径流量为4 321 亿m3，其中汛期平均径流量为3 412 亿m3，占全年径流量的79%左右；三峡水库的蓄水极有可能改变下泄径流水沙的行进规律和时空变化，进而造成径流变化，故三峡工程调蓄前后宜昌站径流变化的多尺度分析亟待进一步深入的研究[20]。本文采用长江流域宜昌水文站1950-2016年实测(年、月、日)径流量资料进行分析。研究区域见图1。

图1 研究区域图Fig.1 Research area

2.2 径流趋势项分析

采用ESMD法分别对宜昌站1950-2016年的多年实测年、月、日平均径流量序列进行模态分解，结果如图2～4所示。年径流序列模态分解最大迭代次数为50，剩余极值点最少个数为6，最佳筛选次数为50次时，分解出了4个模态以及一个趋势余量R；趋势项余量R可表征序列整体趋势变化特征，通过观察ESMD方法分解所得宜昌水文站年平均径流量序列的自适应性趋势余量R(见图2)的变化，可得出过去67年内三峡工程宜昌水文站年平均径流量整体变化趋势在1950-2003年维持平稳，2003年后年平均径流量开始表现出逐年下降的趋势，且趋势较为明显。

图2 年径流序列ESMD分解的模态分量及趋势项Fig.2 Modal components and trend terms of ESMD decomposition of monthly runoff series

同理，月径流最大迭代次数为40，剩余极值点最少个数为10，最佳筛选次数为24次，分解出了6个模态及趋势余量R(见图3)，通过观察分解后所得剩余趋势项余量R得出宜昌水文站月平均实测径流整体为轻微波动的下降趋势。

图3 月径流序列ESMD分解的模态分量及趋势项Fig.3 Modal components and trend terms of ESMD decomposition of annual runoff series

日径流最大迭代次数为40次，最少剩余极值点数为4个时，最佳筛选次数为22次，分解出11个模态分量及一个趋势余量R(见图4)，三峡工程宜昌水文站各日平均径流量的整体变化趋势为小幅度、不显著的下降趋势。

图4 日径流序列ESMD分解的模态分量及趋势项Fig.4 Modal components and trend terms of ESMD decomposition of daily runoff series

利用M-K法对宜昌站年、月、日径流实测资料进行趋势项分析加以验证上述结论，见图5、图6；年径流UFk基本均小于0，且均满足小于临界曲线(-1.96，+1.96)，置信度未通过99%可信度检验，宜昌站年平均径流呈不显著的下降趋势。同样月径流和日径流序列也均表现为不显著的下降趋势，与ESMD法所得结论相一致。

图5 年径流MK法运行结果UF统计量图Fig.5 Annual runoff MK method operation results UF statistics chart

图6 月径流MK法运行结果UF统计量图Fig.6 Monthly runoff MK method operation results UF statistics chart

综上，径流变化趋势不同尺度下在三峡水库蓄水(2003年)前，主要受诸多因素(如气候、流域用水量、流域下垫面条件、水利工程等)影响，具有时间上的随机性，在时间趋势上表现出一定的波动性；而在三峡水库正式调蓄(2003年)后，不同尺度下径流均表现为递减的变化趋势，由于三峡水库的蓄水增加，使水文站观测到的径流量减少。同时不同时间尺度趋势项表现出了一致性，验证了水文研究多时间尺度分析具有嵌套性。

2.3 径流周期分析

ESMD法分解所得的模态分量各自具有独立性及代表性，故每个模态分量各自体现原序列中不同尺度的周期。利用ESMD时频分析周期图法分别估算年、月、日径流各个分解模态的平均周期如表1。

表1 ESMD法不同尺度下径流周期结果表Tab.1 Runoff cycle results table at different scales of ESMD method

同时运用小波变换法对径流进行周期分析得到小波方差图，周期分析结果如表2。小波变换法所得周期结果与ESMD方法所得周期结果相接近，证明所求结果基本合理，但二者结果仍不尽相同，这主要是由于小波分析方法仍存在着如小波基函数的选取、恒定的多分辨率等局限。故与小波变换分析方法相比较，ESMD方法具有更强的自适应性和灵活性，更有助于探索事物的内在规律。

表2 小波变换法周期分析结果汇总表Tab.2 Wavelet transform method periodic analysis result summary table

对比3种尺度下径流的周期结果发现：月径流51.91月及229.43月周期对应年径流5年及17年周期，日径流的368 d及2 879.1 d周期与月径流12.08月及94.47月周期对应相近，周期角度上径流不同时间尺度上呈现出一致性；对于复杂的径流时间序列来说，其同时存在年际、年内等多种尺度的周期变化，并在不同时间尺度上表现出嵌套性。其中周期中表现为3年的年径流周期，其极有可能与ENSO 事件(包括周期为3.5年的厄尔尼若暖事件与周期为5～6年的拉尼娜冷事件)密切相关； 10年的周期与日月交食年现象的平均周期表现为相一致[21]； 17年周期恰好与20世纪30年代以来太阳黑子的活动周期相近似[22]。

2.4 径流突变分析

将分解得到的模态利用直接插值法进行时频分析，绘出其瞬时频率与瞬时振幅的相应时变图如图7，图7中Fi代表各模态频率，Ai代表振幅。图7从模态1的对照关系来看，2003年的点发生了突然的低频、高振幅振荡，表明2003年为突变点，模态3中1984年显示出高频、小振幅振荡也为突变点。同理分析月径流得到1954年7月、1960年1月、1978年5月及1996年9月发生低频高振幅振荡、1972年7月及2003年6月发生高频小振幅振荡为突变点。日径流1954.4.23、1976.10.30、1983.5.23、1972.7.13、1996.9.25、2003.6.21、1963.9.9、1983.5.27等为突变点。

图7 ESMD法年径流数据各模态频率(Fi)与振幅(Ai)时变图Fig.7 ESMD method annual runoff data modal frequency (Fi) and amplitude (Ai) time-varying map

通过希尔伯特-黄变换方法侧面验证三峡工程径流突变性结论，经过运行希尔伯特谱分析相关程序后得到希尔伯特-黄(Hilbert-Huang)谱如图8，根据谱图分布规律，得出宜昌站1950.1.1-2016.12.31期间实测年平均径流量的突变点为1958、1972、1978、1998及2002年；月径流突变点为1954.7、1968.9、1979.7 、1999.2及2003.6；日径流突变点分别为1957.1.21、1973.8.13、1996.9.13、2002.11.16及2008.3.27等，如表3所示。希尔伯特-黄变换法所得突变点中的1973.8.13、1996.9.13分别与ESMD方法的1972.7.13、1996.9.25相接近，证明两种方法所求得结果基本合理，ESMD方法具有更强的适用性。

图8 日径流的希尔伯特-黄谱图Fig.8 Hilbert-yellow spectrum of daily runoff y

突变点年径流/年1958, 1972, 1978, 1998, 2003月径流/(年月)195407, 196809, 1979.7, 199902, 200306日径流/(年月日)19570121, 19730813, 19960913, 20021116, 20080327

综上，不同尺度下径流突变结果如表4所示，其中年月日径流均存在2003年6月左右的突变点，月、日径流均存在1954年7月、1972年7月、1996年9月左右的突变点，年月日径流突变点表现一致性。特别地， 2003年6月的突变点与三峡工程实际正式开始一期蓄水(2003年6月1日)时间吻合，可知三峡工程的蓄水对径流时间序列前后的一致性产生了明显影响，使径流量发生了十分明显的突变，三峡的调蓄作用是引起径流量突变的主要原因。

表4 不同尺度下径流突变结果表Tab.4 Runoff mutation results table at different scales

3 结 语

ESMD方法作为一种新兴的数据分析方法，十分适合非线性非平稳数据的多尺度分析，本文将ESMD法应用于三峡调蓄前后径流的多尺度性变化规律研究上，对三峡工程宜昌水文站1950-2010年径流资料进行分解和分析，进而得出得出以下结论。

(1)三峡工程宜昌水文站径流在不同尺度下均表现为不显著的递减趋势，因此预测未来几年三峡径流量仍有可能继续下降，但并不显著。

(2)三峡径流在不同尺度下的周期表现出相互的包含及嵌套性，其中3年、10年及17年的周期分别与ENSO事件、日月交食年现象及太阳黑子的活动周期可能存在着密切的关系。

(3)突变规律中年月日径流均存在2003年6月左右的突变点，月、日径流均存在1954年7月、1972年7月、1996年9月左右的突变点，3种尺度下的突变时刻相互之间存在一致性，而其中的2003年6月的突变点恰好与三峡工程实际正式开始一期蓄水(2003年6月1日至6月15日)时间段吻合，即三峡工程的蓄水对径流时间序列前后的一致性产生了明显影响。

通过三峡径流的分析研究表明ESMD方法在径流多尺度分析趋势、周期和突变上具有较好的应用效果，其具有更高的灵活性及更强的表现力，十分适合非线性非平稳数据的多尺度分析。但目前该方法在水文特别是径流研究方面的应用还尚在起步阶段，继续开展该方法应用于水文上的研究前景广阔。