（四）Landat_5 TM 遥感影像利用6s模型大气校正以及地表反射率计算

注：本文章只是实验笔记，记录求解过程，其中包含很多错误因素，欢迎读者指正！

校正过程

一、6S模型所需参数输入

1. 几何条件 Geometrical conditions
   igeom [0-7]: 7 其代表的是传感器的类型，由于我的是Landsat影像，故igeom=1
   其中SPOT卫星代表的igeom为6
2. 时间、经度、纬度 (month,day,hour,long,lat)
   通过查询头文件，可以发现以上值：

其中经纬度为影像四个角点的位置，我们求平均得到以下输入格式：
8 2 2.7 107.56525 30.30198

1. 大气模型 Atmospheric model
   idatm [0-8]: 2 其代表的该影像地区的大气类型，由于所在区为中纬度夏季大气模式，故idatm=2

2. 气溶胶类型 Aerosol model(type)
   iaer[0-12]: 3

3. 气溶胶浓度 Aerosol model (concentration)
   大气厚度：23 （单位km）

4. 目标高 The altitude of target
   xps: -0.02 (目标高度＝0.02km)

5. 传感器的高度 The sensor altitude
   xpp: -1000 (卫星高度＝1000km)

6. 波段状况 The spectral conditions
   iwave[25-30]:25 （25即TM band1, 26即TM band2，27即TM band3，28即TM band4，29即TM band5, 30即TM band6）

7. 地面反射率 Ground reflectance (type)
   inhomo: 0 (地面均匀)
   idirec: 0 (地面无方向影响)
   igroun: 0 (反射率不随波长变化)
   ro: 0.2 (地面波段反射率＝0.2)

8. 激活大气订正的方式 Atmospheric correction mode
   rapp: 0
   我们将以上参数传入6s模型，首先window+R输入cmd，在cmd中输入6s模型的路径，然后再键入以下参数：
   
   
   重复上述步骤：

波段	xa、xb、xc
band1	0.00322、0.12016、0.12834
band2	0.00321、0.06567、0.09102
band3	0.00341、0.033703、0.067
band4	0.00481、0.01723、0.04092
band5	0.00676、0.00261、0.01069
band6	0.05805、0.00109、0.00482（不建议使用）

二、开始校正

根据6s模型中的公式，我们可以将其整理得到以下结果（式中：measured radiance 表示辐射度 acr 表示地表反射率）：

xaxbxc=[ 0.00322,0.12016,0.12834;
    0.00321,0.06567,0.09102;
    0.00341,0.03703,0.067;
    0.00481,0.01723,0.04092;
    0.00676,0.00261,0.01069]
measured_radiance={image{1},image{2},image{3},image{4},...
   image{5} };
acr={};
y={};
for oi=1:length(measured_radiance)
    y{oi}=xaxbxc(oi,1).*double(measured_radiance{oi})-xaxbxc(oi,2);
    acr{oi}=y{oi}./(1+xaxbxc(oi,3).*y{oi});
end
%绘制波谱反射曲线
msgbox('请在图像中选取3个点！！！');
figure;
imshow(hend.*3);
[x_b,y_b]=ginput;
hold on
scatter(x_b,y_b,1000,'r','+');
point_data={};
for ya=1:length(Reflectance_value)
     map=double(Reflectance_value{ya});
     yi=map(int16(y_b),int16(x_b));
     point_data{ya}=diag(yi);
end
matrix=cell2mat(point_data);
ld_x=0.3:0.4:2.7;
ld_y=matrix(1,:);
ld_y(ld_y==Inf)=0;
figure;
plot(ld_x,ld_y,'Color','r','LineWidth',1);
xlabel('波长（微米）');
ylabel('大气表观反射率');
title('地物光谱反射率');
hold on
ld_y2=matrix(2,:);
ld_y2(ld_y2==Inf)=0;
plot(ld_x,ld_y2,'Color','b','LineWidth',1);
hold on
ld_y3=matrix(3,:);
ld_y3(ld_y3==Inf)=0;
plot(ld_x,ld_y3,'Color','g','LineWidth',1);
legend('point 1', 'point 2', 'point 3');
hold off;
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
point_data2={};
for yz=1:length(acr)
    map2=double(acr{yz});
    yi2=map2(int16(y_b),int16(x_b));
    point_data2{yz}=diag(yi2);
end
matrix2=cell2mat(point_data2);
cd_x=0.2:0.4:1.8;
cd_y=matrix2(1,:);
cd_y(cd_y==Inf)=0;
figure;
axis([0,1.8,0,1]);
plot(cd_x,cd_y,'Color','r','LineWidth',1);
xlabel('波长（微米）');
ylabel('地表反射率');
title('地物光谱反射率');
hold on
cd_y2=matrix2(2,:);
cd_y2(cd_y2==Inf)=0;
plot(cd_x,cd_y2,'Color','b','LineWidth',1);
hold on
cd_y3=matrix2(3,:);
cd_y3(cd_y3==Inf)=0;
plot(cd_x,cd_y3,'Color','g','LineWidth',1);
legend('point 1', 'point 2', 'point 3');
hold off;