遥感时序分析

引言

遥感时间序列分析把同一地点在不同时刻的观测连成一条曲线，再从这条曲线上读出地表的变化规律。它回答的问题和单期影像完全不同：不是「这块地现在是什么」，而是「这块地的植被什么时候开始返青、什么时候枯黄、过去二十年是变绿还是退化、哪一年被砍伐或火烧」。这类问题服务于物候监测、作物估产、生态退化评估、碳循环建模与灾害恢复追踪。

难点集中在数据侧而非算法侧。光学影像被云污染，一个像元一年可能只有十几次有效观测，且分布不均；不同传感器的波段响应与分辨率不同，直接拼接会引入系统偏差；时序里的突变（火灾、砍伐、洪水）与渐变（气候趋势、缓慢退化）混在一起，需要区分对待。算法本身（谐波拟合、Mann-Kendall、BFAST）都是成熟方法，真正的工程量在于把脏的、稀疏的、多源的观测整理成可建模的干净序列。

本文聚焦「曲线怎么建、怎么建模、怎么读」，与变化检测的边界要讲清楚：变化检测篇关注两期或多期之间的变化图斑与变化类型，回答「哪里变了、变成什么」；本篇关注连续时序的物候、趋势与断点，回答「曲线形态如何、何时突变、长期趋势是升是降」。两者共享时序数据，但目标与输出不同，前者的产物是变化图，后者的产物是物候参数、趋势斜率与断点年份。

阅读前建议先理解光谱指数的构造与异常值处理，见 光谱指数计算 ；时序数据的可用性取决于云掩膜质量，见 云检测与云掩膜工程实践 。

目录

1. 时序分析的问题类型与边界
2. 时序数据的构建与合成
3. 云填补与缺口重建
4. 时序平滑与谐波拟合
5. 物候参数提取
6. 趋势检测：Mann-Kendall 与 Theil-Sen
7. 断点检测：BFAST 与 CCDC
8. 多传感器时序融合与协调
9. 工程实现与质量评估

1. 时序分析的问题类型与边界

先分清几类时序问题的输出差异，它们的建模方式与验证方法都不同。

物候与趋势是连续变化，断点是突变事件。这个区分决定了预处理策略：做趋势分析时要用平滑抑制噪声但保留长期信号；做断点检测时必须保留突变，绝不能用填补数据，否则断点会被抹平。

与变化检测的边界再强调一次。两期影像相减得到变化图，这是变化检测；把十年观测拟合成曲线求斜率，这是趋势分析。前者对配准与辐射一致性极敏感，后者对时序长度与观测密度敏感。工程上常同时使用：先用变化检测找出可疑图斑，再用时序分析确认变化发生的时间与类型。

2. 时序数据的构建与合成

时序数据的第一步是把零散的观测整理成规则时间网格上的序列。以 Sentinel-2 为例，双星重访 5 天，但受云影响，实际有效观测间隔可能长达数周。构建方式是先做云掩膜，再把有效观测按时间排序，形成一个不规则的时间戳数组与对应的指数值数组。

import numpy as np
import pandas as pd

def build_series(obs):
    # obs: [{"date": "2020-03-05", "ndvi": 0.31, "quality": 0.9}, ...]
    df = pd.DataFrame(obs)
    df["date"] = pd.to_datetime(df["date"])
    df = df.sort_values("date").drop_duplicates("date")
    df["ndvi"] = df["ndvi"].clip(-1.0, 1.0)              # 物理范围约束
    df["doy"] = df["date"].dt.dayofyear                    # 年积日
    df["year"] = df["date"].dt.year
    return df.reset_index(drop=True)

两个纪律。第一，指数值要裁剪到物理范围（NDVI 在 -1 到 1），越界值几乎都是云残留或异常像元，不裁会污染后续拟合。第二，要保留质量标记与观测数量，没有这两个信息，下游无法判断某个像元的值可信度。

时序合成是另一条路线。当观测太稀时，先按周期合成再分析：10 天、16 天或月合成，用最大值合成或质量优先合成。最大值合成对植被最常用，因为它抑制云、偏向晴空观测，代价是对负向变化（如落叶、灾害）不敏感。质量优先合成按质量波段排序取最优观测，精度更高但依赖质量波段的可靠性。

