简介：面向遥感地质与矿产勘探从业者及高校科研人员，这套岩性智能识别系统以极端随机树为核心分类器，并引入布谷鸟算法、粒子群优化进行超参数寻优，覆盖多光谱遥感图像处理、高光谱数据特征提取到岩性分类识别与地质填图自动化的完整流程。包内共12个文件，主体为8个Python脚本，对应数据预处理、格式转换、多表合并与多种训练策略等环节，另有说明文件、Markdown文档及附赠资源文档辅助阅读，预训练模型结果以pickle格式保存，总大小约92KB，便于直接复现和二次开发。项目模块划分清晰，从数据读取到模型评估均有独立脚本承接，适合有一定机器学习基础、希望将元启发式算法用于地质遥感分类的开发者改造使用。已有47人学习，在同类岩性智能识别案例中属于轻量而完整的参考样例。
1. 遥感地质学里的岩性智能识别：为什么是“树模型+群体优化”而不是深度网络
一张30米分辨率的多光谱影像摆在GIS里，岩性界线在假彩色合成图上其实看得见，但真要把每条界线都画出来，地质队员在野外跑一个月也不一定画得完。遥感地质学这几年把岩性识别这件事交给机器学习，很多人一上来直接堆深度网络，结果训练样本不够、标签噪声又大，效果反而不如极端随机树模型加认真调参来得稳。这套系统打包的路线就是后者：用极端随机树模型做分类器，用布谷鸟算法和粒子群优化算法去搜超参数，把多光谱遥感图像处理和高光谱数据特征提取在进模型之前先管干净，最终输出适合矿产资源勘探与地质填图的岩性分类图。适合谁用？地质调查院、矿产勘查项目的技术员，以及想把自己算法落到地学场景的机器学习工程师。它解决的不是算法炫技，而是野外填图的效率问题——把目视解译交给模型，人只负责验证。
2. 多光谱与高光谱数据管线：预处理和特征提取决定分类上限
所有岩性分类项目的上限都不在模型，而在进模型之前那批特征。多光谱和高光谱的差别不只是波段数量：多光谱每个波段几十纳米宽，高光谱可以细到几纳米，后者能看见吸收谷，但也带来更强的噪声和波段冗余。下面这套处理管线我几乎每个项目都复用，改的只是波段号和传感器参数。
2.1 数据源选择与预处理标准流程
做岩性识别，数据源选择第一。免费且覆盖好的星载多光谱数据一般用Landsat 8/9 OLI和Sentinel-2，前者30米分辨率、11个波段，后者10米分辨率、12个波段，两者都有短波红外，对识别粘土矿物和碳酸盐很关键。高光谱数据通常来自机载或星载传感器，波段数从几十到几百不等，光谱分辨率高，但覆盖范围小、成本高，一般只在重点矿区或图幅内做局部精细分类。
| 数据源类型 | 典型波段数 | 空间分辨率 | 适用场景 |
|---|---|---|---|
| 多光谱（Landsat OLI） | 11 | 30 m | 区域地质填图、构造解译 |
| 多光谱（Sentinel-2） | 12 | 10 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_estimators=100, max_depth=20, min_samples_leaf=2, n_jobs=-1, random_state=42 ) cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) scores = cross_val_score(base_model, X, y, cv=cv, scoring='f1_macro') print(f'baseline 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_estimators=int(params[0]), max_depth=int(params[1]), min_samples_leaf=int(params[2]), n_jobs=-1, random_state=42 ) return cross_val_score(model, X, y, cv=5, scoring='f1_macro').mean() bounds = np.array([[50, 300], [5, 50], [1, 10]]) def pso_search(X, y, iters=20, n_particles=15): positions = np.random.uniform(bounds[:, 0], bounds[:, 1], size=(n_particles, bounds.shape[0])) velocities = np.random.uniform(-3, 3, size=positions.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, iters=15, n_particles=12) print(f'PSO best: {best_params}, macro-F1: {best_score:.4f}')PSO里每个粒子就是一组参数位置，速度控制参数移动的幅度。速度clip到[-5,5]是关键，不然粒子可能一步飞出边界，后面几次迭代都在边界上反复横跳。位置clip保证搜索始终在预设参数范围内。
布谷鸟搜索接着做精调。把PSO找到的最优位置作为布谷鸟的一个初始巢，再做若干轮莱维飞行。
def levy_flight(pos, bounds, scale=0.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(f'CS 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_size=256): rows, cols = features.shape[0], features.shape[1] out = np.zeros((rows, cols), dtype=np.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), axis=1) pred = np.zeros(h * w, dtype=np.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_weight='balanced'，评估指标只用宏F1和Kappa；对少数类做SMOTE过采样时要小心，别在同一个地理位置附近重复采样，否则验证集和训练集共享相似样本，精度虚高。
6. 把分类结果变成地质图：后处理、野外验证与特征重要性核查
模型输出的原始分类结果是椒盐噪声严重的栅格，直接拿去填图会被地质队员拒收。我一般做三步后处理：用3×3或5×5的众数滤波去掉孤立像元；用连通域分析删除面积小于最小可填图单元（例如0.25平方公里）的碎斑；按地质图规范把相近类别归并成一个谱系类，比如把不同粒度的花岗岩类归并为一个大类。如果项目需要矢量成果，再转矢量并简化边界。
验证不能只看训练分数。我把野外独立检查点分成两份，一份在调参时用，一份留到模型全部定稿才拿出来算混淆矩阵和Kappa。Kappa不低于0.75我才敢把图交出去。另外一定要打印极端随机树的feature_importances_，对照地质常识检查：如果粘土矿物指数排在首位，合理；如果NDVI排首位，说明植被干扰没有清理干净，模型学的是植被分布而不是岩性。
| 验证方式 | 具体操作 | 可接受阈值 |
|---|---|---|
| 独立验证集 | 空间分块留出的野外检查点 | Kappa 不低于 0.75 |
| 特征重要性核查 | 对比排序结果与已知地质规律 | 诊断性特征排前 |
| 野外抽查 | 随机抽点实地核对岩性 | 错误点集中处补样本 |
我早期交过一张“训练集极漂亮、野外一查就翻车”的岩性图，从那以后验证区和训练区永远按空间切分，连调参阶段都只用训练区的子集。后来每个项目我都留一块没参与任何训练的验证区，确认无误再出图。岩性识别这个方向，模型能跑通只是开始，后处理、验证、可解释性这三关过了，才算真正能用。希望帮到你。
本文还有配套的精品资源，点击获取
