近70年黄河入海水沙通量演变特征

王俊杰, 拾 兵, 巴彦斌

(中国海洋大学 工程学院, 山东 青岛 266100)

河流入海水沙通量对河口地貌塑造、岸滩侵蚀以及近海岸泥沙输移有着重要作用[1]；同时，在河流水动力的作用下会使得泥沙颗粒携带大量污染物迁移到河口及近海岸，严重影响河口生态环境[2]；再者，河口及近海岸区域存在的悬浮泥沙颗粒会与水体微生物以及其他有机物相互作用形成悬浮物，悬浮物的大量存在会影响水体透明度，进一步影响水体生物的光合作用和初级生产力[3]。因而关于入海水沙通量的研究是陆海相互作用研究的重要科学问题之一。

目前，关于黄河水沙过程的研究已引起了广泛关注，但研究区域多集中在中上游段，并且对水沙过程的特征分析不够全面。如张金萍等[4]仅研究了黄河潼关站水沙过程的趋势性、周期性，并未对其变异性进行讨论，且方法较为单一。欧阳潮波等[5]对黄河河龙区间的水沙长期趋势、年际变异、以及水沙过程对人类活动的响应进行了讨论，对水沙的周期性并未涉及。姚文艺等[6]利用“水文法”和“水保法”对黄河上中游水沙演变的驱动成因和未来趋势进行了研究，认为未来黄河水沙量仍处于下降趋势。另外，其他众多研究成果[7-10]表明黄河入海水沙通量呈显著下降的趋势，且人类活动(水库建设、引水引沙和水土保持措施等)的加剧以及气候变化的影响是水沙通量减少的主要因素。那么，在多种耦合因素的影响下，最近70 a黄河入海水沙呈现怎样的多尺度特征？径流和输沙是否存在相似的演变规律？以及水沙之间是否具有协同影响和共振周期？为了明晰以上问题，本文以黄河利津站1950—2018年水沙资料为基础，应用滑动平均法对水沙趋势性进行讨论，以Mann-Kendall突变检验法、滑动T检验等对水沙序列进行年际变异性分析，同时利用经验模态分解和Morlet小波分析对水沙过程的多尺度周期进行探讨，另外利用交叉小波分析和小波相干性对径流和输沙的共振周期进行研究。旨在为进一步加强对黄河口及其邻近海域的生态治理、三角洲的侵蚀保护等提供理论依据。

1 研究区概况与数据来源

黄河流域处于东亚海陆季风区的北部，是中国的第二大河，以水少沙多而闻名于世，具有“水沙异源”的显著特点[11]。利津站是黄河流域的入海控制站,在利津水文站至入海口之间的泥沙淤积量无法确定时,可将利津站实测悬移质输沙量近似视为黄河入海泥沙通量[12]。本文所探讨的黄河入海水沙通量数据为1950—2018年利津水文站实测资料，其主要来源分为三部分：一是《全国主要河流水文特征统计》，时间跨度为1950—1979年；二是黄河利津水文站，时间跨度为1980—1999年；三是《中国河流泥沙公报》，时间跨度为2000—2018年。资料均通过水文的“三性”审查，研究区利津水文站地理位置如图1所示。

2 研究方法

趋势性、周期性和突变性是水文时间序列的重要特征[13]，本文采用5 a滑动平均法、Mann-Kendall突变检验法、滑动T检验等方法对水沙序列的趋势性和变异点进行分析，并通过小波变换、经验模态分解等方法对水沙通量的周期性进行研究，利用交叉小波分析和小波相干谱对径流量和输沙量的共振周期与相干性进行研究，其中Mann-Kendall突变检验法、滑动T检验、小波分析等方法为水文分析常用方法，因此不再赘述，经验模态分解方法和交叉小波分析方法如下：

图1 研究区及水文站位置

(1) EMD方法。经验模态分解是Huang等人[14]在1998年提出的一种自适应信号的分解方法，它利用信号内部的特征尺度变化作频率与能量的解析，将非平稳、非线性信号分解为若干本征模态函数(intrinsic mode function,IMF)和一个残余分量。对于给定的原始信号X(t)，其分解过程为：

(1)

② 计算h1(t)=X(t)-m1(t)

(2)

③ 判断h1是否是IMF，若不是，则继续重复① 、②过程，直到均值包络线趋于零，此时，记h1为c1，c1即为筛选出的第一个IMF，它表示原始信号的最高频分量。

④ 将c1从原始信号中分离，得到去除一个高频分量的新信号r1(t)，即：

r1(t)=X(t)-c1(t)

(3)

r1(t)也称为残余信号，再对新信号重复以上步骤，直到第n个残余信号变为单调函数，即可终止分解，最终得到一组IMF和一个残差信号rn(t)。原始信号X(t)可以重构为：

(4)

式中：rn(t)为残差项，代表信号的平均走势。

