MODIS热红外温度数据反演的算法实现及其与地温数据的相关性分析

路 茜，张铁宝

（四川省地震局，成都610041）

0 引言

在地震前地热异常的早期研究中，可以发现不少关于地震、地质构造与地表及浅层地温关系研究的记录[1-3]。近年来，卫星遥感技术在现代光学、航天技术和计算机技术的基础上迅速地发展了起来。在地震预报和构造分析上的应用也取得了初步的进展，由于其具有覆盖范围广、实时监测及较高的空间分辨率等优点，成为了观测、分析地表热辐射的新技术手段。地表温度是地表与大气相互作用过程中的一个重要的物理参数。陆面温度遥感反演一直是许多学者力图解决的一个难题。其难点在于存在其他大气成分、地表层热状况以及热存储释放过程的复杂性和热传递方式等因素的影响。尽管如此，随着遥感技术的不断发展，应用卫星红外资料对陆面温度的研究已取得了很大的发展，陆面 “分裂窗”方法使得陆面温度反演由难变易。

1 MODIS数据前期处理

卫星遥感热红外亮温信息用于地震前兆异常分析研究的难点之一来自云层的干扰来自地表的热辐射信息难以穿透云层被卫星传感器接收到.若在温度反演中引入有云数据，则得到的反演结果存在干扰，其中必然存在云顶温度，即比实际的地表温度低很多。MODIS数据有36个光谱通道，多光谱联合云检测是MODIS数据的应用研究特色之一。 例如， 波段 0.65 μm、 0.85 μm、 1.38 μm、1.6 μm、 8.6 μm 和 11 μm 的组合应用就可以非常准确地得到不同云层的特性。多光谱的研究应用可以将一景数据的云和晴空分为多个不同的类别。云和雪在0.65 μm波段看起来极其相似，而在1.6 μm波段的图像中却大不相同；位于对流层之上的1.38 μm波段图象对于高空的冰晶云特别敏感；8.6 μm和11 μm的组合应用可以区分卷云和低空的可降水云[4]。根据有云数据特征，利用通道1及通道31阈值计算进行云检测，本文中的数据范围为中国大陆西部地区，采用了通道1反射率大于0.2与通道31亮温值小于273 K，该阈值会随着季节的变化作调整。

2 温度数据分裂窗法

为从6个热红外地表亮温波段中获取更为统一的数据源，并探寻多种MODIS数据应用方法，现对热红外亮温数据反演为温度数据的算法实现做初步研究。本文的方法是分裂窗地表温度普遍性算法。我们已知地表光谱反射率，根据多通道海面温度算法修正大气因素就能得到陆地表面温度[5]。这种方法跟其他方法的不同之处在于：①不要求精确的大气廓线；②不需要坐标系中模拟辐射转换；③其精确性取决于表面辐射率。分裂窗法地表温度普遍性算法采用一个简单模式，即使在处理过程中模拟辐射转换不能节省处理时间，处理时间也不长。该方法被作为MODIS地表温度反演算法的首要算法，是基于以下关于大气、地表温度和地表辐射的几点考虑的：① 由于大气（特别是水汽廓线）在垂直、水平方向上的变化幅度大，要计算其他相对湿气廓线，使之精确度达到误差10%以内是很困难的；②地表温度随着地区和时间变化，且与地面空气温度也联系较小，不同地表覆盖类型白天和夜晚温度可能相差10°C以上[6]；③ 多数地表覆盖类型的辐射率在MODIS 31、32通道中是相对稳定的，如图1所示。

图1 各种不同的地表覆盖类型在MODIS 31、32波段的平均辐射率（ε31、ε32分别为31、32波段的辐射率）Fig.1 Average radiation rate of MODIS band 31， 32 with various land cover types （ε31、 ε32are radiation rate of band31， 32）

3 NDVI植被归一化指数理论

