CN104156567A - 一种耦合卫星遥感影像大气校正和地形校正过程的地表反射率获取技术 - Google Patents

Abstract

一种耦合卫星遥感影像大气校正和地形校正过程的地表反射率获取技术，将遥感影像大气校正和地形校正过程有机地结合起来，实现卫星遥感影像地表反射率的有效反演。该方法基于遥感影像大气校正和地形校正的物理机制，是一种物理方法，具有普适性。

Description

一种耦合卫星遥感影像大气校正和地形校正过程的地表反射率获取技术

技术领域

本发明涉及一种卫星遥感影像地表反射率获取的方法，能够应用在农业、林业、气象、生态环境以及国防安全等行业部门。 

背景技术

卫星传感器在获取地表信息过程中会受到多种因素的干扰，从而造成卫星传感器获取信息的失真，大气和地形影响是造成卫星传感器信息失真的最重要的两个因素。为了从卫星遥感数据高精度地获取地表信息(地表反射率)，必须消除或者减弱大气和地形因素对卫星遥感数据的影响，这个过程分别称为大气校正和地形校正。 

随着定量遥感的发展，国内外研究者围绕遥感影像大气校正和地形校正开展了大量工作，提出了不少大气校正和地形校正方法。卫星遥感影像大气校正的主要目的是减少大气吸收和散射对卫星信号的影响。在大气校正方面，比较常用的方法包括：辐射传输模型法、不变目标法、暗目标减法、直方图匹配法等。在山区，地表接收的辐射能量和传感器所接收的信号要受到地形起伏的影响。表现在遥感图像上，阴坡上的像元接收到的入射辐照度较弱、亮度值较低；与此相反，阳坡上的像元接收到的入射辐照度较强、亮度值较高，造成“同物异谱”、“异物同谱”问题，这种光谱信息的失真，严重影响了山区遥感图像的信息提取精度，成为山区遥感图像定量应用的一个障碍。针对地形效应问题，近年来陆续提出了一些地形校正方法，比较常用的方法包括C模型、SCS模型、Minnaert模型等。 

目前关于遥感影像大气校正和地形校正方法的研究存在两个明显的问题：第一，不少的大气校正和地形校正方法是经验性的，普适性差，例如常用的不变目标法、暗目标减法以及C地形校正模型等；第二，大气校正和地形校正算法分别基于不同的理论基础，大气校正算法大多基于大气辐射传输过程，而地形校正算法则基于太阳-地表-传感器三者的几何关系，因此大气校正和地形校正算法分别独立开展研究，大气校正和地形校正被视为两个独立的环节，在实际操作过程中，通常先进行大气校正，然后进行地形校正。然而在卫星传感器获取地表信息过程中，大气影响和地形影响是交织在一起的，这种大气校正和地形校正相互割裂的做法必然对地表信息的获取产生影响。本发明从大气校正和地形校正的理论基础出发，将遥感影像大气校正和地形校正过程有机地结合起来，实现了卫星遥感影像地表反射率的有效反演。该发明充分发挥了遥感数据面特征的优势，而且它基于遥感影像大气校正和地形校正的物理机制，是一种物理方法，具有普适性，适用于所有的陆地卫星遥感影像。 

发明内容

本发明的目的在于提供一种基于物理机制的卫星遥感影像地表反射率获取技术，以克服现有方法的的不足，从而促进相关领域遥感应用的发展。 

为实现上述目的，本发明提出的方法包括以下步骤： 

第一步、对卫星影像进行几何校正； 

第二步、辐射定标，将影像像元亮度值转换为星上辐射亮度； 

L_λ＝Gain·DN+Bias 

L_λ为波段λ的星上辐射亮度，DN为波段λ的像元值，Gain和Bias分别为波段λ的增益和偏置。 

第三步、获取大气程辐射 

将直方图中频数累加和等于总象素个数0.01％的像素DN值作为暗目标像元值DN_min，程辐射L_p的计算式变为： 

L_p＝Gain·DN_min+Bias 

