作物模型狂奔：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)