一种地下水模型输出不确定性分析方法

摘要

本发明公开一种地下水模型输出不确定性分析方法，其将频率分析与敏感性分析方法相结合。其中频率分析过程包括参数估计与假设检验过程，选择7中具有代表性的分布函数作为备选概率密度函数，并进行假设检验，为地下水位序列选择合适的概率密度函数。敏感性分析方法包括逐步回归分析和互熵分析。逐步回归分析能够对输入变量的不确定重要性进行趋势分析，互熵分析能够对不确定性因子进行良好的识别。该方法能够弥补传统不确定性分析方法在研究内容上的不足，分析输出变量的概率分布特征，识别影响输出变量概率分布的关键不确定性因子。从而更好的理解地下水模型不确定性的产生与来源，为地下水模拟的数据收集工作提供反馈，减少模型的不确定性。

说明书

技术领域

本发明涉及一种模型输出不确定性分析方法，具体涉及一种对地下水模型输出不确定性分析方法。

背景技术

随着数值模拟方法及计算机技术的发展，地下水模型已经成为地下水资源管理与规划的基本工具。地下水系统是一个复杂、开放的巨系统，它受到水文、气象、地质条件及人类活动的影响。由于观测资料的限制，当利用地下水模型进行数值模拟时，模拟结果往往与实际观测存在偏差，即地下水模型的不确定性，从而影响地下水模拟的可靠性。因此，为了降低地下水模拟与预测的风险，地下水模型的不确定性分析成为当前国内外研究的热点。

目前，地下水模型不确定性分析的方法主要有：①参数不确定性分析方法，对由模型参数、边界条件等导致的不确定性进行分析，对模型输入参数的后验概率分布进行反演，并对模型的输出进行预测（Beven, K., Freer, J., 2001. Equifinality, data assimilation, and uncertainty estimation in mechanistic modeling of complex environmental systems using the GLUE methodology. J Hydrol, 249(1-4): 11-29; Hassan, A.E., et al, 2008. Uncertainty assessment of a stochastic groundwater flow model using GLUE analysis. J Hydrol, 362(1-2): 89-109）；②概念模型不确定性分析方法，考虑模型结构的不确定性，综合多个概念模型的结果进行不确定性预测（Rojas, R., et al, 2008, Conceptual model uncertainty in groundwater modeling: Combining generalized likelihood uncertainty estimation and Bayesian model averaging, Water Resour Res, 44(12); Ye, M., et al, 2010, A Model-Averaging Method for Assessing Groundwater Conceptual Model Uncertainty, Ground Water, 48(5), 716-728）；③模型输出的敏感性分析方法，对模型的输出（如水位、浓度）进行敏感性分析，识别地下水模型的不确定性因子（Mishra, S., et al, 2009, Global Sensitivity Analysis Techniques for Probabilistic Ground Water Modeling. Ground Water, 47(5): 730-747）。这些方法被广泛应用于地下水数值模拟领域，且受到水文学者们的认同。尽管如此，这些方法不能对地下水模型输出的概率分布特征及其影响因子进行分析。模型输出的概率分布特征作为模型不确定性的最终表现形式，决定了模型模拟及预测的准确性和精确性。并且，模型输出概率分布的影响因素能够更好的理解地下水模拟不确定性的产生及来源。

发明内容

发明目的：本发明的目的在于针对现有技术的不足，提供一种能够更好地理解模型不确定性的来源，识别地下水模型的关键不确定性因子，从而控制不确定性产生的地下水模型输出不确定性分析方法。

技术方案：本发明所述的地下水模型输出不确定性分析方法，包括以下步骤：

频率分析：

（1）在地下水模型输出变量的频率分析过程中，选用正态、对数正态、gamma-2，对数gamma-2、p-III、对数p-III及均匀分布作为备选的概率密度函数；

（2）根据最大熵原理（Principle of Maximum Entropy, POME），对正态、对数正态、gamma-2，对数gamma-2、p-III、对数p-III及均匀分布函数进行参数估计；

（3）选择Chi-squared检验，为地下水模型输出变量选择合适的概率密度函数；

敏感性分析：

（4）选择逐步回归分析方法，对地下水模型输出变量的概率分布特征与模型输入变量进行线性相关分析；

（5）选择互熵分析方法，对下水模型输出变量的概率分布特征与模型输入变量的相关关系进行分析。

上述步骤1）中，选择的7中概率分布函数具有不同的分布特征，具有代表性和广泛性。因此，能够防止由于备选分布函数不够充分而导致对地下水模型输出变量概率分布估计的偏差。

上述步骤4）中，采用基于线性模型的逐步回归分析，能够对变量的相关关系进行趋势分析，并且与互熵分析的结果进行对比，识别模型输出变量概率分布的不确定性因子。

