基于天宫二号数据的北京市地表温度反演研究

白沁灵孟 丹*邢 婧

(1.资源环境与地理信息系统北京市重点实验室，北京100048；2.首都师范大学资源环境与旅游学院，北京100048)

1 引言

地表温度体现地球表面物质的能量变化与交换过程，是地球科学研究的一个重要指标。地表温度的区域分布在气候变化、植被生态、环境监测和城市热岛等研究领域都有着重要的应用价值[1]。

随着遥感卫星技术的发展，越来越多的卫星搭载了波段和分辨率等参数不一的热红外传感器，国内外学者针对不同的传感器提出多种地表温度反演方法，这些算法大致可以分为5类：单通道算法、多通道算法、多角度算法、多时相算法和高光谱反演算法[2]。其中多通道算法用2个或2个以上的热红外波段反演地表温度，减少了一些参数的假设，因此反演精度较高[3]。多通道算法中，针对2个热红外通道数据的劈窗算法计算简单、方法成熟，是目前应用最为广泛的温度反演算法之一。劈窗算法是由McMillin[4]提出的一种根据AVHRR数据的2个热红外波段来反演海面温度的方法，后被推广应用于陆地表面温度反演。该算法利用2个相邻通道对大气吸收作用的不同，通过2个通道辐射亮度的各种线性组合来剔除大气的影响，反演地表温度。近年来，国内外学者基于 ASTER[5]、MODIS[6]、Landsat-TIRS[7]、Himawari 8 AHI[8]等数据开展了一系列的劈窗算法研究和应用。

天宫二号宽波段数据具有2个相邻的热红外波段(T1∶8.125 ～8.825 μm，T2 ∶8.925 ～9.275 μm)。当前劈窗算法多是针对10～14 μm波长区间的数据，很少涉及8～10 μm区间内的热红外通道。孟鹏等[9]证明了将8～10 μm波长区间的数据用于劈窗算法反演地表温度具有一定的可行性。基于以上研究，本文在传统劈窗算法的基础上，调整算法参数，对北京市的地表温度进行反演，并利用MODIS地表温度产品对反演结果进行了精度验证，最后统计分析北京市不同土地利用类型的地表温度情况。

2 研究区与数据

2.1 研究区概况

北京市地处华北平原，位于北纬39°36′到41°03′、东经 115°25′到 117°30′之间，以山地和平原地形为主，地势西北高、东南低。三面环山，东南部是一片缓缓向渤海倾斜的平原，海拔介于10～2157 m之间，总面积达 16 410.54 km2，属于暖温带半湿润大陆性季风气候，夏季高温多雨，冬季寒冷少雨，春、秋短促，年均降水量 595 mm，其中75%集中在夏季，年平均气温11～13℃。

2.2 遥感数据源

天宫二号空间实验室搭载的宽波段成像仪是新一代宽波段、宽视场和图谱合一的光学遥感器，国内外首次在单台仪器上实现了可见光、短波红外和热红外多光谱大视场全推扫成像的组合集成功能，对地表覆盖信息、地表温度等研究具有重要意义[10]。

