CN106960108A - 基于贝叶斯网络的水库上游来水压力分析方法 - Google Patents

Abstract

本发明提供了一种基于贝叶斯网络的水库上游来水压力分析方法，包括：步骤1、建立压力源特定时刻污染物浓度与受体特定时刻污染物浓度之间的函数关系；步骤2、确定受体数目，并建立贝叶斯层次模型；步骤3、对所述贝叶斯层次模型进行参数估计，获得参数估计值；步骤4、基于所述参数估计值，计算各个压力源对有效浓度的贡献比例；步骤5、设定受体污染物的浓度阈值，以超出阈值的百分比表征压力源对受体的压力，并划分不同的预警等级。该方法建立起了上游多个压力源对下游受体的污染物浓度之间的关系，对河流污染物的判断、预警、分析等，提供了可靠的数据依据及方法。

Description

基于贝叶斯网络的水库上游来水压力分析方法

技术领域

本发明涉及一种水体环境监测分析方法，属于水体环境监测和治理领域，尤其涉及一种基于贝叶斯网络的水库水体压力源分析方法。

背景技术

大坝在水力发电、防洪、航行等方面发挥着重要的作用，然而，一些生态和环境问题也随之而来，例如水坝库区环境的富营养化和藻类水华问题，以及可能的来自上游水源的生活污水、轮船废弃物、工业污水污染等问题。

为防治库区水环境问题，尤其是水体富营养化问题，需要对水体各个阶段链条进行串联研究，现有技术中，很多研究者进行了研究，以三峡库区为例。

(1)对于上游来水，集中在对上游断面水质时空变化特征的分析。例如，(曹承进etal.,2008；郑丙辉et al.,2008)以2003-2004年的监测数据为基础，分别探究了库区上游3个主要污染源(长江、嘉陵江、乌江)总磷(TP)、总氮(TN)浓度的变化特征；(Han et al.,2016)分析了长江、嘉陵江、乌江2004-2013年的TP浓度数据，探究了这3个污染源的TP输出特征，发现乌江断面TP浓度发生了激增，并指出在高水位时期大坝导致回水会携带大量可溶性磷进行支流，造成支流可溶性磷的累积并为支流春季水华提供充足营养盐；(Ren etal.,2015)扩展了研究范围，分析了2003-2010年沱江、长江、乌江、嘉陵江、岷江等5个主要污染源TN浓度的时空特征；(Ren et al.,2016)分析了岷江近10年来的TP变化趋势，指出流域污染防治应该点源和非点源并重。

(2)对于库区水质，主要是比较水库蓄水前后干、支流水质的变化情况以及水质评价方法的研究。例如，(郑丙辉et al.,2006)根据模型计算得到富营养化敏感指数，将库区水体分为河流型、过渡型和湖泊型3种，并确定了过渡型和湖泊型水体的营养指标分级标准值，提出了营养化水平的综合评价方法；(刘辉et al.,2010)比较了三峡库区不同蓄水阶段干、支流的水质状况，发现库区干流现状水质较好且稳定，支流水质有所下降且蓄水后部分支流的局部区段发生了水华；(吕怡兵et al.,2007)根据2003年5、6、10月的监测数据对水库水质进行了评价；(胡莲et al.,2013)分析了三峡水库试验性蓄水后小江流域的水质指标的时空变化规律。

(3)由于回水区水华现象频发，干流倒灌和支流中上游对回水区的影响备受研究者关注。例如，(操满et al.,2015)研究了干流倒灌对梅溪河营养盐浓度的影响，发现干流倒灌的影响不仅局限于河湾区域，对于河流中上游亦有影响，干流倒灌为藻类爆发提供了有利条件；(Yang et al.,2010)的研究指出，干流营养盐倒灌会增加支流水体富营养化的风险；(Ye et al.,2007)分析了香溪河湾春季水华的空间分布及限制性营养盐，发现硅是其主要限制性营养盐；(李凤清et al.,2008)的研究表明，香溪河对香溪河湾TN、TP的贡献率分别为68.5％和91.7％；(刘德富et al.,2016)对已有的三峡库区水动力学特征及其对水环境的影响和支流水华发生机理研究进行了总结，提出了对以后研究的期望；(宋林旭etal.,2016)采用SWAT模型对香溪河氮、磷的非点源分布规律进行了探索，指出支流高岚河流域和古夫流域应当作为香溪河流域污染控制的重点。(张磊et al.,2015)对小江流域的研究表明，在水库调控的高水位时期(9月-次年4月)，干流倒灌是回水区高阳平湖监测断面硝氮和溶解性磷的主要来源。

