CN105183989B - 一种Landsat8卫星数据地表反射率反演方法 - Google Patents

Abstract

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

Description

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

技术领域

本发明涉及一种从Landsat 8卫星数据反演地表反射率的方法，能够应用在农业、林业、气象、生态环境以及国防军事等遥感部门。

背景技术

为了从卫星遥感影像反演得到地表反射率，我们需要对遥感影像进行大气校正。暗目标法是一种较常用的卫星影像大气校正方法。该方法的基本原理就是在地表均匀、朗伯面反射、大气性质均一，忽略大气多次散射的前提下，假定待校正的遥感图像上存在反射率为0的黑暗像元区域，由于大气的影响，使得这些像元的反射率实际不为0，可以认为增加的这部分辐射是由于大气程辐射而产生的。利用暗目标像元值计算出程辐射，借助于暗目标法即可实现大气校正。大气透过率是决定暗目标法校正精度的关键因素。在影响大气透过率的四种因素(大气分子、臭氧含量、气溶胶和可降水汽)中，大气分子和臭氧含量比较稳定，在空间上差异不大，而气溶胶和可降水汽在空间上变动较大，因此必须提高气溶胶和可降水汽的估算精度才能提高大气透过率的精度，从而提高最终大气校正的精度。对于传统的Landsat卫星数据(比如Landsat 5)而言，我们无法从卫星数据本身来得到气溶胶和可降水汽，必须依赖外部数据源来获取这两个参数，通常通过气象数据或者与Landsat卫星过境时间接近的Terra/MODIS数据来间接获取，然而不管是利用气象数据还是MODIS数据都具有明显的局限性：气象数据是一种点数据，而遥感数据是一种面数据，气象数据以点代面的方式会导致较大的误差，而且对于偏远地区或者历史存档卫星数据，获取对应的气象数据就非常困难；MODIS数据和Landsat数据在成像时间和空间分辨率上存在较大差异，两种数据之间的几何配准和投影转换也会带来误差。更重要的是，对于中国大部分地区，Landsat数据和MODIS数据之间的地理重叠区域往往非常小(小于三分之一)，甚至找不到与Landsat数据在地理位置上对应的MODIS数据。以上这些缺陷给传统的Landsat卫星数据大气校正处理造成了非常大的困难。

幸运的是，新型Landsat8卫星数据的波段设置给基于Landsat8数据本身进行大气校正带来了可能。与传统的Landsat系列卫星(Landsat5、7)相比，Landsat 8在波段的数量、波段的光谱范围和影像的辐射分辨率上进行了改进。Landsat8携带了两个传感器：1)Operational Land Imager(OLI)和Thermal Infrared Sensor(TIRS)。OLI传感器在可见光、近红外和短波红外区域接收九个光谱波段的数据，新增加了一个蓝波段(波段1)和卷云波段(波段9)。TIRS传感器将原来Landsat5、7的热红外波段一分为二，设置成两个热红外通道(Band10：10.6-11.19μm；Band 11：11.5-12.51μm)。对于Landsat8数据，可以基于第一波段来获取气溶胶光学厚度，而可降水汽可以利用Landsat8的两个热红外通道基于劈窗协方差-方差比算法来反演，这样就可以实现基于Landsat8数据本身来获取地表反射率。该发明对于实现利用Landsat8数据业务化地生产地表反射率产品具有重要的现实意义。

发明内容

本发明的目的在于提供一种Landsat 8卫星数据地表反射率反演方法，该方法主要基于Landsat8影像本身来反演，对外部数据源的需求非常低，容易获取，因此具有较强的实用性。

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

第一步、计算星上辐射亮度；

第二步、基于Landsat 8热红外波段获取可降水汽，计算水蒸汽光学厚度；

第三步、基于Landsat 8第一波段利用查找表方法获取气溶胶光学厚度；

第四步、大气校正和地表反射率反演；

(4-1)计算大气程辐射；

(4-2)计算日地距离；

(4-3)计算天空光漫射到地表面的光谱辐照度；

(4-4)计算瑞利散射光学厚度和臭氧吸收光学厚度；

(4-5)计算大气透过率；

(4-6)计算地表反射率。

附图说明

图1Landsat 8热红外波段大气透过率比值和可降水汽的关系

