ArcGIS与随机森林结合生成不确定性分布图的入门指南

1次阅读

没有评论

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

在传统的地理空间数据分析中，克里金插值等确定性方法占据主导地位。但这些方法存在明显的局限性：

对小样本数据敏感，当训练数据不足时容易产生过拟合
难以量化预测结果的不确定性范围
对非线性关系的表达能力有限

特别是在环境监测、地质灾害评估等场景中，决策者往往不仅需要预测值，更需要了解预测的可信程度。这就需要我们寻找能够同时输出预测值和不确定性评估的方法。

随机森林算法因其独特的优势成为空间不确定性建模的理想选择：

与克里金相比：能自动捕捉非线性特征，不依赖严格的平稳性假设
与贝叶斯最大熵相比：计算效率更高，适合处理高维特征
自带 OOB 误差估计，天然支持不确定性量化

最关键的是，随机森林的树结构差异本身就反映了模型对数据不确定性的响应，这为不确定性分布图的生成提供了理论基础。

使用 ArcPy 进行数据准备工作：

import arcpy
from arcpy.sa import *

# 设置工作空间
arcpy.env.workspace = "input.gdb"
arcpy.env.overwriteOutput = True

# 坐标系转换
arcpy.Project_management("raw_data.shp", "projected_data.shp", 
                         "PROJCS['WGS_1984_UTM_Zone_50N']")

# 异常值处理（以 NDVI 为例）def clean_ndvi(ndvi_layer):
    # 将无效值设为 NoData
    out_con = Con((ndvi_layer >= -1) & (ndvi_layer <= 1), 
                 ndvi_layer)
    # 使用焦点统计填充异常值
    out_focal = FocalStatistics(out_con, NbrRectangle(3,3), "MEAN")
    return out_focal

利用 scikit-learn 实现核心建模流程：

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score
import numpy as np

# 构造空间特征矩阵
X = np.column_stack([arcpy.RasterToNumPyArray("elevation.tif").flatten(),
    arcpy.RasterToNumPyArray("ndvi.tif").flatten(),
    # 添加其他特征...
])

# 移除无效值
valid_mask = ~np.isnan(X).any(axis=1)
X_clean = X[valid_mask]
y_clean = y[valid_mask]

# 初始化模型
rf = RandomForestRegressor(
    n_estimators=500,
    oob_score=True,
    n_jobs=-1  # 启用并行
)

# 交叉验证评估
cv_scores = cross_val_score(rf, X_clean, y_clean, cv=5)
print(f"CV R2 平均得分: {np.mean(cv_scores):.3f}")

# 训练模型
rf.fit(X_clean, y_clean)

采用 OOB 方法计算预测标准差：

# 获取各样本的 OOB 预测
oob_pred = np.stack([tree.predict(X_clean) for tree in rf.estimators_
    if hasattr(tree, 'tree_')
], axis=0)

# 计算标准差作为不确定性指标
oob_std = np.std(oob_pred, axis=0)

将不确定性结果渲染为专题地图：

# 将 numpy 数组转回栅格
def array_to_raster(array, template_raster):
    arcpy.env.outputCoordinateSystem = template_raster
    result = arcpy.NumPyArrayToRaster(
        array,
        arcpy.Point(template_raster.extent.XMin, template_raster.extent.YMin),
        template_raster.meanCellWidth,
        template_raster.meanCellHeight
    )
    return result

# 生成不确定性分布图
uncertainty_raster = array_to_raster(oob_std.reshape(arcpy.RasterToNumPyArray("ndvi.tif").shape),
    "ndvi.tif"
)

# 设置分类渲染
classified = Reclassify(
    uncertainty_raster,
    "VALUE",
    RemapRange([[0,0.1,"1"], [0.1,0.2,"2"], [0.2,1,"3"]])
)

# 导出结果
classified.save("uncertainty_map.tif")

处理大型栅格时采用分块策略：

# 设置处理区块大小
arcpy.env.compression = "LZ77"
arcpy.env.tileSize = "256 256"

# 使用迭代器分块处理
for x, y in tile_coordinates:
    extent = f"{x} {y} {x+256} {y+256}"
    arcpy.env.extent = extent
    # 执行分块处理...

在 ArcGIS Pro 中配置并行处理：

打开 Geoprocessing Options
设置 Parallel Processing Factor 为 70-80%
在 Python 脚本中添加：arcpy.env.parallelProcessingFactor = "75%"

使用 SHAP 分析特征重要性：

import shap

# 计算 SHAP 值
explainer = shap.TreeExplainer(rf)
shap_values = explainer.shap_values(X_clean)

# 可视化
shap.summary_plot(shap_values, X_clean, feature_names=["高程","NDVI"])

本方法可进一步扩展至：