论经验性植被指数是根据叶子的典型光谱反射率特征得到的。由于叶绿素吸收在蓝色（470 nm）和红色（670 nm）波段最敏感，可见光波段的反射能量很低。几乎所有的近红外（NIR）辐射都被散射，很少吸收，且散射程度因叶冠光学特性、结构特性而异。因此红色波段和近红外波段的反差是对植被量很敏感的度量波段。无植被或少植被区反差最小，中等植被区的反差是红色和近红外波段的变化结果，高植被区则只有近红外波段对反差有贡献，红色波段趋于饱和，不再变化。这种反差对比可以用比值（NIR/red）， 差分（NIR-red）， 线性组合（x1×red+x2×NIR）或上述三者的组合来增强。 一般来说，植被指数就是对这种对比的度量，也是叶冠结构 （LAI，LAD）和生物参数 （FAPAR）的函数。Deering[7]提出 “归一化差分植被指数”（NDVI），将比值限定在 [-1，1]范围内：

经过了辐射校正和几何初校正的MODIS数据源，再进一步对每日、每景图像进行几何精校正、除条纹、去云等处理，进而进行NDVI计算及合成。计算公式为：

其中b1、b2为MODIS的第1、2波段。

4 温度计算的实现

4.1 分裂窗法公式参数选取

经过对NDVI原理的初步了解以后，下面介绍运用NDVI极限方法来进行分裂窗法地表温度反演的方法。

分裂窗法的经典公式为[8]：

式（3）中， ε=0.5（ε31+ε32）， Δε=ε31-ε32。

利用NDVI极限方法来求解ε31、ε32的值，其关系为： 如果： NDVI＜NDVImin， 则 ε31=0.979 5-0.056 5 ρ1，ε32=0.981 5-0.027 5 ρ1， ρ1为第 1 波段反射率； 如果： NDVImin≤NDVI≤NDVImax， 则： ε31=0.968+0.021 fV，ε32=0.974-0.015 fV。fV=如果NDVI＞NDVImax,， 则 ε31=ε32=0.989 。

NDVImin、NDVImax是通过大气修正的裸露土壤和叶面积指数 LAI（Leaf Area Index）的 大小极NDVI值，取决于所研究的数据图像。

从前面谈到的植被指数理论可知，植被指数是对地表植被覆盖率的简单、有效和经验的衡量。它会随时间、空间变化。对于每个MODIS数据点来说，不同时间的植被指数就有所不同，因为要较准确的计算地表温度，需要根据NDVI的变化情况来正确地选择参数ε31、ε32。这里我们所采用的方法是建立NDVI最值数据库来存储每个数据点的NDVImin、NDVImax，在读取数据文件时获取每个MODIS数据点的NDVI值与当月NDVImin、NDVImax进行比较并替换，在此基础上查找上月该点所有 数据的最大值和最小值和当前值进行比较从而选择合适的参数ε31、ε32，这样就可以保证计算温度数据文件的实时准确性。

由于接收到的各数据点的观测角度都不尽相同，因此分裂窗法地表温度普遍性算法中公式的系数也不相同（如表1）。

MODIS扫描观测宽度达2 330 km，观测角度达到±58°，接收到的数据与经度线，卫星运行轨道方向与经度存在一定的夹角。接收到的MODIS数据经过定标定位及等经纬度投影以后，再按照所需要的区域截取包含所需热红外地表波段的数据文件，由此得到的亮温数据文件中的每个像元（Pixels，注下文用p表示）数据的观测角度均不同，多年数据的海量像元参与计算就需要一个庞大的观测角度数据库。为了便于计算，本文采用了各观测角度的系数取平均值的办法，由此得到的平均值系数如上表1中所示。

4.2 温度计算结果

综上所述步骤， 各参数值 （ε31、ε32、 A1、 A2、A3、B1、B2、B3）带入分裂窗法公式中进行计算，由此得到各点亮温数据的地表温度LST（Land Surface Temperature）。例如：下图2为由2004年3月29日北纬30°、东经100°为中心，长宽各为2 000 km的去云亮温文件（图 2a）与NDVI文件（图2b）而得到的LST数据文件（图 2c）。

4.3 与MODIS亮温产品比较

经过分裂窗算法计算出的温度文件，首先与MODIS地表亮温产品进行比较。比较时段为2004年3月至2006年11月。在MODIS地表亮温产品包含的6个热红外地表亮温波段，通道21、22、23的数据形态较一致，通道31、32的数据形态较一致，而通道20则与他们都不尽相同。在此给出几个像元点的反演温度与代表性的通道20、通道21、通道31的比较结果（图3）。

