ENVI下基于劈窗算法从MODIS数据中反演海表温度

ENVI下基于劈窗算法从MODIS数据中反演海表温度(2013-12-03 19:42:18)转载▼分类：遥感解决方案标签：劈窗算法modis反演温度教育劈窗算法最初是为反演海面温度开发的，具体地说是针对NOAA/AVHRR的4和5通道设计的，后来也被用来反演地表温度，这种算法较成熟，精度也高。

劈窗算法以地表热辐射传导方程为基础，利用10~13μm 大气窗口内，两个相邻热红外通道(一般为10.5～11.5μm、11.5～12.5μm)对大气吸收作用的不同，通过两个通道测量值的各种组合来剔除大气的影响，进行大气和地表比辐射率的修正。

表达式为：T S= T 4+ A (T 4- T 5) + B其中：T S为地表真实温度，T 4和 T 5分别为AVHRR的4和5通道，A和B为常量。

AVHRR 的通道4(10.15~11.13μm) 和通道5 (11.15~12.15μm) 恰与MODIS 的第31 波段 (10.178~ 11.128μm) 和32 波段 (11.177~ 12.127μm ) 的中心波长相对应, 可将MODIS 的31、32 波段资料, 用于劈窗算法进行地表温度计算。

很多学者对这个算法进行了推演，得到很多新的算法，如覃志豪、毛克彪等人。

本文就是使用其他学者推演的算法。

利用MODIS数据劈窗算法反演海表温度技术流程如下图：图1技术流程图注：（1）按照本流程反演出来的结果是SST。

陆地上的值可以视为无效值，若要得到正确的陆表温度，需要加入海陆分离的步骤，以及城镇和自然表面的比辐射率计算。

（2）MODIS数据下载：nasa官网：/.下面详细介绍处理流程。

第一步：数据打开Open->Open as->EOS->MODIS，选择MOD021KM.A2013185.0245.005.2013185094144.hdf文件，点击OK打开。

分为三个数据集：热红外数据（20-36波段），可见光到短波红外的辐射率数据（1-19、26波段），可见光到短波红外的反射率数据（1-19、26波段）。

第二步：辐射定标由于ENVI默认的设置是对MODIS数据进行自动的辐射定标，所以打开的数据即是经过了辐射定标的数据。

图2打开modis 02级数据第三步：几何校正MODIS1b数据是hdf格式，自带有经纬度坐标信息，可自动进行几何校正，分别对热红外数据集和反射率数据集进行自动几何校正。

反射率数据集合校正：（1）打开工具/Geometric Correction/Georeference by Sensor/Georeference MODIS，选择反射率数据集，点击Spatial Subset，选择大亚湾的大概位置，点击Spectral Subet，选择2和19波段。

点击OK。

图3选择数据同时选择空间子集和光谱子集（2）在Georeference MODIS Parameter面板设置为UTM， WGS84，49带。

图4 Georeference MODIS Parameter面板（3）点击OK，在结果输出面板，设置路径和文件名输出。

热红外数据集几何校正：方法与反射率数据集校正方法一样，选择31和32两个波段校正。

第四步：海表温度反演本文使用的是《高懋芳，覃志豪等，MODIS数据反演地表温度的基本参数估计方法》中分裂窗算法模型进行海表温度反演，旨在学习ENVI中的操作流程。

算法为：T s = A 0 + A 1*T 31 - A 2* T 32 ( 公式1)其中：T s 是地表温度 ( K ) , T 31 和 T 32 分别是M OD I S 第 31 和 32 波段的亮度温度； A 0, A 1 和 A 2 是分裂窗算法的参数，分别定义如下:A 0 = [ D 32( 1 - C31 - D 31) / ( D 32 C31 -D 31 C32) ] a31 - [ D 31( 1 - C32 - D 32) /( D 32 C31 - D31 C32) ] a 32 (公式2)A 1 = 1 + D 31/ ( D32C31 - D 31 C32) + [ D 32( 1 -C31 - D 31) / ( D 32C31 - D31 C32) ]b31 ( 公式3)A 2 = D 31/ ( D 32 C31 - D 31C32) + [ D 31( 1 - C32 - D 32) / ( D 32 C31 - D 31C32) ]b32 ( 公式4)式中， a 31，b31，a 32和 b 32是常量, 根据 MODIS的波段特征确定, 在地表温度 0 ~ 5 e 范围内, 这些常量分别可取 a 31 = -64.60363 , b31 = 0.440817, a 32 = -68.72575, b32 = 0.473453。

上述公式的中间参数分别计算如下:C i =Ɛiτi (ɵ) ( 公式6)D i = [ 1 -τi (ɵ)] [ 1 + ( 1 -Ɛi) τi (ɵ)] ( 公式7)式中: i 是指 MODIS的第31和32波段, 分别为 i= 31 或 32; τi (ɵ)是视角为ɵ的大气透过率; Ɛi是波段 i的地表比辐射率。

