4、时序特征提取：植被指数（NDVI、EVI、LAI）计算、物候特征提取（SOS、EOS、LOS）
好，咱们进入正题。时序特征提取这块，说白了就是把卫星影像里那些像素级别的数字，翻译成我们能理解的生态语言。你拿到一堆时间序列的遥感影像，不能光盯着灰度值看，得算出有物理意义的指标来。
我个人习惯把这一步分成两个阶段：植被指数计算和物候参数提取。前者是基础，后者是应用。咱们一个一个来拆解。
4.1 植被指数计算：NDVI、EVI、LAI
植被指数，本质上是波段间的数学组合。目的是增强植被信号，抑制土壤、大气这些干扰项。我最早做这个的时候，以为公式越复杂越好，后来发现，简单的东西往往最可靠。
4.1.1 NDVI（归一化植被指数）
NDVI 是遥感界的「Hello World」。公式简单到令人发指：
NDVI = (NIR - Red) / (NIR + Red)其中 NIR 是近红外波段，Red 是红波段。对于 Landsat 8 来说，就是 Band 5 和 Band 4。
我在项目中遇到过一个问题：直接用原始 DN 值算 NDVI，结果全是负值。后来才意识到，必须先做辐射定标和大气校正。你想想看，DN 值本身没有物理意义，不校正就是瞎算。
核心要点：
- NDVI 范围：-1 到 1。健康植被通常在 0.3 到 0.8 之间。
- 水体、雪、云：通常为负值或接近 0。
- 裸土、建筑：接近 0。
下面是我常用的 Python 计算片段，用 numpy 向量化操作，速度很快：
import numpy as np
def calculate_ndvi(nir_band, red_band):
"""
计算 NDVI
nir_band, red_band: 2D numpy arrays
"""
# 防止除零
denominator = nir_band + red_band
denominator = np.where(denominator == 0, 0.0001, denominator)
ndvi = (nir_band - red_band) / denominator
# 裁剪到有效范围
ndvi = np.clip(ndvi, -1, 1)
return ndvi避坑指南：我曾经直接用 float32 算，结果在大面积水体区域出现 NaN。后来加了 np.clip 和分母保护，问题解决。另外，输出时建议乘以 10000 存成 int16，能省不少磁盘空间。
4.1.2 EVI（增强型植被指数）
NDVI 在高植被覆盖区容易饱和，说白了就是到了 0.8 以上就分不出差别了。EVI 就是为了解决这个问题。它还加入了蓝波段来校正大气影响。
EVI = 2.5 * (NIR - Red) / (NIR + 6 * Red - 7.5 * Blue + 1)这里的系数 2.5、6、7.5、1 是经验值，来自 MODIS 的算法。我建议你不要随便改，除非你有充分的物理依据。
嗯，这里要注意：EVI 对大气校正的质量要求更高。如果你用的是 Level-1 数据（未做大气校正），算出来的 EVI 可能还不如 NDVI 靠谱。
| 指数 | 优势 | 劣势 | 适用场景 |
|---|---|---|---|
| NDVI | 简单、稳定、计算快 | 高植被区饱和 | 低-中植被覆盖 |
| EVI | 抗饱和、抗大气 | 需要蓝波段、计算复杂 | 高植被覆盖、热带雨林 |
4.1.3 LAI（叶面积指数）
LAI 不是直接算出来的，而是反演出来的。它定义为单位地表面积上叶片总面积的一半。LAI 和 NDVI/EVI 有很强的经验关系，但更精确的做法是用辐射传输模型。
我常用的经验公式是这样的：
LAI = 0.57 * exp(2.33 * NDVI)这个公式来自某个经典文献，适用于草地和农田。对于森林，系数要调整。我在做东北森林项目时，发现这个公式高估了针叶林的 LAI，后来改用基于 EVI 的公式才校准过来。
重要提醒：LAI 反演没有万能公式。不同植被类型、不同传感器、不同大气条件，都需要重新标定。如果你追求精度，建议用 PROSAIL 模型做查找表反演，虽然慢，但靠谱。
4.2 物候特征提取：SOS、EOS、LOS
算完植被指数，我们得到的是每个时间点的「快照」。但生态学家关心的是过程：什么时候开始变绿？什么时候开始变黄？生长季有多长？
这就是物候特征提取要干的事。三个核心参数：
- SOS（Start of Season）：生长季开始日期
- EOS（End of Season）：生长季结束日期
- LOS（Length of Season）：生长季长度，LOS = EOS - SOS
4.2.1 时间序列重构
在提取物候参数之前，必须先做一件事：平滑去噪。原始 NDVI 时间序列充满了噪声——云、阴影、气溶胶都会造成异常值。
我常用的方法有三种：
- S-G 滤波（Savitzky-Golay）：保留细节，适合高时间分辨率数据。
- 双逻辑斯蒂拟合（Double Logistic）：物候学意义明确，参数可解释。
- Whittaker 平滑：速度快，适合大数据量。
我个人习惯用双逻辑斯蒂拟合。为什么？因为它拟合出来的曲线参数直接对应 SOS 和 EOS，不需要再额外做阈值判断。
from scipy.optimize import curve_fit
import numpy as np
def double_logistic(t, a, b, c, d, e, f):
"""
双逻辑斯蒂函数
t: 时间（儒略日）
a: 背景值（冬季NDVI）
b: 振幅（生长季峰值-背景值）
c: SOS相关参数
d: SOS相关参数
e: EOS相关参数
f: EOS相关参数
"""
return a + b * (1/(1+np.exp(-c*(t-d))) - 1/(1+np.exp(-e*(t-f))))
# 拟合示例（伪代码）
# popt, _ = curve_fit(double_logistic, doy, ndvi_smoothed, p0=initial_guess)经验之谈：初始参数 p0 很重要。我一般用冬季 NDVI 均值作为 a，用峰值减去冬季值作为 b。c 和 e 初始设为 0.1，d 设为 150（约 5 月底），f 设为 250（约 9 月初）。这样拟合收敛快，不容易掉进局部最优。
4.2.2 阈值法提取 SOS/EOS
如果你不想用拟合方法，阈值法是最直观的。常用的阈值有：
- 固定阈值：比如 NDVI > 0.2 就算进入生长季。简单粗暴，但不同区域不通用。
- 动态阈值：比如 NDVI 达到当年振幅的 20% 就算 SOS。适应性更好。
动态阈值的计算步骤：
- 计算每年 NDVI 的振幅：Amplitude = max(NDVI) - min(NDVI)
- 设定阈值比例：比如 20%
- SOS = 第一个 NDVI 超过 (min + 0.2 * Amplitude) 的日期
- EOS = 最后一个 NDVI 低于 (min + 0.2 * Amplitude) 的日期
注意：阈值法对噪声敏感。我建议在提取前先做一次中值滤波（窗口大小 3-5 个时间点），能有效减少误判。
4.2.3 物候参数的质量控制
提取出来的 SOS 和 EOS 不一定都合理。我一般会做以下过滤：
- SOS 必须在 1 到 180 天之间（北半球）
- EOS 必须在 180 到 365 天之间
- LOS 必须大于 30 天（否则可能是噪声）
- 峰值 NDVI 必须大于 0.3（否则可能是裸土或水体）
我曾经在青藏高原的项目中，发现有些像元的 SOS 算出来是第 5 天——1 月份就进入生长季了？这显然不合理。后来排查发现是雪盖造成的 NDVI 异常。加上积雪掩膜后，数据才正常。
4.3 本章知识体系
下面这张图概括了本章的核心逻辑。从原始影像到物候参数，每一步都有坑，每一步也都有技巧。
这张图把整个流程串起来了。你从最上面的原始影像开始，一路往下走，每一步都有对应的算法和工具。我个人建议你先把 NDVI 算熟，再碰 EVI 和 LAI。物候提取这块，先从阈值法入手，理解原理后再上拟合方法。
好了，这一章的内容就到这里。代码和公式都给了，剩下的就是动手跑一遍。记住，遥感分析不是玄学，是工程。每一步都要问自己：这个参数合理吗？这个结果能解释吗？