3.2 水文风险分析

3.2.1 水文风险分析基本程序

水文风险分析的程序与风险管理框架中风险识别与评估的步骤基本一致，即确定风险目标、划定目标所处的时空范围、选取风险分析方法、准备风险分析数据资料、计算目标风险概率、估算目标风险影响等，如图3.1所示。

第一步，明确水文风险分析的目标。在进行水文风险分析之前，首先要确定风险分析的对象及风险预期目标，即根据研究或管理需要，明确分析对象是供水风险、需水风险，还是暴雨洪水风险，或者其中的一种或多种风险等，初步判断分析对象可能达到的风险预期目标，多个风险共同研究时还需要考虑目标之间是否有直接或者重要的影响。确定风险目标后，还需要进一步考虑实现目标的具体任务以及技术路线等。

第二步，划定水文风险分析的时空范围。根据研究或管理需求，结合研究区域的风险分析目标和相关资料，界定分析所限定的时间和空间范围，初步判断在该时空范围内目标风险可能在什么地方发生、什么时间发生、发生的概率有多大、其危害有多大范围、持续的时间等，并且需要对一些风险较大的区域和时段进行重点分析。例如，对于供水水文风险，从狭义的角度讲，其风险发生的区域是供水的水源地，可据此查找并列出研究对象区域范围内的水源地范围，即为供水水文风险的空间范围；时间范围则根据其监测资料的系列长短和管理与研究的实际需求确定。

第三步，选取水文风险分析方法。风险分析方法有很多，如贝叶斯网络分析方法、蒙特卡洛法、灰色随机风险分析法、最大熵风险分析法等，每一种方法都有其优缺点及适用范围，因此必须根据研究对象及目标选择适宜的方法。

第四步，准备水文风险分析所需的数据资料。在确定风险分析方法后，根据方法所需要的数据资料及风险目标，在研究区时空范围内搜集与风险分析指标有关的数据信息，如气象信息、水文信息、地形地貌信息、社会经济发展信息等。当同一个风险指标可以选择一个或多个风险表征信息时，需要考虑不同信息之间的关联度，在同等相关程度条件下，要优先考虑能够直接测量获取或资料较为丰富的指标数据，如降水、气温、蒸散发、社会经济用水等。对于非定量的数据信息，还需要采取层次分析法、模糊数学方法等途径，借助于专家经验知识将定性的指标信息按照一定的标准或规则进行定量化描述。例如，降水信息可以采用气象站监测得到的降水系列值，而供水保证率指标在没有详细数据资料的情况下也可以根据熟悉情况的专家进行经验判定。

第五步，计算目标水文风险的概率。在充分的资料收集和准备基础上，采用选定的风险分析计算方法，计算和统计目标风险的频率、强度、持续时间和发生范围等风险指标值。

第六步，估算目标水文风险的影响。根据估算的水文风险概率、强度以及时空变化特征等信息，结合实际调查资料估算该水文风险可能造成的人员、财产损失量及其影响的时空范围。

3.2.2 水文风险变量筛选

风险变量是表征风险大小或严重程度的指示器。水文风险分析中，在明确风险目标后就需要对表征该风险目标的变量进行筛选，具体方法可参考水文风险因子识别一章介绍的定性或定量方法。从通过风险识别确定的主要因素选取适宜的指标，指标的确定要具有代表性和针对性，即要求能够准确反映该因素的状态变化情况，并与风险控制目标相呼应，能够通过该指标描述或解决风险目标对应的问题。一项因素可以选取一个或多个代表性指标，但是在具体操作过程中要注意避免与其他相近指标的重复选择，还要兼顾指标量化的难易程度，以供水水文风险分析为例，其目标是供水水文风险的大小，影响因素包括气象、水文、人工用水、调水工程引水等多个方面，其中气象因素涉及降水、气温、蒸散发、风速、辐射等多个指标，通过敏感性分析及风险识别，可以确定降水、蒸散发是两个最具有代表性和针对性的指标；同时，水文指标一般可以用径流过程、水库蓄水量、水库蓄水水位等指标表征，但是考虑到径流过程较易获取，可以优先考虑选取。最后对选取的指标进行定量化描述，收集整理相关的数据资料，为下一步分析做好准备。

