黄河入海口亚三角洲演变与下游河道关系探究

尤延锋，王奕童，王远见，许月萍，郭玉雪

（1.浙江大学建筑工程学院水利工程学系，杭州310058；2.黄河水利科学研究院，郑州450003）

0 引 言

随着社会经济和科学技术的发展，人们开始利用淤积形成的河口三角洲进行繁衍生息，创造出了大量的经济价值。黄河是中国的第二长河，是中华民族文明最主要的发源地。历史时期黄河下游入海口发生了多次改道，严重影响了下游生态环境，而在新中国成立以后，通过多次人为改道，花费了巨大投资保持现有河道稳定，现行的清水沟入海口河道已经使用了超过40年。同时，黄河沿程修建了超过3 000 座水库，其中大中型水库170 多座，水库调洪拦沙作用以及黄河中游的水土保持措施对于黄河下游河道的水沙变化情况以及入海口的淤积变化情况产生了巨大的影响。

黄河河道和三角洲的淤积问题一直是研究的热点问题。李健等［1］以数值模拟为研究手段，构建了黄河上游段的连续弯曲河道通航因素分析模型，模型提出的整治方法实际效果明显。陈建国等［2］分析了小浪底水库运行十年之后水库淤积及下游河道的再造床过程及特点。申震洲等［3］总结了黄河砒砂岩区的地域特征特性，并提出了根据空间结构分异特征、植被生境和群落分异特征耦合机理对黄河流域河道侵蚀规律问题等进行研究。蔡蓉蓉等［4］利用SOM-K 神经网络均值聚类耦合的方法对黄河中游潼关水文站进行了水沙组合分类研究。赵连军等［5］考虑了黄河下游河道与河口水沙演进传播的特点，研制开发了计算黄河河道演变与河口演变耦合作用的水沙数学模型。刘慰等［6］选取了黄河下游4 个典型断面的实测资料以及相关的水沙数据，求解了不同时期河道断面沉积速率。杜小康等［7］根据不同影响因素作用下的两级波幅、时均波幅提出了西河口12m 水位的改道标准。彭俊等［8］认为河道淤积冲刷表现的分界点应设置为下游河道含沙量18.6 kg/m3。黄河三角洲的冲淤面积变化有诸多提取方法，遥感影像提取是当前十分有效的方法。相比于传统的海岸线测绘方法，通过遥感技术，将各类物理手段、地学分析和数学方法作为基础对目标进行分析，所获得的数据具有范围大、时效新等特点，是测定海岸线冲淤演变的一种有效手段。Mcfeeters［9］于为了求解水陆边界，于1996年提出了使用归一化差异水体指数（NDWI）结合遥感数据的方法进行影像处理和边界提取；常军［10］通过对1976-2000年期间共20景遥感影像进行处理，对黄河三角洲地区海岸线实现了动态监测；徐涵秋［11］对归一化差异指数（MNDWI）进行了改进，使得研究地区的水体边界计算结果更容易区分阴影和水体。王英珍等［12］通过对比1999-2016年黄河下游游荡段汛后卫星遥感影像与实测淤积断面资料，分析了游荡段的主槽摆动特点。Ji［13］通过遥感影像分析三角洲演变过程，并通过校核潮汐站水位数据，得到了黄河三角洲新河嘴生成的时间序列。河口三角洲的演变与上游的来水来沙条件、河道的演变情况存在复杂的相关关系，人们难以简单地探明他们之间的关系。

本研究将利用遥感影像和下游关键河道断面的演变数据以及入海水沙数据，简要分析黄河近年来水沙变化的情况，其次对入海口的遥感数据进行解析，提取了黄河入海口三角洲面积及河长。最后利用黄河下游关键断面的水沙数据以及淤积数据，构建滞后响应模型对利津站的3 000 m3/s 流量下的水位进行了计算，完整的提出了一种对黄河亚三角洲演变过程特征描述的方法，为深入了解黄河三角洲近年来的淤积变化提供一定的参考。

1 研究区域及数据

黄河下游河口河道的多次变迁导致黄河入海口三角洲在历史上曾经存在过多个入海口，自1976年人工改道清水沟以来，黄河河口区域在人为控制下，基本保持稳定。本文研究区域如图1 所示，下游河长约为768 km，占黄河总河长的14%，下游面积占总流域面积约为3.0%［14］。本文选取花园口和利津两个水文站1950-2010年的日尺度实测水文泥沙数据进行分析。卫星数据选用1985-2015年的Landsat卫星遥感影像。

