优化的高光谱解混算法实现

高光谱解混是遥感图像处理中的重要技术，用于从混合像元中提取纯光谱特征(端元)和它们的比例(丰度)。

% 优化的高光谱解混算法
% 包含VCA、FCLS、SUnSAL、CLSUnSAL和基于深度学习的解混方法

clear;
close all;
clc;
warning off;

%% 参数设置
dataType = 'synthetic'; % 'synthetic' 或 'real'
numEndmembers = 5;      % 端元数量
sparsityLevel = 0.1;    % 丰度稀疏性水平(仅对合成数据有效)
SNR = 30;               % 信噪比(dB)

% 算法参数
maxIter = 100;          % 最大迭代次数
lambda = 0.1;           % 稀疏性参数
mu = 0.01;              % 平滑性参数

%% 数据生成/加载
if strcmp(dataType, 'synthetic')
    % 生成合成高光谱数据
    [Y, E, A, rows, cols, bands] = generateSyntheticData(numEndmembers, sparsityLevel, SNR);
else
    % 加载真实高光谱数据(这里使用示例数据，实际应用中需替换为真实数据)
    load('samson_data.mat'); % 假设有samson_data.mat文件
    Y = data; % 数据矩阵[bands x pixels]
    rows = 95; cols = 95; bands = size(Y, 1);
    numEndmembers = 3;
end

% 数据显示
figure('Name', '高光谱数据立方体');
subplot(1,2,1);
imagesc(reshape(Y(50,:), rows, cols)); 
title('波段50图像');
colorbar;

subplot(1,2,2);
plot(Y(:,1:100:end)); 
title('随机像元光谱');
xlabel('波段'); ylabel('反射率');

%% 端元提取 - VCA算法
disp('执行VCA端元提取...');
E_vca = vca(Y, 'Endmembers', numEndmembers, 'SNR', SNR);
fprintf('VCA端元提取完成\n');

% 显示提取的端元
figure('Name', 'VCA提取的端元');
plot(E_vca);
title('VCA提取的端元光谱');
xlabel('波段'); ylabel('反射率');

%% 丰度估计 - 多种算法比较

% 1. 全约束最小二乘法(FCLS)
disp('执行FCLS丰度估计...');
A_fcls = FCLS(E_vca, Y);
fprintf('FCLS丰度估计完成\n');

% 2. 稀疏解混算法(SUnSAL)
disp('执行SUnSAL丰度估计...');
A_sunsal = sunsal(E_vca, Y, 'lambda', lambda, 'AL_iters', maxIter);
fprintf('SUnSAL丰度估计完成\n');

% 3. 协作稀疏解混算法(CLSUnSAL)
disp('执行CLSUnSAL丰度估计...');
A_clsunsal = clsunsal(E_vca, Y, 'lambda', lambda, 'mu', mu, 'iter', maxIter);
fprintf('CLSUnSAL丰度估计完成\n');

% 4. 基于深度学习的解混(需要Deep Learning Toolbox)
try
    disp('执行深度学习解混...');
    A_dl = deep_unmixing(E_vca, Y, rows, cols);
    fprintf('深度学习解混完成\n');
    useDL = true;
catch
    warning('深度学习解混不可用，需要Deep Learning Toolbox');
    useDL = false;
    A_dl = [];
end

%% 结果可视化
% 显示丰度图
figure('Name', '丰度估计结果比较');
for i = 1:numEndmembers
    subplot(4, numEndmembers, i);
    imagesc(reshape(A_fcls(i,:), rows, cols));
    title(sprintf('FCLS 端元%d', i));
    
    subplot(4, numEndmembers, i + numEndmembers);
    imagesc(reshape(A_sunsal(i,:), rows, cols));
    title(sprintf('SUnSAL 端元%d', i));
    
    subplot(4, numEndmembers, i + 2*numEndmembers);
    imagesc(reshape(A_clsunsal(i,:), rows, cols));
    title(sprintf('CLSUnSAL 端元%d', i));
    
    if useDL
        subplot(4, numEndmembers, i + 3*numEndmembers);
        imagesc(reshape(A_dl(i,:), rows, cols));
        title(sprintf('DL 端元%d', i));
    end
end
colormap jet;

% 计算和显示重建误差
Y_recon_fcls = E_vca * A_fcls;
reconstruction_error_fcls = norm(Y - Y_recon_fcls, 'fro') / norm(Y, 'fro');

Y_recon_sunsal = E_vca * A_sunsal;
reconstruction_error_sunsal = norm(Y - Y_recon_sunsal, 'fro') / norm(Y, 'fro');

Y_recon_clsunsal = E_vca * A_clsunsal;
reconstruction_error_clsunsal = norm(Y - Y_recon_clsunsal, 'fro') / norm(Y, 'fro');