3.2.3 水文风险分析计算方法

比较常见的水文风险分析方法有贝叶斯网络分析法、蒙特卡洛法、灰色随机风险分析法、最大熵风险分析法等。其中贝叶斯网络分析法既可以用于识别主要风险因子，也可以对风险的概率进行定量计算，其具体的计算原理在第二章已有详细介绍，这里不再赘述，重点介绍一下其他几种风险分析方法的计算原理和计算公式。

3.2.3.1 蒙特卡洛法

蒙特卡洛法，又称统计试验法或随机抽样法。以概率统计理论为基础，以随机抽样为手段，通过预先制订各影响因素的操作规则和变化模式，人工产生各指标因子的随机数，从大量的统计分析计算结果中找出指标因子的概率分布特征，进而求出相应指标因子的风险率。

蒙特卡洛法的基本思想早在17世纪就被人们发现，并利用事件发生的“频率”来决定事件的“概率”。例如19世纪出现的蒲丰氏投针试验，就是将长为2L的一根针任意投到地面上，用针与一组相间距离为2a（L&lt；a）的平行线相交的频率代替概率P，然后利用P与L和a的关系式推求圆周率π的近似值。为了得到具有一定精度的近似解，需要进行多次试验，由于人工进行大量的试验较为困难，甚至是无法做到的。因此，蒙特卡洛法的基本思想虽然早已提出，但是一直没有得到广泛应用。直到上世纪40年代以来，随着电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验通过计算机来完成，从而给蒙特卡洛法的应用提供了广阔的空间，在各个行业都发挥了重要的作用。

1.基本原理

基本原理：假定函数Y=f（X₁，X₂，…，X_n），其中变量X₁，X₂，…，X_n的概率分布或发生频率已知。

蒙特卡洛法利用一个随机数发生器通过直接或间接抽样取出每一组随机变量X₁，X₂，…，X_n分布的一组值X_1i，X_2i，…，X_ni，然后按Y对于X₁，X₂，…，X_n的关系式确定函数Y的值y_i：

按照上述方法反复独立抽样（模拟）多次（i=1，2，…，n），便可得到函数Y的一批抽样数值y₁，y₂，…，y_n，当模拟次数足够多时，便可以给出与实际情况相近的函数Y的概率分布及其特征数字。

2.收敛性

由蒙特卡洛法的基本原理可知，该方法是由随机变量X的简单子样X₁，X₂，…，X_n算术平均求得：

作为所求解的近似值，由大数定律可知，若X₁，X₂，…，X_n满足独立同分布，且具有有限期望值E（X）&lt；∞，则：

即随机变量X的简单子样的算术平均值，当子样数N充分大时，以概率1收敛于他的期望值E（X）。

3.误差

蒙特卡洛法的近似值与真值的误差问题，可以通过概率论的中心极限定理找出答案。根据该定理，如果随机变量序列X₁，X₂，…，X_n满足独立同分布，且具有有限非零的方差σ²，即：

其中，f（x）是X的概率分布密度函数，则有：

当N充分大时，有如下近似式：

其中，α称为置信度，1-α为置信水平。

由式（3.6）可知，不等式近似地以概率1-α成立，且误差收敛速度的阶为O（N^-1/2）。那么，可将蒙特卡洛法的误差ε定义为：

其中，λ_α与置信度α是一一对应的，根据问题的要求确定出置信水平后，查阅标准正态分布表，就可以确定λ_α值。

需要指出的是，蒙特卡洛法的误差为概率误差，与数值计算方法中的误差概念是有区别的；其次，误差中的均方差σ是未知的，需要采用其估计值代替，σ的估计值计算如下：

给定置信度α，误差ε由均方差σ和试验次数N决定，减小方差σ²或增大试验次数N；若均方差σ固定，把精度提高一个数量级，试验次数N则需要提高两个数量级；相反，如果能减小估计的均方差σ，比如降低一半，则误差ε也减小一半，相当于试验次数增大4倍的效果。

4.优缺点