图1 研究区域Fig.1 Study area

2 分析方法

2.1 遥感影像的提取方法

遥感技术是一种精确、直观，数据获取方便的技术，地理信息系统则在图像处理和空间区域分析中扮演着重要角色，用其研究黄河河口段三角洲的淤积变化具有十分重要的意义。本文选取的遥感数据来源于美国陆地卫星数据Landsat 多时相MSS、TM、ETM+系列数据，其空间分辨率为30 m。本文考虑了研究需求以及相关数据资料特征等，将改进的归一化差异水体指数MNDWI（Modified-NDWI）作为遥感影像的计算工具，通过对水体边界处理后得到的遥感影响进行阈值分割法提取影响，并进行海岸线修正。

式中：Green代表绿光波段（如Landsat TM 影像中2 波段）；MIR为中红外波段（如Landsat TM影像中的5波段）。

考虑Landsat数据获取的影像资料均为地形矫正影像（已经进行了系统辐射、地形处理以及几何矫正等），采用阈值分割算法对影像资料进行预处理由于预处理后图像水陆边界对比度差异显著，水边线可以清晰准确的被分割提取。由于黄河入海口狭窄、潮差小，海岸线被海水影像的时间仅占其裸露时间的极小部分，因此本文中数据处理时没有考虑水位波动。采取以下步骤进行遥感影像的解析以及海岸线的提取：①使用ENVI 5.3 进行遥感影像的波段合成计算，得到MNDWI 影像，确定影像阈值后，对其进行二值化处理，使得影像中水体与陆地的分界更加清晰［图2（a）、图2（b）］；②通过ENVI 5.3 绘制出初步的海岸线边界，将其导入ArcMap 10.2 中，进行边界的提取和细部处理；③进行海岸线边界细部修正，得到最终的海岸线。

图2 遥感影像处理Fig.2 Remote sensing image processing

2.2 河道高程变化确定

主要采用积分方法，通过对各个年份和断面的河道测点数据进行积分计算，河道断面的河床底部平均高程计算方法如下：

2.3 河口滞后响应模型

黄河河口演变趋势有着显著的滞后响应特征［15］。吴保生［16，17］提出的滞后响应模型正是基于河床演变的滞后特征建立的。河床的特征变量y在河床受到外部因素扰动后，特征变量的变化可以用下式表示：

式中：t为时间；β和ye均为常数。

考虑到河床演变特征的滞后性，并不是所有的外部扰动在其发生后都能马上达到相应的稳定状态，在整个滞后响应演变过程中，新的扰动发生将产生新的稳定状态，不断迭代。将新外部扰动发生时刻的特征变量作为新的初始状态，通过递推可得到滞后响应模型的多步递推模式：

式中：y为河道的某一特征变量；ye为该特征变量在河道受到外部扰动之后将要到达的新的平衡状态；n为递推时段；i为时段编号。

郑珊等［15］通过构建黄河河道比降与其水位关系，利用滞后响应模型建立了利津3 000 m3/s 流量下的水位计算方法。将河道比降作为河床的特征量，由（4）式可得。

河道比降的平衡值可以用以下经验关系式进行计算：

式中：W为利津站的年来水量，亿m3；Ws为利津站的年来沙量，亿t；K、a、b均为经验参数。

此外，河道比降还可以利用式（7）计算，即：

式中：Z为利津站的3 000 m3/s流量对应水位；L为利津至入海口出口的河道长度。

基于式（5）、（6）、（7），可以推导得到利津3 000 m3/s 流量水位计算公式：

式 中：ep1=e-β1Δt，ep2=e-β2Δt，并且有F(Ⅰ) =LⅠ，i(1-ep1)。

其中，ZⅠ和ZⅡ分别为利津至改道点河道和改道点河道的3 000 m3/s流量对应水位。

2.4 灰狼优化算法

灰狼优化算法是一种智能优化算法，2014年由Mirjalili［18］提出。在灰狼算法中，严格的等级制度是灰狼算法的核心，狼群被区分为α、β、δ、ω狼，具体的等级制度如图3 所示。其中，α狼是灰狼群体是拥有最强能力和最高优先级的个体，也称之为头狼，在算法中我们将他看作为离最优值最接近的个体；而作为适应度次之的个体，β和δ狼则在捕猎过程中将协助头狼管理群体，制定围猎方案，因此他们也作为α狼的候补者；而狼群剩余的大部分狼群将被认定为是ω狼，主要用于平衡狼群之间的内部关系，并且协助α、β、δ狼进行捕猎行动。

