灞河流域水环境数值模拟研究

陈 军，吴 程，李 爽，王 颖，龚 毅

(1 中国电建集团市政规划设计研究院有限公司，广东 珠海 519000； 2 中国电建集团昆明勘测设计研究院有限公司，云南 昆明 650032；3 西安理工大学 水利水电学院，陕西 西安 710048)

我国水资源时空分布不均，随着全球气候变化的加剧，水资源短缺与水环境恶化问题愈显突出。城市水系周边人口密集、工企业聚集、污染源众多，人口快速增长与城市化进程加快造成了一系列城市水环境问题[1-2]。城市水系承担着调蓄防洪、维持区域生物多样性、气候微调和景观文娱等功能[3]，是城市经济发展和生态平衡维系的重要依托。如何对城市水系水环境进行高效管理是目前研究热点之一[4-5]。

作为水环境问题快速识别与高效管理的重要技术手段之一，水环境模型得到了广泛应用[6-10]。水环境模型按照计算域尺度可分为水域模型与流域模型，其中水域模型最具代表性的有DHI MIKE21/MIKE3、EFDC(environmental fluid dynamics code)和Delft3D等，流域模型最具代表性的有SWAT(soil & water assessment tool)和HSPF(hydrological simulation program-fortran)[11-13]等。目前相关研究以水质指标空间分布与典型控制断面(点位)的指标变化特征以及水环境对各类工程响应特性的研究为主[14-16]，但利用水环境模型对水质数据进行空间扩展后再进行综合水质评价的研究较少。

灞河是西安市第二大河流，下游城区段受乡村、城镇等典型人类活动的影响，在浐灞橡胶群坝工程实施后，河道内形成了较为稳定的水域面积，水体流速减缓，使得河道流态更加接近湖泊，形成了典型的兼具河道与湖泊特征的湖泊型城市河道。前期已有学者利用水环境模型针对灞河做了如下研究，Li等[17]利用SWAT模型研究了灞河流域非点源污染空间分布特征且识别了重点源区；胡胜等[18]基于数字滤波法和SWAT模型，借助空间分析与全局趋势分析方法对灞河流域年均基流空间分布特征进行了解析，为流域地下水资源开发与保护提供了决策依据。综上，目前关于水环境模型在灞河层面的应用研究多以灞河流域水文和非点源污染特征为主，较少涉及到灞河和浐河水域尺度。

本研究以灞河流域下游城区段灞河和浐河水域为对象，利用DHI MIKE21FM模型中HD模块耦合Ecolab模块，选取溶解氧(dissolved oxygen,DO)、化学需氧量(chemical oxygen demand,COD)、叶绿素a(chlorophyll a,Chl-a)、氨氮(ammonia nitrogen,NH3-N)和总磷(total phosphorus,TP)5个水质指标作为状态变量，参考WQ模板中各状态变量的生化反应方程对其修改后进行水环境模型构建，利用监测数据对模型校准后，开展灞河和浐河水环境数值模拟研究，分析灞河口水质指标变化特性。基于数值模拟结果利用BP(back propagation)神经网络算法对研究区时均综合水质类型进行评价，解析基准年、丰水年、平水年与枯水年研究区时均综合水质类型空间分布特征，利用主成分分析(principal component analysis,PCA)法识别了研究区水质类型的主导性指标，以期为湖泊型城市河道水环境数值模拟研究提供参考。

1 研究区域概况

灞河流域位于陕西省中部，西安市东南部，发源于秦岭北麓蓝田县灞源乡九道沟，自东南流向西北，最终于灞桥区三郎村附近汇入渭河。灞河总长104.1 km，河床平均比降6.2‰，流域面积2 581 km2，流域内现有罗李村和马渡王2座水文站[19]。灞河主要支流有浐河、蓝桥河与辋峪河。浐河系灞河下游最大支流，全长64 km，河床平均比降8.9‰，流域面积752.8 km2。灞河流域属温带半湿润大陆性季风气候，多年平均气温13.3 ℃，多年平均降水量720 mm，多年平均蒸发量946.6 mm，降水时空分布不均，降水量由南向北递减，降水主要集中在6-10月。灞河流域内有5区1县，分别为商州区、长安区、雁塔区、灞桥区、浐灞生态区和蓝田县[20]。流域上游地类以林地为主，中游以耕地为主，下游以建设用地为主。灞河是西安市重要的饮用水源地之一[19]，也是西安市经济与生态发展的重要载体。研究区以灞河流域下游城区段内灞河和浐河为主，南起陇海线北至灞河入渭河口，沿水流方向自东南向西北贯穿浐灞生态区核心区。灞河流域及研究区域位置如图1所示。