首先，采用蒙特卡洛法进行风险分析，能够较为逼真地描述具有随机性的事物特征及物理过程，因而该方法能够部分替代物理实验，甚至可以得到物理实验难以得到的结果。其次，蒙特卡洛法受几何条件限制较小，在具有随机性质的问题中，如果考虑的系统几何条件很复杂，采用一般的数值方法难以得到满意的结果，使用蒙特卡洛法一般不会有原则上的困难。第三，蒙特卡洛法的收敛速度与问题的维数无关。前面在误差的定义中提到，在给定的置信水平情况下，蒙特卡洛法的收敛速度为O（N^-1/2），而问题维数的变化只能引起抽样时间及估计量计算时间的变化，并不会影响误差值，即使用蒙特卡洛法时，抽取的子样总数N与维数S无关，维数的增加，除了增加相应的计算量外，不影响问题的误差，这也决定了蒙特卡洛法对多维问题的适应性。第四，蒙特卡洛法具有同时计算多个方案与多个未知量的能力。当遇到多个方案的问题时，一般方法需要逐个计算，而蒙特卡洛法则可以同时计算所有方案，其计算量几乎与计算一个方案的计算量相当。除此之外，蒙特卡罗法还可以同时得到若干个所求量，而不是像常规方法那样逐一计算所求量。第五，蒙特卡洛法的误差容易确定。对于常规计算方法，得到计算结果与真实值的误差并不容易，而蒙特卡洛法可以根据误差公式，在计算所求量的同时计算出误差值，并且该方法还不存在有效位数损失的问题。第六，程序结构简单，利用计算机实现蒙特卡洛法较为简单。

但蒙特卡洛法也存在一些不足。首先，收敛速度较慢，一般不容易得到精确度较高的近似结果，并且对于维数较少的问题（三维以下），其收敛效果也不如其他方法。其次，蒙特卡罗法的误差具有概率性，而不是一般意义下是误差。第三，计算结果与系统大小有关，对于大系统或小概率事件的计算问题上，蒙特卡洛法的计算结果往往低于真实值，而数值方法则比较适用于大系统，因此通常会考虑将蒙特卡洛法与数值方法结合起来，取长补短，解决单纯使用一种方法难以解决的问题。

5.应用

蒙特卡洛法自身的特点，使得它在各行各业的应用越来越广泛。对于水文风险分析来说，由于水文系统涉及的风险因子众多，在运用蒙特卡洛法时可将风险变量进行分解，将分解得到的最终变量作为输入变量，并假设各变量之间是相互独立的，然后对每个变量采用蒙特卡洛法进行分析，求出其风险概率分布。理论上讲，变量分解得越细致，模拟结果的可靠性也就越高。但是分解过细就有可能造成变量之间存在一定的相关性，违背了变量间相互独立这一假设前提。因此，运用此法时应根据实际情况对变量分解的程度以及模拟的次数进行一定的限制。

3.2.3.2 灰色随机风险分析法

由于水文系统的不确定性既有源于系统内在的不确定性，也存在着模型的不确定性、参数的不确定性和获取信息的不足和不精确性，加之对水文系统的认识水平有限，总是会出现水文系统部分信息已知、部分未知的情况，此时可将该系统称为灰色系统，将系统中的变量视为灰变量，应用灰色区间预测方法来度量系统的不确定性。因此，简单的说，水文风险就是指水文系统遭受的外来负荷大于系统本身某一水平的阻抗或承载能力的概率，而灰色随机风险分析方法是在灰色随机风险率的方法基础上，对水文风险率的灰色不确定性进行描述和量化。

1.灰色概率分布理论基础

为了搞清楚灰色随机风险分析方法的基本思路，这里先了解一下灰色概率分布的一些基本概念和理论基础，如灰色概率、灰色概率分布、灰色概率密度、灰色期望和灰色方差的定义及数学表达方式（胡国华，2001）。

（1）灰色概率：称映射P_G：Ψ→P（［0，1］）为灰色概率，若P_G是（Ω，Ψ）上的闭区间集值测度，且1∈P_G（Ω），（Ω，Ψ，P_G）称为灰色概率空间。

（2）灰色概率分布：在样本空间Ω上，取值于实数域的函数ξ，称为是样本空间Ω上的随机变量，并称F_G（x）是随机变量ξ的灰色概率分布函数，简称灰色概率分布，记为：

式中：P_G为实数域上随机变量ξ的灰色概率估计。

