📚 课程目录第2章作物生长模型驱动的产量预测

2. 作物生长模型驱动的产量预测

任务说明

本案例针对特定地块，利用历史气象数据（温度、降水、日照）与土壤数据（有机质、速效氮），通过简化的作物生长模型（基于积温与水分胁迫的干物质积累模型）预测该地块的作物产量。你需要完成数据准备、模型参数设定、产量计算与结果输出，并验证模型对输入变化的响应是否符合农学常识。

数据与模型定义

模型核心假设：作物产量与生育期内有效积温累积呈线性关系，但受水分供应（降水+土壤初始含水量）的胁迫限制。土壤肥力（有机质、速效氮）作为产量修正系数。

输入数据（以冬小麦为例，生育期150天）：

气象数据：逐日最高温、最低温、降水量（单位：℃、mm）
土壤数据：土壤有效含水量（mm）、有机质含量（%）、速效氮含量（mg/kg）

模型公式：

有效积温 GDD = Σ max(0, (Tmax+Tmin)/2 - Tbase)  其中 Tbase=5℃
水分胁迫系数 W = min(1, (降水累计 + 初始土壤有效水) / 作物需水量)
作物需水量 ETc = 累计GDD × 0.25（mm/℃·d）
肥力修正 F = 0.8 + 0.1×有机质(%) + 0.002×速效氮(mg/kg)
产量 Y = 潜在产量系数 K × GDD × W × F   （K=0.8 kg/ha·℃）

代码实现

import json

def read_inputs(weather_file, soil_file):
    """读取气象与土壤JSON数据"""
    with open(weather_file, 'r', encoding='utf-8') as f:
        weather = json.load(f)
    with open(soil_file, 'r', encoding='utf-8') as f:
        soil = json.load(f)
    return weather, soil

def calculate_gdd(weather, tbase=5.0):
    """计算逐日有效积温并返回累计值"""
    total_gdd = 0.0
    for day in weather['daily']:
        tavg = (day['tmax'] + day['tmin']) / 2.0
        gdd = max(0.0, tavg - tbase)
        total_gdd += gdd
    return total_gdd

def water_stress_factor(weather, soil, total_gdd):
    """计算水分胁迫系数"""
    total_rain = sum(day['precip'] for day in weather['daily'])
    available_water = total_rain + soil['initial_soil_water']
    etc = total_gdd * 0.25  # 作物需水量
    w = min(1.0, available_water / etc) if etc > 0 else 1.0
    return w

def fertility_factor(soil):
    """计算土壤肥力修正系数"""
    return 0.8 + 0.1 * soil['organic_matter'] + 0.002 * soil['available_n']

def predict_yield(weather_file, soil_file, k=0.8):
    """主函数：预测产量(kg/ha)"""
    weather, soil = read_inputs(weather_file, soil_file)
    total_gdd = calculate_gdd(weather)
    w = water_stress_factor(weather, soil, total_gdd)
    f = fertility_factor(soil)
    yield_pred = k * total_gdd * w * f
    return {
        'total_gdd': round(total_gdd, 1),
        'water_stress': round(w, 3),
        'fertility_factor': round(f, 3),
        'yield_kg_per_ha': round(yield_pred, 1)
    }

# 示例数据文件：weather.json
weather_data = {
    "daily": [
        {"tmax": 12, "tmin": 4, "precip": 2},
        {"tmax": 15, "tmin": 6, "precip": 0},
        {"tmax": 18, "tmin": 8, "precip": 5},
        {"tmax": 20, "tmin": 10, "precip": 3},
        {"tmax": 22, "tmin": 12, "precip": 0}
    ]
}
with open('weather.json', 'w') as f:
    json.dump(weather_data, f)

# 示例数据文件：soil.json
soil_data = {
    "initial_soil_water": 40,
    "organic_matter": 2.5,
    "available_n": 80
}
with open('soil.json', 'w') as f:
    json.dump(soil_data, f)

# 执行预测
result = predict_yield('weather.json', 'soil.json')
print(json.dumps(result, ensure_ascii=False, indent=2))

运行与验证

运行方式：将上述代码保存为 yield_predict.py，在Python 3环境中执行 python yield_predict.py。代码会自动生成示例气象与土壤数据文件并运行预测。

验证方法：

检查输出结果中 total_gdd 是否等于逐日有效积温之和（示例数据中第1天(12+4)/2-5=3，第2天(15+6)/2-5=5.5，第3天(18+8)/2-5=8，第4天(20+10)/2-5=10，第5天(22+12)/2-5=12，合计38.5℃）。
验证水分胁迫：累计降水10mm + 初始土壤水40mm = 50mm，需水量=38.5×0.25=9.625mm，W=min(1, 50/9.625)=1.0，说明水分充足。
验证肥力修正：F=0.8+0.1×2.5+0.002×80=1.21。
预期产量：0.8×38.5×1.0×1.21=37.3 kg/ha（与代码输出一致）。
修改weather.json中降水为0，重新运行，观察产量是否下降（水分胁迫生效）。

判定标准：程序无报错运行，输出结果与上述手算一致；当降水减少时产量降低，证明模型对水分输入敏感。