基于平移相关法的黄河口溯源冲淤范围分析

王党伟,吉祖稳,邓安军

(中国水利水电科学研究院流域水循环模拟与调控国家重点实验室,北京 100048)

随着气候变化和人类活动对地表过程影响的显著增强,世界上大部分河流的河口形貌都在发生着显著改变[1- 2],入海水沙变化对三角洲演变的影响受到广泛关注[3- 6]。上游水沙变化是河口三角洲及流路演变的重要影响因素,与此同时,入海流路变化对上游河道的水沙输移和形态调整的作用也不可忽视,这一点在入海流路变化剧烈的黄河口表现得尤其明显。黄河是世界上著名的多沙河流,由于大量泥沙在河道及河口区域沉积,造成黄河河道快速抬升,黄河三角洲不断增长,黄河入海流路也随之延伸,多种因素作用导致黄河河口流路多次改道,自1855年黄河于河南铜瓦厢决口夺大清河入渤海至今,河口流路共发生了10次改道[7]。河口流路演变导致黄河下游侵蚀基准面发生变化,河口延伸及流路改道是科学制定黄河下游河道治理与河口治理规划的关键问题之一,关系到治黄战略及策略的确定[8]。黄河下游河道演变既受到上游来水来沙过程的影响,一般称之为沿程冲淤,也受到河口淤积延伸、河口流路改道后造成的侵蚀基准面变化的影响,一般称之为溯源冲淤。由于黄河水沙过程时空分布很不均匀,河口流路淤积延伸,同时流路改道造成流路缩短,2种过程相互影响,确定黄河下游沿程冲淤和溯源冲淤的影响范围十分困难。有关黄河口溯源冲淤影响范围的观点主要分为2种:一种观点认为黄河口是黄河下游的侵蚀基准面,河口冲淤演变会造成整个黄河下游河道平行抬升或下降[9- 10];另外一种观点认为,河口演变的溯源冲淤的影响范围是有限的,不会波及到整个黄河下游河道,但是具体范围又存在巨大分歧,主要的结论有3种,分别为泺口以下河段[11- 12],艾山以下河段[13]和利津以下河段[14- 15]。河口演变对黄河下游河道溯源影响范围的认识仍然存在争议,一方面是由于黄河下游河道演变非常复杂,另一方面由于缺乏具有明确物理意义的研究方法。因此,需要进一步从机理上研究黄河河口延伸引起的溯源淤积的影响范围、作用大小和响应时间,以揭示黄河下游与河口河道演变的相互作用与反馈机制[16]。

本文基于沿程冲淤和溯源冲淤在时间和空间上的传播特征,研究黄河口演变溯源影响范围的定量化判别方法,探讨黄河河口流路演变对黄河下游河道的影响机理和过程。

1 数据与方法

1.1 研究采用的数据

黄河口流路多次改道,入海流路长度发生了较大幅度的变化,黄河下游河道的侵蚀基准面随之改变。新中国成立以来河口发生了3次大的流路改道,如表1所示[7]。1953年以来的入海流路先后为神仙沟、钓口河以及清水沟,流路改道点基本都位于渔洼下游附近,如图1所示。清水沟流路是当前行河的河道,仍然较为稳定。由于历史原因,1950年之前黄河下游水沙资料缺乏完整性,因此,研究的起始时段取为1950年。

表1 黄河入海流路改道情况

图1 黄河下游及入海流路示意

图2 黄河下游年内最大日均流量

黄河下游各站在3 000 m3/s流量下的水位变化与相应河段累计冲淤量变化趋势基本一致[7],可以采用各水文站3 000 m3/s流量下的水位变化来反映其所在河道的冲淤变化。由于1990年之后黄河下游汛期流量减小,导致部分断面多年没有出现3 000 m3/s以上的流量(如图2所示)。为保证数据的连续性,同时考虑到黄河流路大幅变化主要出现在20世纪50年代初到20世纪70年代末,本研究中的水位数据采用花园口、高村、艾山、泺口、利津5个大断面(见图1)上1950—1990年的日均水位,从上到下相邻断面之间的距离依次为174 km、190 km、102 km和168 km。

