简介这份资源提供ISODATA聚类算法的MATLAB实现代码面向从事遥感图像处理、模式识别与无监督分类的学习者和研究人员。ISODATA即迭代自组织数据分析技术是对k-means聚类算法的改进能够克服k-means需预先指定类别数的缺点通过分裂与合并机制自动调整聚类数目在多光谱图像的无监督分类任务中表现尤为有效。资源包共6个文件包含cpp与h源码文件、说明文档、README及若干工程配置文件压缩包约6KB体量轻便便于快速阅读与移植。目前已有485人学习下载。读者可从中获取完整的算法实现逻辑理解分裂、合并、迭代阈值等核心步骤的代码组织方式并可直接将其应用于多光谱遥感影像分类实验或作为进一步改进聚类算法的参考基础。1. 多光谱图像分类为什么绕不开 ISODATA 这步手里有一批多光谱影像波段从可见光到近红外排了七八层但没有任何标注样本也没人告诉你地表到底分几类。这种场景下ISODATAIterative Self-Organizing Data Analysis Technique几乎是绕不开的一步。它和 K-means 最大的区别在于聚类数 K 不是拍脑袋定的而是在迭代中通过分裂和合并自动调整最后收敛到一个相对合理的类别数。这对多光谱数据特别友好因为波段多、地物混杂你很难事先知道该分几类。这篇笔记面向两类人一是刚拿到多光谱数据、想快速跑出分类结果的新手二是用过 K-means 但被固定 K 值坑过、想换一套自适应方案的熟手。我会把 ISODATA 的核心参数、MATLAB 实现步骤、调参逻辑和踩坑记录讲清楚代码可以直接抄参数可以照着改。2. ISODATA 的四个核心机制与 MATLAB 选型理由2.1 分裂、合并、删除三个动作撑起自适应聚类ISODATA 的迭代过程本质上是三个动作的循环分裂、合并、删除。分裂针对的是某个簇的标准差在某个维度上超过阈值说明这个簇内部太散需要拆成两个合并针对的是两个簇中心距离小于阈值说明它们其实是一类需要合到一起删除针对的是某个簇的样本数太少低于最小样本数阈值直接丢掉。这三个动作的触发条件都依赖参数设定。分裂看的是各维度标准差和分裂阈值的关系合并看的是簇间距离和合并阈值的关系删除看的是簇内样本数和最小样本数的关系。理解这三个动作的触发逻辑比死记代码更重要因为调参调的就是这些阈值。多光谱数据的特点是波段间相关性高但不同地物的光谱曲线差异集中在特定波段。比如植被在近红外波段反射率高水体在近红外波段反射率低建筑在可见光波段反射率均匀。ISODATA 在分裂时按维度检查标准差恰好能捕捉到这种波段间的差异把混在一起的地物拆开。2.2 为什么在多光谱场景选 ISODATA 而不是 K-meansK-means 需要预先指定 K 值这在多光谱无监督分类里是个硬伤。你可能会说那我用肘部法或者轮廓系数选 K 不就行了问题是多光谱数据维度高肘部法的拐点往往不明显轮廓系数计算量也大。ISODATA 把选 K 的问题转化为调阈值的问题虽然阈值也不好调但至少它有一个自适应的过程初始 K 可以给一个粗略值让算法自己去分裂和合并。另一个原因是 ISODATA 对初始中心不敏感。K-means 对初始中心很敏感不同的初始化可能得到完全不同的结果。ISODATA 因为允许分裂和合并即使初始中心选得不好迭代过程中也有机会调整过来。当然这不是说初始中心可以随便选只是容错空间更大。从计算量看ISODATA 每轮迭代要做一次分类、一次统计、一次分裂合并判断比 K-means 多了一些开销但在多光谱数据量级下比如几百乘几百像素、七八个波段MATLAB 里跑几十轮迭代也就是几十秒的事完全可以接受。2.3 MATLAB 实现的最小代码框架下面是一个可以直接跑的 ISODATA 核心框架输入是三维多光谱数据立方体输出是分类标签图和最终簇中心。function [labels, centers] isodata_core(data, init_k, max_iter, ... min_samples, split_std, merge_dist, max_k) % data: H x W x B 的多光谱数据 % init_k: 初始聚类数 % max_iter: 最大迭代次数 % min_samples: 每簇最小样本数 % split_std: 分裂标准差阈值 % merge_dist: 合并距离阈值 % max_k: 最大允许簇数 [H, W, B] size(data); X reshape(data, H*W, B); % 展平成 N x B N size(X, 1); % 初始化随机选 init_k 个样本作为初始中心 rng(42); % 固定随机种子保证可复现 idx randperm(N, init_k); centers X(idx, :); for iter 1:max_iter % 步骤1分配每个样本到最近中心 dists pdist2(X, centers); [~, labels] min(dists, [], 2); % 步骤2统计每簇的样本数和标准差 new_centers zeros(size(centers)); cluster_std zeros(size(centers)); valid_clusters []; for k 1:size(centers, 1) member X(labels k, :); if size(member, 1) min_samples new_centers(k, :) mean(member, 1); cluster_std(k, :) std(member, 0, 1); valid_clusters [valid_clusters, k]; end end % 步骤3删除小簇样本数不足的簇直接丢弃 centers new_centers(valid_clusters, :); cluster_std cluster_std(valid_clusters, :); % 步骤4分裂判断 split_flag false; for k 1:size(centers, 1) if max(cluster_std(k, :)) split_std ... size(centers, 1) max_k % 沿标准差最大的维度分裂 [~, split_dim] max(cluster_std(k, :)); offset zeros(1, B); offset(split_dim) split_std * 0.5; centers [centers; centers(k, :) offset]; centers(k, :) centers(k, :) - offset; split_flag true; break; % 每轮最多分裂一次避免震荡 end end % 步骤5合并判断 if size(centers, 1) 1 merge_flag false; for i 1:size(centers, 1)-1 for j i1:size(centers, 1) if norm(centers(i, :) - centers(j, :)) merge_dist % 合并两个簇中心为加权平均 centers(i, :) (centers(i, :) centers(j, :)) / 2; centers(j, :) []; merge_flag true; break; end end if merge_flag, break; end end end % 收敛判断中心变化小于阈值 if ~split_flag ~merge_flag fprintf(ISODATA 在第 %d 轮收敛最终簇数 %d\n, ... iter, size(centers, 1)); break; end end % 最终分配 dists pdist2(X, centers); [~, labels] min(dists, [], 2); labels reshape(labels, H, W); end这段代码的逻辑说明先把三维数据展平成二维矩阵每行是一个像素的光谱向量。初始化用随机选点的方式固定随机种子保证每次跑结果一致。主循环里先做分配再统计每簇的均值和标准差然后依次执行删除、分裂、合并。分裂时沿标准差最大的维度偏移偏移量取分裂阈值的一半这个系数可以调。合并时用两个中心的算术平均作为新中心也可以改成按样本数加权。参数说明init_k一般给 5 到 10给太小分裂次数多给太大合并次数多max_iter给 50 到 100 足够大多数情况二三十轮就收敛min_samples按总像素数的 0.5% 到 1% 给太小会导致碎片簇多太大会把真实的小地物删掉split_std是关键参数多光谱数据归一化后一般给 0.1 到 0.3merge_dist同样在归一化后给 0.2 到 0.5max_k按经验给 15 到 20防止无限分裂。2.4 数据预处理归一化和波段选择多光谱数据不同波段的量纲可能差很多比如某些波段是 0 到 255 的 DN 值某些是 0 到 1 的反射率。不归一化直接跑 ISODATA标准差大的波段会主导距离计算分裂判断也会偏向那个波段。常见做法是逐波段做 min-max 归一化把每个波段拉到 0 到 1 之间。% 逐波段 min-max 归一化 [H, W, B] size(data); data_norm zeros(H, W, B); for b 1:B band data(:, :, b); band_min min(band(:)); band_max max(band(:)); if band_max band_min data_norm(:, :, b) (band - band_min) / (band_max - band_min); else data_norm(:, :, b) band; % 常量波段保持原样 end end归一化之后分裂阈值和合并阈值才有可比性。如果跳过这一步你会发现分裂阈值给 0.1 根本没反应给 10 又疯狂分裂就是因为量纲没统一。波段选择方面如果波段数超过 10 个建议先做一次相关性分析把高度相关的波段去掉几个。ISODATA 的距离计算在高维空间会受冗余波段影响两个本来很近的像素因为冗余波段的噪声被拉远。常见做法是计算波段间相关系数矩阵相关系数大于 0.95 的波段对只保留一个。3. 从数据立方体到分类图完整跑通流程3.1 数据读取与维度检查多光谱数据常见的格式有 ENVI 的 .dat .hdr、GeoTIFF、MAT 文件。MATLAB 里读 GeoTIFF 用readgeoraster读 ENVI 格式需要自己解析头文件或者用第三方函数。不管哪种格式读进来之后第一件事是检查维度。% 以 MAT 文件为例 load(multispectral_data.mat); % 假设变量名是 ms_data [H, W, B] size(ms_data); fprintf(数据尺寸: %d x %d, 波段数: %d\n, H, W, B); % 检查是否有异常值 if any(isnan(ms_data(:))) warning(数据中存在 NaN将用波段均值填充); for b 1:B band ms_data(:, :, b); band_mean mean(band(~isnan(band))); band(isnan(band)) band_mean; ms_data(:, :, b) band; end end维度检查的目的是确认数据是 H x W x B 的排列。有些数据读进来是 B x H x W需要 permute 一下。NaN 值在多光谱数据里很常见比如云遮挡或者传感器异常直接跑 ISODATA 会导致整个簇的统计量变成 NaN所以要先填充。3.2 参数初始化与迭代控制参数初始化不是随便填几个数就完事需要根据数据特点来定。下面是一个参数初始化的参考表参数含义推荐范围调整方向init_k初始聚类数5-10数据复杂给大简单给小max_iter最大迭代次数50-100一般不用改min_samples最小样本数总像素的 0.5%-1%小地物多给小split_std分裂标准差阈值0.1-0.3分类细给小的merge_dist合并距离阈值0.2-0.5类别少给大的max_k最大簇数15-20地物复杂给大迭代控制里有一个容易忽略的点分裂和合并不能在同一轮里同时触发太多次否则会震荡。我一般每轮最多分裂一次、合并一次让算法慢慢调整。如果一轮里分裂三次又合并两次中心位置会剧烈变化收敛不了。% 迭代控制示例 options struct(); options.init_k 8; options.max_iter 80; options.min_samples round(H * W * 0.008); options.split_std 0.15; options.merge_dist 0.35; options.max_k 18; [labels, centers] isodata_core(ms_data_norm, options.init_k, ... options.max_iter, options.min_samples, options.split_std, ... options.merge_dist, options.max_k);3.3 结果可视化与类别数统计跑完之后需要看结果。最直接的是把标签图用伪彩色显示再统计每个类别的像素数和占比。% 可视化分类结果 figure; subplot(1, 2, 1); imagesc(labels); colormap(jet(max(labels(:)))); colorbar; title(ISODATA 分类结果); subplot(1, 2, 2); unique_labels unique(labels); counts histcounts(labels(:), [unique_labels; max(unique_labels)1]); bar(unique_labels, counts); xlabel(类别编号); ylabel(像素数); title(各类别像素分布); % 打印统计信息 fprintf(最终类别数: %d\n, length(unique_labels)); for i 1:length(unique_labels) fprintf(类别 %d: %d 像素 (%.2f%%)\n, ... unique_labels(i), counts(i), ... counts(i) / numel(labels) * 100); end看分类图的时候注意两点一是空间连续性同一类地物应该连成片如果出现大量椒盐噪声说明分裂阈值太小或者波段噪声太大二是类别占比如果某一类占了 80% 以上说明合并阈值太大把不同地物合到一起了。3.4 用混淆矩阵做无监督结果的间接验证无监督分类没有真实标签怎么验证常见做法是找一小块已知区域做目视判读或者用高分辨率影像做参考。如果手头有少量标注样本可以算一下聚类结果和标注的匹配度。% 假设有少量标注样本 ground_truth部分像素 % 计算聚类标签与真实标签的混淆矩阵 valid_idx ground_truth 0; pred_labels labels(valid_idx); true_labels ground_truth(valid_idx); % 用匈牙利算法做标签匹配MATLAB 自带 matchpairs cost_matrix zeros(max(pred_labels), max(true_labels)); for i 1:max(pred_labels) for j 1:max(true_labels) cost_matrix(i, j) -sum(pred_labels i true_labels j); end end matched matchpairs(cost_matrix, 1e6); % 计算匹配后的准确率 correct 0; for i 1:size(matched, 1) correct correct sum(pred_labels matched(i, 1) ... true_labels matched(i, 2)); end accuracy correct / length(true_labels); fprintf(匹配后准确率: %.2f%%\n, accuracy * 100);这段代码的逻辑是无监督聚类的类别编号和真实标签的编号没有对应关系所以先用匈牙利算法做最优匹配再算准确率。matchpairs是 MATLAB 自带的二分图匹配函数cost_matrix用负的匹配数因为matchpairs求的是最小代价。4. 调参踩坑记录分裂阈值、合并阈值和波段噪声4.1 分裂阈值给太小分类图全是椒盐噪声现象跑完 ISODATA 之后分类图上同一片地物里出现大量零散的不同类别像素看起来像撒了盐和胡椒。原因分裂阈值split_std给得太小比如给了 0.05导致每个簇只要有一点标准差就触发分裂最后簇数膨胀到 max_k 上限每个像素几乎自成一类。解决把split_std调回 0.15 到 0.25 之间同时检查数据是否归一化。如果归一化后标准差仍然很大说明波段噪声高可以先做一次 3x3 均值滤波再跑聚类。4.2 合并阈值给太大不同地物被强行合并现象分类结果只有三四个类别明显把植被和建筑合到了一起分类图上看不出细节。原因合并阈值merge_dist给得太大比如给了 0.8导致两个本来有差异的簇中心因为距离小于阈值被合并。解决把merge_dist降到 0.3 左右同时检查归一化方式。如果用的是 z-score 归一化距离尺度会和 min-max 归一化不同阈值需要重新调。我一般统一用 min-max 归一化阈值范围比较稳定。4.3 波段间量纲不统一距离计算被单一波段主导现象分类结果和某个特定波段的空间分布高度相似其他波段好像没起作用。原因没有做逐波段归一化某个波段的数值范围远大于其他波段距离计算时那个波段的贡献占了主导。解决跑之前先做逐波段 min-max 归一化再检查每个波段归一化后的均值和标准差是否接近。如果某个波段归一化后标准差仍然很小说明这个波段信息量低可以考虑去掉。4.4 初始聚类数给太少分裂次数过多导致收敛慢现象迭代了七八十轮还没收敛每次输出都在分裂簇数慢慢往上爬。原因init_k给得太小比如给了 2而实际地物有七八类算法需要不断分裂才能逼近合理类别数。解决把init_k提到 6 到 10 之间让初始中心覆盖大部分地物类型。初始中心可以用 k-means 的方式选比纯随机选更分散。4.5 最小样本数给太大小地物被删除现象分类结果里没有水体类别但影像里明明有河流。原因min_samples给得太大比如给了总像素的 5%河流像素占比不到这个数在删除步骤被丢掉了。解决把min_samples降到 0.5% 到 1% 之间。如果小地物很重要可以进一步降到 0.2%但要注意碎片簇会变多需要配合调大合并阈值。5. 进阶技巧用分裂合并日志反推地物类别数跑 ISODATA 的时候我习惯在每轮迭代打印当前簇数和触发动作这样能看出算法在哪个阶段稳定下来。下面是一个日志输出的示例% 在迭代循环里加日志 if split_flag fprintf(第 %d 轮: 分裂当前簇数 %d\n, iter, size(centers, 1)); elseif merge_flag fprintf(第 %d 轮: 合并当前簇数 %d\n, iter, size(centers, 1)); else fprintf(第 %d 轮: 无分裂无合并当前簇数 %d\n, iter, size(centers, 1)); end日志的价值在于如果簇数在某个值附近稳定了五六轮然后突然又分裂说明分裂阈值卡在临界点上可以微调一下。如果簇数一直在缓慢增长说明分裂阈值偏小。如果簇数先增后减再稳定说明算法在探索合理的类别数最终稳定的值就是比较可信的类别数。另一个技巧是用多次运行取共识。ISODATA 对初始中心有一定敏感性跑十次每次用不同的随机种子然后统计每个像素在十次运行中被分到同一类的频率。频率高的像素分类可信度高频率低的像素可能是边界区域。这个做法计算量乘以十但在数据量不大的情况下可以接受。% 多次运行取共识 n_runs 10; all_labels zeros(H, W, n_runs); for r 1:n_runs rng(r * 100); % 不同随机种子 [labels_r, ~] isodata_core(ms_data_norm, options.init_k, ... options.max_iter, options.min_samples, options.split_std, ... options.merge_dist, options.max_k); all_labels(:, :, r) labels_r; end % 计算共识度每个像素在多次运行中标签一致的频率 % 这里简化处理统计每个像素的标签众数出现频率 consensus zeros(H, W); for i 1:H for j 1:W pixel_labels squeeze(all_labels(i, j, :)); consensus(i, j) max(histcounts(pixel_labels)) / n_runs; end end figure; imagesc(consensus); colorbar; title(分类共识度图);共识度图里亮的地方表示多次运行结果一致暗的地方表示分类不稳定。暗区域通常是地物边界或者混合像素这些地方本来就不该强行分类。我一般会把共识度低于 0.6 的像素标记为“未分类”在后续应用里单独处理。最后一个习惯每次调完参数跑完把参数组合和对应的类别数、目视效果记在一个表格里。ISODATA 调参是个经验活记录多了就能看出规律。比如我自己的记录里归一化后的多光谱数据split_std在 0.15 到 0.2 之间、merge_dist在 0.3 到 0.4 之间大多数场景都能得到比较合理的结果。这个范围不是绝对的但可以作为一个起点。希望帮到你。本文还有配套的精品资源点击获取