大汶河流域泰安段水质变化及其与水量响应关系研究

侯恩光，张升第，王如岩，张春霞，杨会清，段建筑，高振勇，谭秀翠*

（1.泰安市水文局，山东 泰安271018；2.山东农业大学，山东 泰安271018）

0 引 言

【研究意义】受自然因素与人类活动的共同影响，地表河流水质发生了变化[1]，河流水质的时空变化特征研究已经成为水资源管理领域的重要工作[2]。大汶河作为山东省泰安市的母亲河，其水量丰富，但年际和年内分配不均，近些年来，受沿岸工农业快速发展的影响，部分河段水质受到明显影响，危及流域水生态系统。因此，开展大汶河流域泰安段的水质变化及其与流量的响应关系研究势在必行。【研究进展】水质与水量是水资源的双重属性，二者之间有着紧密的关系。一般来说，河流水质在1年内的变化都是由不同时期水量大小不同的差异所引起的[3]。因此，为分析不同时期水量差异所引起的水质变化，常采用季节性Kendall 检验法分析水质的变化趋势[4-5]，但该方法并未将水量与水质建立关系。采用数值模型研究水量水质关系，需要输入水动力与水质的各项参数[6-7]，其直接决定模型的计算精度。为了更真实地反映流域的水质状况，Copula 函数被应用于构建水量与水质的联合分布，分析二者的变化关系。2003年意大利学者De Michele 等[8]将Copula 函数应用于水文水资源领域，之后该方法得到广泛应用。张娜等[9]采用Gumbel-Hougaard Copula 函数构建了年最大日雨量与年最大7日雨量之间的联合分布，并基于联合分布推求设计暴雨值。杜懿等[10]研究显示二元t-Copula 函数最能体现郁江南宁水文站洪水峰量联合分布的特性规律。黄锋等[11]采用Frank Copula 函数建立年最大洪水发生时间和洪峰流量两变量的联合分布，有效挖掘了入库洪水信息。廖显薇等[12]利用Copula 函数构建二维干旱变量的联合分布。由于Copula 函数能够灵活地构造边缘分布为任意分布的水文变量联合分布，张翔等[13]将Copula 函数引入到了水量水质的研究中，构建二者的联合分布函数，并对淮河流域蚌埠闸的水量水质联合概率分布函数进行计算。【切入点】目前对于大汶河流域水质的研究主要以各类水质指标为基础，开展水质评价、纳污能力计算[14]、水生态健康评价[15]等，而对于水量影响下的水质变化规律考虑不够充分。【拟解决的关键问题】因此，以大汶河流域泰安段为研究对象，在分析水质时空变化规律及演变趋势的基础上，采用Copula函数构建二维水量（流量）与水质、三维水量（降水、流量）与水质的联合分布，计算联合概率值，明晰水质与水量的响应关系，研究成果对大汶河流域水资源管理和保护方案的合理制定、实现水资源高效利用具有重要意义。

1 材料与方法

1.1 研究区概况

大汶河是黄河下游最大的一条支流，发源于济南市钢城区，自东向西流经济南、泰安2 地，大汶河流域面积8 945 km2，年平均降水量694.5 mm，流域年平均气温12.9 ℃[16]，流域分布概况见图1。根据2018年《泰安市水资源公报》，大汶河区水资源总量15.48亿m3，占泰安市水资源总量的82%，是地区发展的重要水源。大汶河水质主要受沿岸工矿企业、居民生活区及农业面源污染影响，2018年大汶河水系评价总河长307 km，其中，Ⅲ类水占4.6%，Ⅳ类水占67.7%，劣Ⅴ类水占27.7%。《泰安市2018年度区域水功能区水质监测报告》显示，2018年泰安市主要入河排污口污水入河量为2.03 亿m3/a，其中化学需氧量（COD）占27%，氨氮（NH3-N）占4%，这2 个水质指标为泰安市水功能区限制纳污控制指标。

图1 大汶河流域分布Fig.1 Distribution of the Dawen river basin

为研究大汶河流域泰安段水质变化规律及其与水量变化的响应关系，本文采用2010—2019年大汶河流域泰安段干流及支流分布的7 个水质站点（渐汶河、北望、楼德、大汶口、南李村、障城、戴村坝）的2 个主要污染指标化学需氧量（COD）和氨氮（NH3-N）及1 个国家级水文站（大汶口）的逐月实测资料（降水量、流量）进行分析。

