变化检测与影像时间序列分析

引言

变化检测要回答的问题很朴素：同一地点在两个时间点或一段时间内，地表发生了什么改变。它的应用极广，从城市扩张、森林砍伐、耕地撂荒，到灾后评估、冰川退缩、水体消长。相比单期分类，变化检测的难度不在识别地物，而在排除「假变化」。

假变化的来源几乎都是干扰而非真实地变：两期影像的配准误差、太阳高度角与观测几何差异、大气状况不同、物候阶段不同、云与阴影残留、传感器更换带来的波段差异。任何一项处理不当，都会让变化图斑里混入大量噪声，让结果无法使用。

本文按「经典方法到时间序列」的顺序展开。先厘清变化检测的任务类型，再讲代数差分与阈值、变化向量分析、IR-MAD 这些双期方法，然后进入深度学习变化检测，最后落到影像时间序列的断点检测、谐波回归与 BFAST、CCDC 这类连续监测方法，并给出后处理与验证的工程要点。

指数差分的计算见 植被指数与光谱指数计算 ，分类后比较依赖单期分类精度，配准与投影一致性见 遥感坐标系统与投影 。

目录

1. 变化检测的问题定义与类型
2. 代数差分与阈值法
3. 变化向量分析 CVA
4. IR-MAD 与统计方法
5. 分类后比较法
6. 深度学习变化检测
7. 影像时间序列与断点检测
8. 谐波回归与 BFAST 与 CCDC
9. 变化图斑后处理与验证

1. 变化检测的问题定义与类型

变化检测不是单一任务，先分清类型才能选方法。

二值变化只需判断「变没变」，多类变化还要判断「变成什么」。后者难度高一个数量级，因为要同时保证前后两期的分类精度，误差会叠加。

变化检测还需明确三个前置条件：第一，两期影像必须经过一致的辐射与大气校正，最好都是地表反射率；第二，空间配准误差要控制在一个像元以内，最好半个像元；第三，要区分真实地变与物候、光照这类周期性或系统性变化。这三条不满足，任何算法都救不回来。

变化检测的评价单位也要事先确定：像元级、对象级还是图斑级。像元级评价最严格，但边界像元的类别天然模糊；对象级或图斑级评价更贴近应用，比如「这片林地被识别为变化」比「这 500 个像元里 380 个判对」更有意义。评价单位要与应用需求一致，否则指标好看但不可用。

2. 代数差分与阈值法

最直接的方法是影像差分：两期影像逐像元相减。

dImage = Image_t2 - Image_t1

单波段差分对噪声敏感，实践中更常用指数差分，因为它把多维信息压缩到一维，且对光照鲁棒。植被变化用 dNDVI，燃烧用 dNBR，水体用 dNDWI：

import numpy as np

def dndvi(ndvi_t1, ndvi_t2, valid):
    diff = ndvi_t2 - ndvi_t1
    return np.where(valid, diff, np.nan)             # 无效像元置 NaN

差分的核心难题是阈值选取。阈值太低会把噪声判为变化，太高会漏掉弱变化。常用方法：

• 经验阈值：如 dNDVI 绝对值大于 0.2 判为变化，简单但需本地标定。
• Otsu 大津法：在差分直方图上自动找类间方差最大的阈值，适合双峰分布。
• K-means：把差分值聚成两类或三类，取类边界作阈值。
• 统计阈值：均值加减 k 倍标准差，k 取 1.5 到 2.5。

from skimage.filters import threshold_otsu

valid_diff = diff[~np.isnan(diff)]
thr = threshold_otsu(valid_diff)                     # 自动阈值
change = np.abs(diff) > thr

除了差分，比值法也是常用的代数方法：用两期影像相除而非相减。比值法对乘性的辐射差异不敏感，因为比值会抵消整体的比例变化，适合未做严格校正的场景。

rImage = Image_t2 / Image_t1

差值与比值各有侧重：

阈值法的最大问题是变化像元占比很小时直方图不呈双峰，Otsu 会失效。这时应先做变化检测的粗筛，再在候选区做精细判定，或者改用能处理类别不平衡的统计方法。

差分法的另一局限是无法给出变化类型。它只能回答「变了」，不能回答「从什么变成什么」，多类变化需求要交给后续方法。

3. 变化向量分析 CVA

变化向量分析（CVA）把每个像元的多波段差分视为一个向量，用向量的模长表示变化强度，用方向表示变化类型。

Δv = [B2_t2 - B2_t1, B3_t2 - B3_t1, ..., Bn_t2 - Bn_t1]
magnitude = sqrt(sum(Δv_i^2))

模长大的像元变化剧烈，模长小的是未变化。方向角则编码了变化类型，比如植被到建筑与植被到水体的向量方向不同。

