简介：这份资源提供ISODATA聚类算法的MATLAB实现，面向从事遥感图像处理、模式识别与无监督分类的学习者和研究人员。ISODATA即迭代自组织数据分析技术，是对k-means聚类算法的改进，能够克服k-means需预先指定类别数的局限，通过分裂与合并机制自动调整聚类数目，在多光谱图像的无监督分类任务中表现尤为有效。资源包共6个文件，以cpp与h源码文件为核心，另含txt说明、md文档及gitignore、gitattributes等工程配置文件，整体约6KB，结构精简，便于直接阅读算法实现逻辑或移植到MATLAB环境中调试运行。目前已有485人学习下载，适合希望理解ISODATA聚类原理、对比k-means改进思路，或需要将其应用于多光谱遥感图像分类实验的读者参考，可帮助快速掌握算法核心流程与参数调节方法。

1. 多光谱图像分类为什么先跑 ISODATA 而不是 K-means

拿到一幅多光谱影像，第一反应往往是做无监督分类，看看地物到底能分成几类。很多人顺手就上 K-means，结果跑完发现类别数得自己拍脑袋定，定少了地物混在一起，定多了同一片地被切得七零八落。ISODATA（Iterative Self-Organizing Data Analysis Technique）解决的正是这个痛点：它允许在迭代过程中自动合并过于接近的类、拆散方差过大的类，最终类别数由数据自己"长"出来，而不是你事先硬塞一个 K 值。

这篇笔记面向的是手里有 ENVI 或 MATLAB 环境、想用 ISODATA 对多光谱数据做无监督分类的从业者。我会把 ISODATA 的 MATLAB 实现从原理到代码逐层拆开，包括聚类中心初始化、分裂与合并的判定条件、迭代终止逻辑，以及多光谱波段数据怎么组织成算法能吃的矩阵。读完你能自己写出一份可跑的 ISODATA 脚本，知道每个参数改了会发生什么，也能判断什么场景下 ISODATA 比 K-means 更值得用。

2. ISODATA 的算法骨架：分裂、合并与迭代终止

2.1 从 K-means 到 ISODATA 多了哪三个机制

K-means 的流程很干净：定 K 个中心，把每个样本分给最近的中心，重算中心，循环到中心不动。ISODATA 在这套骨架上加了三件事，让它能自适应调整类别数。

第一是分裂。当一个类的标准差在某个维度上超过阈值，说明这个类内部散得太开，可能混了两种地物，就把它沿方差最大的方向劈成两个。第二是合并。当两个类的中心距离小于阈值，说明它们本质上是同一类，就合并成一个。第三是每次迭代都重新判定分裂与合并，而不是只在开头定一次。这三个机制叠加，导致 ISODATA 的类别数在迭代中是动态变化的，可能先增后减，最后收敛到一个稳定值。

理解这一点很关键：ISODATA 不是"先定类数再聚类"，而是"边聚类边调类数"。所以它的输出类别数不可预测，但通常比拍脑袋定的 K 更贴合数据分布。

2.2 多光谱数据怎么组织成算法输入矩阵

多光谱影像在 MATLAB 里通常是一个三维数组，行×列×波段。ISODATA 处理的是样本向量，所以第一步是把三维数组展平成二维矩阵：每一行是一个像素，每一列是一个波段。

% 假设 img 是 H×W×B 的多光谱影像，数据类型为 double [H, W, B] = size(img); % 展平为 (H*W)×B 的样本矩阵，每行一个像素 data = reshape(img, H*W, B); % 去掉无效像素（如全零或 NaN 的背景） valid_mask = all(~isnan(data), 2) & any(data ~= 0, 2); data = data(valid_mask, :); fprintf('有效像素数：%d，波段数：%d\n', size(data,1), size(data,2));

这里有几个细节值得说。reshape的顺序是按列优先，MATLAB 里三维数组展平后像素顺序和img(:,:,1)的列扫描一致，这对后续把分类结果还原成图像很重要。valid_mask那行同时排除了 NaN 和全零像素，因为多光谱影像边缘常有填充值，不排除会污染聚类中心。如果你的数据已经做过辐射定标和大气校正，数值范围差异大，建议在展平后做一次归一化，否则数值大的波段会主导距离计算。

2.3 距离度量与聚类中心的初始化策略

ISODATA 用欧氏距离衡量像素与中心的接近程度。在多光谱场景下，如果各波段量纲差异大，欧氏距离会被大量纲波段绑架。常见做法是逐波段标准化：

% 逐波段标准化，使每个波段均值为0、标准差为1 data_norm = (data - mean(data)) ./ std(data);