7 个水质站点位于泰安市水功能区划中的工业用水区，其目标水质为Ⅳ类水。

1.2 研究方法

1.2.1 季节性Kendall 检验法

月尺度的水质数据具有很强的季节性特征，因此，采用季节性Kendall 检验法对大汶河流域泰安段的水质进行趋势分析。该方法可以实现对多年同一月份的水质质量浓度进行比对，能有效地避免流量周期性变化带来的影响[17-18]，进而揭示水质随时间的变化规律及发展趋势。

当有n年p月的水质数据，则X表示为：

式中：x11，…，xnp为每月水质质量浓度数据，p≤12。

对各年第i月的水质质量浓度值进行比较，并计算Si，计算式为：

式中：1≤k

当n≥10，S服从正态分布，且标准差Z满足，

对于已知的趋势检验显著性水平α，取值0.1 和0.01，当|Z|≤Zα/2，则满足0 假设条件。

1.2.2 Copula 联合分布

Coupla 函数是定义域为[0，1]均匀分布的多维联合分布函数，可以利用它将多个随机变量的边缘分布联结起来构造联合分布函数[13]，

式中：C为Copula 函数；θ为Copula 参数；F1，F2，…，Fn为各随机变量的边际分布。

常用的Copula函数有多种类型，其中Archimedean型在水文领域最为常用，常用的函数有Clayton Copula函数、Gumbel Copula 函数、Frank Copula 函数3 种，其函数表达式见表1[19]。

表1 3 种Copula 函数表达式Table 1 Three kinds of Copula function expressions

采用下式计算联合分布函数的经验频率Pei。

式中：mi为观测值中满足X≤xi，Y≤yi（二维联合分布）或满足X≤xi，Y≤yi，Z≤zi（三维联合分布）的联合观测值的个数；n为总的数据对个数。

为选择最优的Copula 函数作为联合分布函数，采用离差平方和最小准则（OLS）、AIC信息准则来评价Copula 函数的拟合优度[20]。

OLS的计算式为：

式中：Pei、Pi分别为联合经验频率和联合理论频率；i为样本序号；n为样本容量。OLS值越小，Copula函数的拟合效果越好。

AIC信息准则计算式为：

式中：m为模型参数的个数。AIC值越小，Copula 函数的拟合效果越好。

2 结果与分析

2.1 水质变化趋势

2010—2019年，大汶河流域泰安段COD 与NH3-N 质量浓度均呈下降趋势（图2），分别由2010年的34.31、1.98 mg/L 降低到2019年的23.57、0.57 mg/L，水质呈变好趋势。在空间上，由上游渐汶河站到下游戴村坝站，COD 与NH3-N 的质量浓度呈下降趋势，但质量浓度最高的断面出现在楼德站，其位于大汶河的支流柴汶河上，其COD 与NH3-N 的质量浓度分别为31.24、1.56 mg/L，均超过Ⅳ类水质标准，受其影响，位于其下游的干流水质站大汶口站的COD 与NH3-N 质量浓度也较高。通过比较可以看出，楼德站为水污染防治的重点治理断面。

图2 2010—2019年大汶河流域泰安段水质时空变化Fig.2 Temporal and spatial variation of water quality in Taian Section of Dawen River Basin from 2010 to 2019

采用季节性Kendall 检验方法对大汶河流域泰安段的逐月水质数据进行趋势分析。

表2 为水质变化趋势分析结果，7 个站点2 个指标的14 次评价结果中，10 个评价结果呈高度显著下降趋势，1 个呈显著下降趋势，表明整体水质呈改善趋势。虽然，渐汶河站的NH3-N 质量浓度呈显著上升趋势，但其平均质量浓度较低，为0.62 mg/L，可达到Ⅲ类水质标准，未呈现污染状态。有2 个评价结果为变化趋势不明显，其中渐汶河站COD 质量浓度满足目标水质要求，但楼德站的NH3-N 质量浓度相对较高，在部分月份达到劣Ⅴ类水质，占比19%，是水污染防治的重点治理指标。

表2 水质变化趋势分析Table 2 Results of water quality change trend analysis

2.2 水质水量变化关系

