问题背景 #

我有一个带短时间序列的2D复数图像，希望对相似像素进行聚类/分割。图像包含大致静态基底，叠加了随时间变化（主要是复数角度）的斑点区域，这些斑点在图像模值中也略有体现。

前期尝试 #

• 首次尝试k-means，但仅依据均值聚类（时序和角度信息区分度更高），效果不佳；
• 第二次尝试ICA，选取幅值最大的k个分量，虽识别出部分差异区域，但未找到关注的小目标区域（视觉易识别但区域小）。

当前现状 #

因前两种方法效果不佳，调研认为谱聚类适用，但受限于内存无法直接处理，故考虑对数据分片后分别进行谱聚类。有资料建议先分片聚类再合并，但未给出具体方法，相关论文并非分片处理全部数据，而是剔除远离主成分的向量。

针对具体问题的解答 #

问题1：如何合并各分片的聚类结果？ #

分片聚类后，每个分片的簇编号是独立的，核心思路是先对齐簇的本质特征，再映射全局编号，具体步骤如下：

1. 提取簇的代表性特征：针对你场景中目标区域的时序角度变化核心特征，对每个分片的每个簇，计算以下统计量：
   • 时序角度的均值/标准差
   • 时序模值的均值/标准差
   • 角度时序的自相关系数（体现时间相关性）
     参考Matlab代码片段：

   % 假设cluster_labels_cell是存储各分片簇标签的cell数组
   global_cluster_centers = [];
   for i = 1:length(cluster_labels_cell)
       slab_data = complexrowsTranspose((i-1)*npix+1:min(i*npix, size(complexrowsTranspose,1)), :);
       labels = cluster_labels_cell{i};
       unique_labels = unique(labels);
       for lbl = unique_labels
           idx = labels == lbl;
           % 拆解复数时序数据，计算特征
           angle_seq = angle(reshape(slab_data(idx,:)', [15, sum(idx)]));
           norm_seq = abs(reshape(slab_data(idx,:)', [15, sum(idx)]));
           feat = [mean(angle_seq,1)', std(angle_seq,1)', mean(norm_seq,1)', std(norm_seq,1)'];
           global_cluster_centers = [global_cluster_centers; mean(feat,1)];
       end
   end
   

2. 全局簇对齐：将所有分片的簇中心用k-means聚成你需要的全局K个簇，得到每个分片簇对应的全局簇编号映射表。
3. 像素映射：遍历每个分片的像素，根据其所在分片簇的映射关系，替换为全局簇编号，最终拼接成完整图像的聚类结果。

这种方法完全避开了分片间编号混乱的问题，用簇的本质特征做对齐，更贴合你的场景需求。

问题2：如何处理跨分片的像素亲和性？ #

分片后直接构建全局亲和矩阵确实会内存爆炸，这里有两种落地性强的思路：

1. 锚点式全局亲和性构建：

   • 从每个分片中按簇均匀采样少量代表性像素（比如每个分片选100个）作为全局锚点；
   • 计算所有锚点之间的亲和性（用高斯核：exp(-||x-y||²/(2σ²))），构建小尺寸的全局锚点亲和矩阵和拉普拉斯矩阵；
   • 对每个分片的像素，计算其到所有锚点的亲和性，将像素特征映射到锚点空间，再基于锚点的聚类结果推导像素的簇标签。
     这种方法的内存占用仅由锚点数量决定，能大幅降低压力。

2. 特征向量合并法：

   • 你提到的合并所有分片的谱聚类特征向量是可行的，但要先对每个分片的特征向量做L2归一化，避免分片间的特征尺度差异；
   • 将所有分片的特征向量拼接成一个大矩阵，再做一次全局k-means，得到最终聚类结果。
     对于非对称的拉普拉斯矩阵，建议转换为对称形式：L_sym = D^(-1/2) * (A + A')/2 * D^(-1/2)（D为度矩阵），这样后续特征分解的结果会更稳定。

问题3：内存占用更低的替代方法 #

结合你的小目标、时序特征主导、内存受限的场景，推荐以下几种方法：

1. 稀疏谱聚类：
   不要构建全连接的亲和矩阵，只保留每个像素的K个最近邻（比如K=20），构建稀疏亲和矩阵。Matlab中可以用knnsearch实现，再用稀疏矩阵专用的eigs函数计算特征向量，内存占用会大幅降低：

   % 构建稀疏亲和矩阵
   [idx, dist] = knnsearch(complexrowsTranspose, complexrowsTranspose, 'K', 21); % 包含自身
   idx = idx(:,2:end); % 去掉自身
   dist = dist(:,2:end);
   sigma = median(dist(:));
   A = sparse(repmat(1:size(complexrowsTranspose,1),1,20)', idx(:), exp(-dist(:).^2/(2*sigma^2)), size(complexrowsTranspose,1), size(complexrowsTranspose,1));
   A = (A + A')/2; % 转对称
   D = sparse(1:size(A,1), 1:size(A,1), sum(A,2));
   L = D - A;
   [V, ~] = eigs(L, K, 'smallestabs'); % 取最小的K个特征值对应的特征向量
   

2. 基于密度的聚类（DBSCAN/HDBSCAN）：
   这类方法不需要构建全局矩阵，只计算每个点的局部邻居，非常适合你的小目标场景（小目标区域像素特征相似、密度高）。Matlab中可以直接用dbscan函数（需Statistics and Machine Learning Toolbox），特征选用时序角度/模值的统计量即可。
3. 先降维再聚类：
   先对每个像素的时序数据做降维，比如用PCA把30维（15维实部+15维虚部）降到5-10维（保留95%以上方差）。降维后数据维度大幅降低，后续无论是谱聚类还是k-means，内存和计算时间都会显著减少，而且PCA能很好捕捉到时序变化的核心模式。

补充分片聚类+合并代码示例 #

基于你提供的代码框架，补充完整的分片谱聚类与结果合并流程：

%% 分片谱聚类
cluster_labels_cell = cell(nslabpix + (nrestpix>0), 1);
K_cluster = 3; % 假设每分片分3类
for i = 1:nslabpix
    slab_data = complexrowsTranspose((i-1)*npix+1:i*npix, :);
    % 构建稀疏亲和矩阵
    [idx, dist] = knnsearch(slab_data, slab_data, 'K', 21);
    idx = idx(:,2:end); dist = dist(:,2:end);
    sigma = median(dist(:));
    A_slab = sparse(repmat(1:npix,1,20)', idx(:), exp(-dist(:).^2/(2*sigma^2)), npix, npix);
    A_slab = (A_slab + A_slab')/2;
    D_slab = sparse(1:npix,1:npix, sum(A_slab,2));
    L_slab = D_slab - A_slab;
    % 提取特征向量
    [V_slab, ~] = eigs(L_slab, 5, 'smallestabs');
    % k-means聚类
    [labels, ~] = kmeans(V_slab, K_cluster);
    cluster_labels_cell{i} = labels;
end
% 处理剩余像素
if nrestpix > 0
    slab_data = complexrowsTranspose(nslabpix*npix+1:end, :);
    [idx, dist] = knnsearch(slab_data, slab_data, 'K', 21);
    idx = idx(:,2:end); dist = dist(:,2:end);
    sigma = median(dist(:));
    A_slab = sparse(repmat(1:nrestpix,1,20)', idx(:), exp(-dist(:).^2/(2*sigma^2)), nrestpix, nrestpix);
    A_slab = (A_slab + A_slab')/2;
    D_slab = sparse(1:nrestpix,1:nrestpix, sum(A_slab,2));
    L_slab = D_slab - A_slab;
    [V_slab, ~] = eigs(L_slab, 5, 'smallestabs');
    [labels, ~] = kmeans(V_slab, K_cluster);
    cluster_labels_cell{end} = labels;
end

%% 合并聚类结果
% 提取所有簇中心特征
global_feats = [];
for i = 1:length(cluster_labels_cell)
    slab_data = complexrowsTranspose((i-1)*npix+1:min(i*npix, size(complexrowsTranspose,1)), :);
    labels = cluster_labels_cell{i};
    unique_lbls = unique(labels);
    for j = 1:length(unique_lbls)
        idx = labels == unique_lbls(j);
        angle_seq = angle(reshape(slab_data(idx,:)', [15, sum(idx)]));
        norm_seq = abs(reshape(slab_data(idx,:)', [15, sum(idx)]));
        feat = mean([mean(angle_seq,1), std(angle_seq,1), mean(norm_seq,1), std(norm_seq,1)],2);
        global_feats = [global_feats; feat'];
    end
end
% 全局簇对齐
[global_cluster_map, ~] = kmeans(global_feats, K_cluster);
% 映射每个分片的簇到全局编号
global_labels = zeros(size(complexrowsTranspose,1),1);
ptr = 1;
for i = 1:length(cluster_labels_cell)
    slab_labels = cluster_labels_cell{i};
    unique_lbls = unique(slab_labels);
    slab_map = zeros(max(unique_lbls),1);
    for j = 1:length(unique_lbls)
        feat_idx = (i-1)*length(unique_lbls) + j;
        slab_map(unique_lbls(j)) = global_cluster_map(feat_idx);
    end
    global_labels(ptr:ptr+length(slab_labels)-1) = slab_map(slab_labels);
    ptr = ptr + length(slab_labels);
end
% 转成图像格式
cluster_image = reshape(global_labels, [921,921]);
figure; imshow(cluster_image, []); title('全局聚类结果');

内容的提问来源于stack exchange，提问作者Leo