尾闾河道延伸或缩短与黄河口河道的泥沙淤积或冲刷密切关联[12],因此可以采用尾闾长度作为黄河口冲淤变化的指征,尾闾河道长度来源于文献[10]中黄河山东水文水资源局的实测数据。

1.2 研究方法

沿程冲淤和溯源冲淤在影响传播的途径上存在明显差别,沿程冲淤是从上游往下游发展,而溯源冲淤是从下游往上游发展。沿程冲淤发展相对较快,溯源冲淤发展相对较为缓慢,其影响是渐变式的,两者在时空上都存在一个发展过程,可以将沿程冲淤和溯源冲淤的发展过程分别看作是冲淤扰动波向河道下游和上游的传播过程。洪水波是天然河道中最常见的扰动波之一,以洪水波在河道中的传播过程为例来说明扰动波的传播特征和本文方法的基本原理。河道上、下游典型的洪水过程如图3(a)所示。

由图3(a)可见,洪水波从上游往下游传播过程中波形变化不大,洪峰流量略有减少。下游洪水的发生过程与上游洪水过程之间紧密相关,理论上两者之间应该具有比较高的相关度。采用Pearson相关分析法量化2场洪水之间的关联度,其相关系数计算方法为

(1)

采用式(1)计算得到2场洪水之间Pearson相关系数为0.25,结果表明2个洪水波之间相关度较低,分析结果与实际的物理过程之间存在较大偏差,直接采用相关分析不能准确反映2场洪水之间的关联性。通过观察2个洪水波之间的形态可以发现,上游断面波峰的出现时间与下游波峰的出现时间相差5 d,如果把下游洪水过程向前平移5 d,则得到如图3(b)所示的图形,再采用相关分析法计算得到2场洪水之间的相关系数为0.95。由此可见,造成计算结果出现偏差的原因在于,在相关分析中没有考虑物理过程之间的时间差,因此导致分析结果不能反映真实情况,而通过对数据平移后再做相关分析不仅可以更准确地反映出波形之间的关联性,还可以反向推理出扰动波传播的方向和时间差。

图3 典型洪水传播过程

河道冲淤的扰动波在空间上的发展过程具有双向性,河道上、下游断面冲淤过程的时间差是反映河道冲淤演变向上下游传播的重要特征。根据扰动传播的规律和溯源冲淤发展过程的特征,本文在相关分析的基础上,将其中1组数据平移不同时间长度后计算2个数据序列的相关系数,根据相关系数的变化规律不仅可以看出不同断面上冲淤主要是受上游影响还是受下游影响,还可以得到这种影响的传播规律。这种方法的关键在于通过对数据进行平移然后计算和分析相关系数分布,因此称之为平移相关分析法。以花园口和高村的1950—1990年水位数据为例,其具体计算方法如下:

记花园口水位变化数据为x=(X50,X51,…,X90),高村水位变化数据为y=(Y50,Y51,…,Y90),这2组数据的Pearson相关系数为

(2)

按照平移相关分析法的基本概念,2组数据的平移相关系数分别为

(3)

式中:R0为花园口与高村同一年水位变化的关联度；Ri为花园口当年水位变化与高村后i年水位变化的相关度,i=1,2,…,40；R-i为高村当年水位与花园口后i年水位变化的相关度,i=1,2,…,40；其他相关系数以此类推。角标正值代表y序列受x序列的影响,即上游影响下游(沿程冲淤影响)；角标负值代表x序列受y序列的影响,即下游影响上游(溯源冲淤影响)。通过分析平移相关系数分布来确定黄河口流路演变对黄河下游河道冲淤的影响过程和影响范围。

2 结果及分析

2.1 水位变化过程分析