河流流量具有周期性变化，河流水质组分质量浓度大多受流量周期性变化的影响[21]，而呈现出不同的变化特征。因此，本文分析大汶河流域泰安段7 个水质站点汛期（6—9月）与非汛期（其余月份）的水质变化规律。由图3 可知，大汶河流域泰安段的7 个水质站点，汛期的水质优于非汛期，汛期COD 和NH3-N 的质量浓度均达到Ⅳ类水质标准，其中南李村、障城、戴村坝3 个断面的NH3-N 的质量浓度低于0.5 mg/L，达到II 类水质标准。在非汛期，水质质量浓度相对较高，其中楼德站的非汛期的COD 与NH3-N的质量浓度均不满足水功能区目标水质要求。

图3 不同时期水质变化关系Fig.3 The relationship of water quality change in different periods

水质和水量是河流不可分割的两方面属性，质以量为载体，量的多少直接影响其环境承载力的大小[22]。因此，本文以大汶口站逐月实测降水量（P）、流量（Q）、化学需氧量（COD）、氨氮（NH3-N）为水量水质的分析指标，进行相关分析。由非汛期到汛期，大汶口站月平均降水量由22.26 mm 增大到128.04 mm，月平均流量由8.86 m3/s 增大到52.02 m3/s，COD 与NH3-N 的质量浓度分别下降6%、37%，Ⅳ类水达标月份占比由73%提高到83%。

由大汶口站水量与水质的相关分析图可以看出，见图4，随着P、Q的增大，COD 与NH3-N 质量浓度呈减小趋势，质量浓度峰值分别为45、7.66 mg/L，出现在2012年5月与2011年2月，其对应的降水量与流量均低于月平均值。在非汛期，数据点比较集中，靠近Z轴沿纵向分布，在汛期，数据点比较分散，在空间上比较接近于XY平面。综上可知，大汶河泰安段水量与水质之间存在较为密切的联系。

图4 不同时期水量与水质变化关系Fig.4 The relationship between water quantity and water quality in different periods

2.3 水量水质联合分布

河流的水量与水质是相互联系、相互影响的，为反映水量变化影响下的水质变化特征，有必要进行水量水质的联合分析，采用Copula 函数构建大汶河流域泰安段的水量水质的二维及三维联合分布，描述二者之间的相关性结构。

Copula 函数的参数估计一般采用相关性指标法、适线法、极大似然法[23]，采用相关性指标法确定二维Copula 函数参数，采用极大似然法确定三维Copula函数参数，参数取值见表3。

表3 Copula 函数参数估计Table 3 Parameter estimation of Copula functions

对常用的联合分布函数Clayton Copula 函数、Gumbel Copula 函数、Frank Copula 函数进行拟合优度检验，选择OLS 与AIC 值最小的Copula 函数作为构造水量水质联合分布的函数。由表4 和表5 可以确定，水量与水质的二维联合分布采用Frank Copula 函数，三维联合分布采用Clayton Copula 函数。

表4 二维Copula 函数的拟合优度检验Table 4 The goodness-of-fit test for two-dimensional Copula function

表5 三维Copula 函数的拟合优度检验Table 6 The goodness-of-fit test for three-dimensional Copula function

大汶口站水量与水质的二维、三维经验联合分布概率值与Copula 函数理论联合分布概率值相关关系，见图5、图6。由图5、图6 可知，对于选定的二维Frank Copula 函数和三维Clayton Copula 函数，理论联合分布与经验联合分布的拟合效果较好，R2均在0.85 以上，说明选取的二维、三维Copula 函数是合理的，可用来分析水量水质联合分布问题。

图5 二维Frank Copula 函数的联合分布概率Fig.5 Joint distribution probability of two-dimensional Frank Copula function

图6 三维Clayton 函数的联合分布概率Fig.6 Joint distribution probability of three-dimensional Clayton function

采用Frank Copula函数构建水量与水质二维联合分布，绘制概率分布图及等值线图，如图7、图8 所示。流量与COD、NH3-N 的二维联合分布规律比较类似，即当流量、水质指标质量浓度增大时，二者的联合概率值越大；当联合概率一定时，COD、NH3-N质量浓度随着流量的增大而降低，并趋于稳定；当流量一定，在任意联合概率分布下，COD 质量浓度高于NH3-N 质量浓度；等值线的疏密，表明了水质质量浓度随流量变化速率的大小，联合概率低于0.7 时，等值线较为密集，水质质量浓度变幅较小。根据等值线图，可以判定任意流量情况下，COD、NH3-N 质量浓度的发生概率，其分析结果对于水污染防治、合理开展水环境规划，实现水资源的可持续利用具有重要的实践意义。