在以控制库区富营养化和水华问题为目的进行的“上游来水→库区干流水质→支流水质”研究中，“上游来水→库区干流水质”的研究尚未引起足够重视，现有技术以及学术领域也未提出有效的水质分析方法。尽管从水质状态的角度来看，库区干流水质较好且稳定，但是由于干流倒灌对回水区营养盐具有重要的补给作用，而上游来水是库区干流污染物的重要来源，因此有必要建立上游来水对库区干流水质的响应关系，探究上游来水对干流水质的影响。

发明内容

有鉴于此，本发明实施例提供一种基于贝叶斯网络的水库上游来水压力分析方法，具体而言，本发明提供了以下的具体技术方案：

一种基于贝叶斯网络的水库上游来水压力分析方法，该方法包括：

步骤1、建立压力源特定时刻污染物浓度与受体特定时刻污染物浓度之间的函数关系，所述函数关系如下：

C_i+t＝f(C_i1)

C_i+t＝g(C_i2)

其中，t表示污染物由压力源到受体所需的时间，C_i1表示压力源某种污染物i时刻的浓度，C_i+t和C_i2表示受体在i+t和i时刻该污染物浓度；

步骤2、确定受体数目，并建立贝叶斯层次模型；

步骤3、对所述贝叶斯层次模型进行参数估计，获得参数估计值；

步骤4、基于所述参数估计值，计算各个压力源对有效浓度的贡献比例。

步骤5、设定上游压力源的不同取值，通过贝叶斯网络，求取不同受体污染物的概率分布，并设定受体污染物的浓度阈值，以超出阈值的百分比表征压力源对受体的压力，并划分不同的预警等级。

优选地，所述步骤2中的贝叶斯层次模型具体如下：

TN_ij＝a_j×TN_L_ij+b_j+ε_TN，ij

TP_ij＝c_j×TP_L_ij+d_j+ε_TP，ij

TN_L_ij＝α_j×CJ_TN_ij+β_j×JLJ_TN_ij+γ_j×WJ_TN_ij

TP_L_ij＝α_j×CJ_TP_ij+β_j×JLJ_TP_ij+γ_j×WJ_TP_ij

m_a，m_b，m_c，m_d～N(0，10)；σ_a，σ_b，σ_c，σ_d～U(0，10)

α～U(0，1000)；β～U(0，1000)；γ～U(0，1000)

j＝1，2，3，4，...

其中，j表示受体的序号，i表示某一特定受体监测值的序号，TN和TP分别表示受体浓度，CJ、JLJ和WJ分别表示压力源1、压力源2和压力源3三个压力源，TN_L和TP_L表示三个压力源对受体回归的“有效浓度”；参数a、c和b、d分别表示有效浓度对受体浓度的回归斜率和截距，这些参数服从同一个未知分布，α、β和γ分别表示压力源1、压力源2和压力源3对有效浓度的回归系数，和表示受体浓度的回归残差；m_a、m_b、m_c、m_d和σ_a、σ_b、σ_c、σ_d是层次模型的超参数；N(m，σ²)和U(a，b)分别表示正态分布和均匀分布。

优选地，所述步骤2中的贝叶斯层次模型，对于仅有2个压力源的目标对象，其模型建立具体如下：

TN_i＝a×(α×CJ_TN_i+β×JLJ_TN_i)+b+ε_TN，i

TP_i＝c×(α×CJ_TP_i+β×JLJ_TP_i)+d+ε_TP，i

其中，i表示某一特定受体监测值的序号，TN和TP分别表示受体浓度，CJ、JLJ分别表示压力源1、压力源2两个压力源，TN_L和TP_L表示三个压力源对受体回归的“有效浓度”；ε为残差项；参数a、c和b、d分别表示有效浓度对受体浓度的回归斜率和截距，这些参数服从同一个未知分布，α、β分别表示压力源1、压力源2对有效浓度的回归系数。

优选地，所述步骤4中的有效浓度的贡献比例计算方法为：

其中，分别表示压力源1、压力源2和压力源3中某种污染物的平均值，α、β和γ分别表示压力源1、压力源2和压力源3对有效浓度的回归系数。

优选地，所述步骤3中的参数估计，采用贝叶斯参数估计，具体方法如下：

其中π(θ)和π(θ|y)分别为参数θ的先验分布和后验分布，f(y|θ)表示在特定参数条件下观测值发生的概率，即为似然函数。

优选地，所述步骤5中，设定上游压力源的不同取值，求取不同受体污染物的概率分布，具体方法为：

将自变量浓度分布的5％-95％分位数分为N段，对每一段进行积分，求解平均值；抽样次数为M次，以M的平均值作为因变量的值，所述M、N为正整数。上述的M、N值，可以根据实际需要进行调整，例如，可以设定M＝N＝100。