（3）灰色概率密度：若ξ为随机变量，F_G（x）为它的灰色概率分布，如果存在函数f_G（x），使对于任意的x，有：

则称f_G（x）为F_G（x）的灰色概率密度函数，简称灰色概率密度。为了更为直观地表现灰色概率分布的灰色特性，还可将式（3.10）表达为：

（4）灰色期望：在空间（Ω，Ψ，P_G）中，若ξ是一个随机变量，其灰色概率分布为F_G（x），当满足

时，则称ξ的数学期望（简称灰色期望）存在，且

（5）灰色方差：在空间（Ω，Ψ，P_G）中，若ξ是一个随机变量，其灰色概率分布为F_G（x），称ξ的灰色方差存在，且

2.灰色随机风险率的表达形式

水文系统的风险可以解读为系统不能够实现原定功能的概率，或者系统外部负荷大于系统本身的阻抗或承灾能力的概率。对此，胡国华等（2001）分三种情况对灰色随机风险进行定义和表达，即单一阻抗与单一负荷失效情况、系统阻抗与外来负荷失效情况、系统极限状态情况，下面具体介绍对上述三种情况的定义与数学表达形式。

（1）系统单一阻抗与单一负荷失效情况。在系统风险率由一个系统阻抗与一个系统外来负荷构成情况下，当考虑系统阻抗或承载能力为一个遵从某一经典分布的随机变量（用X表示），而系统外来负荷为一个遵从某一灰色概率分布的随机变量（用Y_G表示）时，则系统失效的风险率定义为：

式中：P_G为灰色概率；R_G为灰色-随机风险率。

如果能求得X的概率密度函数f_X（x）和Y_G的灰色概率分布函数或灰色概率密度函数时，式（3.15）还可以写成：

那么，失效风险率R_G的质的量度就可以采用曲线f_X（x）分别与与的重叠部分来表示。

（2）系统阻抗与外来负荷失效情况。对于水文系统来说，若系统阻抗和系统外来负荷为分别遵从某一灰色概率分布F_G（x）与F_G（y）的随机变量X_G、Y_G，则系统失效的风险率可表示为：

如果能求得X_G、Y_G的灰色概率分布函数、或者灰色概率密度函数、时，则曲线、（x）分别与与的重叠部分即为失效风险率R_G的质的量度。

当考虑R_G为描述失效风险率的最小与最大值区间时，式（3.17）还可以进一步展开如下：

（3）系统极限状态情况。水文系统由许多因素或变量构成，假定采用向量X_G=（X_G1，X_G2，…，X_Gn）来表示水文系统中的这些因素或变量，用g（X_G）表示水文系统性能函数，则系统失效导致水文风险发生的概率即为g（X_G）&lt；0的概率。

假设系统的状态变量X_G服从用灰色概率分布函数或灰色概率密度函数表征的某一灰色概率分布，则可定义水文系统的临界性能要求为g（X_G）=0，即系统的“极限状态”。相应地，可知水文系统的“安全状态”和“失效状态”可分别表示为：

根据几何学原理，极限状态方程g（X_G）=0是一个n维的表面，称之为“失效面”。显然，由于状态变量X_G的灰色特性，存在2ⁿ个失效面。若失效面到原点的最小距离用向量d=（d₁，d₂，…，）表示，则由于失效面的位置可用失效面到原点的最小距离来表示，而相对于原点的失效面的位置又决定系统的风险性。因此，所有最小距离的最小值min（d）与最大值max（d）可近似地用于风险性的度量，即［min（d），max（d）］。因此，上述变量X_G1，X_G2，…，X_Gn的联合概率密度函数为（X_G1，X_G2，…，X_Gn），失效状态的概率为：

简写为

3.灰色随机风险率的求解方法

由于灰色-随机风险率可以分解为2ⁿ个随机风险率的形式来表达，因此可将灰色-随机风险率转换成一般的随机风险率，进而采用其他方法进行计算，如改进一阶二矩法。原则上，R_G应从向量R_G=（R₁，R₂，…，）中通过取最大、最小值得到，即

在实际运用中，容易通过分析相关信息得知各因素x_G变化对R_G的影响。因此，能够构造各因素的一、二阶矩取值的两两组合，并从系列组合中求出最小和最大风险率值，分别记为R_min、R_max，即

