混合像元分解模型误差分析

1.关于混合像元分解模型误差分析背景介绍

2.多波段平均RMSE求解方法如上（最后一个公式)  

  利用matlab求解：RMSE函数构造

function  [RMSE]=R_cal(data, fdata, endm)

% function [R R_ave R_rms] =R_cal(data, fdata, endm)

%data - 真实数据

%fdata - 丰度图（各组分在混合像元所占的比例）

%endm - 端元(各组分在各波段的反射率)

%Rre - 相对误差

%RMS - 均方根误差

[c r v]=size(data);

[c r fv]=size(fdata);

for i=1:c

for j=1:r

R(i,j,:)=fdata(i,j,1)*reshape(endm(:,1),1,1,v)+fdata(i,j,2)*reshape(endm(:,2),1,1,v)+...

fdata(i,j,3)*reshape(endm(:,3),1,1,v)+fdata(i,j,4)*reshape(endm(:,4),1,1,v);  %没个组分比例fdata()乘以各组分在 各波段的反射率，求和就得到(i,j,v)大小的三维矩阵

   end

end

R=R-data;   %残差  (i,j,v)大小的三维矩阵，每一个波段都有一个残差

R_ave=mean(R,3);   %平均残差  对第三维求平均值，就得到平均残差 (i,j,1)大小的二维矩阵

%%%%%%

for i=1:c

   for j=1:r

       RMSE(i,j)=sqrt(sum(R(i,j,:).*R(i,j,:))/v);  %求每个点（i,j)的rmse 各波段残差平方求和，除以波段数开根号就得到(i,j,1)大小的二维矩阵的平均残差

   end  

end

3.两种方法将结果写入到图片
%% ........................................方法1

data=imread('文件路径原始数据');

fdata=imread('文件路径丰度数据');

endm=load('文件路径端元数据');

[RMSE]=R_cal(data, fdata, endm);  %调用函数

imwrite(RMSE(:,:,1),'文件路径输出图片文件','tif'); %将RMSE数据写入图片tiff格式

%利用arcgis加载‘输出图片文件’,可以发现图像为灰度图，因此要还原数据真实值

%首先，将tiff输出为grid格式

%然后，给grid，添加字段value_orig,利用field calculator,计算字段值，公式为round(value/255,2)

%其次，将grid转化为点文件，字段选择value_orig

%最后，将点文件，转化为栅格文件，及还原了RMSE的真实值

%%  .......................................方法2

data=imread('文件路径原始数据');

fdata=imread('文件路径丰度数据');

endm=load('文件路径端元数据');

[RMSE]=R_cal(data, fdata, endm);  %调用函数

rmse=roundn(im2double(RMSE),-2);

dlmwrite('rmse.txt',rmse, 'delimiter',' ');%将rmse矩阵写入txt

%首先，找到该研究区头文件，给rmse.txt添加头文件

%然后，利用ASCII to Raster 工具，将 rmse转化为栅格文件，输出数据类型选浮点型float

%注：如果不知道头文件，自己手动添加

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%头文件格式

ncols         1522

nrows         1012

xllcorner     -5.5555549743902e-011

yllcorner     -5.5567994650119e-011

cellsize      0.013888889

NODATA_value  -9999

结果如下：(没有裁剪）