合成周期与有效观测数（Sentinel-2 双星，5 天重访）
区域        10 天窗口    16 天窗口    30 天窗口
赤道雨林    1~2          2~3          4~7
温带大陆    4~6          6~9          12~18
干旱沙漠    8~10         13~16        25~30

经验规则是每个合成窗口至少要有 3 个有效观测，中值或最大值合成才有统计意义。低于 3 个时应拉长周期，或退回「用邻近时相填补」。

3. 云填补与缺口重建

缺口填补是时序分析绕不开的一步，但也是最容易引入虚假信号的一步。方法分三类。

时间维插值：线性、样条或三次插值，用前后有效观测填补中间缺口。简单但对长缺口不可靠，连续 30 天以上的缺口插出来的直线与真实物候不符。

时序拟合填补：用谐波模型拟合整年曲线，用拟合值填补缺口。平滑、抗噪，但会抹平突变，且假设地物变化是季节性的。

时空邻域填补：结合空间上相邻像元的信息，用 NSPI、GNSPI 这类时空加权方法。能保留空间细节，计算量大。

import numpy as np

def harmonic_fill(doy, values, max_gap_days=45, n_harm=2):
    # 用谐波拟合填补缺口，缺口过长的位置保留 NaN
    w = 2 * np.pi / 365.25
    cols = [np.ones_like(doy, dtype="float64")]
    for k in range(1, n_harm + 1):
        cols += [np.sin(k * w * doy), np.cos(k * w * doy)]
    A = np.stack(cols, axis=1)
    coef, *_ = np.linalg.lstsq(A, values, rcond=None)
    filled = A @ coef
    # 缺口超过阈值的位置不填，避免虚构长段信号
    gap = np.diff(doy, prepend=doy[0])
    long_gap = np.convolve(gap, np.ones(max_gap_days), "same") > max_gap_days
    return np.where(long_gap, np.nan, filled)

工程纪律有三条。第一，限制插值跨度，超过阈值的缺口标记为缺失而不是强行填补。第二，填补结果必须与原观测区分标记，元数据里记录每个像元哪些时刻是填补的。第三，做断点或异常检测时只用原始观测，绝不用填补值，否则火灾、砍伐这类突变会被抹平成渐变，断点检测彻底失效。

4. 时序平滑与谐波拟合

谐波模型是植被时序分析的核心工具。它假设植被指数随时间呈周期性变化，用有限阶的正弦余弦叠加拟合季节曲线。一二阶谐波通常就能捕捉双季作物的双峰，三四阶用于更复杂的物候。

import numpy as np

def harmonic_features(doy, values, n_harm=2):
    w = 2 * np.pi / 365.25
    cols = [np.ones_like(doy, dtype="float64")]
    for k in range(1, n_harm + 1):
        cols += [np.sin(k * w * doy), np.cos(k * w * doy)]
    A = np.stack(cols, axis=1)
    coef, *_ = np.linalg.lstsq(A, values, rcond=None)
    fitted = A @ coef
    resid = values - fitted
    amp = [np.hypot(coef[1 + 2 * (k - 1)], coef[2 + 2 * (k - 1)]) for k in range(1, n_harm + 1)]
    return {"fitted": fitted, "resid": resid, "amplitude": amp, "coef": coef}

拟合得到的谐波系数本身就是特征。一阶振幅反映季节性强度，相位反映峰值出现的时间（即物候早晚），常数项是年均值。这些特征可以送入随机森林做作物分类，比原始时序更抗噪、维度更低。

平滑的另一条路线是滤波。Savitzky-Golay 滤波在保留峰值形状的同时去噪，是 MODIS 时序处理的经典方法；双逻辑斯蒂拟合用两个 S 形函数拼出生长与衰老段，参数直接对应物候指标，在作物时序里常用。

from scipy.signal import savgol_filter

def sg_smooth(values, window=7, order=2):
    # 窗口须为奇数，order 小于 window
    return savgol_filter(values, window_length=window, polyorder=order)

平滑强度要按观测密度调。观测密时可以用小窗口保留细节，观测稀时大窗口更稳但会压低峰值。平滑前一定要先按时间排序并处理重复时刻，否则滤波会把时间上不连续的点当作相邻点处理。

5. 物候参数提取

物候参数是时序分析最有业务价值的输出，描述植被生长季的时间特征。

