(4)
| 束龙仓1，朱元生1，孙庆义2，彭绪民2 |
| (1.河海大学水文水资源及环境学院；2.济宁市水利局) |
摘 要：传统的地下水允许开采量计算大都有是采用确定性模型进行的。而实际上，在其计算中存在着一些不确定性因素，它们直接影响着计算结果的精度，而“不确定性”与“风险”是紧密相联系的。为此，本文以山东省济宁市地下水资源评价为例，在详细介绍了风险分析的基本概念和风险分析方法的基础上，对地下水允许开采量的确定进行了风险分析，所采用的风险分析方法为蒙特卡罗方法，其中的随机数产生采用混合同余法。经风险分析所得结果较传统的水文地质学方法计算的结果，更符合实际情况。为在今后的地下水资源评价工作中，更加合理地、可靠地确定水源地的允许开采量提供了一种切实可行的方法，使得与地下水开采有关决策的失误减小到最低程度。
关键词：地下水；允许开采量；不确定性；蒙特卡罗方法；灵敏度分析；风险分析 地下水是我国城市，尤其是北方城市的重要供水水源。近年来，随着城市规模的扩大、工农业生产的迅速发展和人民生活水平的不断提高，对地下水的需求量日益增加。由于城市的工业过于集中市区，加之长期过量开采地下水，已引起了一系列的环境地质总是如地下水水位的持续下降、地面沉降、地面塌陷、地裂缝、海(咸)水入侵、地下水污染等，这些环境地质问题已严重制约着城市的发展，甚至威胁着城市的存亡。所谓地下水的过量开采是指地下水的实际开采量超过了允许开采量。因此，允许开采量的准确确定是非常重要的，这也是地下水资源评价的主要任务之一。长期以来，地下水资源评价工作大都是采用确定性模型进行的。而实际上，在地下水资源评价中存在着许多不确定性因素，它们直接影响着评价结果的精度。因此，应对地下水允许开采量的确定进行风险分析。
1 风险分析的基本概念 目前有十多种关于风险(Risk)的定义，不可能取得完全统一，它会随着具体研究问题的不同而有一定的差别。简单地说，风险就是系统失效的概率，即工程达不到预期目标的概率。当提出“什么是风险”时?实际上是提出了三个问题，即：“在这项工程中或在系统运行中会出现什么事故?”，“发生事故的可能性有多大?”或“如果发生了，其后果会怎样?”。系统运行中会出现什么事故，可用事故树Si(即事故的集合)来表示，描绘系统会出现什么事故?是如何出现的?与地下水开采有关的主要事故或主要环境问题有：区域地下水水位的持续下降、地下水污染(或水质恶化)、地面沉降(塌陷)、泉水断流、海水入侵等。出现事故的可能性用Li表示；出现事故的损失或后果用Xi表示，包括了各种可能出现的损失的总和，如：社会的、环境的、经济的等。因此可用三位体的全集{<Si,Li,Xi>}c来定量地定义风险R：
R={<Si,Li,Xi>}c |
(1) |
而风险分析(Risk analysis)则相应包括三个相互联系的部分，即：风险辩识、风险估算和风险评价。风险辩识即为分析系统在特定的决策条件下发生什么样的事故，对地下水供水工程来说，就是分析在某一决策条件下运行供水工程会产生什么环境问题；风险估算则着眼于定量地描述事故发生的概率，即计算出事故发生的可能性有多大；风险评价的主要任务是阐述事故将会产生什么样的社会、经济和环境负效应，即对事故所产生的后果进行评价。
2 风险分析的研究方法 目前，用于风险及不确定性分析的方法很多，如：重现期法、直接积分法、蒙特卡罗法、均值一次二阶距法、加强一次二阶距法和二阶距法等。本文所采用的风险分析方法是目前被广泛应用于科研和生产实际工作中的蒙特卡罗方法，该方法的优点在于无需以显式给出风险因子与研究变量之间的关系，而是通过从实测资料的统计分析中，求得其经验关系[1]。蒙特卡罗(Monte-Carlo，简称为M-C)方法，又称统计试验方法，是人工产生和利用随机数方法的总称。它是一类通过对有关的随机变量或随机过程的随机抽样，来求解数学、物理和工程技术问题近似解的数值方法。具体来说，就是对所要求解的问题，构造一种随机变量或随机过程，使其某一数值特征(如数学期望)为所求问题的解，然后对所构造的随机变量或过程进行抽样，并由得到的样本算出相应的参数值，作为所求问题的近似解。
用蒙特卡罗方法模拟产生满足一定分布的水文地质参数，首先要用某种特定的方法产生该种分布的随机数，这一过程称为随机抽样。随机变量的分布有多种(如：均匀分布、正态分布、皮尔逊Ⅲ型、F-分布等)，不同分布对应的随机数序列也不同。但就随机数的产生而言，最基本的随机变量是在区间[0，1]上服从均匀分布的随机变量。服从其它分布随机变量的随机数可由[0，1]上均匀分布的随机数变换产生。
产生随机数的方法大致可分三类，即随机数表方法、物理方法和数学方法。目前，最为常用的是数学方法。用数学方法产生随机数是通过数学递推式运算实现的，而由此产生的数值序列到一定长度以后，或退化为零，或周而复始地出现周期现象。因此，由数学方法产生的随机数并不是真正的随机数。但是，只要其周期充分长，而且能够通过[0，1]上均匀分布变量所需满足的各种统计检验，如分布的均匀性、抽样的随机性、独立性以及前后的一致性等，则从统计模拟的角度来看，就可以把它们当作随机数使用。为了和真正的随机数相区别，通常把用数学方法产生的随机数称为“伪随机数”。用数学方法产生的“伪随机数”的特点是速度快，占用计算机的内存小，对所模拟的问题可以进行复查，一般又有较好的概率统计性质。
用数学方法产生伪随机数的方法有多种，包括迭代取中法、移位法和同余法，而最为常用的是同余法。同余法包括乘同余法、加同余法、混合同余法等[2]。
本次采用混合同余法，其递推公式为：
xi+1=λxi+c |
(2) |
γi+1=xi+1/m (mod m) |
(3) |
其中，x,λ,c和m都取非负整数。当c=0时，上式即为乘同余法的递推公式。
按上述方法由计算机产生的伪随机数序列为[0，1]上的均匀随机数，在实际问题中尚需变换为给定或已知分布的随机样本值，即进行随机抽样。由于本文研究的水文地质参数大都服从均匀分布或正态分布，因此只需介绍正态随机变量的抽样。
正态随机变量的随机抽样方法有：哈斯汀(Hasting)有理逼近方法、统计近似方法和变换方法[3]。本次采用统计近似方法。产生均值为μ、方差为σ2的正态分布随机数Y的计算公式为：
|
(4) |
其中n足够大。通常取n=12，其近似程度已是相当好了，此时有：
|
(5) |
其中RNi为0到1之间的均匀分布的随机数[2]。
3 地下水允许开采量确定的风险分析 对地下水供水水源地的持续运行来说，风险Pf就是指供水工程承担的供水外来负荷L(Loading)大于工程本身的承载能力R(Resistance或Capacity)的概率：Pf=P(L>R) |
(6) |
式中：L与R均可为包括众多水文地质参数的函数。对地下水供水水源地来说，L主要为地下水的开采量，R则为地下水的允许开采量。因此，允许开采量的确定是非常重要，其数值的准确与否直接关系着水源地能否持续运行。允许开采量的计算方法有很多，但每一种都包含着一些不确定性因素。例如在最简单的通过计算地下水的补给量来确定允许开采量的方法中，要涉及到降水入渗补给量、河道渗漏补给量、灌溉水回渗量、其它含水层的越流补给量、地下水的侧向径流补给量等。
本文作者长期在山东省济宁市开展地下水资源评价与此同时管理工作，掌握了大量的地质及地下水的动态资料等，且深感过去的资源评价仅局限于确定性模型研究的不足。为此，选定济宁市作为研究实例。对于济宁市的深层承压含水层来说，地下水的补给量主要为来自浅层含水层的越流补给量和深层含水层本身的侧向径流补给量[4]，即：
Q补=ΔH×b×F+T×J×L |
(7) |
式中：Q补为深层承压含水层的补给量(×104m3/d)；ΔH为深层承压含水层与浅层含水层之间的水头差(m)；b为越流补给系数(1/d)；F为计算区面积(m2)；T为深层承压含水层的导水系数(m2/d)；J为深层承压水的水力坡度；L为计算区的周长(m).
式中的F和L是根据实际情况给定的确定数，而ΔH、b、T和J是反映当地水文地质条件的变量或参数，要经过野外的观测、抽水试验和绘制等水位线图等方法来确定，在这些过程中，都存在不确定性，往往只能根据掌握的观测数据和具体的水文地质条件，合理地确定其取值范围。如对济宁市的深层承压含水层，具体取值范围见表1(1987年的值)，表中各参数单位见上式的说明。
| 表1 含水层的有关参数值 | ||||
| 参数 | ΔH | b | T | J |
| 取值 | 10～15 | (1.3～3)×10-4 | 759～1597 | 0.001～0.0038 |
常规的水文地质计算方法是取上表中各参数的平均值，这样计算的地下水补给量为16.4×104m3/d.依据可开采系数法(本例取0.9)，其允许开采量可用下式直接估算，即：
Q允许=0.9×Q补给=14.76×104m3/d |
(8) |
即1987年济宁市深层地下水的允许开采量为14.76×104m3/d.而1987年的实际开采量为14.35×104m3/d，小于或近似等于允许开采量。实际的情况是1987年深层地下水的水位却仍在下降，其主要原因是允许开采量确定过程中，存在着一些不确定性因素，而传统的水文地质计算中并没有充分考虑这些不确定性，求得的数值可能偏大。正是由于这些不确定性因素的存在，使得在采取一定的开采量条件下，都存在有超出允许值风险。显然实际开采量愈小，则风险愈小。因此，必须对其进行风险分析。
由于灵敏度可度量一种因子的变化对另一因子的影响程度，因此本文采用了灵敏度分析来筛选主风险因子。某一模型因变量对模型输入参数的灵敏度可用因变量对该输入参数的偏导数来表示，即：
|
(9) |
式中的Xi.k为模型的因变量对第i观测点，第k参数的灵敏度系数。由于不同的参数其单位不同，这样不同参数的灵敏度系数就无可比性，因此采用下列标准化无量纲形式[5]，即
|
(10) |
据此式可计算出各因子的灵敏度系数，结果见表2.
| 表2 灵敏度系数计算结果 | ||||||||
| 项目 | ΔH | b | T | L | ||||
| 系数值 | 0.59 | -0.59 | 0.59 | -0.62 | 0.41 | -0.41 | 0.41 | -0.41 |
从表2中的结果可知，四个因子的灵敏度系数值较接近，表明它们对补给量的计算结果的影响都较显著。因此，将其都作为随机变量来考虑。
为此，用蒙特卡罗方法对上述四个因子进行模拟。限于客观资料条件，未能分别对各参数作统计分析，确定各自的概率分布概型，我们假定各参数服从均匀分布，且相互独立，则利用公式(5)计算求得的地下水补给量以及相应的概率见表3.
|
|||||||||||
| 补给量 | >16.4 | >14.76 | >14.35 | >12.12 | >9.89 | >8.32 | >6.65 | ||||
| 概率/% | 50.22 | 56.42 | 59.40 | 70.00 | 80.00 | 90.00 | 100.0 | ||||
从表3的结果可知，补给量大于16.4×104m3/d的概率为50.22%，由于我们假定各参数服从对称的均匀分布，所以，以参数的平均值计算的传统方法求得的补给量的概率约50%，超过或不足的机会基本相同。经风险分析得出：大于传统方法确定的允许开采量14.76×104m3/d的概率也只有56.42%，仅提高了6%.1987年实际开采量为14.35×104m3/d，相应的概率为59.40%，也就是说，该开采量的保证程度只有60%左右。当开采量仅为6.65×104m3/d时，才有接近100%的把握(不考虑其它条件的变化情况下).也就是说，当以1987年的实际开采量14.35×104m3/d开采时，虽然其小于或近似等于传统方法确定的允许开采量，但发生环境问题的可能性仍有40%.
当然，应该考虑到1987年的降水量仅547.4mm，比频率为75%对应的587.0mm还小，表明1987年为枯水年。若遇丰水年，浅层地下水的水位将更高，浅层与深层地下水的水头差更大，则越流补给量加大，使得以上确定的允许开采量保证程度提高。
在济宁市，因过量开采引起的环境问题主要有地下水水位的持续下降、地下水污染和地面沉降。这些环境问题已引起巨大的经济损失(限于篇幅，不再赘述).
4 结语以济宁市为例所进行的地下水允许开采量确定的风险分析，所得结果较传统水文地质学方法计算的结果，更符合实际情况。也为在今后的地下水资源评价工作中，更加合理地、可靠地确定水源地的允许开采量提供了一种切实可行的方法，使与地下水开采有着决策的失误减小到最低程度。研究认为，只有进行了有关的风险分析研究工作，才能提供各种地下水开发利用决策的风险信息，并据此提出实现地下水可持续开采的具体措施，使地下水能够实现预期的可持续开采的目标和产生应有的社会、经济和环境效益，不致造成无法挽回的损失。
同时，本课题的研究为传统的水文地质学开辟了一个析的研究领域，可有效地弥补以往确定性模型的不足。
参考文献:[1] 徐钟济。蒙特卡罗方法[M]。上海：上海科学技术出版社，1985.
[2] 徐士良。FORTRAN常用算法程序集(第二版)[M]。北京：清华大学出版社，1995.
[3] 赵国藩。工程结构可靠性理论与应用[M]。大连：大连理工大学出版社，1996.
[4] 长春地质学院学报编辑室。长春地质学院科学研究论文集(1988)[C]。长春：吉林科学技术出版社，1991，190～213.
[5] Zheng Chunmiao,G D Bennett.Applied contaminant transport modeling—theory and practice[M]。Van Norstrand Reinhold, 1995，287～311.
网页制作
