一种Landsat8卫星数据地表反射率反演方法，该方法主要基于Landsat8影像本身来反演地表反射率，对外部数据源的需求非常低，容易获取，克服了传统Landsat数据地表反射率反演必须依赖较多外部数据源造成的局限，因此该方法具有较强的实用性，对于实现利用Landsat8数据业务化地生产地表反射率产品具有重要的现实意义。
【技术实现步骤摘要】
本专利技术设及一种从Landsat 8卫星数据反演地表反射率的方法，能够应用在农 业、林业、气象、生态环境W及国防军事等遥感部口。
技术介绍
为了从卫星遥感影像反演得到地表反射率，我们需要对遥感影像进行大气校正。 暗目标法是一种较常用的卫星影像大气校正方法。该方法的基本原理就是在地表均匀、朗 伯面反射、大气性质均一，忽略大气多次散射的前提下，假定待校正的遥感图像上存在反射 率为0的黑暗像元区域，由于大气的影响，使得运些像元的反射率实际不为0,可W认为增 加的运部分福射是由于大气程福射而产生的。利用暗目标像元值计算出程福射，借助于暗 目标法即可实现大气校正。大气透过率是决定暗目标法校正精度的关键因素。在影响大气 透过率的四种因素（大气分子、臭氧含量、气溶胶和可降水汽）中，大气分子和臭氧含量比 较稳定，在空间上差异不大，而气溶胶和可降水汽在空间上变动较大，因此必须提高气溶胶 和可降水汽的估算精度才能提高大气透过率的精度，从而提高最终大气校正的精度。对于 传统的Landsat卫星数据（比如Landsat 5)而言，我们无法从卫星数据本身来得到气溶 胶和可降水汽，必须依赖外部数据源来获取运两个参数，通常通过气象数据或者与Landsat 卫星过境时间接近的Terra/M孤IS数据来间接获取，然而不管是利用气象数据还是版)DIS 数据都具有明显的局限性：气象数据是一种点数据，而遥感数据是一种面数据，气象数据W 点代面的方式会导致较大的误差，而且对于偏远地区或者历史存档卫星数据，获取对应的 气象数据就非常困难;M孤IS数据和Landsat数据在成像时间和空间分辨率上存在较大差 异，两种数据之间的几何配准和投影转换也会带来误差。更重要的是，对于中国大部分地 区，Landsat数据和M孤IS数据之间的地理重叠区域往往非常小（小于=分之一），甚至找 不到与Landsat数据在地理位置上对应的版)DIS数据。W上运些缺陷给传统的Landsat卫 星数据大气校正处理造成了非常大的困难。 幸运的是，新型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巧数据业务化地生产地表反射率产品具有重要的 现实意义。
技术实现思路
本专利技术的目的在于提供，该方法主 要基于Landsa巧影像本身来反演，对外部数据源的需求非常低，容易获取，因此具有较强 的实用性。 为实现上述目的，本专利技术提出的方法包括W下步骤： 第一步、计算星上福射亮度； 第二步、基于Landsat8热红外波段获取可降水汽，计算水蒸汽光学厚度； 第=步、基于Landsat8第一波段利用查找表方法获取气溶胶光学厚度； 第四步、大气校正和地表反射率反演； (4-1)计算大气程福射； (4-。计算日地距离； (4-3)计算天空光漫射到地表面的光谱福照度； (4-4)计算瑞利散射光学厚度和臭氧吸收光学厚度； (4-5)计算大气透过率； (4-6)计算地表反射率。【附图说明】 图1Landsat8热红外波段大气透过率比值和可降水汽的关系【具体实施方式】 遥感影像地表反射率反演是W福射传输方程为基础的，在假定地表均匀、朗伯面 反射、大气性质均一，忽略大气多次散射的前提下，星上福射亮度和地表反射率的关系如式 1所示：[001 引（式 1) 式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):(式。 星上福射亮度Lg。,可W由像元亮度值经福射定标得到 [002引Lsat=MLQcal+AL(式 ^[002引其中，Lsat是星上福射亮度，ML为波段的增益，AL为波段的偏置，Q。。1为影像DN值，M式PAL从Landsat8头文件获得。 程福射Lp由式4计算。[002引 Lp=MLQCAldark+AL(式 4)[002引 QCALdgtk是影像中暗目标的亮度值，影像中的暗目标可W选择阴影区域、洁净的水 体或者浓密植被区域。 日地距离d(天文单位）根据式5和式6计算： l/d2= 1. 000110+0. 0：M221cosr+0. 001280sinr+0. 〇〇〇719cos2r+0. 〇〇〇〇77sin 2r(式W r= 2jt化-l)/365 (式 6) dn为儒略日（即Landsa巧影像获取日期距离1月1日的天数），如果遇到罔年， 则用366代替式6中的365。 太阳照射方向和传感器观测方向的大气透过率T,和Ty是大气光学厚度I的函 数。 其中，If、X。、T。和X。分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光 学厚度和大气水蒸汽光学厚度。是传感器观测天顶角，对于Landsa巧数据，0y可W近 似为0。 瑞利散射光学厚度Tf和臭氧吸收光学厚度I。相对稳定，可W表达为波长入和 高程h的函数： 其中，^是Landsa巧影像各波段的中屯、波长（ym)，h是海拔高度（km)。 水蒸汽光学厚度T表达为：(式 11) 其中W是可降水汽（cm)，是水汽吸收系数，各波长的直可W查文献得 到。相对大气量M由下式获得： 1 (式 12) 目Z是太阳天顶角。 计算水蒸汽光学厚度勺关键在于获取可降水汽W。 W基于劈窗协方差-方差比算法来反演，该算法假设在无云条件下，N个相邻像元 区域内（对于Landsat8,N可W取值为20,即窗口大小为20像元*20像元），大气条件和 比福射率不发生改变，仅地表溫度发生改变，W按下式计算： W=a(X/Xi)+b(式 13) 其中，T1为i波段的大气透过率，Ij为j波段的大气透过率，e1为i波段的比 福射率，e,为j波段的比福射率，k表示第k个像元，Ti,k为第k个像元i波段的星上亮度 溫度，T，,k为第k个像元j波段的星上亮度溫度，^为i波段N个像元的平均星上亮度溫度， 为j波段本文档来自技高网...
【技术保护点】
一种Landsat 8卫星数据地表反射率反演方法，其步骤为：第一步、计算星上辐射亮度；Lsat＝MLQcal+AL (式1)其中，Lsat是星上辐射亮度，ML为波段的增益，AL为波段的偏置，Qcal为影像DN值，ML和AL从Landsat 8头文件获得；第二步、计算水蒸汽光学厚度；(2‑1)获取可降水汽w；w＝a(τj/τi)+b (式2)τjτi=ϵiϵjRj,i]]>且Rj,i=Σk=1N(Ti,k-T‾i)(Tj,k-T‾j)Σk=1N(Tj,k-T‾i)2]]> (式3)其中，τi为i波段的大气透过率，τj为j波段的大气透过率，εi为i波段的比辐射率，εj为j波段的比辐射率，k表示第k个像元，Ti，k为第k个像元i波段的星上亮度温度，Tj，k为第k个像元j波段的星上亮度温度，为i波段N个像元的平均星上亮度温度，为j波段N个像元的平均星上亮度温度；对于Landsat8数据，i，j分别为10，11；系数a和b可以通过式4‑5获得：w＝‑18.973(τ11/τ10)+19.13 R2＝0.9663，τ11/τ10＞0.9 (式4)w＝‑13.412(τ11/τ10)+14.158 R2＝0.9366，τ11/τ10＜0.9 (式5)Landsat 8第10波段和第11波段的星上亮度温度按下式计算：T＝K2/ln(1+K1/Lsat) (式6)其中，Lsat是星上辐射亮度，T是星上亮度温度，K1和K2为常数，从Landsat 8头文件获得；比辐射率ε利用NDVI(Normalized Difference Vegetation Index)阈值法来获取：NDVI=DNband5-DNband4DNband5+DNband4]]> (式7)其中DNband5和DNband4分别表示Landsat8第5波段和第4波段影像的DN值；当NDVI＜NDVIs时，ε＝εs，其中NDVIs是纯裸土区域的NDVI，εs是土壤的比辐射率；当NDVI＞NDVIv时，ε＝εv，其中NDVIv是纯植被区域的NDVI，εv是植被的比辐射率；当NDVIs≤NDVI≤NDVIv时，ε＝εs(1‑FVC)+εvFVCFVC是植被覆盖度：FVC=[NDVI-NDVIsNDVIv-NDVIs]2]]> (式8)NDVIs和NDVIv可以从图像上选取均质的裸土区域和植被区域来获取；εs和εv通过MODIS UCSB比辐射率库和Landsat 8TIRS波谱响应函数计算得到；(2‑2)水蒸汽光学厚度τw可表达为：τw=0.2385awλwM(1+20.07awλwM)0.45]]> (式9)其中w是可降水汽(cm)，awλ是水汽吸收系数，相对大气量M由下式获得：M＝[cosθz+0.15(93.885‑θz)‑1.253]‑1 (式10)θz是太阳天顶角；第三步、计算气溶胶光学厚度；(3‑1)计算太阳天顶角和太阳方位角，包括以下步骤：(3‑1‑1)计算太阳时：t=ts+0.170sin(4π(J-80)373)-0.129sin(2π(J-8)355)+12(SM-L)π]]> (式11)其中t是太阳时(小时数，带小数位)，ts是标准时(小时数，带小数位)，SM是该时区标准经线的经度(弧度)，L是该像元点的经度(弧度)，J是儒略日；(3‑1‑2)计算太阳赤纬：δ=0.4093sin(2π(J-81)368)]]> (式12)其中δ是太阳赤纬(弧度)，J是儒略日；(3‑1‑3)计算太阳天顶角：θz=π2-arcsin(sin l sinδ+cos l cosδcosπt12)]]> (式13)其中θz是太阳天顶角(弧度)，l是该像元点的纬度(弧度)，δ是太阳赤纬(弧度)，t是太阳时；(3‑1‑4)计算太阳方位角：(式14)其中是太阳方位角，θz是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬；(3‑2)计算Landsat8第一波段的星上反射率：ρsat＝(MρQcal+Aρ)/cos(θz) (式15)ρsat为星上反射率，Mρ为增益，Ap为偏置，这两个参数从Landsat8头文件获得；Qcal为影像DN值，θz是太阳天顶角；(3‑3)计算气溶胶光学厚度；(式16)其中，ρsat为星上反射率，ρ0为大气的路径辐射项等效反射率，μs和μv分别为太阳天顶角θz与观测天顶角θv的余弦，为相对方位角；r为地表反射率，S为大气整层向下的半球反射率，T(μs)和T(μv)分别为...
【技术特征摘要】
【专利技术属性】
技术研发人员：张兆明，何国金，王猛猛，龙腾飞，王桂周，张晓美，
申请(专利权)人：中国科学院遥感与数字地球研究所，
类型：发明
国别省市：北京;11
还没有人留言评论。发表了对其他浏览者有用的留言会获得科技券。