与现有技术相比，本发明技术方案提供了一种可以有效对水库、水坝中水源上游来水对库区水质影响进行分析的方法，实用性强，并通过对主要污染源和“上游来水→库区干流水质”响应关系的识别，增强对库区上游来水对各个监测断面压力的定量认识，从保护干流水质的角度出发，为上游来水污染防治提供建议；同时完善库区污染防治研究的“上游来水→库区干流水质→支流水质”链条，为进一步构建整个库区的大系统模型提供支持。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本发明实施例的污染物从压力源到受体过程的概念模型；

图2为本发明实施例的研究对象水域实例图；

图3为本发明实施例的寸滩断面的响应关系；

图4为本发明实施例的清溪场断面的响应关系；

图5为本发明实施例的晒网坝断面的响应关系；

图6为本发明实施例的培石断面的响应关系；

图7为本发明实施例的培石断面的响应关系。

具体实施方式

下面结合附图对本发明实施例进行详细描述。应当明确，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本领域技术人员应当知晓，下述具体实施例或具体实施方式，是本发明为进一步解释具体的发明内容而列举的一系列优化的设置方式，而该些设置方式之间均是可以相互结合或者相互关联使用的，除非在本发明明确提出了其中某些或某一具体实施例或实施方式无法与其他的实施例或实施方式进行关联设置或共同使用。同时，下述的具体实施例或实施方式仅作为最优化的设置方式，而不作为限定本发明的保护范围的理解。

研究区域以长江某水段为研究对象，如图2所示，其中寸滩、清溪场、晒网坝、培石和银杏沱为库区的主要监测断面(即为受体)，关注的水质指标为总氮(TN)和总磷(TP)；上游来水主要考虑3个来源，即为以朱沱断面为代表的长江、以北温泉断面为代表的嘉陵江和以麻柳嘴为代表的乌江。收集了2008-2013年的水质监测数据，其中银杏沱断面TN监测数据缺失，故而不分析上游来水TN对银杏沱TN的影响，其它4个库区干流监测站点同时分析TN和TP；由于寸滩断面在乌江断面以上，因此其上游来水压力仅考虑长江和嘉陵江，其它4个库区干流断面的上游来水压力考虑3个来源。

本研究关注的水质指标是TN和TP的浓度，区分水库调控导致的高、低水位对响应关系的影响，同时分析响应关系的沿程特征。本研究的主要内容包括3个部分：

(1)受体主要压力源识别。根据收集到污染源和受体水质监测数据，建立“污染源→受体”的响应模型，根据不同污染源对受体TN、TP浓度的贡献率，筛选不同受体的主要污染源。在建立响应模型时，由于受体均为三峡库区的水质监测站点，因此采用贝叶斯层次模型(Bayesian Hierarchical Model,BHM)对回归系数进行分层，即设定一个全局参数作为各个受体参数的共同先验分布；此外，假设不同压力源的单位污染物浓度对同一受体不同污染物(TN和TP)的相对重要程度相同；压力源对受体的压力需要根据各个压力源对受体的相对重要程度以及污染物的平均浓度确定；由于高、低水位差距悬殊，研究中假设两种水位时期的响应关系完全独立。综上，上游来水压力源对同一受体TN或TP的联系通过不同压力源对受体的相对性体现，受体响应关系的沿程特征通过BHM联系起来，而高、低水位则认为影响巨大而相互独立的。

(2)“主要压力源→受体”响应关系的建立。在识别出主要的压力源后，对于各个受体，分别建立不同水位时期“主要压力源→受体”响应关系，探究水位高、低对响应关系的影响，探究响应关系的沿程特征，并比较TN、TP响应关系沿程特征的异同。“主要压力源→受体”响应关系的建立采用贝叶斯网络(Bayesian Network,BN)方法，建立主要压力源与受体之间的概率响应关系，为下一步进行压力风险分析提供便利。

(3)主要压力源对受体的压力风险分析。风险分析的目的是探究在控制受体营养盐浓度在一定水平时，探究不同的超标率下对应的主要污染源的最大污染物浓度，或者探究不同污染物浓度对应的受体超标率。显然，当设定不同的超标率对应不同的预警等级时，压力风险分析即为上游来水压力对库区受体的预警分析。本研究借助建立的“主要压力源→受体”的BN，通过多次抽样进行风险分析。

