4. 作物识别特征工程

前三节我们完成了影像预处理与样本准备，现在面对一个关键问题：把什么样的数值交给分类器，才能让不同作物在特征空间中真正分开。原始多光谱影像只有几个波段，单靠蓝、绿、红、近红外反射率直接分类，往往出现玉米与水稻混分、早稻与晚稻错判的情况。特征工程的任务，就是把原始波段变换、组合、统计成一组对作物物候和结构敏感的判别量，并从中筛掉冗余、保留有效信息。本节围绕光谱特征、植被指数、纹理特征、时序统计量四个维度展开，并完成一次可复现的特征提取与重要性分析。

四个维度的特征各自解决什么问题

光谱特征是基础层，直接使用各波段反射率或波段间的简单运算，如近红外与红光之差、短波红外与近红外之比。它保留了地物最原始的辐射信息，对区分裸土、水体、建筑与植被有效，但对同属植被的不同作物分辨力有限。植被指数通过波段组合放大植被的生理差异，NDVI 反映绿度与覆盖度，EVI 在高生物量区抗饱和，NDWI 对水分敏感因而能区分水稻田与旱作，SAVI 引入土壤调节因子适合苗期稀疏植被。不同指数对作物物候的响应曲线形状不同，这正是分类的信息来源。纹理特征描述像元邻域的灰度空间结构，常用灰度共生矩阵的均值、方差、对比度、同质性、熵等统计量，能刻画地块内部的垄作方向、种植密度与田块破碎度，对区分果园与农田、大田与菜地有补充作用。时序统计量则把一整季的多期观测压缩成少数几个数：NDVI 时序的最大值、最小值、均值、标准差、振幅、峰值出现时间、积分面积等。作物的播期、抽穗期、成熟期不同，时序曲线的形状和关键节点就不同，这是区分同谱异物的核心手段。

特征提取的具体实现

假设我们使用 Python 与 rasterio、numpy 处理一个县区的 6 期 Sentinel-2 影像，样本为 4 类作物：水稻、玉米、大豆、蔬菜，每类 200 个地块中心像元。提取流程分三步。第一步，逐期计算植被指数并堆叠为多波段栅格。以 NDVI 为例，公式为 (NIR−R)/(NIR+R)，计算时需先把反射率缩放到 0～1 并处理分母为零的情况：

import numpy as np, rasterio

def ndvi(nir, red, eps=1e-6):
    return (nir - red) / (nir + red + eps)

def lswi(nir, swir1, eps=1e-6):
    return (nir - swir1) / (nir + swir1 + eps)

# 逐期读取并计算，堆叠为 (期数, 高, 宽) 数组
ndvi_stack = np.stack([ndvi(nir_i, red_i) for nir_i, red_i in zip(nir_list, red_list)])

第二步，对每个样本点提取像元值，形成“样本×特征”表。光谱特征取 6 期各波段的均值，植被指数取 6 期 NDVI、EVI、NDWI、LSWI 的逐期值，纹理特征在近红外波段上用 5×5 窗口计算 GLCM 的均值、方差、对比度、同质性、熵共 5 个量，时序统计量则对 NDVI 序列计算最大值、最小值、均值、标准差、振幅（最大值减最小值）、峰值出现期次和时序积分面积。这样每个样本得到约 40 维候选特征。

第三步，做共线性初筛。计算特征两两之间的 Pearson 相关系数，若 |r| 大于 0.9，保留物理意义更明确或与作物物候更相关的一个。例如 6 期 NDVI 的均值与 EVI 均值高度相关，可只保留 EVI 均值；红光反射率与 NDVI 强负相关，可去掉红光反射率。初筛后通常剩 20～25 维。

特征重要性分析与结果判读

用随机森林在训练集上训练，并输出基于不纯度减少的重要性，再用验证集上的置换重要性交叉验证。一个合理的示例结果是：NDVI 峰值出现期次、NDVI 振幅、LSWI 均值排在前三，说明物候节律和水分状况是区分水稻与旱作的关键；GLCM 熵排在第五，说明田块内部结构对区分蔬菜与大豆有贡献；而单期红光反射率排在末位，属于可剔除的冗余特征。需要强调，这些排序是示例性假设，实际项目必须用本地样本重新训练得到，不能直接套用。

可直接套用的特征清单模板