简介基于PCA-LDA的光谱数据降维处理MATLAB源码为高光谱图像分析、特征提取与分类识别场景提供了一套可直接运行的实现方案尤其适合需要处理高维光谱数据、降低冗余并提升模型区分度的研究人员与工程师。源码将无监督的主成分分析与有监督的线性判别分析串联使用先借助前者去除噪声、压缩数据维度再借助后者最大化类间差异使投影后的特征更利于后续分类器训练兼顾了运算效率与判别精度。压缩包共8个文件含4个M程序文件和4个ASV自动保存备份文件程序文件包括主程序、数据训练和测试等模块结构清晰便于对照学习算法流程和进行二次修改。整个资源压缩后仅6KB体积轻量已有2211人学习浏览。通过源码包读者可以清晰掌握组合降维的实现步骤、关键函数调用与参数调整方法并可直接迁移到自己的光谱分类或高光谱图像预处理任务中具有较高的参考和复用价值。1. 高光谱图像为什么一定要降维一个晚上只跑完一轮训练的教训接过一个高光谱图像分类的题目时我踩过最重的坑就是“原样上模型”。一块 512×512 像素、200 个波段的 Cube按像元展开是 26 万个样本乘 200 维特征丢进 SVM 做交叉验证一个晚上只跑完一轮训练结果还过拟合。后来我才意识到高光谱维度灾难的核心不是内存不够而是样本量和特征维度的比例完全失衡。你在每个像素上同时拿到的是上百个连续波段的反射率特征之间高度相关噪声、水分吸收带、仪器响应全都叠在里面。PCA-LDA 光谱数据降维的思路就是先用 PCA 做无监督压缩扔掉冗余和噪声再用 LDA 在低维空间里找真正能把类别分开的方向两步串在一起是高光谱图像分类里最实用、解释成本最低的一条路。这篇笔记把整个流程展开从数据读取到 K 值选择再到那些让你彻底翻车的边界问题全部写清楚。不管你是刚进这个方向的研究生还是做农产品无损检测、地物分类的工程师按这套流程走完至少能在一个晚上跑通全流程。2. 数据准备把高光谱立方体拆成“样本×波段”矩阵2.1 为什么几乎所有教程都不告诉你第一步如此关键高光谱图像不是一张普通照片它是一个三维立方体行、列、波段。常规 RGB 图像只有 3 个通道而高光谱图像往往有 100 到 1000 个波段。读取之后你没法直接拿立方体去做 PCA 或 LDA必须把它重排成二维矩阵每一行是一个像元样本每一列是一个波段特征行数可能是几十万列数是波段数。这一步听起来简单但很多新手第一次加载完数据就直接瞎了因为矩阵里充满了 NaN、全零行、饱和像素。这些脏数据不清理PCA 的协方差矩阵会被拉偏后面 LDA 更是直接奇异最后得到的结果连你自己都不敢信。还有一个非常容易被忽略的问题如果做的是监督分类标签必须和样本一一对应标签不齐LDA 算类内散度矩阵的时候会直接报索引越界。所以数据准备是整个 PCA-LDA 降维流程里最不能省的一步我一般会花掉整个项目三分之一的时间在这里但回报是后面的代码几乎不会因为数据问题返工。2.2 读取立方体并转成训练矩阵的 MATLAB 代码常见的数据格式有 ENVI.hdr .dat/.raw、MATLAB 的 .mat、HDF5 等。这里以 .mat 格式为例写一段通用的读取和矩阵化代码% 读取高光谱立方体: cube 是 rows × cols × bands 的三维数组 data load(hsi_cube.mat); cube data.cube; [rows, cols, bands] size(cube); % 把三维立方体重排为二维矩阵: 每行是一个像元的光谱 X reshape(cube, rows * cols, bands); % 清洗异常值: 全零像元、含 NaN/Inf 的像元都会污染协方差矩阵 intensity mean(X, 2); valid isfinite(sum(X, 2)) intensity 0; X X(valid, :); % 生成标签: 这里模拟二分类实际项目中替换为你的标签文件 labels load(hsi_labels.mat); y labels.gt; % ground truth 矩阵 y y(:); % 拉直成向量 y y(valid); % 与清洗后的样本对齐这段代码里有一个很容易犯的错reshape 是按列优先展开的也就是说先固定第一列的所有行再取第二列。如果你后续要做空间验证就必须保留行列坐标否则将来想按区域划分训练集时根本找不到样本在哪。清洗异常值这里我用的是亮度滤波全零像元通常是黑框或背景含 NaN 或 Inf 的像元通常来自传感器坏点或辐射定标异常。有人会问为什么不做最大值滤波但实际场景里饱和像元只在特定波段出现直接删掉整个像元损失信息太多我一般只做全零和 NaN 清洗然后把极端值用中位数填充。如果你需要做辐射校正或白板校正要在这之前完成不能让校正后出现异常值。2.3 训练集测试集划分与标准化这是数据泄漏的重灾区划分训练集和测试集时最常见的做法是按像元随机分层抽样。但注意高光谱图像的空间自相关性极强相邻像元的光谱几乎一样如果随机采样测试集里会有大量训练集像元的“近亲”测出来的精度虚高。所以如果你要做严谨的验证最好像遥感领域常见做法那样按空间块划分训练集把左上角区域做训练右下角区域做测试或者至少做区域划分后再随机抽。这会在第 5 章展开这里先给一份可跑的划分和标准化代码rng(42); % 固定随机种子确保结果可复现 cv cvpartition(y, HoldOut, 0.3); X_train X(training(cv), :); y_train y(training(cv)); X_test X(test(cv), :); y_test y(test(cv)); % 标准化: 只使用训练集的均值和标准差 mu mean(X_train, 1); sigma std(X_train, 0, 1); X_train_norm (X_train - mu) ./ sigma; X_test_norm (X_test - mu) ./ sigma;最后三行是整个预处理里最容易被写错的地方。很多人会图省事把 X_train 和 X_test 拼起来一起做 zscore然后训练精度极高、测试精度也很高但这种精度没有任何实际价值因为测试集的信息已经被“泄漏”到训练过程里了。正确的做法是永远只在训练集上计算 mu 和 sigma然后用同一组参数去标准化测试集。同理后面 PCA 的载荷矩阵也只能从训练集学习测试集必须用训练集的投影矩阵去变换。记住一句话测试集只能在最后被使用一次。3. PCA 无监督降维从协方差矩阵到主成分选择3.1 为什么 PCA 的核心是协方差矩阵的特征分解PCA 的原理在所有资料里都能看到但很多人只记住了“最大方差方向”这句话却不知道为什么一定要算协方差矩阵。其实逻辑很简单你对数据做了中心化后协方差矩阵里每个元素表示两个波段之间的线性相关程度对角线上是每个波段自身的方差。PCA 就是在寻找一组新的正交基使得数据在这组基上的投影方差尽可能大而这个优化问题最终就归结为协方差矩阵的特征值分解特征向量就是新基的方向特征值就是该方向上的方差大小。高光谱波段之间本来就高度相关比如 700nm 和 710nm 两个波段的反射率几乎完全线性相关所以协方差矩阵的秩很低前少数几个特征值就能解释绝大部分总方差——这正是高光谱数据能降维的根本原因。你可能已经发现相关热词里经常出现“pca 原理为什么用协方差矩阵”这类问题。这里顺带说清楚如果你对数据做了标准化协方差矩阵就等于相关矩阵但如果你没有标准化两者的结果会差异很大。高光谱数据的波段值通常是反射率或辐射值量纲一致但动态范围不同比如叶绿素吸收波段数值低近红外波段数值高。不标准化直接做 PCA高方差的波段会强行占据前几个主成分而真正能区分地物的细微光谱差异反而被埋没。所以我的建议是无论量纲是否一致先标准化再 PCA几乎总不会错。3.2 手写 PCA 与 MATLAB 内置 pca 函数的对照既然要理解 PCA 的行为就得知道它内部做了什么。先用矩阵运算手写一遍再和内置函数对照% 手写 PCA: 输入为训练集标准化后的矩阵 X_train_norm [Ns, B] size(X_train_norm); % 计算协方差矩阵 Sigma (X_train_norm * X_train_norm) / (Ns - 1); % 特征值分解 [V, D] eig(Sigma); % eig 返回的特征值按升序排列这里翻转为降序 [evals, idx] sort(diag(D), descend); V V(:, idx); % 取前 K 个主成分 K 20; Vk V(:, 1:K); % 降维后的得分矩阵 score_train X_train_norm * Vk; % 计算各主成分解释方差百分比与累计值 explained evals / sum(evals) * 100; cum_explained cumsum(explained);MATLAB 内置的 pca 函数用的是奇异值分解数值稳定性更好实际项目中我建议直接用 pca手写这段代码是为了让你看到方差解释率的来源。内置函数的等效调用是[coeff, score, ~, ~, explained] pca(X_train_norm); % coeff: 载荷矩阵波段 × 主成分score: 降维后得分这里要特别提醒pca 函数默认对数据进行中心化但不会做标准化。如果你前面已经手动 zscore那没问题如果没有必须通过 Centered 和 Standardize 参数显式控制。我见过太多人在 MATLAB 里直接 pca(X)发现前两个主成分永远是方差最大的波段而非真正有判别力的波段然后反过来怀疑 PCA 没用其实是标准化没做。3.3 K 值怎么选90% 还是 95%碎石图还是玄学这是 PCA 降维里最让人犹豫的一步选择 K 的方式直接影响后续 LDA 的表现。高光谱图像场景下常见做法有两种一种是用累计解释方差阈值一般取 95% 以上因为高光谱波段虽然多但有效信息往往压缩在前十几个主成分里另一种是看碎石图的拐点特征值陡降之后趋于平缓的位置就是 K。我一般会把两个方法结合起来% 画累计解释方差曲线 figure; plot(cum_explained(1:min(50, length(cum_explained))), o-, LineWidth, 1.5); xlabel(主成分数量); ylabel(累计解释方差百分比); grid on; title(选择 K 的碎石图); % 找到累计解释方差达到95%的K K_95 find(cum_explained 95, 1); fprintf(达到95%%解释方差所需主成分数: %d\n, K_95);在标准高光谱数据集上这个 K 通常在 10 到 30 之间。但注意累计解释方差 95% 是一个经验值不是数学上保证分类性能的阈值。有些光谱差异很小但非常有判别力的地物比如两个不同品种的作物它们的信息可能落在特征值相对靠后的主成分上。所以如果后续 LDA 分类精度一直上不去可以试着把 K 调大到 50 甚至 80观察精度变化。K 的选择本质上带一点玄学成分最好的办法是在后面接上分类器做交叉验证直接以分类精度作为 K 的判定标准。这个技巧我会在第 5 章的避坑清单里放一段完整的上层调参代码。4. LDA 监督降维利用标签找判别方向以及奇异矩阵的硬伤4.1 LDA 与 PCA 的本质区别最大化的是“可分离性”PCA 是完全没有使用标签的无监督方法它找的方向是数据方差最大的方向LDA 是监督方法它找的方向是类间散度最大、类内散度最小的方向。换句话说PCA 在回答“数据本身长什么样”LDA 在回答“不同类别之间有什么区别”。对高光谱分类而言LDA 天然更适合作为分类前的降维手段因为它直接面向判别任务。LDA 的核心是两个散度矩阵类内散度矩阵 S_w 和类间散度矩阵 S_b。S_w 是所有类别的样本围绕各自类别中心散布的累积S_b 是各分类中心围绕全局中心的散布并且每个类按样本量加权。它的目标可以写成 Fisher 判别准则的形式最大化 trace((W^T S_w W)^{-1} W^T S_b W)这个问题的解是广义特征值分解 S_b w λ S_w w 的前 C-1 个特征向量其中 C 是类别数。也就是说不管你的原始特征有多少维LDA 最多只能降到 C-1 维。二分类问题降维到 1 维三分类问题降维到 2 维。这是 LDA 的极限也是很多人觉得 LDA 降维后信息量不够的根源。4.2 高光谱场景下 LDA 为什么必须依赖 PCA 先降维直接对标准化后的原始光谱 X_train_norm 做 LDA你会立刻撞上一个硬伤S_w 是奇异的根本无法求逆。原因在于样本数远小于波段数。假设你有 1000 个训练像元、200 个波段S_w 是一个 200×200 的矩阵但它的秩最多只能达到 1000 减去类别数也就是远小于 200。秩亏矩阵没有逆你用 MATLAB 的 eig(S_b, S_w) 得到的结果会是一堆 NaN 或者虚数。这就是为什么标题里写的是“PCA-LDA”而不是单独的 LDAPCA 先把数据压缩到最多不超过样本数的低维空间通常控制 K 在 30 以内此时 S_w 满秩LDA 才能稳定工作。如果不先做 PCA几乎所有高光谱数据都会在这个位置翻车。4.3 在 PCA 压缩后的低维空间上手写 LDA这里给出一个手写 LDA 的函数它接收 PCA 降维后的得分矩阵和训练标签返回投影矩阵function [Wlda] mylda(X, y) % 输入: X 是样本×特征的低维得分矩阵, y 是列向量标签 classes unique(y); C numel(classes); [n, d] size(X); mu_all mean(X, 1); Sw zeros(d); Sb zeros(d); for c 1:C Xc X(y classes(c), :); nc size(Xc, 1); muc mean(Xc, 1); % 类内散度: 各类样本减去类均值后的外积和 Xc_center Xc - muc; Sw Sw Xc_center * Xc_center; % 类间散度: 类中心相对全局中心的偏移, 按样本量加权 diff muc - mu_all; Sb Sb nc * (diff * diff); end % 广义特征值分解: 解 Sb * w lambda * Sw * w [V, D] eig(Sb, Sw); [evals, idx] sort(real(diag(D)), descend); V V(:, idx); % LDA 最多保留 C-1 个判别方向 Wlda V(:, 1:C-1); % 丢弃虚部: 当 Sw 接近奇异时会出现极小的虚部 Wlda real(Wlda); end代码里最需要注意的就是最后两行。高光谱数据经过 PCA 后虽然 Sw 理论上是满秩的但数值上某些特征值非常接近零eig 的广义求解会对这些小特征值极其敏感结果会带出复数虚部。直接拿复数矩阵去算后面的投影分类器会报错。我的处理是排序时取实部投影矩阵再强制 real()这样能在数值上稳定工作。此外这里的 Sw 是散度矩阵不是协方差矩阵所以没有除以样本数。除以不除以不影响特征向量的方向但会影响特征值的大小对后续投影不产生实质影响。4.4 完成降维后接分类器fitcdiscr 与精度验证PCA-LDA 的最终目的不是降维本身而是为分类器提供更好的输入。LDA 投影之后的特征维度非常低二分类只有 1 维此时最自然的搭配是线性判别分类器也就是 fitcdiscr% 用训练集做 PCA 获得投影矩阵和降维得分 [coeff, score_train, ~, ~, ~] pca(X_train_norm); score_train_k score_train(:, 1:K); % 用 LDA 找到判别投影方向 Wlda mylda(score_train_k, y_train); % 训练分类器 X_lda_train score_train_k * Wlda; mdl fitcdiscr(X_lda_train, y_train); % 测试集变换: 先做同样标准化, 再 PCA, 再 LDA score_test_k (X_test_norm * coeff(:, 1:K)); X_lda_test score_test_k * Wlda; % 预测与精度 pred predict(mdl, X_lda_test); acc mean(pred y_test); fprintf(PCA-LDA 分类精度: %.2f%%\n, acc * 100);测试集变换时有一个地方非常关键用 X_test_norm 乘以 coeff 的时候coeff 必须来自训练集的 pca 结果而不是对测试集重新做 pca。如果对测试集重新拟合 PCA相当于测试集的光谱信息被用于构造投影矩阵这会造成与标准化同源的数据泄漏问题。我在实际项目中把 coeff 保存为 .mat 文件推理阶段只加载不下算这样能从流程上防止这种错误。5. PCA-LDA 实战中的避坑与排查从维度报错到精度上不去5.1 LDA 报错“矩阵接近奇异值”或结果全是 NaN在做 LDA 时最常见的一个报错是矩阵接近奇异或者广义特征值分解后结果全是 NaN。原因是直接拿原始 200 维光谱矩阵喂给 mylda类内散度矩阵 Sw 的秩远小于维度是不可逆的。解决方法是必须先用 PCA 降维让 K 小于训练样本数减去类别数。我通常把 K 设置在 20 到 50 之间也就是在累计解释方差 95% 附近这样能将光谱压缩到 LDA 可计算的维度区间。如果你把 K 设到 200LDA 照样报错因为 PCA 只是改变了坐标系的旋转并没有减少特征数量除非你显式截断主成分数量。这是我第一次做高光谱降维时踩过的坑当时以为 PCA 和 LDA 连用就万事大吉结果 PCA 后忘记截断主成分LDA 直接求逆失败。5.2 训练精度 99%测试精度一塌糊涂数据泄漏怎么排查训练精度极高、测试精度低的问题在高光谱图像上最容易被误判为“模型过拟合”但实际上往往是流程上的数据泄漏。最常见的原因是标准化的时候把训练集和测试集拼在一起算了均值和方差导致测试集的信息进入了训练阶段。另一个更隐蔽的原因是 PCA 载荷矩阵用全量数据拟合而不是仅用训练集。排查方式很简单检查 mu、sigma 和 coeff 这三个变量是在哪一步计算出来的如果它们是基于 X_train_norm 还是基于全量 X问题一眼就能看出来。还有一点就是本章 2.3 中提到的空间自相关导致的样本泄漏两个相邻像元的光谱几乎一致随机划分后测试集里混入训练集邻居测试精度虚高。这两种泄漏都要从流程上堵住而不是靠调模型参数解决。5.3 LDA 降维后精度反而不如单独 PCA维度压太狠了在二分类任务中LDA 最多只能输出 1 个判别方向所有信息被投影到一条线上。如果这条线的方向稍微被噪声带偏分类精度就会断崖式下降。所以当你发现 PCA-LDA 的精度低于单独 PCA 加 SVM 时第一反应不应该是怀疑 LDA 没用而是应该检查 K 值是多少。如果 K 只有 5主成分里包含的判别信息可能不够如果 K 是 30信息通常已经足够。我会做一个简单的 K 值扫描来验证而不是直接拍脑袋选一个数% 在 K 的候选值上做 5 折交叉验证以分类精度为准则 K_candidates 5:5:60; cv_acc zeros(length(K_candidates), 1); for i 1:length(K_candidates) k K_candidates(i); % 先用 pca 降维到 k [coeff_k, score_k] pca(X_train_norm, NumComponents, k); % 再 LDA 投影 Wlda mylda(score_k, y_train); score_lda score_k * Wlda; % 用 fitcdiscr 做五折交叉验证 mdl fitcdiscr(score_lda, y_train, CrossVal, on); cv_acc(i) 1 - kfoldLoss(mdl); fprintf(K%2d, CV acc%.4f\n, k, cv_acc(i)); end % 选交叉验证精度最高的 K [~, best_idx] max(cv_acc); best_K K_candidates(best_idx);这种循环跑起来非常快因为每个 K 上只做 PCA 到 k 维再 LDA模型训练开销极低。比直接猜 K 靠谱得多。你会发现最优 K 通常不是一个精确的点而是一个平缓的平台取平台中间的值就好不必追求极致精度。5.4 数据量太少、类别不平衡时LDA 的方向被大类绑架高光谱图像的标签数据通常非常稀缺一块图上千类中只有几片区域有标注如果某类样本只有几十个像元而另一类有上万个像元LDA 的目标函数会因为 S_b 中按样本量加权的项被大类主导。此时 LDA 找出的方向会重点区分大类与全局中心小类的判别信息被严重忽略。解决方式有两个方向一是对各类训练样本做数量均衡常见做法是随机欠采样或 SMOTE 过采样但高光谱数据上我一般只做欠采样另一种是修改 LDA 的类间权重项把按样本量权重改成均等权重。具体改法就是在算 Sb 的时候不乘 nc把每个类别等权对待。代价是牺牲大类精度来换小类召回率但分类任务本来就是宁可整体精度略降也要每个类别都别漏掉这个调整在农产品分类任务里经常能救命。5.5 预处理阶段中文注释乱码导致变量被误删MATLAB 编码环境惹的祸这个问题看起来和算法无关但我在实际项目里被它浪费过一整个下午。MATLAB 2023 及之后的版本默认使用 UTF-8 编码但很多在老版本中以中文本地化编码GBK保存的脚本文件打开后中文注释会变成乱码。乱码本身不影响运行但如果注释文字里恰好包含换行符被破坏代码块的 if/end 配对就会错乱你可能会在删乱码时误删掉关键变量初始化语句。解决办法很简单在 MATLAB 主页设置里把语言编码切换为 UTF-8或者一步到位写代码时变量名用英文注释用英文或尽量少的短中文。数据预处理脚本中的任何一行错误都会被后续 PCA-LDA 逐步放大最终结论完全不可信。6. 进阶用 PCA 载荷和 LDA 判别向量反选波段给模型找回物理意义PCA-LDA 作为降维手段的输出是一组抽象的主成分很多人把它当成一个“黑匣子”算出分类精度就算完事不敢去解释它。其实降维除了让分类器跑起来之外还有另一个高价值的应用就是通过载荷矩阵反推原始波段的重要性挑选出物理上可解释的代表波段。这在噪声检测、地物识别、农产品品质检测场景里比单纯追求分类精度更有工程价值因为你可以告诉别人“我们主要用 680nm 和 720nm 附近的波段就能把这两类分开”而不是说“我们用前 15 个主成分”。做法是直接把 LDA 的判别方向投影回原始波段空间。PCA 的载荷矩阵 coeff 描述的是每个主成分由各原始波段线性组合的权重LDA 的方向则描述各个主成分的判别权重两者相乘得到每个原始波段对最终判别方向的贡献。贡献值绝对值越大的波段说明它对分类越重要% 假设 K 已经通过 5.3 中的交叉验证选好 [coeff, ~, ~, ~, ~] pca(X_train_norm); score_k score_train_norm * coeff(:, 1:best_K); Wlda mylda(score_k, y_train); % 投影回原始波段空间: 原始波段 × 判别方向 proj coeff(:, 1:best_K) * Wlda; % 多类时 Wlda 有 C-1 列合并为每个波段的总体贡献 band_importance sum(abs(proj), 2); % 排序并取出前 15 个最重要的波段 [~, idx_top] sort(band_importance, descend); top_bands idx_top(1:15); % 可视化: 画出平均光谱并标记重要波段 wavelengths 400:10:2500; % 示例波段名称替换为你自己的 figure; plot(wavelengths, mean(X_train_norm), k, LineWidth, 1.2); hold on; stem(wavelengths(top_bands), band_importance(top_bands), r); xlabel(波长 (nm)); ylabel(反射率 / 贡献权重); legend(平均光谱, 重要波段);代码里需要注意一点band_importance 是 PCA 载荷和 LDA 方向乘积绝对值求和的结果它代表的是线性贡献权重不是相关显著性。如果两类样本的判别信息主要存在于红边区域比如 680nm 到 750nm这个权重曲线通常会在那个区间出现一个清晰的峰这和植被遥感里的“红边效应”能够相互印证。遇到这种情况算是运气好模型和物理规律对上了。如果权重峰杂乱无章那就得对 PCA-LDA 流程本身打一个问号K 是否选得太小把关键主成分丢弃了或者训练集里存在异常值没有清洗干净。我自己的习惯是每次选完波段后拿选出的 15 个波段直接训练一个决策树或线性 SVM如果精度和 PCA-LDA 全流程的精度差距在三个百分点以内说明波段选择是有效的如果差距过大说明关键信息被遗漏了需要回头调整 K 或检查数据质量。这个方法把降维模型从一个黑匣子变成可解释的工具也让结果在汇报时更容易被信任。最后分享一个我自己的习惯每次跑完 PCA-LDA我都会把 coeff、Wlda 和标准化参数 mu、sigma 一并存成 .mat 文件并写一行注释记录 K 值和当时的交叉验证精度。这样下次换一批数据进来不用从头调参直接在已有投影方向的基础上做迁移尝试能省下大把时间。高光谱降维这件事真正的壁垒不是你懂多少原理而是你能不能把数据清洗、参数选择、评估方法这套流程固定下来不翻车。希望这篇笔记能帮到你。本文还有配套的精品资源点击获取