上述步骤5）中，互熵分析基于信息熵理论，能够克服线性模型在相关分析中的限制，对复杂多变量的非线性关系进行良好的识别。

本发明通过频率分析方法获得地下水模型输出的概率分布特征，采用逐步回归分析和互熵分析方法识别模型输出概率分布的关键影响因子。模型输出变量的概率分布作为地下水模型不确定性的直接体现，它是不确定性研究的最终目标。变量的概率分布范围、形状和位置由该变量的概率密度函数决定。因此，输出变量概率分布的影响因子能够在根本上控制模型输出不确定的来源及形成过程。地下水模型输出变量概率分布的敏感性分析能够更好的理解模型不确定性的来源，识别地下水模型的关键不确定性因子，从而有利于控制不确定性的产生。

本发明将频率分析方法与敏感性分析方法相结合，建立一种模型输出不确定性分析方法。本发明相对于现有的技术具有以下有益效果：

（1）将频率分析用于地下水模型的不确定性分析中，能够获得地下水输出变量所服从的概率密度函数，对输出变量的概率分布特征进行定量描述；

（2）与常规的地下水模型敏感性分析不同，本次敏感性分析的对象是地下水模型输出变量的概率分布，能够从更深层次理解模型的输入对输出变量的影响；

（3）采用互熵分析方法，能够比传统方法更加有效、准确的识别模型输出变量概率分布的关键不确定性因子；

（4）通过识别地下水模型输出变量概率分布的不确定性因子，能够更好的解释模型输出不确定性的来源，为实际地下水模拟的数据收集工作提供反馈，有利于减少模型的不确定性。

附图说明

图1为本发明实施例中地下水模型示意图。

图2为本发明实施例中地下水位序列均值和方差的逐步回归分析的结果，(a), (b), (c)分别表示输出变量y为观测井1, 2, 3的水位序列的均值，(d), (e), (f)分别表示输出变量y为观测井1, 2, 3的水位序列的方差。

图3为本发明实施例中地下水位序列均值和方差的互熵分析结果，(a), (b), (c)分别表示输出变量y为观测井1, 2, 3的水位序列的均值，(d), (e), (f)分别表示输出变量y为观测井1, 2, 3的水位序列的方差。

具体实施方式

本发明基于一个理想的地下水流模型，对模型输出变量的概率分布进行分析。地下水模型的输出变量为地下水位序列。通过频率分析方法获取地下水位序列的概率分布特征，再对其进行敏感性分析，识别关键的不确定性因子。

理想地下水模型的建立

如图1所示，该理想三维地下水模型在平面上为一个矩形，长3600m，宽1800m，剖分为20 m* 20m的栅格单元。模型在垂向上延伸53m，共分为三层，从顶部至底部依次为：潜水含水层（1-30m），弱透水层（31-33m），承压含水层（34-53m），各含水层产状水平。假设含水介质为非均质，渗透系数场在层内符合平稳分布。通过空间随机函数理论来表达渗透系数K的空间变异性，采用各向同性的指数协方差函数描述lnK在各层内的分布。表1为模型各层内lnK协方差函数的参数。模型各层的给水度（Sy）和储水系数（S）均假设为均质分布。

表1 地下水模型各含水层内渗透系数K空间分布参数

如图1所示，模型的北部和南部均为给定流量边界，分别为入流和出流边界。模型西部为通用水头边界（general head boundary）。东部为一条宽20m的河流，河床厚度为2m，河床底部高程为45m。模型顶部均匀接受降水补给。在模型承压含水层内设置一口抽水井和三个观测点。表2为模型参数的取值区间，模型边界条件的设置如表3，表4所示。

表2 地下水模型参数的取值区间

表3 地下水模型边界参数的取值区间

参数取值区间河床渗透系数 (m/d)0.1-10.0通用水头边界导水性50.0-5000.0

表4 地下水模型边界条件的取值区间和变化范围

边界取值区间变化范围降水率 (mm/month)5.0-50.0[-100%, 100%]抽水率(10⁴m³/month)0.2-2.0[-100%, 100%]入流率 (10⁴m³/month)0.1-1.0[-100%, 100%]处理率 (10⁴m³/month)0.1-1.0[-100%, 100%]河流水位 (m)40.0-44.0[-1.0m, 1.0m]通用边界水头(m)46.0-50.0[-1.0m, 1.0m]

1.2 频率分析方法

频率分析是水文预报中广泛使用的一种分析方法。频率分析的内容是为观测数据系列(x₁, x₂, …, x_n)选择合适的概率密度函数，从而分析变量的分布特征及进行预报。频率分析的主要步骤有：（1）选择合适的概率密度函数；（2）根据最大熵原理，基于观测数据序列，对备选概率密度函数进行参数估计；（3）对变量进行不确定性预报。