初始化聚类中心有两种常见策略。一种是随机选 K 个像素作为初始中心，简单但结果不稳定，跑两次可能不一样。另一种是沿数据主轴等间隔取点，稳定性更好。我一般用第一种加固定随机种子，方便复现：

rng(42); % 固定种子，保证每次运行结果一致 K_init = 8; % 初始类别数，ISODATA 会在此基础上调整 idx = randperm(size(data_norm,1), K_init); centers = data_norm(idx, :);

初始 K 不用太精确，ISODATA 的分裂合并会修正它。但也不能太离谱，比如真实地物有 5 类你初始给 50，迭代前期会疯狂合并，收敛慢。经验值是真实类别数的 1.5 到 2 倍。

3. 用 MATLAB 实现 ISODATA 的核心迭代循环

3.1 分裂判定：标准差阈值与分裂方向

分裂的核心逻辑是：对每个类，计算它在各波段上的标准差，如果最大标准差超过sigma_threshold，并且当前类别数没超过上限，就分裂。

function [centers, labels] = isodata_iterate(data, centers, params) % params 包含：sigma_thr, d_merge, K_max, K_min, max_iter for iter = 1:params.max_iter % 步骤1：分配每个样本到最近的中心 dists = pdist2(data, centers); [~, labels] = min(dists, [], 2); % 步骤2：重算每个类的中心 K = size(centers, 1); new_centers = zeros(K, size(data,2)); for k = 1:K member = data(labels == k, :); if ~isempty(member) new_centers(k, :) = mean(member, 1); else new_centers(k, :) = centers(k, :); % 空类保持原中心 end end centers = new_centers; % 步骤3：分裂判定 if K < params.K_max for k = 1:K member = data(labels == k, :); if size(member,1) < 2, continue; end sigma = std(member, 0, 1); [max_sigma, max_dim] = max(sigma); if max_sigma > params.sigma_thr % 沿方差最大维度分裂 delta = zeros(1, size(data,2)); delta(max_dim) = 0.5 * max_sigma; centers = [centers; centers(k,:) + delta]; centers(k,:) = centers(k,:) - delta; K = K + 1; end end end % 步骤4：合并判定 if K > params.K_min dist_matrix = pdist2(centers, centers); dist_matrix(logical(eye(K))) = inf; % 排除自身 [min_dist, merge_idx] = min(dist_matrix(:)); if min_dist < params.d_merge [i, j] = ind2sub([K, K], merge_idx); % 合并 i 和 j，新中心取加权平均 ni = sum(labels == i); nj = sum(labels == j); merged = (ni * centers(i,:) + nj * centers(j,:)) / (ni + nj); centers(i,:) = merged; centers(j,:) = []; K = K - 1; end end % 步骤5：检查中心是否稳定 if iter > 1 && max(abs(centers(:) - prev_centers(:))) < 1e-4 break; end prev_centers = centers; end end

这段代码把 ISODATA 的五个步骤串在一个循环里。pdist2计算所有样本到所有中心的距离矩阵，min取每行最小值得到标签。分裂时沿标准差最大的维度偏移半个标准差，这是经典做法，偏移量太大会导致新中心跑到别的类地盘，太小则分裂后很快又被合并回去。合并时用加权平均而不是简单平均，因为两个类样本数可能差很多，简单平均会让小类被大类吞掉。

3.2 合并判定：类间距离阈值怎么设

d_merge是最难调的参数。设太小，该合并的不合并，类别数偏多；设太大，不该合并的也合并了，地物混在一起。一个实用的起点是取所有类间距离的中位数的 0.5 到 0.7 倍：

% 自动估算 d_merge 的初始值 dists = pdist(centers); d_merge = 0.6 * median(dists);

但注意，这个估算只在第一次迭代前做一次，后续迭代中类中心在变，固定阈值可能不再合适。更稳的做法是每次迭代重新估算，或者用相对阈值——比如类间距离小于所有类间距离均值的某个比例时才合并。

3.3 迭代终止条件与结果还原成分类图

终止条件有三个：达到最大迭代次数、中心变化小于阈值、或者类别数在上下限之间且中心稳定。实际跑的时候，我一般设max_iter = 50，中心变化阈值1e-4，大多数数据在 20 到 30 次迭代内收敛。

跑完之后要把标签还原成图像：

% 把标签填回有效像素位置，再还原成 H×W full_labels = zeros(H*W, 1); full_labels(valid_mask) = labels; class_map = reshape(full_labels, H, W); % 可视化 figure; imagesc(class_map); colormap(jet(K)); colorbar; title('ISODATA 分类结果');

