1. 从田间到模型:为什么作物模型参数优化是个绕不开的坎
搞作物模型的人都有一个共识:模型框架本身其实大同小异,真正决定模拟精度的,是参数。WOFOST、DSSAT、APSIM这些主流模型,结构上都是把光合、呼吸、蒸腾、分配、衰老这些过程用数学公式串起来,但参数一旦不对,模拟结果和实测值能差出天际。我自己刚开始接触WOFOST的时候,觉得只要把气象数据、土壤数据、管理数据填进去,模型就能跑出合理的结果。实际一跑才发现,出苗日期差三天,最终产量模拟能差20%以上。这就是参数敏感性的问题。
作物模型参数优化,说白了就是用各种算法去“校准”模型里的那些不确定参数,让模拟值尽量逼近实测值。这些参数包括比叶面积、光能利用率、干物质分配系数、生育期积温参数、根系吸水参数等等。它们有的可以直接测量,比如积温参数可以通过生育期观测数据算出来;但更多的参数很难直接测,比如干物质向不同器官的分配系数,你不可能把一株作物拆开去称每一部分的动态变化。所以必须通过优化算法来反推。
为什么这件事现在越来越热?因为随着遥感技术和高光谱设备的普及,我们有了大量中间状态量的观测数据——叶面积指数、生物量、叶绿素含量、冠层覆盖度等等。这些数据可以作为优化过程中的约束条件,让参数校准不再只依赖最终产量这一个端点。这就把“参数优化”和“状态量光谱检测”两个方向绑在了一起。你光谱反演出来的叶面积指数越准,参数优化的目标函数就越可靠,模拟结果就越接近真实。
这个方向适合谁?如果你是做农业遥感的研究生,或者在做作物生长模拟的科研人员,又或者是农业技术推广中需要用到模型预测的从业者,那这套东西你迟早要碰。哪怕你只是用模型做灌溉决策或者产量预报,参数没调好,后面的一切都是空中楼阁。我见过太多人拿着默认参数跑模型,然后抱怨模型不准——问题不在模型,在你没做参数优化。
2. 核心概念拆解:WOFOST、光谱检测与干物质分配到底怎么串起来
2.1 WOFOST模型的基本骨架与参数体系
WOFOST是荷兰瓦赫宁根大学开发的一套作物生长模拟模型,全称是WOrld FOod STudies。它的核心思路是把作物生长分成两个层次:潜在生产水平和水分限制水平。潜在生产水平下,作物生长只受温度、辐射和作物自身遗传特性控制;水分限制水平下,还要考虑土壤水分供应对蒸腾和光合的影响。这个模型最大的优点是结构清晰、参数可解释性强,不像有些机器学习模型那样是个黑箱。
WOFOST里的参数大致可以分成几类。第一类是发育参数,包括各生育阶段的有效积温、光周期敏感性系数、春化需求等。第二类是光合与呼吸参数,比如最大光合速率、光能利用效率、维持呼吸系数。第三类是分配参数,也就是干物质向根、茎、叶、穗各器官的分配比例随发育阶段的变化。第四类是形态参数,比如比叶面积、叶面积指数最大值、根深增长速率。第五类是水分相关参数,比如根系吸水速率、土壤水分传导率等。
这些参数里,发育参数和形态参数相对容易获取,因为可以通过田间观测直接标定。但分配参数和光合参数就很难直接测。举个例子,干物质分配系数在WOFOST里是一组随发育阶段变化的曲线,出苗到开花期间分配系数逐渐从根向叶和茎转移,开花后主要向穗分配。这条曲线你怎么测?你只能通过破坏性取样,在不同生育期分别称量根、茎、叶、穗的干重,然后反推分配比例。但破坏性取样本身就有误差,而且样本量有限,不可能每天都取。所以分配参数的不确定性很大,必须通过优化来降低。
2.2 光谱检测如何提供状态量约束
光谱检测在作物模型里的角色,是提供高频次、无损的状态量观测。传统方法测叶面积指数,你得用冠层分析仪或者破坏性取样,费时费力,而且空间代表性有限。高光谱或多光谱传感器可以快速获取冠层反射率,然后通过植被指数或者辐射传输模型反演出叶面积指数、叶绿素含量、生物量等状态量。这些状态量正好是WOFOST模型每天都会模拟的输出变量。
这就形成了一个闭环:模型模拟出叶面积指数,光谱反演出叶面积指数,两者一对比,差异就是优化算法的目标。你可以用优化算法不断调整参数,让模拟的叶面积指数序列和光谱反演的叶面积指数序列最接近。这比只用最终产量做优化要稳健得多,因为产量只是一个端点,而叶面积指数是一个时间序列,包含了整个生育期的动态信息。
但这里有个坑:光谱反演本身也有误差。植被指数饱和、土壤背景干扰、冠层结构变化、光照条件差异,都会影响反演精度。如果你直接把光谱反演的叶面积指数当作真值去优化参数,那反演误差就会传递到参数里。所以实际操作中,通常会给光谱反演值一个不确定性范围,或者用滤波方法平滑时间序列,再或者把光谱反演和模型模拟做数据同化,而不是简单的点对点匹配。
2.3 干物质分配研究的核心问题
干物质分配是作物模型里最敏感也最不确定的部分之一。它决定了光合产物有多少分配到叶片用于继续生长,有多少分配到茎秆用于支撑,有多少分配到籽粒形成产量。分配系数的微小变化,经过整个生育期的累积,会导致产量模拟的巨大差异。
目前干物质分配的研究主要有几个方向。一是基于生理机制的分配模型,比如根据源库关系、激素调控、糖信号传导来模拟分配。二是基于经验统计的分配模型,比如WOFOST里用的那种随发育阶段变化的分配系数表。三是基于环境驱动的分配模型,考虑水分、氮素、温度对分配的影响。这几种方向各有优劣,机制模型解释性强但参数多且难标定,经验模型简单但外推能力差,环境驱动模型介于两者之间。
从参数优化的角度看,干物质分配参数是最需要优化的对象。因为其他参数比如积温、比叶面积,你还能通过观测直接算出来,但分配系数真的只能靠反推。而且分配系数不是常数,它随发育阶段变化,所以你要么优化一组分段常数,要么优化一条曲线的形状参数。这就把优化问题的维度提高了,对优化算法的要求也更高。
3. 参数优化实操:从数据准备到算法选型的完整流程
3.1 数据准备:气象、土壤、管理、观测四件套
参数优化的第一步是数据准备。没有数据,再好的算法也跑不起来。你需要四类数据:气象数据、土壤数据、管理数据和观测数据。
气象数据包括逐日太阳辐射、最高气温、最低气温、降水量、风速、相对湿度等。这些数据可以从气象站获取,也可以用再分析资料插值。注意,WOFOST对辐射和温度最敏感,这两个变量的精度直接影响模拟结果。如果你用的气象数据本身就有偏差,那参数优化就是在错误的基础上调参,越调越偏。
土壤数据包括土壤质地、容重、田间持水量、凋萎含水量、饱和含水量、饱和导水率等。这些参数决定了水分限制水平下的模拟结果。如果你只做潜在生产水平的模拟,土壤数据可以简化,但如果你要考虑水分胁迫,土壤参数必须准确。
管理数据包括播种日期、播种深度、种植密度、灌溉日期和灌溉量、施肥日期和施肥量等。这些数据通常来自田间记录,相对容易获取。但要注意,WOFOST对播种日期非常敏感,播种日期差几天,整个生育期模拟都会偏移。
观测数据是参数优化的目标。你需要至少一个生育期内多次观测的叶面积指数、生物量、产量等数据。观测频次越高,优化结果越可靠。如果只有最终产量一个观测点,那优化问题就是欠定的,很多组参数都能得到相似的产量,你无法判断哪组参数更合理。
3.2 敏感性分析:先找出哪些参数值得优化
不是所有参数都需要优化。有些参数对模拟结果影响很小,优化它们只会增加计算负担,甚至导致过拟合。所以第一步应该做敏感性分析,找出对目标变量影响最大的参数。
敏感性分析的方法有很多,最简单的是局部敏感性分析,也就是每次只改变一个参数,看模拟结果的变化幅度。这种方法计算量小,但忽略了参数之间的交互作用。更全面的是全局敏感性分析,比如Sobol方法或Morris筛选法,可以同时考虑多个参数的变化及其交互效应。
我在实际项目中的做法是:先用Morris筛选法快速筛一遍,找出敏感性高的参数子集,然后再用Sobol方法对子集做定量分析。这样既节省计算时间,又能得到可靠的敏感性排序。根据我的经验,WOFOST里最敏感的参数通常是:出苗到开花的有效积温、开花到成熟的有效积温、最大光合速率、比叶面积、干物质分配系数。这几个参数优化好了,模拟精度能提升一大截。
3.3 优化算法选型:遗传算法、粒子群与贝叶斯优化的取舍
参数优化本质上是一个高维非凸优化问题。目标函数通常是模拟值与观测值的均方根误差或纳什效率系数。由于WOFOST模型是非线性的,目标函数往往有多个局部最优解,传统的梯度下降法很容易陷入局部最优。所以实践中常用启发式算法。
遗传算法是用的最多的。它通过模拟自然选择和遗传变异来搜索最优解,全局搜索能力强,不容易陷入局部最优。但遗传算法的计算量很大,因为每一代都要评估大量个体,而每次评估都要跑一遍WOFOST模型。如果你的参数维度是10维,种群大小是100,迭代100代,那就是10000次模型运行。每次运行假设1秒钟,那就是将近3个小时。这还只是一个站点、一个品种的优化。
粒子群算法和遗传算法类似,也是全局优化算法,但收敛速度通常更快。它通过模拟鸟群觅食行为来更新粒子位置,参数少、实现简单。但粒子群算法容易早熟收敛,特别是在高维问题上。我的经验是,如果参数维度低于8维,粒子群算法效果不错;如果维度更高,遗传算法更稳健。
贝叶斯优化是近年来比较热门的方法,特别适合目标函数评估代价高的场景。它通过构建目标函数的概率代理模型,用采集函数来决定下一个采样点,从而用最少的模型运行次数找到最优解。如果你跑一次WOFOST要好几秒甚至更久,贝叶斯优化能帮你省下大量计算时间。但贝叶斯优化的实现比遗传算法复杂,需要额外的库支持,比如Python的scikit-optimize或GPyOpt。
3.4 目标函数设计:单一产量还是时间序列
目标函数的设计直接决定了优化结果的质量。最简单的目标函数是只用最终产量,也就是让模拟产量和实测产量的差异最小。但这样做的问题很明显:产量只是一个点,很多组参数都能得到相近的产量,但它们的叶面积指数动态、生物量累积过程可能完全不同。这样的参数外推到其他年份或其他站点,很可能失效。
更好的做法是用时间序列作为目标函数。把叶面积指数、地上生物量、各器官生物量的观测序列都纳入目标函数,让模拟值和观测值的整体差异最小。这样优化出来的参数不仅拟合了最终产量,还拟合了整个生长过程,外推能力更强。
但时间序列目标函数也有问题:不同变量的量纲不同,叶面积指数的范围是0到8,生物量的范围可能是0到20000公斤每公顷,直接相加会导致生物量主导目标函数。所以需要归一化或者加权。常用的做法是计算每个变量的相对误差,然后取平均。或者用多目标优化算法,比如NSGA-II,同时优化多个目标,最后从帕累托前沿里选一组折中解。
4. 光谱检测与数据同化:让模型和遥感真正对话
4.1 光谱反演叶面积指数的实操要点
光谱反演叶面积指数的核心是建立冠层反射率和叶面积指数之间的统计关系或物理关系。统计关系就是直接用植被指数做回归,比如归一化植被指数、增强型植被指数、比值植被指数等。物理关系则是用辐射传输模型,比如PROSAIL模型,模拟不同叶面积指数下的冠层反射率,然后通过查找表或优化算法反演。
统计方法简单快速,但外推能力差。你在一个站点建立的回归关系,换到另一个站点可能就不适用了,因为土壤背景、冠层结构、叶片生化参数都变了。物理方法机理明确,外推能力好,但计算量大,而且需要先验知识来约束反演过程。
我在实际操作中的经验是:如果只是做同一站点的时序监测,统计方法够用,但要注意植被指数饱和问题。当叶面积指数超过4到5时,归一化植被指数就饱和了,反演精度急剧下降。这时候可以用红边植被指数或者短波红外波段来补充信息。如果要做跨站点或跨年份的应用,最好用物理方法,或者至少用统计方法加上站点特异性校正。
还有一个容易被忽略的问题:光谱观测的角度和太阳角度。冠层反射率是方向性的,不同观测角度得到的反射率不同。如果你用的光谱数据没有经过角度校正,反演的叶面积指数就会有系统性偏差。所以尽量用多角度观测数据,或者至少记录观测时的太阳天顶角和观测天顶角,在反演时加以考虑。
4.2 数据同化:把光谱观测融入模型运行
数据同化是把观测数据融入模型状态估计的过程。在作物模型里,常用的同化方法有卡尔曼滤波、粒子滤波和变分同化。卡尔曼滤波适合线性系统,但WOFOST是非线性的,所以扩展卡尔曼滤波或者集合卡尔曼滤波更常用。粒子滤波适合强非线性系统,但计算量大。变分同化则是把同化问题转化为优化问题,通过调整控制变量来最小化观测和模拟的差异。
数据同化和参数优化有什么区别?参数优化是调整模型参数,让模型长期平均行为更接近观测;数据同化是调整模型状态,让模型在特定时刻更接近观测。两者可以结合:先用参数优化校准模型参数,再用数据同化在运行过程中不断更新状态变量。这样既能保证模型整体行为合理,又能吸收实时观测信息。
实际操作中,数据同化的难点在于观测误差和模型误差的估计。观测误差可以从光谱反演的验证数据里估计,但模型误差很难量化。如果观测误差设得太小,同化结果会过度依赖观测,模型自身的物理约束被削弱;如果观测误差设得太大,同化就没什么效果。我的经验是,观测误差设为反演均方根误差的1到2倍,模型误差设为观测误差的5到10倍,然后通过试错调整。
4.3 时间序列滤波:处理光谱数据的不连续性
光谱观测往往不是每天都有,云雨天气会导致数据缺失,不同卫星过境时间也不一样。所以光谱反演出来的叶面积指数时间序列通常是稀疏且不规则的。直接拿这样的数据去优化参数或同化,效果不会好。
时间序列滤波是解决这个问题的常用手段。你可以用Savitzky-Golay滤波器平滑噪声,用双逻辑斯蒂函数拟合整个生育期的叶面积指数动态,或者用高斯过程回归插补缺失值。这些方法各有优劣:Savitzky-Golay简单但需要等间隔数据;双逻辑斯蒂函数符合作物生长规律但假设了固定的曲线形状;高斯过程回归灵活但计算量大。
我比较推荐双逻辑斯蒂函数拟合,因为作物叶面积指数的动态确实符合S型曲线:出苗后缓慢增长,然后快速上升,达到峰值后缓慢下降。用双逻辑斯蒂函数拟合,既能平滑噪声,又能插补缺失值,还能提取关键物候参数比如出苗期、抽穗期、成熟期。这些物候参数反过来又可以用来约束WOFOST的发育参数,形成良性循环。
5. 干物质分配的参数化:从经验系数到环境驱动
5.1 WOFOST中的分配系数表及其局限
WOFOST里干物质分配是用一张表来描述的:每个发育阶段对应一组分配系数,分别给根、茎、叶、穗。这张表通常是基于文献或田间试验得到的,不同品种、不同环境下的分配系数表可能不同。模型运行时,根据当前发育阶段查表得到分配系数,然后把每日新增干物质按比例分配到各器官。
这张表的局限性很明显。首先,它是经验性的,没有考虑环境对分配的影响。实际上,水分胁迫、氮素胁迫、温度胁迫都会改变分配模式。比如干旱条件下,作物会把更多干物质分配到根系以增强吸水能力,分配到地上部分的就少了。但WOFOST的分配系数表是固定的,不会随环境变化。其次,这张表是离散的,发育阶段之间的分配系数是线性插值,但实际分配系数的变化可能不是线性的。第三,这张表通常是针对特定品种的,换品种就要换表,但很多品种根本没有可用的分配系数表。
5.2 分配系数的优化策略:分段常数还是连续曲线
在参数优化中,分配系数怎么参数化是个关键问题。最简单的方式是分段常数:把整个生育期分成几个阶段,每个阶段有一组固定的分配系数。比如出苗到开花、开花到灌浆、灌浆到成熟三个阶段,每个阶段优化四个分配系数,总共12个参数。这种方式参数少,优化起来相对容易,但不够灵活,可能无法捕捉分配系数的连续变化。
更灵活的方式是用连续曲线描述分配系数随发育阶段的变化。比如用多项式、样条函数或者逻辑斯蒂函数来拟合分配系数曲线。这样参数少的时候可以是一条直线,参数多的时候可以是一条复杂曲线。但参数多了,优化难度也上去了,而且容易过拟合。
我的建议是:如果你有高频次的器官生物量观测数据,可以用连续曲线;如果观测数据稀疏,分段常数更稳健。另外,你可以先用分段常数优化,看看分配系数的大致变化趋势,然后再决定是否用连续曲线细化。不要一上来就搞复杂模型,简单模型往往更可靠。
5.3 环境驱动分配模型:引入水分和氮素胁迫因子
如果你想让分配模型更机理化,可以引入环境胁迫因子来调整分配系数。比如定义一个水分胁迫因子,当土壤水分低于临界值时,增加向根系的分配比例,减少向叶片的分配比例。类似地,氮素胁迫因子可以调整向叶片的分配,因为叶片是氮素的主要需求器官。
这种环境驱动分配模型的优点是外推能力强,因为它考虑了环境变化对分配的影响。但缺点是参数更多,而且胁迫因子和分配系数之间的函数关系需要假设。你可以用线性函数、指数函数或者分段函数来描述这种关系,但每种假设都会引入不确定性。
实际操作中,我建议先用固定分配系数跑一遍,看看模拟值和观测值的差异有多大。如果差异主要在水分胁迫年份,那就说明分配系数需要环境驱动。如果差异在所有年份都差不多,那可能是其他参数的问题,不一定是分配系数。不要一上来就怀疑分配模型,先排查其他更简单的因素。
6. 常见问题与排查技巧实录
6.1 参数优化结果不稳定怎么办
参数优化结果不稳定是常见问题。同样的数据、同样的算法,跑两次得到的最优参数不一样。这通常是因为目标函数有多个局部最优解,而优化算法的随机性导致每次收敛到不同的局部最优。
解决办法有几个。一是增加种群大小或迭代次数,让算法有更多机会探索全局最优。二是多次运行优化算法,取目标函数最优的那次结果。三是用不同的优化算法交叉验证,如果不同算法都收敛到相似的参数值,那结果就比较可靠。四是增加观测数据,特别是时间序列数据,让目标函数更明确,减少多解性。
我自己的做法是:先用遗传算法跑5次,每次种群200、迭代200代,然后取最优解。如果5次结果差异很大,说明目标函数多解性严重,需要增加观测数据或者减少优化参数。如果5次结果接近,那就以最优解为准,再用粒子群算法局部搜索一下,看能不能进一步改进。
6.2 光谱反演叶面积指数偏高或偏低怎么调
光谱反演叶面积指数出现系统性偏高或偏低,通常有几个原因。一是植被指数和叶面积指数的回归关系是用特定数据集建立的,换到你的站点就不适用了。这时候需要用你的站点数据重新标定回归系数。二是土壤背景反射率的影响。如果土壤颜色浅、反射率高,植被指数会被拉低,反演的叶面积指数偏低。这时候可以用土壤调整植被指数或者修改土壤线参数。三是冠层结构的影响。如果作物行距大、冠层封闭度低,冠层反射率会包含大量土壤信息,反演精度下降。这时候可以用多角度观测或者方向性反射率模型来校正。
我的经验是:先用地面实测叶面积指数验证光谱反演结果,看看偏差有多大、是系统性偏差还是随机偏差。如果是系统性偏差,就重新标定回归模型或者调整辐射传输模型的参数。如果是随机偏差,就增加观测频次或者用时间序列滤波来平滑。
6.3 干物质分配模拟和观测对不上怎么排查
干物质分配模拟和观测对不上,可能的原因很多。首先要检查发育阶段模拟是否准确。如果发育阶段偏移了,分配系数查表就会查错,分配结果自然不对。所以先看开花期、成熟期的模拟值和观测值差多少,如果差太多,先优化发育参数。
如果发育阶段没问题,那就看总生物量模拟是否准确。如果总生物量偏高或偏低,分配比例再对也没用。总生物量取决于光合和呼吸参数,所以先优化这些参数。
如果总生物量和发育阶段都对了,但分配比例不对,那就是分配系数的问题。这时候可以固定其他参数,只优化分配系数。如果优化后还是对不上,那可能是分配模型的结构问题,比如WOFOST的固定分配系数表无法捕捉环境对分配的影响,需要考虑环境驱动分配模型。
6.4 常见问题速查表
| 问题现象 | 可能原因 | 排查方法 | 解决思路 |
|---|---|---|---|
| 参数优化结果不稳定 | 目标函数多解 | 多次运行算法 | 增加观测数据、减少参数维度 |
| 叶面积指数反演偏高 | 土壤背景干扰 | 检查土壤反射率 | 用土壤调整植被指数 |
| 叶面积指数反演偏低 | 植被指数饱和 | 检查高叶面积指数时的植被指数 | 用红边或短波红外波段 |
| 干物质分配模拟偏差大 | 发育阶段偏移 | 对比开花期和成熟期 | 先优化发育参数 |
| 总生物量模拟偏差大 | 光合参数不准 | 对比生物量累积曲线 | 优化光合和呼吸参数 |
| 产量模拟偏差大 | 分配系数不准 | 对比各器官生物量 | 优化分配系数 |
| 水分胁迫年份模拟差 | 分配未考虑环境 | 对比干旱和正常年份 | 引入环境驱动分配模型 |
7. 实操案例:一个冬小麦参数优化的完整记录
7.1 数据概况与预处理
我拿一个冬小麦站点的数据来做案例。气象数据是2018到2020年三个生长季的逐日数据,来自附近气象站。土壤数据是站点实测的,质地是壤土,田间持水量0.32,凋萎含水量0.12。管理数据包括播种日期、播种量、灌溉和施肥记录。观测数据包括每两周一次的叶面积指数、地上生物量、各器官生物量,以及最终产量。
预处理主要是气象数据的质量控制和缺失值插补。太阳辐射用日照时数换算,最高最低气温检查异常值,降水数据做累积曲线检查。土壤数据用实测值,没有做特殊处理。观测数据做了平滑处理,叶面积指数用双逻辑斯蒂函数拟合,生物量用线性插值补齐。
7.2 敏感性分析结果
用Morris筛选法对WOFOST的20个参数做敏感性分析,目标变量是叶面积指数和产量。结果显示,对叶面积指数最敏感的参数是比叶面积、出苗到开花的有效积温、最大光合速率。对产量最敏感的参数是开花到成熟的有效积温、干物质分配系数、最大光合速率。综合来看,我选了8个参数做优化:出苗到开花积温、开花到成熟积温、比叶面积、最大光合速率、初始分配系数、开花期分配系数、灌浆期分配系数、成熟期分配系数。
7.3 优化过程与结果
优化算法用遗传算法,种群大小200,迭代300代,交叉概率0.8,变异概率0.05。目标函数是叶面积指数、地上生物量和产量的归一化均方根误差的加权平均,权重分别是0.4、0.3、0.3。优化跑了5次,取目标函数最小的那次结果。
优化后的参数和默认参数相比,出苗到开花积温降低了约50度日,开花到成熟积温增加了约30度日,比叶面积增加了约15%,最大光合速率降低了约10%。分配系数方面,开花期向穗的分配比例增加了,灌浆期向穗的分配比例也增加了,说明默认参数低估了穗的库强度。
优化后的模拟结果:叶面积指数的均方根误差从1.2降到0.6,地上生物量的均方根误差从1500公斤每公顷降到800公斤每公顷,产量的相对误差从12%降到5%以内。三个生长季的模拟精度都有提升,说明优化参数没有过拟合。
7.4 光谱数据同化的效果
在参数优化的基础上,我又做了光谱数据同化。用 Sentinel-2 影像反演的叶面积指数,通过集合卡尔曼滤波同化到WOFOST里。同化后,叶面积指数的模拟精度进一步提升,均方根误差降到0.4。产量模拟精度也有改善,相对误差降到3%以内。特别是在水分胁迫年份,同化的效果更明显,因为光谱观测能捕捉到水分胁迫导致的叶面积指数下降,模型通过同化及时调整了状态变量。
8. 工具链与代码实现要点
8.1 Python生态下的WOFOST调用
WOFOST有Python接口,叫PCSE,全称Python Crop Simulation Environment。安装很简单,pip install pcse就行。PCSE提供了WOFOST的完整实现,可以方便地修改参数、运行模型、提取输出。
调用WOFOST的基本流程是:先定义气象数据提供器、土壤数据提供器、作物参数提供器,然后创建WOFOST模型对象,最后运行模型并提取输出。下面是一个简化的代码示例:
import pcse from pcse.models import Wofost72_WLP_FD from pcse.fileinput import YAMLCropDataProvider from pcse.db import NASAPowerWeatherDataProvider # 加载作物参数 cropdata = YAMLCropDataProvider() cropdata.set_active_crop('wheat', 'Winter_wheat_101') # 加载气象数据 weather = NASAPowerWeatherDataProvider(latitude=35, longitude=116) # 定义土壤参数 soil = { 'SMFCF': 0.32, 'SMW': 0.12, 'SM0': 0.45, 'CRAIRC': 0.06, 'K0': 10.0, 'SOPE': 1.5, 'KSUB': 10.0, 'RDMSOL': 120.0 } # 定义管理参数 sowing_date = '2018-10-15' harvest_date = '2019-06-10' # 创建模型 model = Wofost72_WLP_FD( parameter_providers=(cropdata, soil), weather_data=weather, agromanagement={ '2018-10-15': { 'CropCalendar': { 'crop_name': 'wheat', 'variety_name': 'Winter_wheat_101', 'crop_start_date': sowing_date, 'crop_end_date': harvest_date } } } ) # 运行模型 model.run_till_terminate() output = model.get_output()这段代码跑完,output里就有逐日的叶面积指数、生物量、产量等变量。你可以把output转成Pandas DataFrame,方便后续分析和优化。
8.2 优化算法的代码框架
优化算法的代码框架可以这样写:先定义一个目标函数,输入是参数向量,输出是模拟值和观测值的差异。然后在目标函数里调用WOFOST模型,提取模拟值,计算差异。最后用遗传算法或粒子群算法来最小化目标函数。
import numpy as np from scipy.optimize import differential_evolution def objective(params): # 把参数向量映射到WOFOST参数 cropdata.set_override('TSUM1', params[0]) cropdata.set_override('TSUM2', params[1]) cropdata.set_override('SLATB', params[2]) cropdata.set_override('AMAXTB', params[3]) # ... 其他参数 # 运行模型 model = Wofost72_WLP_FD(...) model.run_till_terminate() output = model.get_output() # 提取模拟值 sim_lai = [o['LAI'] for o in output] sim_biomass = [o['TAGP'] for o in output] sim_yield = output[-1]['TWSO'] # 计算目标函数 rmse_lai = np.sqrt(np.mean((np.array(sim_lai) - obs_lai)**2)) rmse_biomass = np.sqrt(np.mean((np.array(sim_biomass) - obs_biomass)**2)) rmse_yield = np.sqrt((sim_yield - obs_yield)**2) return 0.4*rmse_lai + 0.3*rmse_biomass + 0.3*rmse_yield # 定义参数边界 bounds = [(800, 1200), (600, 1000), (0.002, 0.004), (20, 40)] # 运行差分进化算法 result = differential_evolution(objective, bounds, maxiter=300, popsize=20) print(result.x)这个框架用的是差分进化算法,和遗传算法类似但实现更简单。你可以根据实际情况换成其他算法。
8.3 光谱反演与同化的代码要点
光谱反演叶面积指数可以用Python的spectral库或者自己写植被指数计算。同化可以用集合卡尔曼滤波,Python里有filterpy库可以用。
import numpy as np from filterpy.kalman import EnsembleKalmanFilter # 定义状态转移函数(WOFOST模型) def state_transition(state, dt): # 这里简化处理,实际要调用WOFOST return state # 定义观测函数(光谱反演) def observation_model(state): # 这里简化处理,实际要用植被指数反演 return state[0] # 假设状态第一个变量是叶面积指数 # 初始化集合卡尔曼滤波 enkf = EnsembleKalmanFilter( x=np.array([0.5]), # 初始状态 P=np.array([[0.1]]), # 初始协方差 dim_z=1, # 观测维度 dt=1.0, N=50 # 集合大小 ) # 同化循环 for t in range(len(obs_lai)): enkf.predict() if not np.isnan(obs_lai[t]): enkf.update(np.array([obs_lai[t]])) print(f"Time {t}: LAI = {enkf.x[0]:.2f}")这段代码只是示意,实际同化需要把WOFOST的状态变量和观测变量对应起来,还要处理模型误差和观测误差的协方差矩阵。
9. 我踩过的坑与实操心得
第一个坑是气象数据的辐射单位。WOFOST默认的辐射单位是焦耳每平方厘米每天,但很多气象站给的是兆焦耳每平方米每天。这两个单位差100倍,如果不换算,模拟结果会完全不对。我刚开始跑的时候,叶面积指数模拟值只有0.1,查了半天才发现是辐射单位错了。
第二个坑是土壤水分参数的标定。WOFOST里的田间持水量和凋萎含水量不是直接测量的值,而是模型特定的参数。如果你直接用实测的田间持水量,模拟结果可能不对。最好用模型自带的土壤参数数据库,或者通过优化来标定。
第三个坑是分配系数的优化顺序。我一开始把所有参数一起优化,结果分配系数优化得乱七八糟,因为其他参数不准的时候,分配系数会被错误地补偿。后来我改成先优化发育参数和光合参数,固定这些参数后再优化分配系数,结果就稳定多了。
第四个坑是光谱反演的时间匹配。光谱影像的过境时间和模型输出的时间可能不一致,比如影像可能是上午10点,模型输出是每日平均值。如果不做时间匹配,同化效果会打折扣。我的做法是把模型输出插值到影像过境时刻,或者把影像反演值当作日平均值处理。
第五个坑是优化算法的计算时间。跑一次遗传算法优化可能要几个小时甚至几天,特别是参数维度高的时候。我的建议是先用小种群、少迭代快速试跑,看看目标函数的大致范围,然后再用大种群、多迭代精细优化。另外,可以用并行计算来加速,PCSE支持多进程运行,可以把种群评估并行化。
最后分享一个小技巧:在优化之前,先用默认参数跑一遍,看看模拟值和观测值的差异在哪里。如果差异主要在生育后期,那可能是发育参数或分配参数的问题;如果差异在整个生育期都存在,那可能是光合参数或初始参数的问题。这样有针对性地优化,比盲目调参效率高得多。