简介WOFOST源代码压缩包面向农业科研人员、作物模型开发者及农业院校师生可支撑作物生产潜力模拟、农业管理策略评估与气候变化影响预测等研究场景。代码源自荷兰瓦赫宁根大学内部含316个文件核心类型包括wof模型配置与输入参数、forFortran程序源码、ltm/cab辅助配置与数据文件等另有少量dat、inc及bat/shell脚本整体体积仅1.51MB结构清晰便于阅读与移植。目前已有990人浏览学习在农业模型圈内受到一定关注。通过源码可深入理解气候、土壤水文、养分、作物生理、收获与管理六大模块的算法逻辑与数据结构既能调试优化运行效率也能扩展新作物或新环境参数还可结合遥感/GIS数据做敏感性分析与精度提升是学习经典作物模型机制、开展本地化二次开发的一份宝贵基础材料。1. WOFOST 源代码一套把作物生长过程变成可计算模块的经典实现WOFOSTWorld FOod STudies是瓦赫宁根大学开发的作物生长模拟模型源代码从 1980 年代延续至今Fortran 到 Python 的版本都有。很多做农业遥感估产、气候变化影响评估的工程师第一次接触它是想搞懂 CERES 之外另一套机理模型是怎么组织的。和统计模型不同WOFOST 不追求拟合产量而是从光能利用、叶片生长、生物量分配这些生理过程推导出逐日生长量——这意味着源代码里几乎每一个函数都对应一个农学概念。读这份代码的价值在于它把“作物生长”这件看似依赖经验的事拆成了可替换的模块你既可以直接跑通做区域产量模拟也可以替换其中光截获、蒸散、物候等子模块去适配自己的研究场景。本文从代码阅读和运行的角度把 WOFOST 的工程实现路径梳理一遍源码结构、关键算法的代码落点、参数文件怎么改、输出怎么验证。2. WOFOST 源码形态版本选择、仓库结构与数据流2.1 Fortran 原版与 Python 重写的取舍接触 WOFOST 时第一个选择是版本。经典 Fortran 版本如 WOFOST 7.x是模型科学的源头所有算法从这里扩散到 DSSAT、APSIM 等平台。Fortran 版的优点是计算效率高、与历史文献中的参数表完全对应比如 van Diepen 等 1989 年发表的那套参数说明直接以 Fortran 代码逻辑为参照。缺点是它保留了早期结构化编程的风格全局变量较多数据传递路径读起来比较费劲。Python 版的 pcse 库则是瓦赫宁根研究团队官方维护的重新实现。它的核心益处是模块边界清晰把 Fortran 版里的 COMMON BLOCK 拆成单独的类比如Assimilation、Evapotranspiration、Partitioning各自独立遇到问题可以直接断点到某个类里排查。对于刚上手的人我一般建议从 pcse 入手跑通整体流程同时把 Fortran 版的教科书和源码放在手边对照——因为许多参数的原始描述、默认值注释还保留在 Fortran 代码里。2.2 代码仓库里必须认识的关键文件与目录以 pcse 库为例源代码的核心结构大致如下pcse/ ├── core.py # 模型运行框架定义了 run 循环 ├── models/ │ ├── Wofost80.py # WOFOST 8.0 版本的具体组装 │ ├── Wofost71.py # 兼容 7.1 的老版本参数结构 ├── inputs.py # 天气数据、作物参数、土壤参数接入层 ├── outputs.py # 输出变量定义与格式控制 ├── util.py # 通用工具如湿度计算、日期转换 ├── soil.py # 土壤水分平衡模块 ├── assimilation.py # CO2 同化模块 ├── evapotranspiration.py# 蒸散模块 ├── phenology.py # 物候发育模块 ├── partitioning.py # 干物质分配模块 ├── leaf_dynamics.py # 叶片生长与衰老模块 ├── stem_dynamics.py # 茎生长模块 ├── storage_organ.py # 储存器官产量形成模块 ├── root_dynamics.py # 根系生长模块 ├── waterbalance.py # 水分平衡模块 ├── warnings.py # 警告处理 ├── data/ │ ├── cropdata/ # 作物参数文件cab.par, wwhe.par 等 │ ├── soildata/ # 土壤参数文件ec1.par, ec3.par 等 │ └── weatherdata/ # 示例气象数据以上是 Python 版的结构参考每个模块名称在 Fortran 版中也都能找到对应的子程序。阅读时有几条主线一是物候链从phenology.py出发追踪发育阶段怎么算二是碳链从assimilation.py出发追踪 CO2 同化产物怎么变成生物量三是水链从waterbalance.py出发追踪水分胁迫使气孔关闭的反馈路径。这三条线在models/Wofost80.py的__call__方法中汇合这就构成了模型运行的主循环。2.3 主循环的驱动方式与模块间的数据契约WOFOST 源代码的骨架是“逐日推进”的循环——每天的计算依赖前一天的状态变量所以本质上是一个离散时间动态系统。在 pcse 中最外层由时间步进器驱动从出苗/播种期开始逐日调用模型直到达到成熟期或最大物候期。# 简单示例手动跑一个生长季 from pcse.base import WeatherDataContainer from pcse.models import Wofost80 # 生成一个简化的天气数据容器 weather WeatherDataContainer( latitude31.2, longitude121.5, elevation4, day2000 # 这里仅是示意实际需逐日提供 )实际使用中不会这样逐日构造天气而是从pcse.inputs.CABOWeatherDataProvider读取包含逐日数据的文本文件。关键点在于模型的每个模块都接收并返回约定好的变量名比如DVS发育阶段TSUM积温LAI叶面积指数TAGP地上总生物量WSO储存器官重量模块之间的耦合由变量名的约定来保证这就是为什么读源码时看到某个类的方法是带self.kiosk这种字典——它是变量共享总线。理解这一点后再去追具体的算法实现思路会顺很多。3. 核心算法在源代码中的实现与参数作用位置3.1 物候发育模块积温驱动与源代码中的分阶段处理物候模型控制作物什么时期进入什么发育阶段。WOFOST 用DVS表示发育状态出苗为 0开花通常在 1成熟为 2。DVS 的推进靠积温TSUM的累计——但它不是简单的逐日温度求和而是区分营养生长和生殖生长两个阶段每个阶段有各自的积温上限。在phenology.py中核心代码大致逻辑如下class DVS_Phenology: 基于积温的发育阶段推进 def __init__(self, params): self.TSUM1 params[TSUM1] # 出苗到开花的有效积温 self.TSUM2 params[TSUM2] # 开花到成熟的积温 self.DVS 0.0 self.TSUM 0.0 def __call__(self, day, delta_t, weather): # 日平均温度 temp (weather.TMIN weather.TMAX) / 2.0 # 有效积温扣除基点温度 if temp self.TBASE: dtsum temp - self.TBASE else: dtsum 0.0 # 依据当前阶段决定推进增量 if self.DVS 1.0: # 营养生长期积温增量直接累加 self.TSUM dtsum self.DVS self.TSUM / self.TSUM1 else: # 生殖生长期只累加开花后的积温 self.TSUM dtsum self.DVS 1.0 (self.TSUM - self.TSUM1) / self.TSUM2 # 发育阶段必须限制在 0~2 之间 self.DVS min(self.DVS, 2.0) return self.DVS这段代码要注意TSUM1与TSUM2的实际取值在作物参数文件里。比如春小麦的TSUM1约在 900~1200 ℃·d 区间冬小麦则要根据春化处理来确定。很多新手调参时只改产量相关的参数实际上生育期长短由这两个值决定——改错了后续光合作用、分配全部错位。3.2 光能利用与 CO2 同化从光截获到总同化的代码链路这是 WOFOST 机理最强的部分。模型没有直接用“光能利用率”这种经验系数而是先算冠层截获的光合有效辐射再通过光合作用响应曲线求瞬时同化速率最后对一天的光周期积分。assimilation.py中的实现核心是高斯积分——因为同化速率在一天内随太阳高度角变化WOFOST 用三点高斯积分逼近日总同化量。# WOFOST 日同化量的核心逻辑简化自 assimilation.py import math # 大气顶层辐射转化为光合有效辐射的系数 F_PAR 0.5 def total_assimilation(daytime_temp, radiation, lai, params): 计算日总同化量kg CO2 / ha / day AMAX params[AMAX] # 单叶最大光合速率 EFF params[EFF] # 光能初始利用效率 KDIF params[KDIF] # 散射光消光系数 # 将总辐射换算为光合有效辐射 PAR radiation * F_PAR # 依据叶面积指数计算截获比例 if lai 0: f_int 1.0 - math.exp(-KDIF * lai) else: f_int 0.0 # 日同化总量 入射PAR × 截获比例 × 最大同化速率的函数密度依赖 assimilation PAR * f_int * AMAX return assimilation实际源码中会区分阴叶和阳叶、考虑湿度对气孔导度的影响所以远比上面复杂。但读者抓住这条主线就够用来追踪问题如果模拟产量偏低先查 LAI——如果 LAI 上不去光截获就低日同化总量自然低产量不可能高。EFF表示光响应曲线的初始斜率通常在 0.45~0.55 kg CO2 / (MJ / m²) 左右AMAX在 30~70 kg CO2 / ha / h 之间C3 和 C4 作物差异很大KDIF约 0.6 左右。调参时注意AMAX 的调整直接影响产量上限但绝不能盲目调高——它需要与叶片氮含量、温度响应曲线配套调整。3.3 干物质分配DVS 驱动的分配系数表WOFOST 将每天形成的同化产物按阶段分配到根、茎、叶、储存器官。分配系数是 DVS 的分段线性表具体数值存在作物参数文件的FR表里。partitioning.py中有一个关键表DVS_FR [ (0.00, 0.50, 0.35, 0.15, 0.00), (0.50, 0.35, 0.35, 0.30, 0.00), (1.00, 0.15, 0.30, 0.40, 0.15), (2.00, 0.00, 0.10, 0.20, 0.70), ]每行格式是(DVS, 根分配系数, 茎分配系数, 叶分配系数, 储存器官分配系数)。逐日计算时先找到当前 DVS 所在区间用线性插值得到当天的分配比例。这种分段表的特点是“刚性”——某个阶段的分配比例固定不变它假设环境不改变分配策略。实际上高温或干旱会使更多碳流向根部如果要模拟这样的胁迫效应就需要修改partitioning.py的逻辑或引入干旱系数。我一般建议初学者把主要精力放在这几个模块上phenology→assimilation→partitioning先把这三段的因果关系练熟再碰蒸散和水分平衡。4. 跑通一个最小 WOFOST 模拟并定位输出变量4.1 环境准备与最小运行示例实际操作时先用 pcse 库跑一个最小案例。需要注意的是 WOFOST 对天气数据格式有严格要求缺测值用 -999 占位日期必须是连续日序列。可以用示例数据快速验证安装成功pip install pcse # 查看自带的示例数据位置 python -c import pcse; print(pcse.__file__)然后写一个最小脚本跑完整个过程import sys from pathlib import Path import pandas as pd from pcse.inputs import CABOWeatherDataProvider from pcse.fileinput import PCSEFileReader from pcse.models import Wofost80 from pcse.base import ParameterProvider # 配置数据路径需要事先下载 WOFOST 的 crop/soil/weather 示例数据 crop_data PCSEFileReader(cropdata/wwhe.par) soil_data PCSEFileReader(soildata/ec3.par) weather_data CABOWeatherDataProvider(weatherdata/示例气象站) # 组装参数提供器 parameters ParameterProvider(cropdatacrop_data, soildatasoil_data) # 建立模型并运行 wofost Wofost80(parameters, weather_data, soil_data) wofost.run_till_terminate() # 获取输出结果 output wofost.get_output() df pd.DataFrame(output) print(df.tail())这段代码有几个容易踩的坑PCSEFileReader要求参数文件编码为 ASCII含中文注释的文件会报编码错土壤参数文件里的初始含水量字段必须和所选土壤类型匹配不匹配时模型的蒸发过程会异常run_till_terminate()在没有设置最大天数时会一直跑到成熟条件触发。4.2 输出变量的单位换算与核心检查项模型输出的原始变量单位与国际制单位不同这是初读源代码最容易误解的地方变量名含义输出单位换算成常用单位LAI叶面积指数m²/m²无需换算TAGP地上总生物量kg/ha直接读数 × 0.001 t/haWSO储存器官干重kg/ha直接读数 × 0.001 t/haTRA作物实际蒸腾cm/day× 10 mm/daySM根区土壤含水量cm³/cm³无需换算验证模型合理性时按这个顺序来先看物候——DVS是否在预期天数达到 1 和 2再看LAI峰值——一般作物在 4~6 之间过大或过小都说明光截获参数或分配系数有问题最后看WSO与TAGP的比例收获指数是否处于合理范围。4.3 用 Python 批量跑多个站点或年份实际业务中不会只跑单站。借助 Python 脚本可以批量模拟import glob import pandas as pd from pcse.models import Wofost80 from pcse.base import ParameterProvider from pcse.fileinput import PCSEFileReader, CABOWeatherDataProvider results [] weather_files glob.glob(weather/*.csv) crop_params PCSEFileReader(cropdata/wwhe.par) soil_params PCSEFileReader(soildata/ec3.par) for weather_file in weather_files: # 从文件名提取站点名称 station Path(weather_file).stem weather_data CABOWeatherDataProvider(weather_file) parameters ParameterProvider(cropdatacrop_params, soildatasoil_params) wofost Wofost80(parameters, weather_data, soil_data) wofost.run_till_terminate() output pd.DataFrame(wofost.get_output()) # 提取出苗至成熟的终点时刻数据 final output.iloc[-1] results.append({ station: station, maturity_date: final[day], WSO: final[WSO], TAGP: final[TAGP], LAI_max: output[LAI].max(), }) summary pd.DataFrame(results)这段代码比单站版本多了两层一是遍历所有天气文件二是从完整输出里提取成熟期终值。注意运行多个站点时需要为每个站点重新创建Wofost80实例——因为模型的内部状态是有记忆的复用同一个实例会让后一个站点的初始状态承接前一个站点的结束状态。4.4 参数文件的直接修改与敏感性初判WOFOST 的参数文件是文本格式直接编辑即可。以作物参数为例关键字段如下AMAX 45.0 # 单叶最大光合速率 [kg CO2/ha/hr] EFF 0.50 # 光能初始利用效率 [kg CO2/(MJ/m2)] TSUM1 1000.0 # 出苗到开花积温 [degC·d] TSUM2 800.0 # 开花到成熟积温 [degC·d] TBASE 0.0 # 发育基点温度 [degC] KDIF 0.6 # 散射辐射消光系数 RGRLAI 0.008 # 早期叶面积指数相对增长率 [1/day]做区位迁移研究时优先调TSUM1/TSUM2匹配当地观测的生育期做水分胁迫研究时优先检查soil.par的土壤水分特征曲线参数而不是调作物的光合参数。5. 从源代码延伸模型替换、加速运算与气象数据接口5.1 基于源代码替换光合模块的一个实践技巧WOFOST 的光合模块假定每日同化速率只受温度和辐射影响在极端高温38℃场景下模拟值常常偏高。常见做法是修改assimilation.py把高温胁迫引入 AMAX 的日衰减# 自定义的日同化模块示意 def calc_assimilation_with_heat_stress(TMAX, radiation, lai, params): 把极端高温对光合的抑制加入原模型 amax params[AMAX] # 高温抑制曲线35℃ 后线性衰减 if TMAX 35.0: reduction max(0.0, 1.0 - (TMAX - 35.0) * 0.05) amax * reduction # 其余逻辑复用原始代码 ...修改源代码后需要做的验证是跑一遍对照组的旧版本比较产量和 LAI 的差值确认只有高温日的同化量发生变化而不是储入了系统性偏差。5.2 气象数据接口把 NetCDF 转成 WOFOST 可读格式实际操作中常遇到气象数据是 NetCDF 格式而 WOFOST 需要逐日文本。转换时注意单位import xarray as xr import pandas as pd ds xr.open_dataset(era5_daily.nc) lat, lon 31.2, 121.5 # 选择最近的格点 ds_point ds.sel(latitudelat, longitudelon, methodnearest) # 构建 WOFOST 所需的字段 df pd.DataFrame({ DAY: ds_point.time.values, IRRAD: ds_point.ssrd.values / 1e6, # J/m2 - MJ/m2 TMIN: ds_point.t2m.values - 273.15, # K - C TMAX: ds_point.t2m.values - 273.15, VAP: ds_point.d2m.values, WIND: ds_point.wind_speed.values, # 需要单独变量或从 u/v 合成 RAIN: ds_point.tp.values * 1000, # m - mm }) df.to_csv(weather_station.csv, indexFalse)WOFOST 的天气格式中IRRAD是全天天顶辐射MJ/m²/day注意别和光合有效辐射混淆。VAP是水汽压单位 kPa——ERA5 的露点温度需要先换算成水汽压再填充。5.3 加速大规模模拟的两个实用策略面积较大的应用需要模拟成千上万个格点。WOFOST 每个格点就是一次完整的时间循环串行很慢。常见做法多进程并发multiprocessing.Pool按格点并行网格数据适合 CPU 密集场景。相关性裁剪大面积运行时很多格点气候相似可以先把气象数据聚类每类只跑一次再用空间插值恢复全区域。这种方法做省级估产时加速比可以到 20 倍以上。from multiprocessing import Pool def run_one_grid(grid_info): lat, lon grid_info weather load_weather_for_point(lat, lon) wofost Wofost80(parameters, weather, soil) wofost.run_till_terminate() # 取出终点产量 output pd.DataFrame(wofost.get_output()) return {lat: lat, lon: lon, yield: output[WSO].iloc[-1]} if __name__ __main__: grid_cells [(lat, lon) for lat in range(28, 35) for lon in range(118, 123)] with Pool(processes8) as pool: result pool.map(run_one_grid, grid_cells)注意 Windows 系统上需要把Pool调用放在if __name__ __main__:下否则进程递归执行时会报错。5.4 检验你的修改是否偏离原始模型的基准测试策略修改源码后判断是否偏离基准最好用同一套天气、参数、土壤数据跑修改前后两个版本输出逐日的 DVS、LAI、TAGP、WSO 四组变量计算均方根误差from sklearn.metrics import mean_squared_error import numpy as np rmse np.sqrt(mean_squared_error(base_df[WSO], modified_df[WSO])) print(fWSO RMSE {rmse:.2f} kg/ha)WSO 的 RMSE 在营养生长阶段应当接近 0因为此时 WSO 通常为 0到灌浆期后才无意义——因此更合理的做法是把对比区间限定在开花以后。此外检查 LAI 的峰值如果 LAI 在所有时段都系统性偏低说明光合参数修改引入了与辐射无关的偏移如果仅高温日偏离则说明胁迫逻辑起效了。代码修改的边界也要明确不要为了拟合一个站点的产量而把AMAX调到 100——这会让模型对辐射的响应完全失真。每一步代码修改都要有农学依据否则模拟结果只是数字游戏。本文还有配套的精品资源点击获取