valid_mask在这里再次派上用场，因为展平时去掉了无效像素，还原时必须按同样的索引填回去，否则分类图会错位。这是新手最容易翻车的地方之一，跑完发现图像是斜的或者有条纹，八成是索引对不上。

4. 参数调优与多光谱场景的避坑记录

4.1 初始类别数、分裂阈值、合并阈值的联动关系

这三个参数不是独立的。初始 K 大，分裂阈值可以设小一点，因为类已经够多；初始 K 小，分裂阈值要设小才能让类长出来。合并阈值和分裂阈值是反向的：分裂阈值小，类容易分裂，类别数偏多，这时候合并阈值要相应放大来平衡。

我一般按这个顺序调：先固定合并阈值，调分裂阈值让类别数接近预期；再固定分裂阈值，调合并阈值消除明显重复的类；最后微调初始 K。整个过程跑三五次就能找到稳定组合。

4.2 波段选择与降维对分类结果的影响

多光谱数据波段多，不是所有波段都对分类有用。比如可见光波段之间高度相关，全放进去会让距离计算被冗余信息主导。常见做法是先做主成分分析，取前几个主成分跑 ISODATA，再把标签映射回原空间。

% PCA 降维后再聚类 [coeff, score, ~, ~, explained] = pca(data_norm); n_components = find(cumsum(explained) > 95, 1); % 保留95%方差 data_pca = score(:, 1:n_components);

但要注意，PCA 后的主成分是原始波段的线性组合，物理意义不直观。如果后续要做地物识别，最好保留原始波段做验证。我的习惯是 PCA 跑一遍看类别数，再用原始波段跑一遍对比，两者类别数接近就说明波段冗余不严重。

4.3 避坑记录：ISODATA 跑出来的五类典型翻车

现象一：类别数一直增长到上限。原因通常是分裂阈值设得太小，或者数据没标准化导致某些波段方差天然很大。解决方法是先做逐波段标准化，再把sigma_thr调大 0.1 到 0.2。

现象二：类别数收敛到 1。合并阈值设太大了，所有类都被合并。检查d_merge是否超过了类间距离的最大值，如果是，直接减半。

现象三：分类图有条纹或斜纹。展平和还原的索引不一致。检查valid_mask在展平和还原时是否用了同一个逻辑数组，以及reshape的维度顺序是否匹配。

现象四：同一类地物被分成多块不相邻区域。这是无监督分类的固有特性，ISODATA 只看光谱不看空间。如果要做空间一致性，需要在聚类后加一步空间滤波或马尔可夫随机场后处理。

现象五：跑两次结果不一样。初始化用了随机选点且没固定种子。加rng(42)或者改用等间隔初始化。如果数据量特别大，pdist2可能内存溢出，这时候要分块计算距离。

注意：ISODATA 的类别数是数据驱动的，不要强行让它输出你想要的类数。如果跑出来的类数和预期差很多，先检查数据预处理，再调参数，最后才怀疑算法本身。

5. 从分类结果到可用产品的最后一步

跑出分类图只是开始，真正要落地还得做几件事。第一是给每个类赋予地物标签，这需要结合地面真值或者目视解译，ISODATA 只告诉你"这些像素是一伙的"，不告诉你"这伙是水体还是植被"。第二是做精度评估，常用混淆矩阵和 Kappa 系数，但无监督分类没有训练样本，通常只能做定性评估或者和已有分类产品对比。

一个实用的技巧是用 ISODATA 的类别数作为 K-means 的 K 值，再跑一遍 K-means。因为 ISODATA 收敛后的类别数已经反映了数据的自然分组数，用它去初始化 K-means，既保留了 ISODATA 的自适应能力，又得到了 K-means 的稳定输出。我做过几组对比，这种"ISODATA 定 K + K-means 精聚类"的组合，在类别数合理性和运行稳定性上都比单独用其中一个好。

% 用 ISODATA 收敛的类别数作为 K-means 的 K K_final = size(centers, 1); [labels_km, centers_km] = kmeans(data_norm, K_final, 'Replicates', 5);

Replicates设 5 表示跑 5 次取最优，避免局部最优。这一步的计算量比 ISODATA 小，因为 K 已经定了。

最后说个我自己的习惯：每次跑完 ISODATA，先把类别数和每类的像素数打出来看一眼。如果某一类只有几十个像素，大概率是噪声或者边缘效应，可以直接合并到最近的类。如果某一类占了 80% 以上，说明分裂不够，得回去调sigma_thr。这个检查花不了几秒钟，但能省掉很多返工。希望帮到你。

本文还有配套的精品资源，点击获取