植被指数与光谱指数计算

引言

光谱指数是遥感里性价比最高的分析手段之一。用两三个波段的线性或非线性组合，就能把植被长势、水体范围、燃烧痕迹这些地物属性从反射率中分离出来。它的计算成本极低，却常常比复杂模型更稳健，因为物理机理清晰、跨传感器可移植。

工程上的难点不在公式本身，而在前提条件。指数只有在经过辐射定标与大气校正之后的反射率上才有物理意义；直接拿 DN 值或表观反射率算 NDVI，跨时相比较会系统性偏移。饱和、土壤背景、传感器波段差异、异常值（水体、阴影、云边缘）都会让指数在真实场景里失效。

本文按「原理到落地」的顺序展开。先讲指数的数学本质与 NDVI 的推导，再覆盖 SAVI、EVI 这类针对土壤与大气改进的指数，然后扩展到水体、雪与燃烧指数，最后落到反射率前提、批量计算代码、时序合成与生物物理参数回归。

前置的辐射定标与大气校正见 辐射定标与大气校正 ，云污染处理见 云检测与云掩膜 。指数结果进一步用于变化检测与时序分析，是后续专题的输入。

目录

1. 光谱指数的数学本质
2. NDVI 的推导与局限
3. SAVI 与土壤调节
4. EVI 与饱和和大气改进
5. 水体与雪指数
6. 燃烧与灾害指数
7. 反射率前提与异常值
8. 批量计算与内存优化
9. 时序合成与生物物理回归

1. 光谱指数的数学本质

遥感波段记录的是地表在该波段区间的反射辐射。不同地物在不同波段的反射率差异，构成了识别它们的基础。健康植被的叶绿素强吸收红光、叶肉组织强反射近红外，这个「红低近红高」的对比就是植被指数的物理来源。

最朴素的表达是比值 NIR/RED。它确实能放大植被信号，但缺点是无界，值域随光照与传感器变化，阈值无法跨影像复用。归一化差分把它压缩到固定区间，同时抵消了乘性的光照与观测几何差异：

ND = (NIR - RED) / (NIR + RED)

分子体现对比，分母做归一化。只要地物反射率按比例缩放，比如不同太阳高度角下的整体亮度变化，分子分母同比例变化，结果近似不变。这就是归一化差分对光照鲁棒的来源，也是几乎所有指数都采用它的原因。

指数不是万能的分类器，它只把某一维物理属性投影到一个标量上。用之前要先问两件事：这一维属性是否足以区分目标地物，以及这个投影对噪声有多敏感。很多失败案例不是公式写错，而是选错了投影维度。

归一化差分之外的构造

归一化差分不是唯一的构造方式，常见的还有三类：

• 比值型：如 SR = NIR / RED，无界但动态范围大，适合做相对变化，不适合固定阈值。
• 距离型：如 BAI，衡量像元到某个参考点的距离，对特定地物敏感。
• 正交型：如 PVI，把植被信号投影到与土壤线正交的方向，最大程度剥离土壤亮度。

选择构造方式的原则是：目标地物的光谱差异落在哪个方向，就用哪个方向的投影。归一化差分适合差异体现为「此消彼长」的场景，比值型适合差异体现为「整体强弱」的场景。此外还要注意波段的物理宽度，宽波段与窄波段的同名指数不能直接互比。

2. NDVI 的推导与局限

NDVI 是最经典的植被指数，定义为近红外与红光的归一化差分：

NDVI = (NIR - RED) / (NIR + RED)

理论值域 -1 到 1。健康浓密植被约 0.7 到 0.9，稀疏植被 0.2 到 0.4，裸土 0.1 到 0.2，水体为负值，云与雪也常为负。这些经验区间只在反射率数据上成立，用在 DN 上会完全失真。

典型地物的 NDVI 经验区间如下，可用于快速目视判读与质检：

import numpy as np

def ndvi(nir, red, eps=1e-6):
    nir = nir.astype("float32")                      # 避免整数除法截断
    red = red.astype("float32")
    return (nir - red) / (nir + red + eps)           # eps 防止除零

