水文仿真软件：MODFLOW_（8）.灵敏度分析与不确定性评估

最新推荐文章于 2025-05-16 06:37:56 发布

kkchenjj

最新推荐文章于 2025-05-16 06:37:56 发布

阅读量718

点赞数 10

分类专栏：环境仿真文章标签：仿真模拟环境仿真模拟仿真

本文链接：https://blog.csdn.net/weixin_42749425/article/details/147967527

版权

环境仿真专栏收录该内容

475 篇文章

订阅专栏

灵敏度分析与不确定性评估

在水文仿真软件的二次开发中，灵敏度分析和不确定性评估是两个至关重要的步骤。通过这些分析，可以更好地理解模型参数对模拟结果的影响，以及模型预测的可靠性。本节将详细介绍如何在MODFLOW中进行灵敏度分析和不确定性评估，并提供具体的代码和数据样例。

在这里插入图片描述

灵敏度分析

灵敏度分析的概念

灵敏度分析（Sensitivity Analysis, SA）是评估模型参数变化对模型输出影响的过程。通过灵敏度分析，可以确定哪些参数对模型输出最为敏感，从而为模型校准和优化提供指导。MODFLOW中常见的灵敏度分析方法包括局部灵敏度分析和全局灵敏度分析。

局部灵敏度分析

局部灵敏度分析（Local Sensitivity Analysis, LSA）是通过改变单个参数的值来评估其对模型输出的影响。这种方法简单直接，但只能提供参数在特定点的敏感性信息。

实现步骤

定义参数和输出变量
改变参数值
运行模型
记录输出变量的变化
计算灵敏度系数

代码示例

假设我们要对MODFLOW中的一个参数（例如储水系数）进行局部灵敏度分析。以下是一个Python脚本示例，使用Flopy库来实现这一过程：


import flopy

import numpy as np



# 定义模型参数

model_name = 'example_model'

mf = flopy.modflow.Modflow.load(f'{model_name}.nam', version='mf2005', exe_name='mf2005.exe')



# 定义储水系数参数

storage_coeff = mf.lpf.sy.array



# 定义输出变量（例如地下水位）

output_variable = 'head'



# 定义参数变化范围

param_range = np.linspace(0.1, 0.9, 9)



# 存储结果

results = []



# 进行局部灵敏度分析

for param in param_range:

    # 修改储水系数

    mf.lpf.sy.array = param

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 记录结果

    results.append((param, np.mean(head)))



# 计算灵敏度系数

base_param = 0.5  # 基准参数值

base_head = np.mean([head for param, head in results if param == base_param])

sensitivity_coefficients = [(param, (head - base_head) / (param - base_param)) for param, head in results]



# 输出结果

for param, sc in sensitivity_coefficients:

    print(f'储水系数: {param:.2f}, 灵敏度系数: {sc:.2f}')

全局灵敏度分析

全局灵敏度分析（Global Sensitivity Analysis, GSA）是通过同时改变多个参数的值来评估其对模型输出的综合影响。这种方法可以提供更全面的敏感性信息，但计算成本较高。

实现步骤

定义参数和输出变量
生成参数组合
运行模型
记录输出变量的变化
计算灵敏度指标

代码示例

假设我们要对多个参数（例如储水系数和导水率）进行全局灵敏度分析。以下是一个Python脚本示例，使用Flopy库和SALib库来实现这一过程：


import flopy

import numpy as np

import pandas as pd

from SALib.sample import saltelli

from SALib.analyze import sobol



# 定义模型参数

model_name = 'example_model'

mf = flopy.modflow.Modflow.load(f'{model_name}.nam', version='mf2005', exe_name='mf2005.exe')



# 定义参数范围

problem = {

    'num_vars': 2,

    'names': ['sy', 'k'],

    'bounds': [[0.1, 0.9], [10, 100]]

}



# 生成参数组合

param_values = saltelli.sample(problem, 1024)



# 存储结果

results = []



