收透过率和气溶胶散射透过率，两者之积为气溶胶透过率，Fg是气溶 胶散射向下分量的比例，各参数的具体计算参见文献怔irdRE，RiordanC.Simplesolar spectralmodelfordirect曰nddiffuseirr曰di曰neeonhorizont曰1 曰ndtilted planesattheearth'ssurfaceforcloudlessatmospheres.Journalofclimateand appliedmeteorology, 1986, 25 :87-97.]。
[0074] 气溶胶光学厚度可W利用浓密植被法来反演，然而运种方法只适用于影像中存在 浓密植被的情况，对于影像中不存在浓密植被的情况（比如沙漠地区），该方法不适用。本 发明基于Landsa巧第1波段来反演气溶胶光学厚度。
[00巧]气溶胶光学厚度反演基于W下公式：
[0076]
(式 23)
[0077]其中，L,。,为星上福射亮度，y郝yy分别为太阳天顶角0,与观测天顶角0V的 余弦，，为相对方位角；L。为观测方向的路径福射项，r为地表反射率，S为大气整层向下的 半球反射率，T为大气透过率，T(yJ和T(iiy)分别为太阳照射方向和传感器观测方向的 大气透过率，F。为大气层顶太阳光的福射通量密度。
[0078] 利用垂直入射的太阳福射量iigF。对上式进行归一化可得星上反射率Pg。，
[007引
（式 24)
[0080] 其中，P。为大气的路径福射项等效反射率。
[008。 上式中，S、P。和T(yjT(i0运3个参数代表大气的状况，因此可W从中获取大 气气溶胶光学厚度。在LandsatS第一波段，卫星传感器接收到的福射中大气反射的贡献较 大，地表反射的贡献微弱，假定同期的地表反射率不变，将晴好天的地表反射率（本发明利 用Landsa巧卫星过境前的版)DIS8天合成地表反射率产品，M孤09产品（第=波段））代 入上式反演气溶胶光学厚度。在具体反演中是用福射传输模型构建气溶胶光学厚度和太阳 天顶角、观测天顶角和相对方位角之间的查找表，然后利用MODIS的8天合成地表反射率产 品去除地表反射贡献，得到气溶胶光学厚度。
[0082] 由于Landsa巧影像可W看做是星下点成像，整景影像范围内观测天顶角可W近 似为0,运样就无需考虑观测天顶角和观测方位角。下面介绍如何逐像元地计算太阳天顶角 和太阳方位角，包括W下步骤：
[008引 （1)计算太阳时：
[0084]
(式25)
[0085]其中t是太阳时（小时数，带小数位），ts是标准时（小时数，带小数位），SM是该 时区标准经线的经度（弧度），L是该像元点的经度（弧度），J是儒略日。
[0086] 似计算太阳赤缔：
[0087]
(式 26)
[008引其中5是太阳赤缔（弧度），J是儒略日。
[0089] (3)计算太阳天顶角
[0090]
(式 27)
[0091] 其中0Z是太阳天顶角（弧度），1是该像元点的缔度（弧度），5是太阳赤缔（弧 度），t是太阳时。
[0092] (4)计算太阳方位角
[0093]
(式28)
[0094] 其中批是太阳方位角，0虎太阳天顶角，1是该像元点的缔度，5是太阳赤缔。
[0095] 按式（29)计算Landsa巧第一波段的星上反射率：
[009引 Psat=(MpQcal+Ap)/c0S(目Z)(式 29)
[0097] P为星上反射率，MP为增益，AP为偏置，运两个参数从Landsa巧头文件获得； Qca巧影像DN值，目Z是太阳天顶角。
[0098] 查找表利用6S模型设定不同的条件进行福射传输计算得到。设定参数包括：12个 太阳天顶角（0-66度，间隔6度），16个太阳与卫星之间的相对方位角（0-180度，间隔12 度），大气气溶胶模式参数假设为大陆性气溶胶，并设立20个大气气溶胶光学厚度值（0-2， 间隔0. 1)，波段自定义为Landsat8 0LI传感器第一波段的响应函数。
[0099] 根据计算得到的观测几何（太阳天顶角和相对方位角），对查找表进行线性插值， 得到不同光学厚度下的大气参数S、P。和T(yjT(iiy)，代入式（24)，同时将Landsa巧卫 星过境前的M0DIS8天合成地表反射率产品（M0D09产品，第=波段）代入式（24)获得不 同气溶胶光学厚度下的星上反射率；然后利用式29得到的Landsa巧第一波段的星上反射 率进行线性插值，得到气溶胶光学厚度。
【主权项】
1. 一种Landsat 8卫星数据地表反射率反演方法，其步骤为： 第一步、计算星上辐射亮度； Lsat= M LQcal+AL (式 1) 其中，Lsat是星上辐射亮度，M A波段的增益，A A波段的偏置，Q 为影像DN值，M # A1^从Landsat 8头文件获得； 第二步、计算水蒸汽光学厚度； (2-1)获取可降水汽w ; w = a ( τ / τ ;) +b (式 2)其中，TiSi波段的大气透过率，τ ,为」波段的大气透过率，ε 波段的比辐射 率，^为j波段的比辐射率，k表示第k个像元，T 1>k为第k个像元i波段的星上亮度温度， Thk为第k个像元j波段的星上亮度温度，f为i波段N个像元的平均星上亮度温度，f为 j波段N个像元的平均星上亮度温度；对于LandsatS数据，i，j分别为10,11 ; 系数a和b可以通过式4-5获得： w = -18. 973 ( τ J τ 10)+19. 13 R2= 0. 9663, τ J τ 10> 〇. 9 (式 4) w = -13. 412( τ J τ 10)+14. 158 R2= 0. 9366, τ J τ 10< 〇. 9 (式 5) Landsat 8第10波段和第11波段的星上亮度温度按下式计算： T = K2An(HK1Asat)(式 6) 其中，Lsat是星上辐射亮度，T是星上亮度温度，Kl和Κ2为常数，从Landsat 8头文件 获得； 比福射率 ε 利用 NDVI (Normalized Difference Vegetation Index)阈值法来获取：其中DNbandJP DN band4*别表示Landsat8第5波段和第4波段影像的DN值； 当NDVI<NDVIS时，ε = ε s，其中勵￥1；3是纯裸土区域的NDVI，ε 5是土壤的比辐射 率； 当NDVI>NDVIv时，ε = ε ν，其中勵￥1,是纯植被区域的NDVI，ε ￥是植被的比辐射 率； 当 NDVIsS NDVI 彡 NDVIv时，ε = ε s(l-FVC)+evFVC FVC是植被覆盖度：NDVI，NDVIv可以从图像上选取均质的裸土区域和植被区域来获取；ε，ε ν通过 MODIS UCSB比辐射率库和Landsat 8TIRS波谱响应函数计算得到； (2-2)水蒸汽光学厚度Tw可表达为：其中W是可降水汽（cm)，awA是水汽吸收系数，相对大气量M由下式获得： M = [cos Θ z+〇. 15(93. 885- Θ z) L 253] 1 (式 10) Θ 2是太阳天顶角； 第三步、计算气溶胶光学厚度； (3-1)计算太阳天顶角和太阳方位角，包括以下步骤： (3-1-1)计算太阳时：其中t是太阳时（小时数，带小数位），、是标准时（小时数，带小数位），SM是该时区 标准经线的经度（弧度），L是该像元点的经度（弧度），J是儒略日； (3-卜2)计算太阳赤炜：其中δ是太阳赤炜（弧度），J是儒略日； (3-1-3)计算太阳天顶角：其中θζ是太阳天顶角（弧度），1是该像元点的炜度（弧度），δ是太阳赤炜（弧度）， t是太阳时； (3-1-4)计算太阳方位角：其中％是太阳方位角，92是太阳天顶角，1是该像元点的炜度，δ是太阳赤炜； (3-2)计算LandsatS第一波段的星上反射率： P sat= (M PQcai+Ap)/cos( θ ζ)(式 15) P sat为星上反射率，M ρ为增益，A ρ为偏置，这两个参数从Landsat8头文件获得；Q ^为 影像DN值，02是太阳天顶角； (3-3)计算气溶胶光学厚度；I 其中，Psat为星上反射率，P。为大气的路径辐射项等效反射率，μ 3和μ v分别为太阳 天顶角θζ与观测天顶角Θ 余弦，A为相对方位角^为地表反射率，S为大气整层向下 的半球反射率，T(ys)和τ(μν)分别为太阳照射方向和传感器观测方向的大气透过率； 利用6S模型设定不同的条件构建查找表；设定参数包括：12个太阳天顶角（0-66度， 间隔6度），16个太阳与卫星之间的相对方位角（0-180度，间隔12度），大气气溶胶模式参 数假设为大陆性气溶胶，并设立20个大气气溶胶光学厚度值（0-2,间隔0. 1)，波段自定义 为Landsat 8 OLI传感器第一波段的响应函数； 根据计算得到的观测几何（太阳天项角和相对方位角），对查找表进行线性插值，得到 不同光学厚度下的大气参数s、P。和τ(μ 3)τ(μν)，代入式（16)，同时将LandsatS卫星过 境前的M0DIS8天合成地表反射率产品（M0D09产品，第三波段）代入式（16)获得不同气溶 胶光学厚度下的星上反射率；然后利用式15得到的LandsatS第一波段的星上反射率进行 线性插值，得到气溶胶光学厚度； 第四步、大气校正和地表反射率反演； (4-1)计算大气程辐射Lp; Lp=MLQCALdark+AL (式 17) QCALdarii是影像中暗目标的壳度值； (4-2)计算日地距离d; 1/d2= 1. 000110+0. 034221cos Γ +0. 001280sin Γ +0. 000719cos2 Γ +0. 000077sin2 Γ (式 18) Γ = 231 (dn-l)/365 (式 19) dn为儒略曰； (4-3)计算天空光漫射到地表面的光谱辐照度Ed_; Edown= Er+Ea (式 20) 其中，民为瑞利散射分量，E a为气溶胶散射分量；其中，E。是大气层外相应波长的太阳光谱辐照度，d是日地距离，Θ 2是太阳天顶角，Tc^ Tw和1\分别是臭氧吸收透过率、水汽吸收透过率和瑞利散射透过率，T aa和T as分别是气溶 胶吸收透过率和气溶胶散射透过率，Fs是气溶胶散射向下分量的比例； (4-4)计算瑞利散射光学厚度τ ^和臭氧吸收光学厚度τ ^其中，λ是Landsat8影像各波段的中心波长，ym;h是海拔高度，km; (4-5)计算大气透过率； Tz = exp (- τ /cos τ ζ) = exp {(- τ r- τ a- τ 〇- τ w) /cos θ J (式 25) Tv= exp (- τ /cos θ ν) = exp {(- τ r- τ a- τ。- τ w) /cos θ y}(式 25) 其中，TJPTv分别是太阳照射方向和传感器观测方向的大气透过率，τ ρ τ3、^和τ w 分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度， 02是太阳天顶角，Θ v是传感器观测天顶角； (4-6)计算地表反射率；其中，P是地表反射率，Lsat是星上辐射亮度，Lp是程辐射，d是日地距离，1\是传感器 观测方向的大气透过率，E。是大气层外相应波长的太阳光谱辐照度，θ z是太阳天顶角，Tz 是太阳照射方向上的大气透过率，Ed_是天空光漫射到地表面的光谱辐照度。
【专利摘要】一种Landsat8卫星数据地表反射率反演方法，该方法主要基于Landsat8影像本身来反演地表反射率，对外部数据源的需求非常低，容易获取，克服了传统Landsat数据地表反射率反演必须依赖较多外部数据源造成的局限，因此该方法具有较强的实用性，对于实现利用Landsat8数据业务化地生产地表反射率产品具有重要的现实意义。
【IPC分类】G06F17/50
【公开号】CN105183989
【申请号】CN201510563531
【发明人】张兆明, 何国金, 王猛猛, 龙腾飞, 王桂周, 张晓美
【申请人】中国科学院遥感与数字地球研究所
【公开日】2015年12月23日
【申请日】2015年9月8日