图3 灰狼等级示意图Fig.3 Grey Wolf Optimizer level diagram

当目标的位置确定后，进行个体和猎物距离的计算以及个体位置的更新。

式中：D是狼群与猎物的距离；t是当前的迭代次数；X是灰狼的位置；Xp则是猎物的位置；A和C均为系数，a为控制参数，a值较大时候，全域搜索效率更高，a值较小时，则可以迅速收敛。r1和r2分别为在［0，1］之内随机生成的随机数。

在搜索猎物的过程中，将储存的三个最佳的猎物潜在位置，并根据设置的适应度函数fitness 对潜在位置进行筛选判别后，进行自身更新。如图4 所示，当候选解ω狼在猎物周围时，说明其位置更好，立刻更新优势狼的位置。当灰狼的搜寻次数和搜寻条件达到要求后，求解得到的α狼就是所求的最优解。

图4 灰狼算法位置更新Fig.4 Location update of Grey Wolf Optimizer

本文将利用灰狼优化算法求解滞后响应模型，对利津水位变化特征进行计算分析。

3 计算结果与分析

3.1 花园口段和利津段的水沙变化过程

花园口水文站是黄河下游河道的第一个重要的水文测站，利津水文站是黄河下游的最后一个控制水文站，通过对这两个水文站的水沙分析对黄河下游河道的水沙特征变化情况进行了简要分析。图5所示为黄河下游来水来沙量变化过程。

受到人类活动和气候变化的影响，黄河下游的泥沙和径流量变化较大。从图5（c）看出，花园口和利津站的累计输沙量存在两个转折点，分别为1979年和1999年，现在有观点认为分别是由于1979年实施小流域综合管理措施和1999年小浪底水库投入使用造成；故结合此前的研究成果，将花园口径流量和泥沙量分为3 个时期，分别是1950-1979，1979-1999，2000-2011。在1950-1979年间，花园口平均年径流量达到454.89 亿m3，平均年沙量为12.85 亿t，1980-1999年间，花园口平均年径流量为334.93 亿m3，平均年沙量为7.31 亿t，年流量较前一时间段变化-29.4%，年沙量较前一时间段变化-34.7%；2000-2011年，花园口的平均年径流量为233.78 亿m3，平均年沙量为0.99 亿t，年流量较前一阶段变化-33.1%，年沙量较前一阶段变化-88.6%。

图5 黄河下游水沙变化过程Fig.5 Change process of water and sediment in the lower Yellow River

利津站为黄河入海口的关键控制水文站，可将其径流量和泥沙量与花园口站相似的分为1950-1979，1979-1999，2000-2015；在1950-1979年，利津站的年平均流量为431.17 亿m3，年平均输沙量为11.01 亿t。1980-1999年间，利津站的年平均流量为213.34 亿m3，年平均输沙量为5.14亿t，年平均流量较上一时期变化-50.0%，年平均输沙量较上一时期变化-47.7%；2000年-2015年间，利津的年平均流量为159.62 亿m3，年平均输沙量为1.31 亿t，年平均流量较上一时期变化-34.11%，年平均输沙量较上一时期变化-47.7%。

3.2 河口三角洲面积及河长变化

河口三角洲的发展过程受到水沙条件和海洋动力的叠加影响，随着大量上游来沙输送至黄河河口，河口海岸发生淤积外延、河道增长、比降变缓等。而伴随着侵蚀基准面的抬高，输沙能力降低，河口的淤积段将不断抬升，并逐渐向上发展，进而抬高洪水水位，造成严重的洪水灾害。

为了控制三角洲演变导致的灾变问题，黄河入海口三角洲的面积变化及河长演变问题显得尤为重要。本文使用了MND⁃WI 方法对Landsat 卫星遥感数据进行了计算处理，通过计算1996年后河长的变化情况和新淤积的河口三角洲面积，如图6所示，得到了近年来黄河入海口的三角洲面积变化趋势以及河长变化，所选取的河长计算起始断面为清8断面。

图6 面积计算区域Fig.6 Calculation area

表1 不同时间段花园口和利津水文站的径流与输沙量变化Tab.1 Runoff and sediment transport at Huayuankou and Lijin stations in different time periods