if useDL
    Y_recon_dl = E_vca * A_dl;
    reconstruction_error_dl = norm(Y - Y_recon_dl, 'fro') / norm(Y, 'fro');
end

fprintf('\n重建误差:\n');
fprintf('FCLS: %.4f\n', reconstruction_error_fcls);
fprintf('SUnSAL: %.4f\n', reconstruction_error_sunsal);
fprintf('CLSUnSAL: %.4f\n', reconstruction_error_clsunsal);
if useDL
    fprintf('深度学习: %.4f\n', reconstruction_error_dl);
end

%% 端元提取函数 - VCA
function E = vca(Y, varargin)
    % VCA算法实现
    % 参数解析
    p = inputParser;
    addParameter(p, 'Endmembers', 0, @isnumeric);
    addParameter(p, 'SNR', 0, @isnumeric);
    parse(p, varargin{:});
    
    numEndmembers = p.Results.Endmembers;
    SNR = p.Results.SNR;
    
    [bands, pixels] = size(Y);
    
    % 估计噪声方差
    if SNR == 0
        % 自动估计SNR
        Y_mean = mean(Y, 2);
        Y_centered = Y - Y_mean;
        sigma2 = mean(var(Y_centered, 0, 2));
        SNR_est = 10 * log10(mean(Y_mean.^2) / sigma2);
        SNR = SNR_est;
    end
    
    % 数据投影到信号子空间
    if SNR > 15 + 10 * log10(numEndmembers)
        % 高SNR情况
        [U, S, ~] = svd(Y * Y' / pixels);
        Up = U(:, 1:numEndmembers);
        Y_proj = Up' * Y;
    else
        % 低SNR情况
        Y_mean = mean(Y, 2);
        Y_centered = Y - Y_mean;
        [U, S, ~] = svd(Y_centered * Y_centered' / pixels);
        Up = U(:, 1:numEndmembers-1);
        Y_proj = Up' * Y;
        Y_proj = [Y_proj; ones(1, pixels)];
    end
    
    % 初始化端元矩阵
    E = zeros(numEndmembers, numEndmembers);
    A = Y_proj;
    
    % 寻找第一个端元(投影数据中最大范数的向量)
    [~, idx] = max(sum(A.^2, 1));
    E(:, 1) = A(:, idx);
    
    % 初始化投影矩阵
    for i = 2:numEndmembers
        % 计算与已有端元正交的投影矩阵
        if i > 2
            P = eye(numEndmembers) - E(:, 1:i-1) * pinv(E(:, 1:i-1));
        else
            u = E(:, 1);
            P = eye(numEndmembers) - u * u' / (u' * u);
        end
        
        % 投影所有向量
        A_proj = P * A;
        
        % 找到投影后范数最大的向量
        [~, idx] = max(sum(A_proj.^2, 1));
        E(:, i) = A(:, idx);
    end
    
    % 如果使用了信号子空间投影，需要将端元映射回原始空间
    if SNR > 15 + 10 * log10(numEndmembers)
        E = Up * E;
    else
        E = Up * E(1:end-1, :) + Y_mean;
    end
end

%% 丰度估计函数 - FCLS
function A = FCLS(E, Y)
    % 全约束最小二乘丰度估计
    [bands, pixels] = size(Y);
    numEndmembers = size(E, 2);
    
    % 构建扩展矩阵以包含和为一约束
    delta = 1e-5; % 正则化参数
    E_ext = [E; delta * ones(1, numEndmembers)];
    Y_ext = [Y; delta * ones(1, pixels)];
    
    % 使用非负最小二乘
    A = zeros(numEndmembers, pixels);
    for i = 1:pixels
        A(:, i) = lsqnonneg(E_ext, Y_ext(:, i));
    end
    
    % 归一化丰度(和为1)
    A = A ./ sum(A, 1);
end

%% 丰度估计函数 - SUnSAL
function A = sunsal(E, Y, varargin)
    % 稀疏解混算法(SUnSAL)
    % 参数解析
    p = inputParser;
    addParameter(p, 'lambda', 0.1, @isnumeric);
    addParameter(p, 'AL_iters', 100, @isnumeric);
    parse(p, varargin{:});
    
    lambda = p.Results.lambda;
    maxIter = p.Results.AL_iters;
    
    [bands, pixels] = size(Y);
    numEndmembers = size(E, 2);
    
    % 使用交替方向乘子法(ADMM)求解
    % 问题形式: min (1/2)||Y - EA||_F^2 + lambda||A||_1 s.t. A >= 0, 1^T A = 1
    
    % 初始化
    A = zeros(numEndmembers, pixels);
    Z = A;
    U = zeros(numEndmembers, pixels);
    
    % 预计算矩阵
    EtE = E' * E;
    EtY = E' * Y;
    inv_EtE = inv(EtE + eye(numEndmembers));
    
    % ADMM迭代
    for iter = 1:maxIter
        % A更新步骤
        A = inv_EtE * (EtY + Z - U);
        
        % Z更新步骤(带有约束的投影)
        Z_prev = Z;
        Z = soft_threshold(A + U, lambda);
        
        % 应用非负约束和和为一约束
        Z = max(Z, 0);
        Z = Z ./ sum(Z, 1);
        
        % 对偶变量更新
        U = U + A - Z;
        
        % 检查收敛性
        if norm(Z - Z_prev, 'fro') < 1e-6
            break;
        end
    end
    
    A = Z;
end

%% 丰度估计函数 - CLSUnSAL
function A = clsunsal(E, Y, varargin)
    % 协作稀疏解混算法(CLSUnSAL)
    % 参数解析
    p = inputParser;
    addParameter(p, 'lambda', 0.1, @isnumeric);
    addParameter(p, 'mu', 0.01, @isnumeric);
    addParameter(p, 'iter', 100, @isnumeric);
    parse(p, varargin{:});
    
    lambda = p.Results.lambda;
    mu = p.Results.mu;
    maxIter = p.Results.iter;
    
    [bands, pixels] = size(Y);
    numEndmembers = size(E, 2);
    
    % 使用ADMM求解协作稀疏问题
    % 问题形式: min (1/2)||Y - EA||_F^2 + lambda||A||_{2,1} + mu TV(A) s.t. A >= 0, 1^T A = 1
    
    % 初始化
    A = zeros(numEndmembers, pixels);
    Z = A;
    U = zeros(numEndmembers, pixels);
    V = zeros(numEndmembers, pixels);
    
    % 预计算矩阵
    EtE = E' * E;
    EtY = E' * Y;
    inv_EtE = inv(EtE + eye(numEndmembers));
    
    % 定义TV正则化项(各向异性TV)
    Dx = diff_matrix(rows);
    Dy = diff_matrix(cols);
    
    % ADMM迭代
    for iter = 1:maxIter
        % A更新步骤
        A_prev = A;
        A = inv_EtE * (EtY + Z - U + mu * (Dx' * V(:, 1:end-1) + Dy' * V(:, 1:end-1)));
        
        % Z更新步骤(协作软阈值)
        Z_prev = Z;
        Z = group_soft_threshold(A + U, lambda);
        
        % 应用非负约束和和为一约束
        Z = max(Z, 0);
        Z = Z ./ sum(Z, 1);
        
        % V更新步骤(TV项)
        V = soft_threshold(Dx * A + Dy * A, 1);
        
        % 对偶变量更新
        U = U + A - Z;
        
        % 检查收敛性
        if norm(Z - Z_prev, 'fro') < 1e-6 && norm(A - A_prev, 'fro') < 1e-6
            break;
        end
    end
    
    A = Z;
end

%% 深度学习解混函数
function A = deep_unmixing(E, Y, rows, cols)
    % 基于深度学习的解混方法
    % 需要Deep Learning Toolbox
    
    [bands, pixels] = size(Y);
    numEndmembers = size(E, 2);
    
    % 准备数据
    X = reshape(Y', rows, cols, bands); % 转换为图像格式
    X = normalize(X, 1, 'range'); % 归一化到[0,1]
    
    % 定义网络架构
    layers = [
        imageInputLayer([rows cols bands], 'Normalization', 'none', 'Name', 'input')
        
        convolution2dLayer(3, 32, 'Padding', 'same', 'Name', 'conv1')
        batchNormalizationLayer('Name', 'bn1')
        reluLayer('Name', 'relu1')
        
        convolution2dLayer(3, 64, 'Padding', 'same', 'Name', 'conv2')
        batchNormalizationLayer('Name', 'bn2')
        reluLayer('Name', 'relu2')
        
        convolution2dLayer(3, 128, 'Padding', 'same', 'Name', 'conv3')
        batchNormalizationLayer('Name', 'bn3')
        reluLayer('Name', 'relu3')
        
        convolution2dLayer(3, numEndmembers, 'Padding', 'same', 'Name', 'conv4')
        softmaxLayer('Name', 'softmax')
        regressionLayer('Name', 'output')
    ];
    
    % 训练选项
    options = trainingOptions('adam', ...
        'MaxEpochs', 50, ...
        'MiniBatchSize', 16, ...
        'InitialLearnRate', 1e-3, ...
        'LearnRateSchedule', 'piecewise', ...
        'LearnRateDropFactor', 0.5, ...
        'LearnRateDropPeriod', 20, ...
        'Shuffle', 'every-epoch', ...
        'Verbose', false);
    
    % 生成训练目标(使用FCLS作为伪标签)
    A_target = FCLS(E, Y);
    Y_target = reshape(A_target', rows, cols, numEndmembers);
    
    % 训练网络
    net = trainNetwork(X, Y_target, layers, options);
    
    % 预测丰度
    A_pred = predict(net, X);
    A = reshape(A_pred, pixels, numEndmembers)';
end

%% 辅助函数 - 软阈值
function X = soft_threshold(Z, lambda)
    % 软阈值函数
    X = sign(Z) .* max(abs(Z) - lambda, 0);
end

%% 辅助函数 - 组软阈值
function X = group_soft_threshold(Z, lambda)
    % 组软阈值函数(用于协作稀疏)
    norms = sqrt(sum(Z.^2, 1));
    scale = max(1 - lambda ./ (norms + eps), 0);
    X = Z .* scale;
end

%% 辅助函数 - 差分矩阵
function D = diff_matrix(n)
    % 创建一维差分矩阵
    D = diag(-ones(n,1)) + diag(ones(n-1,1), 1);
    D = D(1:end-1, :);
end

%% 辅助函数 - 生成合成数据
function [Y, E, A, rows, cols, bands] = generateSyntheticData(numEndmembers, sparsity, SNR)
    % 生成合成高光谱数据
    
    % 设置参数
    rows = 64;
    cols = 64;
    pixels = rows * cols;
    bands = 200;
    
    % 生成端元光谱(使用USGS光谱库或随机生成)
    E = rand(bands, numEndmembers);
    
    % 生成稀疏丰度图
    A = zeros(numEndmembers, pixels);
    for i = 1:numEndmembers
        % 创建空间平滑的丰度图
        abundance = smooth_abundance(rows, cols, sparsity);
        A(i, :) = abundance(:);
    end
    
    % 确保丰度和为1
    A = A ./ sum(A, 1);
    
    % 生成高光谱数据
    Y = E * A;
    
    % 添加噪声
    signal_power = norm(Y, 'fro')^2 / (bands * pixels);
    noise_power = signal_power / (10^(SNR/10));
    noise = sqrt(noise_power) * randn(size(Y));
    Y = Y + noise;
end

%% 辅助函数 - 生成平滑丰度图
function A = smooth_abundance(rows, cols, sparsity)
    % 生成空间平滑的丰度图
    
    % 创建随机点作为高丰度区域中心
    num_points = max(1, round(sparsity * rows * cols));
    points = rand(num_points, 2);
    points(:,1) = round(points(:,1) * rows);
    points(:,2) = round(points(:,2) * cols);
    
    % 初始化丰度图
    A = zeros(rows, cols);
    
    % 为每个点创建高斯分布
    for i = 1:num_points
        r = points(i, 1);
        c = points(i, 2);
        
        % 创建高斯核
        [X, Y] = meshgrid(1:cols, 1:rows);
        kernel = exp(-((X - c).^2 + (Y - r).^2) / (2 * (min(rows, cols)/10)^2));
        
        A = A + kernel;
    end
    
    % 归一化到[0,1]
    A = A / max(A(:));
end

说明

这个高光谱解混算法实现包含以下主要部分：

1. 数据生成/加载

• 支持合成数据和真实数据
• 合成数据生成考虑了空间平滑性和稀疏性

2. 端元提取

• 实现了VCA(Vertex Component Analysis)算法
• 自动估计信号子空间维度
• 适应不同信噪比条件

3. 丰度估计算法

• FCLS(Fully Constrained Least Squares)：全约束最小二乘法
• SUnSAL(Sparse Unmixing by Variable Splitting and Augmented Lagrangian)：稀疏解混算法
• CLSUnSAL(Collaborative SUnSAL)：协作稀疏解混算法，结合空间上下文信息
• 深度学习解混：基于卷积神经网络的端到端解混方法

4. 优化技术

• 使用ADMM(Alternating Direction Method of Multipliers)求解优化问题
• 包含稀疏约束(L1正则化)
• 包含空间平滑约束(TV正则化)
• 协作稀疏约束(L2,1正则化)

5. 结果评估与可视化

• 丰度图可视化
• 重建误差计算
• 算法性能比较

使用

1. 设置参数(dataType, numEndmembers, sparsityLevel, SNR)

2. 运行代码，算法将自动执行以下步骤：

   • 数据生成/加载
   • VCA端元提取
   • 多种丰度估计算法
   • 结果可视化和性能评估

3. 对于真实数据应用，需要：

   • 准备真实高光谱数据文件
   • 调整参数以适应实际数据特性
   • 可能需要调整正则化参数(λ, μ)以获得最佳结果

推荐代码 优化的高光谱解混算法 www.youwenfan.com/contentcnl/54929.html

算法特点

1. VCA端元提取：基于几何原理，能够有效识别数据凸包顶点
2. SUnSAL：利用稀疏性先验，适合端元数量未知的情况
3. CLSUnSAL：结合空间上下文信息，提高解混精度
4. 深度学习方法：端到端学习，能够捕捉复杂非线性关系