高光谱解混是遥感图像处理中的重要技术，用于从混合像元中提取纯光谱特征(端元)和它们的比例(丰度)。
% 优化的高光谱解混算法
% 包含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)
trydisp('执行深度学习解混...');A_dl = deep_unmixing(E_vca, Y, rows, cols);fprintf('深度学习解混完成\n');useDL = true;
catchwarning('深度学习解混不可用，需要Deep Learning Toolbox');useDL = false;A_dl = [];
end%% 结果可视化
% 显示丰度图
figure('Name', '丰度估计结果比较');
for i = 1:numEndmemberssubplot(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 useDLsubplot(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 useDLY_recon_dl = E_vca * A_dl;reconstruction_error_dl = norm(Y - Y_recon_dl, 'fro') / norm(Y, 'fro');
endfprintf('\n重建误差:\n');
fprintf('FCLS: %.4f\n', reconstruction_error_fcls);
fprintf('SUnSAL: %.4f\n', reconstruction_error_sunsal);
fprintf('CLSUnSAL: %.4f\n', reconstruction_error_clsunsal);
if useDLfprintf('深度学习: %.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% 自动估计SNRY_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 > 2P = eye(numEndmembers) - E(:, 1:i-1) * pinv(E(:, 1:i-1));elseu = 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;elseE = 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:pixelsA(:, 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-6break;endendA = 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-6break;endendA = 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% 确保丰度和为1A = 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_pointsr = 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. 结果评估与可视化
- 丰度图可视化
- 重建误差计算
- 算法性能比较
使用
-
设置参数(
dataType,numEndmembers,sparsityLevel,SNR) -
运行代码，算法将自动执行以下步骤：
- 数据生成/加载
- VCA端元提取
- 多种丰度估计算法
- 结果可视化和性能评估
-
对于真实数据应用，需要：
- 准备真实高光谱数据文件
- 调整参数以适应实际数据特性
- 可能需要调整正则化参数(λ, μ)以获得最佳结果
推荐代码 优化的高光谱解混算法 www.youwenfan.com/contentcnl/54929.html
算法特点
- VCA端元提取：基于几何原理，能够有效识别数据凸包顶点
- SUnSAL：利用稀疏性先验，适合端元数量未知的情况
- CLSUnSAL：结合空间上下文信息，提高解混精度
- 深度学习方法：端到端学习，能够捕捉复杂非线性关系
