一种Landsat8卫星数据地表反射率反演方法

摘要

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

主权项

一种Landsat 8卫星数据地表反射率反演方法，其步骤为：第一步、计算星上辐射亮度；L_sat＝M_LQ_cal+A_L   (式1)其中，L_sat是星上辐射亮度，M_L为波段的增益，A_L为波段的偏置，Q_cal为影像DN值，M_L和A_L从Landsat 8头文件获得；第二步、计算水蒸汽光学厚度；(2‑1)获取可降水汽w；w＝a(τ_j/τ_i)+b   (式2)$\frac{{\tau }_j}{{\tau }_i}=\frac{{ϵ}_i}{{ϵ}_j}{R}_{j,i}$且$R_{j,i}=\frac{\underset{k=1}{\overset{N}{\Sigma }}(T_{i,k}-{\bar{T}}_i)(T_{j,k}-{\bar{T}}_j)}{\underset{k=1}{\overset{N}{\Sigma }}{(T_{j,k}-{\bar{T}}_i)}^2}$   (式3)其中，τ_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；系数a和b可以通过式4‑5获得：w＝‑18.973(τ₁₁/τ₁₀)+19.13 R²＝0.9663，τ₁₁/τ₁₀＞0.9   (式4)w＝‑13.412(τ₁₁/τ₁₀)+14.158 R²＝0.9366，τ₁₁/τ₁₀＜0.9   (式5)Landsat 8第10波段和第11波段的星上亮度温度按下式计算：T＝K₂/ln(1+K₁/L_sat)   (式6)其中，L_sat是星上辐射亮度，T是星上亮度温度，K1和K2为常数，从Landsat 8头文件获得；比辐射率ε利用NDVI(Normalized Difference Vegetation Index)阈值法来获取：$NDVI=\frac{{\mathrm{DN}}_{band5}-{\mathrm{DN}}_{band4}}{{\mathrm{DN}}_{band5}+{\mathrm{DN}}_{band4}}$   (式7)其中DN_band5和DN_band4分别表示Landsat8第5波段和第4波段影像的DN值；当NDVI＜NDVI_s时，ε＝ε_s，其中NDVI_s是纯裸土区域的NDVI，ε_s是土壤的比辐射率；当NDVI＞NDVI_v时，ε＝ε_v，其中NDVI_v是纯植被区域的NDVI，ε_v是植被的比辐射率；当NDVI_s≤NDVI≤NDVI_v时，ε＝ε_s(1‑FVC)+ε_vFVCFVC是植被覆盖度：$FVC={\left[\frac{NDVI-{\mathrm{NDVI}}_s}{{\mathrm{NDVI}}_v-{\mathrm{NDVI}}_s}\right]}^2$   (式8)NDVI_s和NDVI_v可以从图像上选取均质的裸土区域和植被区域来获取；ε_s和ε_v通过MODIS UCSB比辐射率库和Landsat 8TIRS波谱响应函数计算得到；(2‑2)水蒸汽光学厚度τ_w可表达为：${\tau }_w=\frac{0.2385a_{w\lambda }wM}{{(1+20.07a_{w\lambda }wM)}^{0.45}}$   (式9)其中w是可降水汽(cm)，a_wλ是水汽吸收系数，相对大气量M由下式获得：M＝[cosθ_z+0.15(93.885‑θ_z)^‑1.253]^‑1   (式10)θ_z是太阳天顶角；第三步、计算气溶胶光学厚度；(3‑1)计算太阳天顶角和太阳方位角，包括以下步骤：(3‑1‑1)计算太阳时：$t=t_s+0.170\sin \left(\frac{4\pi (J-80)}{373}\right)-0.129\sin \left(\frac{2\pi (J-8)}{355}\right)+\frac{12(SM-L)}{\pi }$   (式11)其中t是太阳时(小时数，带小数位)，t_s是标准时(小时数，带小数位)，SM是该时区标准经线的经度(弧度)，L是该像元点的经度(弧度)，J是儒略日；(3‑1‑2)计算太阳赤纬：$\delta =0.4093sin\left(\frac{2\pi (J-81)}{368}\right)$   (式12)其中δ是太阳赤纬(弧度)，J是儒略日；(3‑1‑3)计算太阳天顶角：${\theta }_z=\frac{\pi }{2}-\arcsin (\sin lsin\delta +\cos lcos\delta cos\frac{\pi t}{12})$   (式13)其中θ_z是太阳天顶角(弧度)，l是该像元点的纬度(弧度)，δ是太阳赤纬(弧度)，t是太阳时；(3‑1‑4)计算太阳方位角：(式14)其中是太阳方位角，θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬；(3‑2)计算Landsat8第一波段的星上反射率：ρ_sat＝(M_ρQ_cal+A_ρ)/cos(θ_z)   (式15)ρ_sat为星上反射率，M_ρ为增益，A_p为偏置，这两个参数从Landsat8头文件获得；Q_cal为影像DN值，θ_z是太阳天顶角；(3‑3)计算气溶胶光学厚度；(式16)其中，ρ_sat为星上反射率，ρ₀为大气的路径辐射项等效反射率，μ_s和μ_v分别为太阳天顶角θ_z与观测天顶角θ_v的余弦，为相对方位角；r为地表反射率，S为大气整层向下的半球反射率，T(μ_s)和T(μ_v)分别为太阳照射方向和传感器观测方向的大气透过率；利用6S模型设定不同的条件构建查找表；设定参数包括：12个太阳天顶角(0‑66度，间隔6度)，16个太阳与卫星之间的相对方位角(0‑180度，间隔12度)，大气气溶胶模式参数假设为大陆性气溶胶，并设立20个大气气溶胶光学厚度值(0‑2，间隔0.1)，波段自定义为Landsat 8 OLI传感器第一波段的响应函数；根据计算得到的观测几何(太阳天项角和相对方位角)，对查找表进行线性插值，得到不同光学厚度下的大气参数S、ρ₀和T(μ_s)T(μ_v)，代入式(16)，同时将Landsat8卫星过境前的MODIS8天合成地表反射率产品(MOD09产品，第三波段)代入式(16)获得不同气溶胶光学厚度下的星上反射率；然后利用式15得到的Landsat8第一波段的星上反射率进行线性插值，得到气溶胶光学厚度；第四步、大气校正和地表反射率反演；(4‑1)计算大气程辐射L_p；L_p＝M_LQCAL_dark+A_L   (式17)QCAL_dark是影像中暗目标的亮度值；(4‑2)计算日地距离d；1/d²＝1.000110+0.034221cosΓ+0.001280sinΓ+0.000719cos2Γ+0.000077sin2Γ   (式18)Γ＝2π(d_n‑1)/365   (式19)d_n为儒略日；(4‑3)计算天空光漫射到地表面的光谱辐照度E_down；E_down＝E_r+E_a   (式20)其中，E_r为瑞利散射分量，E_a为气溶胶散射分量；$E_r=E_0/d^{2}cos\left({\theta }_2\right)T_o{T}_w{T}_{aa}(1-T_r^{0.95})0.5$   (式21)$E_a=E_0/d^{2}cos\left({\theta }_z\right)T_o{T}_w{T}_{aa}{T}_r^{1.5}(1-T_{as})F_s$   (式22)其中，E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离，θ_z是太阳天顶角，T_o、T_w和T_r分别是臭氧吸收透过率、水汽吸收透过率和瑞利散射透过率，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，F_s是气溶胶散射向下分量的比例；(4‑4)计算瑞利散射光学厚度τ_r和臭氧吸收光学厚度τ_oτ_r＝[0.0088λ^{(‑4.15+0.2λ)}][exp(‑0.1188h‑0.00116h²)]   (式23)${\tau }_o=0.03(1-\frac{1.0183}{1+0.0183\exp (h/5)})\exp [-277{(\lambda -0.6)}^2]$   (式24)其中，λ是Landsat8影像各波段的中心波长，μm；h是海拔高度，km；(4‑5)计算大气透过率；T_z＝exp(‑τ/cosτ_z)＝exp{(‑τ_r‑τ_a‑τ_o‑τ_w)/cosθ_z}   (式25)T_v＝exp(‑τ/cosθ_v)＝exp{(‑τ_r‑τ_a‑τ_o‑τ_w)/cosθ_y}   (式25)其中，T_z和T_v分别是太阳照射方向和传感器观测方向的大气透过率，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度，θ_z是太阳天顶角，θ_v是传感器观测天顶角；(4‑6)计算地表反射率；$\rho =\frac{\pi (L_{sat}-L_p)d^2}{T_v\left(E_{0}cos\right({\theta }_z)T_z+E_{down})}$   (式27)其中，ρ是地表反射率，L_sat是星上辐射亮度，L_p是程辐射，d是日地距离，T_v是传感器观测方向的大气透过率，E₀是大气层外相应波长的太阳光谱辐照度，θ_z是太阳天顶角，T_z是太阳照射方向上的大气透过率，E_down是天空光漫射到地表面的光谱辐照度。