作物模型狂奔：WOFOST(PCSE) 数据同化思路

去B吧，这里没图

整体思路：PCSE -》 敏感性分析 -》调参 -》同化

0、准备工作

0.0 电脑环境

我用的Win10啦，Linux、Mac可能得自己再去微调一下。

0.1 Python IDE

我用的Pycharm，个人感觉最好使的IDE，没有之一。

Python 解释器随便装个咯，我用的 Python 3.10 版本。

对于 Python，我其实也是个小白，只会写些小脚本，哈哈哈哈哈。

1、核心依赖包

1.1 PCSE

Wofost 模型Python版本

1.2 SALib

敏感性分析用的包

1.3 SPOTPY

参数优化包，算法很多，我用的SA

2、编码思路

2.1 敏感性分析

目标变量 地上部生物量:TAGP；最大叶面积指数:LAIMAX；器官重:TWSO；
target_variable = run_details.target_variable

敏感度分析方法 efast; sobol

sa_method = run_details.sa_method

模型类别 限水:WLP; 潜在:PP

run_type = run_details.run_type

模型类别 作物:crop; 作物&水:cropAndWater

params_type = run_details.params_type

为每个标量参数确定一个合理的区间

problems_parameters = pd.read_excel(params_dir / f"sensitivity_analysis_parameters_{params_type}.xlsx")
problem = {
‘num_vars’: len(problems_parameters),
‘names’: problems_parameters[‘name’].tolist(),
‘bounds’: problems_parameters[[‘low’, ‘up’]].values.tolist()
}

随机数种子

seed = 2000

二阶

calc_second_order = True

样本数

nsamples = 256
print(“\n================ 参数抽样开始 ================\n”)
st = datetime.datetime.now().timestamp()

生成抽样参数集

paramsets = None

efast

nsamples = 65 * len(problems_parameters)
paramsets = fast_sampler.sample(problem, nsamples, seed=seed)

保存参数集

paramsets_df = pd.DataFrame(paramsets)
paramsets_df.columns = problem[‘names’]
et = datetime.datetime.now().timestamp()
print(f"\n================ 生成参数集结束 {et-st}s ================\n")

开并行狂奔

target_results = []
with tqdm(total=len(paramsets)) as pbar:
# cpu核数进程池
with mp.Pool(mp.cpu_count()) as pool:
# 并行执行
for result in pool.imap(run_wofost_partial, paramsets):
target_results.append(result)
pbar.update()

et = datetime.datetime.now().timestamp()
print(f"\n================ 执行模拟结束 {et-st}s ================\n")

敏感性分析

Si = fast.analyze(problem, target_results, seed=seed)

2.2 调参

话不多说，直接模拟退火

初始值用的参数默认值，然后上下波动优化

初始化代价函数计算器

objfunc_calculator = ObjectiveFunctionCalculator()

待优化参数边界

lowers, uppers, steps, defaults = [], [], [], []
defaults = [22.8, 0.00406, 525, 962, 24.64, 0.455]
for key, default, fluctuate in tuning_parameters:
lowers.append(default * (1 - fluctuate))
uppers.append(default * (1 + fluctuate))
steps.append(default * 2 * fluctuate / step_times)

参数初始位置

firstguess = defaults
sa = SA(func=objfunc_calculator, x0=firstguess, T_max=100, T_min=1e-7, L=1000, max_stay_counter=5, lb=lowers, ub=uppers)
best_x, best_y = sa.run()

2.3 数据同化

这里用的ENKF同化叶面积指数。

定义观测数据集变量

observations_for_DA = [(row[‘day’].date(), {‘LAI’: (row[‘LAI’], std_lai[index])}) for index, row in df_observation_LAI.iterrows()]

初始化了 WOFOST 模型的集合。

for i in range(ensemble_size):
for par, distr in override_parameters.items():
p.set_override(par, distr[i])
member = Wofost72_WLP_FD(p, weatherdataprovider, agromanagement)
ensemble.append(member)

每个可用的观测值重复同化步骤

for i in range(0, len(observations_for_DA)):
# 模型狂奔
da_enkf_single(i)

最后，我们可以使用 run_till_terminate() 调用运行 WOFOST 集合直到生长季节结束。

for member in ensemble:
member.run_till_terminate()

欢迎各位道友关注、留言、私聊、交流病情。

去B站讨论吧，平常不登录CSDN

挂个小广告不会太过份吧~

需要指导的话，那就得让我挣点零花咯，嘿嘿。PS：不说虚的，程序狂奔才是最重要的。

【闲鱼】https://m.tb.cn/h.5ttgPfa?tk=ue5dW9B3RMm HU9046 「我在闲鱼发布了【wofost python版本 pcse代码 指导】」