表1 不同观测角度与分裂窗算法系数的对应表Table 1 Corresponding relation between different observational angle and split window algorithm

图2 温度文件的反演过程示例图Fig.2 The example of assimilation procedure of the temperature files

图3 反演温度与MODIS地表温度产品的曲线图（4p）Fig.3 The sequence diagram of assimilation of LST and MODIS thermal infrared band（4p）

图3显示，较高纬度地区（a点）的反演温度及各数据产品之间相关性较中纬度地区（b、c、d点）高，可能因为a点位于新疆、青海交界地区，植被覆盖少，干扰较少，年变清晰。而b、c、d点分别位于缅甸热带雨林地区、四川甘孜州理塘地区、四川和贵州交界地区，植被覆盖较高，或云多，地表温度实测存在干扰因素较多。

经过多点的反复比较（表2）可以看出，由分裂窗算法计算出的反演温度与21和31通道的MODIS数据存在较好的相关性，其中与31通道相关系数接近1。而20通道的数据则存在较大的波动性，相关性较差。

表2 反演温度与MODIS亮温产品之间的相关系数Table 2 The correlation coefficient between assimilation of LST and MODIS thermal infrared band

5 MODIS反演温度数据与地温数据的初步相关性分析

5.1 气象地温数据和反演温度数据收集情况

本文所采用的气象数据地温资料为地面气象站观测所得，这些观测站是位于我省鲜水河、安宁河、龙门山等主要断裂带区域内的24个气象台站观测记录到的地表温度（0 cm）以及地下浅层（包括5～320 cm）温度，由于与反演所得到的地表层温度做比较，因此本文采用地表温度（0 cm）和地下浅层（40 cm）的日均值数据及其平均所得到的月均值数据。大部分台站的资料起止时间为2004年1月至2006年7月。

为了配合气象台站观测点地表温度数据的时间， MODIS反演温度数据也选定了上述相同时间数据作为分析数据。MODIS热红外波段数据的空间分辨率为1 km，而地面气象台观测地温是用以代表一个区域的地温。因此本文根据气象观测台的经纬度分别选取了MODIS反演温度数据的不同大小区域的均值进行相关性对比分析。反演温度数据经过多天数据拼接减少云数据覆盖和去云处理。

5.2 反演温度（20×20 p）与各深度层（0 cm、40 cm）地温日均值的相关性

利用前面分裂窗法计算所得的多年MODIS反演温度与地温的日均值数据进行对比分析。在此给出甘孜、宜宾、马尔康、康定四个地温观测点与反演温度的对比分析。本研究利用线性相关系数计算公式：

计算反演温度与地温的相关性。图4给出了反演温度曲线及四个地温观测点日均值曲线。

图4可以看出，地温日均值曲线较反演温度曲线稳定，波动小。这可能是由卫星传感器接收信号受到大气颗粒、太阳辐射等多方面因素的影响造成的。反演温度与地温都出现了在夏季时波动较大、冬季时较为稳定的现象。0 cm地温与反演温度的相似性比40 cm地温与反演温度的相似性好。表3给出了反演温度与4个地区地温的相关系数。

图4 反演温度与甘孜（a）、宜宾（b）、马尔康（c）、康定（d）0 m、 0.4 m 地温的日均值曲线图Fig.4 The daily mean profiles of assimilation of land surface temperature （LST） and geotemperature of Ganzi， Yibin， Maerkang，Kangding area

表3 反演温度（20×20p）与地温（日均值）之间的相关系数Table 3 The correlation coefficient between assimilation of LST （20×20p） and daily mean geotemperature

5.4 反演温度（20×20 p）与各深度层（0 cm、40 cm）地温月均值的相关性

用日均值数据进行月平均处理，可得到月均值时序图。由于篇幅的原因，在此仅给出月均值相关性分析的结果。从月均值的相关系数来看，比日均值有了明显的提高（图4）。

5.3 不同大小区域范围反演温度与地温日均值相关性的比较分析