图7 流量与COD 的二维联合分布Fig.7 Two-dimensional joint distribution of flow and COD

图8 流量与NH3-N 的二维联合分布Fig.8 Two-dimensional joint distribution of flow and NH3-N

采用Clayton Copula 函数构建水量与水质三维联合分布图，结果如图9、图10 所示。

图9 P、Q 与COD 的三维联合分布概率Fig.9 Three-dimensional joint distribution probability diagram of precipitation,flow and COD

图10 P、Q 与NH3-N 的三维联合分布概率Fig.10 Three-dimensional joint distribution probability diagram of precipitation,flow and NH3-N

可以看出，COD 与NH3-N 的三维联合概率分布规律并不一致。当联合概率值小于0.6 时，随着P、Q、水质指标质量浓度的增大，三者联合概率值增大，其在XY 轴的投影形态比较类似，且主要沿着流量增大的方向变化。当联合概率值大于0.6 时，COD 与NH3-N 的三维联合概率的演变趋势发生变化，其对应的流量基本在200 m3/s 左右，但对应的降水区间有较大差异。

3 讨论

郝守宁等[24]研究显示，水质随河流方向有改善的趋势，大汶河流域泰安段的水质也呈现了相同的变化规律，但不同河流水质的沿程变化特征会受到沿岸城市分布、人类活动的剧烈程度影响，陈善荣等[25]对长江干流59 个水质监测断面进行水质变化特征研究，结果显示上游水质好于中下游。

蔡帅[26]研究表明，受来水量影响，大连市蔡房身大桥断面水质季节变化较为明显，丰水期（8月）水质优于平水期（10月）和枯水期（4月）水质。这与大汶河流域泰安段7 个水质站点的汛期与非汛期水质变化规律一致。冯卫等[2]研究表明，受人为因素影响，河流沿岸点源的污水排放在丰水期、平水期和枯水期存在较大差异，会导致水质变化的时间差异。

综上可以看出，河流水质的时空变化规律受自然因素、人为因素综合影响，在不同的流域会呈现不同的变化特征。

本文分别采用Frank Copula、Clayton Copula 函数构建了大汶口站水量水质的二维、三维联合分布，而张翔等[13]对淮河流域蚌埠闸的水量水质的联合分布研究采用的是Gumbel-Hougard Copula 函数，主要是因为Copula 函数本身不具有水文物理基础，不能从物理意义上确定各随机变量的理论联合分布模型，通常通过拟合优度检验，选择与实测资料配合最优的Copula 函数[27]。刘章君等[27]综述了Copula 函数在水文领域的研究进展，证实Copula 函数是一种灵活构造多变量联合分布和处理多变量问题的有效工具，具有良好的适用性和不可替代的优越性。对于水量水质的研究，不仅用于构建二者的联合分布模型，在水量水质联合风险的综合评价[28]中也体现出较大的潜力。本文对于水量水质关系的分析主要采用大汶河干流大汶口水文站的实测资料，受资料限值，缺少对其上游、下游及支流的水量水质关系分析，需要在未来的研究中做进一步分析讨论。在构建水量与水质二维、三维联合分布的基础上，应继续开展应用研究，开展水量水质风险分析，为水污染防治、水量水质综合监测管理提供参考。

4 结论

1）在时间上与空间上，大汶河流域泰安段COD与NH3-N 的质量浓度均呈下降趋势，其中楼德站为水污染防治的重点治理断面。

2）季节性Kendall 检验方法结果表明，71%的评价结果为高度显著下降趋势。

3）大汶口站实测水量水质相关分析表明，随着降水量与流量的增大，COD 与NH3-N 质量浓度下降。非汛期到汛期，COD 与NH3-N 的质量浓度分别下降6%、37%，Ⅳ类水达标月份占比提高10%。

4）通过拟合优度检验，选择Frank Copula 函数构建大汶口站的水量与水质的二维联合分布，COD与NH3-N 的概率分布规律基本一致。选择Clayton Copula 函数构建水量与水质的三维联合分布，当联合概率值大于0.6 时，COD 与NH3-N 的概率分布规律发生变化。