def cva(bands_t1, bands_t2):
    dv = bands_t2 - bands_t1                         # 形状 (n_bands, H, W)
    mag = np.sqrt((dv ** 2).sum(axis=0))             # 变化强度
    angle = np.arctan2(dv[1], dv[0])                 # 二维方向示例
    return mag, angle

CVA 的改进是加入角度约束，只有模长与方向都显著偏离未变化分布的像元才判为变化，能显著降低虚警。

CVA 的前提是波段数一致且已做一致校正。跨传感器做 CVA 要先做波段对齐与交叉标定，否则向量方向会被传感器差异污染。

阈值化 CVA 的判定逻辑是把模长与方向联合起来：

thr_mag = np.nanmean(mag) + 2.0 * np.nanstd(mag)     # 模长阈值
thr_ang = np.deg2rad(15.0)                           # 方向容差
change = (mag > thr_mag) & (np.abs(angle) > thr_ang) # 模长大且方向偏离

联合判定能显著降低虚警，因为单纯模长大可能来自光照差异，而方向偏离才更可能对应真实的地物类型改变。

4. IR-MAD 与统计方法

多元变化检测（MAD）基于典型相关分析，找出两期影像间相关性最弱的分量，这些分量承载了变化信息。

设两期影像 X 与 Y，MAD 求一组线性组合 a、b，使得 a’X 与 b’Y 的方差最大、相关最小。得到的 MAD 分量按变化信息量排序，前几个分量通常就包含了绝大部分变化。

IR-MAD 是 MAD 的迭代加权版本，它迭代地为每个像元分配权重，未变化像元权重高、变化像元权重低，从而在存在大量变化像元时仍能稳健估计未变化的背景分布。

def ir_mad_weights(X, Y, iters=10):
    n = X.shape[0]
    w = np.ones(n)                                   # 初始权重全为 1
    for _ in range(iters):
        Xw, Yw = X * w[:, None], Y * w[:, None]      # 在加权数据上做典型相关，得到 MAD 分量
        mad = compute_mad(Xw, Yw)
        stat = (mad ** 2).sum(axis=1)                # 卡方统计量
        w = 1 - chi2_cdf(stat, df=mad.shape[1])      # 更新权重，变化像元权重下降
    return w, stat

判定变化用卡方检验：统计量超过给定置信度对应的卡方分位数即判为变化。IR-MAD 的优势是无需人工设阈值、能自动处理辐射差异、对多波段联合利用充分，缺点是需要足够多的未变化像元来稳定估计背景，变化面积占比过高时会退化。

实践上，IR-MAD 常作为粗筛，输出变化概率图，再用阈值或聚类得到二值变化，最后叠加分类做变化类型判定。

MAD 分量的解释性也很有价值。第一个 MAD 分量往往对应最主要的变化方向，通过分析各波段在该分量上的载荷，可以推断变化的主导类型。例如载荷集中在近红外与短波红外，通常指向植被变化。

mad_vecs, mad_vars = compute_mad(X, Y)               # 特征向量与各分量方差
loadings = mad_vecs[:, 0]                            # 第一分量各波段权重
dominant = np.argmax(np.abs(loadings))               # 主导波段

5. 分类后比较法

分类后比较法先对两期影像各自分类，再逐像元比较类别标签，标签不同即判为变化。

def post_classification(c1, c2, ignore=0):
    valid = (c1 != ignore) & (c2 != ignore)
    change = (c1 != c2) & valid
    from_to = np.where(change, c1 * 100 + c2, 0)     # 编码变化类型
    return change, from_to

变化类型用类别编码相乘得到，便于统计各转移方向：

有了转移矩阵，就能得到各类之间的转换面积，这是土地利用变化分析的核心产出。

它的最大优点是直接给出变化类型，无需额外设计。最大缺点是误差叠加：两期分类误差会相乘，如果每期精度 85%，那么「未变化」的判定正确率理论上会降到约 72%。

因此分类后比较对单期分类精度要求极高，通常要求每期总体精度在 90% 以上，且类别体系严格一致。提升手段包括：用同一模型同一特征训练两期分类、做类别后处理平滑、对变化图斑做人工复核。

单期分类的方法选择与精度控制见 遥感影像分类 。分类后比较与直接差分可以互补：差分定位「哪里变了」，分类后比较回答「变成什么」。

6. 深度学习变化检测

深度变化检测把两期影像作为输入，端到端输出变化图。主流结构是孪生网络：两期共享权重的编码器分别提特征，再在特征层做差分或拼接，最后解码成变化图。

import torch
import torch.nn as nn