具体实施方式

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

(式1)

式1中：L_sat是星上辐射亮度，L_p是程辐射，F_d是地表接收到的辐照度，T_v是传感器观测方向的大气透过率，s是大气下界的半球反射率，ρ是地表反射率，d是日地距离。其中F_d＝E_b+E_down，E_down是由天空光漫射到地表面的光谱辐照度，E_b是太阳直射辐照度，E_b＝E₀cos(θ_z)T_z，E₀是大气层外相应波长的太阳光谱辐照度，可由探测器响应函数计算得到。T_z是太阳照射方向上的大气透过率。θ_z是太阳天顶角。因为s值很小，通常可以忽略，所以由式1可以得到地表反射率的计算公式(式2)：

(式2)

星上辐射亮度L_sat可以由像元亮度值经辐射定标得到

L_sat＝M_LQ_cal+A_L (式3)

其中，L_sat是星上辐射亮度，M_L为波段的增益，A_L为波段的偏置，Q_cal为影像DN值，M_L和A_L从Landsat 8头文件获得。

程辐射L_p由式4计算。

L_p＝M_LQCAL_dark+A_L (式4)

QCAL_dark是影像中暗目标的亮度值，影像中的暗目标可以选择阴影区域、洁净的水体或者浓密植被区域。

日地距离d(天文单位)根据式5和式6计算：

1/d²＝1.000110+0.034221cosΓ+0.001280sinΓ+0.000719cos2Γ+0.000077sin2Γ (式5)

Γ＝2π(d_n-1)/365 (式6)

d_n为儒略日(即Landsat8影像获取日期距离1月1日的天数)，如果遇到闰年，则用366代替式6中的365。

太阳照射方向和传感器观测方向的大气透过率T_z和T_v是大气光学厚度τ的函数。

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

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

其中，τ_r、τ_a、τ₀和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度。θ_v是传感器观测天顶角，对于Landsat8数据，θ_v可以近似为0。

瑞利散射光学厚度τ_r和臭氧吸收光学厚度τ_o相对稳定，可以表达为波长λ和高程h的函数：

τ_r＝[0.0088λ^{(-4.15+0.2λ)}][exp(-0.1188h-0.00116h²)] (式9)

(式10)

其中，λ是Landsat8影像各波段的中心波长(μm)，h是海拔高度(km)。

水蒸汽光学厚度τ_w可表达为：

(式11)

其中w是可降水汽(cm)，a_wλ是水汽吸收系数，各波长的a_wλ值可以查文献得到[BirdR E，Riordan C.Simple solar spectral model for direct and diffuse irradianceon horizontal and tilted planes at the earth′s surface for cloudlessatmospheres.Journal of climate and applied meteorology，1986，25：87-97.]。相对大气量M由下式获得：

M＝[cosθ_z+0.15(93.885-θ_z)^-1.253]^-1 (式12)

θ_z是太阳天顶角。

计算水蒸汽光学厚度τ_w的关键在于获取可降水汽w。

w基于劈窗协方差-方差比算法来反演，该算法假设在无云条件下，N个相邻像元区域内(对于Landsat 8，N可以取值为20，即窗口大小为20像元*20像元)，大气条件和比辐射率不发生改变，仅地表温度发生改变，w按下式计算：

w＝a(τ_j/τ_i)+b (式13)

且 (式14)

其中，τ_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。

针对Landsat8 TIRS数据，采用式(13)和(14)反演w，需要确定系数a和b，系数a和b可以通过大气辐射传输模型模拟可降水汽与热红外波段大气透过率比值的关系求解得到。

利用MODTRAN4.0大气辐射传输模型和TIGR(Thermodynamic Initial GuessRetrieval，TIGR)数据库来模拟可降水汽w与Landsat8热红外波段大气透过率比值τ₁₁/τ₁₀之间的关系。TIGR数据库是一个由2311条大气剖面组成的气象数据库；其中每条剖面数据都包含了从地表到大气层顶部每层的气压、气温、水汽含量和臭氧含量。TIGR数据库中包括了872条热带大气剖面，742条中纬度大气剖面和697条高纬度大气剖面。将TIGR数据库作为MODTRAN4.0模型的输入来模拟可降水汽w与热红外波段大气透过率比值τ₁₁/τ₁₀之间的关系。图1表示基于2311条TIGR大气剖面和MODTRAN4.0大气辐射传输模型得到的Landsat 8热红外波段大气透过率比值和可降水汽的关系。

