偏最小二乘法与支持向量机在高光谱分类中的应用实践

简介面向高光谱图像分类应用的Matlab程序包集成偏最小二乘法PLS与支持向量机SVM两条技术路线并兼顾BP神经网络等对比方法主要服务于遥感、农业和环境监测等领域的研究人员帮助解决高光谱数据降维、特征提取与地物分类问题。压缩包共含22个文件以mat数据文件14个为主配合6个m脚本和2个asv备份文件既有训练/测试数据集也有可直接运行的分类主程序与自定义核函数实现整体约100KB小巧但结构清晰。目前已吸引213人学习适合初入高光谱分类或希望对比多种算法效果的学习者。通过运行程序可完整体验从数据预处理、特征提取到分类器构建与结果验证的流程尤其能直观理解PLS降维后结合SVM分类的常见协作模式为实际项目中的地物识别提供可复用的代码基础与调参思路。1. 高光谱分类为什么把偏最小二乘法和支持向量机放在一起拿到一份高光谱数据第一件事往往不是训练分类器而是面对一个尴尬的现实波段数动辄几百样本量却常常只有几十到几百个。这种“高维小样本”结构让很多常规分类算法直接失效——特征维数远大于样本数时协方差矩阵不可逆距离度量失真模型很容易在训练集上完美分类、在测试集上彻底崩溃。高光谱分类的核心矛盾从来不是“用什么分类器”而是“怎么在几百个高度相关的波段里提取出真正有判别力的信息”。偏最小二乘法PLS和支持向量机SVM恰好是从两个不同方向解决这个问题的经典组合。PLS走的是“降维回归/分类”路线它把光谱矩阵和类别标签同时投影到低维潜变量空间找的是最能解释类别差异的方向SVM则走“间隔最大化”路线通过核函数把样本映射到高维空间再找分类超平面天然擅长处理小样本问题。前者是线性方法计算快、可解释性强后者可以非线性上限高、鲁棒性好。实际项目里把两者的预测结果对比着看或者用PLS做特征提取再喂给SVM是出现频率最高的落地套路。这篇博文就按这个思路把数据预处理、PLS-DA分类、SVM分类、参数调优和常见坑过一遍所有代码都能直接跑。2. 高光谱分类的预处理与样本划分反射率转换和训练集构造是第一步2.1 高光谱DN值到反射率的转换高光谱传感器记录的原始数据通常是DN值Digital Number受光照条件、传感器响应、大气吸收等因素影响不同时间、不同批次获取的DN值之间没有可比性。分类模型学习的是光谱形态差异如果直接用DN值建模光照变化会被模型误当成地物类别差异训练集和测试集来自不同拍摄条件时准确率会大幅下降。高光谱如何转反射率常见做法有两种有定标参数时做辐射定标和大气校正把DN值转为地表反射率没有定标条件时至少做相对归一化。实操中比较稳妥的处理是如果数据集中有白板或参考板用目标像元DN值除以参考板DN值得到反射率近似值如果没有参考板则对每条光谱做向量归一化消除整体亮度差异。下面给出一段将高光谱DN值转为反射率并进行归一化的Python代码import numpy as np def dn_to_reflectance(dn_img, white_refNone, norm_methodl2): # dn_img: 形状为 (rows, cols, bands) 的原始DN值数据 # white_ref: 可选白板参考光谱形状为 (bands,) if white_ref is not None: # 用白板做相对反射率转换 reflectance dn_img / (white_ref 1e-6) else: # 没有白板时先做最小值偏移再按照波段最大值缩放 min_val dn_img.min(axis(0, 1), keepdimsTrue) max_val dn_img.max(axis(0, 1), keepdimsTrue) reflectance (dn_img - min_val) / (max_val - min_val 1e-6) # 对每条光谱做归一化消除亮度差异 if norm_method l2: norm np.linalg.norm(reflectance, axis-1, keepdimsTrue) reflectance reflectance / (norm 1e-6) elif norm_method mean: mean_val reflectance.mean(axis-1, keepdimsTrue) reflectance reflectance / (mean_val 1e-6) return reflectance这段代码的核心逻辑是如果有白板参考光谱就做逐波段的比值运算得到近似反射率没有白板时做Min-Max拉伸保证每个波段的数值范围在0到1之间。最后的归一化非常重要——高光谱分类对光照强度差异极其敏感L2归一化把每条光谱变成单位向量只保留形状信息丢掉强度信息这能让模型更关注光谱形态而非整体亮度。实际使用中我一般倾向L2归一化它对植被、水体、土壤这类光谱形态差异明显的地物效果更好如果分析的是矿物等高反照率目标均值归一化会更稳定。2.2 训练集和测试集的划分策略高光谱分类中样本的划分直接影响模型评估的可信度。常见的错误做法是随机把所有像素点打乱后划分训练集和测试集这会导致同一地物区域的空间相邻像素同时出现在训练和测试中模型“记住”了空间位置而不是光谱特征评估结果虚高。更严谨的做法是依据空间区域划分——把图像按区域分块部分区域做训练其余区域做测试。以下是一个按空间位置进行分块划分的代码实现def spatial_split(data, labels, train_ratio0.5, split_axisx): # data: 形状为 (rows, cols, bands) # labels: 形状为 (rows, cols) 的地物类别标签 rows, cols data.shape[:2] train_mask np.zeros((rows, cols), dtypebool) test_mask np.zeros((rows, cols), dtypebool) if split_axis x: split_idx int(cols * train_ratio) train_mask[:, :split_idx] True test_mask[:, split_idx:] True elif split_axis y: split_idx int(rows * train_ratio) train_mask[:split_idx, :] True test_mask[split_idx:, :] True else: # 按区域网格划分每隔一个网格取一个区域 block_size 16 for r in range(0, rows, block_size): for c in range(0, cols, block_size): if ((r // block_size) (c // block_size)) % 2 0: train_mask[r:rblock_size, c:cblock_size] True else: test_mask[r:rblock_size, c:cblock_size] True train_indices np.where(train_mask) test_indices np.where(test_mask) X_train data[train_indices] y_train labels[train_indices] X_test data[test_indices] y_test labels[test_indices] # 去除标签为0的背景像素 valid_train y_train 0 valid_test y_test 0 return X_train[valid_train], y_train[valid_train], X_test[valid_test], y_test[valid_test]按空间分块划分的关键意义在于高光谱相邻像素的高度相关性决定了如果训练集和测试集存在空间重叠会导致模型评估结果偏乐观分块划分能够显著降低这种光谱泄漏风险。分块大小是个需要权衡的参数块太大时训练区域和测试区域的类别分布可能不均衡类别较少的小地物会被切掉块太小时又接近随机划分。网格划分的块尺寸建议设置在8到32像素之间同时观察每个区域内的类别分布确保每类样本在两个集合中都有足够数量。另外要注意背景或无效像素的空间分布通常有一定规律要确保它们不会被全部划入训练集否则会影响类别先验概率。3. 用偏最小二乘法做高光谱分类PLS-DA的最小可复现程序3.1 为什么偏最小二乘法适合处理高维光谱数据偏最小二乘法把光谱矩阵X样本×波段和类别标签矩阵Y样本×类别同时分解为潜变量空间中在提取X主成分时要求与Y的协方差最大。这个机制决定了它和主成分分析PCA有本质区别PCA只找X方差最大的方向这些方向很可能与类别判别无关而PLS找的是与类别标签相关性最强的方向在分类任务中信息利用效率更高。高光谱数据的波段之间存在严重的多重共线性相邻波段的反射率高度相关。PLS的潜变量投影在数学上是各波段的线性组合其权重向量体现了哪些波段对分类贡献最大这在高光谱分类中具有重要的物理意义——通过检查权重或VIPVariable Importance in Projection得分可以定位出与地物属性密切相关的特征波段。相比直接使用全部波段PLS-DA输出的类别概率也可以通过Softmax变换获得方便后续与SVM概率输出做对比。3.2 PLS-DA的Python实现与参数设定下面给出使用scikit-learn实现PLS-DA进行高光谱分类的完整代码import numpy as np from sklearn.cross_decomposition import PLSRegression from sklearn.preprocessing import LabelEncoder, StandardScaler from sklearn.metrics import classification_report, confusion_matrix from sklearn.model_selection import cross_val_score def pls_da_classify(X_train, y_train, X_test, y_test, n_components10): # 标签转为数值编码 le LabelEncoder() y_train_enc le.fit_transform(y_train) y_test_enc le.transform(y_test) n_classes len(le.classes_) # 构造二值化的类别矩阵用于PLS回归 Y_train np.zeros((len(y_train_enc), n_classes)) Y_train[np.arange(len(y_train_enc)), y_train_enc] 1 # 对光谱数据做标准化每个波段均值0方差1 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_test_scaled scaler.transform(X_test) # 创建PLS回归模型分类时本质是PLS回归到类别指示矩阵 pls PLSRegression(n_componentsn_components, scaleFalse) pls.fit(X_train_scaled, Y_train) # 预测并取概率最大的类别 Y_pred pls.predict(X_test_scaled) y_pred_enc np.argmax(Y_pred, axis1) y_pred le.inverse_transform(y_pred_enc) # 输出分类报告 print(classification_report(y_test, y_pred)) return pls, y_pred, le这段代码把PLS-DA拆解为回归任务类别标签被编码为one-hot矩阵PLS回归预测每个样本属于各类的得分取最大值对应的类别作为分类结果。有几个参数值得注意。n_components是PLS的潜变量数量也是这个模型最重要的超参数设置过小时欠拟合设置过大时会过度拟合训练集的噪声通常通过交叉验证来确定。scaleFalse的选择基于前面已经做了标准化避免重复计算。训练前用StandardScaler对光谱做标准化在这里几乎是必需的因为不同波段的反射率绝对值差异较大如果不标准化高反射率波段会在潜变量投影中占据过大的权重弱化低反射率但具有判别力的波段贡献。3.3 用交叉验证确定最佳潜变量数选择一个合适的潜变量数量最常见的做法是K折交叉验证观察不同成分数下模型在验证集上的平均准确率和标准差from sklearn.model_selection import StratifiedKFold def tune_pls_components(X_train, y_train, max_components30): # 标签转为数值编码 le LabelEncoder() y_enc le.fit_transform(y_train) n_classes len(le.classes_) Y_train np.zeros((len(y_enc), n_classes)) Y_train[np.arange(len(y_enc)), y_enc] 1 scaler StandardScaler() X_scaled scaler.fit_transform(X_train) # 使用分层K折确保每个折中类别比例一致 cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) best_n, best_score 2, -1 scores_history [] for n in range(2, max_components 1): pls PLSRegression(n_componentsn) scores cross_val_score(pls, X_scaled, Y_train, cvcv, scoringaccuracy) mean_score scores.mean() scores_history.append((n, mean_score, scores.std())) if mean_score best_score: best_n, best_score n, mean_score # 打印前10个成分数的验证结果 for n, mean_s, std_s in scores_history[:10]: print(fcomponents{n}: accuracy{mean_s:.4f} (/- {std_s:.4f})) return best_n, scores_history在实验代码中常常出现一个误区直接选用验证集准确率最高的成分数忽略了标准差的影响。高光谱数据噪声较大某个成分数可能在某个折上表现突出换一个折就明显下降泛化性并不好。一般建议选择在“准确率开始进入平台期”的最小成分数这时模型复杂度最低、泛化性最好。如果PLOT中发现准确率随成分数一直上升且没有平台期多半是训练集和测试集之间存在空间相关性样本划分逻辑需要重新检查。4. 支持向量机高光谱分类实现核函数选择与参数网格搜索4.1 SVM在光谱分类中的适用条件SVM在处理小样本、高维数据时的优势在于它不依赖数据分布的统计假设通过间隔最大化原则寻找最优分类超平面泛化误差上界只与间隔大小有关与特征维数无直接关系。这就解释了为什么在波段数几百、样本数几十的条件下SVM仍能训练出有效的分类器而逻辑回归或线性判别分析在同样的数据规模下往往因为协方差矩阵估计不稳定而失效。SVM的决定边界由支持向量决定也就是距离分类超平面最近的少数样本点。这意味着高光谱数据中大量“容易分类”的像素不会影响模型只有处于类别边界附近、容易混淆的像素才真正决定分类结果。这种稀疏性让SVM对高光谱数据的噪声相对鲁棒但也带来了一个短板对类别不平衡比较敏感。如果某一类像素数量极少支持向量可能全部落在多数类一侧少数的地物类别会被直接忽略。实际使用中需要预先处理类别不平衡比如调整类别权重或者对少数类做过采样。4.2 RBF核SVM的网格搜索参数RBF核函数有两个关键参数C正则化参数和gamma核函数宽度。C控制对误分类样本的惩罚力度越大越容易过拟合gamma控制单个训练样本的影响范围越大则决策边界越复杂。对光谱数据gamma的含义要结合输入范围理解——标准化后的光谱数值通常在0到1之间此时gamma在0.001到0.1这个区间比较常见太大会导致每个样本只影响自己附近极小的范围决策边界变得极其破碎。以下是基于网格搜索的SVM调参和高光谱分类实现import numpy as np from sklearn.svm import SVC from sklearn.model_selection import GridSearchCV, StratifiedKFold from sklearn.preprocessing import StandardScaler from sklearn.metrics import accuracy_score, cohen_kappa_score def svm_grid_search(X_train, y_train, X_test, y_test): # 标准化光谱数据 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_test_scaled scaler.transform(X_test) # 定义参数网格C和gamma按数量级搜索 param_grid { C: [0.5, 1, 10, 50, 100], gamma: [0.001, 0.005, 0.01, 0.05, 0.1], class_weight: [balanced, None] } # RBF核SVM svm SVC(kernelrbf, probabilityTrue, random_state42) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) grid GridSearchCV(svm, param_grid, cvcv, scoringf1_macro, n_jobs-1, verbose1) grid.fit(X_train_scaled, y_train) print(Best params:, grid.best_params_) best_model grid.best_estimator_ y_pred best_model.predict(X_test_scaled) acc accuracy_score(y_test, y_pred) kappa cohen_kappa_score(y_test, y_pred) # 输出总体准确率和Kappa系数 print(fOverall Accuracy: {acc:.4f}) print(fKappa Coefficient: {kappa:.4f}) return best_model, y_pred, scaler网格搜索让模型自动寻找最优的参数组合但实际使用时要注意几个细节class_weightbalanced在类别不平衡时能显著提升小众类别的识别率。高光谱数据中的常见场景是水体、植被这类大块地物样本量充足而某种人工地物或稀有矿物只有几十个像素此时不给少数类加权的话模型的精确率和召回率会严重失衡。使用f1_macro作为评分指标比accuracy更稳妥它给每个类别相同的权重少数类的表现下降会直接反映在得分上。probabilityTrue会启用Platt缩放用交叉验证来估计概率代价是训练时间增加如果不需要概率输出可以关掉以加快训练速度。4.3 PLS降维加SVM的组合策略高光谱分类中一个经常出现的问题是数据中的冗余信息、波段之间存在的高度相关性干扰了SVM的核函数计算。一种有效的折中方案是先用PLS对光谱数据进行降维再把降维后的潜变量作为SVM的输入。这样做的好处是通过PLS降维剔除了部分噪声和非线性干扰SVM的高精度得到保留训练速度也有明显提升。def pls_svm_pipeline(X_train, y_train, X_test, y_test, n_components15): from sklearn.cross_decomposition import PLSRegression # 标签映射为one-hot le LabelEncoder() y_train_enc le.fit_transform(y_train) y_test_enc le.transform(y_test) Y_train np.zeros((len(y_train_enc), len(le.classes_))) Y_train[np.arange(len(y_train_enc)), y_train_enc] 1 # PLS降维 pls PLSRegression(n_componentsn_components) X_train_pls pls.fit_transform(X_train, Y_train)[0] X_test_pls pls.transform(X_test) # 对降维后的特征标准化并训练SVM scaler StandardScaler() X_train_pls_scaled scaler.fit_transform(X_train_pls) X_test_pls_scaled scaler.transform(X_test_pls) svm SVC(kernelrbf, C10, gamma0.05, class_weightbalanced) svm.fit(X_train_pls_scaled, y_train_enc) y_pred_enc svm.predict(X_test_pls_scaled) y_pred le.inverse_transform(y_pred_enc) print(PLS-SVM Accuracy:, accuracy_score(y_test, y_pred)) return svm, pls, y_predPLS-SVM组合的效果并不总优于直接使用原始光谱训练SVM——当PLS降维后的前几个潜变量保留了大部分判别信息时分类准确率通常没有明显下降但训练时间大幅缩短如果类别间差异主要集中在少数波段PLS的线性投影可能把这些信息稀释掉。建议在项目中同时跑“原始光谱SVM”和“PLS-SVM”两组实验对比结果再决定是否采用降维方案。4.4 高光谱分类评估指标与混淆矩阵分析高光谱分类评估最常遇到的困境是总体准确率Overall Accuracy, OA看着很高但某些类别完全没分出来。比如一个区域90%的面积是农田剩下10%是道路和建筑模型把所有像素都预测为农田OA也能达到90%但道路和建筑的类别完全失败。只看OA是不够的至少要结合以下指标判断指标计算方式高光谱场景下的典型问题总体准确率OA正确分类像素数/总像素数大类别占主导时容易虚高平均准确率AA各类别分类准确率的算术平均小类别表现差时能反映出来类别精确率该类预测正确数/预测为该类总数用于判断是否产生大量误报类别召回率该类预测正确数/该类实际总数用于判断该类是否被漏检Kappa系数基于混淆矩阵计算的一致性指标剔除了随机分类的干扰混淆矩阵逐类别的预测与实际对应表定位具体哪两类容易混淆在代码中输出这些指标的实际操作是在预测完成后from sklearn.metrics import classification_report, cohen_kappa_score, confusion_matrix # y_test为真实标签y_pred为模型预测 report classification_report(y_test, y_pred, digits3) print(report) print(Kappa:, cohen_kappa_score(y_test, y_pred)) # 输出混淆矩阵便于定位易混淆的类别对 cm confusion_matrix(y_test, y_pred) for i, row in enumerate(cm): print(fTrue class {y_classes[i]}: {row})面对高光谱分类结果中混淆矩阵显示出的模式需要系统性排查如果两个植被类别频繁混淆优先检查波段范围是否包含红边区域如果某个类别的召回率极低需要查看训练集样本量是否满足基本要求考虑补充训练样本或调整类别权重如果混淆集中在光谱形态相近的地物常见做法是把纹理特征加入SVM输入。分类评估在模型迭代中扮演着定位系统问题的参谋角色不是最终目的。目标是把混淆矩阵中体现的错误模式还原到光谱特征层面理解然后指导数据预处理和参数调整。5. 高光谱分类的5个实战技巧光谱特征、样本量与模型选择一起决定上限5.1 技巧一用主成分分析结果辅助划定感兴趣区域高光谱分类的第一步通常是确定感兴趣区域Region of InterestROI圈定训练样本。手动圈选ROI时容易凭借假彩色合成图的主观经验忽视了光谱数据中的细微差异。实际操作中先在所有光谱上做主成分分析观察前几个主成分分量不同地物通常在特定主成分分量上呈现明显对比度差异。比如PCA第二分量突出显示土壤成分差异第三分量突出显示植被长势差异利用这些差异来圈定ROI比直接看图例选择更准确。5.2 技巧二光谱重采样减少噪声波段干扰高光谱传感器的相邻波段相关性极高很多波段实际是由光谱仪内部插值造成的冗余信息。对这些波段直接建模不仅计算量大还会放大探测器噪声的影响。常见做法是先通过VIP得分或SVM-RFE递归特征消除筛选出重要性较高的波段子集再做分类。RFE-SVM在实验中的波段筛选效果通常好于随机选择但计算时间较长。高光谱400到2500纳米范围内一般保留15到30个代表性波段就可以保留大部分分类能力。5.3 技巧三少数类的样本扩增策略高光谱数据中类别不平衡是常态。当某一类的训练样本只有几个到几十个时SVM几乎无法学习该类别的有效边界。常用的解决办法是对少数类样本做小幅度的光谱扰动扩增在原光谱上叠加标准差为1%到2%的噪声生成新的样本。这种方式在实验中被证实能有效提升SVM对少数类的分类性能但要注意噪声幅度不能超过光谱本身的类内方差否则会让扩增样本偏离真实分布。5.4 技巧四多模型集成避免单一分类器偏好PLS-DA优点在于处理线性分类和快速特征提取SVM擅长处理非线性边界。把两者的预测概率做加权平均或采用投票法组合实际操作中得到的整体分类精度一般比单一分类器提升2到5个百分点。集成方式需要注意各类别在两个模型上的表现差异比如PLS-DA在判别土壤类别上表现好SVM在判别植被亚类上更准确那在投票时给对应的类别预测更高的权重会更合理。5.5 技巧五把分类结果保存为地理编码栅格高光谱分类的最终产出是一张分类图方便后续做面积统计和专题图制图。把像素级的预测结果按照原始空间位置填回去保存为GeoTIFF格式from osgeo import gdal def save_classification_map(pred_labels, rows, cols, geo_transform, proj, output_path): # pred_labels: 一维预测结果长度等于 rows*cols driver gdal.GetDriverByName(GTiff) ds driver.Create(output_path, cols, rows, 1, gdal.GDT_Byte) ds.SetGeoTransform(geo_transform) ds.SetProjection(proj) band ds.GetRasterBand(1) # 将一维预测结果重塑为二维图像 band.WriteArray(pred_labels.reshape(rows, cols)) # 设置NoData值 band.SetNoDataValue(0) ds.FlushCache() ds None print(fClassification map saved to {output_path})保存分类结果时注意把背景像素重新填为0与标准的地物类别编码保持一致。后续在GIS软件中加载分类图时类别编码的可读性直接决定分析效率建议在输出时同时保存一个类名与编码的对应关系表为后续修正和验证提供依据。高光谱分类程序调试到最后拼的不是模型有多复杂而是把这类数据处理环环相扣的工程细节做扎实。本文还有配套的精品资源点击获取