本研究使用天宫二号宽波段成像光谱仪的可见光近红外图像和热红外图像对地表温度进行反演，数据来源于载人航天空间应用数据推广平台(http：//www.msadc.cn/)。本次研究选取2018年4月27日经过系统辐射和几何校正后的两景天宫二号宽波段L2级产品数据，图像质量良好，无云。地表温度反演结果验证借助MODIS地表温度数据产品 MOD11 A(https：//modis.gsfc.nasa.gov/)。土地利用数据来源于中国科学院资源环境科学数据中心(http：//www.resdc.cn/)。

3 研究方法

本文在传统劈窗算法的基础上，针对天宫二号热红外通道的特点，重新定义算法参数。

首先对辐射传输方程中的Planck函数进行线性简化，建立2个热红外通道的亮度温度与辐射亮度之间的关系；接着利用天宫二号数据的第2和第4波段反演大气水汽数据，并计算其大气透过率；然后根据ASTER光谱数据库提供的地物波谱数据以及利用可见光数据计算的归一化植被指数NDVI(Normalized Difference Vegetation Index)计算地表比辐射率；最后联立2个热红外波段的辐射传输方程，即可求得地表温度。技术流程图如图1所示。

图1 技术流程图Fig.1 Technology flowchart

3.1 图像预处理

获取的天宫二号数据是经过系统辐射校正、几何校正后的(入瞳处)辐亮度产品，DN(Digital Number)值是经过放大的辐亮度，经过缩小还原后得到天宫二号数据的辐射亮度。采用Landsat 8影像作为基准参考影像，对天宫二号影像数据进行几何精纠正，校正精度控制在1个像元内。采用FLAASH大气校正模块对天宫二号数据进行大气校正，设置天宫二号的轨道参数，并导入其波谱响应函数。最后利用行政边界数据对影像进行裁切，得到天宫二号宽波段成像仪的北京市遥感影像(图2)。

图2 2018年4月27日北京市天宫二号影像图(红：V9，绿：V10，蓝：V12)Fig.2 TG-2 image of Beijing on April 27， 2018(R：V9， G：V10， B：V12)

3.2 劈窗算法的推导

地表温度反演通过建立辐射传输方程来实现，热辐射传输方程是热红外遥感和地表温度反演的基础[5]，地表辐射经过大气到达传感器的过程可以简单的描述为式(1)：

式中Bi(Ti)是遥感器所接收的辐射亮度；Ts是地表温度；τi(θ)是通道i在遥感器视角为θ下从地面到遥感器的大气透过率；εi是比辐射率；Bi(Ts)是在地表温度为Ts时的黑体辐射亮度；Ti是星上亮度温度；是大气下行辐射；是大气上行辐射。覃志豪等[11]对有详细的推导，有近似解见式(2)：

式中Ta为大气向上平均温度，为向下平均作用温度。

热辐射传输方程涉及Planck函数的计算，由于Planck函数十分复杂，所以在进行地表温度反演研究时，通常对其进行线性简化。覃志豪[12]研究AVHRR数据反演地表温度和 Offer[7]研究Landsat 8数据反演地表温度研究时，均对Planck函数进行泰勒展开，见式(3)：

将线性简化的 Planck函数式(3)带入式(1)，对于天宫二号的2个热红外波段辐射传输方程简化为式(4)：

为了方便计算，定义系数为式(5)：

这样，式(4)可简写为式(6)：

联立2个方程组，求解出地表温度为式(7)：

ai、aj、bi、bj为线性简化的 Planck 函数中的系数。这样，只需确定Planck函数参数、亮度温度、大气透过率、地表比辐射率即可得到地表温度Ts。

3.3 Planck函数线性简化

本文利用Planck函数分别对天宫二号宽波段数据的2个热红外波段的辐射亮度在温度273～322 K区间[13]内的变化关系进行计算，如图3所示，辐射亮度随温度的变化近似线性关系。进行线性回归得到方程式(8)。

图3 天宫二号2个热红外波段的辐射亮度与温度关系图Fig.3 Relationship between radiation intensity and temperature in two thermal infrared bands of TG-2

3.4 亮度温度的计算

天宫二号宽波段数据为经辐射校正和系统几何校正后的产品，利用Planck函数求解亮度温度[14]，计算公式如式(9)所示：

式中Ii为辐射亮度，Ki，1＝2hc2/λ5i，Ki，2＝hc/(kλi)，h为普朗克常数，约为 6.62606896×10-34J˙s，c为光速，为 2.99792×108m˙s-2，k 为玻尔兹曼常数，为 1.3806505×10-23J˙K-1，λi为第i通道内的中心波长。

天宫二号宽波段数据在热红外数据T1波段的中心波长λ1＝8.475 μm，在热红外数据T2波段的中心波长λ2＝9.1 μm，根据Planck公式得出天宫二号宽波段数据的Ki，1和Ki，2。 对于热红外波段T1，K1，1＝2724.1 W˙m-2˙sr-1˙μm-1，K1，2＝1697.3 K，热红外波段T2，K2，1＝1908.6 W˙m-2˙sr-1˙μm-1，K2，2＝1580.7 K。

3.5 大气水汽含量与大气透过率估算

大气水汽含量与大气透过率是地表温度反演过程中的关键参数。在MODIS的大气水汽反演中利用MODIS数据的第2波段(0.841～0.876 μm)和第 19 波段(0.915～0.965 μm)计算其大气水汽含量[15]。天宫二号宽波段数据的第4波段(0.845～0.885 μm)和第 2波段(0.930～0.950 μm)的波长范围与MODIS数据的第2波段和第19波段波长范围相近，故天宫二号宽波段数据反演大气水汽含量的估算公式如式(10)所示：

式中ω是大气水汽含量(g/cm-2)，取α＝0.02和β＝0.6321[15]；ρ2和ρ4是天宫二号数据第2和第4波段的地面反射率。

张晓等[16]模拟出的8～9.3 μm波长范围内2个热红外波段的夏季大气透过率与水汽含量的关系，利用该关系计算大气透过率。计算公式如式(11)所示：

式中τ1、τ2为天宫数据2个热红外波段的大气透过率。

3.6 比辐射率的估算

地表比辐射率主要由地物表面的结构和波长范围决定，天宫二号热红外数据的波段范围为8～9.3 μm。在劈窗算法地表温度反演中，通常假定地表辐射率已知，估算方法主要有经验值法和混合像元法。本文采用覃志豪等[17]提出的混合像元法对地表比辐射进行估计。混合像元法中，将遥感影像分为水体、自然表面和城镇表面，自然表面看作植被和裸土的混合，城镇表面看作城镇和植被的混合，其计算公式如式(12)所示：

式中εsurface为自然表面的比辐射率，εbuilding为城镇表面的比辐射率，Pv为植被覆盖度，εiv、εis、εim为植被、裸土、建筑在热红外波段内的比辐射率，Rv、Rs、Rm为植被、裸土、建筑的温度比率；dε为比辐射率的修正项，在地表相对平整的情况下，这一项的值就很小[18]，研究区域属于平原地区，因此取 dε＝0；

混合像元法中提出了植被、裸土和建筑表面的温度比率Rv、Rs、Rm的估算公式如式(13)所示：

其中Pv为植被覆盖度，由公式Pv＝NDVINDVIs/NDVIv-NDVIs求出，NDVIv和NDVIs分别表示纯植被像元和纯裸土像元的归一化植被指数NDVI值。当NDVI≤NDVIs时，植被覆盖度Pv＝0；当NDVI≥NDVIv时，植被覆盖度Pv＝1。 由于受到季节、地表湿度、大气条件等因素的影响，NDVIv和NDVIs并不是一个固定值[19]，本文通过统计天宫二号数据影像NDVI直方图，取5%和95%的累积百分比作为置信区间，取NDVI的上下限阈值分别为研究区的NDVIv和NDVIs来进行植被覆盖度的近似估计。

根据ASTER光谱数据库提供的典型地物波谱数据计算典型地物的比辐射率，水体、植被、裸地和建筑在8～9.3 μm波段范围内的反射率变化如图4所示。由于该波谱库的数据为反射率数据，根据基尔霍夫定律，处于某一温度的物质比辐射率与吸收率值相同。除水体外，对于一般地物都可不考虑透射，即比辐射率＝吸收率＝1-反射率，水体的比辐射率在热红外波段内变化较小，取0.995[17]。选择几种典型的地物计算其平均比辐射率，典型地物水体、植被、裸地和建筑在天宫二号数据的2个热红外波段T1和T2范围内的比辐射率分别为：0.995、0.9802、0.9457、0.9378 和0.995、0.9761、0.9432、0.9408。

综上，自然表面和城镇表面在天宫二号的热红外波段范围内的比辐射率由公式(14)进行估算。

4 结果分析

4.1 地表温度反演结果验证

图4 典型地物在8～9.4 μm波段范围内的反射率Fig.4 Reflectivity of typical ground objects in the range of 8～ 9.4 μm

采用上述劈窗算法，得到天宫二号地表温度反演结果(图5)，选取2018年4月27日的MODIS地表温度产品(图6)与天宫二号宽波段数据反演结果进行一致性分析。天宫二号热红外波段的分辨率为400 m，而MODIS地表温度产品空间分辨率为1 km，将MODIS温度产品空间分辨率重采样为400 m，进而与同分辨率的反演结果进行对比。从图5、图6可以看出，天宫二号宽波段数据反演得出的北京市地表温度最低温度为292.2 K，低温部分主要出现在密云水库，最高温度为310.4 K，高温区出现在中心城区和延庆地区；MODIS温度产品的最低温度为286.8 K，最高温度为309.0 K，高温区和低温区与天宫二号宽波段数据反演的结果一致。从空间上看，天宫二号宽波段数据反演得出的地表温度结果与MODIS地表温度产品的空间分布情况基本一致，北京地区的地表温度空间分布趋势表现出明显的城市热岛效应，城镇温度高于自然地表的温度，水体的地表温度最低，从图5可以看出密云水库明显边界。

图5 基于天宫二号数据反演的北京市地表温度Fig.5 Beijing LST retrieval from TG-2 data

图6 基于MODIS数据反演的北京市地表温度Fig.6 Beijing LST retrieval from MODIS data

图7 西南-东北方向地表温度剖面图Fig.7 LST profile in SW-NE direction

选取研究区地表温度起伏较多的西南至东北方向对天宫二号和MODIS反演的地表温度做剖面分析，以剖面线的西南端起点为原点，剖面线的位置见图5、图6，分析结果如图7所示。剖面线上，天宫二号数据反演的地表温度整体比MODIS地表温度产品数据高；地表温度剖面起伏情况基本一致；因为天宫数据的原始分辨率较高，其地表温度剖面的起伏较多；高温中心均出现在市区和城镇地区，从市区和城镇地区到郊区地表温度逐渐降低；人类活动频繁地带温度较高，自然表面的温度较低，水体的温度最低。

从研究区随机选取1000个采样点，将天宫二号宽波段数据反演的地表温度结果与MODIS地表温度产品进行对比分析，二者的误差均值为3.5 K，误差的标准差为1.4 K；对其进行线性回归分析，可以看出，二者有着较高的相关性，其相关系数达到0.79，均方根误差为1.28 K(图8)。由于天宫二号数据与MODIS数据过境时间、空间分别率等存在差别，导致两者反演的地表温度存在差异，具体原因还有待下一步工作的研究。

图8 基于MODIS与天宫二号数据反演的地表温度散点图Fig.8 LST scatter plot based on MODIS and TG-2 data inversion

4.2 北京市不同土地利用类型地表温度分析

结合北京市2018年的土地利用图(图9)，统计各土地利用类型地表温度的均值、中位数、最大值、最小值、标准差，如表1所示。

不同土地的利用类型存在差异。其中，平均地表温度较高的是建设用地，为305.40 K，这是由于地表温度的高低主要取决于地面物体的热学性质，如热容量、热惯量、热传导率以及地表潜热通量[20]。城市人工构筑物吸热快而热容量小，升温快，同时，城市人为的热排放也是造成建设用地的地表温度较高的原因之一。水域的平均地表温度最低为301.1 K，原因是水体的热容量大，升温慢，且受热后会以水分蒸发的形式把热量散发出来，降低自身的温度。耕地的平均地表温度仅次于建设用地，原因是其植被覆盖度较低，植被蒸腾作用弱。不同土地利用类型的地表温度的平均值从高到低的顺序为：建设用地＞耕地＞草地＞林地＞未利用土地＞水域。

利用不同土地利用类型所对应的地表温度绘制箱线图，如图10所示。箱线图可以看出反演的温度值是否稳定，而中位数可以对比不同类型数据的趋势。从箱线图可以看出：水域地表温度数据值的箱体较长，说明该地表温度数据变化较大，较不稳定，可能是因为天宫二号宽波段数据的热红外波段的400 m的分辨率容易形成混合像元，且由于河流两岸多为耕地、草地、建设用地等，吸收更多的太阳辐射，混合像元会表现出更高的表面温度，因此会导致河流表面温度过高；而其他土地利用类型的箱体相对较短，说明其地表温度较为稳定。

图9 2018年北京市土地利用图Fig.9 Land use map of Beijing in 2018

表1 不同土地利用类型地表温度情况统计(单位：K)Table 1 Statistics of LST for different Land use(Unit：K)

5 结论

1)天宫二号宽波段数据波段设置可以满足地表温度反演的参数需求。其中的第2波段和第4波段，可以反演出大气水汽含量，进而计算出大气透过率。

图10 北京市不同土地利用类型地表温度箱线图Fig.10 LST Box Map of different Land use types in Beijing

2)基于天宫二号数据反演出的结果与MODIS地表温度产品具有较高的一致性，该参数可推广应用于热红外波段在8～10 μm波长的传感器。验证了天宫二号宽波段数据可以利用改进参数后的劈窗算法进行地表温度反演，天宫二号为地表温度反演提供了新的数据源。

3)北京地区不同土地利用类型地表温度不同，地表温度情况从高到低的顺序为：建设用地＞耕地＞草地＞林地＞未利用土地＞水域。

致谢：感谢载人航天工程提供天宫二号宽波段成像仪数据产品。