图7 显示了1985-2015年间黄河三角洲面积变化情况，由于考虑到人工清8改汊工程在1996年实施，2003年小浪底水库到达预期水位，调水调沙试验正式启动，将研究区域段分为三个时期，分别为1985-1996年、1996-2003年，2003-2015年。在1985-1988年期间，黄河来水来沙量的锐减导致黄河下游河口段海岸普遍遭受侵蚀，河口段三角洲面积减少57 km2，年均减少19 km2，但自1988年起始，河口段三角洲因冲积后近似楔形状，黄河尾闾受大堤影响，河道顺直，河口段一直处于淤积造陆阶段，淤积面积逐年增加，1988-1992年淤积面积增加116 km2，年平均增加29 km2。自1992年后，黄河来水来沙量减少，海岸线变化幅度不大，黄河三角洲淤积处于波动增加的状态，至1995年增加了淤积面积62 km2，年平均增加淤积面积20 km2。1996年黄河清8 改汊工程将黄河改道从北汊入海，重新生成了新河嘴，由于来水来沙量降低，新河嘴淤积缓慢，而旧河嘴不断受到海水侵蚀，黄河河口三角洲总体淤积面积处于小幅波动淤积上升趋势。1999年小浪底水库投入使用之后，入海泥沙量骤减，新河嘴遭受海水侵蚀，面积减少［20］。2003年后，受到调水调沙试验以及黄河中上游以及黄土高原水土保持措施的影响下，入海水沙量逐年减少，新河嘴开始逐步缓慢淤积，而旧河嘴不断退蚀，直至2015年河口段三角洲面积增加25.6 km2。

图7 黄河三角洲面积变化Fig.7 Change of Yellow River Delta area

新河口三角洲的河长及淤积面积计算范围如图8所示。图9 和图10 分别为黄河口改道后新三角洲面积变化和新河长变化情况。

图8 新河口三角洲河长及其淤积面积计算范围Fig.8 Study area of the length of the new estuary delta and its siltation area

图9 黄河口改道后新生成三角洲面积变化情况Fig.9 Area change of the newly formed delta after the diversion of the Yellow River estuary

图10 黄河河口改道后新生成河长变化情况Fig.10 Change of the length of newly formed river after the diversion of the Estuary of the Yellow River

黄河河口三角洲的面积和黄河新河口三角洲河长变化与时间的拟合曲线分别为：

其中S和L分别表示河口三角洲的面积和河长；从图9 得知，黄河河口三角洲的面积变化与时间的相关关系呈抛物线型，面积增长逐渐加快。河长变化呈现指数型函数，增长逐渐放缓。从遥感数据解析中也可得知，黄河河口三角洲的淤积发展随着时间的推进，淤积的面积越来越大；图10也表明在改道初期，河口三角洲将会经历一段河长快速增长的过程，但是河长的增长随着新河床的稳定逐渐放缓。图11 表明了黄河河口三角洲淤积面积与河长的相关关系，在新三角洲发育的后期，黄河入海口的河长增长引起的面积变化逐渐增大，由于新三角洲在改道初期快速发育，三角洲的发育情况以河道延伸为主。虽然三角洲河长发育变缓，但两侧开始逐渐冲积出新的陆地，新三角洲的面积增长速率加快。

图11 黄河新河口三角洲面积与河长关系Fig.11 Relationship between delta area and river length of the new estuary of the Yellow River

经过前述计算，新生成的河嘴面积变化情况如图11 所示；前文计算的河口三角洲情况根据式（14），可以得出在河口改道之后，原三角洲自然退蚀的变化情况（见图12）。

图12 改道后原河口侵蚀变化情况Fig.12 Changes of erosion in former estuaries after diversion

式中：SO为剔除了新河嘴生成增加的面积部分后剩余的三角洲面积，这一部分因为只受到海洋动力作用因而将其用于分析海洋动力作用对三角洲演变的影响；ST为图5 所计算的黄河三角洲面积；SN为三角洲在改道后新淤积生成的面积。

旧河嘴因人工改道废弃使用后，其原先淤积的三角洲滩地受到海洋冲刷作用而逐渐发生退蚀。在改道之后的19年里，消失的陆地面积为77 km2，年均侵蚀面积4.05 km2。

3.3 河口三角洲与入海水沙关系探究

河口三角洲的演变与上游的来水来沙情况、河床的变化情况有着非常重要的关系。探究河口三角洲的演变参数，并构建滞后响应模型建立相关关系，为定量研究黄河入海口的演变变化提供参考。

3.3.1 利津断面平均高程