以下结合图1，对本发明的方法进行详细的说明。上游来水压力源由上游达到下游需要一定的时间，称为时滞时间，从理论上看，在同一时刻压力源和受体的同一污染物之间不存在直接的因果关系，直接建立压力源污染物浓度(Ci1)与受体污染物浓度(Ci2)之间的统计模型似乎缺乏机理过程的支撑。图1给出了污染物从压力源到受体过程的概念模型(假设只有一个压力源和一个受体)。图中，t表示污染物由压力源到受体所需的时间；左侧实心点表示污染源，Ci1和Ci-t分别表示某种污染物某一时刻(i时刻)和i-t时刻的浓度；右侧实心点表示受体Ci+t和Ci2在i+t和i时刻该污染物浓度。一方面，源在i时刻的浓度Ci1由于受到扩散稀释、分解、沉降和再悬浮等因素的影响，在i时刻到达受体时，其浓度变为Ci+t；显然Ci1与Ci+t之间存在一定的函数关系，不妨设为C_i+t＝f(C_i1)。另一方面，受体污染物浓度在不同时刻存在自相关，不妨设C_i+t＝g(C_i2)。因此，在i时刻压力源污染物浓度Ci1与受体污染物浓度Ci2可以通过t时刻受体污染物浓度Ci+t联系起来；这种关系不是机理过程的因果关系，而是统计学上的相关关系，这种关系符合机理过程，虽然不能表征因果关系，但对于建立源与受体之间的关系却是有效的。

之所以不直接建立C_i+t＝f(C_i1)的响应关系，是出于以下两个方面的考虑：(1)在讨论高、低水位的响应关系时，连续的时间序列被打破，因此难以通过统计模型计算得到高、低水位时的时滞时间t；(2)由于不同月份流速不一致，可能导致不同月份时滞时间数值不一致。因而，直接建立C_i+t＝f(C_i1)的响应关系难度较大。

根据以上理论，建立Ci1与Ci2之间的响应关系的显著与否取决于以下两个方面因素：(1)函数C_i+t＝f(C_i1)的显著性。该影响的显著性又受到两个过程的影响，一是由源到受体过程的干扰(使浓度增加或者减少的因素)，使得源浓度与“剩余浓度”(源浓度在Ci1经过t个时刻之后到达受体时的浓度)之间的关系不再显著，例如，若途中有很大的负荷输入或者负荷大部分被分解，则C_i+t＝f(C_i1)的关系必不显著；二是“剩余浓度”与受体“初始浓度”之间的相对大小，如果“剩余浓度”所占的负荷很小，则不能对混合后的受体浓度产生影响，则C_i+t＝f(C_i1)的关系亦不显著。(2)函数C_i+t＝g(C_i2)的显著性，即经过t个时刻之后受体污染物的自相关关系是否显著，若不显著，则最终建立的Ci1与Ci2之间的响应关系亦不显著。可见，对于本研究受体的沿程特征，当受体距离压力源足够远时使得C_i+t＝f(C_i1)或者C_i+t＝g(C_i2)不显著时，得到的响应关系不显著，此距离即为源对受体影响的“有效距离”(L)。需要指出的是，以上理论是本研究建立模型和对结果进行解释的重要理论基础。

由于本研究未收集到关于模型参数的有效信息，因此模型的先验分布一律采用无信息的先验分布形式，即给定一个很大的方差。参数估计的依据是贝叶斯公式：其中π(θ)和π(θ|y)分别为参数θ的先验分布和后验分布，f(y|θ)表示在特定参数条件下观测值发生的概率，即为似然函数。

本研究采用的贝叶斯层次模型中“分层”的含义是指在参数估计过程中对参数进行分层，由于贝叶斯层次模型的超参数采用模糊扁平先验分布时，其结果与随机效应模型相近，此时贝叶斯层次模型可以用R软件(请对R软件的具体名称或定义进行说明)“lme4”包中的lmer()函数简捷地实现。

本研究中，采用的贝叶斯层次模型形式如下：

TN_ij＝a_j×TN_L_ij+b_j+ε_TN,ij

TP_ij＝c_j×TP_L_ij+d_j+ε_TP,ij

TN_L_ij＝α_j×CJ_TN_ij+β_j×JLJ_TN_ij+γ_j×WJ_TN_ij

TP_L_ij＝α_j×CJ_TP_ij+β_j×JLJ_TP_ij+γ_j×WJ_TP_ij

m_a,m_b,m_c,m_d～N(0,10)；σ_a,σ_b,σ_c,σ_d～U(0,10)

α～U(0,1000)；β～U(0,1000)；γ～U(0,1000)

j＝1,2,3,4

其中，j表示受体的序号(1-4分别表示清溪场、晒网坝、培石和银杏沱4个受体)，i表示某一特定受体监测值的序号，TN和TP分别表示受体浓度，CJ、JLJ和WJ分别表示长江、嘉陵江和乌江3个压力源，TN_L和TP_L表示3个压力源对受体回归的“有效浓度”；参数a、c和b、d分别表示有效浓度对受体浓度的回归斜率和截距，这些参数服从同一个未知分布，α、β和γ分别表示长江、嘉陵江和乌江对有效浓度的回归系数，和表示受体浓度的回归残差；m_a、m_b、m_c、m_d和σ_a、σ_b、σ_c、σ_d是层次模型的超参数；N(m,σ²)和U(a,b)分别表示正态分布和均匀分布；ε为残差项。由于寸滩仅有2个压力源，因此无法加入到上述贝叶斯层次模型中，而对单独处理，方程为：