图1 Landsat 8热红外波段大气透过率比值和可降水汽的关系

如图1所示，Landsat8数据11波段和10波段大气透过率比值和可降水汽有很好的相关性。从图1中可以看出，在透过率比值为0.9处存在一个拐点，为了更好地拟合大气透过率比值和可降水汽之间的关系式，以0.9为分界点将大气透过率比值分成两段进行拟合，得到大气透过率比值和可降水汽之间的关系式：

w＝-18.973(τ₁₁/τ₁₀)+19.13 R²＝0.9663，τ₁₁/τ₁₀＞0.9 (式15)

w＝-13.412(τ₁₁/τ₁₀)+14.158 R²＝0.9366，τ₁₁/τ₁₀＜0.9 (式16)

由式15和式16可以得到系数a和b。

Landsat 8第10波段和第11波段的星上亮度温度按下式计算：

T＝K₂/ln(1+K₁/L_sat) (式17)

其中，L_sat是星上辐射亮度，T是星上亮度温度，K1和K2为常数，从Landsat 8头文件获得。

比辐射率ε利用NDVI(Normalized Difference Vegetation Index)阈值法来获取：

(式18)

其中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)+ε_vFVC

FVC是植被覆盖度：

(式19)

NDVI_s和NDVI_v可以从图像上选取均质的裸土区域和植被区域来获取。ε_s和ε_v通过MODIS UCSB比辐射率库和Landsat 8 TIRS波谱响应函数计算得到。

如果不考虑地表和大气间的多次反射，天空光漫射到地表面的光谱辐照度E_down可以表示为：

E_down＝E_r+E_a (式20)

其中，E_r为瑞利散射分量，E_a为气溶胶散射分量。

(式21)

(式22)

其中T_o、T_w和T_r分别是臭氧吸收透过率、水汽吸收透过率和瑞利散射透过率，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，两者之积为气溶胶透过率，F_s是气溶胶散射向下分量的比例，各参数的具体计算参见文献[BirdR E，Riordan C.Simple solarspectral model for direct and diffuse irradiance on horizontal and tiltedplanes at the earth′s surface for cloudless atmospheres.Journal of climateand applied meteorology，1986，25：87-97.]。

气溶胶光学厚度可以利用浓密植被法来反演，然而这种方法只适用于影像中存在浓密植被的情况，对于影像中不存在浓密植被的情况(比如沙漠地区)，该方法不适用。本发明基于Landsat8第1波段来反演气溶胶光学厚度。

气溶胶光学厚度反演基于以下公式：

(式23)

其中，L_sat为星上辐射亮度，μ_s和μ_v分别为太阳天顶角θ_z与观测天顶角θ_v的余弦，为相对方位角；L₀为观测方向的路径辐射项，r为地表反射率，S为大气整层向下的半球反射率，T为大气透过率，T(μ_s)和T(μ_v)分别为太阳照射方向和传感器观测方向的大气透过率，F₀为大气层顶太阳光的辐射通量密度。

利用垂直入射的太阳辐射量μ_sF₀对上式进行归一化可得星上反射率ρ_sat

(式24)

其中，ρ₀为大气的路径辐射项等效反射率。

上式中，S、ρ₀和T(μ_s)T(μ_v)这3个参数代表大气的状况，因此可以从中获取大气气溶胶光学厚度。在Landsat8第一波段，卫星传感器接收到的辐射中大气反射的贡献较大，地表反射的贡献微弱，假定同期的地表反射率不变，将晴好天的地表反射率(本发明利用Landsat8卫星过境前的MODIS 8天合成地表反射率产品，MOD09产品(第三波段))代入上式反演气溶胶光学厚度。在具体反演中是用辐射传输模型构建气溶胶光学厚度和太阳天顶角、观测天顶角和相对方位角之间的查找表，然后利用MODIS的8天合成地表反射率产品去除地表反射贡献，得到气溶胶光学厚度。

