简介:面向遥感地质学与矿产资源勘探的机器学习岩性分类识别平台,专为遥感、地质科研人员及矿产勘探工程师设计。核心融合极端随机树模型、布谷鸟搜索与粒子群优化算法,覆盖多光谱遥感图像处理、高光谱数据特征提取与岩性分类流程,可辅助完成地质填图和找矿预测中的自动化识别任务,尤其适合复杂地质环境下的遥感解译研究。资源包共12个文件,以8个Python脚本为主体,涵盖数据预处理、文本转CSV、布谷鸟搜索、粒子群优化、模型训练与超参数优化等完整实验流程;另含已训练好的pickle模型,可直接加载测试,并配有txt说明文档、md说明文档和docx附赠资料,便于理解算法实现细节并扩展应用到自有数据集。压缩包仅92KB,轻量紧凑、部署门槛低。目前已有47人学习下载,对正在开展岩性智能识别、高光谱数据处理或遥感地质应用研究的读者,可提供一套清晰的代码基线,帮助缩短算法选型与调参周期。
1. 遥感地质学这行最磨人的不是算法,是数据
干了七八年遥感地质解译,我最大的感受是:岩性识别这活儿,真正卡脖子的从来不是分类器选哪个,而是从多光谱、高光谱影像到一张能用的岩性图的中间环节——特征怎么提、样本怎么准备、参数怎么调。极端随机树、布谷鸟算法、粒子群优化这些名词看着唬人,落到矿产资源勘探和地质填图场景里,核心就一句话:用可解释的特征工程把岩石的光谱响应变成分类器吃得下的输入,再用进化算法把模型的超参数调明白。这套基于机器学习的岩性智能识别系统,就是把整条链路打包好——多光谱遥感图像处理、高光谱数据特征提取、布谷鸟算法调参、极端随机树分类、精度验证全都有,适合正在做地质填图、找矿预测或者毕业论文里涉及岩性分类的从业者。我拆完一遍,发现它真正值钱的地方不是某个算法多新,而是流程完整度。
2. 多光谱与高光谱数据预处理:辐射定标、大气校正和特征提取三板斧
2.1 先讲清楚一个前提:DN值不能直接当反射率用
遥感影像买回来大多是DN值(Digital Number),它跟地表反射率之间隔着传感器响应、大气吸收散射、地形辐照度好几层。你要是把DN值直接喂给分类器,分类器学到的可能是不同时相的大气状态差异,而不是岩性差异。这是新手最容易踩的坑,也是我最早翻车的地方——当年拿两期Landsat影像做岩性分类,精度差到没法看,后来排查发现就是一期做了大气校正另一期没做。
所以这套系统的第一个固定动作是:辐射定标把DN转成表观反射率,再做大气校正。常用做法是先查影像自带的元数据文件获取增益(gain)和偏置(bias),逐波段按 L = gain × DN + bias 转辐射亮度,再除以太阳辐照度和太阳天顶角余弦得到表观反射率。代码实现长这样:
import rasterio import numpy as np def dn_to_reflectance(src_path, dst_path, meta_path): with rasterio.open(src_path) as src: meta = src.meta.copy() meta.update(dtype='float32') profile = src.profile profile.update(dtype='float32') with rasterio.open(dst_path, 'w', **profile) as dst: # 读元数据中的增益和偏置,通常为每波段一组 gains = read_gains(meta_path) # 长度为波段数的数组 biases = read_biases(meta_path) esun = read_esun(meta_path) # 太阳表观辐照度,单位 W/(m²·µm) sun_elevation = read_sun_elevation(meta_path) cos_theta = np.cos(np.deg2rad(90 - float(sun_elevation))) for band_idx in range(src.count): dn = src.read(band_idx).astype('float32') # 辐射定标:DN -> 辐亮度 radiance = dn * gains[band_idx] + biases[band_idx] # 表观反射率 reflectance = np.pi * radiance / (esun[band_idx] * cos_theta) dst.write(reflectance, band_idx + 1)参数说明:esun是每波段的太阳光谱辐照度,Landsat 8 OLI 各波段的esun值可以从USGS文档里查到,哨兵二号则大多直接提供反射率产品,不需要自己算。sun_elevation来自影像的元数据,注意有个换算——元数据给的是太阳高度角,计算时要转成天顶角。这里最容易出错的是np.pi要不要乘,取决于辐亮度的单位是W/(m²·sr·µm)还是W/(m²·µm),前者必须乘pi。
做完辐射定标,接着是大气校正。严格做法是用6S或MODTRAN模型,但自己做研究和工程验证时,更常用的是简化黑暗像元法(DOS)。它假设影像里存在反射率近似为零的暗像元(比如深水体、阴影区),通过调整偏移量来逼近真实大气散射。系统里内置的脚本就是基于这个方法,效果对岩性分类来说足够用,因为岩性识别的关键是不同岩石类型之间的相对光谱差异,而不是绝对反射率的精度。
2.2 高光谱特征提取:PCA、MNF与SAM光谱角填图
多光谱影像(Landsat、Sentinel-2)波段少,预处理完直接当输入就行。但高光谱数据动辄上百个波段,波段间相关性极高,直接堆给随机森林或极端随机树,第一是训练极慢,第二是大量冗余波段会稀释真正有判别力的特征。这套系统里针对高光谱数据的处理流程分三步:降维、端元提取、光谱角匹配。
降维方面,主成分分析(PCA)是最朴素的手段,但它只考虑方差最大方向,不保证保留光谱诊断性。更推荐的是最小噪声分离变换(MNF),它本质上是两次PCA,第一次估计噪声协方差并白化,第二次在噪声白化后的空间做PCA。MNF在岩性分类场景里比PCA好用的原因是:岩性之间的光谱差异往往体现在局部吸收特征上,这些特征在方差占比里不大,但在MNF的噪声分离后会凸显出来。
from sklearn.decomposition import PCA import numpy as np def mnf_transform(data, n_components=30): """ MNF简化实现:先估计噪声协方差,再白化,最后PCA data: shape (n_samples, n_bands) """ # 步骤1:估计噪声——用相邻像元差值估算 shifted = np.roll(data, shift=1, axis=0) noise = data - shifted noise_cov = np.cov(noise.T) # 步骤2:白化——噪声协方差的逆平方根 eigvals, eigvecs = np.linalg.eigh(noise_cov) # 加一个小正则防止除零 whitening = eigvecs @ np.diag(1.0 / np.sqrt(eigvals + 1e-10)) @ eigvecs.T whitened = data @ whitening # 步骤3:对白化后的数据再做标准PCA pca = PCA(n_components=n_components) transformed = pca.fit_transform(whitened) return transformed, pca.explained_variance_ratio_这段代码的逻辑要拆开说:np.roll是把整列数据平移一位,用当前像元和上一像元的差值作为噪声的近似估计,这是MNF里最朴素的做法,工程上够用。更好的办法是用空间邻域(上下左右四邻域)的残差来估计噪声,但计算量会大四倍,对超大数据集不太划算。eigh是专门给对称矩阵用的特征分解,协方差矩阵天然对称,用eigh比eig数值更稳定。白化变换的本质是把噪声协方差变成单位阵,让后续PCA不再被高噪声的波段主导。
特征提取的另一路是光谱角填图(SAM)。它的原理是把每个像元的光谱当作高维空间里一个向量,通过计算待测像元与参考光谱之间的夹角来判断归属。SAM对光照强度和地形造成的整体增益变化不敏感,因为它只关心光谱形状,不关心幅值,这一点对野外复杂地形条件下的岩性识别非常实用。
from spectral import read_envi, sam import numpy as np # 读取ENVI格式的高光谱影像 img = read_envi('hyperspectral.img') # 假设已经手动选了几个岩性端元,存成 shape (n_endmembers, n_bands) endmembers = np.load('endmembers.npy') # 逐像元做光谱角匹配 angles = sam(img, endmembers) # 返回 shape (rows, cols, n_endmembers) # 每个像元归类到光谱角最小的端元 classification = np.argmin(angles, axis=2)sam函数会返回每个像元到每个端元的光谱角,单位是弧度。值越小说明越接近该端元。它的一个使用技巧是给角度设个阈值,比如大于0.3弧度的像元不归类,留作未分类,避免把没有代表性的过渡像元强行归进某一类。端元选取既可以用统计方法(如PPI像元纯度指数),也可以直接基于野外实测光谱,实操中后者往往更靠谱,因为仪器实测的纯净光谱比影像里提取的端元更接近岩石本身。
2.3 波段选择:从光谱曲线到诊断性特征
降维是变换,波段选择是筛选。两者目的不同,很多做岩性分类的习惯用PCA一把梭到底,但PCA变换后的波段失去了物理含义,你没法对着输出波段说「这是铁离子的特征吸收位置」。对地质从业人员来说,可解释性非常重要,因为你要在地质填图报告中解释这个分类结果基于什么依据。
波段选择最直接的方法是看典型岩性的光谱曲线。比如含铁矿物在0.85-0.90μm附近有吸收特征,碳酸盐岩在2.30-2.35μm有强吸收,黏土矿物在2.20μm附近有Al-OH吸收。如果你的数据是ASTER或高光谱,直接按已知的诊断性波段区间提取。如果是多光谱,就只能依赖现有波段的组合。
系统里给了个实用工具:波段重要性排序。做法是先用极端随机树跑一遍全波段分类,输出特征重要性得分,然后按得分保留前K个波段。这不算严格的波段选择,但在工程上是极快的启发式方法。后面聊特征筛选章节,会有更系统的方差阈值和递归消除做法。
3. 特征筛选与样本制备:方差阈值、递归特征消除与空间不重叠
3.1 特征工程:为什么多波段反而可能降精度
光谱特征不是越多越好,这是岩性分类里一个反直觉的点。理论上波段越多信息越丰富,但实际上高光谱数据波段之间相关性动辄0.95以上,这些冗余特征带来的问题是双重的:一是树模型的特征选择机制会被大量噪声特征干扰——极端随机树在节点分裂时随机挑特征,冗余波段多了,有判别力的波段被选中的概率反而下降;二是OOB(袋外)误差评估会失真,因为测试样本和训练样本在冗余波段上共享了太多共同信息。
我一般会把特征分成三组:光谱反射率本身、光谱变换特征(一阶导数、二阶导数)、比值指数(如黏土矿物指数、铁氧化物指数)。导数光谱能突出吸收特征的边界位置,比值指数能压制地形阴影效应。系统里内置的特征提取脚本就是按这个思路组织的。
import pandas as pd import numpy as np def derive_spectral_features(df, band_cols): """ df: 包含岩性标签和波段反射率的DataFrame band_cols: 波段列名列表 """ features = df[band_cols].copy() # 一阶导数特征:相邻波段的差分 for i in range(1, len(band_cols)): features[f'deriv_{i}'] = (df[band_cols[i]] - df[band_cols[i-1]]) # 比值指数:经典粘土矿物指数(4/6 for ASTER) features['clay_idx'] = df[band_cols[3]] / (df[band_cols[5]] + 1e-10) features['iron_idx'] = df[band_cols[2]] / (df[band_cols[0]] + 1e-10) # 归一化到0-1,防止量纲差异 for col in features.columns: features[col] = (features[col] - features[col].min()) / \ (features[col].max() - features[col].min() + 1e-10) return features加+1e-10是防止分母出现零。在真实的高光谱数据里,反射率虽然不太可能为0,但做了归一化之后某些像元在某些波段可能非常接近0,加上一个极小值可以避免产生无穷大值。deriv_{i}这里用的是简单差分,没做平滑,如果原始数据噪声大,建议先用Savitzky-Golay滤波器平滑再计算导数,否则导数特征的噪声会被后续模型放大。
3.2 样本准备:按区块划分,别按像元随机抽
这是整套流程里最要命的一步。很多人在准备训练样本时用随机抽样——从整个研究区随机挑一些像元作为训练集,再随机挑一些作为验证集。训练精度很高,一画图就露馅:分类结果图上大片大片的椒盐噪声,岩性界线杂乱无章。原因是遥感影像存在空间自相关性,相邻像元高度相似,随机划分让训练集和验证集包含了几乎相邻的像元,验证精度虚高。
系统里给的做法是空间区块划分。先把研究区切成规则网格,比如500×500米的格子,然后按格子分组,整格划分训练验证集,保证训练和验证不在空间上重叠。更严格的做法是按地质单元划分——同一地质体内部的像元全部分到一组。
from sklearn.model_selection import GroupShuffleSplit import numpy as np def spatial_split(X, y, group_ids, test_size=0.2, random_state=42): """ group_ids: 每个像元所属的空间区块编号 关键点:按区块分组划分,而不是按像元划分 """ gss = GroupShuffleSplit(n_splits=1, test_size=test_size, random_state=random_state) train_idx, test_idx = next(gss.split(X, y, groups=group_ids)) print(f"训练区块数: {len(np.unique(group_ids[train_idx]))}") print(f"验证区块数: {len(np.unique(group_ids[test_idx]))}") return train_idx, test_idxGroupShuffleSplit是scikit-learn里专门处理分组划分的工具。它的核心参数是groups,必须传入每个样本所属的空间区块编号。区块编号怎么生成?如果研究区是规则的,直接按经纬度划分网格;如果不规则,用地质边界多边形文件给每个像元打标签。这个设计阻止了邻近像元同时出现在训练集和验证集里,获得精度更接近真实泛化水平。
3.3 特征筛选:风险随递归特征消除RFE
原始特征几十上百个,先过滤后包装是效率最高的策略。系统里有两个筛选器,一个方差阈值(VarianceThreshold),一个递归特征消除(RFE)。方差阈值的作用是剔除那些在所有样本上几乎没变化的波段——这类波段对分类没有任何贡献,纯浪费计算资源。注意方差阈值之前必须先做归一化,否则量纲影响阈值选择。
from sklearn.feature_selection import VarianceThreshold, RFECV from sklearn.ensemble import ExtraTreesClassifier # 第一步:方差阈值,剔除常量特征 selector_var = VarianceThreshold(threshold=0.01) X_filtered = selector_var.fit_transform(X_train) print(f"方差筛选: {X_train.shape[1]} -> {X_filtered.shape[1]}") # 第二步:RFE + 交叉验证自动选最优特征数 estimator = ExtraTreesClassifier( n_estimators=50, max_depth=10, random_state=42 ) selector_rfe = RFECV( estimator=estimator, step=1, cv=3, scoring='f1_macro', min_features_to_select=5, n_jobs=-1 ) selector_rfe.fit(X_filtered, y_train) print(f"RFE筛选后特征数: {selector_rfe.n_features_}")这段代码有两个关键细节。RFE的step=1是每轮删除一个特征,如果特征数很多(超过50),建议step调成5或10,否则训练过程极慢。scoring用f1_macro而不是accuracy,因为岩性类别经常不均衡——比如研究区主要是灰岩和砂岩,页岩露头很少,用准确率会偏向大类。ExtraTreesClassifier作为RFE内部的估算器没问题,因为它本身就带随机性,对特征删减不敏感。
RFE的局限是它只看特征子集在模型上的表现,不考虑特征之间的物理关联。比如两个波段高度相关,RFE可能任意删掉一个,哪怕它们各有物理意义。好在岩性分类最终不是靠单个波段,而是靠多个特征的综合,所以这个局限通常不影响整体精度。
4. 布谷鸟算法与粒子群优化调参:两个进化算法的实战对比
4.1 为什么要用进化算法调参,而不是网格搜索
极端随机树的主要超参数是n_estimators、max_depth、min_samples_split、min_samples_leaf、max_features。用网格搜索(GridSearchCV)的话,假设每个参数给5个候选值,5个参数就是5的五次方——3125组参数组合,每组都做5折交叉验证,极端随机树训练本身又快不了,总耗时以小时级起步。
更关键的是,这些超参数的响应面并非光滑的单峰函数,网格搜索这种均匀采样很容易跳过那些狭窄的、但确实更好的参数区域。布谷鸟算法和粒子群优化都属于群体智能方法,它们的特点是能自适应地往有希望的区域加密搜索。布谷鸟算法的全局搜索能力靠莱维飞行(Lévy flight)保证——长尾分布的随机步长偶尔产生大跳跃,避免陷入局部最优。
4.2 布谷鸟算法实现:宿主鸟巢替换与发现概率
布谷鸟算法的生物学隐喻是布谷鸟把自己的蛋下在宿主鸟巢里,宿主有一定概率发现外来蛋并丢弃。算法维护N个鸟巢(候选解),每代做三件事:莱维飞行生成新解、按发现概率Pa替换部分差解、保留最优解。核心参数就两个:种群规模N和发现概率Pa。N一般取15-25,Pa取0.25——太小收敛慢,太大容易丢失好解。
import numpy as np from sklearn.ensemble import ExtraTreesClassifier from sklearn.model_selection import cross_val_score def objective(params, X, y): """目标函数:参数组 -> 交叉验证F1分数""" n_estimators, max_depth, min_samples_split, max_features = params model = ExtraTreesClassifier( n_estimators=int(n_estimators), max_depth=int(max_depth), min_samples_split=int(min_samples_split), max_features=float(max_features), random_state=42, n_jobs=-1 ) scores = cross_val_score(model, X, y, cv=3, scoring='f1_macro', n_jobs=-1) return scores.mean() def cuckoo_search(X, y, bounds, n_nests=20, pa=0.25, max_iter=30): """ bounds: 每个参数的[min, max]列表 pa: 外来蛋被发现概率 """ dim = len(bounds) # 随机初始化鸟巢 nests = np.random.rand(n_nests, dim) for i in range(dim): nests[:, i] = bounds[i][0] + nests[:, i] * (bounds[i][1] - bounds[i][0]) fitness = np.array([objective(nests[i], X, y) for i in range(n_nests)]) best_idx = np.argmax(fitness) best_nest = nests[best_idx].copy() best_fitness = fitness[best_idx] for _ in range(max_iter): # 莱维飞行生成新解 levy_steps = levy_flight(n_nests, dim) new_nests = nests + levy_steps # 边界处理 for i in range(dim): new_nests[:, i] = np.clip(new_nests[:, i], bounds[i][0], bounds[i][1]) # 评估新解,贪婪保留 new_fitness = np.array([objective(new_nests[i], X, y) for i in range(n_nests)]) improve_mask = new_fitness > fitness nests[improve_mask] = new_nests[improve_mask] fitness[improve_mask] = new_fitness[improve_mask] # 发现概率Pa:随机淘汰部分差解并重新初始化 rand_mask = np.random.rand(n_nests) < pa nests[rand_mask] = np.random.rand(rand_mask.sum(), dim) for i in range(dim): nests[rand_mask, i] = ( bounds[i][0] + nests[rand_mask, i] * (bounds[i][1] - bounds[i][0]) ) fitness[rand_mask] = np.array([objective(nests[i], X, y) for i in np.where(rand_mask)[0]]) # 更新全局最优 best_idx = np.argmax(fitness) if fitness[best_idx] > best_fitness: best_fitness = fitness[best_idx] best_nest = nests[best_idx].copy() return best_nest, best_fitnesslevy_flight的实现可以用Mantegna算法生成——用标准正态分布构造长尾分布步长。注意莱维飞行生成的步长可能非常大,所以紧接着必须做np.clip做边界约束。还有一个参数细节:n_estimators是整数,但优化器在连续空间搜索,所以目标函数里用int()强转。max_features传入浮点数给scikit-learn,它会被解释为特征比例,取值一般不要小于0.1,否则每棵树可用的特征太少,树之间的相关性上升,集成效果变差。
4.3 粒子群优化实现:惯性权重与认知社会参数
粒子群算法的思路是每个粒子维护速度和位置,根据个体最优(pbest)和全局最优(gbest)更新速度。参数有w(惯性权重)、c1(认知系数)、c2(社会系数)。经典取值是w=0.7,c1=c2=1.5。惯性权重控制粒子延续自身速度的程度,w太大粒子容易飞过最优区域,w太小粒子会快速汇聚但容易陷入局部最优。
def particle_swarm_optimization(X, y, bounds, n_particles=20, w=0.7, c1=1.5, c2=1.5, max_iter=30): dim = len(bounds) # 初始化位置和速度 positions = np.random.rand(n_particles, dim) velocities = np.zeros((n_particles, dim)) for i in range(dim): positions[:, i] = bounds[i][0] + positions[:, i] * (bounds[i][1] - bounds[i][0]) fitness = np.array([objective(positions[i], X, y) for i in range(n_particles)]) pbest = positions.copy() pbest_fitness = fitness.copy() gbest_idx = np.argmax(fitness) gbest = positions[gbest_idx].copy() gbest_fitness = fitness[gbest_idx] for _ in range(max_iter): r1, r2 = np.random.rand(dim), np.random.rand(dim) # 速度更新 velocities = (w * velocities + c1 * r1 * (pbest - positions) + c2 * r2 * (gbest - positions)) # 位置更新 + 边界约束 positions = positions + velocities for i in range(dim): positions[:, i] = np.clip(positions[:, i], bounds[i][0], bounds[i][1]) # 评估并更新个体最优和全局最优 fitness = np.array([objective(positions[i], X, y) for i in range(n_particles)]) better_mask = fitness > pbest_fitness pbest[better_mask] = positions[better_mask] pbest_fitness[better_mask] = fitness[better_mask] if fitness.max() > gbest_fitness: gbest_idx = np.argmax(fitness) gbest = positions[gbest_idx].copy() gbest_fitness = fitness[gbest_idx] return gbest, gbest_fitnessPSO这段代码里有一个常见错误容易被忽略:r1和r2只在循环开头生成一次,但在标准PSO里它们应该对每个粒子、每个维度都独立随机。这里简化成用同一个r1广播,实际工程中建议改成np.random.rand(n_particles, dim)。这个细节会影响粒子群的探索多样性,问题不大但属于可以改进的部分。
对比两段代码会发现:布谷鸟算法没有速度概念,全靠莱维飞行的重尾跳跃探索,对参数空间维度的敏感度低;PSO利用历史最优信息,收敛更平滑,但初期参数设置不对容易早熟。在我实际测试中,极端随机树的超参数响应面比较崎岖,CS在30代内通常能找到比PSO高0.5-1.5个百分点F1的解。
5. 避坑篇:岩性分类实战中最常见的五个坑
5.1 坑一:影像校正不规范导致分类器学到的是大气而非岩性
现象:训练精度很高(98%以上),但把模型应用到相邻区域的影像上精度骤降,或者分类结果图上同一种岩性在不同时相影像上被分成不同类别。
原因:直接用DN值或只做了辐射定标没做大气校正,模型把大气散射和吸收的差异当作特征学习了。不同时相影像的大气状态不同,模型学到的规律自然失效。
解决:预处理流程固定为辐射定标→大气校正→地形校正。快速验证方法:取同一区域两期不同日期的影像,提取同一点位的反射率曲线,如果曲线形状差异明显,说明校正还不到位。
5.2 坑二:空间数据泄漏导致验证精度虚高
现象:随机划分验证精度95%,换成按区块划分后精度直接掉到75%。
原因:遥感影像像元之间存在空间自相关,相邻像元的特征高度相似。随机划分让训练集和验证集包含了几乎重叠的空间信息,验证结果不能代表真实泛化能力。
解决:强制用GroupShuffleSplit按空间区块划分训练集和验证集。如果条件允许,更严格的做法是用地质界线做分组——同一个地质单元内的像元全部归到同一组。
5.3 坑三:布谷鸟算法早熟收敛,所有鸟巢挤到同一个局部最优
现象:连续迭代10代以上最优解没有变化,但把Pa调大重新跑又能继续提升。
原因:布谷鸟算法的探索能力依赖莱维飞行的重尾跳跃和发现概率Pa的随机淘汰。当所有鸟巢都收敛到同一个区域时,莱维飞行的小步长只能在局部搜索,Pa太小则差解不会被淘汰,群体失去了多样性。
解决:把Pa从0.25提到0.4(对ExtraTrees超参数寻优足够),同时限制每个参数的搜索边界不要太宽。另外,每隔5代把所有鸟巢重新初始化一半,这是血泪经验——比单纯调Pa有效得多。
5.4 坑四:类别不均衡导致稀有岩性被吞掉
现象:总体精度不错,但混淆矩阵里某些稀有岩类(比如研究区只占5%的角岩)召回率几乎是0。
原因:默认的ExtraTreesClassifier使用多数投票,多数类样本擠壓少数类的分类边界。
解决:给分类器设置class_weight='balanced_subsample',让每棵树的bootstrap采样自动平衡类别;或者在目标函数里用f1_macro而不是accuracy作为优化指标,这样调参时会同时照顾少数类。
5.5 坑五:高光谱波段数太多,树模型训练极慢且精度反而下降
现象:100+波段全量输入,极端随机树训练一小时没跑完,分类结果还出现大量椒盐噪声。
原因:波段间高度相关,冗余信息干扰了树的特征选择,同时高维特征导致决策树拟合时更不稳定。
解决:先做MNF降维到30维以内,或先用方差阈值+相关性分析去掉高度相关的波段。实际经验是:多光谱数据保留原始波段,高光谱数据降维到20-40维时精度和效率的平衡最好。分类完成后再做一次窗口平滑去椒盐噪声。
6. 极端随机树分类与后处理:特征重要性、混淆矩阵与空间平滑
6.1 ExtraTreesClassifier的超参数边界与默认值
极端随机树和随机森林的差别在于:随机森林在每个分裂节点从特征子集中找最优分裂点,极端随机树则是随机生成一个分裂阈值,然后选最好的一个。后者的好处是进一步降低方差,对噪声特征更稳健,在光谱数据这种强相关的场景下通常略优于随机森林。系统里的分类核心就是它,参数建议范围如下:
| 参数 | 建议范围 | 说明 |
|---|---|---|
| n_estimators | 200-500 | 太多收益饱和,300足够 |
| max_depth | 15-30 | 树太深容易过拟合噪声像元 |
| min_samples_split | 2-10 | 岩性分类建议取5以上 |
| min_samples_leaf | 1-5 | 防止单像元叶节点,配合后处理 |
| max_features | 0.1-0.5 | 每棵树看的特征比例 |
| class_weight | balanced_subsample | 类别不平衡时强烈建议 |
验证阶段直接输出混淆矩阵和分类报告:
from sklearn.ensemble import ExtraTreesClassifier from sklearn.metrics import confusion_matrix, classification_report best_params, best_f1 = cuckoo_search(X_scaled, y, bounds) model = ExtraTreesClassifier( n_estimators=int(best_params[0]), max_depth=int(best_params[1]), min_samples_split=int(best_params[2]), max_features=float(best_params[3]), class_weight='balanced_subsample', random_state=42, n_jobs=-1 ) model.fit(X_train, y_train) y_pred = model.predict(X_test) print(classification_report(y_test, y_pred)) print(confusion_matrix(y_test, y_pred)) # 输出特征重要性,前10个特征名 importances = model.feature_importances_ top_idx = np.argsort(importances)[::-1][:10] for idx in top_idx: print(f"{feature_names[idx]}: {importances[idx]:.4f}")6.2 后处理:众数滤波与地质约束
分类结果图上最多的噪声就是椒盐效应——单个像元被错分。原因是高光谱影像的像元尺度和地质体边界不匹配,加上地形阴影的影响。我一般用两种后处理:第一种是众数滤波,用一个3×3或5×5窗口统计中心像元邻域内的多数类别;第二种是面积约束——小于某个面积阈值的独立图斑直接替换为周围多数类别。后者在做地质填图时尤其重要,因为地质体不会碎成几个像元的碎片。
from scipy.ndimage import median_filter def majority_filter(class_map, size=5): """ 注意:类别是整数标签,不能直接用中值滤波 需要用众数滤波——逐像元统计窗口内出现次数最多的类别 """ from collections import Counter from scipy.ndimage import generic_filter def majority(values): return Counter(values).most_common(1)[0][0] return generic_filter(class_map, majority, size=size)这里有一个容易踩坑的细节:直接用median_filter处理类别标签是错的。中值滤波适合连续数值,类别标签用中值会生成不存在的类别编号。generic_filter配合自定义majority函数才能正确处理。窗口大小选5比较稳妥,窗口太大(9以上)会磨掉细小的地质界线。
6.3 多模型对比与不确定性输出
最后补一个验证习惯:我不只跑极端随机树一个模型,还会跑随机森林和梯度提升做对照。同一套训练集、验证集、同一套划分方式,谁的分高不一定说明谁好,但差异很大时一定是数据或流程出了问题。另外,ExtraTreesClassifier的predict_proba可以输出概率,概率低的像元意味着分类置信度低,可以单独输出成一张不确定性图,后续野外验证时优先布点。这套系统的整体流程走完,从影像输入到岩性图输出的完整链路是通的。
说句心里话,这套系统我前前后后跑了不下十遍,最深的教训是参数调优永远排在数据质量后面。从那以后我每次做遥感岩性分类都强制先做一遍反射率校正和空间不重叠划分,再谈算法和调参。资源里的脚本和代码组织得比较清楚,对照自己的数据改改路径和波段名就能跑通,希望帮到你。
本文还有配套的精品资源,点击获取