备选概率密度函数的参数估计过程如下：

（1）均匀分布

概率密度函数为： (1)

参数估计为： (2)

θ₁，θ₂表示均匀分布的上、下边界，(x₁, …, x_n)为数据系列，n为样本数量，X_max，X_min分别为数据系列的最大值与最小值；

（2）正态分布

概率密度函数为： (3)

参数估计为：，μ表示样本的均值，σ²表示样本的方差；(x₁, …, x_n)为数据系列，n为样本数量；

（3）gamma-2分布

概率密度函数为： (4)

基于最大熵原理（Principle of Maximum Entropy, POME），gamma-2参数估计为：

 (5)

α，β分别表示gamma-2分布函数的形状参数和尺度参数；

（4）P-III分布

概率密度函数为： (6)

基于POME，P-III分布的参数估计为：

(7)

α，β，c分别表示P-III分布函数的形状参数、尺度参数和位置参数，σ_x为样本标准方差，(x₁, …, x_n)为数据系列，n为样本数量，C_s为偏态系数，C_v为离差系数；

（E）对数正态、对数gamma-2、对数p-III分布

首先对数据系列进行对数转换，再对对数正态、对数gamma-2及对数P-III分布进行相应的参数估计。

对备选的概率密度函数进行Chi-squared检验，其步骤如下：

（1）用k-1个数将数轴分成k个区间(-∞, t₁], (t₁, t₂], …, (t_k-2, t_k-1], (t_k-1, +∞]，k ≈ 1.87(n-1)^0.4，n为样本数量；

（2）计算样本序列(x₁, x₂, …, x_n)落入各个区间内的数量n_i, i=1, 2, … , k，并计算备择概率密度函数f₀(x)在各个区间的概率：

 (8)

（3）计算Chi-squared统计量χ²：

 (9)

（4）指定置信水平α，若p (χ² ≥ χ2 1-α) ≥ α，认为样本序列服从备择概率密度函数，否则拒绝。

敏感性分析方法

逐步回归分析是一种简单、方便的线性回归方法。它假设输出变量与输入变量满足线性关系：

 (10)

y′为拟合的输出变量，x_j为输入变量，b_j为回归系数，j = 1,.., k，表示回归变量的数量。逐步回归分析将回归变量逐个引入，每一步都要进行测试确保每个回归系数都具有显著性，在一定的置信水平下，进行t分布测试。回归模型的建立可以分为两个步骤：（1）不断将具有最大偏相关系数的变量导入回归模型；（2）对模型中所有变量进行显著性检验，排除偏相关系数不显著的变量，直到模型内所有变量均具有显著重要性，且未进入模型的变量的偏相关系数均不显著。

互熵分析是一种基于信息熵理论的敏感性分析方法。通过建立列联表（contingency tables）来确定数据(X, Y)之间的相关关系。将输入变量x的取值范围分为i个区间，输出变量y的取值范围分为j个区间。建立i行j列的列联表，统计样本数据落入每个表格单元（i, j）的数量N_ij, p_ij表示样本数据落入区间（x_i, y_j）的概率，p_ij = N_ij/N, N为样本数量。p_i., p_.j分别表示样本落入区间x_i, y_j的概率，p_i. = N_i./N, p_.j = N_.j/N, N_i., N_.j分别表示样本落入区间x_i, y_j的数量。

基于信息论，变量x, y的信息熵H(x), H(y)为：

 (11)

变量x, y的联立信息熵H(x, y)为：

  (12)

在信息论中，变量间的互熵I(x, y)表示它们的相互依赖程度。x, y的互熵表示由于x（或y）的信息而导致y（或x）不确定性的减少。即：

 (13)

在互熵分析中，输入变量x对与输出变量y的不确定重要性（uncertainty importance）有以下两种指标表示：

（1）不确定性系数U(x, y)

 (14)

（2）R统计量

 (15)

这两个指标均在[0, 1]间变化，表示x与y独立或完全相关。

模拟方案设计

基于Monte Carlo模拟技术，对获取的地下水位序列进行频率分析，并对其概率分布进行敏感性分析，识别关键不确定性因子。Monte Carlo模拟步骤如下：

（1）对地下水模型进行空间离散和时间离散。模拟期设为100年，一个月为一个应力期，每个应力期分为6个时间段。输出水位为每个应力期内6个模拟水位的平均值。1-100年内某个特定月份的水位平均值序列作为研究对象。

