锐客网

  1. 首页 > it技术 > >

【RS学习笔记】基于MODIS劈窗算法地表温度反演

目的：
1）掌握基于遥感的陆表温度反演方面的理论知识。
2）能够根据劈窗算法反演地表温度。
步骤：
1 理论知识
陆地表面温度(Land Surface Temperature ,LST) 是一个重要的地球物理参数，用卫星遥感数据提取海洋温度(Sea Surface Temperature ,SST)已较为成熟，可以在全球范围内达到1K精度。由于陆地表面比海洋表面复杂得多,导致陆地表面温度反演的精度较低，陆地表面温度反演成了一个亟待解决的科学难题。
现有的地表温度反演算法大致有以下四种：大气校正法、单通道算法、分裂窗算法和多波段算法。大气校正法和单通道算法需要大气实时剖面数据，单通道算法适用于只有一个热红外波段的数据，如 Landsat TM /ETM数据；劈窗算法适合于两个热红外波段的数据，如 MOAAAVHRR 和MODIS；多波段算法适合于多个热红外波段的数据，所需参数多，运算复杂且需要白天晚上两景数据，反演难度较大，就成熟程度而言，多波段算法还在发展之中。到目前为止，劈窗算法是目前发展比较成熟的地表温度遥感反演方法。这一算法需要两个彼此相邻的热红外波段遥感数据来进行地表温度的反演。劈窗算法主要是针对 NOAA-AVHRRR 的热红外通道4和5数据来推导。在 MODIS 的8个热红外波段中，第31和第32波段最接近于 AVHRR 通道4和5的波段范围，因而最适用于分裂窗算法。
本文主要学习基于劈窗算法，利用MODIS 数据反演海水温度，涉及MODIS数据的几何校正、裁剪、和 Band math 工具使用。基于 ENVI5.3软件下操作完成。实验流程如图所示
【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片

2 预处理
2.1 打开数据
打开ENVI 5.3软件，点击【Open】/【Open as】/【EOS】/【MODIS】，选择MOD021KM.A2013185.0245.005.2013185094144.hdf
文件，点击OK打开。
可以看到数据分为三个数据集：热红外数据（20-36波段），可见光到短波红外的辐射率数据（1-19、26波段），可见光到短波红外的反射率数据（1-19、26 波段）。
【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片

2 辐射定标
由于ENVI默认设置是对MODIS数据进行自动辐射定标，所以打开的数据即是已经过辐射定标的数据。
2.3 几何校正
MODIS1b数据为hdf格式，自带有经纬度坐标信息，可自动进行几何校正，分别对热红外 数据集和反射率数据集进行自动几何校正。
【【RS学习笔记】基于MODIS劈窗算法地表温度反演】（1）反射率数据几何校正。打开工具/【Geometric Correction】/【Georeference by Sensor】/【Georeference
MODIS】，选择反射率数据集，点击 【Spatial Subset】，选择大亚湾的大概位置，点击【Spectral Subet】，选择2和19波段
【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片

（2）热红外数据集几何校正。打开工具/【Geometric Correction】/【Georeference by Sensor】/【Georeference
MODIS】，选择反射率数据集，点击 【Spatial Subset】，选择大亚湾的大概位置，点击【SpectralSubet】，选择31和32波段。
【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片

