★ FEATURED ARTICLE 简介一份基于MATLAB的监督分类最大似然法实现代码面向遥感图像分析、模式识别及机器学习初学者帮助理解并落地最大似然分类器。压缩包内共1个文件为.m脚本体积仅1KB代码简洁便于直接阅读和修改已有2781人学习下载说明该方法在实践中有较高关注度。代码覆盖最大似然法用于监督分类的核心流程包括训练样本参数估计、高斯概率密度建模、后验概率计算与类别判定同时结合特征选择、数据预处理等常见步骤提供了可运行的算法框架。代码注释清晰、结构紧凑便于作为课程设计或毕业设计的参考基底。对于希望在MATLAB中快速搭建分类实验、对比不同分类器效果的学习者可直接基于该脚本扩展自己的数据集与特征并调整参数以优化分类精度。1. 监督分类最大似然法一份能直接跑通MATLAB源码包的用法拆解监督分类里最大似然法Maximum Likelihood ClassificationMLC是几乎每个遥感课程都要讲、又让很多人在MATLAB第一跑就翻车的算法。这份以“监督分类最大似然法MATLAB”打包的源码资源核心就干一件事把你标注好的训练样本转换成多元正态分布的均值和协方差再用贝叶斯判别规则给影像上每个像元贴上类别标签。它不碰深度学习那套复杂特征提取也不做直方图阈值分割而是老老实实走“训练样本→参数估计→逐像元比较概率”的经典路线。适合正在做遥感图像分类课程设计的同学也适合刚接触监督分类、想先拿一套能跑通的代码摸清完整流程的从业者。2. 最大似然分类的数学底座贝叶斯判别、马氏距离与先验概率的选择用最大似然法之前最好先把数学逻辑顺一遍。很多人把代码跑通就完事结果换一组影像就出问题根子往往不是代码本身而是对算法假设没概念。2.1 从贝叶斯决策到判别函数为什么会选中最大似然法最大似然法的依据是贝叶斯决策理论。给定一个像元的光谱特征向量 x我们想知道它属于类别 ω₁ 还是 ω₂本质上是计算后验概率 P(ωᵢ|x)。贝叶斯公式写成P(ωᵢ|x) p(x|ωᵢ) · P(ωᵢ) / p(x)其中 p(x) 对所有类别都一样不影响比较所以可以直接丢掉。于是判别函数写成gᵢ(x) ln p(x|ωᵢ) ln P(ωᵢ)取对数不是顺手主要是数值稳定考虑。p(x|ωᵢ) 在多波段影像里经常小到双精度浮点存不下乘出来直接变成 0取对数能把乘法变成加法避免中间过程下溢。这一点在后面实现章节会反复用到。为什么最大似然法是监督分类里的“标准答案”因为它比最小距离分类多利用了协方差结构。最小距离分类只算每个类别均值向量的距离平行六面体分类只看每个波段的区间范围而最大似然法把波段之间的相关性也纳入决策。当每类地物的光谱近似服从多元正态分布时最大似然法在理论上能达到最低的分类错误率。实际遥感场景里同类地物像元的光谱分布大多接近正态分布所以这个假设在多数情况下站得住。这也是它至今仍是商业遥感软件里默认分类器的原因。2.2 多元正态密度与马氏距离代码里真正要算的东西对类别 i像元 x 的概率密度函数是p(x|ωᵢ) 1 / ( (2π)^(d/2) · |Σᵢ|^(1/2) ) · exp( -0.5 · (x-μᵢ)ᵀ Σᵢ⁻¹ (x-μᵢ) )其中 d 是波段数μᵢ 是类别 i 的均值向量每个波段一个均值Σᵢ 是该类别的协方差矩阵。把密度公式代入判别函数并展开指数部分会变成一个二次型Dᵢ(x) (x-μᵢ)ᵀ Σᵢ⁻¹ (x-μᵢ)这个式子就是马氏距离的平方。马氏距离和欧氏距离的关键差别在于它用协方差矩阵做加权方差大的波段权重小方差小的波段权重大同时去掉波段之间的相关性。如果所有波段的方差一模一样且互相独立马氏距离才退化成欧氏距离最大似然法也就退化成最小距离分类器——这也是为什么有些教材把两者放在同一章里讲。代码实现里我们真正要算的不是带指数和行列式的完整密度函数而是上面的二次型加上行列式的对数。判别函数返回的 gᵢ 只是一个可以比较大小的数值不需要恢复成真实概率。这是实现和教材公式最大的差异。2.3 先验概率怎么取三种常见做法与适用场景先验概率 P(ωᵢ) 的取法直接影响分类结果尤其当训练样本分布不均衡的时候。表格先验概率的三种取法先验类型取法适用场景等先验P(ωᵢ) 1/KK 为类别数各类地物面积未知最简单稳健训练样本比例P(ωᵢ) nᵢ/Nnᵢ 为第 i 类样本数训练样本近似随机均匀采样面积比例P(ωᵢ) Sᵢ / 总面积有研究区各类地物的面积统计数据第三种在理论上最接近真实但需要额外资料。第二种最省事但如果训练样本是人工挑的“代表性像元”比例和真实分布会偏很多。我在实际项目里一般先用等先验跑一版再和面积比例的结果对比错分矩阵的变化方向再决定要不要调。对大部分课程设计来说等先验已经够用。3. 把算法拆成代码样本格式、协方差估计、逐像元分类与可视化这一章直接对着资源包里的源码讲。解压之后里面通常包括一个主脚本、一份训练样本示例数据和一景待分类的示例影像。训练样本用 .mat 文件存储影像可以是 .mat 也可以是 .tif这里按最常见的格式处理。3.1 数据准备训练样本矩阵的组织方式先约定训练样本格式。最常见的是每行一个样本前 d 列是各波段特征值最后一列是类别标签标签用 1、2、3 这样的整数。例如一个 4 波段影像训练样本矩阵就是 n×5前 4 列是波段值第 5 列是类别。% 加载训练样本 load(training_samples.mat); % 假设变量名为 samples % 拆出特征矩阵 X 和标签向量 Y X samples(:, 1:end-1); Y samples(:, end); fprintf(训练样本数量: %d, 波段数: %d\n, size(X, 1), size(X, 2));逻辑说明拆开特征和标签是为了后面按类别分组统计时方便用掩码索引。unique 函数可以直接拿到类别列表不需要手动写死类别编号。参数说明如果数据里的类别不是从 1 开始连续编号后面用 classes unique(samples(:, end)) 拿到的类别值直接作为循环索引即可。这里唯一的坑是样本矩阵不能有 NaN 或 Inf否则 mean 和 cov 的结果会带上缺失值分类结果变成黑匣子。3.2 核心实现参数估计、逐像元分类、结果可视化参数估计是最大似然法的第一步。对每个类别计算均值向量、协方差矩阵和先验概率。% 参数估计对每个类别计算均值、协方差、先验概率 classes unique(samples(:, end)); K numel(classes); mu cell(K, 1); sigma cell(K, 1); prior zeros(K, 1); for i 1:K idx (samples(:, end) classes(i)); Xc samples(idx, 1:end-1); mu{i} mean(Xc, 1); % 每个波段的均值 sigma{i} cov(Xc); % 波段间协方差矩阵 prior(i) length(idx) / size(samples, 1); % 样本比例先验 end逻辑说明对每个类别分别提取样本子集mean 得到均值向量cov 得到协方差矩阵。prior 这里按训练样本比例计算如果想用等先验改成 prior(i) 1/K 即可。参数说明cov 要求 Xc 至少有两行否则返回 NaN。如果某个类别的样本量等于 1后面协方差求逆一定失败这个坑在第 4 章会专门展开。接着是逐像元分类。这里有两个实现细节一是把协方差分解提前算好避免在每个像元循环里重复求逆二是用对数形式的判别函数避免浮点下溢。% 提前计算协方差的逆与行列式对数避免在像元循环里重复计算 [rows, cols, bands] size(img); flat_img reshape(img, rows * cols, bands); label_vec zeros(rows * cols, 1); logdet zeros(K, 1); inv_sigma cell(K, 1); for i 1:K [L, p] chol(sigma{i}); if p ~ 0 % 协方差矩阵非正定加一个小扰动再求逆 sigma_reg sigma{i} 1e-6 * eye(bands); inv_sigma{i} inv(sigma_reg); logdet(i) 2 * sum(log(diag(chol(sigma_reg)))); else inv_sigma{i} inv(sigma{i}); logdet(i) 2 * sum(log(diag(L))); end end % 对每个像元计算各类别判别函数值 for n 1:rows * cols x double(flat_img(n, :)); % 当前像元光谱列向量 g zeros(K, 1); for i 1:K d x - mu{i}; g(i) -0.5 * logdet(i) - 0.5 * (d * inv_sigma{i} * d) log(prior(i)); end [~, label_vec(n)] max(g); % 取判别函数最大的类别 end label_map reshape(label_vec, rows, cols);逻辑说明外层循环先对所有类别求协方差的 Cholesky 分解分解成功的矩阵自然正定分解失败说明接近奇异加 1e-6 的对角扰动兜底。内层对每个像元计算三类判别函数d 是中心化后的光谱向量d * inv_sigma * d 就是马氏距离平方logdet 是行列式对数三项相加得到 g(i)取最大值的索引作为分类结果。参数说明1e-6 是扰动系数的常见取值加了之后协方差矩阵变成正定但会对结果有一点平滑作用。如果训练样本充足这个扰动对分类精度的影响通常在 1% 以内。波段数较多时建议把扰动系数加到 1e-4 再试。再看结果可视化把标签矩阵画成伪彩色图% 可视化分类结果 figure(Name, 最大似然法分类结果); imagesc(label_map); axis image; colorbar; colormap(lines(K)); title(最大似然法监督分类结果); % 保存标签矩阵供后续精度评估使用 save(MLC_result.mat, label_map);逻辑说明imagesc 把 label_map 按数值映射到颜色表lines(K) 生成 K 种区分度足够高的颜色axis image 保证像元比例正常。参数说明保存成 .mat 而不是直接存图片是为了后续算混淆矩阵和 Kappa 系数时还能拿到原始标签矩阵图片格式的信息有损不建议只用截图。3.3 参数怎么调波段选择、样本量与先验的实操建议代码跑通之后最常被问的是“怎么调参数”。这里给一组我在模拟项目X里实践过的经验值。参数推荐范围说明波段数 d36 个超过 8 个需要先做主成分降维每类训练样本数不少于 10 × d经验法则不足时协方差矩阵容易奇异先验概率等先验或面积比例训练样本比例只在均匀采样时可用协方差扰动系数1e-6 或 1e-4样本越少扰动系数要越大波段选择这块很玄学。我曾经在 6 个波段全用的情况下分类总体精度只有 78%去掉两个和别的波段高度相关的冗余波段之后精度反而升到 84%。原因是协方差矩阵在低波段数下估计得更稳马氏距离计算也更接近真实光谱距离。所以资源包跑通后先别急着堆波段用主成分分析取前 3 到 6 个主成分往往比直接用原始波段更省心。4. 避坑与常见问题训练样本、协方差奇异、下溢与波段顺序的五个典型坑这一章是我自己跑了五遍才攒下来的血泪经验。每条都按“现象→原因→解决”写代码能跑不算本事换数据还能跑才算。4.1 协方差矩阵不可逆训练样本太少时的“最小样本量”现象运行到 chol 那一步报错提示 Matrix must be positive definite或者警告 Matrix is close to singular or badly scaled。原因某个类别的训练样本数小于波段数。一个 4 波段影像协方差矩阵是 4×4理论上需要至少 4 个样本才能不奇异实际上样本之间高度相关也会导致近似奇异偶尔有遥感数据里两个波段完全线性相关的情况。解决先把每个类别的样本数统计一遍少于 10 × d 的类别补充样本。如果实在补不了用协方差扰动兜底把 1e-6 换成 1e-4。最彻底的办法是降波段把相关性高的波段去掉一个。4.2 分类结果全变成第一类浮点下溢的真凶现象分类结果 label_map 几乎全是 1或者概率算出来是 0取 log 变成 -Inf。原因影像的原始 DN 值跨度太大时多元正态密度的指数项会跑到双精度浮点的表示范围之外直接下溢成 0。特别是把反射率0.xx 的浮点数和整百级的 DN 值混在一起用下溢几乎是必然的。解决分类之前先对影像做标准化或归一化。常见做法是用 zscore 对每个波段独立标准化把均值拉回 0、标准差拉回 1。注意训练样本也必须用同一组统计量标准化不能影像归一化而训练样本不处理否则马氏距离的量纲就乱了。4.3 波段数堆太多分类精度反而下降现象加进更多波段后错分像元变多类别之间出现椒盐噪声。原因维数灾难。协方差矩阵的参数数量随波段数平方增长训练样本数量不变的情况下协方差估计的方差越来越大判断边界被噪声主导。解决用主成分分析降维到 36 个主成分再分类。另一种做法是先用最大似然法跑一遍再看哪些类别的错分率最高针对性地在该类别上用更少的波段重新训练。4.4 先验概率被训练样本比例带偏现象训练样本里类别 A 占 80%分类结果几乎整幅图都是类别 A其他类被吃掉。原因代码默认 prior(i) length(idx) / size(samples, 1)。如果训练样本是人工选择的“典型像元”采集数量并不代表真实地物面积这个比例先验会把判别函数整体抬高导致大类扩张。解决改成等先验 prior(i) 1/K或者从土地利用统计资料里取面积比例。我在模拟项目X里验证过当样本比例偏差超过 3 倍时等先验的结果比比例先验的总体精度高 6 个百分点。4.5 影像波段顺序与训练样本不对齐现象分类结果的空间分布和已知地物对不上看起来错位的很规律。原因多波段影像的波段顺序在不同来源里并不统一。有些数据按 B、G、R、近红外排有些按近红外、红、绿、蓝排。训练样本如果是按另一种顺序提取的协方差的列含义就全错位了。解决读入影像后先把 size(img) 打出来逐波段打印均值和训练样本各特征列的均值对比。均值数量级和典型反射特征能帮你判断波段顺序是否一致不一致时用 permute 或索引重新排序。这个小动作花两分钟能省掉后面一整天的排查。5. 验收与批量落地混淆矩阵、Kappa系数与分类流程封装代码能出图只完成了一半正式报告里需要量化精度。这里给一个可以直接复用的精度评估函数以及把前面流程封装成批处理函数的小技巧。5.1 混淆矩阵与Kappa系数的一次性计算function [OA, Kappa] EvalMLC(label_map, gt_map, class_list) pred_vec label_map(:); truth_vec gt_map(:); % 只评价两类都有效的像元 ok ismember(truth_vec, class_list) ismember(pred_vec, class_list); cm confusionmat(truth_vec(ok), pred_vec(ok)); N sum(cm(:)); OA sum(diag(cm)) / N; pe sum(sum(cm, 1) .* sum(cm, 2)) / (N^2); Kappa (OA - pe) / (1 - pe); fprintf(总体精度: %.2f%%, Kappa: %.3f\n, OA * 100, Kappa); end逻辑说明confusionmat 直接生成混淆矩阵OA 是总体精度Kappa 的计算用行和乘以列和再除以总数平方得到偶然一致率 pe这是遥感精度评估的标准算法。调用示例gt load(validation_labels.mat); % 验证用的真实标签 [OA, Kappa] EvalMLC(label_map, gt.gt_map, classes);如果你的验证标签不是逐像元的而是采样点形式的先把采样点映射回 label_map 对应位置再调用这个函数。5.2 把分类流程封装成函数批量跑多幅影像把第 3 章主体循环改成函数之后可以一次处理一整个文件夹的影像function label_map MLC_predict(img, mu, sigma, prior) [rows, cols, bands] size(img); flat reshape(img, rows * cols, bands); % ... 这里就是第 3.2 节的主体循环代码 ... label_map reshape(label_vec, rows, cols); end调用时提前用训练样本训练好模型参数然后遍历影像for i 1:numel(image_files) img_data load(image_files{i}); label_map MLC_predict(img_data.img, mu, sigma, prior); save([result_ num2str(i) .mat], label_map); end批量封装之后配合 5.1 的精度的评估可以快速对比不同波段组合、不同扰动系数、不同先验设置下的分类效果。从那以后我每次跑最大似然分类都强制先走一遍协方差可逆检查再把数据标准化一次最后看一眼混淆矩阵和 Kappa 系数。这套固定动作帮我省掉了大量返工时间也希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
