基于遥感影像的漠河市地表温度反演

Land Surface Temperature Retrieval of Mohe City Based on Remote Sensing Image

摘要: 地表温度不仅对全球性的气候变化有重要意义，而且在各种灾害的监测预测方面，以及现代城市建设中也起到了重要作用。随着遥感卫星技术的发展，卫星遥感影像数据被广泛应用于众多领域之中，Landsat系列卫星影像的广泛应用对推动遥感技术的发展起到了重要作用。本文以黑龙江省漠河市为研究区域，利用2018年8月23日和2018年8月30日的Landsat8卫星影像，对地表温度进行反演，对反演结果进行分析。研究结果表明，漠河市最高温度为36.5909℃，最低温度为14.2906℃，平均温度为23.2028℃，标准差为1.5777，全区气温偏低，低温区位于水体及植被覆盖度高的地区，高温区位于城镇人口密度较大的地区，对研究本地城市热岛效应和森林火灾预测监控具有重要意义。

Abstract: Land surface temperature is not only of great significance to global climate change, but also plays an important role in the monitoring and prediction of various disasters, as well as in the construction of modern cities. With the wide application of satellite remote sensing image data in various fields and the rapid development of remote sensing satellite technology, the wide application of Landsat series satellite images plays an important role in promoting the development of remote sensing technology. This paper takes Mohe city, Heilongjiang province as the research area, the land surface temperature (LST) retrieval was conducted from remote sensing Landsat8 images acquired on August 23, 2018 and August 30, 2018 based on radioactive transfer equation algorithm. The highest temperature in Mohe city was 36.5909˚C, the lowest temperature is 14.2906˚C and the average temperature was 23.2028˚C. The standard deviation is 1.5777. The temperature is low in the whole area. Areas with water and high vegetation coverage are low temperature areas, and the densely populated area of the city is the high temperature zone. This experiment is of great significance to study the local heat island effect and forecast and monitor forest fires.

1. 绪论

地表温度(Land surface temperature)是平衡地球系统水量与能量之间关系的一个重要因素，体现了地球表面物质能量变化与交换过程，是地–气系统相互作用过程中一个极其重要的参数 [1]，不仅对于全球性的气候变化来讲有着重要意义，而且在各种灾害的监测、预防方面，以及现代城市建设中也都起到很重要的作用 [2]。另外陆面温度也被广泛应用于日常生活中，主要常见于农业生产方面和研究环境之中，都需要测量地表温度作为必要的研究要素。要获取全球或区域尺度的地表温度的时空分布，只能依靠热红外遥感。伴随卫星遥感影像数据被广泛应用于众多领域之中，再加上近年来遥感卫星技术和热红外遥感技术的不断进步并迅速发展，实现了很方便的大范围地实时监测地表温度 [3]，在其中对遥感技术发展起到重要推动作用的要数Landsat系列卫星影像的广泛使用 [4]。当前，在众多研究学者不断推导以及大量的实验论证下，他们提供了多种不同的算法反演地表温度 [5]。其中多波段算法需要采用2个或2个以上的热红外波段来反演地表温度，可以减少一些假设的参数的使用，因此多波段算法反演精度较高 [6]，然而多波段算法相对要复杂些，以大量地气象资料作为支撑，但单波段算法依赖于一个热红外波段就足够可以反演地表温度，操作比较简单。本研究以黑龙江省漠河市为研究目标，借助两组卫星遥感数据，使用辐射传导方程法反演了漠河市的地表温度并分析了结果。

2. 研究区概况