本文首先通过公式（2）计算利津断面的平均高程，如图13所示。黄河入海河道的三次大型改道分别出现在1953年、1964年、1976年，这与图13 中观察的陡降趋势相应对。利津河道平均高程在每次改道之后都会迅速下降，新河道的河长短，比降大，过流的流速变快，水位低，水流的冲刷携沙能力增强，河道遭受明显冲刷；但在改道之后经过1～2年时间，河道高程又将有一个迅速回升过程，发生这种现象的原因主要是由于黄河入海口的淤积增长速率降低，河长开始增长，河道比降趋于平稳，河道泥沙的淤积量大于冲刷量。20 世纪80年代，清水沟河口严重淤积，为了油田的安全，通过疏浚工程将河口三角洲和河道里的泥沙进行清理，导致利津河道高程明显下降。

图13 利津断面平均高程Fig.13 Average elevation of Lijin

3.3.2 滞后响应模型计算

黄河河口河长与利津的累计来沙量之间通常有较好的相关关系，河口河长与累计泥沙淤积量呈线性相关［20，21］，我们将使用线性方程来拟合他们之间的相关关系。

黄河尾闾河道1996年进行人工改道，原河路从清水沟改道汊8 入海，而在2007年黄河尾闾又发生自然出汊，河道发生改道。图14显示了通过遥感数据提取的1996年和2007年两次河道改道河长数据与利津累计来沙量之间的相关关系。

图14 累计输沙量与新河长相关关系Fig.14 Relationship between cumulative sediment transport and new river length

利用灰狼优化算法，求解得到相关参数，选择选择1980年到2000年数据为训练集，2001-2012年的数据为测试集。灰狼优化算法参数初值设定参考了前人的研究内容［16］，参数初值范围设置分别为（0，1）、（-1，0）、（0，1）、（0，1）。设置种群大小为100，迭代次数为500 次，适应度函数选择纳什系数NSE。表2为利津站水位灰狼算法拟合参数表，而图15则是适应度函数的变化曲线图。将计算河长公式代入式（8）计算得到利津3 000 m3/s流量水位，得到的结果如图16。

图15 利津站水位适应度函数的变化曲线Fig.15 The variation curve of water level fitness function at Lijin Station

图16 3 000 m3/s流量对应的水位计算结果Fig.16 Calculation result of flow level of Lijin 3 000 m3/s

表2 利津站水位灰狼算法拟合参数Tab.2 Fitting parameters of water level grey Wolf algorithm in Lijin Station

滞后响应模型对利津的3 000 m3/s 流量水位进行计算有着较好的拟合效果，在趋势上基本保持一致，相关系数R2=0.83，并且很好的拟合了水位的变化趋势，可以较好的拟合实测结果。图17为利津和西河口的河道比降计算结果。

图17 河道比降计算结果Fig.17 Calculation result of channel slope

通过式（7），得到河道的比降数值，模拟结果与实测结果R2=0.83，模拟效果很好。可以明显发现，利津河道比降的比降最低值为0.11，并在未来的几年时间内保持稳定。西河口河道比降在改道后同样也会逐年降低，并且于0.145 出达到最低值，保持在最低值波动，而随着再次改道后，又将重复循环这一过程。该结果与郑珊提出的改道点以下流路平衡比降约在0.15～0.16的结果十分接近［13］。如果能够合理的预测未来水沙数据以及河道特征演变情况，本文所提出模型将很好的预测未来的水位和河道比降变化，进而对河道的特征变化及改道风险做出可靠判断。

4 结 语

（1）河口三角洲的海岸线变化频繁，发生了多次的淤进、侵蚀等现象。黄河的来水来沙条件是黄河三角洲演变的物质基础，通过本文的分析可以得出：黄河三角洲的水沙关系特征变化近年来十分明显，自20 世纪80年代起，随着退耕还林、水库兴建等措施，进入河口的水沙量大幅度削减，2000-2015年利津站平均输沙量相对1950-1970年减少79.1%。

（2）通过遥感数据对三角洲的淤积和侵蚀进行了分析，发现黄河入海口的冲淤演变同时受到了来水来沙条件和河口改道变迁的影响，在河道改道初期，新生成的河道河长与经历快速增长阶段和平稳增长阶段，随着河长的增长到一定程度，三角洲面积增长也将放缓。年均发生造成侵蚀大约为4 km2。

（3）河口三角洲的演变主要是由来水来沙条件变化以及河口改道造成的变化所影响，考虑到他们之间复杂的物理作用机制，采用滞后响应模型可以较好地拟合利津站3 000 m3/s 流量对应的水位以及相应的河道比降。为定量研究黄河冲淤演变提供参考和建议。 □