Python 中的 NDVI 双 Logistic 曲线拟合

Question

我正在尝试用 Python 进行 NDVI 双逻辑曲线拟合。这种双逻辑曲线拟合由 Beck 等人发表。 2006 年。有一个 R 包 greenbrown 已经完成了一年多的工作，到目前为止我还没有得到很好的效果或者没有曲线拟合。我想在 R 中进行与此代码相同的曲线拟合：https://github.com/r-forge/greenbrown/blob/master/pkg/greenbrown/R/FitDoubleLogBeck.R

贝克等人。 2006

我每年每个地块有 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()

此代码的结果：

Answer 1

我认为这里发生了两个问题：

1. 你的数学模型和你的代码不匹配。
2. 优化器在优化此函数时采取了大得离谱的步骤。

你的数学模型和你的代码不匹配，因为模型是将两条 sigmoid 曲线相加，但你却将它们相乘。另外，模型使用 mNDVI - wNDVI，但代码使用 mNDVI。

我使用数学模型重新实现了你的曲线，并将其分成几行，以便更容易理解。

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

您应该始终测试曲线函数，以确保您正在优化您认为正在优化的函数。

我发现的另一个问题是优化器对 S 和 A 参数采取了大得离谱的步骤。

这些是优化器尝试的前 8 个 S、A 值。

S, A=(150.0, 150.0)
S, A=(150.0, 150.0)
S, A=(150.0, 150.0)
S, A=(150.00000223517418, 150.0)
S, A=(150.0, 150.00000223517418)
S, A=(150.0, 150.0)
S, A=(150.0, 150.0)
S, A=(15721.812769000722, 15722.840546778361)

发生这种情况时，我发现确保 X 和 Y 具有相似的大小通常很有帮助。为此，我重新参数化了您的 X 变量，使其从 [0, 1] 而不是 [0, 365] 变化。

t_9_obs = np.array([30, 60, 90, 120, 150, 180, 210, 240, 270]) / 365

这样做之后，我得到了更合理的配合。

完整代码如下：

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):
    sigmoid1 = 1 / (1 + np.exp(-mS * (t - S)))
    sigmoid2 = 1 / (1 + np.exp( mA * (t - A)))
    seasonal_term = sigmoid1 + sigmoid2 - 1
    return wNDVI + (mNDVI - wNDVI) * seasonal_term

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]) / 365

##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 * 365, 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 * 365, 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()

警告：如果您正在解释此回归的系数，请务必小心。该曲线有多个等效拟合。例如，如果您有系数 [wNDVI, mNDVI, S, A, mS, mA]，则 [wNDVI, mNDVI, A, S, -mA, -mS] 是给出相同曲线的等效解。