# 进行全局灵敏度分析

for params in param_values:

    # 修改储水系数和导水率

    mf.lpf.sy.array = params[0]

    mf.lpf.k.array = params[1]

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 记录结果

    results.append(np.mean(head))



# 计算灵敏度指标

Si = sobol.analyze(problem, np.array(results))



# 输出结果

print('S1 (First-order indices):', Si['S1'])

print('ST (Total-order indices):', Si['ST'])

不确定性评估

不确定性评估的概念

不确定性评估（Uncertainty Analysis, UA）是评估模型参数的不确定性和模型结构的不确定性对模型输出的影响。通过不确定性评估，可以了解模型预测的可靠性，并为决策提供科学依据。

方法

常见的不确定性评估方法包括蒙特卡洛模拟（Monte Carlo Simulation, MCS）和贝叶斯方法（Bayesian Methods）。

蒙特卡洛模拟

蒙特卡洛模拟是通过随机抽样参数值来评估模型输出的不确定性。这种方法简单易行，但需要大量的模拟次数。

实现步骤

定义参数和输出变量
生成随机参数值
运行模型
记录输出变量的变化
统计分析结果

代码示例

假设我们要对储水系数和导水率进行蒙特卡洛模拟。以下是一个Python脚本示例，使用Flopy库和NumPy库来实现这一过程：


import flopy

import numpy as np



# 定义模型参数

model_name = 'example_model'

mf = flopy.modflow.Modflow.load(f'{model_name}.nam', version='mf2005', exe_name='mf2005.exe')



# 定义参数范围和分布

param_distributions = {

    'sy': np.random.uniform(0.1, 0.9, 1000),

    'k': np.random.uniform(10, 100, 1000)

}



# 存储结果

results = []



# 进行蒙特卡洛模拟

for sy, k in zip(param_distributions['sy'], param_distributions['k']):

    # 修改储水系数和导水率

    mf.lpf.sy.array = sy

    mf.lpf.k.array = k

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 记录结果

    results.append(np.mean(head))



# 统计分析结果

mean_head = np.mean(results)

std_head = np.std(results)



# 输出结果

print(f'平均地下水位: {mean_head:.2f}')

print(f'地下水位标准差: {std_head:.2f}')

贝叶斯方法

贝叶斯方法是一种基于贝叶斯定理的不确定性评估方法，通过先验分布和后验分布来评估参数的不确定性。

实现步骤

定义参数的先验分布
生成参数值
运行模型
记录输出变量的变化
计算后验分布

代码示例

假设我们要使用贝叶斯方法评估储水系数和导水率的不确定性。以下是一个Python脚本示例，使用Flopy库和PyMC3库来实现这一过程：


import flopy

import numpy as np

import pandas as pd

import pymc3 as pm



# 定义模型参数

model_name = 'example_model'

mf = flopy.modflow.Modflow.load(f'{model_name}.nam', version='mf2005', exe_name='mf2005.exe')



# 定义参数的先验分布

with pm.Model() as model:

    sy = pm.Uniform('sy', lower=0.1, upper=0.9)

    k = pm.Uniform('k', lower=10, upper=100)

    

    # 定义模型输出变量

    head = pm.Normal('head', mu=0, sigma=1, observed=np.array([0.0]))  # 假设观测数据为0.0



# 生成参数值

with model:

    trace = pm.sample(1000, tune=1000, chains=2)



# 存储结果

results = []



# 进行贝叶斯模拟

for sy_val, k_val in zip(trace['sy'], trace['k']):

    # 修改储水系数和导水率

    mf.lpf.sy.array = sy_val

    mf.lpf.k.array = k_val

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 记录结果

    results.append(np.mean(head))



# 输出后验分布

pm.traceplot(trace)

pm.summary(trace)

结合灵敏度分析和不确定性评估

重要性分析

结合灵敏度分析和不确定性评估，可以通过重要性分析（Importance Analysis）来确定哪些参数对模型输出的不确定性贡献最大。这有助于模型的进一步优化和校准。