3 海表面温度反演
本文使用《MODIS数据反演地表温度的基本参数估计方法》中分裂窗算法模型进行海表温度反演，旨在学习
ENVI 中的操作流程。算法为：
Ts= A 0 + A 1T31 - A2 T32 （1）
其中：Ts是地表温度（K）,T31和T32分别是MODIS第31和32波段的亮度温度；A0, A 1和A2是分裂窗算法的参数，分别定义如下:
A0 = [D 32(1- C31 - D31) / (D32 C31 -D31C32)] a31 - [D31(1 - C32 - D32)/(D32 C31 - D31C32)] a32 （2）
A1 = 1+D31/(D32C31-D31C32)+[D32(1-C31- D31) / ( D32C31 - D31C32)] b31 （3）
A2 = D31/(D32C31 - D31C32) + [D31(1- C32 - D32) / (D32 C31 -D31C32)]b32 （4）
式中，a31，b31，a32和b32是常量,根据MODIS的波段特征确定,
在地表温度0 ~ 5 e范围内,这些常量分别可取a31 = -64.60363 , b31= 0.440817, a32 = -68.72575, b32 = 0.473453。
上述公式的中间参数分别计算如下:
Ci= ?iτi(?) （5）
Di = [1 -τi (?)] [1 + (1- ?i) τi (?)] （6）
式中: i是指MODIS 的第31和32波段,分别为i= 31或32; τi (?)是视角为 ?的大气透过率; ?i 是波段i的地表比辐射率。
由以上公式可以看出, 该算法要求卫星遥感器31和32 波段数据来计算星上亮度温度, 同时还要求已知大气透过率和地表比辐射率,才能进行地表温度的反演。
3.1 大气透过率计算
大气透过率τi (?)是计算地表温度的基本参数, 通常是通过大气水汽含量来估计。经前人研究，可以用MODIS 第2和19波段来反演大气水分含量，然后再根据大气水分含量与大气透过率之间的关系来估计大气透过率。对于MODIS 图像中的任何一个像元，其可能的大气水分含量用下式估计：
【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片
（7）
式中: ω是大气水分含量(g*cm-2) ，α和β常量，分别α= 0.02 和 β=0.6321；ρ19 和ρ2分别是MODIS 第19和第2波段的地面反射率。
在Envi中，点击工具箱/【Band Algebra】/【Band Math】工具计算大气水分含量，输入：((0.02-alog(b19/b2))/0.6321)^2
B19：第19波段反射率，B2：第2波段反射率。
大气透过率的计算中，水汽是最主要的考虑因素，毛克彪等将MODTRAN 等大气模型模拟出来的两者的关系，应用到MODIS数据中，提高了地表温度反演的精度和实时性，本文采用模拟效果较好的指数关系模拟方程，拟合度达到了0.99 以上，公式为：
τ31= 2.89798-1.88366*exp{-[ω/（-21.22704）]} （8）
τ32= -3.59289+4.60414*exp{-[ω/（-32.70639）]} （9）
式中，ω是水汽含量。点击工具箱/【Band Algebra】/【Band Math】工具计算大气透过率：
31 波段大气透过率表达式：2.89798-1.88366*exp(b1/21.22704)
32 波段大气透过率表达式：-3.59289+4.60414*exp(b1/32.70639)
B1：大气水分含量
3.2 地表比辐射率计算
地表比辐射率主要取决于地表的物质结构，对MODIS来说，大致分为水面、城镇和自然表面。对于反演来说，利用混合像元分解的方法，根据植被覆盖率来计算自然表面和城镇的比辐射率，水体的可以用常量：?31水体=0.996，?32水体=0.992。
有了这些参数，我们就可以计算 C31、C32、D31、D32中间参数，在工具箱/【Band Algebra】/【Band Math】工具中分别输入：
C31=0.996*b31
C32=0.992*b32
D31=(1-b31)*(1+(1-0.996)*b31)
D32=(1-b32)*(1+(1-0.996)*b32)
其中B31：31波段大气透过率，B32：32波段大气透过率温计算
将图像DN值定标为热辐射强度之后, 可用Planck函数求解出星上亮度温度, 计算公式如下：
Ti= Ki2 / ln (1+Ki1 /Ii )
式中, Ki1和Ki2是常量，对于第i = 31 波段, 分别为K31,1 =729 .541636 W?m-2?sr-1?um-1，K31,2 = 1304.413871K ; 对于第i = 32 波段, 为 K32,1= 474 . 6 84780 W?m-2?sr-1?um-1，, K 31,2 = 1196 . 978785 K。点击工具箱/【Band Algebra】/【Band Math】分别输计算31 和 32 的亮温。
31 波段亮温：1304.413871/alog(1+729.541636/b31)
32 波段亮温：1196.978785
/alog(1+474.684780/b32)
B31: 31波段辐射亮度值；B32：32波段辐射亮度值
3.4 A0、A1、A2参数计算
在【Band Math】中依次输入：
A0=b4*(1-b1-b3)/(b4b1-b3b2)(-64.60363)-b3(1-b2-b4)/(b4b1-b3b2)*(-68.72575)
A1=1+b3/(b4b1-b3b2)+b4*(1-b1-b3)/( b4b1-b3b2)* 0.440817
A2=b3/(b4b1-b3b2)+b3*(1-b2-b4)/(b4b1-b3b2)*0.473453
其中，b1：C31，B2：C32 ，B3：D31，B4：D32
3.5 温度计算
把这些参数带入公式1中计算温度值，在【Band math】中输入：
Ts=b0+b1b31-b2b32-273
其中，b0：A0参数，B1：A1参数，B2：A2参数，B31：B31 亮温值，B32：B32亮温值。
点击【Statistics】/【Compute Statistics】工具统计看到该结果会有一些数量很少的异常值（图3.14），几十个像素，这些异常值大多是处于影像的边缘：小于-40的像元有57个，大于32度的34个像元。可以将这些像元处理，在【Band math】中输入：-40>b1<32，得到最终反演结果。
【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片

【RS学习笔记】基于MODIS劈窗算法地表温度反演
文章图片

    推荐阅读


    上一篇：canvas|canvas制作一个简易时钟

    下一篇：js scroll案例,当到一定程度时,元素由绝对定位改为固定定位window.pageXOffset