地温数据来自与气象台站的测定，而反演温度则是依据气象台站的准确经纬度来确定的数据点位的单个像元或多个像元的均值。卫星数据受到大气、辐射、地表形态等多种因素的影响，单个像元的数据可用率较低，可靠性也较低，通过多像元的数据求均值可以弱化干扰，更能反映区域的实际情况。本文就这个问题，采用不同大小区域范围反演温度与地温进行对比、相关性分析。表5给出了结果，由此可以看出除了马尔康地区以外，其他三个地区经过多像元数据平均后，一定程度地弱化了干扰因素，相关系数有所提高。甘孜、马尔康、康定属于高原地区，植被覆盖率较低、地表覆盖类型较单一，这可能是其实测地温与各大小区域范围反演温度相关系数相差不大的原因，而宜宾位于四川盆地，长江流域，地貌以中低山地和丘陵为主，从表5中可以看出最大区域范围比单个像元与反演温度的相关系数相差0.2。

表4 反演温度（20×20 p）与地温（月均值）之间的相关系数Table 4 The correlation coefficient between assimilation of LST （20×20 p） and monthly mean

表5 不同大小区域范围反演温度与地温（日均值）之间的相关系数Table 5 The correlation coefficient between assimilation of LST of different regional ranges and monthly mean geotemperature

6 讨论

（1）分裂窗法地表温度算法实现中存在的误差因素：公式中的7个系数的确定存在一定的误差，本文在确定的时候是采取的取均值的方法，但实际观测角度不同时，该计算系数是有一定误差的；还有的选取也不是非常准确，由于的季节性变化，其值也是不稳定的，值域范围也不容易确定，该如何来确定其一定时间范围的最大、最小值，使得参数更加准确适合，还需要反复的试验。

（2）据上文分析得出，地表温度与遥感数据反演温度的日值相关性系数仅在0.6左右，而月值相关性却在0.8左右。这可能受到两方面原因的影响：①地温数据是由气象台站在每天多个时段测量所得温度值的平均值，而遥感卫星数据则是每天某个特定的时刻接收到的，本文所用的TERRA卫星数据是在上午10点至12点的实测数据。这两者之间必定存在差异。② 卫星遥感属于 “面”观测，而地面温度测点属于 “点”观测，通过卫星遥感获得的地表温度代表一定空间区域内的平均温度，通过地表实测的方式获得的地表温度代表测点的温度。遥感与实测之间，空间尺度存在不确定性[9]，这可能为二者带来差异。③ 遥感热红外数据的每日数值随机波动远大于地表测定温度，可能是因为遥感数据受到大气吸收、植被覆盖、人类活动等干扰引起的，而进行月均值处理以后这些干扰明显被弱化了。将遥感卫星热红外数据应用于地震监测预报工作时，应尽量采取一些去噪声的处理，如小波变换、滑动滤波等，以便于在应用中提取地壳内部变化引起的热红外信息。

致谢：本文所采用的地表及浅层地表温度数据来自四川省气象局，在此对他们的大力支持表示感谢。

参考文献：

[1]尤传侠.浅层地温在滇西南地震预报中的应用 [J].地震研究，1989，12（3）：254-259.

[2]张治洮.浅层地温变化与强震、旱涝的关系[J].中国减灾,1992，2（3）：50-50.

[3]刘小凤.青藏高原北部浅层地温异常特征及中短期地震预测[J].高原地震，2002，14（2）：1-7.

[4]刘玉洁，杨忠东.MODIS遥感信息处理原理与算法[M].北京：科学出版社，2001.

[5]Wan Z.J.Dozier.A.Generalized split-window algorithm for retrieving land surface temperature measurement from space[J].IEEE Trans Geosci Remote sensing，1996（34）：892-905.

[6] Betts A K，Ball J H，Beljaars A C M，et al.The land surface-atmosphere interaction：A review based on observation and global modeling perspectives[J].J.Geophys.Res，1996，101（D3）：7209-7225.

[7] Deering D W.Rangeland Reflectance Characteristics measured by aircraft and Spacecraft Sensor[D].Ph.D.Diss,Texas A&M University，1978.

[8]Zheng M W.MODIS Land-Surface Temperature Algorithm Theoretical Basis Document （LST ATBD）.Institute for Computational Earth System Science [D].University of California，Santa Barbara，1999.

[9] 陈顺云，刘培洵，刘力强，等.遥感与实测地表温度的对比分析及在地震研究中的意义[J].地球物理学报，2011，54（3）：747-755.