class SiameseChangeNet(nn.Module):
    def __init__(self, encoder, n_classes=2):
        super().__init__()
        self.encoder = encoder                       # 两期共享权重
        self.head = nn.Conv2d(encoder.out_channels * 2, n_classes, 1)

    def forward(self, t1, t2):
        f1 = self.encoder(t1)
        f2 = self.encoder(t2)
        feat = torch.cat([f1, f2, torch.abs(f1 - f2)], dim=1)   # 拼接加绝对差
        return self.head(feat)

几个关键设计点：

• 特征融合方式：拼接、绝对差、乘性交互各有侧重，绝对差对变化敏感，拼接保留更多信息。
• 时序注意力：用注意力让网络关注真正变化的区域，抑制物候引起的系统性差异。
• 多尺度监督：在多个分辨率上计算损失，兼顾大图斑与细碎变化。
• 变化先验：把差分图作为额外输入通道，给网络一个强先验。

Transformer 变化检测（如 ChangeFormer）用自注意力建模两期的长距离关系，在复杂场景下表现更好，代价是数据与算力需求更高。

深度方法的共同挑战是样本。变化样本天然稀疏，标注成本高。对策包括：用双期分类结果自动生成弱标签再人工修正、用合成变化增强、用预训练编码器降低数据需求。评价时除 mIoU 外还要看变化类的 F1，因为变化类占比低，mIoU 会被未变化类主导。

评价深度变化检测时要注意几点：变化类 F1 必须单独报告，因为未变化类占比常超过九成；边界像元可做容差处理，配准误差会让边界处天然不一致；模型要在跨区域、跨季节的数据上测试，避免只在同源数据上刷高分。

7. 影像时间序列与断点检测

双期方法只能看到两个快照，密集时间序列能捕捉变化的时刻与过程。Landsat 自 1984 年、Sentinel-2 自 2015 年起提供了可用的密集时序，是长时序监测的基础。

断点检测的核心是：给一条像元的时间序列，找出均值或趋势发生突变的时刻。

import numpy as np

def simple_breakpoint(series, dates, min_size=5):
    best_t, best_gain = None, 0.0
    total_var = np.nanvar(series) * len(series)
    for i in range(min_size, len(series) - min_size):
        left, right = series[:i], series[i:]
        var = np.nanvar(left) * len(left) + np.nanvar(right) * len(right)
        gain = total_var - var                       # 分割后方差下降越多越好
        if gain > best_gain:
            best_gain, best_t = gain, dates[i]
    return best_t, best_gain

断点检测的工程要点：

• 时间序列要先去噪与填补，云与缺失期不能直接参与统计。
• 要区分突变（火烧、砍伐）与渐变（退化、恢复），前者是阶跃，后者是趋势改变。
• 单个像元的序列噪声大，通常先用空间邻域或时序滤波稳定，再检测断点。
• 检测到的断点要做显著性检验，避免把噪声波动误判为变化。

时序去噪与填补是断点检测的前置步骤，常用线性插值填补缺失，用中值或 Savitzky-Golay 抑制噪声：

import pandas as pd

s = pd.Series(values, index=pd.to_datetime(dates))
s = s.interpolate(method="time", limit=3)             # 按时间插值，最多补 3 期
s = s.rolling(window=3, center=True, min_periods=1).median()   # 中值滤波去尖峰

时序数据量巨大，逐像元处理需要并行与分块。任务编排与调度可参考 数据管道编排 ，把像元块作为独立任务分发，是处理全国乃至全球尺度时序的标准做法。

8. 谐波回归与 BFAST 与 CCDC

真实地表的时间序列既有季节性周期，又有趋势与突变。谐波回归把序列分解为周期项与趋势项：

y(t) = a0 + a1 * t + sum_k [ c_k * cos(2*pi*k*t/T) + s_k * sin(2*pi*k*t/T) ] + e(t)

周期项刻画物候，趋势项刻画缓慢变化，残差刻画异常。检测变化时，看的是残差是否持续偏离，或者趋势项是否发生改变。

BFAST 把断点检测嵌入谐波模型：先拟合季节加趋势模型，再检测趋势或季节分量的断点，最后迭代精化。它输出的断点带时间与类型（趋势改变、季节改变），适合植被退化与恢复监测。

CCDC 用滑动窗口的谐波模型逐时相判断像元是否偏离模型预测，偏离超过阈值即标记为变化，然后用新数据重建模型。它天生适合连续监测，能给出变化的时刻而不只是「变了」。

这些方法共同的难点是参数多、对噪声敏感、计算量大。工程上要在精度与成本间权衡，比如用较粗的空间分辨率做区域筛查，再对候选区做精细分析。

CCDC 的窗口长度是核心参数：窗口太短模型不稳，太长对变化反应迟钝，常用 1.5 到 2 年的观测。判定阈值通常用残差标准差的倍数，或直接对残差做卡方检验。CCDC 会为每个像元维护一套模型与变化历史，存储与计算开销都不小，工程上要按需裁剪空间范围。

