Python中NDVI时间序列的LWR实现代码及适配方案咨询

一、现有LWR脚本的复用性判断 #

大部分通用LWR（局部加权回归）脚本的核心是加权拟合局部样本，这个逻辑完全适配时间序列场景——只要将时间维度作为自变量、NDVI作为因变量即可。无需完全推翻现有实现，重点是做时间序列特有的针对性修改。

二、核心修改要点 #

1. 时间轴的数值化转换 #

• 将日期转换为连续数值（例如从序列起始日开始的天数：days_since_start = (date - start_date).dt.days），因为LWR要求输入连续型自变量，无法直接处理日期字符串或对象。
• 确保局部窗口内包含足够样本，适配Sentinel（约5天间隔）与Landsat（16天间隔）融合后的不规则采样特性。

2. 局部窗口的时间自适应设置 #

• 把通用LWR的“固定样本数窗口”改成固定时间跨度窗口（例如滑动覆盖30天），避免因采样间隔不均导致窗口内样本数量波动过大。
• 示例调整：将窗口参数从n_neighbors=10替换为time_window=pd.Timedelta(days=30)，通过时间差筛选窗口内的有效样本。

3. 权重函数的时间相关性优化 #

• 将通用LWR的空间距离权重替换为时间差权重：离目标时间越近的样本权重越高，比如采用高斯权重公式：weights = np.exp(-(time_diff ** 2) / (2 * sigma ** 2))，其中sigma为时间带宽参数，可根据数据采样频率调整（如以Sentinel为主时设为10天）。

4. 平滑与日插值的衔接处理 #

• 遵循先平滑、后插值的顺序：对每个目标日期运行LWR拟合得到平滑NDVI值，再补全所有日度时间节点的结果，避免直接对原始噪声数据插值导致误差传递。

三、代码修改示例（Python） #

假设现有通用LWR脚本用于拟合x-y数据，修改为时间序列适配版本：

import pandas as pd
import numpy as np

def lwr_time_series(ts_data, target_date, time_window=30, sigma=10):
    # ts_data: 包含'date'（日期列）和'ndvi'（NDVI值列）的DataFrame
    start_date = ts_data['date'].min()
    # 转换日期为连续天数
    ts_data['days'] = (ts_data['date'] - start_date).dt.days
    target_days = (target_date - start_date).days
    
    # 筛选时间窗口内的样本
    window_mask = np.abs(ts_data['days'] - target_days) <= time_window
    window_data = ts_data[window_mask]
    
    # 样本数不足时返回最近原始值
    if len(window_data) < 3:
        return ts_data.iloc[(ts_data['days'] - target_days).abs().idxmin()]['ndvi']
    
    # 计算时间差权重
    time_diff = window_data['days'] - target_days
    weights = np.exp(-(time_diff ** 2) / (2 * sigma ** 2))
    
    # 加权线性回归（LWR核心逻辑）
    X = np.column_stack((np.ones(len(window_data)), window_data['days']))
    y = window_data['ndvi']
    W = np.diag(weights)
    beta = np.linalg.inv(X.T @ W @ X) @ X.T @ W @ y
    
    # 预测目标日期的平滑NDVI值
    return beta[0] + beta[1] * target_days

# 生成日度平滑序列
ndvi_fused = pd.read_csv('fused_ndvi.csv', parse_dates=['date'])
date_range = pd.date_range(start=ndvi_fused['date'].min(), end=ndvi_fused['date'].max(), freq='D')
smoothed_daily = []
for date in date_range:
    smoothed_val = lwr_time_series(ndvi_fused, date)
    smoothed_daily.append({'date': date, 'ndvi_smoothed': smoothed_val})

smoothed_daily_df = pd.DataFrame(smoothed_daily)

四、无需完全替换的场景 #

如果现有脚本已支持自定义权重和自变量输入，只需完成以下3步即可：

• 将时间维度作为唯一自变量传入
• 替换权重计算逻辑为时间相关函数
• 循环遍历每个目标日运行LWR拟合

内容的提问来源于stack exchange，提问作者Estefanía Pizarro Arias