漠河市在东经120˚12'到127˚00'和北纬50˚11'到53˚33'之间，是中国最北、纬度最高的县，坐落在黑龙江省的西北部，属于大兴安岭山脉北麓地段。寒温带大陆性季风气候是当地的主要气候类型，大陆和海洋季风相互影响，使得小范围内形成的气候千变万化，局部地区气候相差较大，冬季气候严寒干燥，温度常在零下30℃~40℃，最低温度曾达到零下52.3℃，冬日持续时间较长，主要受极地大陆气团控制；在副热带海洋气团的影响下，夏季降水丰富，温暖潮湿，但持续时间较短，只有半个月左右，而且昼长夜短，白天最长可以达到19个小时以上；由于冬夏季风交替控制形成的季风气候，春、秋两季气候变化较大，春季降水少，地表干旱，且经常发生大风天气，所以春季容易发生森林火灾，秋季降温迅速，霜降经常提早发生，使农作物受到冻害。漠河市处于我国的高纬度地区，是多年冻土区，由于温度低且持续时间久，冻土层最厚可达100 m以下，是我国多年冻土区人口较为密集的城市之一。冻土层一定程度上影响了当地野生动植物的生存生长，更适合耐寒，喜低温的动植物生存。漠河市存在热岛效应，而由冻土引发的工程和环境建设问题十分突出，且漠河市植被覆盖率较高，大面积的林区使得森林火灾的预防和监控问题尤为重要，解决这些问题都需要研究地表温度。

3. 数据与数据预处理

采用美国地质调查局(USGS)上下载的2018年8月23日行列号为(122, 23)标识为(LC81220232018235 LGN00)和8月30日行列号为(123, 23)标识为(LC81230232018242LGN00)的Landsat8卫星遥感影像数据，影像无云，质量较好。数据预处理包括辐射定标和大气校正。使用ENVI5.1的辐射定标工具Radiometric Calibration将两张影像的热红外数据和多光谱数据进行辐射定标，通过辐射定标后的影像DN值可靠。之后选用Flaash校正法，使用Atmospheric Correction工具对辐射定标后的数据进行大气校正，基本消除了对地物反射以及大气分子、气溶胶散射的影响 [7]。

4. 地表温度反演原理

传统获取地表温度数据的方法主要是通过收集研究区域地面的各个监测站点的监测数据，虽然精度相对较高，但是成本也相对高，收集数据耗费时间长，且由于监测站点分散，收集到的数据为离散点数据，要想计算得出大范围的连续时间内的地表温度只能使用离散点数据结果进行空间插值这一方法，但是这样计算出的地表温度值和实际结果相差比较大，且具有不确定性。

辐射传导方程法还有单窗算法和劈窗算法 [8] 等算法都可用于热红外遥感来反演地表温度。劈窗算法适合于两个热红外波段的数据，如MOAA-AVHRR和MODIS，大气校正法和单窗算法都需要大气实时剖面数据，而单窗算法适用于只有一个热红外波段的数据。

4.1. 单窗算法

覃志豪于2001年提出了单窗算法(MW算法)，是针对TM数据只有一个热红外波段来反演地表温度的一种算法。自此以后，又有许多的学者研究并证明，单窗算法不仅繁衍精度高，而且除了TM数据，同样可以运用在ETM+和Landsat8数据中。计算公式为：

$T_S=\left(a_6\left(1-C_6-D_6\right)+\left(b_6\left(1-C_6-D_6\right)+C_6+D_6\right)T_{sensor}-D_6{T}_6\right)/C_6$ (1)

公式中的T_S代表地表温度(K)，T_sensor表示式传感器上的亮度温度(K)，而T_a为大气平均温度(K)；a、b是参考系数，若地表温度处于0℃~70℃时，a = −67.355351，b = 0.458606；C、D是中间变量，公式如下所列：

$C=\epsilon \tau$ (2)

$D=\left(1-\tau \right)\left[1+\left(1-\epsilon \right)\tau \right]$ (3)

式中，ε为地表比辐射率，τ为地面到传感器的大气总透射率。因此单窗算法反演地表温度的关键是通过计算得到亮度温度T_sensor、大气平均作用温度T_a、地表比辐射率ε、大气透射率τ。

4.2. 劈窗算法

