Python实现NDVI双逻辑斯蒂曲线拟合（对标R包greenbrown）

我需要在Python中实现Beck等人2006年提出的NDVI双逻辑斯蒂曲线拟合方法，该方法已在R语言的greenbrown包中通过FitDoubleLogBeck.R实现，但我当前的Python拟合效果未达预期，希望完全复刻该R脚本的拟合逻辑。每个样地每年有9个NDVI观测值，以下是我的尝试代码及拟合结果：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def double_logistic_function(t, wNDVI, mNDVI, S, A, mS, mA):
    return wNDVI + (mNDVI / (1 + np.exp(-mS * (t - S)))) * (1 - 1 / (1 + np.exp(-mA * (t - A))))

def weight_function(t, S, A, r):
    tr = 100 * (t - S) / (A - S)
    tr = np.clip(tr, 0, 100)
    return np.exp(-np.abs(r) / (1 + tr / 10))

def fit_curve(t, ndvi_observed):
    initial_guess = [np.min(ndvi_observed), np.max(ndvi_observed), np.mean(t), np.mean(t), 1, 1]
    params, _ = curve_fit(double_logistic_function, t, ndvi_observed, p0=initial_guess)
    residuals = ndvi_observed - double_logistic_function(t, *params)
    weights = weight_function(t, params[2], params[3], residuals)
    params, _ = curve_fit(double_logistic_function, t, ndvi_observed, p0=initial_guess, sigma=weights)
    return params

##Example usage with 9 observations over the year
t_9_obs = np.array([30, 60, 90, 120, 150, 180, 210, 240, 270])  

##ndvi_9_observed = np.array([0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2])
ndvi_9_observed = np.array([0.58, 0.583, 0.713, 0.807, 0.832, 0.878, 0.886, 0.863, 0.717])


##Fit the curve
parameters = fit_curve(t_9_obs, ndvi_9_observed)

##Plot the observed NDVI values
plt.scatter(t_9_obs, ndvi_9_observed, label='Observed NDVI')

##Generate points for the fitted curve
t_fit = np.linspace(min(t_9_obs), max(t_9_obs), 1000)
ndvi_fit = double_logistic_function(t_fit, *parameters)

##Plot the fitted curve
plt.plot(t_fit, ndvi_fit, label='Fitted Curve', color='red')

plt.xlabel('Day of the Year')
plt.ylabel('NDVI')
plt.legend()
plt.title('Double Logistic Curve Fitting for 9 NDVI Observations')
plt.show()

拟合结果对比 #

• 我的代码运行拟合结果：
  

• Beck等人2006年提出的双逻辑斯蒂曲线模型示意图：
  

内容的提问来源于stack exchange，提问作者Damiaan van Harteveld