土壤湿度反演IEM模型代码

以下是IEM模型的Python代码：

import numpy as np
from scipy.optimize import minimize

def iem_model(theta, t, P, R):
    """
    IEM模型函数
    参数：
        theta: 模型参数，[a, b, c, d]
        t: 时间序列
        P: 降雨量序列
        R: 壤中含水量序列
    返回值：
        S: 土壤湿度序列
    """
    a, b, c, d = theta
    S = np.zeros_like(R)
    for i in range(1, len(R)):
        dR = R[i] - R[i-1]
        if dR >= 0:
            S[i] = S[i-1] + P[i] - dR/1000.0
        else:
            S[i] = S[i-1] + P[i] - dR*(1.0-a*S[i-1]**b-c*P[i]**d)/1000.0
        if S[i] < 0:
            S[i] = 0
    return S

def cost(theta, t, P, R, S_obs):
    """
    代价函数，使用平方和误差
    """
    S = iem_model(theta, t, P, R)
    return np.sum((S - S_obs)**2)

def optimize_theta(theta0, t, P, R, S_obs):
    """
    使用scipy.optimize.minimize求解最优参数
    """
    res = minimize(cost, theta0, args=(t, P, R, S_obs), method='Nelder-Mead', options={'maxiter':1000})
    return res.x

该代码实现了IEM模型及其参数优化。其中，iem_model函数为IEM模型的主体部分，接受时间序列t、降雨量序列P和壤中含水量序列R，并返回土壤湿度序列S。cost函数为代价函数，使用平方和误差。optimize_theta函数使用scipy.optimize.minimize求解最优参数。