利津水文站是最靠近黄河入海口的水文大断面,该站水位极值随河口流路淤积延伸逐渐上升,与试验中观测到的受河口溯源冲淤影响的水位过程类似[17],尤其是在河口流路发生改道等剧烈变化时,利津水文站的水位都会有相应的变化,两者之间存在着滞后响应关系[18],如图4所示,因此,可以采用利津水文站的水位变化代表河口流路演变和侵蚀基准面的变化。以各断面1950年3 000 m3/s水位为基准,以每年相应流量下的水位与1950年水位的差值来反映河道水位变化过程,结果如图5所示。

由图5可见,黄河下游各站的水位变化趋势总体上是一致的,各断面上水位变化过程的波形基本相同,都呈现出震荡上升的趋势,符合扰动波的传播特征。1950—1964年所有断面上的水位都是先上升后下降的趋势,利津水位从1953年河口改道后就开始显著下降,泺口、艾山水位从1954年开始明显下降,而花园口和高村在1960年之前水位都呈上升趋势。随着新流路逐渐淤积延伸,利津、泺口、艾山水位分别在1955年、1956年和1957年开始回升。黄河下游各站水位都在1960年达到阶段性高点,花园、高村水位相对1950年上涨约1 m,其余三站水位上升幅度不足0.5 m。1960年三门峡开始运行,进入下游水沙量大幅减少,与此同时河口流路出汊导致流路长度缩短(见图4),从花园口到利津的水位都开始显著下降,直到1964年各站水位都达到了阶段性低点。1964—1973年三门峡水库大量排沙导致进入下游沙量大幅增加,1964年河口流路改道后新流路又开始淤积延伸,因此1964年之后各站水位又开始快速上升,1975年各断面水位都比1950年上涨了2 m左右。1975—1977年黄河下游连续发生了大洪水(见图2),高村以上河段主糟大量冲刷,加上1975年河口流路改道河长大幅减小,在此期间黄河下游各站水位都有较大幅度下降,花园口和利津站水位下降幅度最大。1978年之后黄河下游各站水位快速上升,直到1980年之后各站水位又开始下降,上游水位变化与1980—1982年汛期大洪水有关[19],而河口水位下降则与河口流路水面比降增加关系紧密[20]。1985年之后黄河下游各断面上水位经历了一次上涨过程,至1990年花园口、高村、艾山、泺口和利津的水位相对于1950年分别上升了1.78 m、2.41 m、2.43 m、2.44 m和1.90 m。从水位变化可以看出黄河下游河道整体是淤积抬升的。

图4 利津水位与黄河口入海流路的变化过程

图5 黄河下游不同断面上水位随时间变化过程(Q=3 000 m3/s)

虽然各断面上水位变化趋势基本一致,但是相同时间点上各断面间水位变化幅度却会存在较大差异。如1958年从花园口到利津的水位变幅分别为0.99 m、1.38 m、0.44 m、-0.01 m、-0.33 m,靠近河口河段冲刷,而其上游则是淤积的,冲淤变化相差较大；而1972年从花园口到利津的水位变幅分别为1.57 m、1.76 m、1.59 m、1.71 m、1.54 m,所有断面上水位基本是平行抬升的。这种复杂的水位变化现象也是造成黄河口演变对上游溯源影响研究争论较大的原因之一。

2.2 冲淤影响分析

以黄河下游花园口、高村、艾山、泺口和利津5个大断面上的水位数据为基础,采用平移相关分析法计算得到不同断面上水位的相关系数,如图6所示。图中每一个时段代表1 a。横坐标0处的相关系数即为R0,横坐标1处的相关系数为R1,其他以此类推。

图6 黄河下游不同断面间水位变化的平移相关系数

从图6可以看出,花园口和高村之间的平移相关系数最大值为R1=0.934,说明花园口和高村的水位之间存在明显的相关性,且花园口的水位变化先于高村,领先时间约为1 a。以R1为中心,随着平移时段数的增加,花园口和高村之间水位的相关性开始显著降低,平移相关系数的分布类似于开口向下的抛物线形态。高村与艾山、艾山与泺口的水位之间的平移相关系数也类似于开口向下的抛物线形态,相关系数的最大值分别为R0=0.915、R0=0.975,说明高村与艾山、艾山与泺口的水位之间是高度相关的,但在年尺度条件下难以判定其变化的先后性。泺口与利津的水位之间的平移相关系数的分布与其他几个断面之间有比较明显的差别,不再是类似于抛物线形态,而是呈现左侧平移相关系数显著高于右侧的分布形态,R-8到R0值均在0.85左右,从R-9往左侧平移相关系数才开始显著减小,从R0往右侧平移相关系数快速下降,最大值为R-1=0.874,说明泺口水位变化滞后于利津,即泺口在河口冲淤演变的溯源影响范围内,且泺口当年的水位变化与利津前8 a的水位变化之间的关系都是比较紧密的,这也符合溯源冲淤影响具有滞后性的特点[18]。

