基于顶点成分分析的多光谱图像解混方法-CSDN博客

在原有高光谱混合像元分解代码基础上，增加端元提取过程的完整MATLAB实现。该方案采用无监督端元提取算法（VCA），无需已知端元光谱，直接从混合像元邻域数据中自动提取端元，并结合全约束最小二乘法（FCLS）进行丰度反演。

代码：

%% 步骤1: 模拟数据生成（含邻域数据）
numBands = 8; % 光谱波段数
numEndmembers = 3; % 真实端元数量

% 生成三个端元光谱（行向量）
endmembers_true = [
    0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1;  % 植被
    0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8;  % 水体
    0.4, 0.3, 0.2, 0.1, 0.2, 0.3, 0.4, 0.5   % 裸土
];

% 生成混合像元（真实丰度：植被30%，水体50%，裸土20%）
abundances_true = [0.3; 0.5; 0.2];
mixed_pixel = endmembers_true' * abundances_true;

% 模拟邻域数据（100个像元，以混合像元为中心添加噪声）
rng(0);
neighborhood = repmat(mixed_pixel', 100, 1) + randn(100, numBands) * 0.05;

%% 步骤2: 端元提取（VCA算法）
% 估计端元数量（PCA累计贡献率>95%）
[coeff, score, latent] = pca(neighborhood);
explained = cumsum(latent) / sum(latent);
m_est = find(explained > 0.95, 1); % 自动估计端元数量

% 使用VCA提取端元
[endmembers_est] = myVCA(neighborhood, m_est);
endmembers_est = endmembers_est'; % 转置为 m_est × numBands

%% 步骤3: 丰度反演（FCLS）
% 构建端元矩阵（注意维度：波段×端元数）
M_est = endmembers_est'; 

% 全约束最小二乘求解（非负+归一化）
x = lsqnonneg(M_est, mixed_pixel);
abundance_est = x / sum(x);

% 重建光谱
reconstruction = M_est * abundance_est;
RMSE = sqrt(mean((mixed_pixel - reconstruction).^2));

%% 步骤4: 结果可视化
figure('Position', [100, 100, 1200, 500]);

% 1. 丰度对比图
subplot(1, 3, 1);
bar([abundances_true; abundance_est]', 'grouped');
legend({'真实丰度', '估计丰度'}, 'Location', 'northoutside');
title('丰度对比');
set(gca, 'XTickLabel', {'植被', '水体', '裸土'});
ylabel('比例');
ylim([0, 1]);

% 2. 端元光谱对比
subplot(1, 3, 2);
hold on;
colors = ['r', 'g', 'b'];
for i = 1:numEndmembers
    plot(1:numBands, endmembers_true(i, :), [colors(i) '-'], 'LineWidth', 1.5);
    plot(1:numBands, endmembers_est(i, :), [colors(i) '--'], 'LineWidth', 1.5);
end
title('端元光谱对比');
xlabel('波段');
ylabel('反射率');
legend('植被(真)', '植被(估)', '水体(真)', '水体(估)', '裸土(真)', '裸土(估)');
grid on;

% 3. 混合光谱重建效果
subplot(1, 3, 3);
plot(1:numBands, mixed_pixel, 'ko-', 'LineWidth', 2, 'DisplayName', '原始混合光谱');
hold on;
plot(1:numBands, reconstruction, 'r--', 'LineWidth', 1.5, 'DisplayName', '重建光谱');
title(['光谱重建 (RMSE=' num2str(RMSE, '%.4f') ')']);
xlabel('波段');
ylabel('反射率');
legend;
grid on;


function [endmembers] = myVCA(data, num_endmembers)
    % 输入: data(m×n矩阵，m为像元数，n为波段数), num_endmembers(端元数量)
    % 输出: endmembers(端元光谱矩阵，n×num_endmembers)
    
    % 步骤1: 数据预处理（去均值化）
    mean_data = mean(data, 1);
    data_centered = data - mean_data;
    
    % 步骤2: PCA降维到(num_endmembers-1)维
    [coeff, ~, ~] = pca(data_centered);
    proj_matrix = coeff(:, 1:num_endmembers-1);
    data_projected = data_centered * proj_matrix;
    
    % 步骤3: 迭代寻找端元
    endmembers = zeros(num_endmembers, size(data,2));
    u = mean(data_projected);  % 初始投影方向
    for i = 1:num_endmembers
        % 计算投影值并寻找最远点
        projections = abs(data_projected * u');
        [~, idx] = max(projections);
        endmembers(i, :) = data(idx, :);
        
        % 更新投影方向（正交化）
        v = data_projected(idx, :) - u;
        u = v / norm(v);
        data_projected = data_projected - (data_projected * u') * u;
    end
end