遥感新手必看：用Python+ENVI快速识别植被、水体、裸土（附光谱曲线对比图）

遥感影像地物分类实战：从ENVI光谱分析到Python自动化处理

当第一次接触遥感影像时，面对那些五彩斑斓的像素点，很多初学者会感到无从下手。实际上，每一类地物都有其独特的光谱"指纹"，就像超市商品的条形码一样，通过分析这些特征，我们就能让计算机自动识别出植被、水体、裸土等地表覆盖类型。本文将带你用ENVI和Python这两个工具，从零开始掌握遥感地物分类的实用技能。

1. 准备工作与环境配置

在开始实际操作前，我们需要准备好软硬件环境和数据资源。对于遥感处理来说，合适的数据和工具是成功的一半。

1.1 数据获取与预处理

推荐使用Landsat 8/9或Sentinel-2卫星数据，这些数据免费且容易获取。美国地质调查局(USGS)的EarthExplorer和欧空局的Copernicus Open Access Hub是两个主要的数据源平台。

下载数据时需要注意：

• 选择云量低于10%的影像
• 优先选择生长季(植被茂盛时期)的数据
• 确保研究区域没有大气污染等干扰因素

数据下载后通常需要进行以下预处理步骤：

1. 辐射定标：将DN值转换为辐射亮度或反射率
2. 大气校正：消除大气散射和吸收的影响
3. 影像裁剪：只保留研究区域，减少计算量

1.2 软件安装与配置

ENVI是遥感领域最常用的专业软件之一，其光谱分析功能尤为强大。同时，我们需要配置Python环境来处理自动化任务：

# 创建Python虚拟环境 python -m venv rs_env source rs_env/bin/activate # Linux/Mac rs_env\Scripts\activate # Windows # 安装必要的Python库 pip install numpy rasterio scikit-learn matplotlib jupyterlab

对于ENVI，建议安装最新版本，并确保Python接口配置正确。ENVI的Python API可以让我们在Python环境中调用ENVI的功能，实现更灵活的自动化处理。

2. ENVI中的光谱特征分析

理解不同地物的光谱特征是准确分类的基础。ENVI提供了强大的光谱分析工具，让我们能够直观地观察和比较各类地物的反射特性。

2.1 典型地物的光谱曲线特征

在ENVI中加载影像后，可以通过以下步骤查看光谱曲线：

1. 打开Display窗口并加载影像
2. 使用Spectral Profile工具
3. 在影像上点击感兴趣的区域

不同地物的典型光谱特征如下表所示：

提示：在实际分析时，建议先选择纯净的"训练样本"区域，避免混合像元影响光谱特征的典型性。

2.2 光谱特征的可视化比较

ENVI的Spectral Library功能允许我们将不同地物的光谱曲线叠加显示，便于直观比较。例如，健康植被和水体的光谱曲线在近红外波段差异非常明显：

# 在Python中使用rasterio绘制光谱曲线示例 import rasterio import matplotlib.pyplot as plt with rasterio.open('landsat.tif') as src: # 读取不同地物的像元值 vegetation = src.read()[:, 100, 200] # 植被样本点 water = src.read()[:, 150, 300] # 水体样本点 bands = range(1, src.count+1) plt.plot(bands, vegetation, 'g-', label='Vegetation') plt.plot(bands, water, 'b-', label='Water') plt.xlabel('Band Number') plt.ylabel('Reflectance') plt.legend() plt.show()

这种可视化对比能帮助我们快速理解各类地物的光谱差异，为后续的分类算法选择合适的光谱特征。

3. Python自动化分类实现

有了对光谱特征的理解后，我们可以用Python编写自动化脚本实现地物分类。这种方法比手动操作更高效，且可重复使用。

3.1 基于阈值的简单分类

最简单的分类方法是设定各波段的阈值范围：

import numpy as np import rasterio def simple_threshold_classification(input_img, output_img): with rasterio.open(input_img) as src: # 读取近红外和红光波段 nir = src.read(5) # Landsat 8近红外波段 red = src.read(4) # Landsat 8红光波段 # 计算NDVI(归一化植被指数) ndvi = (nir - red) / (nir + red + 1e-10) # 定义分类规则 classified = np.zeros_like(red, dtype=np.uint8) classified[ndvi > 0.6] = 1 # 植被 classified[(nir < 0.1) & (red < 0.1)] = 2 # 水体 classified[(ndvi <= 0.6) & (ndvi > 0.2)] = 3 # 裸土 # 保存结果 profile = src.profile profile.update(dtype=rasterio.uint8, count=1) with rasterio.open(output_img, 'w', **profile) as dst: dst.write(classified, 1) # 使用示例 simple_threshold_classification('landsat.tif', 'classified.tif')

这种方法虽然简单，但对于初步分类和快速评估非常有效。主要优点是计算速度快，易于理解和调整。

3.2 机器学习分类方法

对于更复杂的场景，可以使用机器学习算法提高分类精度。scikit-learn提供了多种分类算法：

from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split import numpy as np def machine_learning_classification(img_path, training_data): """ img_path: 输入影像路径 training_data: 训练样本，格式为[(x,y,class), ...] """ with rasterio.open(img_path) as src: # 准备训练数据 X = [] y = [] for x, y, cls in training_data: spectra = src.sample((x, y)) X.append(spectra[0]) y.append(cls) X = np.array(X) y = np.array(y) # 划分训练集和测试集 X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.3, random_state=42) # 训练随机森林分类器 clf = RandomForestClassifier(n_estimators=100, random_state=42) clf.fit(X_train, y_train) # 评估模型 accuracy = clf.score(X_test, y_test) print(f"模型准确率: {accuracy:.2f}") # 对整个影像进行分类 data = src.read() height, width = data.shape[1], data.shape[2] reshaped_data = data.reshape(data.shape[0], -1).T predicted = clf.predict(reshaped_data) classified = predicted.reshape(height, width) # 保存结果 profile = src.profile profile.update(dtype=rasterio.uint8, count=1) with rasterio.open('ml_classified.tif', 'w', **profile) as dst: dst.write(classified, 1) # 示例训练数据 (x坐标, y坐标, 类别) # 类别: 1=植被, 2=水体, 3=裸土 training_samples = [ (100, 200, 1), (150, 300, 2), (200, 150, 3), (120, 180, 1), (160, 320, 2), (180, 140, 3) ] machine_learning_classification('landsat.tif', training_samples)