NDVI 有三个绕不开的局限，理解它们才能选对替代指数。

第一，饱和。当叶面积指数 LAI 超过 3 到 4，红光吸收接近饱和，NDVI 增长趋缓，对高生物量区不敏感，浓密森林内部区分度低。

第二，土壤背景。稀疏植被下冠层未完全覆盖，传感器同时看到植被与土壤，NDVI 随土壤亮度变化而漂移，同一植被在不同土壤上读数不同。

第三，大气与几何。气溶胶、水汽、观测天顶角都会引入偏差，在跨时相比较时尤为明显，这也是长时序监测必须先做大气校正的原因。

饱和与土壤问题催生了后续的改进指数，大气问题则要靠校正或专门指数缓解。

3. SAVI 与土壤调节

土壤调节植被指数 SAVI 引入土壤亮度调节因子 L，压低土壤背景的干扰：

SAVI = (NIR - RED) / (NIR + RED + L) * (1 + L)

L 的取值与植被覆盖度相关，经验配置如下：

当 L 取 0 时 SAVI 退化为 NDVI。实践中 L 难以先验已知，于是有了 OSAVI，固定 L 为 0.16，牺牲一点精度换取无需调参：

OSAVI = (NIR - RED) / (NIR + RED + 0.16)

MSAVI 更进一步，让 L 随 NDVI 自适应变化，无需人工设定。选择上，如果研究区覆盖度变化大，MSAVI 更稳；如果只是做长时序一致性监测，OSAVI 的固定参数更利于跨期可比。SAVI 族的共性是牺牲部分绝对精度，换取对土壤背景的鲁棒，因此在稀疏植被区（干旱区、早期作物）收益最大，在浓密植被区与 NDVI 差别不大。

def savi(nir, red, L=0.5, eps=1e-6):
    nir = nir.astype("float32")
    red = red.astype("float32")
    return (nir - red) / (nir + red + L) * (1.0 + L)     # L 按覆盖度设定

4. EVI 与饱和和大气改进

增强植被指数 EVI 同时针对饱和与大气做了改进：

EVI = G * (NIR - RED) / (NIR + C1 * RED - C2 * BLUE + L)

标准参数为 G=2.5、C1=6.0、C2=7.5、L=1.0。它的三个设计点值得理解：

• 分母加入 BLUE 波段，用蓝光估计气溶胶散射，做部分大气自校正，因此对烟霾更鲁棒。
• 分母中的 RED 项系数放大，配合 L 项，让高植被区不再过早饱和，动态范围更大。
• G 是增益系数，把 EVI 的值域重新拉回接近 NDVI 的尺度，便于沿用经验阈值。

代价是 EVI 需要蓝光波段，对蓝光波段的信噪比要求高；在蓝光噪声大或缺失的传感器上不可用。而且 EVI 对残留云与云边缘更敏感，容易出现异常高值，必须配合严格的云掩膜。

def evi(blue, red, nir, G=2.5, C1=6.0, C2=7.5, L=1.0, eps=1e-6):
    denom = nir + C1 * red - C2 * blue + L
    out = G * (nir - red) / (denom + eps)
    return np.clip(out, -1.0, 1.0)                   # 抑制云边缘产生的异常值

EVI 与 NDVI 的关系是互补：NDVI 在低覆盖区更敏感，EVI 在高覆盖区更敏感。做长时序分析时，常见做法是两个都算，用 EVI 看浓密植被，用 NDVI 做交叉验证，两者出现系统性背离时往往是云掩膜或定标出了问题。

实际项目里不必二选一。把 NDVI、EVI、SAVI 一起算出来，作为特征堆叠输入下游模型，往往比纠结单一指数更有效。

5. 水体与雪指数

水体在近红外强吸收，反射率极低，这让「绿光高、近红外低」成为稳定的水体特征。McFeeters 提出的 NDWI 用绿光与近红外：

NDWI = (GREEN - NIR) / (GREEN + NIR)

