嫩江下游洪水演进及对洪泛区植被影响分析

毛思媛，贾艳红，假冬冬，郑祥民，范宝山，张佳欣

(1.华东师范大学地理科学学院地理信息科学教育部重点实验室，上海 200241；2.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏 南京 210029； 3.中水东北勘测设计研究有限责任公司，吉林 长春 130061；4.水利部寒区工程技术研究中心，吉林 长春 130061)

随着全球极端气候事件的增加，洪水频率也相应增加[1]，这不仅对生命财产安全构成了重大威胁，也对河岸生态系统产生了重大影响。Berthelot等[2]研究表明，洪水频率的增加将导致不同洪水重现期内物种的组成、多样性及群落结构出现显著差异，特别是在洪水频发区，植被的密度、更新率与物种多样性会显著减少。为有效减轻洪涝灾害造成的巨大损失，促进河流生态建设，模拟流域洪水演进并分析洪水淹没对河岸生物群落的影响对恢复河流基本功能、改善河流水生态环境及促进河流生态治理具有重要意义。

近年来，随着洪水管理及生态河道治理政策的推行，如长江流域“共抓大保护，不搞大开发，生态优先，绿色发展”战略、黄河流域“积极实施生态保护和高质量发展”战略，松花江流域的洪水管理与生态环境保护也日渐受到人们的重视。目前已在嫩江洪水管理与洪水演进方面展开了一些研究，如孙亚楠等[3]利用MIKE11水动力学方法模拟了不同洪水组合下嫩江受第二松花江洪水顶托的影响,杨浩淼[4]基于HydroInfo二维水动力模型模拟了嫩江湾湿地拟建道路对河道防洪的影响,李洋[5]采用HEC-RAS一二维水动力耦合模型对嫩江左岸齐齐哈尔以上防洪保护区进行了洪水风险图模型的构建与验证。在洪水对河岸生态系统的影响方面国内也有部分研究成果，如饶品增等[6]利用Landsat 7/8数据和MNDWI方法对2013年黑龙江省洪水对植被的影响进行了评估；余晓等[7]通过MIKE11一维非恒定流水动力学模型模拟了额尔古纳河洪水淹没范围并分析了湿地植被变化情况。目前关于不同频率洪水运动特征及其对生态系统影响的研究相对较少。

不同重现期的洪水淹没程度显著影响着河岸植被的群落分布。本文在利用水动力学方法建立嫩江流域河道洪水演进模型的基础上，进一步模拟了典型河段4种不同重现期的洪水运动特征，并选取MODIS卫星的MOD13Q1植被指数数据说明不同重现期洪水运动特征及淹没范围具有的生态学意义，旨在一定程度上弥补当前河道管理与生态河道建设缺乏规范与定量性研究的不足，同时也为河道实现近自然治理的目标提供基础数据。

1 研究区概况

嫩江是松花江的最大支流，位于黑龙江省西北部，地处温带大陆性季风气候区，夏季高温多雨，降雨集中，上游河网密度大，支流多，河道全长1 370 km，流域面积29.7万 km2，是我国水旱灾害频发的地区[8]。嫩江江桥—大赉河段位于嫩江下游三岔河口以上平原地带(图1)，河长约225 km[9]，主要流经黑龙江省齐齐哈尔市、大庆市与吉林省白城市，主河道布设有江桥、大赉两个水文测站，以及位于两水文测站间的白沙滩水位站(123°53′E, 46°18′N)。该河段为典型的平原性河道，地形复杂，主流宽度小，河滩宽度大，平均坡降约0.0045‰～0.005‰[9]。河道水流平缓，河槽蜿蜒曲折，河道内沙洲、汊河密布，多呈网状，且沿江两岸滩地上广泛分布着大片湿地、泡沼群及牛轭湖。大洪水期间，江道与泡沼相连，两岸滩地洪水漫溢，最宽处可达20.5 km，最大水深达7.4 m，河槽调蓄量大，水位涨落缓慢，洪水历时长，使得汛期两岸堤防的防洪任务较重[10]，也显著影响着沿岸地区的植被生长。

图1 研究区概况