第四步、根据陆地卫星遥感影像的经纬度和影像成像时间，获取同一地理范围、同步过境的MODIS气溶胶光学厚度(0.47μm和0.66μm)、总大气水蒸汽含量、地表气压和臭氧浓度影像，并将各影像做几何校正，转换为和陆地卫星遥感影像一致的投影方式和分辨率，且影像行列数一致。 

第五步、由总大气水蒸汽含量W计算大气水蒸汽光学厚度； 

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$

其中，τ_w是大气水蒸汽光学厚度，a_wλ是水汽吸收系数，w是可降水汽，可降水汽w与总大气水蒸汽含量W在数值上相等，M是相对大气量； 

第六步、计算臭氧吸收光学厚度 

τ_o＝C_ozone*A_ozone(λ) 

τ_o是臭氧吸收光学厚度，C_ozone为MODIS臭氧浓度(单位为Dobson)，A_ozone(λ)为波段λ的臭氧吸收系数； 

第七步、由MODIS 0.47μm和0.66μm的气溶胶光学厚度，利用下式计算波长指数(α)及大气浑浊度系数(β) 

$\mathrm{\α}=-\frac{{1\mathrm{n\τ}}_a\left({\mathrm{\λ}}_1\right)-{1\mathrm{n\τ}}_a\left({\mathrm{\λ}}_2\right)}{{1\mathrm{n\λ}}_1-{1\mathrm{n\λ}}_2}$

β＝τ_a(λ₁)λ₁ ^α＝τ_a(λ₂)λ₂ ^α

λ₁和λ₂分别为0.47μm和0.66μm，τ_a(λ₁)和τ_a(λ₂)分别为0.47μm和0.66μm的气溶胶光学厚度； 

在得到α和β后，可以利用下式得到任意波长λ的气溶胶光学厚度τ_a(λ)： 

τ_a(λ)＝βλ^-α

其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度； 

第八步、计算瑞利散射光学厚度； 

${\mathrm{\τ}}_r=\frac{p}{1013.25}0.008569{\mathrm{\λ}}^{-4}(1+0.0133{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})$

其中，τ_r为瑞利散射光学厚度，λ是影像各波段的中心波长，μm；p为MODIS的地表气压产品； 

第九步、计算太阳照射方向和传感器观测方向的大气透过率； 

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_a-τ_w)/cosθ_z} 

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_a-τ_w)/cosθ_v} 

其中，T_z和T_v分别是太阳照射方向和传感器观测方向的大气透过率，τ是大气光学厚度，θ_z是太阳天顶角，θ_v是传感器观测天顶角，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度； 

第十步、计算天空散射辐射 

E_f＝(E_r+E_a)(3+cos(2S))/4 

其中，E_f为天空散射辐射，E_r为瑞利散射分量，E_a为气溶胶散射分量，S为倾斜坡面的坡度角； 

E_r＝0.5×E₀/d²cosθ_zexp{(-τ_o-τ_w)/cosθ_z}T_aa(1-exp{(-0.95τ_r)/cosθ_z}) 

E_a＝E₀/d²cosθ_zexp{(-τ_o-τ_w)/cosθ_z}T_aaexp{(-1.5τ_r)/cosθ_z}(1-T_as)F_s

其中，E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离，θ_z是太阳天顶角，τ_o、τ_w和τ_r分别是臭氧吸收光学厚度、水汽吸收光学厚度和瑞利散射光学厚度，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，F_s是气溶胶散射向下分量的比例； 

日地距离d²(天文单位)按下式计算： 

$\mathrm{\Γ}=\frac{360}{365}(\mathrm{dn}-1)$

$d^2=\frac{1}{1.00011+0.034221\cos \mathrm{\Γ}+0.00128\sin \mathrm{\Γ}+0.000719\cos \mathrm{\Γ}+0.00077\sin 2\mathrm{\Γ}}$

其中，dn是儒略日，即影像获取日期距离1月1日的天数； 

第十一步、计算地形反射辐射； 

E_s＝0.5×(E₀cosθ_zT_z+E_r+E_a)×ρ_mean×(1-cos S) 

其中E_s为坡面像元接收的周围地形的反射辐射，S为倾斜坡面的坡度角；ρ_mean表示相邻物体的平均反射率，取消除大气影响后的反射率影像的平均值来表示，具体计算式如下： 