TN_i＝a×(α×CJ_TN_i+β×JLJ_TN_i)+b+ε_TN,i

TP_i＝c×(α×CJ_TP_i+β×JLJ_TP_i)+d+ε_TP,i

变量和参数含义以及参数的先验分布同上。

贝叶斯层次模型参数后验分布的获得需要很高的计算量，常用的软件包括WinBUGS、JAGS和STAN等，这些软件都可以通过R软件很方便地调用。本研究采用R软件的“rstan”软件包进行贝叶斯(层次)模型的参数估计和推断，总共迭代100000次，后50000次用于推断参数后验分布(取1000次有效样本)，以R-hat≈1.0作为模型收敛的标识。在得到参数估计值后，计算各个压力源对有效浓度的贡献比例，计算方法(以长江为例)为：

其中表示长江某种污染物的平均值。主要污染源的判定准则为：对有效浓度贡献率超过70％的污染源。

本研究采用主观方法确定贝叶斯网络的结构，采用R软件中的“bnlearn”包进行贝叶斯网络参数学习(连续贝叶斯网络)；设定上游压力源的不同取值，探究不同受体污染物的概率分布(将自变量浓度分布的5％-95％分位数分为100段，对每一段进行积分，求解平均值；抽样次数为100次，以100的平均值作为因变量的值)；设定受体污染物的浓度阈值，以超出阈值的百分比表征源对受体的压力，并划分不同的预警等级。

下面结合具体的水域实例，对本发明的技术方案进行进一步详细的解释。

为方便表述，我们将函数C′_t+k＝f(C_t)的显著性称为“源显著性”，将函数C′_t+k＝g(C′_t)的显著性称为“受体显著性”。根据上述的提出的技术方案原理，采用贝叶斯层次模型对5个受体断面TN、TP浓度在高、低水位的主要压力源进行识别，取对受体营养盐回归“有效浓度”的总贡献率超过70％的压力源作为受体的主要污染源。除了依据统计模型的结果外，在最终确定主要压力源时，对于某一特定的营养盐，还需要结合压力源对受体影响的沿程特征进行分析，尤其是对同时受3个压力源影响的清溪场、晒网坝、培石和银杏沱断面：(1)在相同流速条件下，下游受体的源显著性和受体显著性均弱于上游受体，因此上游受体的模型拟合效果要优于下游受体；(2)若某一压力源不是上游受体的显著压力源，则亦不应该称为下游受体的显著压力源；(3)高水位时的源显著性弱于低水位时的源显著性，而受体显著性则相反，在高、低水位时源对受体的影响程度不一定具有特定的规律性。此外，不同的源与受体的距离不同，距离越大，两种显著性越弱，源对受体的影响也就越弱。

1、寸滩断面主要污染源识别

由表1可知，无论在高、低水位，寸滩断面TN浓度的主要压力源均为嘉陵江，且长江TN浓度对寸滩TN浓度没有显著影响；寸滩断面TP浓度的主要压力源为长江，在低水位时嘉陵江的贡献率占29％，而在高水位时嘉陵江影响不显著。根据实测数据可知，长江TN浓度远低于嘉陵江，而TP浓度嘉陵江低于长江，因此寸滩断面TN、TP的主要压力源可由2个压力源的浓度大小值来解释。对于TP浓度，由于嘉陵江TP浓度很低，但仍然在高水位时占29％的贡献率，可以推断主要是由于流速减缓导致源显著性减低造成的。

表1：各个压力源对寸滩断面有效浓度的平均贡献率

(注：表中的0表示回归系数不显著,下同)

2、其它4个受体主要压力源识别

除寸滩外的4个受体的主要压力源如表2所示。对于TN，在高水位时，清溪场的主要压力源为乌江，而晒网坝断面则不存在显著压力源，到了培石断面则嘉陵江被识别为主要压力源。显然，这种现象与源对受体影响的沿程特征是不相符合的，需要在下一步分析时予以特别重视。在低水位时，清溪场和晒网坝的主要压力源均为嘉陵江和乌江，与在高水位时相比，多了一个嘉陵江作为主要压力源，原因应为：低水位流速增加，使得嘉陵江与这两个受体之间的源显著性增强；这也表现了水位对主要压力源的影响，说明了对水位进行高、低水位划分的必要性和合理性。低水位时由培石到银杏沱主要压力源的改变也应在分析时给予特别重视。对于TP，在高水位时，清溪场和晒网坝的主要压力源均为乌江，这与乌江TP浓度的平均值和波动性均为3个压力源的最大值相符合；但是到了培石断面，将长江识别为主要压力源，值得警惕。在低水位时，清溪场和晒网坝均以乌江为主要压力源，而到了培石出现了嘉陵江作为主要压力源，需要在后续分析中特别注意。