（2）设置模型参数，包含渗透系数、给水度、储水系数、河床渗透系数、通用水头边界的导水率。根据模型含水层的渗透系数空间变异特征（表1），渗透系数随机场由直接傅里叶变化方法（direct Fourier transform）生成。其余的模型输入参数在相应的取值范围内（表2，表3）均匀随机抽样生成。

（3）设置模型边界条件，包含抽水率、边界流入量、边界流出量、河水水位、通用边界水头。由模型边界设置的范围（表4），边界条件的赋值分为两个步骤：①在相应边界条件的取值区间内，为每月（1-12月）均匀随机抽取一个均值；②基于该均值，在相应的变化范围均与抽样100次，作为模拟期内该月份的边界条件值。

（4）基于Modflow-2005建立地下水数值模拟模型，运行模型并对各观测点的地下水位序列进行频率分析。基于地下水位序列，为每个备择概率密度函数进行参数估计，并进行Chi-Square检验。设置显著水位α为0.05，若所有备择函数均未能通过Chi-Square检验，则将该水位序列标记为未知概率分布。

（5）重复步骤2—步骤4，直到服从每个备选概率密度函数的水位序列所占的比例保持稳定。

（6）收集Monte Carlo模拟的数据样本，包含水位序列、水位序列所服从的概率密度函数及其参数、水位序列所对应的模型参数与边界条件。

（7）对水位序列概率密度函数的参数与模型输入变量（模型参数与边界条件）进行逐步回归分析和互熵分析。

地下水位序列概率分布的敏感性结果与讨论

通过对地下水位序列的频率分析，83.83%的水位序列服从正态分布，服从对数正态、gamma-2、对数gamma-2、P-III、对数P-III、均匀分布及未知分布的水位序列所占的比例分别为： 83.72%，83.48%，83.28%，1.23%，1.03%，9.44%，和13.05%。因此，大多数的水位序列服从正态分布（或对数正态分布）。并且，当gamma-2分布的参数α（shaper parameter）为一个正整数n时，可以将其看做n个指数分布的总和。由中心极限定律，当n较大时（在本例中约为20000），gamma-2分布将趋近于正态分布。因此，这两个分布所占的比例相近。

选择服从正态分布的水位序列作为分析对象，将水位序列的均值与方差作为输出变量，地下水模型的输入参数与边界条件作为输入变量（如表5所示），分别进行逐步回归分析和互熵分析。

表5 地下水位序列概率分布敏感性分析的输入变量及其编号

变量编号变量编号潜水层K均值1入流率均值11弱透水层K均值2入流量方差12承压层K均值3出流率均值13潜水层给水度4出流量方差14潜水层储水系数5抽水率均值15承压层储水系数6抽水率方差16河床渗透系数7河水位均值17通用水头边界导水率8河水位方差18降水率均值9通用边界水头均值19降水率方差10通用边界水头方差20

图2为地下水位序列均值和方差的逐步回归分析结果。(a), (b), (c)分别表示输出变量y为观测井1, 2, 3的水位序列的均值，(d), (e), (f)分别表示输出变量y为观测井1, 2, 3的水位序列的方差。对于水位序列的均值，前3位不确定性因子为变量1, 3, 19。对于水位序列的方差，前3位不确定性因子为变量1, 4, 7。回归模型中，变量的回归系数的绝对值一般较小，对于均值回归模型，最大值约为0.3，平均值约为0.12。对于方差模型，最大值约为0.4，平均值约为0.18。这主要是由于这些水位序列的均值和方差与输入变量之间并非简单的线性关系，而是更为复杂的非线性关系。

图3为地下水位序列均值和方差的互熵分析结果。(a), (b), (c)分别表示输出变量y为观测井1, 2, 3的水位序列的均值，(d), (e), (f)分别表示输出变量y为观测井1, 2, 3的水位序列的方差。水位序列的均值与方差的前3位不确定性因子分别为变量1, 3, 19，和变量1, 4, 7，这与逐步回归分析的结果相同。但是，回归系数较小的变量与互熵分析具有不同的不确定重要性识别结果。因此，逐步回归分析能够对关键的不确定性因子进行识别，而不能对普通的不确定性因子进行分析。

结论

    针对传统地下水模型不确定性分析方法在研究内容上的缺陷，本发明将频率分析与敏感性分析方法相结合，对地下水模型输出变量的概率分布进行敏感性分析，识别关键的不确定性因子。通过一个理想地下水模型，将观测点地下水位作为模型输出变量，选择7中常见的分布作为备选概率密度函数，并进行假设检验。对服从正态分布的水位序列的均值与方差进行逐步回归分析和互熵分析。该方法能够获取输出变量的概率分布特征，有效地识别输出变量的不确定性因子。从而为不确定性分析的数据收集工作提供反馈，减少模型的不确定性。

如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。