[发明专利]一种Landsat8卫星数据地表反射率反演方法

摘要

一种Landsat 8卫星数据地表反射率反演方法，该方法主要基于Landsat8影像本身来反演地表反射率，对外部数据源的需求非常低，容易获取，克服了传统Landsat数据地表反射率反演必须依赖较多外部数据源造成的局限，因此该方法具有较强的实用性，对于实现利用Landsat8数据业务化地生产地表反射率产品具有重要的现实意义。

主权项

1.一种Landsat 8卫星数据地表反射率反演方法，其步骤为：/n第一步、计算星上辐射亮度；/nL_sat＝M_LQ_cal+A_L (式1)/n其中，L_sat是星上辐射亮度，M_L为波段的增益，A_L为波段的偏置，Q_cal为影像DN值，M_L和A_L从Landsat 8头文件获得；/n第二步、计算水蒸汽光学厚度；/n(2-1)获取可降水汽w；/nw＝a(τ_j/τ_i)+b (式2)/n 且/n其中，τ_i为i波段的大气透过率，τ_j为j波段的大气透过率，ε_i为i波段的比辐射率，ε_j为j波段的比辐射率，k表示第k个像元，T_i，k为第k个像元i波段的星上亮度温度，T_j，k为第k个像元j波段的星上亮度温度，为i波段N个像元的平均星上亮度温度，为j波段N个像元的平均星上亮度温度；对于Landsat8数据，i，j分别为10，11；/n系数a和b可以通过式4-5获得，式4-5中的R²为决定系数；/nw＝-18.973(τ₁₁/τ₁₀)+19.13，R²＝0.9663，τ₁₁/τ₁₀＞0.9 (式4)/nw＝-13.412(τ₁₁/τ₁₀)+14.158，R²＝0.9366，τ₁₁/τ₁₀＜0.9 (式5)/nLandsat 8第10波段和第11波段的星上亮度温度按下式计算：/nT＝K₂/ln(1+K₁/L_sat) (式6)/n其中，L_sat是星上辐射亮度，T是星上亮度温度，K₁和K₂为常数，从Landsat 8头文件获得；/n比辐射率ε利用NDVI即Normalized Difference Vegetation Index阈值法来获取：/n /n其中DN_band5和DN_band4分别表示Landsat8第5波段和第4波段影像的DN值；/n当NDVI＜NDVI_s时，ε＝ε_s，其中NDVI_s是纯裸土区域的NDVI，ε_s是土壤的比辐射率；/n当NDVI＞NDVI_v时，ε＝ε_v，其中NDVI_v是纯植被区域的NDVI，ε_v是植被的比辐射率；/n当NDVI_s≤NDVI≤NDVI_v时，ε＝ε_s(1-FVC)+ε_vFVC/nFVC是植被覆盖度：/n /nNDVI_s和NDVI_v可以从图像上选取均质的裸土区域和植被区域来获取；ε_s和ε_v通过MODISUCSB比辐射率库和Landsat 8 TIRS波谱响应函数计算得到；/n(2-2)水蒸汽光学厚度τ_w可表达为：/n /n其中w是可降水汽，单位为cm，a_wλ是水汽吸收系数，相对大气量M由下式获得：/nM＝[cosθ_z+0.15(93.885-θ_z)^-1.253]^-1 (式10)/nθ_z是太阳天顶角；/n第三步、计算气溶胶光学厚度；/n(3-1)计算太阳天顶角和太阳方位角，包括以下步骤：/n(3-1-1)计算太阳时：/n /n其中t是太阳时，单位为小时数，带小数位，t_s是标准时，单位为小时数，带小数位，SM是该标准时对应时区标准经线的经度，L是该像元点的经度，J是儒略日；/n(3-1-2)计算太阳赤纬：/n /n其中δ是太阳赤纬，J是儒略日；/n(3-1-3)计算太阳天顶角：/n /n其中θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬，t是太阳时；/n(3-1-4)计算太阳方位角：/n /n其中是太阳方位角，θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬；/n(3-2)计算Landsat8第一波段的星上反射率：/nρ_sat＝(M_ρQ_cal+A_ρ)/cos(θ_z) (式15)/nρ_sat为星上反射率，M_ρ为增益，A_ρ为偏置，这两个参数从Landsat8头文件获得；Q_cal为影像DN值，θ_z是太阳天顶角；/n(3-3)计算气溶胶光学厚度；/n /n其中，ρ_sat为星上反射率，ρ₀为大气的路径辐射项等效反射率，μ_s和μ_v分别为太阳天顶角θ_z与观测天顶角θ_v的余弦，为相对方位角；r为地表反射率，S为大气整层向下的半球反射率，T(μ_s)和T(μ_v)分别为太阳照射方向和传感器观测方向的大气透过率；/n利用6S模型设定不同的条件构建查找表；设定参数包括：12个太阳天顶角，取值范围为0-66度，间隔6度，16个太阳与卫星之间的相对方位角，取值范围为0-180度，间隔12度，大气气溶胶模式参数假设为大陆性气溶胶，并设立20个大气气溶胶光学厚度值，取值范围为0-2，间隔0.1，波段自定义为Landsat 8 OLI传感器第一波段的响应函数；/n根据计算得到的太阳天顶角和相对方位角，对查找表进行线性插值，得到不同光学厚度下的大气参数S、ρ₀和T(μ_s)T(μ_v)，代入式(16)，同时将Landsat8卫星过境前的MODIS 8天合成地表反射率产品，即MOD09产品，取第三波段，代入式(16)获得不同气溶胶光学厚度下的星上反射率；然后利用式15得到的Landsat8第一波段的星上反射率进行线性插值，得到气溶胶光学厚度；/n第四步、大气校正和地表反射率反演；/n(4-1)计算大气程辐射L_p；/nL_p＝M_LQCAL_dark+A_L (式17)/nQCAL_dark是影像中暗目标的亮度值；/n(4-2)计算日地距离d；/n1/d²＝1.000110+0.034221cosΓ+0.001280sinΓ+0.000719cos2Γ+0.000077sin2Γ(式18)/nΓ＝2π(d_n-1)/365 (式19)/nd_n为儒略日；/n(4-3)计算天空光漫射到地表面的光谱辐照度E_down；/nE_down＝E_r+E_a (式20)/n其中，E_r为瑞利散射分量，E_a为气溶胶散射分量；/nE_r＝E₀/d²cos(θ_z)T_oT_wT_aa(1-T_r^0.95)0.5 (式21)/n /n其中，E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离，θ_z是太阳天顶角，T_o、T_w和T_r分别是臭氧吸收透过率、水汽吸收透过率和瑞利散射透过率，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，F_s是气溶胶散射向下分量的比例；/n(4-4)计算瑞利散射光学厚度τ_r和臭氧吸收光学厚度τ_o/nτ_r＝[0.0088λ^{(-4.15+0.2λ)}][exp(-0.1188h-0.00116h²)] (式23)/n /n其中，λ是Landsat8影像各波段的中心波长，单位为μm；h是海拔高度，单位为km；/n(4-5)计算大气透过率；/nT_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z} (式25)/nT_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v} (式26)/n其中，T_z和T_v分别是太阳照射方向和传感器观测方向的大气透过率，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度，θ_z是太阳天顶角，θ_v是传感器观测天顶角；/n(4-6)计算地表反射率；/n /n其中，ρ是地表反射率，L_sat是星上辐射亮度，L_p是程辐射，d是日地距离，T_v是传感器观测方向的大气透过率，E₀是大气层外相应波长的太阳光谱辐照度，θ_z是太阳天顶角，T_z是太阳照射方向上的大气透过率，E_down是天空光漫射到地表面的光谱辐照度。/n