常用物候指标
SOS  生长季开始（Start of Season）
EOS  生长季结束（End of Season）
LOS  生长季长度 = EOS - SOS
POS  峰值出现时间（Peak of Season）
POP  峰值大小（Peak value）
AUC  曲线下面积，近似累积生产力

提取方法有三类。

阈值法：设定一个基准值（如振幅的 20% 或年最小值的某个比例），曲线第一次跨过该阈值的时间即 SOS，最后一次跨过即 EOS。简单、可复现，但阈值选择影响结果，且对噪声敏感。

导数法：找曲线的最大上升速率点作为 SOS，最大下降速率点作为 EOS。对平滑后的曲线效果稳定，是较主流的方法。

拟合参数法：用双逻辑斯蒂或谐波拟合，从拟合参数直接解析出物候时间。抗噪最好，但假设曲线形状符合模型，对不规则曲线会失真。

import numpy as np

def phenology_threshold(doy, values, base_frac=0.2):
    base = np.nanmin(values)
    amp = np.nanmax(values) - base
    thr = base + base_frac * amp
    above = values > thr
    if not above.any():
        return {"SOS": np.nan, "EOS": np.nan, "POS": np.nan, "LOS": np.nan}
    idx = np.where(above)[0]
    sos, eos = doy[idx[0]], doy[idx[-1]]
    pos = doy[np.nanargmax(values)]
    return {"SOS": sos, "EOS": eos, "POS": pos, "LOS": (eos - sos) % 365}

物候提取的前提是曲线完整。观测太稀或缺口太多时，SOS 与 EOS 的估计误差可达数周。因此物候产品通常附带一个质量标记，标注该像元的观测数与拟合残差，下游按质量过滤。

物候的验证依赖地面观测网络，如物候相机与人工物候记录。遥感 SOS 与地面观测之间存在系统性偏差（遥感看的是冠层绿度，地面看的是展叶或开花），跨区域比较时要确认物候定义一致。

6. 趋势检测：Mann-Kendall 与 Theil-Sen

趋势分析的目的是判断某个像元的植被指数在多年尺度上是升是降，以及这个趋势是否显著。

最小二乘线性回归直接给出斜率，但它对离群值敏感，一年火灾就能把斜率拉偏。更稳健的组合是非参数方法：Theil-Sen 斜率估计配合 Mann-Kendall 趋势检验。

Theil-Sen 斜率取所有点对斜率的中位数，对离群值鲁棒。Mann-Kendall 检验是非参数检验，不假设数据分布，检验序列是否存在单调趋势。

import numpy as np
from scipy.stats import kendalltau

def theil_sen_slope(years, values):
    n = len(values)
    slopes = []
    for i in range(n):
        for j in range(i + 1, n):
            if years[j] != years[i]:
                slopes.append((values[j] - values[i]) / (years[j] - years[i]))
    return float(np.median(slopes)) if slopes else np.nan

def mann_kendall(values):
    # 返回 tau 与双侧 p 值，p 小表示趋势显著
    tau, p = kendalltau(np.arange(len(values)), values)
    return float(tau), float(p)

两者配合使用：Theil-Sen 给趋势大小，Mann-Kendall 给显著性。业务上通常只保留「显著且斜率绝对值超过阈值」的像元，避免把噪声波动当趋势。

趋势分析有几个陷阱。第一是时序长度，短于 10 年的序列趋势估计不稳，且容易把周期性波动误判为趋势。第二是观测数不均，早期观测少、近期观测多的序列，趋势会被观测密度差异放大。第三是突变污染，一次火灾会在灾后留下数年恢复期，把整个序列的斜率拉负，这种像元应先做断点检测再分段算趋势。

季节 Mann-Kendall 是对标准 MK 的改进，先扣除季节分量再检验趋势，适合植被这类强季节性序列，能显著降低假阳性。

7. 断点检测：BFAST 与 CCDC

断点检测要回答「这条曲线在第几年发生了结构性突变」。突变分两类：渐变型（气候驱动的缓慢变化）与突变型（火灾、砍伐、洪水、城市化）。BFAST 与 CCDC 是两个主流框架。

BFAST（Breaks For Additive Season and Trend）把序列分解为趋势、季节与残差三部分，在每个分量上检测断点，再把断点合并成最终的突变时刻。它对季节性序列效果好，能区分季节内的异常与结构性的突变。

