一个基于线性混合模型(Linear Mixing Model, LMM)的最简光谱解混 MATLAB 代码示例,包含模拟数据生成、端元提取、丰度反演及结果评估的全流程。代码基于全约束最小二乘法(Fully Constrained Least Squares, FCLS)实现,确保丰度满足非负性(abundances ≥ 0)和归一性(Σ abundances = 1)约束。
代码:
%% 步骤1: 模拟数据生成(3个端元,100个波段)
% 定义端元光谱矩阵(行向量为端元,列向量为波段)
endmembers = [0.8, 0.7, 0.6, 0.5, 0.4; % 植被(绿色曲线)
0.1, 0.2, 0.3, 0.2, 0.1; % 水体(蓝色曲线)
0.3, 0.4, 0.5, 0.6, 0.7]; % 裸土(棕色曲线)
[numEndmembers, numBands] = size(endmembers);
% 生成混合像元(端元按比例混合)
abundances_true = [0.3; 0.5; 0.2]; % 真实丰度向量(和=1)
mixed_pixel = endmembers' * abundances_true; % 线性混合模型:Y = M*A
% 添加高斯噪声(信噪比SNR=30dB)
SNR = 30;
noise_power = norm(mixed_pixel)^2 / (length(mixed_pixel) * 10^(SNR/10));
noise = sqrt(noise_power) * randn(numBands, 1);
mixed_pixel_noisy = mixed_pixel + noise; % 含噪混合像元
%% 步骤2: 端元提取(假设已知端元光谱,实际可用VCA/PPI提取[6](@ref))
% 此处直接使用模拟的端元矩阵作为已知端元
M_known = endmembers'; % 端元矩阵转置为波段×端元数
%% 步骤3: 丰度反演(FCLS:全约束最小二乘)
A_est = lsqnonneg(M_known, mixed_pixel_noisy); % 非负约束
A_est = A_est / sum(A_est); % 丰度归一化
%% 步骤4: 结果评估与可视化
% 计算重建误差(RMSE)
reconstruction = M_known * A_est;
RMSE = sqrt(mean((mixed_pixel_noisy - reconstruction).^2));
fprintf('重建RMSE: %.4f\n', RMSE);
% 绘制真实丰度与估计丰度对比
figure;
subplot(2,1,1);
bar([abundances_true, A_est]);
legend('真实丰度', '估计丰度');
title('丰度对比');
ylabel('比例');
set(gca, 'XTickLabel', {'植被', '水体', '裸土'});
% 绘制光谱曲线对比
subplot(2,1,2);
plot(1:numBands, mixed_pixel_noisy, 'ko-', 'LineWidth', 1.5, 'DisplayName', '混合像元(含噪)');
hold on;
plot(1:numBands, reconstruction, 'r--', 'LineWidth', 1.5, 'DisplayName', '重建光谱');
title('光谱重建结果');
xlabel('波段');
ylabel('反射率');
legend;
grid on;
公式:
其中,y为混合像元光谱,为端元矩阵,a为丰度向量,n为噪声。
#丰度反演
求解非负约束最小二乘问题:
归一化。