表2：受体主要压力源识别结果

总体上，对于TN浓度，乌江和嘉陵江由于浓度较高，而被识别为受体的主要压力源；对于TP浓度，乌江由于浓度较高且距离各个压力源最近，而被识别为主要压力源。仅从对主要压力源的识别来看，水位对TN影响更大。

根据受体压力源识别的结果，对不同受体的不同营养盐，分别构建在高、低水位时的贝叶斯网络，表征“源→受体”的响应关系。通过分析响应关系，可以(1)探究主要压力源对受体的影响是否足够大，(2)比较高、低水位的响应关系是否存在差异，并分析这种差异可能的原因。在识别响应主要压力源时，依据的是压力源对回归有效浓度的贡献比例，贡献比例小的压力源对受体营养盐浓度的影响是微不足道的，而贡献比例大的压力源对受体浓度的影响也不一定大。因而，进行响应关系的分析，一方面是对压力源识别的必要补充，另一方面是下一步压力风险分析的基础。根据贝叶斯网络得到的响应关系如下。

(1)寸滩断面响应关系分析

图3左侧图为高、低水位时，长江TP浓度对寸滩TP浓度的响应关系(横坐标为长江TP浓度，纵坐标为寸滩TP浓度)；右侧图为高、低水位时，嘉陵江TN浓度对寸滩TN浓度的响应关系(横坐标为嘉陵江TN浓度，纵坐标为寸滩TN浓度)。对于TP，高、低水位对于响应关系的影响较大，对于TN影响较小。总体而言，在高水位时，同样源浓度对应的受体浓度较高。

(2)清溪场响应关系分析

结合图4，左侧图为在高、低水位时，乌江TP浓度对清溪场TP浓度的响应关系；右侧图黑色实线为在高水位时，乌江TN浓度对清溪场TN浓度的响应关系，较长的黑色虚线表示在低水位时，控制嘉陵江TN浓度在其平均值时，乌江TN浓度对清溪场TN浓度的响应关系，较短的黑色虚线表示在低水位时，控制乌江TN浓度在其平均浓度是，嘉陵江TN浓度对清溪场TN浓度的响应关系。总体而言，在高水位时，同样源浓度对应的受体浓度较高。此外，对比寸滩和清溪场断面营养盐浓度可知，由于乌江的营养盐输入使得清溪场断面的营养盐浓度发生较大幅度升高，而清溪场营养盐浓度受乌江营养盐浓度影响的主导作用亦体现出来。相对而言，长江和嘉陵江断面由于距离清溪场断面较远，使得与清溪场断面的源显著性和受体显著性均减弱，因而与受体浓度不具有显著关系；采用本研究所用的方法，难以将长江和嘉陵江对受体的影响定量化，需要进一步进行数据收集，提升研究的尺度(例如以年数据建立响应关系)。按照营养盐浓度沿程的变化特点，长江和嘉陵江对下游受体影响的显著程度不应强于清溪场断面。

(3)晒网坝响应关系分析

结合图5，左侧图为在高、低水位时，乌江TP浓度对晒网坝TP浓度的响应关系；右侧较长的黑色虚线表示在低水位时，控制嘉陵江TN浓度在其平均值时，乌江TN浓度对晒网坝TN浓度的响应关系，较短的黑色虚线表示在低水位时，控制乌江TN浓度在其平均浓度时，嘉陵江TN浓度对晒网坝TN浓度的响应关系，在高水位时没有影响显著的压力源。对于TP，在低水位时同样源浓度对应的受体浓度较高，这主要是由于随着距离增加，高水位相对于低水位的受体显著性减弱导致的。对于TN，虽然根据压力源识别的结果，在低水位时存在主要压力源乌江和嘉陵江，但是这2个压力源对受体的TN浓度影响均很小，这也与高水位时无显著压力源的结果相一致，也就是说到了晒网坝断面，无论是高、低水位，所有压力源对其营养盐浓度波动性的影响均已经不显著。

(4)培石响应关系分析

结合图6，左侧图中，黑色实线代表高水位时长江TP浓度对培石TP浓度的影响，黑色虚线表示低水位时乌江TP浓度对培石TP浓度的影响；右侧图为在高、低水位时嘉陵江TN浓度对培石TN浓度的影响。对于TP，在高水位时，出现了长江作为主要压力源的情况，从其响应关系来看对TP浓度的影响也很大；然而进一步的分析发现，培石的TP浓度高于上游断面清溪场的TP浓度，尤其是在后半段时间，据此推断在2010年前后培石断面有新的且足够对培石TP浓度造成较大影响的新压力源进入。因此，培石断面TP的主要压力源的统计分析结果不能让人信服的，其响应关系的显著性仅仅是由于数据的协同性导致的，而缺乏科学合理的因果关系推理。因此，在实践中，根据现有的数据分析3个压力源对培石营养TP的压力风险也是没有意义的。而对TN的分析结果则进一步印证了我们在3.2.3中的推断，即3个上游压力源的源显著性和受体显著性已经很低，对于培石断面TN的波动性已经不能造成很大的影响。同样，根据现有的数据分析3个压力源对培石营养TN的压力风险也是没有意义的

