4、植被指数计算：NDVI、EVI、SAVI的原理与计算

各位同学，今天我们来聊聊植被指数。说实话，这玩意儿是我入行遥感时最先接触的东西，也是我后来做大宗商品预测时用得最多的特征之一。你想想看，农作物长得好不好，直接关系到期货价格走势，而植被指数就是那个能提前告诉你「地里啥情况」的指标。

我个人习惯把植被指数理解为「卫星给植物做的体检报告」。NDVI、EVI、SAVI这三个，就是最常用的三项检查指标。咱们一个一个来拆解。

4.1 NDVI：最经典的植被指数

NDVI全称是归一化植被指数。公式很简单：

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

其中NIR是近红外波段，Red是红波段。为什么这么算？因为健康植物对红光的吸收很强（用于光合作用），对近红外的反射很强（细胞结构散射）。一减一除，就把植被信息放大了。

我在项目中遇到过一个问题：用Landsat 8数据算NDVI时，发现城市区域的值偏高。后来一查，原来是建筑材料的反射特性跟植被有点像。嗯，这里要注意，NDVI在高植被覆盖区容易饱和，说白了就是到了0.8以上就分不清是森林还是灌木了。

4.2 EVI：增强型植被指数

EVI就是为了解决NDVI的饱和问题而生的。它加入了蓝光波段来校正大气影响，还调整了土壤背景的干扰。

EVI = G * (NIR - Red) / (NIR + C1 * Red - C2 * Blue + L)

参数说明：

• G = 2.5（增益系数）
• C1 = 6.0（红光校正系数）
• C2 = 7.5（蓝光校正系数）
• L = 1.0（土壤调节系数）

你看这个公式，比NDVI复杂了不少。但效果确实好，尤其是在高植被区和大气条件不好的时候。我记得有一次做巴西大豆产区的监测，雨季时云多雾重，NDVI数据全是噪声，换成EVI后信号就清晰多了。

4.3 SAVI：土壤调整植被指数

SAVI是专门为半干旱地区设计的。它的核心思想是：当植被覆盖度低时，土壤背景对NDVI的影响很大，需要引入一个土壤调节因子L。

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

L的取值规则：

植被覆盖度	L值	适用场景
高（>80%）	0	等同于NDVI
中等（50%-80%）	0.5	一般农田
低（<50%）	1.0	荒漠、稀疏草原

我曾经在内蒙古草原做项目时，发现NDVI值普遍偏低，但实地考察发现草长得还行。换成SAVI（L=0.5）后，数值就合理多了。说白了，SAVI就是给NDVI加了个「土壤滤镜」。

4.4 用numpy和rasterio实现

理论讲完了，咱们直接上代码。我个人习惯用rasterio读数据，用numpy做计算，又快又稳。

import numpy as np
import rasterio

def calculate_vegetation_indices(red_path, nir_path, blue_path=None, index_type='NDVI'):
    """
    计算植被指数
    red_path: 红波段文件路径
    nir_path: 近红外波段文件路径
    blue_path: 蓝波段文件路径（EVI需要）
    index_type: 'NDVI', 'EVI', 'SAVI'
    """
    with rasterio.open(red_path) as src_red:
        red = src_red.read(1).astype(np.float32)
        profile = src_red.profile
    
    with rasterio.open(nir_path) as src_nir:
        nir = src_nir.read(1).astype(np.float32)
    
    # 防止除零
    mask = (red + nir) == 0
    
    if index_type == 'NDVI':
        ndvi = (nir - red) / (nir + red + 1e-10)
        ndvi[mask] = -999
        return ndvi
    
    elif index_type == 'SAVI':
        L = 0.5  # 中等植被覆盖度
        savi = (nir - red) / (nir + red + L) * (1 + L)
        savi[mask] = -999
        return savi
    
    elif index_type == 'EVI':
        if blue_path is None:
            raise ValueError("EVI需要蓝波段数据")
        with rasterio.open(blue_path) as src_blue:
            blue = src_blue.read(1).astype(np.float32)
        
        G, C1, C2, L = 2.5, 6.0, 7.5, 1.0
        evi = G * (nir - red) / (nir + C1 * red - C2 * blue + L + 1e-10)
        evi[mask] = -999
        return evi

# 使用示例
ndvi_result = calculate_vegetation_indices('red.tif', 'nir.tif', index_type='NDVI')

# 保存结果
with rasterio.open('ndvi_output.tif', 'w', **profile) as dst:
    dst.write(ndvi_result, 1)

4.5 时间序列植被指数提取

做大宗商品预测，单张影像远远不够。我们需要的是时间序列——比如过去5年每16天的NDVI数据。

我的做法是三步走：

1. 批量读取：用glob或os.listdir遍历所有影像文件
2. 逐张计算：对每张影像调用上面的函数
3. 按区域聚合：提取研究区（比如某个县）的平均值

import glob
import pandas as pd
from shapely.geometry import box
import geopandas as gpd

def extract_ts_ndvi(image_folder, shapefile_path):
    """
    提取时间序列NDVI
    image_folder: 存放多时相影像的文件夹
    shapefile_path: 研究区边界矢量文件
    """
    # 读取研究区边界
    study_area = gpd.read_file(shapefile_path)
    
    # 获取所有影像文件（假设命名格式：YYYYMMDD_red.tif）
    red_files = sorted(glob.glob(f"{image_folder}/*_red.tif"))
    
    ts_data = []
    
    for red_file in red_files:
        # 从文件名提取日期
        date_str = red_file.split('/')[-1].split('_')[0]
        nir_file = red_file.replace('_red.tif', '_nir.tif')
        
        # 计算NDVI
        ndvi = calculate_vegetation_indices(red_file, nir_file, index_type='NDVI')
        
        # 用rasterio.mask提取研究区内的像元
        with rasterio.open(red_file) as src:
            out_image, out_transform = rasterio.mask.mask(src, study_area.geometry, crop=True)
            # 裁剪NDVI
            ndvi_cropped = ndvi[src.window_tl[0]:src.window_tl[0]+out_image.shape[1],
                                src.window_tl[1]:src.window_tl[1]+out_image.shape[2]]
        
        # 计算平均值（排除无效值）
        valid_pixels = ndvi_cropped[ndvi_cropped > -999]
        mean_ndvi = np.mean(valid_pixels) if len(valid_pixels) > 0 else np.nan
        
        ts_data.append({'date': date_str, 'ndvi': mean_ndvi})
    
    # 转为DataFrame
    df = pd.DataFrame(ts_data)
    df['date'] = pd.to_datetime(df['date'], format='%Y%m%d')
    df = df.sort_values('date')
    
    return df

# 使用示例
ts_df = extract_ts_ndvi('./landsat_images/', './study_area.shp')
print(ts_df.head())

4.6 知识体系总览

为了让大家更直观地理解本章的知识结构，我画了一张流程图：

这张图把整个流程串起来了：从原始影像到三种指数计算，再到时间序列提取，最后成为预测模型的输入特征。你想想看，每一步都有坑，但每一步也都有优化空间。

好了，植被指数这块就讲到这里。代码我已经在实际项目中跑过很多次了，你们直接拿去用就行。记住一点：没有最好的指数，只有最合适的指数。根据你的研究区和数据源，灵活选择才是王道。