CCDC（Continuous Change Detection and Classification）用滑动窗口做逐时相预测，当连续多次观测偏离预测模型的置信区间时判定为变化，并重新拟合模型。它适合近实时的变化监测，能给出变化发生的时间与前后类别。

import numpy as np

def ccdc_lite(doy, values, n_harm=3, consec=3, k=3.0):
    # 简化版 CCDC：滑动拟合谐波模型，连续多次超限即报断点
    w = 2 * np.pi / 365.25
    cols = [np.ones_like(doy, dtype="float64")]
    for h in range(1, n_harm + 1):
        cols += [np.sin(h * w * doy), np.cos(h * w * doy)]
    A = np.stack(cols, axis=1)
    coef, *_ = np.linalg.lstsq(A, values, rcond=None)
    pred = A @ coef
    resid = values - pred
    sigma = np.median(np.abs(resid - np.median(resid))) * 1.4826 + 1e-6
    out = np.abs(resid) > k * sigma
    # 连续超限才确认，避免单点噪声触发
    breaks = []
    run = 0
    for i, flag in enumerate(out):
        run = run + 1 if flag else 0
        if run == consec:
            breaks.append(doy[i - consec + 1])
    return breaks, sigma

用中位数绝对偏差（MAD）而非标准差估计噪声水平，是因为突变本身是离群值，标准差会被它拉大从而降低灵敏度，1.4826 是使 MAD 与标准差在正态分布下可比的一致性常数。

断点检测的验证最难。地面事件记录（采伐许可、火灾报告）往往不完整，且时间粒度粗。实践中常用交叉验证：用一部分观测拟合模型，预测剩余观测，看断点前后的预测误差是否显著变化。

8. 多传感器时序融合与协调

单一传感器的时序常不够长或不够密。Landsat 从 1984 年至今提供 30 米观测但重访 16 天，Sentinel-2 从 2015 年起提供 10 米观测重访 5 天。把两者融合能同时得到长历史与高频率，但必须先解决协调问题。

主要差异与处理：

协调的核心是交叉定标：选取两传感器同期过境的影像，统计同一地物在两传感器上的反射率关系，拟合线性或二次变换，把一个传感器归一化到另一个的基准。不做这步直接拼接，序列会在传感器切换处出现台阶，趋势分析会把这个台阶误判为突变。

import numpy as np

def cross_calibrate(ref, tgt):
    # 用同期观测拟合 tgt -> ref 的线性关系
    mask = np.isfinite(ref) & np.isfinite(tgt)
    a, b = np.polyfit(tgt[mask], ref[mask], 1)
    return a, b        # ref_hat = a * tgt + b

融合后的序列还要做一致性检查：统计传感器切换点前后一定窗口内的均值差异，差异超过阈值的像元标记为可疑。这一步能抓出定标异常或配准错位。

SAR 是光学时序的重要补充。Sentinel-1 的 C 波段后向散射不受云影响，在常年多云区域能维持观测连续性。但后向散射与植被指数的物理含义不同，不能直接拼进同一条 NDVI 序列，只能作为独立通道或用于填补缺口时的辅助信息，SAR 的处理见 SAR 与 InSAR 处理 。

9. 工程实现与质量评估

时序分析的工程实现要考虑数据量与计算效率。一个中等区域（如一个省）十年 Sentinel-2 时序，按 10 米分辨率可能有数十亿像元乘上百个时刻，无法一次读入。做法是按块处理：把区域切成瓦片，逐瓦片读入时序立方体，逐像元做拟合与检测，输出物候、趋势与断点的栅格。

import numpy as np
import xarray as xr

def analyze_chunk(cube, doy, n_harm=2):
    # cube: (t, y, x)，返回物候、趋势、断点栅格
    t, h, w = cube.shape
    sos = np.full((h, w), np.nan, dtype="float32")
    slope = np.full((h, w), np.nan, dtype="float32")
    for y in range(h):
        for x in range(w):
            v = cube[:, y, x]
            if np.isfinite(v).sum() < 8:            # 观测太少，跳过
                continue
            valid = np.isfinite(v)
            ph = phenology_threshold(doy[valid], v[valid])
            sos[y, x] = ph["SOS"]
            slope[y, x] = theil_sen_slope(doy[valid] / 365.25, v[valid])
    return sos, slope

