简介面向遥感图像岩性识别任务提供基于极端随机树模型的完整Python实现并结合布谷鸟搜索与粒子群优化算法对模型参数进行自动调优适合计算机、人工智能、遥感、自动化等专业的在校生、教师或工程师用于课程设计、毕业设计或项目初期验证。压缩包共10个文件仅63KB主体为8个Python脚本覆盖数据预处理、txt/csv转换、模型训练网格搜索、布谷鸟、粒子群以及结果合并等环节另包含1个训练好的RF模型pickle文件和1份README说明文档结构清晰、上手门槛低。资源当前已有159人学习代码经测试可正常运行可用于复现遥感岩性识别流程也可在现有逻辑上修改参数或扩展识别目标作为高分毕设或课设的参考资料较为合适。下载后按README指引即可快速理解各脚本作用与运行顺序。1. 极端随机树做岩性识别为什么不用随机森林岩性识别是遥感图像分类里典型的“样本少、类别杂、边界糊”任务一张多光谱影像对应几种岩石单元标注样本少则几百像素、多则几万像素类别比例经常差出十倍。用极端随机树做这个事是我在做课程项目和地质解译时反复验证过的一条稳路——它在特征维数不高、标注噪声不小的场景下往往比随机森林训练更快、方差更低也更不容易被少数类带偏。这篇笔记从特征矩阵构建、模型选型、参数边界到布谷鸟和粒子群两种优化算法调参再到验证集设计和常见坑位给你一套可以直接照着跑的 Python 流程。适合正在赶岩性识别课程作业的同学也适合拿群智能算法调参却总得不到稳定结果的人。2. 特征准备把遥感影像变成 ExtraTrees 能吃的特征矩阵2.1 光谱特征和纹理怎么搭岩性识别到底喂什么给模型极端随机树本身不关心你喂的是原始波段还是花哨的衍生特征它只关心矩阵里每一列有没有区分度。岩性识别里最常用的组合是“原始波段反射率 邻域窗口统计量”。不同岩性的光谱曲线在可见光到近红外的斜率、吸收深度上有差异这是第一个信号来源但纯光谱对灰岩和白云岩这类“光谱长得像”的类别往往分不开此时表面粗糙度、风化纹理、沟槽密度会提供第二个信号也就是纹理。纹理特征我不太建议一上来就上 GLCM 或者小波那一套成本高、参数多、容易过拟合。实用做法是取 3×3 和 5×5 窗口的均值、方差配合原始波段一起喂给模型。窗口均值相当于低通滤波抹掉传感器噪声窗口方差直接反映邻域辐射变化程度岩性表面的均匀性差异会被它捕捉。特征维度控制在几十维以内就够ExtraTrees 对高维噪声有韧性但样本不够时特征太多反而会稀释每个分裂节点的信息。2.2 训练样本与遥感图像标注三种真值来源和清洗方法岩性识别的真值来源通常有三种野外调查的 GPS 点位、已有地质图的矢量化、以及直接在影像上人工勾绘的逐像素标注。第三种在课程作业里最常见代价是“遥感图像标注”的质量非常不稳定——岩性界线本身是渐变过渡带不同人勾出来的边界能差出几十个像素过渡带像素混入两个类别是家常便饭。我一般会在标注后做两步清洗。第一步把岩性边界向外收缩 5 到 10 个像元把过渡带样本直接扔掉第二步检查每类样本的光谱直方图如果某一类出现明显双峰说明标注里混进了异常地物需要回到影像上看是不是把阴影、植被覆盖区或者云影子圈进来了。样本量上每类至少要有 500 到 2000 个像素类间比例超过 10 倍就必须靠类别权重兜底否则再好的树也会偏心。2.3 用 Python 把 tif 拆成特征矩阵最小实现代码环境上只需要基础的地理栅格和机器学习库先确认装好了再往下走。pip install rasterio scikit-learn scipy numpyimport numpy as np import rasterio from scipy.ndimage import uniform_filter def build_feature_stack(tif_path, out_path, band_count5): 从单景多光谱影像构建逐像素特征矩阵。 band_count: 假设前5个波段参与建模一般对应可见光近红外。 with rasterio.open(tif_path) as src: data src.read().astype(np.float32) # 原始DN/反射率转float32省内存 bands data.shape[0] feat [data] # 第1组原始波段值 # 两种窗口统计量3x3保留细节纹理5x5反映更大范围的岩性均匀度 for win in [3, 5]: smooth uniform_filter(data, size(1, win, win), modereflect) smooth_sq uniform_filter(data ** 2, size(1, win, win), modereflect) feat.append(smooth) feat.append(smooth_sq - smooth ** 2) # Var[X] E[X^2] - E[X]^2 features np.concatenate(feat, axis0) # 按波段方向堆叠 features features.transpose(1, 2, 0).reshape(-1, features.shape[0]) np.save(out_path, features) return features这段代码的逻辑分三层第一层读原始影像并保持波段顺序第二层用uniform_filter计算窗口均值和方差注意size(1, win, win)里的 1 表示不动波段轴只平滑空间维度第三层把原始波段、窗口均值、窗口方差按通道方向拼起来展平成(像素数, 特征数)的矩阵存盘。uniform_filter计算均值用的是递归方式比逐像素滑窗快一个量级reflect边界模式能避免影像边缘出现黑边伪特征。拿到特征矩阵后下一步是把它和标签对齐并且这里有一个关键选择训练集和验证集不能随机打散。from scipy import ndimage def build_samples(feat_path, label_path, out_npz, val_ratio0.2, seed42): X np.load(feat_path) with rasterio.open(label_path) as src: labels src.read(1) # 标签tif0背景1..n岩性类别 flat labels.flatten() mask flat 0 X_all, y_all X[mask], flat[mask] # 按8连通图斑划分保证同一块岩性整体进入训练或验证集 obj, n_obj ndimage.label((labels 0).astype(np.uint8)) obj_ids list(range(1, n_obj 1)) rng np.random.default_rng(seed) val_obj set(rng.choice(obj_ids, max(1, int(val_ratio * n_obj)), replaceFalse)) val_mask np.isin(obj, list(val_obj)) (flat 0) np.savez(out_npz, X_trX_all[~val_mask], y_try_all[~val_mask], X_vaX_all[val_mask], y_vay_all[val_mask])这段代码把标签影像按连通域聚成一个个图斑再以图斑为单位抽样。这样做的原因是遥感影像存在强空间自相关——同一图斑里的像素光谱几乎一致如果随机打散模型训练时见过目标图斑的相似样本验证集精度会虚高到失真。按图斑划分后验证集才真正代表“模型去预测一片没见过的区域”。注意如果某一类岩性只分布在两三个大图斑里按图斑划分会导致验证集缺类跑之前先打印一遍np.unique(y_tr)确认每类都在。3. 极端随机树建模关键参数和先跑通的第一版基线3.1 ExtraTrees 与随机森林的两点本质差异极端随机树Extremely Randomized Trees常被当成随机森林的变体但两者有两点本质区别。第一随机森林每棵树用 bootstrap 重采样得到近似一样的训练子集而 ExtraTrees 直接用全部训练样本不重抽样。第二随机森林在分裂节点时遍历候选特征的所有取值找最优切分阈值ExtraTrees 对每个候选特征只随机生成少量阈值、挑其中最好的一个就用。换句话说ExtraTrees 把“找最佳切分点”这一步的随机性拉满单棵树的偏差略高但森林整体的方差更低。放在岩性识别场景里这很划算。遥感影像存在饱和像元、坏线和标注噪声随机森林会把一些噪声处的阈值当成固定结构反复利用ExtraTrees 的随机切分等于自带了一层正则化。训练时间上也有明显优势省掉了每棵树的阈值排序所以树数可以放心加大。我跑岩性数据时ExtraTrees 300 棵树的训练时间大概只有随机森林的 60% 到 70%。3.2 n_estimators、min_samples_leaf、max_features岩性场景下的参数范围岩性识别里真正影响结果的超参数就三个树数、叶节点最小样本数、每节点候选特征数。类别不平衡问题则交给class_weight解决不需要扔进优化器里搜。参数作用岩性场景建议给优化器的搜索范围n_estimators树的数量决定方差下降程度300 左右开始出现精度平台200 ~ 600min_samples_leaf叶节点最少样本数控制过拟合3 ~ 5别用 1岩性标签本身有噪声1 ~ 10max_features每次分裂考察的特征比例特征维度小时从 sqrt 起步再放宽0.3 ~ 0.8class_weight类别权重直接设 balanced_subsample不参与搜索min_samples_leaf是三个参数里最容易被低估的一个。岩性标签的边界地带像素光谱是混合的如果允许叶节点落到单像素模型会把这种无法复现的混合像素当成规则记住预测图会出现大量椒盐点。max_features在 ExtraTrees 里同样生效传 float 时 sklearn 会按特征总数乘比例取整所以优化器里可以直接搜连续值。3.3 先跑一个基线ExtraTreesClassifier 的最短可用代码调参数之前先跑一个手设的基线这个基线分数就是你判断优化算法有没有用的标尺。from sklearn.ensemble import ExtraTreesClassifier from sklearn.model_selection import StratifiedKFold, cross_val_score data np.load(samples.npz) X, y data[X_tr], data[y_tr] et ExtraTreesClassifier( n_estimators300, min_samples_leaf3, max_featuressqrt, class_weightbalanced_subsample, n_jobs-1, # 所有CPU核心并行树多也不怕 random_state42 ) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) score cross_val_score(et, X, y, cvcv, scoringf1_macro) print(fbaseline f1_macro: {score.mean():.4f} /- {score.std():.4f})这里用f1_macro而不是 accuracy是因为岩性类别不平衡时accuracy 会被大类别带跑偏哪怕你完全忽略灰岩也能拿高分。f1_macro对每个类别算 F1 再取平均少数类表现直接反映在分数里。class_weightbalanced_subsample和普通balanced的区别在于它会根据每棵树的样本子集动态调整权重和不做 bootstrap 的 ExtraTrees 配合更自然。如果这个基线的 f1_macro 低于 0.6先别急着调参回头检查特征对齐和标签质量调参治不了数据病。4. 布谷鸟和粒子群调参给 ExtraTrees 洗超参数的正确姿势4.1 别急着网格搜索为什么岩性识别调参会翻车很多人拿到调参任务第一反应是GridSearchCV在岩性识别场景里这通常会在两个地方翻车。第一是计算量爆炸三个参数各取 5 档就是 125 组组合每组组合做 5 折交叉验证、每折训练 300 棵树普通笔记本要跑好几个小时。第二是网格搜索只能覆盖离散格点而最优参数很可能落在格点之间的缝隙里参数之间还有交互效应树数少时合适的max_features和树数多时完全不一样固定网格很难描出这种曲面。布谷鸟和粒子群这类群智能算法把超参数当成连续空间里的位置向量让几十个个体在边界内朝更优区域聚集迭代固定次数后取历史最优。它不会漏掉连续区域里的好点而且每次迭代都可以并行评估时间预算更容易控制。调参问题的核心从“选哪些档位”变成了“设哪些边界”后者经验上更好判断。4.2 布谷鸟搜索Levy 飞行和发现概率 pa 是怎么回事布谷鸟搜索模拟的是布谷鸟借巢孵蛋的行为算法里只有三条规则每只布谷鸟一次只产一个蛋一个蛋对应一个候选解宿主鸟巢里最好的蛋会被保留到下一代剩余的蛋以概率pa被宿主发现并抛弃布谷鸟需要重新生成新蛋。真正起作用的是第二步的移动方式——Levy 飞行。它生成的重尾分布步长里既有大量的小步抖动也有小概率的远距离跳跃这让鸟巢既能局部精修又能跳出局部最优。在参数空间里这表现为前期能跨过无价值的参数平原后期在好参数附近细搜。pa控制的是全局探索强度pa越大越多的巢被重建探索越激进pa太小则容易早熟收敛。岩性调参里我一般把pa设在 0.15 到 0.5alpha缩放系数 0.01 起步。4.3 粒子群优化速度更新与惯性权重的设置粒子群的逻辑比布谷鸟更直观每个粒子记住自己的历史最优位置pbest同时知道整个群体的全局最优gbest下一时刻的速度由三部分合成——惯性项延续当前运动方向认知项把自己拉回个人历史最优社会项把群体拉向全局最优。三个系数里c1和c2通常取 2.0真正需要调的是惯性权重w。岩性调参里推荐的做法是让w从 0.9 线性递减到 0.4前期大步探索参数空间后期精细收敛。另一个容易被忽略的细节是速度上限vmax不设上限的话粒子可能一步从n_estimators200串到600之外在边界上来回震荡浪费迭代。把vmax限制在参数范围的 20% 左右收敛会稳很多。4.4 两个优化器调 ExtraTrees 的 Python 实现两个算法共用同一个目标函数和同一份参数边界。目标函数做 5 折分层交叉验证返回f1_macro优化器负责最大化这个值。import numpy as np from sklearn.model_selection import StratifiedKFold, cross_val_score def etree_objective(params, X_tr, y_tr, folds5): n_estimators int(round(params[0])) min_samples_leaf max(1, int(round(params[1]))) max_features float(np.clip(params[2], 0.1, 1.0)) model ExtraTreesClassifier( n_estimatorsn_estimators, min_samples_leafmin_samples_leaf, max_featuresmax_features, class_weightbalanced_subsample, n_jobs-1, random_state42 ) # 固定折划分和随机种子保证每次比较的是参数差异而不是抽样噪声 cv StratifiedKFold(n_splitsfolds, shuffleTrue, random_state42) return cross_val_score(model, X_tr, y_tr, cvcv, scoringf1_macro).mean()这里有一个必须注意的点目标函数里的StratifiedKFold一定要固定random_state。如果每次调用都重新打乱优化器看到的目标值波动会大于参数变化带来的差异它会在噪声里乱撞最后搜出来的组合比手设基线还差。固定折之后同样的参数每次评估得到相同分数优化器才有方向可循。布谷鸟主体如下Levy 步长用标准 Mantegna 算法生成。import math def levy_flight(beta1.5): sigma (math.gamma(1 beta) * math.sin(math.pi * beta / 2) / (math.gamma((1 beta) / 2) * beta * 2 ** ((beta - 1) / 2))) ** (1 / beta) u np.random.normal(0, sigma) v np.random.normal(0, 1) return u / (abs(v) ** (1 / beta)) # 重尾步长偶尔远跳、多数短移 def cuckoo_search(objective, bounds, n_nest20, pa0.25, max_iter15, alpha0.01, seed42): rng np.random.default_rng(seed) dim len(bounds) lo np.array([b[0] for b in bounds], dtypefloat) hi np.array([b[1] for b in bounds], dtypefloat) nests rng.uniform(lo, hi, size(n_nest, dim)) fit np.array([objective(p) for p in nests]) best_i int(fit.argmax()) best_p, best_f nests[best_i].copy(), fit[best_i] for _ in range(max_iter): # 第一阶段Levy飞行更新所有鸟巢方向偏向当前最优解 for i in range(n_nest): step alpha * levy_flight(1.5) * (nests[i] - best_p) new_nest np.clip(nests[i] step, lo, hi) new_f objective(new_nest) if new_f fit[i]: nests[i], fit[i] new_nest, new_f # 第二阶段宿主以概率pa抛弃旧巢随机重建保证全局探索 for i in range(n_nest): if rng.random() pa: new_nest rng.uniform(lo, hi) new_f objective(new_nest) if new_f fit[i]: nests[i], fit[i] new_nest, new_f cur_best int(fit.argmax()) if fit[cur_best] best_f: best_p, best_f nests[cur_best].copy(), fit[cur_best] return best_p, best_falpha是步长缩放系数。参数边界跨度很大时比如n_estimators从 200 到 600alpha0.01会让位移整体偏小收敛慢这时可以放大到 0.05。pa0.25是个中庸值岩性数据类别越不平衡我越倾向调高到 0.35用更多随机重建抵抗局部陷阱。粒子群主体写法略有不同多了一个速度限制。def particle_swarm(objective, bounds, n_particle20, max_iter15, w_start0.9, w_end0.4, c12.0, c22.0, seed42): rng np.random.default_rng(seed) dim len(bounds) lo np.array([b[0] for b in bounds], dtypefloat) hi np.array([b[1] for b in bounds], dtypefloat) pos rng.uniform(lo, hi, size(n_particle, dim)) vel np.zeros((n_particle, dim)) vmax 0.2 * (hi - lo) # 限速防止粒子在边界震荡 pbest_pos pos.copy() pbest_fit np.array([objective(p) for p in pos]) gbest_idx int(pbest_fit.argmax()) gbest_pos, gbest_fit pbest_pos[gbest_idx].copy(), pbest_fit[gbest_idx] for it in range(max_iter): w w_start - (w_start - w_end) * it / (max_iter - 1) # 惯性线性衰减 r1, r2 rng.random(dim), rng.random(dim) vel w * vel c1 * r1 * (pbest_pos - pos) c2 * r2 * (gbest_pos - pos) vel np.clip(vel, -vmax, vmax) pos np.clip(pos vel, lo, hi) fit np.array([objective(p) for p in pos]) better fit pbest_fit pbest_pos[better], pbest_fit[better] pos[better], fit[better] cur_best int(fit.argmax()) if fit[cur_best] gbest_fit: gbest_pos, gbest_fit pos[cur_best].copy(), fit[cur_best] return gbest_pos, gbest_fit两个优化器的调用方式完全一致参数边界按第 3.2 节的表格设置。data np.load(samples.npz) X_tr, y_tr data[X_tr], data[y_tr] # 闭包目标函数把训练矩阵绑定进去 objective_fn lambda p: etree_objective(p, X_tr, y_tr) bounds [(200, 600), (1, 10), (0.3, 0.8)] best_cs, f_cs cuckoo_search(objective_fn, bounds, n_nest20, pa0.25, max_iter15, seed1) best_pso, f_pso particle_swarm(objective_fn, bounds, n_particle20, max_iter15, seed1) print(CS best:, best_cs, f1_macro:, f_cs) print(PSO best:, best_pso, f1_macro:, f_pso)提示目标函数每次调用都要训练 300 棵树并做 5 折交叉验证两个算法各 15 次迭代加起来约 600 次目标函数调用普通笔记本可能要跑数小时。强烈建议先从训练矩阵里按类别各抽 2000 个像素组成小样本矩阵来调参搜到组合后再用全量样本训练最终模型。布谷鸟和粒子群在这个场景下的表现差异也有规律。布谷鸟的 Levy 飞行带有明显的长尾跳跃参数面越崎岖、存在多个局部优区时越占优粒子群前期收敛快但一旦所有粒子聚到同一个局部峰值附近就不容易逃出来。两者放一起跑还有个隐藏好处如果两个算法搜出的最佳参数完全不同但 f1_macro 接近说明目标函数面很平参数选哪套都行挑训练最快的即可。5. 避坑岩性识别里反复出现的五个调参事故现场5.1 类别不平衡少数岩性在预测图上消失现象训练时打印的 f1_macro 有 0.88但预测全图后灰岩分布区几乎全是背景少数类像素被成片吞掉。原因灰岩样本只有板岩的 5%树节点分裂时按基尼不纯度选阈值多数类的区分收益天然更大少数类分裂点很难被选中。解决给 ExtraTreesClassifier 设class_weightbalanced_subsample让每棵树的样本内部按类别权重重新平衡如果仍然不行把少数类样本复制 3 到 5 倍再做训练这不是做数据增强纯粹是给树更多可分裂的少数类样本。5.2 空间自相关精度虚高的真凶现象train_test_split随机划分时交叉验证 f10.94模型一上全图预测就明显不对换成按连通域划分后 f1 掉到 0.71预测图反而和真实地质图对得上。原因同一图斑相邻像素的光谱几乎一样随机划分把模型的“邻居”塞进了验证集模型相当于在背答案不是在做预测。解决用第 2.3 节的ndimage.label按图斑划分训练和验证集写文档时明确说明自己是按空间对象划分而不是随机划分这个细节在地学领域比算法选型更受认可。5.3 特征矩阵内存爆炸现象一景 10000×10000 的 5 波段影像跑build_feature_stack时进程直接被系统杀掉。原因25 维特征矩阵如果按 float64 存全图内存需求超过 25 GB即使只存 float32 也要 12 GB 以上再加窗口统计计算的中间数组轻松翻倍。解决特征矩阵统一用float32并且只在有标签的区域保留特征背景像素不参与训练也不需要保存如果仍然太大把影像裁成若干块分块提取特征和分块预测最后拼接预测结果图不要一口气全图展开。5.4 调参结果反而比手设基线差现象布谷鸟跑了 30 次迭代搜出的参数组合在验证集上的 f1_macro 比第 3.3 节的手设基线还低 0.02。原因有两个一是目标函数里每次交叉验证随机打乱优化器在追噪声二是迭代次数太多个体在某一组验证折上过拟合分数虚高但换折就垮。解决目标函数固定StratifiedKFold的random_state保证同参数同分数打印迭代过程中的best_f曲线连续 5 轮不涨就提前终止min_samples_leaf的搜索下限设成 2别让优化器把单像素叶节点当成宝贝。5.5 分类图椒盐噪声严重现象预测图整体格局正确但岩性界线附近布满黑白椒盐点图斑破碎得没法直接出图。原因逐像素分类只看了每个像素自己的特征没有考虑邻域一致性而岩性在空间上本来就是连续成片的。解决对预测图做一次窗口多数滤波比如用scipy.ndimage的generic_filter配合mode窗口选 5×5 或 7×7更稳妥的做法是先做影像分割对每个分割单元内的预测结果做众数投票把单元内零散像素归一到主要岩性。后处理不改变模型本身但能让最终图件的地质可读性大幅提升。6. 验证与进阶用混淆矩阵和特征重要性给高分说明加分6.1 混淆矩阵比总精度重要得多调参结束后的第一件事不是看 f1_macro而是打印混淆矩阵看哪些类别互相混淆。from sklearn.metrics import confusion_matrix, f1_score final_model ExtraTreesClassifier( n_estimatorsint(best_cs[0]), min_samples_leafint(best_cs[1]), max_featuresbest_cs[2], class_weightbalanced_subsample, n_jobs-1, random_state42 ) final_model.fit(X_tr, y_tr) y_pred final_model.predict(X_va) class_names [灰岩, 板岩, 白云岩, 砂岩] cm confusion_matrix(y_va, y_pred) print(cm) print(per-class f1:, f1_score(y_va, y_pred, averageNone)) print(overall accuracy:, (y_pred y_va).mean())混淆矩阵的意义在于它能指出“哪两类被分混了”。如果灰岩大量被分到白云岩而且反向错分很少说明这两类的光谱特征有重叠但你的样本定义里可能混入了过渡带像素。岩性识别报告里把混淆矩阵的每一行解释清楚比贴一个总精度数字更有说服力这也是课程作业和毕业设计里拉开差距的关键。6.2 特征重要性可视化给文档说明加分的技巧调完模型后feature_importances_是一份现成的可解释材料。特征顺序要和第 2.3 节构建矩阵时的排列顺序严格对应先原始波段再 3×3 均值、方差再 5×5 均值、方差。import pandas as pd feat_names ([fband{b}_raw for b in range(1, 6)] [fband{b}_w3_mean for b in range(1, 6)] [fband{b}_w3_var for b in range(1, 6)] [fband{b}_w5_mean for b in range(1, 6)] [fband{b}_w5_var for b in range(1, 6)]) importance pd.Series(final_model.feature_importances_, indexfeat_names) print(importance.sort_values(ascendingFalse).head(10))如果排名靠前的特征集中在近红外波段的 5×5 方差说明模型在依赖大窗口纹理区分岩性文档里就可以写“风化面粗糙度差异是本次识别的主要信号来源”如果以原始波段为主那说明光谱差异本身足够纹理是辅助。特征重要性和混淆矩阵结合起来正好能回答“为什么选极端随机树”“模型学到了什么”两个问题——把这两页写清楚配上前面的参数搜索迭代曲线整份说明文档的结构就完整了。我现在的习惯是接到任何岩性识别任务第一件事永远是打开影像和标签的直方图而不是急着训练模型。先确认类别比例、样本分布和数据范围再决定要不要上优化算法——调参只是锦上添花样本标注和验证集划分才是决定成败的那一半。这个习惯帮我躲过很多次“模型分数漂亮、出图就翻车”的尴尬希望帮到你。本文还有配套的精品资源点击获取
极端随机树做岩性识别：布谷鸟与粒子群调参实战指南
简介面向遥感图像岩性识别任务提供基于极端随机树模型的完整Python实现并结合布谷鸟搜索与粒子群优化算法对模型参数进行自动调优适合计算机、人工智能、遥感、自动化等专业的在校生、教师或工程师用于课程设计、毕业设计或项目初期验证。压缩包共10个文件仅63KB主体为8个Python脚本覆盖数据预处理、txt/csv转换、模型训练网格搜索、布谷鸟、粒子群以及结果合并等环节另包含1个训练好的RF模型pickle文件和1份README说明文档结构清晰、上手门槛低。资源当前已有159人学习代码经测试可正常运行可用于复现遥感岩性识别流程也可在现有逻辑上修改参数或扩展识别目标作为高分毕设或课设的参考资料较为合适。下载后按README指引即可快速理解各脚本作用与运行顺序。1. 极端随机树做岩性识别为什么不用随机森林岩性识别是遥感图像分类里典型的“样本少、类别杂、边界糊”任务一张多光谱影像对应几种岩石单元标注样本少则几百像素、多则几万像素类别比例经常差出十倍。用极端随机树做这个事是我在做课程项目和地质解译时反复验证过的一条稳路——它在特征维数不高、标注噪声不小的场景下往往比随机森林训练更快、方差更低也更不容易被少数类带偏。这篇笔记从特征矩阵构建、模型选型、参数边界到布谷鸟和粒子群两种优化算法调参再到验证集设计和常见坑位给你一套可以直接照着跑的 Python 流程。适合正在赶岩性识别课程作业的同学也适合拿群智能算法调参却总得不到稳定结果的人。2. 特征准备把遥感影像变成 ExtraTrees 能吃的特征矩阵2.1 光谱特征和纹理怎么搭岩性识别到底喂什么给模型极端随机树本身不关心你喂的是原始波段还是花哨的衍生特征它只关心矩阵里每一列有没有区分度。岩性识别里最常用的组合是“原始波段反射率 邻域窗口统计量”。不同岩性的光谱曲线在可见光到近红外的斜率、吸收深度上有差异这是第一个信号来源但纯光谱对灰岩和白云岩这类“光谱长得像”的类别往往分不开此时表面粗糙度、风化纹理、沟槽密度会提供第二个信号也就是纹理。纹理特征我不太建议一上来就上 GLCM 或者小波那一套成本高、参数多、容易过拟合。实用做法是取 3×3 和 5×5 窗口的均值、方差配合原始波段一起喂给模型。窗口均值相当于低通滤波抹掉传感器噪声窗口方差直接反映邻域辐射变化程度岩性表面的均匀性差异会被它捕捉。特征维度控制在几十维以内就够ExtraTrees 对高维噪声有韧性但样本不够时特征太多反而会稀释每个分裂节点的信息。2.2 训练样本与遥感图像标注三种真值来源和清洗方法岩性识别的真值来源通常有三种野外调查的 GPS 点位、已有地质图的矢量化、以及直接在影像上人工勾绘的逐像素标注。第三种在课程作业里最常见代价是“遥感图像标注”的质量非常不稳定——岩性界线本身是渐变过渡带不同人勾出来的边界能差出几十个像素过渡带像素混入两个类别是家常便饭。我一般会在标注后做两步清洗。第一步把岩性边界向外收缩 5 到 10 个像元把过渡带样本直接扔掉第二步检查每类样本的光谱直方图如果某一类出现明显双峰说明标注里混进了异常地物需要回到影像上看是不是把阴影、植被覆盖区或者云影子圈进来了。样本量上每类至少要有 500 到 2000 个像素类间比例超过 10 倍就必须靠类别权重兜底否则再好的树也会偏心。2.3 用 Python 把 tif 拆成特征矩阵最小实现代码环境上只需要基础的地理栅格和机器学习库先确认装好了再往下走。pip install rasterio scikit-learn scipy numpyimport numpy as np import rasterio from scipy.ndimage import uniform_filter def build_feature_stack(tif_path, out_path, band_count5): 从单景多光谱影像构建逐像素特征矩阵。 band_count: 假设前5个波段参与建模一般对应可见光近红外。 with rasterio.open(tif_path) as src: data src.read().astype(np.float32) # 原始DN/反射率转float32省内存 bands data.shape[0] feat [data] # 第1组原始波段值 # 两种窗口统计量3x3保留细节纹理5x5反映更大范围的岩性均匀度 for win in [3, 5]: smooth uniform_filter(data, size(1, win, win), modereflect) smooth_sq uniform_filter(data ** 2, size(1, win, win), modereflect) feat.append(smooth) feat.append(smooth_sq - smooth ** 2) # Var[X] E[X^2] - E[X]^2 features np.concatenate(feat, axis0) # 按波段方向堆叠 features features.transpose(1, 2, 0).reshape(-1, features.shape[0]) np.save(out_path, features) return features这段代码的逻辑分三层第一层读原始影像并保持波段顺序第二层用uniform_filter计算窗口均值和方差注意size(1, win, win)里的 1 表示不动波段轴只平滑空间维度第三层把原始波段、窗口均值、窗口方差按通道方向拼起来展平成(像素数, 特征数)的矩阵存盘。uniform_filter计算均值用的是递归方式比逐像素滑窗快一个量级reflect边界模式能避免影像边缘出现黑边伪特征。拿到特征矩阵后下一步是把它和标签对齐并且这里有一个关键选择训练集和验证集不能随机打散。from scipy import ndimage def build_samples(feat_path, label_path, out_npz, val_ratio0.2, seed42): X np.load(feat_path) with rasterio.open(label_path) as src: labels src.read(1) # 标签tif0背景1..n岩性类别 flat labels.flatten() mask flat 0 X_all, y_all X[mask], flat[mask] # 按8连通图斑划分保证同一块岩性整体进入训练或验证集 obj, n_obj ndimage.label((labels 0).astype(np.uint8)) obj_ids list(range(1, n_obj 1)) rng np.random.default_rng(seed) val_obj set(rng.choice(obj_ids, max(1, int(val_ratio * n_obj)), replaceFalse)) val_mask np.isin(obj, list(val_obj)) (flat 0) np.savez(out_npz, X_trX_all[~val_mask], y_try_all[~val_mask], X_vaX_all[val_mask], y_vay_all[val_mask])这段代码把标签影像按连通域聚成一个个图斑再以图斑为单位抽样。这样做的原因是遥感影像存在强空间自相关——同一图斑里的像素光谱几乎一致如果随机打散模型训练时见过目标图斑的相似样本验证集精度会虚高到失真。按图斑划分后验证集才真正代表“模型去预测一片没见过的区域”。注意如果某一类岩性只分布在两三个大图斑里按图斑划分会导致验证集缺类跑之前先打印一遍np.unique(y_tr)确认每类都在。3. 极端随机树建模关键参数和先跑通的第一版基线3.1 ExtraTrees 与随机森林的两点本质差异极端随机树Extremely Randomized Trees常被当成随机森林的变体但两者有两点本质区别。第一随机森林每棵树用 bootstrap 重采样得到近似一样的训练子集而 ExtraTrees 直接用全部训练样本不重抽样。第二随机森林在分裂节点时遍历候选特征的所有取值找最优切分阈值ExtraTrees 对每个候选特征只随机生成少量阈值、挑其中最好的一个就用。换句话说ExtraTrees 把“找最佳切分点”这一步的随机性拉满单棵树的偏差略高但森林整体的方差更低。放在岩性识别场景里这很划算。遥感影像存在饱和像元、坏线和标注噪声随机森林会把一些噪声处的阈值当成固定结构反复利用ExtraTrees 的随机切分等于自带了一层正则化。训练时间上也有明显优势省掉了每棵树的阈值排序所以树数可以放心加大。我跑岩性数据时ExtraTrees 300 棵树的训练时间大概只有随机森林的 60% 到 70%。3.2 n_estimators、min_samples_leaf、max_features岩性场景下的参数范围岩性识别里真正影响结果的超参数就三个树数、叶节点最小样本数、每节点候选特征数。类别不平衡问题则交给class_weight解决不需要扔进优化器里搜。参数作用岩性场景建议给优化器的搜索范围n_estimators树的数量决定方差下降程度300 左右开始出现精度平台200 ~ 600min_samples_leaf叶节点最少样本数控制过拟合3 ~ 5别用 1岩性标签本身有噪声1 ~ 10max_features每次分裂考察的特征比例特征维度小时从 sqrt 起步再放宽0.3 ~ 0.8class_weight类别权重直接设 balanced_subsample不参与搜索min_samples_leaf是三个参数里最容易被低估的一个。岩性标签的边界地带像素光谱是混合的如果允许叶节点落到单像素模型会把这种无法复现的混合像素当成规则记住预测图会出现大量椒盐点。max_features在 ExtraTrees 里同样生效传 float 时 sklearn 会按特征总数乘比例取整所以优化器里可以直接搜连续值。3.3 先跑一个基线ExtraTreesClassifier 的最短可用代码调参数之前先跑一个手设的基线这个基线分数就是你判断优化算法有没有用的标尺。from sklearn.ensemble import ExtraTreesClassifier from sklearn.model_selection import StratifiedKFold, cross_val_score data np.load(samples.npz) X, y data[X_tr], data[y_tr] et ExtraTreesClassifier( n_estimators300, min_samples_leaf3, max_featuressqrt, class_weightbalanced_subsample, n_jobs-1, # 所有CPU核心并行树多也不怕 random_state42 ) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) score cross_val_score(et, X, y, cvcv, scoringf1_macro) print(fbaseline f1_macro: {score.mean():.4f} /- {score.std():.4f})这里用f1_macro而不是 accuracy是因为岩性类别不平衡时accuracy 会被大类别带跑偏哪怕你完全忽略灰岩也能拿高分。f1_macro对每个类别算 F1 再取平均少数类表现直接反映在分数里。class_weightbalanced_subsample和普通balanced的区别在于它会根据每棵树的样本子集动态调整权重和不做 bootstrap 的 ExtraTrees 配合更自然。如果这个基线的 f1_macro 低于 0.6先别急着调参回头检查特征对齐和标签质量调参治不了数据病。4. 布谷鸟和粒子群调参给 ExtraTrees 洗超参数的正确姿势4.1 别急着网格搜索为什么岩性识别调参会翻车很多人拿到调参任务第一反应是GridSearchCV在岩性识别场景里这通常会在两个地方翻车。第一是计算量爆炸三个参数各取 5 档就是 125 组组合每组组合做 5 折交叉验证、每折训练 300 棵树普通笔记本要跑好几个小时。第二是网格搜索只能覆盖离散格点而最优参数很可能落在格点之间的缝隙里参数之间还有交互效应树数少时合适的max_features和树数多时完全不一样固定网格很难描出这种曲面。布谷鸟和粒子群这类群智能算法把超参数当成连续空间里的位置向量让几十个个体在边界内朝更优区域聚集迭代固定次数后取历史最优。它不会漏掉连续区域里的好点而且每次迭代都可以并行评估时间预算更容易控制。调参问题的核心从“选哪些档位”变成了“设哪些边界”后者经验上更好判断。4.2 布谷鸟搜索Levy 飞行和发现概率 pa 是怎么回事布谷鸟搜索模拟的是布谷鸟借巢孵蛋的行为算法里只有三条规则每只布谷鸟一次只产一个蛋一个蛋对应一个候选解宿主鸟巢里最好的蛋会被保留到下一代剩余的蛋以概率pa被宿主发现并抛弃布谷鸟需要重新生成新蛋。真正起作用的是第二步的移动方式——Levy 飞行。它生成的重尾分布步长里既有大量的小步抖动也有小概率的远距离跳跃这让鸟巢既能局部精修又能跳出局部最优。在参数空间里这表现为前期能跨过无价值的参数平原后期在好参数附近细搜。pa控制的是全局探索强度pa越大越多的巢被重建探索越激进pa太小则容易早熟收敛。岩性调参里我一般把pa设在 0.15 到 0.5alpha缩放系数 0.01 起步。4.3 粒子群优化速度更新与惯性权重的设置粒子群的逻辑比布谷鸟更直观每个粒子记住自己的历史最优位置pbest同时知道整个群体的全局最优gbest下一时刻的速度由三部分合成——惯性项延续当前运动方向认知项把自己拉回个人历史最优社会项把群体拉向全局最优。三个系数里c1和c2通常取 2.0真正需要调的是惯性权重w。岩性调参里推荐的做法是让w从 0.9 线性递减到 0.4前期大步探索参数空间后期精细收敛。另一个容易被忽略的细节是速度上限vmax不设上限的话粒子可能一步从n_estimators200串到600之外在边界上来回震荡浪费迭代。把vmax限制在参数范围的 20% 左右收敛会稳很多。4.4 两个优化器调 ExtraTrees 的 Python 实现两个算法共用同一个目标函数和同一份参数边界。目标函数做 5 折分层交叉验证返回f1_macro优化器负责最大化这个值。import numpy as np from sklearn.model_selection import StratifiedKFold, cross_val_score def etree_objective(params, X_tr, y_tr, folds5): n_estimators int(round(params[0])) min_samples_leaf max(1, int(round(params[1]))) max_features float(np.clip(params[2], 0.1, 1.0)) model ExtraTreesClassifier( n_estimatorsn_estimators, min_samples_leafmin_samples_leaf, max_featuresmax_features, class_weightbalanced_subsample, n_jobs-1, random_state42 ) # 固定折划分和随机种子保证每次比较的是参数差异而不是抽样噪声 cv StratifiedKFold(n_splitsfolds, shuffleTrue, random_state42) return cross_val_score(model, X_tr, y_tr, cvcv, scoringf1_macro).mean()这里有一个必须注意的点目标函数里的StratifiedKFold一定要固定random_state。如果每次调用都重新打乱优化器看到的目标值波动会大于参数变化带来的差异它会在噪声里乱撞最后搜出来的组合比手设基线还差。固定折之后同样的参数每次评估得到相同分数优化器才有方向可循。布谷鸟主体如下Levy 步长用标准 Mantegna 算法生成。import math def levy_flight(beta1.5): sigma (math.gamma(1 beta) * math.sin(math.pi * beta / 2) / (math.gamma((1 beta) / 2) * beta * 2 ** ((beta - 1) / 2))) ** (1 / beta) u np.random.normal(0, sigma) v np.random.normal(0, 1) return u / (abs(v) ** (1 / beta)) # 重尾步长偶尔远跳、多数短移 def cuckoo_search(objective, bounds, n_nest20, pa0.25, max_iter15, alpha0.01, seed42): rng np.random.default_rng(seed) dim len(bounds) lo np.array([b[0] for b in bounds], dtypefloat) hi np.array([b[1] for b in bounds], dtypefloat) nests rng.uniform(lo, hi, size(n_nest, dim)) fit np.array([objective(p) for p in nests]) best_i int(fit.argmax()) best_p, best_f nests[best_i].copy(), fit[best_i] for _ in range(max_iter): # 第一阶段Levy飞行更新所有鸟巢方向偏向当前最优解 for i in range(n_nest): step alpha * levy_flight(1.5) * (nests[i] - best_p) new_nest np.clip(nests[i] step, lo, hi) new_f objective(new_nest) if new_f fit[i]: nests[i], fit[i] new_nest, new_f # 第二阶段宿主以概率pa抛弃旧巢随机重建保证全局探索 for i in range(n_nest): if rng.random() pa: new_nest rng.uniform(lo, hi) new_f objective(new_nest) if new_f fit[i]: nests[i], fit[i] new_nest, new_f cur_best int(fit.argmax()) if fit[cur_best] best_f: best_p, best_f nests[cur_best].copy(), fit[cur_best] return best_p, best_falpha是步长缩放系数。参数边界跨度很大时比如n_estimators从 200 到 600alpha0.01会让位移整体偏小收敛慢这时可以放大到 0.05。pa0.25是个中庸值岩性数据类别越不平衡我越倾向调高到 0.35用更多随机重建抵抗局部陷阱。粒子群主体写法略有不同多了一个速度限制。def particle_swarm(objective, bounds, n_particle20, max_iter15, w_start0.9, w_end0.4, c12.0, c22.0, seed42): rng np.random.default_rng(seed) dim len(bounds) lo np.array([b[0] for b in bounds], dtypefloat) hi np.array([b[1] for b in bounds], dtypefloat) pos rng.uniform(lo, hi, size(n_particle, dim)) vel np.zeros((n_particle, dim)) vmax 0.2 * (hi - lo) # 限速防止粒子在边界震荡 pbest_pos pos.copy() pbest_fit np.array([objective(p) for p in pos]) gbest_idx int(pbest_fit.argmax()) gbest_pos, gbest_fit pbest_pos[gbest_idx].copy(), pbest_fit[gbest_idx] for it in range(max_iter): w w_start - (w_start - w_end) * it / (max_iter - 1) # 惯性线性衰减 r1, r2 rng.random(dim), rng.random(dim) vel w * vel c1 * r1 * (pbest_pos - pos) c2 * r2 * (gbest_pos - pos) vel np.clip(vel, -vmax, vmax) pos np.clip(pos vel, lo, hi) fit np.array([objective(p) for p in pos]) better fit pbest_fit pbest_pos[better], pbest_fit[better] pos[better], fit[better] cur_best int(fit.argmax()) if fit[cur_best] gbest_fit: gbest_pos, gbest_fit pos[cur_best].copy(), fit[cur_best] return gbest_pos, gbest_fit两个优化器的调用方式完全一致参数边界按第 3.2 节的表格设置。data np.load(samples.npz) X_tr, y_tr data[X_tr], data[y_tr] # 闭包目标函数把训练矩阵绑定进去 objective_fn lambda p: etree_objective(p, X_tr, y_tr) bounds [(200, 600), (1, 10), (0.3, 0.8)] best_cs, f_cs cuckoo_search(objective_fn, bounds, n_nest20, pa0.25, max_iter15, seed1) best_pso, f_pso particle_swarm(objective_fn, bounds, n_particle20, max_iter15, seed1) print(CS best:, best_cs, f1_macro:, f_cs) print(PSO best:, best_pso, f1_macro:, f_pso)提示目标函数每次调用都要训练 300 棵树并做 5 折交叉验证两个算法各 15 次迭代加起来约 600 次目标函数调用普通笔记本可能要跑数小时。强烈建议先从训练矩阵里按类别各抽 2000 个像素组成小样本矩阵来调参搜到组合后再用全量样本训练最终模型。布谷鸟和粒子群在这个场景下的表现差异也有规律。布谷鸟的 Levy 飞行带有明显的长尾跳跃参数面越崎岖、存在多个局部优区时越占优粒子群前期收敛快但一旦所有粒子聚到同一个局部峰值附近就不容易逃出来。两者放一起跑还有个隐藏好处如果两个算法搜出的最佳参数完全不同但 f1_macro 接近说明目标函数面很平参数选哪套都行挑训练最快的即可。5. 避坑岩性识别里反复出现的五个调参事故现场5.1 类别不平衡少数岩性在预测图上消失现象训练时打印的 f1_macro 有 0.88但预测全图后灰岩分布区几乎全是背景少数类像素被成片吞掉。原因灰岩样本只有板岩的 5%树节点分裂时按基尼不纯度选阈值多数类的区分收益天然更大少数类分裂点很难被选中。解决给 ExtraTreesClassifier 设class_weightbalanced_subsample让每棵树的样本内部按类别权重重新平衡如果仍然不行把少数类样本复制 3 到 5 倍再做训练这不是做数据增强纯粹是给树更多可分裂的少数类样本。5.2 空间自相关精度虚高的真凶现象train_test_split随机划分时交叉验证 f10.94模型一上全图预测就明显不对换成按连通域划分后 f1 掉到 0.71预测图反而和真实地质图对得上。原因同一图斑相邻像素的光谱几乎一样随机划分把模型的“邻居”塞进了验证集模型相当于在背答案不是在做预测。解决用第 2.3 节的ndimage.label按图斑划分训练和验证集写文档时明确说明自己是按空间对象划分而不是随机划分这个细节在地学领域比算法选型更受认可。5.3 特征矩阵内存爆炸现象一景 10000×10000 的 5 波段影像跑build_feature_stack时进程直接被系统杀掉。原因25 维特征矩阵如果按 float64 存全图内存需求超过 25 GB即使只存 float32 也要 12 GB 以上再加窗口统计计算的中间数组轻松翻倍。解决特征矩阵统一用float32并且只在有标签的区域保留特征背景像素不参与训练也不需要保存如果仍然太大把影像裁成若干块分块提取特征和分块预测最后拼接预测结果图不要一口气全图展开。5.4 调参结果反而比手设基线差现象布谷鸟跑了 30 次迭代搜出的参数组合在验证集上的 f1_macro 比第 3.3 节的手设基线还低 0.02。原因有两个一是目标函数里每次交叉验证随机打乱优化器在追噪声二是迭代次数太多个体在某一组验证折上过拟合分数虚高但换折就垮。解决目标函数固定StratifiedKFold的random_state保证同参数同分数打印迭代过程中的best_f曲线连续 5 轮不涨就提前终止min_samples_leaf的搜索下限设成 2别让优化器把单像素叶节点当成宝贝。5.5 分类图椒盐噪声严重现象预测图整体格局正确但岩性界线附近布满黑白椒盐点图斑破碎得没法直接出图。原因逐像素分类只看了每个像素自己的特征没有考虑邻域一致性而岩性在空间上本来就是连续成片的。解决对预测图做一次窗口多数滤波比如用scipy.ndimage的generic_filter配合mode窗口选 5×5 或 7×7更稳妥的做法是先做影像分割对每个分割单元内的预测结果做众数投票把单元内零散像素归一到主要岩性。后处理不改变模型本身但能让最终图件的地质可读性大幅提升。6. 验证与进阶用混淆矩阵和特征重要性给高分说明加分6.1 混淆矩阵比总精度重要得多调参结束后的第一件事不是看 f1_macro而是打印混淆矩阵看哪些类别互相混淆。from sklearn.metrics import confusion_matrix, f1_score final_model ExtraTreesClassifier( n_estimatorsint(best_cs[0]), min_samples_leafint(best_cs[1]), max_featuresbest_cs[2], class_weightbalanced_subsample, n_jobs-1, random_state42 ) final_model.fit(X_tr, y_tr) y_pred final_model.predict(X_va) class_names [灰岩, 板岩, 白云岩, 砂岩] cm confusion_matrix(y_va, y_pred) print(cm) print(per-class f1:, f1_score(y_va, y_pred, averageNone)) print(overall accuracy:, (y_pred y_va).mean())混淆矩阵的意义在于它能指出“哪两类被分混了”。如果灰岩大量被分到白云岩而且反向错分很少说明这两类的光谱特征有重叠但你的样本定义里可能混入了过渡带像素。岩性识别报告里把混淆矩阵的每一行解释清楚比贴一个总精度数字更有说服力这也是课程作业和毕业设计里拉开差距的关键。6.2 特征重要性可视化给文档说明加分的技巧调完模型后feature_importances_是一份现成的可解释材料。特征顺序要和第 2.3 节构建矩阵时的排列顺序严格对应先原始波段再 3×3 均值、方差再 5×5 均值、方差。import pandas as pd feat_names ([fband{b}_raw for b in range(1, 6)] [fband{b}_w3_mean for b in range(1, 6)] [fband{b}_w3_var for b in range(1, 6)] [fband{b}_w5_mean for b in range(1, 6)] [fband{b}_w5_var for b in range(1, 6)]) importance pd.Series(final_model.feature_importances_, indexfeat_names) print(importance.sort_values(ascendingFalse).head(10))如果排名靠前的特征集中在近红外波段的 5×5 方差说明模型在依赖大窗口纹理区分岩性文档里就可以写“风化面粗糙度差异是本次识别的主要信号来源”如果以原始波段为主那说明光谱差异本身足够纹理是辅助。特征重要性和混淆矩阵结合起来正好能回答“为什么选极端随机树”“模型学到了什么”两个问题——把这两页写清楚配上前面的参数搜索迭代曲线整份说明文档的结构就完整了。我现在的习惯是接到任何岩性识别任务第一件事永远是打开影像和标签的直方图而不是急着训练模型。先确认类别比例、样本分布和数据范围再决定要不要上优化算法——调参只是锦上添花样本标注和验证集划分才是决定成败的那一半。这个习惯帮我躲过很多次“模型分数漂亮、出图就翻车”的尴尬希望帮到你。本文还有配套的精品资源点击获取