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

摘要

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

主权项

一种耦合卫星遥感影像大气校正和地形校正过程的地表反射率获取技术，其步骤为：第一步、对卫星影像进行几何校正；第二步、辐射定标，将影像像元亮度值转换为星上辐射亮度；L_λ＝Gain·DN+BiasL_λ为波段λ的星上辐射亮度，DN为波段λ的像元值，Gain和Bias分别为波段λ的增益和偏置；第三步、获取大气程辐射将直方图中频数累加和等于总象素个数0.01％的像素DN值作为暗目标像元值DN_min，程辐射L_p的计算式变为：L_p＝Gain·DN_min+Bias第四步、根据陆地卫星遥感影像的经纬度和影像成像时间，获取同一地理范围、同步过境的MODIS气溶胶光学厚度(0.47μm和0.66μm)、总大气水蒸汽含量、地表气压和臭氧浓度影像，并将各影像做几何校正，转换为和陆地卫星遥感影像一致的投影方式和分辨率，且影像行列数一致；第五步、由总大气水蒸汽含量W计算大气水蒸汽光学厚度；${\tau }_w=\frac{0.2385a_{\mathrm{w\lambda }}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\lambda }}\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的气溶胶光学厚度，利用下式计算波长指数(α)及大气浑浊度系数(β)$\alpha =-\frac{\ln {\tau }_a\left({\lambda }_1\right)\ln {\tau }_a\left({\lambda }_2\right)}{\ln {\lambda }_1-{\ln \lambda }_2}$β＝τ_a(λ₁)λ₁^α＝τ_a(λ₂)λ₂^αλ₁和λ₂分别为0.47μm和0.66μm，τ_a(λ₁)和τ_a(λ₂)分别为0.47μm和0.66μm的气溶胶光学厚度；在得到α和β后，可以利用下式得到任意波长λ的气溶胶光学厚度τ_a(λ)：τ_a(λ)＝βλ^‑a其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度；第八步、计算瑞利散射光学厚度；${\tau }_r=\frac{p}{1013.25}{0.008569\lambda }^{-4}(1+{0.0133\lambda }^{-2}+{0.00013\lambda }^{-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²(天文单位)按下式计算：$\Gamma =\frac{360}{365}(\mathrm{dn}-1)$$d^2=\frac{1}{1.00011+0.034221\cos \Gamma +0.00128\sin \Gamma +0.000719\cos \Gamma +0.00077\sin 2\Gamma }$其中，dn是儒略日，即影像获取日期距离1月1日的天数；第十一步、计算地形反射辐射；E_s＝0.5×(E₀cosθ_ZT_Z+E_r+E_a)×ρ_mean×(1‑cosS)其中E_s为坡面像元接收的周围地形的反射辐射，S为倾斜坡面的坡度角；ρ_mean表示相邻物体的平均反射率，取消除大气影响后的反射率影像的平均值来表示，具体计算式如下：${\rho }_{\mathrm{mean}}=\mathrm{average}\left(\frac{\pi (L_{\lambda }-L_p)d^2}{T_v(E_0\cos {\theta }_z{T}_z+E_r+E_a)}\right)$第十二步、利用Dymond公式考虑地形对地表反射的影响，地表反射率的计算公式可以表示为：$\rho =\frac{\pi (L_{\lambda }-L_p)d^2}{T_v(\frac{\Theta 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，则最终地表反射率的计算公式为：$\rho =\frac{\pi (L_{\lambda }-L_p)d^2}{T_v(\frac{\Theta E_0{T}_z(\cos {\theta }_z+1)\cos i}{\cos i+\cos S}+E_f+E_s)}$第十三步、对地表反射率反演结果进行验证。