(5)银杏沱响应关系分析

结合图7，乌江TP浓度与银杏沱TP浓度的响应关系，由于在培石断面TP有其它压力源的输入，因此对于位于培石下游的银杏沱来说，这种显著的响应关系也缺乏因果关系的支撑，不能够用于进一步分析压力影响。

综合以上分析可见：库区监测站点受上游3个源的影响，尽管长江断面占据了绝对优势的负荷比例，但是由于距离受体较远(时滞时间较长)，源显著性和受体显著性均较弱，其对受体浓度的影响并不显著；但这并不说明长江的污染负荷对受体没有影响，本研究认为朱沱较大的负荷主要构成了营养盐TN和TP的平均水平，而波动则主要受到乌江和嘉陵江的影响。虽然乌江流量最小，但是其TP浓度高，且距离受体较近，对受体TP浓度的波动起到主要作用；嘉陵江断面TN浓度较高，但是由于其距离下游受体较远，影响范围较小，在高水位时仅能影响到清溪场，在低水位时可以影响到晒网坝；乌江TP的影响则可以至少到达晒网坝，由于培石有TP负荷的输入，因此不能判别是否对培石和银杏沱具有显著影响。综合以上分析，对寸滩、清溪场和晒网坝进行上游来水压力分析时，培石和银杏沱的意义不大，因此不予分析。

下面，结合上述技术方案的分析，对压力风险进行分析

(1)预警分级

表3：预警等级划分

采用贝叶斯网络分析压力源对受体的压力，贝叶斯网络由于输出结果为概率分布，因此适用于对受体超过某一阈值的概率进行分析。本研究将受体超过某一阈值的概率水平划分为3级预警，如表3，当压力源的浓度使得受体浓度超过阈值的概率为大于60％，但是小于等于75％时，定义为1级预警，以此类推。对于TN，选择所有受体浓度的50％分位数作为阈值，即为1.8mg/L；对于TP，寸滩断面选择0.10mg/L，晒网坝断面选择0.15mg/L。

(2)各地点压力风险分析

表4为寸滩断面压力风险分析的结果。在进行模拟时，设定了TN的最高浓度为4.5mg/L，TP的最高浓度为0.8mg/L，因此表格中“-”的表示需要的压力源浓度超过了上述值。表格中的浓度值为对应的主要压力源在受体达到一定预警等级时需要的浓度值。单位均为mg/L。表5和表6的含义同上。

表4：寸滩断面不同预警等级对应的主要压力源浓度

清溪场的TN浓度在低水位时有两个影响显著的压力源，分别为乌江和嘉陵江，参见表5，对于TN低水位时的压力源值，是控制其中一个压力源浓度值为平均值时得到的。

表5：清溪场断面不同预警等级对应的主要压力源浓度

晒网坝的TP浓度主要压力源为乌江；对于TN浓度，在高水位时没有影响显著的压力源，在低水位时的主要压力源为嘉陵江和乌江，不同预警等级的浓度值也是控制另外一个压力源浓度值为平均值时得到的，如表6。

表6：晒网坝断面不同预警等级对应的主要压力源浓度

结合上述结果可见，同一受体在高、低水位条件下，对应的主要压力源有所差异，即使对于相同的压力源在数值上也存在差异。需要说明的是，同一受体的不同预警等级对应的主要压力源浓度差异不是很大，这与预警等级的划分方法有关。本研究采用的方法是超过某一特定阈值的概率，即超标的风险，该方法实际上是对超标的可能性做一评价，且采取的分位数分别为60％、75％和90％，差距不是特别大。若采用分位数分别为50％、75％和95％，则不同预警等级对应的压力源浓度差距将会变大。此外，还可设定一定的超标概率，设定不同的阈值来确定预警等级。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random AccessMemory，RAM)等。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (6)

1.一种基于贝叶斯网络的水库上游来水压力分析方法，其特征在于，该方法包括：

步骤1、建立压力源特定时刻污染物浓度与受体特定时刻污染物浓度之间的函数关系，所述函数关系如下：

C_i+t＝f(C_i1)

C_i+t＝g(C_i2)

其中，t表示污染物由压力源到受体所需的时间，C_i1表示压力源某种污染物i时刻的浓度，C_i+t和C_i2表示受体在i+t和i时刻该污染物浓度；

步骤2、确定受体数目，并建立贝叶斯层次模型；

步骤3、对所述贝叶斯层次模型进行参数估计，获得参数估计值；