2 一、二维水动力耦合模型构建

MIKE FLOOD一、二维耦合模型是把一维模型和二维模型连接在一起进行动态耦合的模型系统。耦合模型既利用了一维模型和二维模型的优点，又避免了采用单一模型遇到的网格精度和准确性问题[11]。基于此，本文采用丹麦水资源及水环境研究所(DHI)研发的MIKE FLOOD进行河道洪水演进模拟，模型以非恒定流运动的圣维南方程组(Saint-Venant)作为河道洪水模拟计算的基本方程[12]。

2.1 MIKE11一维水动力模型的构建

MIKE11一维水动力模型[13]的建立主要由河网数据、断面数据、边界条件及水动力学参数等组成，具体参数设置如下：

a.河网数据：根据2003年绘制的嫩江流域1∶5 000地形图，本文对江桥—大赉河段河流走向线进行描绘，并将其概化生成河网文件。

b.断面数据：利用ArcGIS对高精度的DEM数据进行断面提取，由于嫩江江桥—大赉河段河长超200 km，考虑到模型运行时的稳定性及准确性，本文共选取65处典型断面，各断面间距约1 000～4 000 m，并对各断面最低点与左右岸进行标注。

c.边界条件：选取江桥站1998年、2013年与不同设计频率(1%、2%、5%、20%)的洪水流量过程作为模型的入流边界条件，选取大赉站同期的水位过程及水位-流量关系作为模型的出流边界条件。由于研究河段长、洪水传播时间长，本文以各年份汛期(6—10月)作为计算时段，历时约123 d，各测站均采用日水位-流量资料。

d.水动力学参数：根据《松花江流域防洪规划》和实地查勘情况，本文将初始水位设为130.65 m，曼宁系数初始设置为0.035，模拟时间步长为5 s，输出结果保存为.res11文件。

2.2 MIKE21 二维水动力模型的构建

MIKE21 FM 水动力模块的建立主要包含网格剖分、地形构建、糙率分区及水动力参数(干湿边界、涡黏系数)的设定等。

a.网格剖分：本文采用非结构三角网格对建模区域进行划分，并根据研究区域地形起伏灵活确定网格大小，如对河堤处进行网格加密。本次模型共划分网格96 886个，节点52 573个。

b.地形构建：将高程散点数据导入网格中进行插值，生成地形文件。