由于Landsat8影像可以看做是星下点成像，整景影像范围内观测天顶角可以近似为0，这样就无需考虑观测天顶角和观测方位角。下面介绍如何逐像元地计算太阳天顶角和太阳方位角，包括以下步骤：

(1)计算太阳时：

(式25)

其中t是太阳时(小时数，带小数位)，t_s是标准时(小时数，带小数位)，SM是该时区标准经线的经度(弧度)，L是该像元点的经度(弧度)，J是儒略日。

(2)计算太阳赤纬：

(式26)

其中δ是太阳赤纬(弧度)，J是儒略日。

(3)计算太阳天顶角

(式27)

其中θ_z是太阳天顶角(弧度)，l是该像元点的纬度(弧度)，δ是太阳赤纬(弧度)，t是太阳时。

(4)计算太阳方位角

(式28)

其中是太阳方位角，θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬。

按式(29)计算Landsat8第一波段的星上反射率：

ρ_sat＝(M_ρQ_cal+A_ρ)/cos(θ_z) (式29)

ρ_sat为星上反射率，M_ρ为增益，A_ρ为偏置，这两个参数从Landsat8头文件获得；Q_cal为影像DN值，θ_z是太阳天顶角。

查找表利用6S模型设定不同的条件进行辐射传输计算得到。设定参数包括：12个太阳天顶角(0-66度，间隔6度)，16个太阳与卫星之间的相对方位角(0-180度，间隔12度)，大气气溶胶模式参数假设为大陆性气溶胶，并设立20个大气气溶胶光学厚度值(0-2，间隔0.1)，波段自定义为Landsat 8 OLI传感器第一波段的响应函数。

根据计算得到的观测几何(太阳天顶角和相对方位角)，对查找表进行线性插值，得到不同光学厚度下的大气参数S、ρ₀和T(μ_s)T(μ_v)，代入式(24)，同时将Landsat8卫星过境前的MODIS 8天合成地表反射率产品(MOD09产品，第三波段)代入式(24)获得不同气溶胶光学厚度下的星上反射率；然后利用式29得到的Landsat8第一波段的星上反射率进行线性插值，得到气溶胶光学厚度。

Claims (1)

1.一种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)

且

其中，τ_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获得，式4-5中的R²为决定系数；

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是星上亮度温度，K₁和K₂为常数，从Landsat 8头文件获得；

比辐射率ε利用NDVI即Normalized Difference Vegetation Index阈值法来获取：

其中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)+ε_vFVC

FVC是植被覆盖度：

NDVI_s和NDVI_v可以从图像上选取均质的裸土区域和植被区域来获取；ε_s和ε_v通过MODISUCSB比辐射率库和Landsat 8 TIRS波谱响应函数计算得到；

(2-2)水蒸汽光学厚度τ_w可表达为：

其中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是标准时，单位为小时数，带小数位，SM是该标准时对应时区标准经线的经度，L是该像元点的经度，J是儒略日；

(3-1-2)计算太阳赤纬：

其中δ是太阳赤纬，J是儒略日；

(3-1-3)计算太阳天顶角：

其中θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬，t是太阳时；

(3-1-4)计算太阳方位角：

其中是太阳方位角，θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬；

(3-2)计算Landsat8第一波段的星上反射率：

ρ_sat＝(M_ρQ_cal+A_ρ)/cos(θ_z) (式15)

ρ_sat为星上反射率，M_ρ为增益，A_ρ为偏置，这两个参数从Landsat8头文件获得；Q_cal为影像DN值，θ_z是太阳天顶角；

(3-3)计算气溶胶光学厚度；

其中，ρ_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卫星过境前的MODIS 8天合成地表反射率产品，即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₀/d²cos(θ_z)T_oT_wT_aa(1-T_r ^0.95)0.5 (式21)

其中，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)

其中，λ是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θ_v} (式26)

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

(4-6)计算地表反射率；

其中，ρ是地表反射率，L_sat是星上辐射亮度，L_p是程辐射，d是日地距离，T_v是传感器观测方向的大气透过率，E₀是大气层外相应波长的太阳光谱辐照度，θ_z是太阳天顶角，T_z是太阳照射方向上的大气透过率，E_down是天空光漫射到地表面的光谱辐照度。