但 NDWI 在建成区容易误判，因为建筑在绿光与近红外的对比与水体相似。Xu 提出的 MNDWI 用中红外替换近红外，显著抑制建筑干扰：

MNDWI = (GREEN - SWIR1) / (GREEN + SWIR1)

雪与云在可见光都很亮，区分靠短波红外：雪在 SWIR 强吸收，云在 SWIR 仍较亮。归一化雪指数 NDSI：

NDSI = (GREEN - SWIR1) / (GREEN + SWIR1)

NDSI 大于 0.4 通常判为雪，但要注意与水体区分，水体 NDSI 也可能偏正。工程上常把 NDSI 与 NDVI 组合，先排除植被区，再叠加地形阴影掩膜。

水体指数的阈值有强地域性。浑浊河流与清澈湖泊的最优阈值不同，工程上应抽取本地样本重新标定，而不是全局照搬 0。

阈值不是固定的。同一景影像里，清澈深水的 NDWI 可达 0.5，而浑浊浅水可能只有 0.1。工程做法是先用 Otsu 或双峰直方图自动确定初始阈值，再用少量人工样本微调，最后固化成本地参数。

水体提取后还要做后处理：去掉面积小于 8 个像元的小斑块，填补内部孔洞，平滑边界。这一步能显著改善面积统计的稳定性，避免碎斑导致的面积高估。掩膜边界像元常处于混合像元状态，边界处的指数介于水陆之间，最好用亚像元或软分类方法处理。

6. 燃烧与灾害指数

归一化燃烧比 NBR 用近红外与短波红外，对燃烧迹地敏感：

NBR = (NIR - SWIR2) / (NIR + SWIR2)

火后植被受损，NIR 下降、SWIR 上升，NBR 显著降低。前后两期相减得到差分燃烧指数：

dNBR = NBR_pre - NBR_post

dNBR 的分级被 USGS 用于烧伤严重度制图，经验阈值大致为：小于 0.1 为未燃烧，0.1 到 0.27 为低，0.27 到 0.66 为中低，更高为高严重度。这些阈值需要在本地标定，不能直接照搬。

另一个常用的是 BAI，用红光与近红外的距离形式，对近期火烧迹地更敏感：

BAI = 1 / ((0.1 - RED)^2 + (0.06 - NIR)^2)

燃烧指数的关键是前后两期的可比性：必须使用同季节、同传感器、经过一致校正的反射率，否则物候差异会淹没燃烧信号。火后恢复监测通常按 NBR 的时序轨迹建模，这与时间序列断点检测是同一套思路。

火后恢复可以用 NBR 的时序轨迹刻画：火烧当年骤降，之后逐年回升，通常用恢复年限或 dNBR 的衰减曲线来量化。

dnbr = nbr_pre - nbr_post
severity = np.digitize(dnbr, [0.1, 0.27, 0.66])          # 分四级严重度

分级结果要与实地调查的烧伤程度交叉验证，因为不同植被类型对 NBR 下降的响应幅度不同，同一 dNBR 在针叶林与草地的含义并不相同。

7. 反射率前提与异常值

指数有效性的第一前提是输入为地表反射率。数据层级大致分三级：

从 DN 到反射率需要辐射定标系数。以 Sentinel-2 L2A 为例，2022 年 1 月起基线 04.00 引入了 1000 的偏移量，处理时必须先减掉，否则所有指数整体偏移：

boa = (dn.astype("float32") - 1000.0) / 10000.0      # L2A 基线 04.00 的定标

异常值处理同样关键。水体、阴影、云边缘会让分母趋近零，指数爆到极端值。工程上常用三招：

• 分母加小量 eps，避免除零。
• 对结果做裁剪，把值域限制在物理合理区间。
• 用有效像元掩膜，把云、水、无效值排除后再统计。

valid = (nir > 0) & (red > 0) & (~cloud_mask)       # 有效像元掩膜
ndvi_clean = np.where(valid, ndvi_arr, np.nan)      # 无效置 NaN

统计时用 np.nanmean 而不是 np.mean，否则 NaN 会污染整幅统计。这一点在区域统计与长时序聚合时最容易踩坑。

