2. 基于蒸散发的灌溉需水量日计算

任务边界与计算口径

本例要交付的是一块目标地块的逐日灌溉需水深度，单位统一为 mm/d。计算链为：气象数据 → 参考蒸散发 ET₀ → 乘以作物系数 Kc 得到作物需水量 ETc → 扣除有效降雨 Pe 得到净灌溉需水量 IRn。所有深度均指水层深度，不涉及流量与时间换算，后者属于控制器执行环节。

采用 FAO-56 推荐的 Penman-Monteith 日尺度公式作为 ET₀ 计算核心：

ET₀ = [0.408·Δ·(Rn − G) + γ·(900/(T+273))·u₂·(es − ea)] / [Δ + γ·(1 + 0.34·u₂)]

其中 Rn 为净辐射（MJ/m²/d），G 为土壤热通量（日尺度取 0），T 为日均气温（℃），u₂ 为 2 m 高处风速（m/s），es 为饱和水汽压（kPa），ea 为实际水汽压（kPa），Δ 为饱和水汽压曲线斜率（kPa/℃），γ 为湿度计常数（kPa/℃）。

计算链与数据流

逐日计算实现

下面给出一个可直接运行的 Python 实现，输入为单日气象记录与作物系数，输出 ET₀、ETc 与 IRn。所有中间量按 FAO-56 定义计算，避免用近似常数替代。

import math

def saturation_vapor_pressure(t):
    # Tetens 公式，t 单位 ℃，返回 kPa
    return 0.6108 * math.exp(17.27 * t / (t + 237.3))

def et0_daily(tmax, tmin, rhmax, rhmin, u2, rs, lat_deg, doy):
    tmean = (tmax + tmin) / 2.0
    # 饱和水汽压曲线斜率 Δ
    es_t = saturation_vapor_pressure(tmean)
    delta = 4098 * es_t / (tmean + 237.3) ** 2
    # 饱和水汽压 es 与实际水汽压 ea
    es = (saturation_vapor_pressure(tmax) + saturation_vapor_pressure(tmin)) / 2.0
    ea = (saturation_vapor_pressure(tmin) * rhmax / 100.0
          + saturation_vapor_pressure(tmax) * rhmin / 100.0) / 2.0
    # 湿度计常数 γ（标准大气压 101.3 kPa）
    gamma = 0.000665 * 101.3
    # 净短波辐射 Rns，草地参考作物反照率 α = 0.23
    rns = (1 - 0.23) * rs
    # 净长波辐射 Rnl
    sigma = 4.903e-9  # MJ/(K^4·m^2·d)
    tk_max, tk_min = tmax + 273.16, tmin + 273.16
    rnl = (sigma * ((tk_max ** 4 + tk_min ** 4) / 2.0)
           * (0.34 - 0.14 * math.sqrt(ea))
           * (1.35 * rs / (0.75 * extraterrestrial_radiation(lat_deg, doy)) - 0.35))
    rn = rns - rnl
    # Penman-Monteith，日尺度 G = 0
    num = 0.408 * delta * rn + gamma * (900 / (tmean + 273)) * u2 * (es - ea)
    den = delta + gamma * (1 + 0.34 * u2)
    return num / den

def extraterrestrial_radiation(lat_deg, doy):
    # 简化 Ra 计算，用于 Rnl 中的 Rs/Rso 比值
    phi = math.radians(lat_deg)
    dr = 1 + 0.033 * math.cos(2 * math.pi * doy / 365)
    dec = 0.409 * math.sin(2 * math.pi * doy / 365 - 1.39)
    ws = math.acos(-math.tan(phi) * math.tan(dec))
    ra = (24 * 60 / math.pi) * 0.0820 * dr * (
        ws * math.sin(phi) * math.sin(dec)
        + math.cos(phi) * math.cos(dec) * math.sin(ws))
    return ra

def effective_rain(p):
    if p <= 25:
        return 0.8 * p
    return 20 + 0.7 * (p - 25)

def daily_irrigation_need(met, kc):
    et0 = et0_daily(met["tmax"], met["tmin"], met["rhmax"], met["rhmin"],
                    met["u2"], met["rs"], met["lat"], met["doy"])
    etc = kc * et0
    pe = effective_rain(met["p"])
    irn = max(0.0, etc - pe)
    return {"ET0": round(et0, 2), "ETc": round(etc, 2),
            "Pe": round(pe, 2), "IRn": round(irn, 2)}

示例数据与结果判读

取华北平原某地块 6 月 15 日（DOY=166，纬度 37°N）的示例气象记录，玉米处于拔节期，Kc 取 0.95。以下数值为按上述公式直接计算得到的示例结果，用于说明判读方法，不代表实测。

判读时注意三点：一是 ET₀ 与 ETc 的差异完全由 Kc 决定，若 Kc 取错生育阶段，需水量会系统性偏移；二是降雨日必须走有效降雨折算，不能直接用 P 扣减；三是 IRn 为 0 表示当日无需灌溉，而不是“需水量为负”。

逐日批量计算与输出

实际使用中通常一次处理一段日期。将每日气象记录组织为字典列表，按日期顺序调用 daily_irrigation_need，即可得到逐日需水深度序列，供后续灌溉排程使用。

records = [
    {"date": "2024-06-15", "tmax": 32.0, "tmin": 20.0, "rhmax": 85, "rhmin": 40,
     "u2": 2.2, "rs": 22.0, "p": 0.0, "lat": 37.0, "doy": 166, "kc": 0.95},
    {"date": "2024-06-16", "tmax": 30.0, "tmin": 19.0, "rhmax": 90, "rhmin": 45,
     "u2": 1.8, "rs": 18.0, "p": 8.0, "lat": 37.0, "doy": 167, "kc": 0.97},
]

for r in records:
    met = {k: r[k] for k in ("tmax", "tmin", "rhmax", "rhmin", "u2", "rs", "p", "lat", "doy")}
    out = daily_irrigation_need(met, r["kc"])
    print(r["date"], out)

验证方法：对同一记录手工核算 Δ、es、ea 与 Rn 四个中间量，若与代码输出一致，则 ET₀ 结果可信；再检查 IRn 是否满足 IRn = max(0, ETc − Pe) 的截断规则。连续多日结果中，若出现 ET₀ 明显高于当地气候上限（如超过 9 mm/d），应优先排查 Rs 单位是否误用为 W/m² 而非 MJ/m²/d。