${\mathrm{\ρ}}_{\mathrm{mean}}=\mathrm{average}\left(\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(E_0\cos {\mathrm{\θ}}_z{T}_z+E_r+E_a)}\right)$

第十二步、利用Dymond公式考虑地形对地表反射的影响，地表反射率的计算公式可以表示为： 

$\mathrm{\ρ}=\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(\frac{\mathrm{\Θ}{E}_0{T}_z(\cos i_h+{\cos e}_h)\cos i}{\cos i+\cos e}+E_f+E_s)}$

i和e分别为直射辐射在坡面上的入射角和出射角，i_h和e_h分别为直射辐射在水平面上的入射角和出射角，入射角i_h与太阳天顶角θ_z相等，则cosi_h＝cosθ_z； 

S为倾斜坡面的坡度角，为太阳方位角，A为倾斜坡面的坡向，θ_z为太阳天顶角； 

假设卫星传感器是星下点观测，则e_h＝0，cose_h＝1，且由几何关系可以推得：e＝S，则最终地表反射率的计算公式为： 

$\mathrm{\ρ}=\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(\frac{\mathrm{\Θ}{E}_0{T}_z(\cos {\mathrm{\θ}}_z+1)\cos i}{\cos i+\cos S}{+E}_f+E_s)}$

第十三步、对地表反射率反演结果进行验证。 

具体实施方式

遥感地表反射率反演是以辐射传输方程为基础的，在假定地表均匀、朗伯面反射、大气性质均一，忽略大气多次散射的前提下，星上辐射亮度和地表反射率的关系如式1所示： 

$L_{\mathrm{\λ}}=L_p+\frac{\mathrm{\ρ}{F}_d{T}_v}{\mathrm{\π}{(1-\mathrm{s\ρ})d}^2}$        (式1) 

式1中，L_λ是星上辐射亮度，L_p是程辐射，F_d是地表接收到的辐照度，T_v是传感器观测方向的大气透过率，ρ是地表反射率，d是日地距离，s是大气下界的半球反射率，因为s值很小，通常可以忽略，所以由式1可以得到地表反射率的计算公式(式2)： 

$\mathrm{\ρ}=\frac{{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d}^2}{T_v{F}_d}$         (式2) 

L_λ由传感器辐射定标得到 

L_λ＝Gain·DN+Bias         (式3) 

L_λ为波段λ的星上辐射亮度，DN为波段λ的像元值，Gain和Bias分别为波段λ的增益和偏置。 

将直方图中频数累加和等于总象素个数0.01％的像素DN值作为暗目标像元值(DN_min)，程辐射L_p的计算式为： 

L_P＝Gain·DN_min+Bias          (式4) 

传感器观测方向的大气透过率T_v按式5获得： 

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v}   (式5) 

其中，τ是大气光学厚度，θ_v是传感器观测天顶角，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度； 

${\mathrm{\τ}}_r=\frac{p}{1013.25}0.008569{\mathrm{\λ}}^{-4}(1+0.0133{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})$    (式6) 

其中，τ_r为瑞利散射光学厚度，λ是影像各波段的中心波长，μm；p为MODIS的地表气压产品； 

由MODIS 0.47μm和0.66μm的气溶胶光学厚度，利用式7-8计算波长指数(α)及大气浑浊度系数(β) 

$\mathrm{\α}=-\frac{{1\mathrm{n\τ}}_a\left({\mathrm{\λ}}_1\right)-{1\mathrm{n\τ}}_a\left({\mathrm{\λ}}_2\right)}{{1\mathrm{n\λ}}_1-{1\mathrm{n\λ}}_2}$    (式7) 

β＝τ_a(λ₁)λ₁ ^α＝τ_a(λ₂)λ₂ ^α    (式8) 

λ₁和λ₂分别为0.47μm和0.66μm，τ_a(λ₁)和τ_a(λ₂)分别为0.47μm和0.66μm的气溶胶光学 厚度； 

在得到α和β后，可以利用式9得到任意波长λ的气溶胶光学厚度τ_a(λ)： 