无效像元掩膜的来源

有效像元掩膜通常来自三层叠加：

1. 传感器自带的质量波段，如 Sentinel-2 的 SCL 场景分类层，直接给出云、阴影、雪、水的标签。
2. 指数自身的物理约束，如 NIR 与 RED 同时大于 0，排除填充值与越界值。
3. 外部产品，如 Landsat 的 QA_PIXEL 位掩码，需要按位解析。

scl = src.read(1)
bad = np.isin(scl, [0, 1, 3, 8, 9, 10, 11])             # 无效、云影、云、卷云、雪
valid = (~bad) & (nir > 0) & (red > 0)

掩膜宁严勿宽。多掩掉一些有效像元，代价是统计样本减少；少掩掉一个云边缘，代价是整幅统计被异常值拉偏，后者更难事后发现。

8. 批量计算与内存优化

真实作业面对的是成百上千景影像，逐像元 Python 循环不可行。rasterio 的分块窗口读取配合 NumPy 向量化是标准做法。

import rasterio
from rasterio.windows import Window

def calc_ndvi_blockwise(src_nir, src_red, dst, block=1024):
    profile = src_nir.profile
    profile.update(dtype="float32", count=1, nodata=np.nan, compress="deflate")
    with rasterio.open(dst, "w", **profile) as out:
        for row in range(0, src_nir.height, block):
            for col in range(0, src_nir.width, block):
                win = Window(col, row, block, block)
                nir = src_nir.read(1, window=win)
                red = src_red.read(1, window=win)
                out.write(ndvi(nir, red), 1, window=win)   # 逐块写出，常驻内存恒定

几个工程要点：

• 分块大小要匹配磁盘与内存，1024×1024 的 float32 块约 4 MB，适合多数场景。
• 输入波段必须对齐，同一景的 NIR 与 RED 网格一致，跨景则要先重采样到统一网格。
• 输出用 float32 而非 float64，指数精度足够且体积减半。
• 用 deflate 压缩，指数影像空间相关性高，压缩率通常可观。

更大规模时，把每景的指数计算包装成独立任务，用队列并行调度，把 IO 密集与计算密集分开。云掩膜要在指数计算之前完成，否则无效像元会参与统计并拉偏结果。

重采样与网格对齐

跨传感器或跨景计算时，先把所有输入重采样到统一网格。推荐用最近邻处理分类与掩膜，用双线性或三次卷积处理连续的反射率，避免在反射率上引入插值伪影。重采样后的网格应与目标分析网格一致，否则后续像元级运算会错位。

gdalwarp -tr 10 10 -t_srs EPSG:32650 -r bilinear \
  -of GTiff -co COMPRESS=DEFLATE s2_b04.tif b04_10m.tif   # 统一到 10 米与 UTM 50N

并行调度上，按景切分任务，每个任务只读自己需要的窗口，输出独立的指数文件，最后用虚拟栅格或 mosaic 拼接。这样单任务内存可控，失败可重试，整体吞吐由并发度而非单机内存决定。

9. 时序合成与生物物理回归

单期指数受云与观测条件影响大，时序分析前要先做合成。最常用的是最大值合成法 MVC：在合成窗口内逐像元取最大值。

stack = np.stack([ndvi_t1, ndvi_t2, ndvi_t3], axis=0)   # 时间维在最前
composite = np.nanmax(stack, axis=0)                    # 逐像元取最大

MVC 的依据是云与阴影会拉低 NDVI，取最大即取最晴空、最接近真实植被的观测。窗口通常取 8 天或 16 天。缺点是会系统性偏高，且对下降趋势不敏感，作物收割这类快速下降会被最大值掩盖。

指数最终常要回归到生物物理参数。以 LAI 为例，经验关系多为指数或幂函数：

lai = -np.log((0.69 - ndvi) / 0.59) / 0.91              # 经验回归，仅作示例