小波相干性(wavelet coherence，WTC)用来反映两个时间序列在时频空间的相干程度，其定义为：

(5)

小波相干性显著性检验采用以红噪声为标准谱的Monte Carlo方法。

3 结果与分析

3.1 水沙通量趋势及变异分析

为分析黄河入海水沙通量的趋势，分别对1950—2018年利津站水沙时间序列做趋势分析和五年滑动平均处理。近70 a黄河入海年径流量呈现波动下降的趋势见图2A。在1950—1970年，径流量平均值较高，处在波动的偏高期，1970年之后整体下降较快；近70 a下降的趋向率为64.3亿m3/10 a，年均径流量为293.7亿m3。输沙量在近70 a也呈现出与径流量相似的趋势特征见图2B。近70 a泥沙下降的趋向率为2.14亿t/10 a，年输沙量平均值为6.62亿t，两者的相关性表明，在α=0.01的置信度水平上呈显著正相关，Pearson相关系数为0.853。

近70 a，在人类活动(水库修建、水土保持措施、干流引水等)和气候变化双重因素的作用下[16]，黄河入海水沙通量发生了深刻变化。1970年之前，黄河流域大型水利枢纽工程较少，河流连通性顺畅，且流域降雨充沛，入海水沙量主要受气候变化影响；70年代以后，黄河干支流陆续修建了刘家峡、龙羊峡、青铜峡、小浪底等水利枢纽，黄河流域的径流逐渐被干支流的水库控制，其中，黄河干流中、上游水库除三门峡外共有7座，总库容高达314.12亿m3[11]，改变了黄河天然的径流状态，河流挟沙能力下降，1986年之后的年平均入海泥沙仅为1986年之前的23.2%。

在分析1950—2018年黄河入海水沙趋势的基础上，通过M-K非参数检验法、滑动T检验、累积距平法对径流和输沙的突变年份进行确定。在入海径流量M-K统计量曲线中，在0.05置信度水平内，UF和UB曲线在1985年出现交点见图3A；入海泥沙量则在1995年，1996年两个年份出现交点见图3B。因此，对应的时间节点可能为水沙时间序列的突变年份。

图2 黄河入海径流量、输沙量变化趋势

图3 黄河入海径流量、输沙量M-K检验

为了进一步阐明近70 a入海输沙量和径流量的突变情况，结合水沙序列的累积距平曲线(图4)，发现1950—2018年水沙时间序列的累积距平曲线大致呈倒“V”型，同时其拐点出现在1985年，说明在1950—1985年，入海水沙量呈现增长的趋势，在1985年之后开始出现下降，整体上经历了由“丰”到“枯”的变化；再进一步结合滑动T检验发现，在步长为3，0.01显著水平下，径流量UF和UB曲线相交的1985年通过滑动T检验；对于输沙量，只有1996年通过了步长为4，显著水平为0.01的滑动T检验，因此，在结合多种方法的基础上，认为入海径流量在1985年发生变异，入海输沙量在1996年发生变异，尽管径流量和输沙量的时间序列具有较强的相关性(相关系数为0.853)，但两者的变异年份仍有一定时差，这与黄河流域“水沙异源”特性以及人类活动干扰等因素有关[12]。

图4 径流和输沙的累积距平曲线

3.2 水沙通量多时间尺度特征

本文采用经验模态分解(Empirical Mode Decomposition，EMD)与Morlet小波分析的方法对1950—2018年黄河入海水沙时间序列进行多时间尺度特征分析。经验模态分解是一种数据驱动算法，具有较好的灵活性和自适应性，能够根据数据固有特征进行分解[17]，图5—6分别为径流量和输沙量的经验模态分解结果，可以看出1950—2018年黄河入海径流量分成了3个本征模态函数(IMF)和1个趋势项(Res)，输沙量被分解成了4个IMF和1个趋势项，分别对径流量和输沙量的IMF和趋势项进行重构，发现与原始序列完全重合。可以认定EMD方法分解结果正确可靠，同时通过对径流量和输沙量的趋势项(Res)的大致走向分析，发现径流量在1950—2000年处于减少的趋势，在2000年左右以后，出现了小幅度上升态势，而输沙量一直处于减少的趋势，这与上文的趋势分析中的结果基本一致。

为了分析径流和输沙的不同时间尺度特征，对径流和输沙经过经验模态分解的各个分量利用FFT求周期的方法进行分析[18]，结果见表1。

图5 径流量EMD分解

图6 输沙量EMD分解

表1 周期识别结果

项目IMF1IMF2IMF3IMF4径流3～59～1119～21—输沙3～45～69～1119～21