τ_a(λ)＝βλ^-α    (式9) 

其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度； 

臭氧吸收光学厚度τ_o按式10获得： 

τ_o＝C_ozone*A_ozone(λ)     (式10) 

τ_o是臭氧吸收光学厚度，C_ozone为MODIS臭氧浓度(单位为Dobson)，A_ozone(λ)为波段λ的臭氧吸收系数，可以根据文献Aerosol Optical Depth Value-Added Product(Koontz A et al.2013)得到。 

https：//www.arm.gov/publications/tech_reports/doe-sc-arm-tr-129.pdf 

大气水蒸汽光学厚度τ_w由总大气水蒸汽含量W计算： 

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$       (式11) 

其中，τ_w是大气水蒸汽光学厚度，a_wλ是水汽吸收系数，各波长的a_wλ值可以查文献得到[Bird R E，Riordan C.Simple solar spectral model for direct and diffuse irradiance on horizontal and tilted planes at the earth′s surface for cloudless atmospheres.Journal of climate and applied meteorology，1986，25：87-97.]。w是可降水汽，cm，可降水汽w与总大气水蒸汽含量W在数值上相等，M是相对大气量，由式12获得： 

M=[cosθ_z+0.15(93.885-θ_z)^-1.253]^-1      

θ_z是太阳天顶角，由影像头文件得到。 

日地距离d²(天文单位)按式13-14计算： 

$\mathrm{\Γ}=\frac{360}{365}(\mathrm{dn}-1)$       (式13) 

$d^2=\frac{1}{1.00011+0.034221\cos \mathrm{\Γ}+0.00128\sin \mathrm{\Γ}+0.000719\cos \mathrm{\Γ}+0.00077\sin 2\mathrm{\Γ}}$    (式14) 

其中，dn是儒略日，即影像获取日期距离1月1日的天数； 

地表接收到的辐照度F_d由三部分组成：太阳直射辐射(E_b)、天空散射辐射(E_f)和周围地形的反射辐射(E_s)，即F_d＝E_b+E_f+E_s，E_b＝ΘE₀cosθ_zT_zcosi/cosθ_z＝ΘE₀T_zcosi，i为坡面太阳入射角，Θ为地形阴影系数，若坡面为阴影区(cosi＜0)，Θ为0；若坡面为非阴影区(cosi≥0)，Θ为1。E₀是大气层外相应波长的太阳光谱辐照度，可由探测器响应函数计算得到。θ_z是太阳天顶角，由影像头文件获得。T_z是太阳照射方向上的大气透过率。 

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z}  (式15) 

cosi由式16得到： 

        (式16) 

S为倾斜坡面的坡度角，为太阳方位角，A为倾斜坡面的坡向，θ_z为太阳天顶角。太阳天顶角和太阳方位角可由卫星图像头文件得到，坡面坡度角和坡向可由DEM数据得到。 

坡面像元接收的天空散射辐射(E_f)由水平面像元接收的天空散射辐射(E_fh)转换得到。E_fh包括两部分：E_r和E_a，E_r为瑞利散射分量，E_a为气溶胶散射分量，假定瑞利散射向下的比例是50％，与太阳天顶角无关；同时假定气溶胶散射向下的比例是F_s，F_s是太阳天顶角的函数。则E_r和E_a的表达式如下： 

E_r＝0.5×E₀/d²cosθ_zexp{(-τ_o-τ_w)/cosθ_z}T_aa(1-exp{(-0.95τ_r)/cosθ_z})  (式17) 

E_a＝E₀/d²cosθ_zexp{(-τ_o-τ_w)/cosθ_z}T_aaexp{(-1.5τ_r)/cosθ_z}(1-T_as)F_s   (式18) 

其中，E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离，θ_z是太阳天顶角，τ_o、τ_w和τ_r分别是臭氧吸收光学厚度、水汽吸收光学厚度和瑞利散射光学厚度，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，F_s是气溶胶散射向下分量的比例，后三个参数的具体计算参见文献[Bird R E，Riordan C.Simple solar spectral model for direct and diffuse irradiance on horizontal and tilted planes at the earth′s surface for cloudless atmospheres.Journal of climate and applied meteorology，1986，25：87-97.]。 

