时间:
一、引言
随着近几十年计算机软硬件和数值计算方法的发展,地下水建模成为人类深入理解开放、复杂的地下水系统的有效手段。但任何数学模型或数值模拟必然都只是真实地下水系统的简化,不确定性不可避免地由各种模型输入因素引入,导致模型的输出或预测结果也带有不确定性。为了更好地利用有限的资金和资源去最优化地降低预测结果的不确定性,需要测算各个模型输入因素的重要性,找出最重要即预测结果对之变化最敏感的不确定模型输入因素,需要引入敏感性分析。
从不确定性角度来理解,敏感性分析是指比较或计算不同的模型不确定性输入因素对输出结果不确定性的贡献的方法。现有的在地下水领域应用的敏感性分析方法大致分为 2 种:局部敏感性分析和全局敏感性分析。全局敏感性方法相较于局部敏感性方法虽然更复杂、计算量更大,但它能够得到适用于模型参数全部取值范围的结论,并考虑输入因素之间的相互作用,因此在最近数十年得到了重视和应用。
尽管全局敏感性分析已经被广泛应用于地下水建模,但多数的敏感性分析研究仅针对地下水模型参数的不确定性,而忽略了另外 2 类地下水建模中所遇到的主要模型输入不确定性:地下水模型本身的不确定性和情景的不确定性。在地下水建模中,情景不确定性被定义为一种影响地下水系统,并存在于未来的多种可能的自然或者人为(比如降雨、抽水等)外部状态导致的不确定性。模型的不确定性是指由于缺乏模拟区域或系统的数据或者对其的理解,多种可能的概念或数学模型并存导致的不确定性。
DAI 等提出了先情景、后模型,最后参数的不确定性多层级关系,并将方差分析法与多层级的不确定性分析框架相结合,推导出新的全局敏感度系数来评估地下水模型中参数的敏感性。但该研究仍然没有针对情景和模型的不确定性进行评价。除此之外,前人的多层级不确定性分析框架只能评估系统中模型的敏感性,无法对模型中各种过程的概念或数学模型的敏感性进行分析。在地下水建模中,通常存在多个概念或数学模型用于刻画同一个地下水流动过程,目前少有研究者对模型中各个过程的概念模型和数学模型进行敏感性分析。
为此,笔者提出了一种改进的层级制敏感性分析方法,量化地下水模型中情景、模型和参数的不确定性,并利用参数分组的方法在层级制不确定性分析框架中评估地下水模型涉及的各个过程及其概念或数学模型的不确定性对模型输出不确定性的贡献,并将提出的新方法应用于理想的地下水污染运移模拟模型以说明该方法的运用。
二、改进的层级制全局敏感性分析方法
2.1 基于方差分析的层级制敏感性分析
将不同来源的不确定性置于多层级框架的不同层,然后利用方差分析将输出结果的不确定性按照层级制框架分配到不同的不确定性输入上,计算各部分的方差在总方差中的占比,并把这种百分比定义为敏感度系数来精确量化不同来源的不确定性的重要性。
按照情景 - 模型 - 参数的多层级框架,利用总方差定律对模型输出结果的总方差进行分解。总方差共被分解为 4 项:情景贡献的方差、模型贡献的方差、参数贡献的方差和观测数据贡献的方差。由于观测数据贡献的方差独立于情景、模型和参数贡献的方差,故本研究不考虑观测数据的不确定性。使用其余 3 部分的方差与总方差的比值定义情景的敏感度系数、模型的敏感度系数和参数的敏感度系数,3 种敏感度系数在数学上的和为 “1”,该条件也是该方法在应用时的检验标准。
2.2 过程子模型的敏感度系数
一个复杂的地下水系统通常由多个过程耦合组成,例如:地下水流动过程、溶质运移过程、能量传输过程等。研究者们通常根据有限的观测信息利用数学方程简化水流运动、污染(溶质)运移等真实过程,这就造成每种过程的数学模型并不唯一,即各种过程的概念或数学模型存在不确定性。目前地下水领域少有研究针对此类不确定性。本研究利用参数分组的方法扩展参数的敏感度系数,在多层级不确定性分析框架中评估参数组的敏感性,从而评估模型中各种过程的概念或数学模型的不确定性对输出结果不确定性的贡献。
将描述某一过程的概念或数学模型包含的所有参数分为一组作为条件对参数贡献的方差进一步分解,等式右边第 1 项代表该参数组,即该数学模型或概念模型(称为过程子模型)的不确定性贡献的方差,右边第 2 项代表除该参数组以外其余参数产生的方差。定义过程子模型的敏感度系数为该参数组不确定性贡献的方差与参数贡献的方差的比值。
三、理想的地下水污染运移模拟模型介绍
通过理想的地下水污染运移模拟模型来说明上述方法的运用。宽度 L=10000 m,地下水流和溶质运移考虑为一维稳定流溶质运移情况。均匀稳定入渗的条件下,考虑 3 种情景,以情景 1 作为基准值,其降雨强度 P 为 1524 mm/a,较干燥和湿润的情景下,降雨强度分别是基准值的 80% 和 180%。两侧河流设定为定水头边界,h₁=330 m,河水位 h₂受到融雪过程的影响。在潜水含水层剖面的中间位置 x₀=5000 m 处存在一个污染源,污染物为乙烯(ETH)。
对于降雨入渗补给过程,考虑 2 种数学模型来描述(R₁和 R₂),入渗补给数学模型 R₁中考虑参数取值的不确定性,参数 a 设定为服从均值为 16.88、标准差为 1 的正态分布;入渗补给数学模型 R₂中参数 b 设定为服从均匀分布的随机参数,分布区间为 [0.1,0.2]。
潜水含水层的地质情况共考虑 2 种数学模型来描述,模型 G₁代表整个河间地块为均质,使用统一的渗透系数 K 来描述;模型 G₂代表河间地块为非均质,以 x=7000 m 为分界处,左侧区域渗透系数为 K₁,右侧渗透系数为 K₂,渗透系数均设定为服从正态分布的随机参数。
河流 2 的水头 h₂受融雪径流过程的影响,河流 2 水位(h₂)和流量 Q 的关系如下,假设河流径流量 Q(m³/s)受上游融雪量控制,计算为 Q=C₁×C_sn×S×SVC×A,使用 2 种数学模型来评估日融雪量 S,即度日模型一般形式(SN₁)和融入辐射变量后的度日模型(SN₂),度日因子 f₁和 f₂分别设定为服从正态分布的随机参数,数学模型 SN₂中太阳短波辐射或者净辐射 Rₙ=80 W/m²,系数 r 设定为服从正态分布的随机参数。
综上,整个潜水含水层一维稳定流溶质运移模型由 3 部分组成:①降雨入渗补给过程可由 2 种数学模型刻画(过程子模型 R₁和 R₂);②地质过程考虑均质和非均质 2 种过程子模型(G₁和 G₂);③融雪模型可由 2 种数学模型刻画(过程子模型 SN₁和 SN₂)。通过组合以上 3 种过程的所有子模型得到 8 种模型用于描述该系统,为每个模型分配相等的先验权重,所考虑的输出(A)是 5100~5700 m 的所有位置(每 100 m 设置 1 个计算点)的乙烯质量浓度,在 MATLAB 软件中采用蒙特卡罗模拟进行编程计算。
四、结果与讨论
4.1 5400 m 处乙烯质量浓度预测结果的敏感性分析
4.1.1 预测结果的不确定性
不同情景下不同模型的预测结果不同,表明预测存在较大的不确定性。其中,包含均质地质过程的模型(均质模型)与包含非均质地质过程的模型(非均质模型)预测结果差异最为明显,非均质模型预测结果的均值明显小于均质模型的均值,且非均质模型的不确定性范围明显大于均质模型的不确定性范围。不同情景之间的预测结果差异并不明显,即在此案例中,情景对结果预测的影响程度并不大。
4.1.2 基于方差分析的层级制敏感性分析结果
模型的敏感度系数(Sₘ)显著大于情景和参数的敏感度系数(S_sc、S_pa),即模型的不确定性是输出结果不确定性最重要的来源;而情景的敏感度系数最小,意味着情景的不确定性可以忽略不计。仅执行 1000 次运算,敏感度系数就已经全部收敛,相对误差小于 0.01,模型运行时长仅 2.88 s,相较于前人对该案例的参数敏感性分析,计算效率大幅提升。
将总方差分解为情景间方差和情景内方差时,情景内方差几乎是结果方差的全部来源,情景间方差可以忽略不计。将情景内方差进一步分解为模型间方差和模型内方差时,模型间方差是情景内方差的主要来源,也是总方差的主要来源。模型内方差是 3 种情景下模型内方差的平均结果,非均质模型中所有参数的不确定性是模型内方差的主要来源。
基于方差分析的层级制敏感性分析结果显示,模型间的不确定性是输出结果的不确定性最重要的来源,这提醒建模者为了显著降低系统的不确定性,应该把更多的精力和资源放在模型结构上。模型内的不确定性即所有参数贡献的不确定性明显小于模型间的不确定性,这也提醒建模者针对该案例,如果花费过多的精力和成本在参数反演上是不明智的选择。另外,尽管非均质模型的方差对模型内方差的贡献更大,但这并不代表均质模型可以在建模的时候被忽略。虽然均质模型的方差贡献小,但值得注意的是,若将该系统的模型确定为均质模型中的一个,系统的总方差将明显减小,不确定性显著降低,这意味着均质模型在这样的系统中应该慎重选择,合理取舍。
4.2 5100~5700 m 范围内所有位置预测结果的敏感性分析
4.2.1 预测结果的不确定性
3 种情景下,各个模型预测的均值不尽相同,表明预测存在不确定性。比较不同情景下所有模型的结果可以发现,随着降雨强度的变化,模型的预测结果并未发生明显的变化,这意味着在 5100~5700 m 的空间范围中,情景间的不确定性很小。但所有模型的预测均值曲线在 5300 m 以后分离程度明显,这表明模型间的预测结果差异很大,即模型间的不确定性很大。
在该空间中大多数位置,均质模型的浓度均值明显高于非均质模型的均值。这是因为非均质模型的等效渗透系数相较于均质模型的渗透系数更小,模拟时间相同时,乙烯在均质模型中运移的更远,污染范围更大,故在靠近污染源的位置,非均质模型的预测浓度均值更大;而在远离污染源的位置,均质模型预测浓度均值更大。
4.2.2 基于方差分析的层级制敏感性分析结果
在 5100~5700 m 的空间范围内,情景的敏感度系数(S_sc)始终远小于其他 2 种系数且几乎接近于 0,这意味着情景的变化对输出结果的不确定性并不重要,即可以忽略情景的不确定性。在 5300 m 后的所有位置,对输出结果的方差贡献最大的是模型间的不确定性(Sₘ)。在 5200 m 处,模型的敏感度系数(Sₘ)明显小于参数的敏感度系数(S_pa),这意味着在该位置系统总方差的主要来源是模型内方差,即所有参数的不确定性贡献的方差。
将总方差分解为情景内方差和情景间方差,所有位置的结果均显示,95% 以上的总方差来自于情景内方差。将情景内方差进一步分解为模型间方差和模型内方差,可知模型间方差始终是情景内方差的主要来源。各部分方差总体上随着距离的增大先增大后减小,在 5400 m 处达到峰值。将 3 种情景下的模型内方差进一步分解为 8 个模型的方差,在 5400 m 及其之前的所有位置,非均质模型预测结果的方差是模型内方差的主要贡献者,而在 5400 m 以后的位置,均质模型预测结果的方差成为模型内方差的主要来源。
4.3 过程子模型的敏感性分析
为了加深对模型结构的理解,使用本研究提出的过程子模型的敏感度系数精确评估了不同地质模型(G₁,G₂)、入渗补给模型(R₁,R₂)和融雪模型(SN₁,SN₂)对输出结果不确定性的贡献。结果显示,所评估的位置不同,敏感度系数的排名不同。总体而言,地质模型始终是最重要的模型,入渗补给模型的敏感性几乎接近于 0。这是因为地质模型中的渗透系数直接影响着污染物乙烯的运移速度,而入渗补给模型涉及的参数通过影响地下水位,间接影响乙烯的运移速度。在 5400 m 及其之前的所有位置,地质模型中的非均质子模型 G₂是所有过程子模型中最重要的模型,其敏感度系数显著大于其余模型的敏感度系数;而在 5500 m 及其以后的位置,均质子模型 G₁是最重要的模型。
综上所述,建议将模型校准的主要精力和资源集中在地质过程的数学模型上。为了提高建模中参数反演效率,可以考虑将入渗补给模型和融雪模型中的不确定性参数固定为某个经验值,优先反演地质过程的数学模型中包含的参数;或者也可以舍弃不敏感的过程子模型,选取更重要的数学模型来描述过程。
五、结语
改进的层级制敏感性分析方法能够量化并比较情景、模型、参数 3 种不同类型的输入不确定性对输出结果不确定性的贡献,同时评估地下水模型中各个过程的概念或数学模型的不确定性对输出结果不确定性的影响程度,为地下水模型提供精确全面的敏感性分析。相较于以往的参数敏感性分析,改进的层级制敏感性分析方法的计算效率大幅提升,并且能够提供更有用的敏感性信息。
本研究使用一个理想的地下水污染运移模拟模型展示了该方法的运用,结果表明,模型的不确定性是该案例中乙烯浓度预测不确定性的主要来源。对所有过程子模型进行分析后发现,地质模型是模型输出结果不确定性最重要的来源,其中非均质子模型相较于均质子模型对结果的不确定性贡献更大。这些信息有助于指导建模者利用有限资源去优先减少该案例的模型不确定性,尤其是描述地质过程的数学模型的不确定性,从而最优地减少预测浓度结果的不确定性。
刘玉姣;戴恒;李跃东;崔节波;文章,中国地质大学(武汉);长江流域环境水科学湖北省重点实验室;生物地质与环境地质国家重点实验室;生态环境部华南环境科学研究所,202405