c.糙率分区：嫩江流域共有耕地、林地、草地、湿地、水体、人工表面及裸地7类地物，地物类别数据主要来源于全国地理信息资源目录服务系统GlobeLand30数据集 (http://www.webmap.cn/commres.do?method=globeIndex)。糙率的选取主要参考《松花江流域防洪规划》与陈晓更等[14]的研究，并结合现场实地查勘与工程经验，在模型率定过程中进行相应的调整。

d.水动力参数: 由于计算区域中存在随水位涨落而变化的动边界，为保证模型计算的连续性，采用“干湿水深处理方式”，将干、湿水深及淹没水深分别设置为0.005m、0.01m与0.05m，涡黏系数根据Smagorinsky公式设置为0.28。

2.3 一、二维水动力模型的耦合

本文采用MIKE FLOOD将MIKE11和MIKE21 FM进行耦合，耦合时采用侧向连接来模拟洪水从河道漫流至洪泛区的运动过程。连接方式允许MIKE21 FM的一系列网格单元从侧面连接到MIKE11的部分河段甚至整个河段[11]。连接方法采用“Cell to Cell”，即对结构物中每1个节点都进行水动量方程计算，之后所有经过计算的水流都被重新分配到MIKE11和MIKE21 FM网格单元中[15]。

3 模型糙率的率定与验证

3.1 率定

3.1.1一维模型糙率的率定

选取1998年6—10月汛期实测水位及沿程洪痕数据用于河道糙率的率定，根据文献[16]及靳宏伟等[17-18]的松花江流域中小洪水糙率，并结合河道特性与现场查勘情况，确定模型参数，得出河道下游出流过程，然后提取河段沿程各断面水位峰值与实测洪痕值进行对比，并参考中间站白沙滩的水位拟合情况进行综合分析。经率定，确定模型计算步长为5 s，河段综合糙率为0.042 7，该糙率符合《松花江流域防洪规划》。

模型模拟结果与实测结果对比如图2所示。由图2(a)可知，嫩江江桥—大赉河段沿程各断面的水位峰值误差均小于0.2 m，该结果符合GB/T22482—2008《水文情报预报规范》[19]评定要求。中间站白沙滩的纳什效率系数(NSE)[20]、决定系数(R2)、一致性指数(d)与均方根误差(RMSE)[21-23]分别为0.98、0.98、0.99和0.174，说明模型拟合结果较好，模拟水位变化过程与实测水位变化过程基本一致，如图2(b)所示。

图2 1998年江桥—大赉河段洪痕率定

3.1.2一、二维耦合模型糙率的率定

由于嫩江流域洪水实际淹没范围矢量图难以获取，本文基于USGS官网(https://earthexplorer.usgs.gov/)下载了洪水时期研究区的Landsat影像，包括Landsat 7和Landsat 8。然后利用ENVI遥感图像处理方法中的监督分类选取洪水淹没区、未淹没区与其他样本，并通过支持向量机分类器进行分类，然后结合遥感目视解译提取的洪水实际淹没范围与模拟的淹没范围进行对比。考虑到汛期遥感影像图受云层影响较大，本文选取2013年8月11日(Landsat 8,云量低于20%)作为一、二维耦合模型糙率的率定时段。率定过程中，模型上边界选用江桥站2013年8月11日前后的实测流量过程，下边界选取大赉站同期的水位过程。经率定，嫩江流域耕地、林地、草地、湿地、水体、人工表面及裸地糙率分别为0.06、0.065、0.04、0.035、0.021、0.07和0.035，时间步长为30 s。模拟淹没范围与实际淹没范围对比结果如图3所示。由图3(a)可知，2013年8月11日洪水实际淹没范围为764.50 km2，模型模拟淹没范围为761.72 km2，模拟与实际相交的淹没范围为700.83 km2，模型相对误差为0.4%，满足《洪水风险图编制导则》[24]的要求。

图3 嫩江江桥—大赉河段率定与验证结果

3.2 糙率的验证

3.2.1一维模型糙率的验证

为保证江桥—大赉河段一维水动力模型的计算精度及普适性，选取2013年典型洪水年的实测水文数据进行验证，验证结果如图4所示。

图4 2013年白沙滩水位模拟与实测结果对比

由图4可知，2013年白沙滩实测水位峰值为136.49 m，模拟水位峰值为136.483 m，绝对误差为0.007 m，相对误差为0.01%，NSE、R2、d与RMSE分别为0.97、0.98、0.99和0.192。由此可见，2013年实测水位峰值与模拟水位峰值的绝对误差与相对误差均满足《水文情报预报规范》的峰值误差精度要求；纳什效率系数、测定系数、一致性指数与均方根误差离散程度低、拟合效果好，即模拟水位变化过程与实测过程基本一致，其中纳什效率系数符合《水文情报预报规范》中的甲级预报精度标准。因此，该模型选用的综合糙率模拟精度高，可为嫩江流域的洪水预报、防洪决策及生态河道建设提供技术支撑。

3.2.2一、二维耦合模型糙率的验证

考虑到汛期遥感影像图受云层影响较大，选取2013年9月4日(Landsat 7，云量低于20%)作为一、二维耦合模型糙率的验证时段。验证过程中，模型上边界选用江桥站2013年9月4日前后的实测流量过程，下边界选取大赉站同期的水位过程，模型模拟的淹没范围与实际淹没范围对比如图3(b)所示。由图3(b)可知，2013年9月4日提取的实际淹没范围为942.70 km2，模拟淹没范围为934.49 km2，模拟与实际相交的淹没范围为895.83 km2，模型相对误差为 1%，满足《洪水风险图编制导则》的要求。

综上所述，本文在一维模型和一、二维耦合模型计算中，选用江桥—大赉河段多个实测断面的洪痕值进行模型糙率的率定与验证，且率定与验证的持续时间较长(时间范围为6—10月)，这样可有效避免选取单一水文站进行糙率率定与验证的偶然性，同时也采用了中间站白沙滩的实测水位过程进行拟合对比，有效增加了糙率率定与验证的合理性。

4 不同重现期洪水演进及其生态学意义

4.1 不同重现期洪水淹没范围数值模拟

为了更好地促进嫩江流域生态河道建设及防洪决策的制定，本文就不同重现期洪水淹没范围进行模拟。其中上边界的流量过程主要参考《松花江流域防洪规划》中的洪峰流量进行放大，下边界的控制条件选取大赉站的水位-流量关系。利用MIKE FLOOD一、二维水动力耦合模型进行不同频率洪水淹没范围数值模拟，结果如图5所示。图5(a)为不同重现期洪水淹没范围，图5(b)为典型断面示意图(断面起点位置为124°3.740′E、45°43.906′N，终点位置为124°8.349′E、45°47.197′N)。

图5 嫩江江桥—大赉河段不同重现期洪水淹没范围

从模型计算水位分析，在5年一遇洪水模拟工况下，洪水演进至8月6日后(白沙滩附近洪水由主河槽溢出滩地)，江桥—大赉河段总淹没面积为58.50 km2，最大水位高程为138.09 m。洪水演进至8月22日后，江桥到大赉河段淹没范围最大，此时最大淹没范围为677.14 km2，最大水深为12.09 m。

在20年一遇洪水模拟工况下，洪水演进至8月6日后，江桥—大赉河段总淹没面积为192.69 km2,最大水位高程为139.75 m。洪水演进至8月22日后，江桥到大赉河段淹没范围最大，此时最大淹没范围为896.06 km2，最大水深为13.86 m。

在50年一遇洪水模拟工况下，洪水演进至8月6日后，江桥—大赉河段总淹没面积为230.46 km2,最大水位高程为140.48 m。洪水演进至8月22日后，江桥到大赉河段淹没范围最大，此时最大淹没范围为1 044.26 km2，最大水深为15.05 m。

在100年一遇洪水模拟工况下，洪水演进至8月6日后，江桥—大赉河段总淹没面积为261.58 km2,最大水位高程为141.16 m。洪水演进至8月22日后，江桥到大赉河段淹没范围最大，此时最大淹没范围为1 234.43 km2，最大水深为15.53 m。

由此可知，随着洪水频率的增加，嫩江江桥—大赉河段的洪水淹没范围也不断增加(图5)，其中20年一遇较5年一遇洪泛增加面积为218.92 km2，50年一遇较20年一遇洪泛增加面积为148.20 km2，100年一遇较50年一遇洪泛增加面积为190.17 km2。

4.2 不同重现期洪水扰动对洪泛区植被指数的影响

洪泛区位于河流与陆地环境的过渡地带，其植被的生长需要一定范围的洪水淹没来提供充足的水分[7]，但洪水淹没的频繁程度又将在一定程度上影响植被的正常生长。由于植被具有不同的光谱特征规律，目前在遥感应用领域将可见光与近红外遥感光谱观测通道进行组合运算, 便可得到衡量地表植被状况的简单、有效和经验度量的植被指数，其中归一化植被指数[25](normalized difference vegetation index，NDVI)是由红光、近红外2个波段计算生成的，其基本公式为

(1)

式中：ρNIR为近红外波段的表观反射率或方向放射率；ρRED为红光波段的表观反射率或方向放射率。

增强型植被指数(enhanced vegetation index，EVI) 是由Red、NIR、Blue(红光、近红外、蓝光)3个波段计算生成的，其基本公式为

(2)

式中:ρBLUE为蓝光波段的表观反射率或者方向放射率；G为增益系数;C1、C2为大气调节参数；L为土壤调节参数。在EVI计算中，通常取L=1，C1=6，C2=7.5，G=2.5。上述两种植被指数作为常用的植被指数，已广泛用来定性和定量评价植被覆盖及其生长活力[26]，对于指示植被的生长状态与空间分布具有重要意义。

4.2.1洪泛区植被指数季节变化分析

为深入了解嫩江流域植被指数的季节变化规律，本文利用ArcGIS对不同洪水频率淹没范围内2000—2019年的MOD13Q1 植被指数数据(时间分辨率16 d，空间分辨率250 m，全年共23景)进行掩膜[7]，并对每一景影像求像元平均值，分析结果如图6所示。由图6可见，2000—2020年不同重现期内的植被指数(NDVI、EVI)具有较强的季节性变化规律，6—8月夏季植被生长最为旺盛，7月底至8月初植被指数达到最大值，随后受洪水影响植被指数从峰值急剧下降，8—9月降低幅度最大，同时受物候条件影响植被在10月底降至完全枯黄，植被指数进入“低谷阶段”，且不同洪水频率下的NDVI和EVI的变化一致性较好。该结论与饶品增等[6]的2013年黑龙江省洪水对植被的影响评估结果一致。

图6 2000—2020年植被指数变化

图7 2012—2016年植被指数变化

4.2.2典型洪水年前后植被指数变化分析

选取2013典型洪水年前后的植被指数进行分析(图7)，该年份洪水重现期达10年一遇。由图7可知，在2013典型洪水年前后，不同重现期内的植被指数(NDVI、EVI)存在显著差异，2013年与2012年相比，2013年不同重现期内的NDVI、EVI峰值均明显降低，可见洪水对该区域的植被造成了显著的影响，其中5年一遇植被指数降低幅度最大，20年一遇较5年一遇洪泛增加区的植被指数也呈现一定的下降趋势；与2013年对比，2014年不同重现期内的NDVI和EVI的峰值均有较明显的提高，其中5年一遇的植被指数增加最多，但与2012年的植被指数峰值还存在一定差距，可见，2014年由于洪水干扰强度降低，洪泛区的植被指数均得到了明显的回升，但回升幅度尚未完全达到洪水前状况，而2015年洪泛区内的植被指数回升则逐渐接近洪水前状况。由此可见，不同重现期洪水对植被具有明显的破坏作用，但这种破坏作用在嫩江流域周期性相对较短，表明该地区植被的可恢复性较强，植被类型主要以1～2 a生为主，与饶品增等[6]的研究结果一致。

4.2.3不同重现期洪泛区植被指数变化分析

为探究不同重现期洪泛区内植被指数的变化规律，本文统计了2000—2020年植被指数的年平均变化值，如图8所示。不同重现期洪泛区内的植被指数存在显著差距，其中频繁淹没区(5年一遇以下)的植被指数最小，中等频繁淹没区(50年一遇较5年一遇洪泛增加区)的植被指数显著增大，50年一遇较20年一遇洪泛增加区的植被指数增加达到最大，而较低频繁淹没区(100年一遇较50年一遇洪泛增加区)的植被指数出现下降趋势。这说明不同重现期洪水扰动对洪泛区植被指数产生了一定程度的影响，该结果与加拿大魁北克南部不同洪水频率对物种多样性的影响结果趋于一致，同时也符合Connell[27]的中度干扰假说。

图8 2000—2020年植被指数的年际变化

根据上述洪水演进数值模拟及植被指数数据分析发现，利用不同重现期的洪水淹没范围可初步对嫩江流域洪泛区植被群落进行分区，这不仅为相关的野外植被调查提供了前期准备，也为嫩江流域防洪规划的制定及生态河道建设提供了基础数据与技术支撑。

5 结 论

a.嫩江下游一、二维水动力耦合模型在率定期与验证期的模拟效果符合甲级预报精度，模型模拟精度较高。

b.不同重现期内的植被指数具有较强的季节变化规律，夏季植被生长的最为旺盛，且NDVI与EVI的变化基本一致。同时洪水对植被具有明显的破坏作用，但这种破坏作用在嫩江流域周期性相对较短，表明该地区植被的可恢复性较强，植被类型主要以1～2 a生为主。

c.受洪水频繁淹没程度影响，频繁淹没区的植被指数最小，中等频繁淹没区的植被指数显著增大，50年一遇较20年一遇洪泛增加区内的植被指数达到最大，而较低频繁淹没区的植被指数则呈现下降趋势，该结论符合中度干扰假说。