ArcGIS与随机森林分类：遥感影像智能解译实战指南

1次阅读

共计 2928 个字符，预计需要花费 8 分钟才能阅读完成。

在遥感影像分类领域，最大似然法（Maximum Likelihood Classification）等传统方法长期占据主导地位。然而随着数据复杂度提升，这些方法暴露出明显短板：

光谱混淆问题：城区建筑与裸土、不同植被类型的光谱特征高度相似时，分类错误率显著上升
人工干预多：需要手动设置概率密度函数和阈值，经验依赖性极强
特征利用单一：仅依赖光谱信息，忽略纹理、时序等多维特征

实测数据显示，在 10 类地物场景中，最大似然法的总体精度通常徘徊在 65%-75% 之间，Kappa 系数低于 0.7。

我们对比了三种主流算法在某湿地保护区分类任务中的表现（Landsat 8 影像，30m 分辨率）：

算法	总体精度	Kappa 系数	训练时间(s)
决策树	82.3%	0.79	28
SVM（RBF 核）	85.1%	0.83	142
随机森林	89.7%	0.88	63

随机森林的 ROC 曲线下面积（AUC）达到 0.93，尤其在以下场景优势明显：

水体与阴影区分（用户精度提升 19%）
混交林类型识别（生产者精度提升 27%）
建筑材质分类（F1-score 提升 33%）

# 执行环境：ArcGIS Pro 3.1 + Python 3.9
import arcpy
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
import numpy as np
import gdal

波段组合策略直接影响模型表现，推荐使用指数增强特征：

计算 NDVI（归一化植被指数）：(NIR-Red)/(NIR+Red)
添加 NDWI（归一化水体指数）：(Green-NIR)/(Green+NIR)
纹理特征：通过 GLCM（灰度共生矩阵）提取对比度

def calculate_ndvi(red_band, nir_band):
    return (nir_band - red_band) / (nir_band + red_band + 1e-10)

# 使用 GDAL 读取训练样本
def load_training_data(tiff_path, sample_points):
    dataset = gdal.Open(tiff_path)
    bands = [dataset.GetRasterBand(i+1).ReadAsArray() for i in range(dataset.RasterCount)]

    X = []
    y = []
    for point in sample_points:
        x_coord, y_coord, class_id = point
        # 转换坐标到像素位置
        px = int((x_coord - gt[0]) / gt[1])
        py = int((y_coord - gt[3]) / gt[5])

        if 0 <= px < dataset.RasterXSize and 0 <= py < dataset.RasterYSize:
            pixel_values = [band[py, px] for band in bands]
            X.append(pixel_values)
            y.append(class_id)

    return np.array(X), np.array(y)

from typing import Tuple, List
import numpy.typing as npt

def train_rf_model(X: npt.NDArray[np.float32], 
    y: npt.NDArray[np.int32],
    n_estimators: int = 100,
    max_depth: int = 10
) -> RandomForestClassifier:
    """
    训练随机森林分类器

    参数:
        X: 特征矩阵 (n_samples, n_features)
        y: 标签向量 (n_samples,)

    返回:
        训练好的随机森林模型
    """
    try:
        clf = RandomForestClassifier(
            n_estimators=n_estimators,
            max_depth=max_depth,
            class_weight='balanced',
            n_jobs=-1  # 启用多核并行
        )
        return clf.fit(X, y)
    except Exception as e:
        arcpy.AddError(f"模型训练失败: {str(e)}")
        raise

def block_process(input_raster: str, block_size: int = 1024):
    """处理超出内存限制的大型影像"""
    ds = gdal.Open(input_raster)
    width = ds.RasterXSize
    height = ds.RasterYSize

    for y in range(0, height, block_size):
        for x in range(0, width, block_size):
            # 计算当前块的实际尺寸
            xsize = min(block_size, width - x)
            ysize = min(block_size, height - y)

            # 读取数据块
            block_data = ds.ReadAsArray(x, y, xsize, ysize)

            # 在此处添加处理逻辑
            # ...