图7所示为利津前8 a水位变化平均值、利津当年水位变化分别与泺口当年水位变化的关系。由图可见,利津前8 a水位平均值与泺口水位在拟合线两侧分布较为紧密,证明这2种物理量之间的密切程度显然高于利津与泺口当年水位之间的关联程度。相关分析结果也显示,利津前2 a、4 a、6 a、8 a和10 a水位变化平均值与泺口当年水位变化值的相关系数分别为0.887、0.901、0.920、0.952和0.926,利津前8 a水位变化平均值与泺口当年水位变化值的相关系数最高,且明显大于利津与泺口当年水位变化的相关系数0.857。由此可见,泺口水位变化主要受包括利津当年在内的前8 a水位变化的累积影响。

根据平移相关分析的原理,图6横坐标0左侧的平移相关系数代表溯源冲淤的影响程度,横坐标0右侧的平移相关系数代表沿程冲淤的影响程度,左右两侧的平移相关系数平均值的对比可以看出溯源冲淤和沿程冲淤的影响程度。左右两侧的平移相关系数均值分别记为RL和RR,结果如图8所示。

图7 利津与泺口水位变化关系

图8 不同河段沿程冲淤与溯源冲淤影响的对比关系

从图8中可以看到,反映沿程冲淤影响的平移相关系数值从上游往下游先增后减,说明沿程冲淤的影响从上游往下游总体是减小的；反映溯源冲淤影响的平移相关系数值从下游往上游逐渐减小,说明沿程冲淤的影响从下游往上游逐渐减弱。平移相关系数的分布符合黄河下游河道的冲淤特征。花园口至艾山河段RR>RL,泺口至利津河段RR>RL,两条曲线的交汇点位于艾山和泺口之间,距离花园口断面约380 km,距离2019年的黄河入海口约350 km。在本文的研究时段内,该点以上以沿程冲淤为主,以下则受溯源冲淤影响较大。

为了避免两侧取平均时人为因素对结果可靠性的影响,分别计算了1 a到10 a平均的RR和RL值,得到了RR和RL交叉点的位置,结果如表2所示。由表可见,RR和RL交叉点的位置在距离花园口359～380 km范围内,取平均所用的年数对计算结果影响不大。

由此可见,黄河口流路演变对黄河下游河道的影响范围主要在艾山以下,大范围上与已有研究结果基本一致[13],说明了本文研究方法和研究成果的合理性。同时从平移相关系数分布对溯源范围进一步细化,得出相对更为确切的溯源影响范围。

表2 取平均年数对RR和RL交叉位置的影响

3 结 论

以1950—1990年黄河下游河道的水位数据为基础,采用平移相关分析法确定了黄河口流路演变对黄河下游河道的溯源影响过程和溯源范围,得到了以下主要结论:

(1) 受入口水沙和尾闾河道演变共同作用,黄河下游不同断面上水位变化过程存在一定的差异性,但不同断面上的水位变化过程总体上是相似的,具有扰动波的传播特征。

(2) 花园口、高村、艾山之间冲淤变化时间相差在1 a左右,泺口断面的冲淤变化与利津断面前8 a的冲淤演变都有比较紧密的关联。

(3) 1950—1990年期间黄河口流路演变溯源影响范围主要集中在艾山以下,距离黄河入海口约350 km范围内的河道冲淤演变主要受河口流路演变的影响,该河段以上河道主要受沿程冲淤变化的影响。