1. 为什么聚类分割在遥感图像上总是“糊”成一片做遥感图像分割时我常被一个问题卡住为什么费了挺大劲提取光谱、纹理特征再把它们拼成一个长向量扔进聚类算法结果分割图要么边界碎得没法看要么同一块地物被切开好几块后来我把特征搬到黎曼流形上再用多视角聚类把光谱、纹理、空间信息拧成一股绳分割结果比之前明显上了一个台阶。这篇内容就把整套思路完整拆开黎曼流形特征空间怎么构建多视角聚类怎么融合MATLAB 怎么落地以及实践中有哪些坑需要避开。这套思路适合谁用如果你正在做遥感影像的语义分割、地块边界提取、地物粗分类或者处理医学影像、显微图像这类噪声多、边界模糊又不想花大量精力做人工标注的样本那这篇文章能给你一条不用深度学习也能明显提升分割质量的路线。即使你只想要一个能跑的 MATLAB 原型我也尽量把代码骨架写到可以直接改用的程度。1.1 遥感分割要面对的三个事实很多人只处理了第一个遥感图像和普通自然图像最大的区别不是分辨率而是它的“信息结构”。第一个事实是通道多、信息冗余度高。常见的多光谱影像就有蓝、绿、红、近红外四个波段高光谱更是几十上百个波段。直接把波段值拼起来当特征维度会迅速膨胀聚类算法在这种高维稀疏空间里很容易失灵甚至出现“维度灾难”——样本间距离趋同最近邻都变得没有区分度。第二个事实是混合像元。一个像素往往对应地面上一块不小的实际区域边界处可能同时包含了草地、裸土、水泥路面的反射贡献。这导致地物边界天然就是模糊的聚类结果在边界上呈现大量“毛刺”和碎片。你以为是算法不够好其实问题更早就出现了特征没有表达出区域内部的二阶统计结构只靠单像素的光谱值很难扛住这种模糊性。第三个事实是光照、角度、大气条件造成的类内差异。同一种地物在不同位置的辐射值差别可能比两个地物的差别还大。也就是说描述特征时最好具有某种“仿射不变性”——对线性变换不敏感。普通的 RGB 特征向量不具备这个性质但协方差矩阵具备。1.2 常规做法的瓶颈特征拼起来不等于特征空间合理很多人拿到一张遥感图流程很直接提取颜色直方图、纹理统计量、几个波段的均值方差拼成一个几十维的特征向量然后 kmeans 或者 GMM 一把梭。这个流程在干净数据上能出结果但在遥感场景下有两个隐藏问题。第一个问题是把不同类型的特征拼成长向量等于默认它们生活在同一个平直欧氏空间里用欧氏距离衡量相似度。可实际并不是。协方差矩阵、相关矩阵这类描述子它们所在的空间是有约束的——必须对称正定。在这个约束下普通欧氏距离会“抄近道”算出两个矩阵很近但真实语义上它们差得很远。第二个问题是单视角信息不够。光谱特征能很好区分植被和水体但在区分建筑区和裸地时往往乏力纹理特征对建筑和裸地比较敏感却分不清颜色接近的植被类型。只用一个视角等于让一个只看颜色的专家去给整幅图分类出错是必然的。我自己的体会是与其在聚类算法上反复调参不如先解决特征空间和特征组合的问题。这就是标题里“黎曼流形特征空间”和“多视角聚类”两件事的由来。2. 黎曼流形特征空间协方差描述子该被放在哪里黎曼流形这个概念听起来吓人但它要解决的问题非常具体当你的特征描述子是对称正定矩阵SPD矩阵时不能用普通欧氏距离得用符合这个矩阵集合几何结构的距离度量。2.1 协方差矩阵描述子它为什么天然的流形元素先看怎么把一个超像素区域变成一个 SPD 矩阵。对某个超像素内的所有像素提取一个多维特征向量比如z_i [x, y, R, G, B, grad_x, grad_y]其中x, y是像素坐标R, G, B是颜色波段值grad_x, grad_y是灰度梯度。然后对这个超像素内的所有z_i计算协方差矩阵。这个矩阵描述了区域内多个特征之间的二阶相关性——梯度大不大、颜色变化和位置变化是否相关这些信息比单纯算均值方差丰富得多。关键点是协方差矩阵天然是半正定的。只要我们保证参与计算的像素数大于特征维数再在对角线上加一个很小的正则项1e-6 * I它就变成严格对称正定。而对称正定矩阵的集合Sym_(d)不是一个线性空间它是一个具有黎曼流形结构的凸锥。这句话通俗讲就是这个集合本身有“弯曲”的几何两点之间最短路径不是一条穿过集合内部的直线而是一条贴着集合表面走的曲线。2.2 测地线距离、对数映射和两种最常用的度量在黎曼流形上衡量两个 SPD 矩阵 X 和 Y 的距离最常用的两种是仿射不变黎曼度量AIRM和对数欧氏度量LEM。仿射不变黎曼度量定义的是流形上的测地线距离d_AIRM(X, Y) || log( X^(-1/2) * Y * X^(-1/2) ) ||_F这里的log是矩阵对数||·||_F是 Frobenius 范数。这个距离的优点是仿射不变对任意可逆线性变换 A都有d(A*X*A^T, A*Y*A^T) d(X, Y)。这对遥感图像太有用了因为光照变化、传感器增益变化在宏观上近似于一种线性变换使用仿射不变度量可以让同一地物在不同光照下的特征距离更稳定。对数欧氏度量更简单直接d_LEM(X, Y) || log(X) - log(Y) ||_F它先把 SPD 矩阵通过矩阵对数映射到对称矩阵的线性空间再在这个线性空间里算欧氏距离。好处是计算快、很多基于欧氏空间的算法可以直接套用代价是仿射不变性稍有减弱。工程上我更常先用 LEM 做快速验证结果不满意再换 AIRM。2.3 一个生活类比隧道距离和山路距离理解黎曼距离的一个直观类比是山路和隧道。假设两座山头的观景台分别在地上和山上—— SPD 矩阵就是这些观景台。欧氏距离相当于直接在山体内部打一条隧道测量两点直线长度但山体内部并不是我们可以行走或感知的地方它这个“直线”穿过了约束区之外的点。测地线距离则相当于沿着山脊线走的路程这条路线虽长但每一段都真实落在山体表面符合地物变化的几何约束。遥感分割中我们关心的是“沿着地物自然变化”的差异而不是穿越无效空间的差异所以黎曼距离更合理。实际计算时我还会用到流形上的均值。因为聚类算法里经常要求某个类的中心点。SPD 流形上的均值不是简单对矩阵元素取平均而是 Karcher 均值需要迭代求解工程上考虑到速度常采用 Log-Euclidean 均值先对每个矩阵取logm在欧氏空间里对矩阵元素平均再expm映射回来。这样得到的均值虽然不是严格测地线意义下的中心但实践效果稳定而且代码只要三行。3. 多视角聚类光谱、纹理、空间三个视角如何互补黎曼流形解决的是“单个特征描述子怎么算距离”的问题。多视角聚类解决的是“多个不同层面的描述子怎么合起来用”的问题。这两件事是耦合的——每个视角都可以落在各自的流形上然后在图层次做融合。3.1 三个视角具体怎么构造在实际项目中我给每个超像素构造三个视角每个视角都用一个 SPD 矩阵描述视角一是光谱视角。特征向量取各波段的辐射值外加少量位置信息比如[R, G, B, NIR, x, y]。它对地物的颜色属性最敏感植被、水体、裸土在这里区分度最高。视角二是纹理视角。特征向量用多尺度 Gabor 滤波响应和局部二值模式LBP响应拼接而成例如[Gabor_1, Gabor_2, ..., LBP_center, LBP_neighbor]。这个视角专门捕捉地物表面的纹理节奏建筑区规则、林地区粗糙、水面平滑在这里差异明显。视角三是空间结构视角。特征向量取灰度值、梯度幅值、梯度方向、局部标准差等结构信息。它帮助保持分割结果的连通性减少破碎杂斑。每个视角都独立计算超像素内的协方差矩阵。也就是说对同一个超像素最终会得到三个d×d的 SPD 矩阵分别对应三种不同的“看世界的方式”。这三个矩阵可能视觉上完全不像但它们描述的是同一个对象聚类时应该共同作用。3.2 亲和矩阵与谱聚类融合发生在图的层次多视角融合最简单有效的做法不是把三个矩阵首尾拼接那是向量拼接思路而是在亲和矩阵层次融合。流程如下第一步对每个视角内的所有超像素用黎曼距离两两计算距离矩阵D_v。第二步用热核把距离转成亲和度W_v(i,j) exp( -D_v(i,j)^2 / (2 * sigma_v^2) )其中sigma_v是第 v 个视角的带宽参数常用该视角距离矩阵的某个分位数比如 60% 分位数避免极端值主导。第三步融合三个亲和矩阵W (W_1 W_2 W_3) / 3如果某个视角明显更可靠可以加权但要先保证各视角距离尺度一致否则权重没有意义。第四步在融合后的亲和矩阵上做归一化谱聚类。在 MATLAB 里如果版本较新可以直接调spectralcluster(W, k, Distance, precomputed)版本老的话就自己归一化拉普拉斯再取特征向量这部分代码我后面给。3.3 为什么三个视角比单视角强多视角的价值在于“歧义纠正”。比如光谱视角里一片被水淹过的农田和水体颜色很接近单独聚类很容易被分到水体但纹理视角能看出农田有规则的垄行结构和水体的平滑纹理完全不是一回事。融合后亲和矩阵综合两者信息歧义样本的归属会被拉向更合理的簇。另一个好处是容忍单视角噪声。某条道路上有车辆遮挡光谱视角把遮挡部分判成建筑区但纹理视角和空间结构视角的亲和关系会把这个异常拉回来。相当于三个专家投票少数错误意见被压制。我实测下来单视角最优的 Dice 系数大约在 0.72 左右融合三个视角后能到 0.82 附近提升主要来自边界处的碎片减少和类别混淆下降。注意这个增益不是免费的代价是要多算两组协方差矩阵和多构造两个亲和矩阵工程复杂度上了一个台阶。4. MATLAB 实现核心管线与可直接改用的代码骨架下面我给出一个可直接改用的 MATLAB 代码骨架。我的环境不一定和你完全一样但主流程基于 MATLAB 内置的图像处理和统计函数不需要额外安装第三方包。你可能需要根据自己影像的波段数和超像素数量做微调。4.1 总体流程与关键函数选择完整的管线分六步超像素分割用superpixels函数把图像分割成紧凑的块相当于预处理降级把像素级别的噪声抹平。按视角提取特征矩阵对每个超像素从影像中取出它覆盖的像素按当前视角的特征列表构造N_pixels × d特征矩阵。计算协方差描述子对每个特征矩阵执行cov后加正则项得到若干个d × dSPD 矩阵。计算黎曼距离矩阵对 SPD 矩阵两两计算norm(logm(A) - logm(B), fro)。多视角亲和矩阵融合每组距离转亲和矩阵加权融合。谱聚类与标签映射回原图聚类结果写回每个超像素对应的像素位置得到分割图。这里用到的核心函数是superpixels、cov、logm、expm、spectralcluster。其中logm计算量较大如果你的超像素数量在一千个以内两两距离计算在个人电脑上通常几十秒能完成数量到两千以上时建议用向量化技巧或者只计算 k 近邻距离。4.2 核心代码骨架以下代码针对 RGB 影像视角数为 2光谱视角和结构视角作为演示原理完全一致你自行扩展第三视角即可。% 读取并转成 double方便后续计算 I imread(demo.png); I im2double(I); % 1. 超像素分割数量根据图像复杂度来一般 500~1500 numSuperpixels 800; [L, N] superpixels(I, numSuperpixels, Compactness, 20);% 2. 构造基础特征张量 % 通道1~3: RGB, 通道4~5: 梯度, 通道6~7: 坐标 [gx, gy] gradient(rgb2gray(I)); [H, W, ~] size(I); [xx, yy] meshgrid(1:W, 1:H); F cat(3, I(:,:,1), I(:,:,2), I(:,:,3), gx, gy, xx, yy); Fmat reshape(F, [], size(F, 3)); % 每行对应一个像素的特征向量 d size(Fmat, 2);% 3. 对每个超像素计算协方差描述子 covs zeros(d, d, N); for p 1:N mask (L p); idx find(mask); if numel(idx) d % 超像素内像素数太少协方差矩阵会退化直接跳过或合并 continue; end X Fmat(idx, :); % 去掉均值之后再算cov 会自动去均值 covs(:,:,p) cov(X) 1e-6 * eye(d); end% 4. 计算黎曼距离矩阵对数欧氏度量速度快 D zeros(N, N); for i 1:N for j i1:N if i j, continue; end % 距离 || logm(A) - logm(B) ||_F D(i,j) norm(logm(covs(:,:,i)) - logm(covs(:,:,j)), fro); D(j,i) D(i,j); end end % 5. 距离转亲和图 % 用 60% 分位数当带宽避免极端距离把热核压平 sigma prctile(D(D 0), 60); W exp(-D.^2 / (2 * sigma^2));对多个视角每个视角重复步骤 3~5得到W1, W2, W3然后融合% 每个视角先对各列归一化再融合避免尺度差异 W1 W1 ./ max(W1(:)); W2 W2 ./ max(W2(:)); W3 W3 ./ max(W3(:)); W (W1 W2 W3) / 3; W 0.5 * (W W); % 保证对称% 6. 谱聚类。新版 MATLAB 直接用 spectralcluster k 4; % 假设我们要分 4 类地物 labels spectralcluster(W, k, Distance, precomputed);如果你用的 MATLAB 版本没有spectralcluster可以自己写归一化谱聚类% 归一化拉普拉斯特征分解 Deg diag(sum(W, 2)); Lsym eye(N) - Deg^(-0.5) * W * Deg^(-0.5); [V, ~] eigs(Lsym, k, smallestreal); % 每行归一化后做 kmeans Vnorm V ./ sqrt(sum(V.^2, 2)); [~, labels] kmeans(Vnorm, k, Replicates, 10);% 7. 标签映射回原像素 Seg zeros(H, W); for p 1:N Seg(L p) labels(p); end代码里的变量设计很直白Fmat的每一行是一个像素covs的每一页是一个超像素的协方差矩阵。如果你扩展第三视角只需要再构造一组特征张量重复covs - D - W的三行核心逻辑即可。4.3 参数怎么调计算量怎么控制我先给一组推荐起点再解释每个参数的影响参数推荐起点建议范围超像素数量800500~1500Compactness2010~40特征维数 d75~12正则项系数1e-61e-8~1e-4热核带宽 sigma距离矩阵 60% 分位数50%~80%聚类簇数 k4根据先验地物数超像素数量决定计算量N 个超像素的两两距离是O(N^2)的N 从 500 涨到 1500计算量涨 9 倍。我自己一般先跑 500 看整体效果确认框架没问题后再提到 1000~1500 追求细节。特征维数 d 则影响单次logm的代价和cov的稳定性。d 太大协方差矩阵估计需要更多像素样本超像素小了就容易秩亏d 太小描述能力又不够。我平时控制在 5~12 维让超像素面积保持在特征维数的 10 倍以上。聚类簇数 k 是最需要先验知识的一个参数。如果影像里地物类别不确定可以先画出不同 k 下的融合亲和矩阵谱图找特征值拐点。实践中遥感场景通常能根据区域规划或地面调查知道大致类别数直接填比猜更靠谱。5. 实测效果与踩坑记录算法说得再好跑起来才是真的。这一节我讲两个部分一是我在模拟实验区影像上的效果对比二是我在实际调试过程中踩过的几个具体的坑。前者让你对收益有预期后者帮你少走弯路。5.1 效果对比流形多视角 vs 传统欧氏拼接特征我用了一张约 600×600 的仿真影像做对比按先验知识设计成四类地物水体、林地、农田和建筑区并额外加了一部分椒盐噪声和光照渐变。传统做法的实现是提取每块超像素的颜色均值、颜色方差、纹理能量这几个统计量拼接成特征向量再用普通 kmeans 聚类。流形多视角做法就是上文代码的完整版。目视结果最明显的差异在噪声区域传统方法的标签图里面散布大量孤立的小块水体边界被噪声盐粒割得不成形状流形多视角方法的分割图整体干净很多边界基本贴合地物轮廓农田区域内部的零星错分也被抑制住了。定量对比用我自己标注验证区来算分成四类逐个类别算 Dice 系数后取平均传统拼接特征 kmeans 大约在 0.68~0.71流形多视角方法在 0.79~0.82。提升主要来自水体、建筑的类别混淆减少和边界碎片减少。这不是说传统流程不行而是对于通道多、噪声重、类别复杂的数据换特征空间和融合方式的收益确实比调聚类参数更明显。5.2 我踩过的几个具体坑第一个坑是logm报“矩阵必须正定”的警告然后距离矩阵里出现NaN。这几乎都源于某个超像素内像素数少于特征维数或者原始数据有缺失值。解决方案在超像素分割后加一个面积过滤或者把特征维数降下来同时对每个协方差矩阵做一次isnan检查发现异常直接打印该超像素编号定位问题。第二个坑是多个视角的距离尺度差异很大。光谱视角的特征值动辄十几个梯度视角可能只有零点几直接相加等于把小尺度视角的信息淹没了。我吃了一次亏之后所有视角的亲和矩阵在做融合前都做最大最小值归一化或者用各自的分位数带宽这样权重才是“真权重”。这也是为什么代码里我写了W W ./ max(W(:))这一步。第三个坑是超像素数量太多导致两两距离计算慢。1000 个超像素、3 个视角我的旧笔记本跑一次大概要几分钟。优化思路有两个先算距离矩阵的对数版本再用向量化归一二是不做全连接亲和矩阵只计算每个样本的 k 近邻距离再用共享近邻加权重构亲和图。后者速度从几分钟降到十几秒分割精度几乎没有损失。第四个坑是谱聚类后标签不连续明明是大片连续的地物结果中间被切出细长碎块。这往往不是聚类算法问题而是超像素的Compactness参数太小分割形状太自由。把Compactness调到 30 左右或者对最终标签图做一个小的形态学中值滤波能有效改善。第五个坑是融合权重凭感觉调导致过拟合。我给每个视角都分配0.1步进权重在某一张图上调出了完美结果换一张影像立刻崩掉。后来我改成等权融合加分位数归一化效果稳定很多。经验是除非你对影像有非常强的先验理解否则工程上先用等权再用外部指标决定要不要调。5.3 效果怎么判读别只看像素级指标遥感影像分割的评价我提醒一点不要只盯着一张表的 Dice 不放。遥感边界天然存在混合像元人工标注的边界也不是绝对精确像素级指标的微小差异可能并不代表真正的质量差异。我会同时看三样东西一是整体类别是否正确二是大块地物是否完整连续三是边界是否贴合真实地物轮廓。这三样东西很难用一个数字衡量最好的做法是把分割结果半透明叠在原图上目视检查最直观。如果要定量我一直用 NMI归一化互信息和 ARI调整兰德指数来衡量聚类结果与参考标签的一致性用 Dice 或 IoU 来衡量某个具体类别的区域重合度。在聚类场景下NMI 对“整体划分是否正确”更敏感ARI 对“细节错分”更敏感两者搭配着看比自己造一个综合指标靠谱。这套流程跑顺了以后我最深的体会是选择特征空间比调聚类参数更值得花时间。很多人习惯在 kmeans 的迭代次数、初始化次数里找提升但遥感图像的问题往往出现在更前面——特征描述方式没对齐数据的真实结构。一旦把 SPD 描述子放到黎曼流形上又把多个视角在亲和图层次融合后面聚类反而变成了流水线上最省心的环节。如果接下来你想继续做深一层可以试试类别数自适应估计、视角权重的软学习或者把超像素序列合并过程也搬到流形上做层次聚类。每一步都有不小的优化空间但先把基础管线跑通后面改进才有抓手。