在得到E_fh后，利用Badescu模型(Badescu V.3D isotropic approximation for solar diffuse irradiance on tilted surfaces.Renewable Energy，2002，26：221-233.)将E_fh转换为E_f

E_f＝E_fh(3+cos(2S))/4        (式19) 

周围地形的反射辐射按式20计算： 

E_s＝0.5×(E₀cosθ_zT_z+E_r+E_a)×ρ_mean×(1-cosS)   (式20) 

其中E_s为坡面像元接收的周围地形的反射辐射，S为倾斜坡面的坡度角，其它参数的意义及计算方法见上文，ρ_mean表示相邻物体的平均反射率，取消除大气影响后的反射率影像的平均值来表示，具体计算式如下： 

${\mathrm{\ρ}}_{\mathrm{mean}}=\mathrm{average}\left(\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(E_0\cos {\mathrm{\θ}}_z{T}_z+E_r+E_a)}\right)$        (式21) 

考虑地形对不同地表辐射反射的影响，并假定坡面对天空散射辐射和周围地形的反射辐 射的反射率相等，则式2可以进一步表示为： 

T_v{ρ_bE_b+ρ_f(E_f+E_s)}＝π(L_λ-L_p)d²   (式22) 

其中，E_b、E_f和E_s分别代表坡面像元接收的太阳直射辐射、天空散射辐射和周围地形的反射辐射。ρ_b为坡面对太阳直射辐射的反射率，ρ_f为坡面对天空散射辐射的反射率，Dymond等[Dymond J R，Shepherd J D.Correction of the topographic effect in remote sensing.IEEE Transactions on Geoscience and Remote Sensing，1999，37(5)：2618-2620.]提出了以下物理模型(式23)来考虑地形起伏对太阳直射辐射反射率的影响，本发明借鉴了Dymond模型。 

$\frac{{\mathrm{\ρ}}_b^h}{{\mathrm{\ρ}}_b}=\frac{\cos \left(i\right)+\cos \left(e\right)}{\cos \left(i_h\right)+\cos \left(e_h\right)}$     (式23) 

ρ_b为坡面对直射辐射的反射率，为水平面对直射辐射的反射率，i和e分别为直射辐射在坡面上的入射角和出射角，i_h和e_h分别为直射辐射在水平面上的入射角和出射角。 

由于卫星传感器所接收的信息大部分来源于地表对太阳直射辐射的反射，天空散射辐射和周围地形的反射辐射相对较小，地形对它们反射率的影响也较小。为了简化计算，假设  ${\mathrm{\ρ}}_f={\mathrm{\ρ}}_b^h,$ 于是： 

$T_v=\{{\mathrm{\ρ}}_b^h{E}_b/\mathrm{\δ}+{\mathrm{\ρ}}_b^h(E_f+E_s)\}=\mathrm{\π}{(L_{\mathrm{\λ}}-L_p)d}^2$    (式24) 

$\mathrm{\δ}=\frac{\cos \left(i\right)+\cos \left(e\right)}{\cos \left(i_h\right)+\cos \left(e_h\right)}$     (式25) 

由几何关系可以得到i_h＝θ_z，经过整理得到地表反射率的计算公式为： 

${\mathrm{\ρ}}_b^h=\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(\frac{\mathrm{\Θ}{E}_0{T}_z(\cos {\mathrm{\θ}}_z+\cos e_h)\cos i}{\cos i+\cos e}{+E}_f+E_s)}$      (式26) 

假设卫星传感器是星下点观测，则e_h＝0，cose_h＝1，且由几何关系可以推得：e＝S，则最终地表反射率的计算公式为： 

$\mathrm{\ρ}=\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(\frac{\mathrm{\Θ}{E}_0{T}_z(\cos {\mathrm{\θ}}_z+1)\cos i}{\cos i+\cos S}+E_f+E_s)}$     (式27) 

利用本发明的方法对2008年3月24日福建永安地区的IRS P6 LISS 3影像进行大气校正和地形校正，获取地表反射率。然后分别从大气效应校正效果和地形效应校正效果两个方面对本发明方法的有效性进行检验。 

