简介面向遥感地质与矿产勘探从业者及高校科研人员这套岩性智能识别系统以极端随机树为核心分类器并引入布谷鸟算法、粒子群优化进行超参数寻优覆盖多光谱遥感图像处理、高光谱数据特征提取到岩性分类识别与地质填图自动化的完整流程。包内共12个文件主体为8个Python脚本对应数据预处理、格式转换、多表合并与多种训练策略等环节另有说明文件、Markdown文档及附赠资源文档辅助阅读预训练模型结果以pickle格式保存总大小约92KB便于直接复现和二次开发。项目模块划分清晰从数据读取到模型评估均有独立脚本承接适合有一定机器学习基础、希望将元启发式算法用于地质遥感分类的开发者改造使用。已有47人学习在同类岩性智能识别案例中属于轻量而完整的参考样例。1. 遥感地质学里的岩性智能识别为什么是“树模型群体优化”而不是深度网络一张30米分辨率的多光谱影像摆在GIS里岩性界线在假彩色合成图上其实看得见但真要把每条界线都画出来地质队员在野外跑一个月也不一定画得完。遥感地质学这几年把岩性识别这件事交给机器学习很多人一上来直接堆深度网络结果训练样本不够、标签噪声又大效果反而不如极端随机树模型加认真调参来得稳。这套系统打包的路线就是后者用极端随机树模型做分类器用布谷鸟算法和粒子群优化算法去搜超参数把多光谱遥感图像处理和高光谱数据特征提取在进模型之前先管干净最终输出适合矿产资源勘探与地质填图的岩性分类图。适合谁用地质调查院、矿产勘查项目的技术员以及想把自己算法落到地学场景的机器学习工程师。它解决的不是算法炫技而是野外填图的效率问题——把目视解译交给模型人只负责验证。2. 多光谱与高光谱数据管线预处理和特征提取决定分类上限所有岩性分类项目的上限都不在模型而在进模型之前那批特征。多光谱和高光谱的差别不只是波段数量多光谱每个波段几十纳米宽高光谱可以细到几纳米后者能看见吸收谷但也带来更强的噪声和波段冗余。下面这套处理管线我几乎每个项目都复用改的只是波段号和传感器参数。2.1 数据源选择与预处理标准流程做岩性识别数据源选择第一。免费且覆盖好的星载多光谱数据一般用Landsat 8/9 OLI和Sentinel-2前者30米分辨率、11个波段后者10米分辨率、12个波段两者都有短波红外对识别粘土矿物和碳酸盐很关键。高光谱数据通常来自机载或星载传感器波段数从几十到几百不等光谱分辨率高但覆盖范围小、成本高一般只在重点矿区或图幅内做局部精细分类。数据源类型典型波段数空间分辨率适用场景多光谱Landsat OLI1130 m区域地质填图、构造解译多光谱Sentinel-21210 m中小尺度岩性单元识别高光谱机载/星载100 以上数米至数千米矿区蚀变矿物填图、精细岩性分类预处理顺序是固定的辐射定标 → 大气校正 → 几何配准 → 云和阴影掩膜 → 重采样到统一像元尺寸。辐射定标把DN值转成传感器处表观反射率大气校正把表观反射率进一步转成地表反射率。大气校正常见做法是FLAASH、6S这类物理模型或者DOS暗像元减法这种快速做法。项目初期验证算法时用DOS节省时间正式出图用物理模型更稳。几何配准这一步容易被忽略不同来源的影像之间如果存在半个像元的偏移训练时模型会把“错位”当成特征学进去到野外实测阶段就会翻车。云和阴影用QA波段或蓝波段阈值做掩膜掩掉后不要填0填NaN后续建模时单独处理。2.2 光谱指数、PCA/MNF与吸收特征三类特征怎么组合特征提取的目标是让岩性的光谱差异在数值上更明显。多光谱方面短波红外波段对粘土矿物和碳酸盐矿物最敏感这类矿物在2.1到2.4微米附近有特征吸收所以用SWIR1/SWIR2比值做粘土矿物指数是常规操作。铁氧化物在红波段和蓝波段之间有吸收差用红波段除以蓝波段能突出铁染信息。植被覆盖区还得加一个NDVI作为辅助特征不是只用来做掩膜而是让模型知道“这个像元有植被干扰光谱可能被混合”。import numpy as np def build_multispectral_features(band_stack, blue_idx, red_idx, nir_idx, swir1_idx, swir2_idx): blue band_stack[..., blue_idx].astype(np.float32) red band_stack[..., red_idx].astype(np.float32) nir band_stack[..., nir_idx].astype(np.float32) swir1 band_stack[..., swir1_idx].astype(np.float32) swir2 band_stack[..., swir2_idx].astype(np.float32) ndvi (nir - red) / (nir red 1e-8) clay_ratio swir1 / (swir2 1e-8) iron_ratio red / (blue 1e-8) return np.stack([blue, red, nir, swir1, swir2, ndvi, clay_ratio, iron_ratio], axis-1)这里每个比值都加了一个极小常数防止除零。不同传感器的波段索引不一样换数据源时最容易出错的就是把Landsat的波段索引套到Sentinel-2上后面避坑章节会专门讲。这些指数计算完后和原始反射率波段一起叠进特征矩阵。高光谱的特征提取思路不同。波段太密集原始反射率直接作为特征会让维度上千训练样本往往只有几千个像元模型很容易把噪声波段的抖动当成规律。MNF最小噪声分离是高光谱降维的标准做法它和PCA的区别在于MNF先估计噪声协方差矩阵把数据白化后再做主成分变换降维后的分量按信噪比排序前20到30个分量基本包含了绝大部分岩性信息。PCA在高光谱上容易被噪声波段主导所以我不太推荐。连续统去除包络线去除是另一个常用方法把光谱曲线归一化到0到1的范围突出吸收峰的相对深度这东西在不同影像之间可比性比原始反射率强适合做矿物识别。2.3 特征矩阵长什么样从像元光谱到训练样本监督分类最终需要的是“像元-标签”对。把所有特征沿通道方向叠起来一个像元就是一个特征向量。标签来源有三个野外路线调查点、已有地质图矢量化、高分辨率影像目视解译。野外点质量最高但数量少地质图矢量化效率高但会遇到成因单元与光谱单元不对应的问题目视解译居中。我的习惯是三者混用但只把野外点作为最终验证集。训练样本不能只取小块连续区域。相邻像元之间空间自相关很强如果训练区和验证区挨在一起模型实际上记住了位置而不是岩性光谱交叉验证分数会虚高一换到新图幅就现原形。正确做法是把训练区和验证区按空间分块比如一个图幅内切出两块不相邻的区域一块训练、一块验证。这一点在第6章还会展开。3. 极端随机树与两种优化算法为什么是它们而不是随机森林和网格搜索模型选型这里标题给了三个关键词极端随机树模型、布谷鸟算法、粒子群优化算法。先理解三者各自的定位再谈组合。3.1 极端随机树的分裂逻辑与随机森林的区别随机森林在节点分裂时会对每个候选特征的每个可能切分点做评估找出信息增益最大的那个切分位置。极端随机树Extra Trees不一样它对每个特征随机生成一个切分阈值再在这些随机阈值里挑增益最大的那个而且每棵树训练时使用全部样本不做bootstrap抽样。这两个差异带来两个工程结果训练速度更快因为不需要对特征值排序寻找最佳切分点模型方差更小因为随机性更强多棵树平均后整体方差被进一步压低。高光谱特征维度上千时极端随机树这种“不考虑最优切分点”的做法反而抑制过拟合。深层网络的参数空间太大在几千个训练像元上很容易把噪声学进去极端随机树在样本量小得多的条件下就能稳定上分而且对标签噪声的容忍度更高。遥感标签几乎不可能干净野外打点错位、地质图边界概括都会引入噪声树模型在这个场景下比深度网络更实用。3.2 布谷鸟算法莱维飞行与寄生繁殖的调参直觉布谷鸟搜索Cuckoo Search的核心是三条理想化规则每只布谷鸟一次产一个蛋随机放进一个宿主巢质量好的蛋解会保留到下一代宿主巢数量固定外来蛋被宿主发现的概率是pa发现后宿主丢弃或重建巢。新解通过莱维飞行产生莱维飞行的步长服从重尾分布意味着大多数时候是小步探索偶尔来一次大步跳跃。这个特性对超参数搜索非常值钱。小步保证在好参数附近精细搜索大步帮助模型跳出局部最优。超参数空间是连续的n_estimators、max_depth、min_samples_leaf这些参数之间存在低相关性常规网格搜索要试几百组参数布谷鸟搜索通常几十次评估就能收敛到可用的参数区间。pa默认取0.25左右巢数n一般取10到25具体搜索次数看模型单次评估成本。3.3 粒子群优化速度-位置更新与CS-PSO混合策略粒子群优化PSO出现得更早每个粒子代表一组候选参数有位置和速度两个属性。速度更新三部分组成惯性项w乘当前速度、个体认知项飞向自己历史最优位置、社会认知项飞向群体最优位置。惯性权重w取0.6到0.9学习因子c1、c2通常取1.5左右。那为什么标题同时给布谷鸟和粒子群因为单独用任何一个都会翻车。PSO收敛快但群体容易被拉到同一个局部最优附近这就是常说的“早熟”布谷鸟随机性强后期收敛慢在评估次数有限时可能浪费大量计算。工程上常见的混合策略是先用PSO快速圈定几个有希望的区域再用布谷鸟算法的莱维飞行在周边做精细搜索或者每迭代几轮把PSO的全局最优位置替换成一个布谷鸟随机跳跃产生的新位置。前者用PSO的开发能力后者用布谷鸟的探索能力互补性很强。4. 从特征到岩性图训练、CS-PSO调参与全图推理的完整流程这一章直接给可复现的流程。假设已经完成了预处理和特征提取得到了一个形状为n_samples, n_features的特征矩阵X和对应的标签向量y。如果是从头开始做先跑通基准模型再做CS-PSO调参最后全图推理。4.1 搭建训练管线数据划分与模型评估指标训练管线第一步是划分数据。用分层抽样的方式划分训练集和测试集保证每个岩性类别在两边都有足够样本。评估指标用宏平均F1而不是准确率因为岩性类别天然不平衡准确率会被大的岩性单元主导宏F1对每个类别的表现一视同仁。import numpy as np from sklearn.ensemble import ExtraTreesClassifier from sklearn.model_selection import StratifiedKFold, cross_val_score X np.load(feature_matrix.npy) # (n_samples, n_features) y np.load(labels.npy) # 整数编码从0开始 base_model ExtraTreesClassifier( n_estimators100, max_depth20, min_samples_leaf2, n_jobs-1, random_state42 ) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(base_model, X, y, cvcv, scoringf1_macro) print(fbaseline macro-F1: {scores.mean():.4f} /- {scores.std():.4f})这里固定random_state保证结果可复现n_jobs-1让模型并行训练。min_samples_leaf设置成2而不是默认的1是因为遥感标签存在空间噪声叶子节点太细容易把噪声像元画成独立图斑。这个基准分数很重要后面调出来的参数如果比基准还差说明搜索过程出了问题。4.2 用CS-PSO搜索极端随机树最优参数参数搜索范围先定好。n_estimators在50到300之间取整数太大训练慢且边际收益递减max_depth在5到50之间树太深会过拟合min_samples_leaf在1到10之间遥感数据最低给到2比较稳。def evaluate_params(X, y, params): model ExtraTreesClassifier( n_estimatorsint(params[0]), max_depthint(params[1]), min_samples_leafint(params[2]), n_jobs-1, random_state42 ) return cross_val_score(model, X, y, cv5, scoringf1_macro).mean() bounds np.array([[50, 300], [5, 50], [1, 10]]) def pso_search(X, y, iters20, n_particles15): positions np.random.uniform(bounds[:, 0], bounds[:, 1], size(n_particles, bounds.shape[0])) velocities np.random.uniform(-3, 3, sizepositions.shape) pbest positions.copy() pbest_score np.array([evaluate_params(X, y, p) for p in positions]) gbest pbest[np.argmax(pbest_score)] gbest_score pbest_score.max() for _ in range(iters): w, c1, c2 0.7, 1.5, 1.5 r1 np.random.random(positions.shape) r2 np.random.random(positions.shape) velocities (w * velocities c1 * r1 * (pbest - positions) c2 * r2 * (gbest - positions)) velocities np.clip(velocities, -5, 5) positions np.clip(positions velocities, bounds[:, 0], bounds[:, 1]) scores np.array([evaluate_params(X, y, p) for p in positions]) improve scores pbest_score pbest[improve] positions[improve] pbest_score[improve] scores[improve] if pbest_score.max() gbest_score: gbest pbest[np.argmax(pbest_score)] gbest_score pbest_score.max() return gbest, gbest_score best_params, best_score pso_search(X, y, iters15, n_particles12) print(fPSO best: {best_params}, macro-F1: {best_score:.4f})PSO里每个粒子就是一组参数位置速度控制参数移动的幅度。速度clip到[-5,5]是关键不然粒子可能一步飞出边界后面几次迭代都在边界上反复横跳。位置clip保证搜索始终在预设参数范围内。布谷鸟搜索接着做精调。把PSO找到的最优位置作为布谷鸟的一个初始巢再做若干轮莱维飞行。def levy_flight(pos, bounds, scale0.3): step np.random.randn(*pos.shape) * scale return np.clip(pos step, bounds[:, 0], bounds[:, 1]) nests np.random.uniform(bounds[:, 0], bounds[:, 1], size(10, bounds.shape[0])) nests[0] best_params nests_score np.array([evaluate_params(X, y, n) for n in nests]) pa 0.25 for _ in range(15): new_nests np.array([levy_flight(n, bounds) for n in nests]) new_scores np.array([evaluate_params(X, y, n) for n in new_nests]) better new_scores nests_score nests[better] new_nests[better] nests_score[better] new_scores[better] # 宿主发现外来蛋随机丢弃部分较差巢 worst_idx np.argsort(nests_score)[:int(pa * len(nests))] for idx in worst_idx: nests[idx] np.random.uniform(bounds[:, 0], bounds[:, 1]) nests_score[idx] evaluate_params(X, y, nests[idx]) best_idx np.argmax(nests_score) print(fCS fine-tune best: {nests[best_idx]}, macro-F1: {nests_score[best_idx]:.4f})布谷鸟的莱维飞行尺度参数scale控制步长0.3意味着大部分新位置离旧位置不会太远偶尔会有跳跃。pa取0.25意味着每轮有约四分之一较差的巢被随机重置这是为了保持种群多样性。提示调参阶段先用2000到5000个像元做子集。全图几十万像元时一次交叉验证要几分钟跑完整轮CS-PSO会等到怀疑人生。子集调出的参数往往在全集上依然可用幅度不大。4.3 全图推理从训练结果到岩性分类图模型调好后就到了全图推理。这一步最常见的问题是整景影像一次性喂进模型导致内存溢出。处理方法是分块预测import rasterio def predict_full_scene(model, features, profile, block_size256): rows, cols features.shape[0], features.shape[1] out np.zeros((rows, cols), dtypenp.uint8) for i in range(0, rows, block_size): for j in range(0, cols, block_size): block features[i:i block_size, j:j block_size] h, w block.shape[0], block.shape[1] flat block.reshape(-1, block.shape[2]) valid np.any(np.isfinite(flat), axis1) pred np.zeros(h * w, dtypenp.uint8) pred[valid] model.predict(flat[valid]) out[i:i h, j:j w] pred.reshape(h, w) with rasterio.open(lithology_map.tif, w, **profile) as dst: dst.write(out, 1) return outblock_size取256每次处理65536个像元配合n_features通常是几十到一两百内存占用完全可以接受。valid掩膜非常重要预处理阶段掩膜掉的云、阴影、水体像元在特征矩阵里是NaN不筛掉直接predict会报错或者产生垃圾类别。5. 岩性识别落地避坑五个常见问题与排查记录这套流程走下来踩坑记录能写满两页纸。把最常见的五个记在这里都是实际项目中死磕过的。5.1 特征矩阵拼接报错波段顺序与像元网格对不上现象代码运行时报维度不一致或者分类结果出图后岩性界线呈现明显的条纹状错位。原因Landsat和Sentinel-2的波段顺序不同直接按习惯的波段索引取值取到的根本不是同一个波长位置另外不同来源影像重采样到不同像元大小时坐标网格没有对齐。解决写一个波段映射表明确记录每个波段在数组里的索引和中心波长所有数据源统一重采样到同一个像元尺寸和投影拼接前打印每个特征矩阵的shape交叉检查。5.2 野外采样点与影像像元错位GPS误差引发的标签噪声现象训练集交叉验证F1高达0.9野外抽查准确率只有一半左右。原因手持GPS精度通常在3到10米在30米分辨率的Landsat影像上还算能接受在米级分辨率的高光谱影像上一个打点直接落到相邻像元上就是错误标签。解决打点位置不要直接转成单像元标签用半径5到10米的缓冲区做多数投票或者野外采样时拍照记录像元中心位置。这是我的血泪经验早期项目吃过很大亏。5.3 高光谱维度灾难上千波段让极端随机树过拟合现象只用原始高光谱反射率不做降维训练集分数接近1.0全图预测结果却是椒盐噪声一片。原因上千个波段里噪声波段占比高训练样本量相对少极端随机树即便随机性再强也会把噪声波段的抖动当成岩性规律。解决先做MNF或PCA降到20到30维再做分类验证时用按空间分块而不是随机划分的数据集不然虚高分数会掩盖过拟合。5.4 调参搜索跑飞PSO速度失控与布谷鸟巢数爆炸现象PSO迭代几轮后所有粒子挤到边界上布谷鸟搜索换了十几次巢都回不到之前PSO找到的分数。原因粒子速度没有限制位置越界后又被clip导致粒子在边界附近反复横跳布谷鸟的pa开太大好巢频繁被随机重置搜索退化成随机采样。解决速度每次更新后clip到参数范围幅度的10%到20%pa用0.2到0.3之间每次迭代记录最优分数变化曲线如果分数停滞超过5轮就检查是不是搜索范围太宽或者速度设置过大。5.5 类别不平衡被忽略小岩性单元总是被吞掉现象岩脉、蚀变带这类面积小的岩性在分类图上几乎消失后处理阶段又被滤波抹掉。原因样本数量天然悬殊模型训练时被大类主导只要把背景类预测对了准确率就很高宏F1偏低但没人看。解决训练时开启class_weightbalanced评估指标只用宏F1和Kappa对少数类做SMOTE过采样时要小心别在同一个地理位置附近重复采样否则验证集和训练集共享相似样本精度虚高。6. 把分类结果变成地质图后处理、野外验证与特征重要性核查模型输出的原始分类结果是椒盐噪声严重的栅格直接拿去填图会被地质队员拒收。我一般做三步后处理用3×3或5×5的众数滤波去掉孤立像元用连通域分析删除面积小于最小可填图单元例如0.25平方公里的碎斑按地质图规范把相近类别归并成一个谱系类比如把不同粒度的花岗岩类归并为一个大类。如果项目需要矢量成果再转矢量并简化边界。验证不能只看训练分数。我把野外独立检查点分成两份一份在调参时用一份留到模型全部定稿才拿出来算混淆矩阵和Kappa。Kappa不低于0.75我才敢把图交出去。另外一定要打印极端随机树的feature_importances_对照地质常识检查如果粘土矿物指数排在首位合理如果NDVI排首位说明植被干扰没有清理干净模型学的是植被分布而不是岩性。验证方式具体操作可接受阈值独立验证集空间分块留出的野外检查点Kappa 不低于 0.75特征重要性核查对比排序结果与已知地质规律诊断性特征排前野外抽查随机抽点实地核对岩性错误点集中处补样本我早期交过一张“训练集极漂亮、野外一查就翻车”的岩性图从那以后验证区和训练区永远按空间切分连调参阶段都只用训练区的子集。后来每个项目我都留一块没参与任何训练的验证区确认无误再出图。岩性识别这个方向模型能跑通只是开始后处理、验证、可解释性这三关过了才算真正能用。希望帮到你。本文还有配套的精品资源点击获取