机器学习方法的优势在于能够自动学习复杂的光谱特征组合，适应不同类型的地物和场景变化。

4. 分类结果验证与优化

分类完成后，我们需要评估结果的准确性并找出可能的改进方向。

4.1 精度评估方法

常用的精度评估指标包括：

• 混淆矩阵：显示各类别被正确和错误分类的情况
• 总体精度：所有正确分类的像元占总像元的比例
• Kappa系数：考虑随机因素的分类精度指标
• 生产者精度：某类别被正确分类的比例
• 用户精度：分类结果中某类别实际正确的比例

from sklearn.metrics import confusion_matrix, classification_report # 假设我们有验证数据 y_true = [1,1,2,2,3,3,1,2,3] # 实际类别 y_pred = [1,1,2,3,3,2,1,2,3] # 预测类别 print("混淆矩阵:") print(confusion_matrix(y_true, y_pred)) print("\n分类报告:") print(classification_report(y_true, y_pred))

4.2 常见问题与解决方案

在实际分类过程中，经常会遇到一些典型问题：

1. 混合像元问题：

   • 现象：一个像元包含多种地物类型
   • 解决方案：使用更高分辨率的影像或亚像元分类方法

2. 同物异谱/同谱异物：

   • 现象：同类地物光谱不同或不同类地物光谱相似
   • 解决方案：增加辅助数据(如DEM、纹理特征)或使用更复杂的分类器

3. 季节性变化影响：

   • 现象：同一地物在不同季节光谱特征不同
   • 解决方案：使用多时相数据或季节性特征库

4. 大气条件干扰：

   • 现象：气溶胶、水汽等影响光谱特征
   • 解决方案：严格的大气校正处理

4.3 分类结果后处理

原始分类结果往往存在噪声和孤立像元，可以通过后处理改善：

from scipy.ndimage import median_filter def postprocess_classification(input_path, output_path): with rasterio.open(input_path) as src: classified = src.read(1) # 中值滤波去除噪声 filtered = median_filter(classified, size=3) # 保存结果 profile = src.profile with rasterio.open(output_path, 'w', **profile) as dst: dst.write(filtered, 1) postprocess_classification('classified.tif', 'filtered.tif')

后处理可以显著改善分类结果的可视化效果和使用价值，但要注意不要过度平滑而丢失重要细节。

5. 实际应用案例与进阶技巧

掌握了基本分类方法后，我们可以将这些技术应用到各种实际问题中，并根据特定需求进行优化。

5.1 植被健康监测

通过多时相分类结果比较，可以监测植被变化：

def detect_vegetation_change(img1_path, img2_path, output_path): with rasterio.open(img1_path) as src1, rasterio.open(img2_path) as src2: # 计算两期影像的NDVI nir1, red1 = src1.read(5), src1.read(4) ndvi1 = (nir1 - red1) / (nir1 + red1 + 1e-10) nir2, red2 = src2.read(5), src2.read(4) ndvi2 = (nir2 - red2) / (nir2 + red2 + 1e-10) # 检测显著变化区域 change = np.zeros_like(ndvi1, dtype=np.int8) change[ndvi1 - ndvi2 > 0.2] = 1 # 植被减少 change[ndvi2 - ndvi1 > 0.2] = 2 # 植被增加 # 保存结果 profile = src1.profile profile.update(dtype=rasterio.int8, count=1) with rasterio.open(output_path, 'w', **profile) as dst: dst.write(change, 1)

5.2 水体提取优化

针对水体提取，可以结合多个水体指数提高精度：

def enhanced_water_extraction(img_path, output_path): with rasterio.open(img_path) as src: green = src.read(3) # 绿光波段 nir = src.read(5) # 近红外波段 swir1 = src.read(6) # 短波红外1 # 计算多种水体指数 ndwi = (green - nir) / (green + nir + 1e-10) # 归一化水体指数 mndwi = (green - swir1) / (green + swir1 + 1e-10) # 改进的归一化水体指数 # 组合指数提取水体 water = ((ndwi > 0) & (mndwi > 0) & (nir < 0.1)).astype(np.uint8) # 保存结果 profile = src.profile profile.update(dtype=rasterio.uint8, count=1) with rasterio.open(output_path, 'w', **profile) as dst: dst.write(water, 1)

5.3 分类流程自动化

对于需要定期处理的任务，可以将整个流程封装成自动化脚本：

import argparse def main(): parser = argparse.ArgumentParser(description='遥感地物分类自动化脚本') parser.add_argument('input', help='输入影像路径') parser.add_argument('output', help='输出分类结果路径') parser.add_argument('--method', choices=['threshold', 'randomforest'], default='threshold', help='分类方法') args = parser.parse_args() if args.method == 'threshold': simple_threshold_classification(args.input, args.output) else: # 这里可以加载预训练的模型或训练数据 training_data = load_training_data() machine_learning_classification(args.input, training_data, args.output) if __name__ == '__main__': main()

这种自动化脚本可以集成到生产环境中，定期处理新获取的遥感数据。