代码示例

以下是一个Python脚本示例，结合局部灵敏度分析和蒙特卡洛模拟来进行重要性分析：


import flopy

import numpy as np

from SALib.sample import saltelli

from SALib.analyze import sobol



# 定义模型参数

model_name = 'example_model'

mf = flopy.modflow.Modflow.load(f'{model_name}.nam', version='mf2005', exe_name='mf2005.exe')



# 定义参数范围

problem = {

    'num_vars': 2,

    'names': ['sy', 'k'],

    'bounds': [[0.1, 0.9], [10, 100]]

}



# 生成参数组合

param_values = saltelli.sample(problem, 1024)



# 存储结果

results = []



# 进行全局灵敏度分析

for params in param_values:

    # 修改储水系数和导水率

    mf.lpf.sy.array = params[0]

    mf.lpf.k.array = params[1]

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 记录结果

    results.append(np.mean(head))



# 计算灵敏度指标

Si = sobol.analyze(problem, np.array(results))



# 输出灵敏度指标

print('S1 (First-order indices):', Si['S1'])

print('ST (Total-order indices):', Si['ST'])



# 进行蒙特卡洛模拟

param_distributions = {

    'sy': np.random.uniform(0.1, 0.9, 1000),

    'k': np.random.uniform(10, 100, 1000)

}



# 存储结果

mc_results = []



# 进行蒙特卡洛模拟

for sy, k in zip(param_distributions['sy'], param_distributions['k']):

    # 修改储水系数和导水率

    mf.lpf.sy.array = sy

    mf.lpf.k.array = k

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 记录结果

    mc_results.append(np.mean(head))



# 统计分析结果

mean_head = np.mean(mc_results)

std_head = np.std(mc_results)



# 输出结果

print(f'平均地下水位: {mean_head:.2f}')

print(f'地下水位标准差: {std_head:.2f}')



# 重要性分析

import matplotlib.pyplot as plt



# 绘制参数与输出变量的关系图

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

axs[0].scatter(param_distributions['sy'], mc_results)

axs[0].set_xlabel('储水系数')

axs[0].set_ylabel('地下水位')

axs[1].scatter(param_distributions['k'], mc_results)

axs[1].set_xlabel('导水率')

axs[1].set_ylabel('地下水位')

plt.show()

模型校准

结合灵敏度分析和不确定性评估的结果，可以进行模型校准，以提高模型的预测精度。模型校准通常涉及参数优化和模型验证。

代码示例

以下是一个Python脚本示例，使用Flopy库和Scipy库进行模型校准：


import flopy

import numpy as np

from scipy.optimize import minimize



# 定义模型参数

model_name = 'example_model'

mf = flopy.modflow.Modflow.load(f'{model_name}.nam', version='mf2005', exe_name='mf2005.exe')



# 定义观测数据

observed_head = np.array([20.0, 25.0, 30.0, 35.0])



# 定义目标函数

def objective_function(params):

    mf.lpf.sy.array = params[0]

    mf.lpf.k.array = params[1]

    

    # 保存模型文件

    mf.write_input()

    

    # 运行模型

    success, buff = mf.run_model(silent=True)

    

    # 读取输出变量

    head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

    

    # 计算残差平方和

    residuals = observed_head - np.mean(head)

    return np.sum(residuals**2)



# 初始参数值

initial_params = [0.5, 50.0]



# 进行优化

result = minimize(objective_function, initial_params, method='BFGS')



# 输出优化结果

print('优化后的储水系数:', result.x[0])

print('优化后的导水率:', result.x[1])



# 验证模型

mf.lpf.sy.array = result.x[0]

mf.lpf.k.array = result.x[1]

mf.write_input()

success, buff = mf.run_model(silent=True)

head = flopy.utils.HeadFile(f'{model_name}.hds').get_data()

print('校准后的平均地下水位:', np.mean(head))