劈窗算法主要是针对NOAA/AVHRR的4和5通道即(10.3~11.3 μm)、(11.5~12.5 μm)设计的，它是利用方程组来反演地表温度，主旨是借助两个通道对水汽吸收和比辐射率存在的差异来完成，MODIS数据的31 (10.78~11.28 μm)、32 (11.77~12.27 μm)波段与NOAA/AVHRR这两个通道十分相似，所以劈窗算法也能够用在MODIS数据中，另外Landsat8遥感影像有着两个热红外波段，其波段信息与MODIS的31、32波段比较接近，因此劈窗算法也可应用于Landsat8数据地表温度反演，公式为：

$T_S=A_0+A_1{T}_i+A_2{T}_j$ (4)

其中，T_S为地表温度(K)，T_i和T_j分别是热通道i和j的亮度温度；A₀、A₁、A₂是分裂窗系数，由大气透过率和地表比辐射率等因子确定。

4.3. 辐射传导方程法

本次研究采用的方法为辐射传导方程法(radioactive transfer equation, RTE)，又称大气校正法，辐射传导方程法所需气象资料比较少，且只需要一个热红外波段，操作简单。

应用原理：卫星热红外传感器接收到的热红外辐射亮度值L_λ由大气上行辐射亮度(大气直接热辐射)L↑、大气下行辐射亮度(大气向地面的热辐射)L↓、地表的真实辐射亮度经过大气衰减后被卫星传感器接收的热辐射热量 [9] 这三部分组成。卫星传感器接收到的热红外辐射亮度值L_λ的表达式可以写为(辐射传输方程)：

$L_{\lambda }=\left[\epsilon \cdot B\left(T_S\right)+\left(1-\epsilon \right)L\downarrow \right]\cdot \tau +L\uparrow$ (5)

式子中，ε表示地表比辐射率、T_S代表地表真实温度(K)、B(T_S)是指用Planck函数所表示的黑体热辐射强度(W∙m⁻²∙sr⁻¹∙μm⁻¹)，τ代表大气在热红外波段透过率。B(T_S)代表温度为T的黑体在热红外波段的辐射亮度，可由公式(1)进行计算：

$B\left(T_S\right)=\left[L_{\lambda }-L\uparrow -\tau \cdot \left(1-\epsilon \right)L\downarrow \right]/\tau \epsilon$ (6)

