一种Landsat8卫星数据地表反射率反演方法
【技术领域】
[0001] 本发明设及一种从Landsat 8卫星数据反演地表反射率的方法，能够应用在农 业、林业、气象、生态环境W及国防军事等遥感部口。
【背景技术】
[0002] 为了从卫星遥感影像反演得到地表反射率，我们需要对遥感影像进行大气校正。 暗目标法是一种较常用的卫星影像大气校正方法。该方法的基本原理就是在地表均匀、朗 伯面反射、大气性质均一，忽略大气多次散射的前提下，假定待校正的遥感图像上存在反射 率为0的黑暗像元区域，由于大气的影响，使得运些像元的反射率实际不为0,可W认为增 加的运部分福射是由于大气程福射而产生的。利用暗目标像元值计算出程福射，借助于暗 目标法即可实现大气校正。大气透过率是决定暗目标法校正精度的关键因素。在影响大气 透过率的四种因素（大气分子、臭氧含量、气溶胶和可降水汽）中，大气分子和臭氧含量比 较稳定，在空间上差异不大，而气溶胶和可降水汽在空间上变动较大，因此必须提高气溶胶 和可降水汽的估算精度才能提高大气透过率的精度，从而提高最终大气校正的精度。对于 传统的Landsat卫星数据（比如Landsat 5)而言，我们无法从卫星数据本身来得到气溶 胶和可降水汽，必须依赖外部数据源来获取运两个参数，通常通过气象数据或者与Landsat 卫星过境时间接近的Terra/M孤IS数据来间接获取，然而不管是利用气象数据还是版)DIS 数据都具有明显的局限性：气象数据是一种点数据，而遥感数据是一种面数据，气象数据W 点代面的方式会导致较大的误差，而且对于偏远地区或者历史存档卫星数据，获取对应的 气象数据就非常困难;M孤IS数据和Landsat数据在成像时间和空间分辨率上存在较大差 异，两种数据之间的几何配准和投影转换也会带来误差。更重要的是，对于中国大部分地 区，Landsat数据和M孤IS数据之间的地理重叠区域往往非常小（小于=分之一），甚至找 不到与Landsat数据在地理位置上对应的版)DIS数据。W上运些缺陷给传统的Landsat卫 星数据大气校正处理造成了非常大的困难。
[0003] 幸运的是，新型Landsa巧卫星数据的波段设置给基于Landsa巧数据本身进行大 气校正带来了可能。与传统的Landsat系列卫星化andsat5、7)相比，Landsat8在波段的 数量、波段的光谱范围和影像的福射分辨率上进行了改进。Landsa巧携带了两个传感器： 1)OperationalLandImager(0LI)和HiermalInfraredSensor(TIR巧。0LI传感器在可 见光、近红外和短波红外区域接收九个光谱波段的数据，新增加了一个蓝波段（波段1)和 卷云波段（波段9)。TIRS传感器将原来Landsat5、7的热红外波段一分为二，设置成两个 热红外通道度andlO:10. 6-11. 19ym;Band11 :11. 5-12. 51ym)。对于Landsat数据，可 W基于第一波段来获取气溶胶光学厚度，而可降水汽可W利用LandsatS的两个热红外通 道基于劈窗协方差-方差比算法来反演，运样就可W实现基于Landsa巧数据本身来获取地 表反射率。该发明对于实现利用Landsa巧数据业务化地生产地表反射率产品具有重要的 现实意义。
【发明内容】
[0004] 本发明的目的在于提供一种Landsat8卫星数据地表反射率反演方法，该方法主 要基于Landsa巧影像本身来反演，对外部数据源的需求非常低，容易获取，因此具有较强 的实用性。
[0005] 为实现上述目的，本发明提出的方法包括W下步骤：
[0006] 第一步、计算星上福射亮度；
[0007] 第二步、基于Landsat8热红外波段获取可降水汽，计算水蒸汽光学厚度；
[0008] 第=步、基于Landsat8第一波段利用查找表方法获取气溶胶光学厚度；
[0009] 第四步、大气校正和地表反射率反演；
[0010] (4-1)计算大气程福射；
[0011] (4-。计算日地距离；
[0012] (4-3)计算天空光漫射到地表面的光谱福照度；
[0013] (4-4)计算瑞利散射光学厚度和臭氧吸收光学厚度；
[0014] (4-5)计算大气透过率；
[0015] (4-6)计算地表反射率。
【附图说明】
[0016] 图1Landsat8热红外波段大气透过率比值和可降水汽的关系
【具体实施方式】
[0017] 遥感影像地表反射率反演是W福射传输方程为基础的，在假定地表均匀、朗伯面 反射、大气性质均一，忽略大气多次散射的前提下，星上福射亮度和地表反射率的关系如式 1所示：
[001 引
（式 1)
[0019] 式1中：是星上福射亮度，Lp是程福射，Fd是地表接收到的福照度，Ty是传感器 观测方向的大气透过率，S是大气下界的半球反射率，P是地表反射率，d是日地距离。其 中Fd=Eb+Ed。?，Ed。?是由天空光漫射到地表面的光谱福照度，Eb是太阳直射福照度，Eb= Eecos( 0 ,)T,，E。是大气层外相应波长的太阳光谱福照度，可由探测器响应函数计算得到。T, 是太阳照射方向上的大气透过率。0,是太阳天顶角。因为S值很小，通常可W忽略，所W 由式1可W得到地表反射率的计算公式（式2):
[0020]
(式。
[0021] 星上福射亮度Lg。,可W由像元亮度值经福射定标得到 [002引Lsat=MLQcal+AL(式 ^
[002引其中，Lsat是星上福射亮度，ML为波段的增益，AL为波段的偏置，Q。。1为影像DN值，M式PAL从Landsat8头文件获得。
[0024] 程福射Lp由式4计算。
[002引 Lp=MLQCAldark+AL(式 4)
[002引 QCALdgtk是影像中暗目标的亮度值，影像中的暗目标可W选择阴影区域、洁净的水 体或者浓密植被区域。
[0027] 日地距离d(天文单位）根据式5和式6计算：
[0028] l/d2= 1. 000110+0. 0：M221cosr+0. 001280sinr+0. 〇〇〇719cos2r+0. 〇〇〇〇77sin 2r(式W
[0029] r= 2jt化-l)/365 (式 6)
[0030] dn为儒略日（即Landsa巧影像获取日期距离1月1日的天数），如果遇到罔年， 则用366代替式6中的365。
[0031] 太阳照射方向和传感器观测方向的大气透过率T,和Ty是大气光学厚度I的函 数。
[0034] 其中，If、X。、T。和X。分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光 学厚度和大气水蒸汽光学厚度。是传感器观测天顶角，对于Landsa巧数据，0y可W近 似为0。
[0035] 瑞利散射光学厚度Tf和臭氧吸收光学厚度I。相对稳定，可W表达为波长入和 高程h的函数：
[0038] 其中，^是Landsa巧影像各波段的中屯、波长（ym)，h是海拔高度（km)。
[0039] 水蒸汽光学厚度T表达为：
[0040]
(式 11)
[0041] 其中W是可降水汽（cm)，是水汽吸收系数，各波长的直可W查文献得 到[BirdRE，RiordanC.Simplesolarspectralmodelfordirectanddiffuse irr曰di曰neeonhorizont曰1 曰ndtiltedpi曰nes曰tthee曰rth'ssurf曰cefor cloudlessatmospheres.Journalofclimateandappliedmeteorology，1986,25 : 87-97.]。相对大气量M由下式获得：
[004引M=[cos0Z+0. 15 巧3. 885- 0Z) 1.253] 1 (式 12)
[0043] 目Z是太阳天顶角。
[0044] 计算水蒸汽光学厚度勺关键在于获取可降水汽W。
[0045] W基于劈窗协方差-方差比算法来反演，该算法假设在无云条件下，N个相邻像元 区域内（对于Landsat8,N可W取值为20,即窗口大小为20像元*20像元），大气条件和 比福射率不发生改变，仅地表溫度发生改变，W按下式计算：
[0046] W=a(X/Xi)+b(式 13)
[0047]
[0048] 其中，T1为i波段的大气透过率，Ij为j波段的大气透过率，e1为i波段的比 福射率，e,为j波段的比福射率，k表示第k个像元，Ti,k为第k个像元i波段的星上亮度 溫度，T，,k为第k个像元j波段的星上亮度溫度，^为i波段N个像元的平均星上亮度溫度， 为j波段N个像元的平均星上亮度溫度。对于Landsa巧数据，i，j分别为10,11。
[0049] 针对LandsatSTIRS数据，采用式（13)和（14)反演W，需要确定系数a和b，系数 a和b可W通过大气福射传输模型模拟可降水汽与热红外波段大气透过率比值的关系求解 得到。
[0050] 利用M0DTRAN4. 0 大气福射传输模型和TIGR(ThermodynamicInitialGuess Retrieval,TIGR)数据库来模拟可降水汽w与Landsa巧热红外波段大气透过率比值Tii/ T1。之间的关系。TIGR数据库是一个由2311条大气剖面组成的气象数据库；其中每条剖面 数据都包含了从地表到大气层顶部每层的气压、气溫、水汽含量和臭氧含量。TIGR数据库 中包括了 872条热带大气剖面，742条中缔度大气剖面和697条高缔度大气剖面。将TIGR 数据库作为M0DTRAN4. 0模型的输入来模拟可降水汽W与热红外波段大气透过率比值T。/ T1。之间的关系。图1表示基于2311条TIGR大气剖面和M0DTRAN4. 0大气福射传输模型 得到的Landsat8热红外波段大气透过率比值和可降水汽的关系。
[0051] 图1Landsat8热红外波段大气透过率比值和可降水汽的关系
[0052] 如图1所示，Landsa巧数据11波段和10波段大气透过率比值和可降水汽有很好 的相关性。从图1中可W看出，在透过率比值为0.9处存在一个拐点，为了更好地拟合大气 透过率比值和可降水汽之间的关系式，W0. 9为分界点将大气透过率比值分成两段进行拟 合，得到大气透过率比值和可降水汽之间的关系式：
[0055] 由式15和式16可W得到系数a和b。
[0056] Landsat8第10波段和第11波段的星上亮度溫度按下式计算：
[0057] T=Kz/lnQ+Ki/Lsat)(式 17)
[0058] 其中，是星上福射亮度，T是星上亮度溫度，K1和K2为常数，从Landsat8头 文件获得。
[0059] 比福射率e利用NDVI(NormalizedDifferenceVegetationIndex)阔值法来获 取：
[0060]
[006。 其中DNbg"d5和DNbg"d4分别表示Landsa巧第5波段和第4波段影像的DN值。
[0062] 当NDVI<NDVIs时，e=eS，其中NDVIs是纯裸±区域的NDVI，eS是上壤的比 福射率；
[0063]当NDVI>NDVIv时，e=eV，其中NDVIv是纯植被区域的NDVI，eV是植被的比 福射率；
[0064] 当NDVIs《NDVI《NDVIV时，e=es(l-FVC)+evFVC
[0065] FVC是植被覆盖度：
[0066]
(式 19)
[0067] NDVI郝NDVIv可W从图像上选取均质的裸上区域和植被区域来获取。es和eV 通过版)DISUCSB比福射率库和Landsat8TIRS波谱响应函数计算得到。
[0068]如果不考虑地表和大气间的多次反射，天空光漫射到地表面的光谱福照度Ed。?可 W表示为：
[006引Edwn=E^a(式 20)
[0070] 其中，Ef为瑞利散射分量，E。为气溶胶散射分量。
[007引其中T。、1；和Tf分别是臭氧吸收透过率、水汽吸收透过率和瑞利散射透过率，T。。 和T。.分别是气溶胶吸