为检验大气效应校正效果，在IRS P6 LISS 3影像上选取典型植被，然后分别提取表观反 射率和本发明方法校正后获取的地表反射率，得到图1。 

图1是本发明方法得到的地表反射率与表观反射率的对比图，系列1代表本发明方法得到的地表反射率，系列2代表表观反射率，横轴代表IRS P6 LISS 3的4个波段 

从图1可以看出，两条曲线在波段1(绿波段)和波段3(近红外波段)有较大差异，经本发明方法校正后，绿波段的大气散射效应和近红外波段的大气吸收效应得到抑制，绿波段和近红外波段的两个植被反射峰得以显示，这也表明了本发明方法取得了明显的大气校正效果。 

对地形效应校正效果的评价，借鉴文献中一个常用的定量评价指标，即“地表反射率和对应太阳入射角余弦的相关性”。校正前，像素的地表反射率与太阳入射角余弦(cosi)的相关性明显，像素地表反射率明显地随着太阳入射角的变化而变化。经过地形校正后，两者的相关性减弱。为了计算两者的相关性，从表观反射率影像、本发明方法得到的地表反射率影像和cosi影像上随机选取采样点(共103个)，则这些采样点线性回归方程的相关系数表示地表反射率和cosi的相关性。相关系数越大，相关性越强。地形校正前后各反射率与cosi的相关性如图2和图3所示(以LISS 3第四波段为例)。 

图2是表观反射率与cosi的相关性图，y轴表示表观反射率，x轴表示cosi 

图3是本发明方法得到的地表反射率与cosi的相关性图，y轴表示地表反射率，x轴表示cosi 

从图2和图3可以看出，利用本发明方法进行校正后，地表反射率与cosi之间的相关性(相关系数R²)明显减弱，表明了本发明方法的有效性。 

Claims (1)

1.一种耦合卫星遥感影像大气校正和地形校正过程的地表反射率获取技术，其步骤为：

第一步、对卫星影像进行几何校正；

第二步、辐射定标，将影像像元亮度值转换为星上辐射亮度；

L_λ＝Gain·DN+Bias

L_λ为波段λ的星上辐射亮度，DN为波段λ的像元值，Gain和Bias分别为波段λ的增益和偏置；

第三步、获取大气程辐射

将直方图中频数累加和等于总象素个数0.01％的像素DN值作为暗目标像元值DN_min，程辐射L_p的计算式变为：

L_p＝Gain·DN_min+Bias

第四步、根据陆地卫星遥感影像的经纬度和影像成像时间，获取同一地理范围、同步过境的MODIS气溶胶光学厚度(0.47μm和0.66μm)、总大气水蒸汽含量、地表气压和臭氧浓度影像，并将各影像做几何校正，转换为和陆地卫星遥感影像一致的投影方式和分辨率，且影像行列数一致；

第五步、由总大气水蒸汽含量W计算大气水蒸汽光学厚度；

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$

其中，τ_w是大气水蒸汽光学厚度，a_wλ是水汽吸收系数，w是可降水汽，可降水汽w与总大气水蒸汽含量W在数值上相等，M是相对大气量；

第六步、计算臭氧吸收光学厚度

τ_o＝C_ozone*A_ozone(λ)

τ_o是臭氧吸收光学厚度，C_ozone为MODIS臭氧浓度(单位为Dobson)，A_ozone(λ)为波段λ的臭氧吸收系数；

第七步、由MODIS0.47μm和0.66μm的气溶胶光学厚度，利用下式计算波长指数(α)及大气浑浊度系数(β)

$\mathrm{\α}=-\frac{\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_1\right)\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_2\right)}{\ln {\mathrm{\λ}}_1-{\ln \mathrm{\λ}}_2}$

β＝τ_a(λ₁)λ₁ ^α＝τ_a(λ₂)λ₂ ^α

λ₁和λ₂分别为0.47μm和0.66μm，τ_a(λ₁)和τ_a(λ₂)分别为0.47μm和0.66μm的气溶胶光学厚度；

在得到α和β后，可以利用下式得到任意波长λ的气溶胶光学厚度τ_a(λ)：