从上可知此类算法需要τ、L↑、L↓三个参数，而在NASA提供的网站(http://aatmcorr.gsfc.nasa.gov/)中输入实时探测到的数据便可以计算生成。此实验应用的Landsat8卫星影像为格林尼治标准时间2018年8月23日2时43分生成的LC81220232018235LGN00这一影像，影像中心纬度为53.098 N，中心经度为123.198E，计算得到的B10波段的大气透射率τ为0.80，L↑为1.50 W∙m⁻²∙sr⁻¹∙μm⁻¹，L↓为2.51 W∙m⁻²∙sr⁻¹∙μm⁻¹；LC81230232018242LGN00成像时间为格林尼治标准时间2018年8月30日2时49分，影像中心纬度为53.098N，中心经度为121.653E，计算得到的B10波段的大气透射率τ为0.87，L↑为1.02 W∙m⁻²∙sr⁻¹∙μm⁻¹，L↓为1.72 W∙m⁻²∙sr⁻¹∙μm⁻¹。

地表温度T_S可由普朗克定律反函数获取：

$T_S=K_2/\ln \left(K_1/B\left(T_S\right)+1\right)$ (7)

对于TIRS Band10，K₁ = 774.8853，K₂ = 1321.0789，以上公式均使用BandMath工具计算实现。

处理流程如图1，将影像校正后做辐射定标和地表比辐射率计算，分别获得辐射亮度图和地表比辐射率图像，二者通过公式计算得到黑体辐射亮度图，最后由普朗克定律反函数得到地表温度。其中地表比辐射率的计算，由于地表物质的组成结构不同，不同地区的地表比辐射率不同，而地表结构比较复杂，对地表的分类方法不同会影响结果的精度。

5. 基于遥感影像的地表温度反演

5.1. 辐射定标

打开Landsat8数据，ENVI会按照波长分为五个数据集。计算辐射亮度值需要计算热红外波段的辐射定标，通过ENVI选用Radiometric Correction/Radiometric Calibration，之后设置定标类型为辐射亮度值(radiance)，输入各种参数执行定标处理，对第10波段进行辐射定标，获得辐射亮度图像。

5.2. 地表比辐射率计算

地表比辐射率是指在同一温度下，地表发射的辐射量比上黑体辐射量，地表比辐射率的大小与地表

Figure 1. Flow chart

图1. 流程图

特征属性有关。目前，测量地表比辐射率主要根据可见光和近红外光谱信息来估计，采用Van的经验公式法(JYGS)或混合像元分离法(HHXY)来计算地表比辐射率 [10]，虽然地表结构比较复杂，但仍然可以在大体上分成植被覆盖区域与裸地区域。本文采用混合像元法(HHXY)来计算。

已有研究发现，水体、城镇、植被和土壤的比辐射率依次是0.995、0.970、0.986、0.972 [11]，要计算这四种典型地物的地表比辐射率时，将地球表面依照归一化差分植被指数NDVI以及中国1:100,000的土地资源分类系统大概划分为水体、城镇、自然表面三种。水陆之间有着明确的界线，而且像元具有单一性，提取过程会较容易些，它的比辐射率定为0.995；城镇及植共同组成城镇用地，植被以及土壤共同组成自然地表 [12]，因此研究这两者是分别运用以下公式计算相对应地表比辐射率：

${\epsilon }_t=0.9608420+0.0860322Pv-0.0671580P{v}^2$ (8)

${\epsilon }_n=0.9643744+0.0614704Pv-0.0461286P{v}^2$ (9)

式子中，Pv为植被覆盖度，通常用以下公式计算：

$Pv=\left(NDVI-NDV{I}_{Soil}\right)/\left(NDV{I}_{Veg}-NDV{I}_{Soil}\right)$ (10)

上述公式中，NDVI是归一化植被指数，可应用于检测植被生长状态，植被覆盖度和消除部分辐射误差等，它能反映出植物冠层的背景影响，如土壤、潮湿地面、雪、枯叶、粗糙度等，且与植被覆盖有关。多种卫星遥感数据反演NDVI产品是地理国情检测云平台推出的生态环境类系列数据产品之一。如果计算植被覆盖度的话，需要利用下列公式从ENVI中提取归一化植被指数(NDVI)，即：

$NDVI=\left(NIR-R\right)/\left(NIR+R\right)$ (11)

上述公式中R代表可见光红光波段，是Landsat8中的第4波段，而NIR代表近红外波段，是Landsat8中的第5波段 [13]。

NDVI_Soil表示全部为裸土或者完全没有植被覆盖区域的NDVI值，而NDVI_Veg是指全部都有植被覆盖区域的NDVI值，纯植被像元的NDVI值。分别取经验值NDVI_Veg = 0.70、NDVI_Soil = 0.05，也就是说当某个像元的NDVI > 0.70时，Pv的取值是1；当NDVI < 0.05，Pv取值是0。

由上述可知，利用NDVI把地表分别看作水体、城镇、自然表面来计算归一化植被指数(NDVI)和植被覆盖度(Pv)也是可以的。当NDVI ≤ 0时，将地表看作为水体，地表比辐射率是0.995；当0 < NDVI < 0.7时，将地表看作为城镇，那么它的地表比辐射率是ε_t，运用公式(4)来计算；当NDVI ≥ 0.7时，则把地表看作为自然表面，此时它的地表比辐射率是ε_n，可用公式(5)进行计算。本次实验0 < NDVI < 0.7，即地表比辐射率使用公式(4)计算。在ENVI计算工具中输入公式：(b1 gt 0.7)*1 + (b1 lt 0.05)*0 + (b1 ge 0.05 and b1 le 0.7)*((b1 − 0.05)/(0.7 − 0.05))，表达式中的b1就代表NDVI，进行计算便可获得植被覆盖度图像。接着在计算工具中输入表达式：0.9608420 + 0.0860322b1 − 0.0671580b12，b1为刚刚获得的植被覆盖率，计算之后可以得到地表比辐射率。

5.3. 黑体辐射亮度计算

在上文中提到的NASA网站查询大气剖面数据，输入相关参数可得到大气剖面信息，LC81220232018235LGN00B10波段的大气透射率τ为0.80，L↑为1.50 W∙m⁻²∙sr⁻¹∙μm⁻¹，L↓为2.51 W∙m⁻²∙sr⁻¹∙μm⁻¹；LC81230232018242LGN00B10波段的大气透射率τ为0.87，L↑为1.02 W∙m⁻²∙sr⁻¹∙μm⁻¹，L↓为1.72 W∙m⁻²∙sr⁻¹∙μm⁻¹。根据公式(6)以及已求出的各个数据，在Bandmath中输入表达式(b2 − 1.50 − 0.80*(1 − b1)2.51)/(0.80b1)和(b2 − 1.02 − 0.87*(1 − b1)1.72)/(0.87b1)，其中b1指的是上文中求得的地表比辐射率，b2指的是最开始求得的第10波段的辐射亮度值，分别计算之后可以得到各自的黑体辐射亮度影像。

5.4. 地表温度反演结果与分析

根据公式(7)以及之前所求的黑体辐射亮度，在Bandmath输入表达式：(1321.08)/alog(774.89/b1 + 1) − 273,其中b1是黑体辐射亮度，计算即可得到地表摄氏温度，将两幅地表温度反演结果影像使用Mosaicking-Seamless Mosaic工具拼接在一起，然后使用ArcMap将漠河市边界矢量数据导入ENVI中，用Subset Data from ROIs来裁剪，得到漠河市地表温度数据图，如图2。

Figure 2. Surface temperature map of Mohe city

图2. 漠河市地表温度图

对本次Landsat8数据辐射传导方程法反演漠河市2018年8月地表温度的结果进行分析，最高温度为36.5909℃，最低温度为14.2905℃，平均温度为23.2028℃，标准差为1.5777，温度主要集中在20℃~23℃，与中国气象数据网的数据对比，反演精度较为良好。结合温度反演结果图可以看出，由于漠河市纬度较高，且地处林区，周边植被覆盖率较高，所以虽然是夏季八月份，地表温度并不高。温度较低的区域是河流水体，温度中等的区域也是面积最大的区域是林区，温度较高的区域是居民地，人类活动对我国高纬度多年冻土区的生态环境造成了显著影响，漠河市存在热岛效应，城镇区域地表温度明显高于全区平均值。

对于热岛效应的研究以及对森林火灾的预防监控等有待进一步研究。

6. 结语

本文以黑龙江省漠河市为研究对象，借助两幅覆盖了漠河市的Landsat8遥感影像，使用ENVI和ARCGIS软件，利用GIS以及遥感技术，对漠河市地表温度进行反演。使用ENVI5.3的工具计算出了研究区域的地表比辐射率，黑体辐射亮度等相关数值以及大气上行辐射亮度、大气下行辐射亮度和大气透射率等相关参数。通过与中国气象数据网的数据相比较，结果表明实验地表温度反演结果精度较好。对漠河市地表温度进行分析发现该研究地区存在热岛效应。为改善热岛效应，应该倡导居民低碳环保。地表温度对监测预防森林火灾，尤其是雷电火、自燃等自然火的发生有重要研究意义。此次实验采用的数据采集时间为夏季，存在水汽等影响，还需要在水汽影响较低的地区进行试验分析。这次地表温度反演的过程中还存在一些问题，比如地表比辐射率的计算存在一定误差，在未来的研究中应该更仔细地考虑地表比辐射率的精度问题。