由表1可知，径流和输沙在不同时间尺度上的周期基本一致，其中输沙量的IMF1和IMF2的周期尺度相近，这可能由经验模态分解的局部筛选性质制约而造成“模态混叠”现象有关[19]，但不影响对水沙时间序列周期的识别。黄河入海径流量具有3个模态，第一阶模态IMF1频率最高，周期最小，IMF1和IMF2的分别具有3～5 a，9～11 a的尺度周期，反映了径流量的年际尺度振荡，而IMF3的周期为19～21 a，是入海径流量年代际尺度周期；黄河入海泥沙具有4个模态，其中，IMF1是频率最高，周期最小的模态，与径流的IMF1对应，具有3～4 a的尺度周期，IMF2，IMF3周期分别为5～6 a，9～11 a，以上3个模态均为输沙量的年际尺度周期，IMF4周期为19～21 a，是入海输沙量的年代际周期，黄河入海水沙呈现出的多尺度特征，表明黄河流域受到人类活动、气候、下垫面变化等因素的综合作用，具有复杂的演变过程。

在采用经验模态分解对黄河入海水沙进行多尺度特征分析的基础上，亦利用Morlet小波分析进行对比，图7分别为1950—2018年黄河径流量和输沙量的小波系数等值线图，可以反映不同时间段下的周期振荡强弱程度，以及水沙时间序列的丰、枯交替的变化现象，其中，实线为丰水丰沙阶段，虚线为枯水枯沙阶段。黄河入海径流量分别存在3～5 a，9～11 a以及20～22 a的尺度周期，与经验模态分解的结果较为吻合，并且入海径流量在年代际尺度周期上大致经历了“丰—枯—丰—枯—丰”5个阶段，未来一段时间内仍为丰水阶段，且20～22 a的尺度周期具有全局性，存在径流的整个阶段，3～5 a，9～11 a的年际尺度周期在1985年之前振荡较强，在1985年之后逐渐减弱，这与上文M-K检验、滑动T检验认为径流在1985年发生变异的结论具有一致性；同时，通过分析黄河入海泥沙量的小波系数等值线图发现，1950—2018年入海泥沙量大致存在3个不同的时间尺度周期，分别为2～5 a，9～11 a，22～24 a，输沙量小波分析的结论与经验模态分解的结果具有一致性，输沙量在22～24 a的年代际尺度周期上大致经历了“丰—枯—丰—枯”4个阶段，并且未来几年仍为枯沙阶段，输沙量的年代际周期在1970年之后振荡逐渐减弱，9～11 a尺度周期在整个时间段内较为稳定，具有全局性。

图7 径流量、输沙量小波系数等值线

综合经验模态分解与小波分析，发现黄河入海水沙的长期演变过程具有显著的年际以及年代际尺度周期变化特征，同时，入海径流量和输沙量的趋势性以及周期尺度特征基本吻合。并且，在1985年以后，黄河入海水沙的各种尺度周期振荡强度逐渐减弱，这与黄河流域大规模人类活动有关，黄河入海水沙进入了新的阶段。

3.3 水沙通量的共振周期

对1950—2018年黄河入海水沙序列进行交叉小波分析，图8径流量和输沙量的交叉小波谱和小波相干谱；其中，箭头方向表示两者之间的相位关系，由左向右表示为同相位关系，由右向左表示反相位，黑色粗实线表示两者之间达到了95%的红噪声检验，黑色的细实线为小波影响锥线(COI)，该曲线以外区域由于受到边缘效应而不予考虑[20]。由径流和输沙的交叉小波谱可以发现，径流和输沙共同的高能量周期为1～6 a，6～15 a,并且在1960—1972年1～6 a的周期尺度通过了95%显著性检验，且二者呈正相位关系；由径流和输沙的小波相干谱可知，黄河入海径流和输沙在整个时频空间低能量区内存在显著的共振周期，且低能量区的显著相关性远大于高能量区，其中，8～16 a的时间尺度在整个时间段内均通过95%的显著性检验(两端时间受边界影响)，且二者以正相位为主，说明黄河入海径流和输沙演变特征具有极强的一致性。

图8 径流和输沙的交叉小波普与小波相干谱

4 结 论

(1) 近70 a黄河入海水沙呈显著下降趋势，黄河入海径流量在1985年发生变异，输沙量在1996年发生变异，且在1970年之后，黄河入海水沙通量的各种尺度信号出现逐渐减弱的特征。

(2) 1950—2018年黄河入海水沙具有显著的周期变化规律，其中，径流量具有3～5 a，9～11 a，19～21 a年的显著周期，输沙量具有3～4 a，5～6 a，9～11 a，19～21 a的显著周期变化规律。

(3) 交叉小波分析表明1950—2018年黄河入海水沙在1～6 a的时间尺度上具有共同的振荡周期，显著性区域集中在1960—1972年，表明黄河入海水沙之间的相关性在该时频空间下最为显著。小波相干谱分析表明，输沙和径流以正相位关系为主，且黄河入海水沙之间的低能量区的显著性远大于高能量区，反映出黄河入海径流和输沙的演变特征具有一致性。