τ_a(λ)＝βλ^-a

其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度；

第八步、计算瑞利散射光学厚度；

${\mathrm{\τ}}_r=\frac{p}{1013.25}{0.008569\mathrm{\λ}}^{-4}(1+{0.0133\mathrm{\λ}}^{-2}+{0.00013\mathrm{\λ}}^{-4})$

其中，τ_r为瑞利散射光学厚度，λ是影像各波段的中心波长，μm；p为MODIS的地表气压产品；

第九步、计算太阳照射方向和传感器观测方向的大气透过率；

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z}

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v}

其中，T_z和T_v分别是太阳照射方向和传感器观测方向的大气透过率，τ是大气光学厚度，θ_z是太阳天顶角，θ_v是传感器观测天顶角，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度；

第十步、计算天空散射辐射

E_f＝(E_r+E_a)(3+cos(2S))/4

其中，E_f为天空散射辐射，E_r为瑞利散射分量，E_a为气溶胶散射分量，S为倾斜坡面的坡度角；假定瑞利散射向下的比例是50％，与太阳天顶角无关；同时假定气溶胶散射向下的比例是F_s，F_s是太阳天顶角的函数，则E_r和E_a的表达式如下：

E_r＝0.5×E₀/d²cosθ_zexp{(-τ_o-τ_w)/cosθ_z}T_aa(1-exp{(-0.95τ_r)/cosθ_z})

E_a＝E₀/d²cosθ_zexp{(-τ_o-τ_w)/cosθ_z}T_aaexp{(-1.5τ_r)/cosθ_z}(1-T_as)F_s

其中，E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离，θ_z是太阳天顶角，τ_o、τ_w和τ_r分别是臭氧吸收光学厚度、水汽吸收光学厚度和瑞利散射光学厚度，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，F_s是气溶胶散射向下分量的比例；

日地距离d²(天文单位)按下式计算：

$\mathrm{\Γ}=\frac{360}{365}(\mathrm{dn}-1)$

$d^2=\frac{1}{1.00011+0.034221\cos \mathrm{\Γ}+0.00128\sin \mathrm{\Γ}+0.000719\cos \mathrm{\Γ}+0.00077\sin 2\mathrm{\Γ}}$

其中，dn是儒略日，即影像获取日期距离1月1日的天数；

第十一步、计算地形反射辐射；

E_s＝0.5×(E₀cosθ_ZT_Z+E_r+E_a)×ρ_mean×(1-cosS)

其中E_s为坡面像元接收的周围地形的反射辐射，S为倾斜坡面的坡度角；ρ_mean表示相邻物体的平均反射率，取消除大气影响后的反射率影像的平均值来表示，具体计算式如下：

${\mathrm{\ρ}}_{\mathrm{mean}}=\mathrm{average}\left(\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(E_0\cos {\mathrm{\θ}}_z{T}_z+E_r+E_a)}\right)$

第十二步、利用Dymond公式考虑地形对地表反射的影响，地表反射率的计算公式可以表示为：

$\mathrm{\ρ}=\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(\frac{\mathrm{\Θ}{E}_0{T}_z(\cos i_h+\cos e_h)\cos i}{\cos i+\cos e}+E_f+E_s)}$

i和e分别为直射辐射在坡面上的入射角和出射角，i_h和e_h分别为直射辐射在水平面上的入射角和出射角，入射角i_h与太阳天顶角θ_z相等，则cosi_h＝cosθ_z；

S为倾斜坡面的坡度角，为太阳方位角，A为倾斜坡面的坡向，θ_z为太阳天顶角；

假设卫星传感器是星下点观测，则e_h＝0，cose_h＝1，且由几何关系可以推得：e＝S，则最终地表反射率的计算公式为：

$\mathrm{\ρ}=\frac{\mathrm{\π}(L_{\mathrm{\λ}}-L_p)d^2}{T_v(\frac{\mathrm{\Θ}{E}_0{T}_z(\cos {\mathrm{\θ}}_z+1)\cos i}{\cos i+\cos S}+E_f+E_s)}$

第十三步、对地表反射率反演结果进行验证。