步骤4、基于所述参数估计值，计算各个压力源对有效浓度的贡献比例。

步骤5、设定上游压力源的不同取值，通过贝叶斯网络，求取不同受体污染物的概率分布，并设定受体污染物的浓度阈值，以超出阈值的百分比表征压力源对受体的压力，并划分不同的预警等级。

2.根据权利要求1所述的方法，其特征在于，所述步骤2中的贝叶斯层次模型具体如下：

TN_ij＝a_j×TN_L_ij+b_j+ε_TN，ij

TP_ij＝c_j×TP_L_ij+d_j+ε_TP，ij

TN_L_ij＝α_j×CJ_TN_ij+β_j×JLJ_TN_ij+γ_j×WJ_TN_ij

TP_L_ij＝α_j×CJ_TP_ij+β_j×JLJ_TP_ij+γ_j×WJ_TP_ij

$a_j~N(m_a,{\mathrm{\σ}}_a^2);b_j~N(m_b,{\mathrm{\σ}}_b^2);c_j~N(m_c,{\mathrm{\σ}}_c^2);d_j~N(m_d,{\mathrm{\σ}}_d^2)$

m_a，m_b，m_c，m_d～N(0，10)；σ_a，σ_b，σ_c，σ_d～U(0，10)

${\mathrm{\ϵ}}_{TN,ij}~N(0,{\mathrm{\σ}}_{TN}^2);{\mathrm{\ϵ}}_{TP,ij}~N(0,{\mathrm{\σ}}_{TP}^2)$

α～U(0，1000)，β～U(0，1000)，γ～U(0，1000)

j＝1，2，3，4，...

其中，j表示受体的序号，i表示某一特定受体监测值的序号，TN和TP分别表示受体浓度，CJ、JLJ和WJ分别表示压力源1、压力源2和压力源3三个压力源，TN_L和TP-L表示三个压力源对受体回归的“有效浓度”；参数a、c和b、d分别表示有效浓度对受体浓度的回归斜率和截距，这些参数服从同一个未知分布，α、β和γ分别表示压力源1、压力源2和压力源3对有效浓度的回归系数，和表示受体浓度的回归残差；m_a、m_b、m_c、m_d和σ_a、σ_b、σ_c、σ_d是层次模型的超参数；N(m，σ²)和U(a，b)分别表示正态分布和均匀分布。

3.根据权利要求1所述的方法，其特征在于，所述步骤2中的贝叶斯层次模型，对于仅有2个压力源的目标对象，其模型建立具体如下：

TN_i＝a×(α×CJ_TN_i+β×JLJ_TN_i)+b+ε_TN，i

TP_i＝c×(α×CJ_TP_i+β×JLJ_TP_i)+d+ε_TP，i

其中，i表示某一特定受体监测值的序号，TN和TP分别表示受体浓度，CJ、JLJ分别表示压力源1、压力源2两个压力源，TN_L和TP_L表示三个压力源对受体回归的“有效浓度”；ε为残差项；参数a、c和b、d分别表示有效浓度对受体浓度的回归斜率和截距，这些参数服从同一个未知分布，α、β分别表示压力源1、压力源2对有效浓度的回归系数。

4.根据权利要求1所述的方法，其特征在于，所述步骤4中的有效浓度的贡献比例计算方法为：

$L=\mathrm{\α}\×\stackrel{\‾}{CJ}+\mathrm{\β}\×\stackrel{\‾}{JLJ}+\mathrm{\γ}\×\stackrel{\‾}{WJ}$

$R_{CJ}=\frac{\mathrm{\α}\×\stackrel{\‾}{CJ}}{L}$

其中，分别表示压力源1、压力源2和压力源3中某种污染物的平均值，α、β和γ分别表示压力源1、压力源2和压力源3对有效浓度的回归系数。

5.根据权利要求1所述的方法，其特征在于，所述步骤3中的参数估计，采用贝叶斯参数估计，具体方法如下：

$\mathrm{\π}\left(\mathrm{\θ}\right|y)=\frac{\mathrm{\π}\left(\mathrm{\θ}\right)f\left(y\right|\mathrm{\θ})}{{\∫}_{\mathrm{\θ}}\mathrm{\π}\left(\mathrm{\θ}\right)f\left(y\right|\mathrm{\θ})d\mathrm{\θ}}$

其中π(θ)和π(θ|y)分别为参数θ的先验分布和后验分布，f(y|θ)表示在特定参数条件下观测值发生的概率，即为似然函数。

6.根据权利要求1所述的方法，其特征在于，所述步骤5中，设定上游压力源的不同取值，求取不同受体污染物的概率分布，具体方法为：

将自变量浓度分布的5％-95％分位数分为N段，对每一段进行积分，求解平均值；抽样次数为M次，以M的平均值作为因变量的值，所述M、N为正整数。