CCDC 单像元流程
初始化窗口拟合谐波模型 -> 逐期预测并检验残差 -> 连续 N 期超阈值则确认变化 -> 用新窗口重建模型

9. 变化图斑后处理与验证

原始变化检测结果通常充满碎斑与噪声，必须后处理。

from scipy import ndimage

labeled, n = ndimage.label(change)                   # 连通域标记
sizes = ndimage.sum(change, labeled, range(1, n + 1))
keep = np.isin(labeled, np.where(sizes >= 9)[0] + 1)  # 去掉小于 9 像元的碎斑
closed = ndimage.binary_closing(keep, iterations=1)   # 闭运算填补小洞

后处理要点：

• 最小图斑面积按应用设定，土地利用变化常取 0.5 到 1 公顷。
• 形态学闭运算填补内部孔洞，开运算去掉毛刺。
• 变化图斑往往沿地物边界出现细线状虚警，是配准误差的典型表现，可用边界掩膜抑制。
• 变化类型合并，把语义相近的类别归并，减少碎片。

验证要用独立的参考数据，抽样方法影响很大。推荐分层随机采样：按变化与未变化分层，变化区多采，未变化区少采，再按面积加权还原。指标除总体精度外，必须报告变化类的召回率、精确率与 F1，因为变化类占比低，总体精度会严重虚高。

验证样本要空间独立于训练与调参数据，且最好来自高分辨率影像或野外调查。把变化面积作为时间序列监控指标，可以发现系统性的假变化。

验证样本量的估算可以按目标精度反推：若希望变化类召回率的标准误不超过 3%，在该类占比较低时需要数百个变化样本，这在变化稀疏的区域意味着要抽样大量图斑。分层抽样能显著降低所需样本量，因为它在变化与未变化层内分别估计，避免被多数层稀释。

权衡取舍

• 双期 vs 时序：双期简单快速，只给两个快照；时序能定位变化时刻但数据与算力需求高。
• 差分 vs CVA：单波段差分简单但噪声敏感，CVA 用多波段联合更稳但需一致校正。
• IR-MAD vs 阈值法：IR-MAD 自动稳健但需足够未变化样本，阈值法简单但需人工标定。
• 分类后比较 vs 直接检测：前者给变化类型但误差叠加，后者定位准但无类型。
• 深度模型 vs 经典方法：深度精度上限高但样本需求大，经典方法在样本少时更实用。
• 阈值松紧：松则虚警多、复核成本高，紧则漏检多、变化被低估，需按业务代价定。

常见坑清单

• 配准误差未控制：现象是变化图斑沿边界呈细线，原因是两期错位，规避方法是配准到亚像元并做边界掩膜。
• 未做一致辐射校正：现象是整幅出现系统性变化，原因是大气与光照差异，规避方法是统一到地表反射率。
• 物候差异误判：现象是农田每年重复出现变化，原因是生育期不同，规避方法是同季比较或用时序模型。
• 阈值跨区照搬：现象是虚警或漏检，原因是阈值有地域性，规避方法是用本地样本标定。
• Otsu 在小变化占比下失效：现象是阈值落在噪声上，原因是直方图非双峰，规避方法是先粗筛再细判。
• 分类后比较精度虚高：现象是变化图斑满地，原因是两期分类误差叠加，规避方法是提升单期精度。
• 忽略云与阴影：现象是云区被判为变化，原因是未掩膜，规避方法是严格云掩膜。
• 最小图斑过小：现象是碎斑淹没真实变化，原因是未做连通域过滤，规避方法是按应用设最小面积。
• 验证样本不独立：现象是精度异常高，原因是样本泄漏，规避方法是空间分块并独立采样。
• 只看总体精度：现象是变化检测看似很准实则漏检，原因是变化类占比低，规避方法是报告变化类 F1。

小结

变化检测的成败，八成取决于前置数据处理，两成取决于算法。一致的辐射校正、亚像元配准、严格的云掩膜、物候对齐，这些准备工作做得越扎实，后面的算法越简单也能得到好结果。反之，再先进的模型也无法从被干扰污染的数据里还原真实地变。

方法选择上，双期任务用差分、CVA 或 IR-MAD，需要变化类型时叠加分类后比较；连续监测任务用谐波回归、BFAST 或 CCDC。深度方法在有足够标注时精度上限最高，但样本与算力是现实约束。无论哪种方法，后处理与独立验证都是不可省略的环节。

下一步建议把变化检测与单期分类、时序指数结合起来，形成「分类定类别、时序定时刻、变化检测定图斑」的组合流程；同时把验证与监控固化进流水线，让变化面积成为可持续观测的指标。