计算量大时用 Dask 或 xarray 的并行能力，把瓦片任务分发到多核或多机。拟合与检测都是逐像元独立操作，天然适合并行，几乎线性加速。

质量评估要报三类指标。第一是覆盖度：有多少像元能提取出物候，多少因观测不足被标记缺失。第二是拟合质量：残差的标准差与 R 方，反映模型是否合适。第三是验证精度：物候估计与地面观测的偏差、趋势检测与已知事件的吻合率。

输出产品必须附带元数据：传感器与产品级别、时间范围、合成或填补方法、模型参数、缺失与填补的标记规则。没有这些信息，下游无法复现，也无法判断某像元值的可信度。元数据里还应记录处理链的版本号，方便在同一区域重跑时对齐口径。

权衡取舍

• 原始观测 vs 合成序列：原始观测保留突变、适合断点检测；合成序列更连续、适合物候与趋势，按分析目标选择。
• 插值填补 vs 标记缺失：插值让曲线连续但会抹平突变，突变敏感的分析必须保留缺失。
• 阈值法 vs 拟合参数法提取物候：阈值法简单可复现但敏感于阈值，拟合参数法抗噪但假设曲线形状。
• 线性回归 vs Theil-Sen：线性回归直观但受离群值影响，Theil-Sen 稳健但计算量大，长时间序列推荐后者。
• 标准 MK vs 季节 MK：强季节性序列用标准 MK 假阳性高，季节 MK 更准但实现更复杂。
• BFAST vs CCDC：BFAST 适合离线历史分析，CCDC 适合近实时监测与变化分类。
• 多传感器融合 vs 单传感器：融合得到更长更密的序列，但引入定标与分辨率协调的复杂度。
• 逐像元 vs 逐对象：逐像元简单通用，逐对象能利用空间上下文但依赖分割质量。

常见坑清单

• 用填补数据做断点检测：现象是火灾与砍伐的突变被抹平，原因是插值假设变化平滑，规避方法是断点检测只用原始观测。
• 指数不裁剪物理范围：现象是拟合出现异常峰值，原因是云残留像元产生越界值，规避方法是把 NDVI 裁剪到 -1 到 1。
• 长缺口强行插值：现象是物候起始日算错数周，原因是长缺口插值偏离真实曲线，规避方法是限制插值跨度并标记缺失。
• 不做交叉定标直接拼接：现象是传感器切换处出现台阶被误判为趋势，原因是波段响应差异，规避方法是同期观测做交叉定标。
• 短序列做趋势：现象是趋势不显著或方向反复，原因是序列短于 10 年，规避方法是延长序列或降低置信要求。
• 忽略观测数不均：现象是趋势被观测密度差异放大，原因是早期观测稀疏，规避方法是按观测数分层或加权。
• 物候阈值选得太高或太低：现象是 SOS 系统性偏早或偏晚，原因是阈值未按振幅自适应，规避方法是用振幅比例而非固定值。
• 平滑过度压平峰值：现象是 POP 被低估、POS 偏移，原因是滤波窗口过大，规避方法是按观测密度调窗口。
• 用标准差估计时序噪声：现象是突变检测灵敏度低，原因是突变拉大了标准差，规避方法是改用 MAD。
• 输出缺元数据：现象是结果无法复现、可信度无法判断，原因是未记录方法与填补标记，规避方法是把处理链版本与标记规则写入元数据。

小结

遥感时序分析的价值在于把「某一刻的状态」升级为「随时间演化的规律」。技术上没有特别新的算法，难的是把被云污染、观测稀疏、多源异构的数据整理成可建模的干净序列，并在物候、趋势、断点三类问题之间选对建模与预处理策略。填补与平滑方便了趋势与物候分析，却会破坏断点检测，这条边界必须在流水线设计时就划清。

落地建议从一个小区域、一个明确问题开始：先建好云掩膜与时序立方体，跑通物候提取并用地面对照验证；再加趋势检测，用季节 MK 配合 Theil-Sen；最后上断点检测，用已知的火灾或采伐记录验证。每一步都要输出质量标记，让下游知道哪些像元可信、哪些是填补的。

下一步可以对照 变化检测工程实践 理解变化图与时序曲线两种视角的分工，也可以结合 光谱指数计算 补充指数构造与异常值处理的细节。