图1 灞河流域位置及各监测点位的分布Fig.1 Distribution of Ba River basin and monitoring point

2 研究数据来源

(1)水文数据。1958-2012年水文数据由西安理工大学图书馆馆藏《黄河流域水文资料》(纸质版)整理得到。(2)气象数据。气象数据由中国气象数据网(http://data.cma.cn)获取。(3)水陆边界数据。利用2016年Sentinel-2两景影像(影像数据来源于美国地质勘探局USGS网站https://earthexplorer.usgs.gov，影像ID：S2A_OPER_MSI_L1C_TL_SGS__20160429T034347_20160429T084621_A004448_T49SCT_N02_01_01、S2A_OPER_MSI_L1C_TL_SGS__20160429T034347_20160429T084621_A004448_T49SCU_N02_01_01)，采用开源GIS工具中最大似然法解译后再根据Google卫片调整后获取。(4)水下地形数据。将现场实测水深、水位数据导入GIS软件后，将水深与水位数据转换成河底高程数据，再与河道岸线高程数据合并后，得到整个研究区水下地形数据。(5)水工建筑物数据。来源于文献[21]。(6)水质数据。于2012年7-10月开展了5期水质监测(2012-07-10，2012-07-25，2012-08-24，2012-09-25和2012-10-15)，其中DO与水温2项指标由便携式多参数水质监测仪进行检测，COD、NH3-N、TP和Chl-a 4项指标均参照《水和废水监测分析方法(第4版)》[22]进行检测。(7)点源污染数据。研究区共有6个污染点源，其流量(于研究点源处多次测量10 L容器注满点源废水所需时间，对时间取算术平均值后换算得到流量数据)和水质(DO、水温、COD、NH3-N、TP和Chl-a)与水质数据(6)同期检测。(8)水深数据。由便携式超深波水深仪测定，与水质数据(6)同期测定。

3 水环境模型构建与数据处理

3.1 水环境模型构建

利用DHI MIKE21FM模型中的HD模块耦合Ecolab模块，选取DO、COD、Chl-a、NH3-N和TP共5个水质指标作为状态变量，参考WQ模板中各状态变量的生化反应方程并对其修改后进行水环境模型构建，修改后各状态变量的生化反应控制方程为[23]：

(1)

(2)

(3)

(4)

(5)

式中：C1表示DO的质量浓度；A表示大气富氧率；S表示光合作用富氧率；B表示呼吸作用耗氧率；CSC2表示COD降解过程中耗氧率；CSC4表示NH3-N降解过程中耗氧率；C2表示COD的质量浓度；DC2表示COD降解率；C3表示Chl-a的质量浓度；PC3表示Chl-a净产率；DEC3表示Chl-a死亡率；C4表示NH3-N的质量浓度；PC2-C4表示COD降解过程中NH3-N产率；DC4表示NH3-N降解率；C5表示TP的质量浓度；PC2-C5表示COD降解过程中TP产率；DC5表示TP降解率。

3.2 水环境模型设置

研究区计算域2个起点分别为灞河6号坝下游100 m处和浐河4号坝下游50 m处，终点为灞河1号坝下游10 m处(图1)。采用非结构化网格进行计算域网格划分，经网格无关性分析后，本次数值模拟最优网格数量为9 254个。模拟时间设置为2012-07-10-10-15，时间步长取10 s。上游水动力边界统一采用流量边界，下游采用水位边界。上游水质边界均采用各指标实测浓度输入，下游水质边界采用第二类边界条件(诺依曼边界条件)。模型中点源均以简单源进行概化，水动力输入条件统一采用流量输入，水质输入条件统一采用浓度输入。由于监测数据期数限制，故本研究未将实测数据分为率定与验证两个阶段，而是将率定与验证阶段进行合并。

3.3 计算工况设置

本研究共设置基准年、丰水年、平水年与枯水年4个计算工况，基准年水动力边界条件采用2012年《黄河流域水文资料》(纸质版)中灞河和浐河的流量与水位数据，水质边界条件采用2012年DO、COD、Chl-a、NH3-N和TP质量浓度的实测数据，丰水年、平水年与枯水年的水质边界条件与基准年保持一致，水动力边界条件分别为丰水年、平水年与枯水年相应流量。丰水年、平水年和枯水年相应流量由以下步骤计算：对马渡王水文站1958-2012年的年平均流量进行P-Ⅲ曲线拟合后，选取Cs=2Cv(其中Cs表示偏态系数，Cv表示变差系数)曲线确定丰水年(保证率P=10%)、平水年(保证率P=50%)和枯水年(保证率P=90%)年平均流量，选取典型年后对年平均流量进行年内分配，分别得到丰水年、平水年与枯水年的水动力边界条件。

3.4 BP神经网络训练

为对研究区水质类型进行准确评估，引入BP神经网络法算法，依据地表水环境质量标准(GB 3838-2002)选取DO、COD、NH3-N和TP的标准限值，利用python3语言numpy包中的random模块下的uniform函数生成4×300的数组作为BP神经网络输入层，并将1～6的十进制数值进行二进制转化后生成1×300的数组，作为BP神经网络输出层，进行BP神经网络训练，其中隐藏层层数采用黄金分割法[24]最终确定为13，输入层与隐藏层之间的激励函数采用tansig函数(双曲正切S型传输函数)，隐藏层与输出层之间的激励函数采用purelin函数(线性传输函数)，样本训练过程中DO中Ⅰ类水标准限值定为0 ℃下的水体饱和溶解氧质量浓度(14.64 mg/L)。

3.5 数据处理

数据统计分析、图形可视化表达以及BP神经网络训练均采用python3语言。计算结果云图采用DHI MIKE自带的Plot Composer进行处理。

4 结果与分析

4.1 模型精度验证

由于监测数据完整性的限制，因此选取图1中的S1、S4点位的水深及S1、S2、S3和S4点位的DO、COD、NH3-N、TP与Chl-a数据对水环境模型模拟值与实测值进行对比分析，为方便后续统计，本研究将研究时段2012-07-10-10-15转化为以天计数，即2012-07-10为第1天，共98 d，水环境模型模拟值与实测值的对比结果见图2。从图2可以观察到，各指标的实测值基本分布于模拟值两侧，且模拟值的变化趋势与实测值基本保持一致。

图2 灞河流域S1、S2、S3和S4监测点位水环境模型模拟值与实测值的对比Fig.2 Comparison of simulated and measured values of water environment model at S1,S2,S3 and S4 monitoring points of Ba River basin

为了对灞河流域水环境模型精度进行量化分析，引入百分比偏差系数(PBIAS)[25]。当PBIAS<0时表示模拟值较实测值偏大，PBIAS>0时表示模拟值较实测值偏小。PBIAS的计算结果如表1所示。

表1 灞河流域水环境模型百分比偏差系数(PBIAS)的计算结果Table 1 Percent bias of water environment model at Ba River basin %

表1显示，除了S1点位的NH3-N(-25.41%)和S3点位的TP(25.06%)外，各项指标的PBIAS绝对值均小于25%，这表明模型模拟值与实测值拟合趋势良好[26]，已建立的模型较好地反映了研究区水质指标的变化情况，可用于该区域不同工况下水质指标的模拟计算。

4.2 灞河口水质分析

根据《陕西省水功能区划》[27]规定，研究区水质目标为地表水Ⅳ类(GB 3838-2002)，相对渭河而言，灞河口水质需满足Ⅳ类目标。本研究计算域未包含灞河1号坝下游至灞河口段，考虑到未纳入计算域的河段长度仅约为1.7 km，且该河段基本无外源污染扰动，因此选取该区域出口断面水质代表灞河口水质进行分析。不同工况下灞河口水质指标变化趋势如图3所示。

图3 不同工况下灞河口断面DO、COD、NH3-N和TP质量浓度的变化Fig.3 Variation of DO,COD,NH3-N and TP concentrations at the Bahekou section for different scenarios

图3表明，TP与DO在研究时段内均满足Ⅳ类水质标准，COD与NH3-N部分时段未达Ⅳ类水质标准，这与文献[28]的论述一致。COD和NH3-N的超标时段主要集中在1～31 d(7月至8月初)。

对灞河口断面COD与NH3-N 2项指标超标情况进行统计，结果见表2。

表2 不同工况下灞河口断面COD和NH3-N质量浓度超标情况统计结果Table 2 Statistics of COD and NH3-N concentrations exceeding standard at the Bahekou section for different scenarios

从表2可以观察到，2012-07-10-10-15共计98 d内，丰水年与平水年COD超标时间少于基准年与枯水年。丰水年与平水年NH3-N超标时间分别为22和27 d，也小于基准年与枯水年。丰水年COD与NH3-N最大超标率大于平水年与基准年，这是因为在水质边界条件相同时，由于丰水年水动力边界的流量在部分时段小于基准年与平水年，且该时段点源流量相对较入口流量大。基准年NH3-N的超标时间明显大于COD。分别对基准年计算时段内研究区域的COD与NH3-N取算术平均值，由这2项指标的均值明确的水质类型与西安市生态环境局公布的《2012年第三季度环境质量监测分析》中的水质类型一致。

4.3 水质类型空间特性分析

研究水质空间分布特性是识别水域重点保护区的重要途径之一。鉴于上述原因，提取水环境模型计算结果中DO、COD、NH3-N和TP共4项指标后，对各指标在时间尺度(2012-07-10-10-15)上取算术平均值，然后利用前期训练好的BP神经网络对研究区综合水质类型进行计算，该计算结果称为时均综合水质类型。不同工况下研究区时均综合水质类型的空间分布如图4所示。从图4可以看出，浐河和灞河在交汇前，灞河水质整体优于浐河水质，灞河水质类型主要为Ⅱ类，部分区域为Ⅲ类；浐河水质类型主要为Ⅳ类，部分区域为Ⅲ类。灞河3号坝下游水质类型除丰水年为Ⅲ类外，其余工况下均为Ⅳ类。不同工况下，尽管2012年7-10月灞河4号坝下游部分区域时均综合水质类型超过Ⅳ类，鉴于其影响范围较小，且该范围内水质下降的主要原因在于河道左岸点源污染的影响，西安浐灞生态区水生态系统保护与修试点工作中的污水治理与中水回用工程实施后，已消除了该范围内的点源影响，因此认为研究区时均综合水质类型在空间范围内于基准年、丰水年、平水年和枯水年基本都满足Ⅳ类水质目标要求。

图4 不同工况下灞河时均综合水质类型的空间分布Fig.4 Spatial distribution of time-average comprehensive water quality at Ba River

由图4可知，受浐河水体扰动，灞河水体在与浐河水体交汇后，交汇区水质类型出现平面梯度分区，水质类型梯度分区从左岸到右岸依次为Ⅳ类-Ⅲ类-Ⅱ类。交汇区水质类型分布特征主要受灞河与浐河交汇前的动量比率以及两水体各项水质指标影响。交汇区水质类型分布特征在不同工况下表现出明显差异，灞河入交汇区前，丰水年动量大于基准年，因此丰水年条件下交汇区Ⅲ类水的带状面积较基准年条件下有所增加，且Ⅲ类水的条带明显向左岸推进。枯水年条件下灞河入交汇区前动量比率小于基准年条件下，交汇区Ⅲ类水条带明显向右岸推进。

为分析不同工况下研究区内水质类型的主导性指标，提取研究区内DO、COD、NH3-N和TP共4项指标，对各指标在时间尺度(2012-07-10-10-15)上取算术平均值，分别对不同工况进行PCA分析。由于KMO(Kaiser-Meyer-Olkin)值均大于0.5且Bartlett显著性概率均为0.000，满足PCA分析条件。根据特征值λi>1的原则提取的1个主成分反映了原始4个指标57%以上的信息。该结果表明，不同工况下COD与NH3-N对第一主成分的贡献率为63%～64%，可以认为不同工况下研究区时均综合水质类型的主导性指标主要为COD与NH3-N。

5 结 论

将DHI MIKE21FM中HD模块与Ecolab模块耦合后构建了适用于灞河流域的水环境模型，研究了灞河口水质指标的变化情况，基于BP神经网络算法结合模拟结果对研究区时均综合水质类型进行了评价，借助PCA识别了研究区时均综合水质类型的主导性指标，得到的研究结论如下：

1)灞河口水质超标因子主要为COD与NH3-N，2012年7-10月NH3-N的超标时间与最大超标率明显大于COD。丰水年与平水年COD与NH3-N的超标时间较基准年明显下降，但最大超标率增加。

2)灞河和浐河交汇前，灞河水质整体优于浐河，灞河水质类型主要为Ⅱ类，部分区域为Ⅲ类；浐河水质类型主要为Ⅳ类，部分区域为Ⅲ类。除丰水年灞河3号坝下游水质类型为Ⅲ类外，基准年、平水年与枯水年该区域水质类型均为Ⅳ类。研究区时均综合水质类型基本满足Ⅳ类水质目标，水质类型主要受COD与NH3-N 2项指标影响。

3)所构建的水环境模型，经验证后具有较高的精度，可适用于类似区域的水质预测及评价。

需要指出的是，由于水质监测数据时长限制，基准年仅计算了7-10月水质指标的变化情况，计算时段并未覆盖1个水文年。因此，本研究得出的结论主要适用于灞河流域雨季的水环境特征，时均综合水质类型并不能代表全年平均的综合水质类型。水环境数值模拟是一个多情境与多维度的过程，本研究在水质模型构建过程中，以实测数据表征了上游的点源与非点源的污染特征，但同时也忽略了下游非点源输入，后续还需进一步研究。