这类回归必须用地面实测数据本地标定，且要注意饱和区间的反演误差放大：NDVI 在 0.8 以上的微小噪声会导致 LAI 反演剧烈波动。更严格的路线是辐射传输模型反演，如 PROSAIL，但计算成本高得多，通常只用于关键时相或验证样本。

时序平滑与谐波拟合

MVC 之后的序列仍有残留噪声，常用 Savitzky-Golay 滤波平滑，它在去噪的同时保留物候曲线的峰谷形态：

from scipy.signal import savgol_filter

smooth = savgol_filter(ndvi_series, window_length=7, polyorder=2)   # 7 期窗口，二次多项式

更结构化的做法是谐波拟合，把一年的 NDVI 序列分解为均值和若干阶正弦分量：

NDVI(t) = a0 + sum_k [ a_k * cos(2*pi*k*t/T) + b_k * sin(2*pi*k*t/T) ]

通常取 2 到 3 阶即可捕捉一年两熟的双峰或一年一熟的单峰。拟合得到的振幅、相位、残差分别对应植被强度、物候时间与异常扰动，是后续断点检测的输入特征。拟合前要保证时间采样近似均匀，缺失期需要插值补足，否则谐波系数会有偏。

权衡取舍

• 指数选择：NDVI 通用但饱和，EVI 抗饱和但要蓝光，SAVI 抗土壤但要调 L，按覆盖度与数据条件选。
• 数据层级：TOA 计算快、无需大气产品，但跨期不可比；BOA 才是指数的正确输入，代价是校正链路更长。
• 计算精度：float32 足够，float64 只浪费内存与带宽，对最终统计无实质改善。
• 合成窗口：窗口越长云污染越少，但对快速变化越迟钝，需按物候节奏权衡。
• 异常值策略：裁剪会掩盖真实极端，掩膜更严谨但需要可靠的云与水体产品。
• 经验回归与物理模型：回归实现简单但泛化差，模型反演更准但需要参数与算力。

常见坑清单

• 用 DN 直接算 NDVI：现象是值域远超 -1 到 1，原因是未定标，规避方法是先转反射率。
• 忘记 Sentinel-2 基线偏移：现象是植被指数整体偏低，原因是未减 1000 偏移，规避方法是按基线版本处理。
• 分母除零未防护：现象是出现 Inf 或 NaN 大片，原因是水体与阴影使分母趋零，规避方法是加 eps。
• 用 np.mean 统计：现象是整幅均值变 NaN，原因是无效像元置了 NaN，规避方法是改用 np.nanmean。
• 忽略云掩膜：现象是 EVI 出现异常高值，原因是云边缘散射，规避方法是严格掩膜后再统计。
• 跨传感器直接比 NDVI：现象是同一地块两传感器差异明显，原因是波段响应函数不同，规避方法是做交叉标定。
• MVC 窗口过长：现象是作物收割期信号滞后，原因是最大值压制了下降，规避方法是缩短窗口或改用时序模型。
• SAVI 的 L 照搬：现象是稀疏区反而更差，原因是 L 与覆盖度不匹配，规避方法是按覆盖度选或改用 OSAVI。
• 输出用 float64：现象是文件体积翻倍、IO 变慢，原因是精度过剩，规避方法是统一 float32。
• 阈值跨区照搬：现象是水体或雪判错，原因是阈值有地域性，规避方法是用本地样本重新标定。

小结

光谱指数的价值在于把多维反射率压缩成可解释、可比较的标量，代价是每一维投影都伴随信息损失与前提假设。理解归一化差分的鲁棒性来源、饱和与土壤背景的物理成因、以及反射率层级对可比性的影响，是正确使用指数的前提。

工程落地时，把定标、掩膜、异常值处理、分块计算、时序合成串成一条稳定流水线，比追求某个「最优指数」更重要。指数只是中间产物，最终要服务于分类、变化检测或生物物理反演这些下游任务。

下一步建议从指数的时序形态入手，理解断点检测与谐波回归如何把 NDVI 序列变成物候与扰动的证据，见 变化检测与影像时间序列分析 ；也可对照 遥感影像分类 看指数如何作为特征进入分类器。