由以上公式可以看出, 该算法要求卫星遥感器的31 和 32 波段数据来计算星上亮度温度, 同时还要求已知大气透过率和地表比辐射率, 才能进行地表温度的反演。

下面是详细操作步骤：（1）大气透过率计算大气透过率τi (ɵ)是计算地表温度的基本参数, 通常是通过大气水汽含量来估计。

经过前人研究，可以用MODIS第 2 和 19 波段来反演大气水分含量，然后再根据大气水分含量与大气透过率之间的关系来估计大气透过率。

对于MODIS图像中的任何一个像元，其可能的大气水分含量用下式估计（公式8）式中: ω是大气水分含量( g *cm-2) , α和β常量，分别α= 0.02 和β=0.6321；ρ19和ρ2分别是MODIS第19和2波段的地面反射率。

使用bandmath工具计算大气水分含量：表达式：((0.02-alog(b19/b2))/0.6321)^2B19：第19波段反射率B2：第2波段反射率大气透过率的计算中，水汽是最主要的考虑因素，毛克彪等将MODTRAN等大气模型模拟出来的两者的关系，应用到MODIS数据中，提高了地表温度反演的精度和实时性，本文采用模拟效果较好的指数关系模拟方程，拟合度达到了0.99以上，公式为：τ31= 2.89798-1.88366*exp{-[ω/（-21.22704）]} （公式9）τ32= -3.59289+4.60414*exp{-[ω/（-32.70639）]} （公式10）式中，ω是水汽含量。

使用bandmath工具计算大气透过率：31波段大气透过率=2.89798-1.88366*exp(b1/21.22704)32波段大气透过率=-3.59289+4.60414*exp(b1/32.70639)B1：大气水分含量。

（2）地表比辐射率的估算地表比辐射率主要取决于地表的物质结构，对MODIS来说，大致分为水面、城镇和自然表面。

对于反演来说，利用混合像元分解的方法，根据植被覆盖率来计算自然表面和城镇的比辐射率，水体的可以用常量：Ɛ=0.996，Ɛ32水体=0.992。

31水体有了这些参数，我们就可以计算C31、C32、D31、D32中间参数，BandMath表示式分别为：C31=0.996*b31C32=0.992*b32D31=(1-b31)*(1+(1-0.996)*b31)D32=(1-b32)*(1+(1-0.992)*b32)其中 B31：31波段大气透过率B32：32波段大气透过率（3）亮度温度的计算将图像DN值定标维热辐射强度之后, 可用Planck函数求解出星上亮度温度, 计算公式如下：T i = K i 2 / l n ( 1 + K i 1 /I i )式中, K i 1和 K i 2是常量，对于第 i = 31 波段, 分别为 K 31 , 1 = 729 .541636 W•m-2•sr-1•um-1，K 31 , 2 = 1304.413871K ; 对于第 i = 32 波段, 为 K 32 , 1 = 474 . 6 84780 W•m-2•sr-1•um-1，, K 31 , 2= 1196 . 978785 K。

使用bandmath工具计算31和32的亮温。

31波段亮温=1304.413871/alog(1+729.541636/b31)B31: 31波段辐射亮度值32波段亮温=1196.978785 /alog(1+474.684780/b32)B32：32波段辐射亮度值（4） A 0, A 1和A 2参数计算接下来我们计算A 0, A 1 和 A 2参数，Bandmath表达式分别为：A 0=b4*(1-b1-b3)/(b4*b1-b3*b2)*(-64.60363)-b3*(1-b2-b4)/(b4*b1-b3*b2)*(-68.72575)A 1=1+b3/( b4*b1-b3*b2)+b4*(1-b1-b3)/( b4*b1-b3*b2)* 0.440817A 2=b3/(b4*b1-b3*b2)+b3*(1-b2-b4)/(b4*b1-b3*b2)* 0.473453其中，b1：C31B2：C32B3：D31B4：D32（5）温度计算把这些参数带入公式1中计算温度值，Bandmath表达式为：T s=b0+b1*b31-b2*b32-273其中：b0：A0参数B1：A1参数B2：A2参数B31：B31亮温值B32：B32亮温值如下为得到的最终地表温度反演结果，单位为摄氏度。

通过/Statistics/Compute Statistics 工具统计看到最初的结果会有一些数量很少的异常值，几十个像素，这些异常值大多是处于影像的边缘。

如下图为统计的结果，小于-40的像元有57个，大于32度的34个像元。

可以将这些像元处理，如用以下bandmath处理：-40>b1<32。

图5初始结果统计图6最终反演结果。