然后采用改进的一阶二矩法进行求解得到结果。

3.2.3.3 最大熵风险分析法

熵是一种表征系统动力学中不能做功的能量的总和，是经典热力学中的一个重要物理概念，最早由德国物理学家克劳修斯于1865年提出。随后，熵的概念逐渐推广至非热力学领域，在信息学、概率论、天体物理、生态科学、生命科学等领域都有着广泛的应用。1957年，Jaynes基于熵理论提出了最大熵原理（Principle of Maximum Entropy，简称POME），并以之成功地解决了信息科学中大量存在的不确定性问题，由此开创了最大熵原理研究的先河。最大熵原理认为，在所有的可行解中，应该选择其熵最大的一个。熵值最大意味着对因为信息数据不足的人为假定（人为添加信息）最小，以此获取的解偏差最小。

在水文风险分析中，最大熵风险分析方法中的“熵”即信息论中所定的“信息熵”，是指水文系统所面临的各种不确定性信息的度量。由于实测资料、所得信息的限制，往往很难通过统计分析的途径确定所有水文风险因子的后验概率分布，而不得不根据所掌握的一些粗略或模糊的认识，如分布的上下限、期望值、分布的形状等约束条件，对各个水文风险因子的先验概率分布做出人为判断或假定。这样所带来的问题是：在推断各个水文风险因子时，很有可能由于不适当地增添了人为额外信息而使得风险分析成果的不确定性增大，影响了水文风险分析的准确性。而最大熵原理将信息熵函数视为风险因子先验概率分布的泛函，由给定的不同约束条件（已知信息）来确定风险因子的先验概率分布。可将信息分为离散信源和连续信源两种情况，以风险分析中几种常用的概率分布和参数估计为例，就能够从数学角度解释了如何在不同已知信息下，由最大熵原理确定风险因子相应的满足约束的最小偏差先验概率分布。或者从某种意义上说，这样所确定的分布是无参数，或者说参数是已知的，即在确定出分布的同时，已经得到了该分布的有关参数，这就避免了上述缺陷。

最大熵风险分析的依据是风险变量的概率特性，其实质是将问题转化为信息处理和寻优问题。在水文系统中，许多风险因子的随机特性都无先验样本而只能获得一些数学特征。最大熵风险分析方法就是利用最大熵准则从这些数学特征的无穷多个概率分布中优选出最合理的分布。这正是水文风险分析中的关键内容之一。

1.最大熵模型

设水文风险随机变量指标x的概率密度分布函数为f（x），则其熵定义为：

相应的最大熵原理表示为：

上诸式中，x=G（x₁，x₂，…，x_n）；Ω为水文风险随机变量指标x的所在集合；m为x的原点矩的阶数；M_i为随机变量x的第i阶原点矩；f（x）为x的概率密度分布函数；n为影响x的风险因子数；b是保证x有意义的变量。

上述数学模型表示在满足已知信息约束下，以熵最大为准则求得水文风险随机变量的概率密度函数。

2.模型求解方法

最大熵模型的求解本质上是一个泛函条件极值问题。为得到f（x）的表达式，可以根据变分法引入拉格朗日乘子λ₀，λ₁，λ₂，…，λ_m，然后令

则问题就转化为使得L这个泛函达到极值。由变分法，得到相应的欧拉方程为：

求解得到：

这就是最大熵概率密度函数的解析形式，只要确定其中的参数λ₀，λ₁，λ₂，…，λ_m，就可以完全确定f（x）。利用式（3.26）、式（3.27），通过数学推导，得到求解λ₀，λ₁，λ₂，…，λ_m应该满足的联立方程组：

将式（3.32）中的方程组看成是有关λ₀，λ₁，λ₂，…，λ_m的m个方程。为了便于数值求解，可将其改写为：

为了便于建立优化函数表达式，令：

上两式中：r_i为残差，可用数值计算方法使其趋于零。

可利用非线性规划等方法求解R的最小值。当R&lt；ε（ε为规定的允许误差）或所有的|r_i｜&lt；ε时即认为式（3.33）收敛，从而求解出λ₁，λ₂，…，λ_m，再按照式（3.32）求得λ₀。将全部的待估计参数代入式（3.31），即可求出风险因子的概率密度函数的解析表达式。