基于DTW的HJ-1B卫星地表温度产品质量检验

张琪，李大成,*，历华，段四波，张霄羽，韩启金，高海霞

1. 太原理工大学，太原 030024 2. 中国科学院 空天信息创新研究院，北京 100094 3. 中国农业科学研究院，北京 100081 4. 山西大学，太原 030006 5. 中国资源卫星应用中心，北京 100094

1 引言

地表温度(land surface temperature，LST)作为衡量地球表面能量平衡状态的重要参数之一,对地气循环系统研究、气候预报、生态环境、林业养殖以及地球化学等领域有深远意义。因此，对地表温度的空间和时间上的精确测量对于气候变化、生态环境监测等各个方面的应用都有重要意义。随着遥感技术的快速发展，热红外遥感为高精度地表温度的获取提供了有效的途径[1]。

为了获得具有较高时空分辨率的地表温度，国内外学者针对不同的遥感数据提出了多种地表温度反演算法[2-4]。单通道算法是针对一个热红外通道提出的地表温度反演算法[5-6]，辐射传输方程法[7]是常用的一种算法。该算法通过输入大气廓线数据并计算大气参数来计算表面温度。由于实时的大气参数一般是利用大气辐射传输模型进行模拟得到的，获取困难，因此适用性不高。为了解决这一问题，文献[8]对辐射传输方程进行简化，提出了一种针对Landsat-TM5数据的单窗算法(MW算法)，该算法不需要大气廓线数据，只需要根据大气透过率、大气水分含量和平均大气温度的经验线性关系对大气进行校正，就能计算得到地表温度。文献[9]提出了一种只需要输入大气水汽含量就可以计算地表温度的普适性单通道算法(JM&S算法)，该方法通过分析大气水汽含量与大气透过率、大气上行辐射和大气下行辐射之间的经验关系来消除大气影响，计算简便。

中国于2008年发射了环境一号A/B卫星(HJ-1A/B)，其中环境一号B星(HJ-1B)搭载两台宽幅多光谱(CCD)相机和一台红外相机(IRS)，重访周期为4 d。CCD相机可提供空间分辨率为30 m的4个波段，IRS相机可提供空间分辨率为300 m的热红外波段[10-11]。自卫星发射以来，国内学者针对环境小卫星的单个热红外波段的数据特点展开了一系列研究工作[12-18]。文献[12]针对环境一号卫星B星进行数据计算，对MW算法和JM&S算法进行对比，结果表明两种算法总体接近，JM&S算法精度更高。文献[13]使用3种单通道算法对HJ-1B数据进行地表温度反演，并对精度进行分析比较，结果表明JM&S算法精度最高。

在获得卫星影像地表温度反演结果后，需要对其真实性进行验证[19]。对于LST产品的真实性检验有4种方法：基于地面的验证、基于卫星LST产品的验证、基于辐射的验证和基于时间序列的比较。基于地面的验证即将反演结果与地面测量的温度进行直接比较，简单易行，常用于验证MODIS的地表温度产品[20]。基于卫星LST产品的验证方法是将两种卫星得到的LST产品进行比较[21]。该方法适用于面积较大的区域，但是会受到卫星传感器的过境时间、观测视角以及不同卫星数据空间分辨率等因素的影响[22]。基于辐射的验证不需要原位测量，而是需要准确估算通道的比表面发射率值以及卫星过境时间的大气水分含量和大气温度[23]。基于时间序列的比较方法通过选取较长时段的遥感地表温度,分析其自身表现的时序特征来发现传感器在运行中可能存在的问题，对地表温度数据进行异常值处理[24]。该方法仅仅能够得到相对精度,对于真实性的评价存在明显不足。

随着中国遥感技术的发展，针对国产卫星的数据反演及真实性检验工作也不断深入。大多数研究对于反演结果的验证较少，通常只是利用同一时刻的卫星LST产品进行验证，缺少对反演结果的定量评价，无法精确应用于其他研究。因此，在HJ-1B卫星方面地表温度的反演及产品的真实性检验尚需要完善。本文以黑河流域多个类型的下垫面为研究区，采用普适性单通道算法对覆盖研究区的2012-2018年清晰无云的HJ-1B遥感影像进行地表温度反演，并利用对应时刻的MOD11A1地表温度和地面站点数据对反演结果进行验证。由于验证数据在时序上比反演数据多，并不完全匹配，为了充分利用验证数据，引入时序相似性检测思想，提出利用动态时间规整(dynamic time warping，DTW)将时间序列进行拉伸或收缩匹配来计算两个时序之间的相似性。

2 研究区及数据预处理

研究区位于甘肃省黑河流域中上游，东经97.1°～102.0°，北纬37.7°～42.7°，该地区气候干旱，降水较少。本研究区内用于验证的地面站点有阿柔超级站(ArouCJZ)、黄藏寺站(HZS)、垭口站(YK)、花寨子荒漠站(HZZ)、神沙窝沙漠站(SSW)、张掖湿地站(SD)和黑河遥感站(HHYG)，站点的地表覆盖类型主要包括高寒草地、戈壁、湿地、荒漠、沙漠、农田等，具体情况如表1所示。

表1 研究区站点情况

本文所选用的遥感数据是HJ-1B数据与MODIS数据，其中HJ-1B数据包括多光谱数据和热红外数据，时间2012-2018年，共60个时相，均可覆盖研究区并清晰无云。MODIS数据包括用于对反演结果进行验证的MOD11A1地表温度产品、用于计算大气水汽含量的MOD021KM以及用于对MOD021KM进行几何校正的MOD03数据。

在进行地表温度反演之前首先对HJ-1B数据进行几何校正，以消除影像因传感器的高度、姿态等因素影响造成的几何畸变，并以Landsat-8遥感影像对其进行影像配准。随后利用定标系数对HJ-1B卫星的CCD影像进行辐射定标，定标公式如下：

式中：L为定标后的辐亮度，单位为W·m-2·sr-1·μm-1；a为绝对定标系数增益；L0为绝对定标系数偏移量，具体值可以在影像的头文件中获取。

卫星传感器在成像过程中会受到大气、光照等因素的影响，因此需要对其进行大气校正，本文采用基于MORTRAN4的FLAASH(fast line-of-sight atmospheric analysis of spectral hypercubes)大气校正模型对定标后的HJ-1B卫星的CCD数据进行大气校正。

为了进行下一步计算及验证，对于MOD021KM数据通过MOD03数据进行几何校正，随后将MOD021KM和MOD11A1数据的投影坐标转换成研究区所在的投影，并将像元尺度重采样到30 m。

本文所研究的地面站点数据来源于黑河流域遥感观测自动气象站，站点每30 s观测一次，每10 min输出一次结果，该结果为这10 min内的平均值。这些站点的观测设备不完全相同，主要有CNR1四分量辐射仪和SI-111红外辐射计两种观测设备，两种设备的计算公式分别为：

(1)

B(TS)=[B(Tr)-(1-ε)Lsky]/ε

(2)

式中：TS为地表温度；B为普朗克公式；F↑为上行短波辐射强度；F↓为下行长波辐射强度；σ为斯蒂芬玻尔兹曼常数，σ=5.67×10-8W/(m2·K4)；εb为宽波段发射率，根据ASTER GED的窄波段发射率或波普库数据可以计算得到；Tr为SI-111红外辐射计观测得到的站点实测温度；Lsky为大气下行辐射强度；ε为SI-111红外辐射计实测的发射率数据。

3 研究方法

本文采用普适性单通道算法(JM&S算法)对HJ-1B卫星数据进行地表温度反演，并利用黑河流域遥感观测站实测数据和MODIS地表温度产品对反演结果进行验证，基于动态时间规整(dynamic time warping，DTW)对反演结果和验证数据的时序相似性进行评价，具体流程如图1所示。

图1 研究技术路线Fig.1 Research technology roadmap

3.1 基于JM&S的HJ-1B LST产品反演

普适性单通道算法(generalized single-channel method，SC)又称JM&S算法，是由Jiménez-Muňoz和Sobrino在2003年针对只有一个热红外通道的遥感数据提出的地表温度反演算法，该算法是通过对普朗克函数做一阶泰勒级数展开所得到的，具体计算公式如下：

TS=γ[ε-1(ψ1Lsensor+ψ2)+ψ3]+δ

δ=-γLsensor+Tsensor

式中：Lsensor为卫星传感器测得的辐射强度，Tsensor为亮温，单位为K；λ为有效波长；C1和C2为普朗克函数常量，C1=1.191 04×108W·μm4·m-2·sr-1，C2=1.438 77×104μm·K；大气函数ψ1、ψ2和ψ3是关于大气水分含量ω的函数。对于大气函数的取值，文献[12]根据HJ-1B卫星的特征对函数进行重新拟合，得到的计算公式如下：

ψ1=0.041 2ω2+0.093 6ω+0.985 6

ψ2=-0.717 4ω2-0.881 2ω+0.394 1

ψ3=0.263 9ω2+0.649 9ω+0.470 3

地物比辐射率ε和大气水分含量ω是JM&S算法的关键参数。大气水分含量反映了大气水分在空间上的不同分布，在地表温度反演中有重要作用。考虑到数据获取的便捷，MODIS数据的通道比值法在计算大气水分含量方面应用广泛。因此本文采用与HJ-1B卫星同时相的MODIS数据的第2波段(大气窗口波段)和第19波段(水汽强烈吸收波段)计算大气水分含量，计算公式如下：

式中：ρ2和ρ19分别为MODIS第2波段和第19波段的表观反射率；α和β为常量，其值依据地表不同而有所不同。本研究区内地表覆盖类型包含植被、裸土及其他混合地表，因此本文采用混合地表类型，对于混合型地表，α=0.02，β=0.651。

按照地物类型的区别，影像分类一般包括植被、水体、裸土及建筑物表面这4类。在估算地物比辐射率方面，常用NDVI阈值法，该方法对地物在多光谱影像不同波段的光谱信息结合，通过NDVI的大小来区分不同地物。针对环境卫星，根据地物的波谱特征，在裸土地区，地物比辐射率取值0.972，纯植被地区的地物比辐射率取为0.99，由于水体接近于黑体，水体的比辐射率取为0.995[25]。混合像元区域地物类型复杂，可由如下公式计算地物比辐射率：

ε=εvPv+εs(1-Pv)+dε

dε=(1-εs)(1-Pv)Fεv

式中：Pv为植被覆盖度，dε为地表几何像元内部的影响；NDVI为归一化植被指数。对于环境卫星，ρ3和ρ4分别为环境一号卫星第3和第4波段的反射率；εv为完全植被覆盖时的比辐射率；εs为完全裸土覆盖时的比辐射率；F为地形因子，F=0.55。

3.2 基于DTW的星地时序非对应LST产品检验

现有的对地表温度反演结果的验证都是基于对应实点的验证，无论是直接验证还是间接验证都无法对时序进行评定。用于验证的数据往往在时序上连续且比反演数据的时序更多，当验证数据与反演数据时间不匹配时，为了充分利用地面站点数据，引入时序相似性思想，通过DTW对两种数据的时序相似性进行质量评价。

DTW通过将时间序列拉伸或收缩匹配来计算两个时间序列之间的相似性，该方法相较于传统欧式距离测度检验方法，克服了时间扭曲问题，在两个时间序列之间寻找一个灵活的时间匹配路径使两个序列对齐，具有更加广泛的实用性。

计算DTW首先需要构建时间序列。分别构建反演结果时序{H}={H1,H2,…,Hm}和检验数据时序{L}={L1,L2,…,Ln}，m和n分别为两个时序的长度。其次，计算两组时序之间任意两个元素之间的差值，生成距离矩阵D(Li,Hj)=(Li-Hj)2，其中1≤i，j≤m。然后根据距离矩阵寻找最短路径，即最佳规整路径。该过程是基于动态规划的思想，假如当前节点是D(Li,Hj)，那么下一个节点必须是在D(Li+1,Hj)，D(Li,Hj+1)，D(Li+1,Hj+1)三者之间选择，并且路径必须是最短的。最后，通过如下递归式计算DTW累计距离：

DTW(i,j)=D(Li,Hj)+

4 结果与分析

本研究选取覆盖黑河流域研究区的2012-2018年清晰无云的HJ-1B卫星影像，采用JM&S的普适性单通道算法对其热红外波段进行地表温度反演，并采用基于地面和基于卫星地表温度产品的两种方法对结果进行精度验证和分析。分别统计各个站点处的实测温度与MODIS地表温度，将反演结果与这两者进行对比，计算平均偏差Bias、标准差STD、均方根误差RMSE、动态时间规整DTW和决定系数R2。

图2为反演温度在7个站点处与地面温度及站点处MODIS地表温度的散点图，图2(a)(c)(e)(g)(i)(k)(m)分别为7个站点的反演温度和地面温度的相关性散点图，图2(b)(d)(f)(h)(j)(l)(n)分别为7个站点的反演温度和MOD11A1地表温度的相关性散点图，图上决定系数R2是在对反演温度和验证温度做线性回归模拟时计算得到的评价系数，无单位。图3为7个站点的反演温度与2种验证温度之间偏差值的折线图。表2为两种验证方式在各个站点处的评价指标值。

对于下垫面为沙漠、荒漠的神沙窝沙漠站和花寨子荒漠站，反演精度最高。图3(a)(b)分别为花寨子荒漠站和神沙窝沙漠站的偏差值曲线，可以看出这两个站点的反演温度和地面温度之间的差值较小，与MODIS地表温度的偏差值稍大些但比其他站点小。从时序上来看没有明显受到时间的影响，可能是由于两个站点的下垫面为荒漠和沙漠，空间异质性较小，地表辐射率随时间变化较小，因此反演精度较高。

对于下垫面为湿地、雪地的张掖湿地站和哑口站，反演精度相对其他站点较低。垭口站的Bias相对其他站点较高，与地面站点相比平均低估3.33 K，与MODIS地表温度相比平均低估4.13 K。张掖湿地站的反演结果与地面站点相比平均低估2.24 K，与MODIS地表温度相比平均低估3.12 K。但这两个站点与验证数据的决定系数R2均在0.9以上，表明反演结果与验证数据的相关性并不差。造成这一结果的原因可能是对于湿地和雪地区域，地表覆盖类型复杂，NDVI阈值法不能准确地计算出该区域的地物比辐射率。

对于下垫面为植被或农作物的3个站点，从图3(c)(d)(e)可看出，阿柔超级站在7～10月反演结果与验证数据差别较大，可能是因为阿柔超级站地表覆盖类型为高寒草地，在夏季温度上升植被生长，对传感器的观测造成影响。黄藏寺站点在3～5月的反演结果与验证数据相差较大，在其他日期相差较小。这可能是因为该站点的地表覆盖类型为小麦，冬小麦的生长季节在春季，成长季地表辐射率变化较大，因此反演结果的精度相对较低。黑河遥感站在夏季的反演结果与验证数据相差比其他季节大，该站点的地表覆盖类型为玉米，生长季节在6～10月，植物生长对地物比辐射率的计算造成影响，因此在植被生长季时反演效果不如其他季节好。

对于7个不同下垫面的站点研究区，反演得到的LST与温度的偏差值均小于其与MODIS地表温度的偏差值，RMSE和STD的表现与Bias一致。图3的曲线图也可看出偏差值大多为负值，表明反演结果大多低于验证数据。基于地面温度的验证显示，除了垭口站的平均偏差最大为-3.33 K外，其他6个站点的与实测温度的平均偏差均在2 K以内。反演结果与地面温度的RMSE最大为垭口站的3.37 K，除垭口站之外其他站点的RMSE最大都不超过3 K。所有站点的反演结果与地面温度的STD都在2 K以内。

与MODIS地表温度相比，除了垭口站和张掖湿地站的平均偏差稍高外，其他站点的偏差值都在3 K以内，RMSE指数最小可达1.26 K(花寨子荒漠站)。STD最大为阿柔超级站的2.04 K，最小为花寨子荒漠站的0.67 K。反演结果与MODIS地表温度的差异可能是由于HJ-1B卫星在本研究区的过境时间并不固定，而MODIS卫星传感器的一景影像的覆盖区域较大，两种卫星的过境时间在不同日期相差不同，因此HJ-1B地表温度的反演结果与MODIS的地表温度数据差别并不稳定。另外，HJ-1B CCD卫星的原始空间分辨率为300 m，MODIS地表温度产品的原始分辨率为1 000 m，由于最终使用的是重采样得到的结果，因此插值计算也可能对结果造成影响。

反演结果与地面温度相比，花寨子荒漠站和神沙窝沙漠站的DTW最小，为0.001 7，表明这两个站点的反演温度与地面温度的相似性最高，这与Bias、RMSE和STD表现一致。从表2其他站点的验证指标也可看出，当Bias、RMSE和STD较小时，DTW时序检验结果也相对较小，表明反演温度与验证温度的相似性较高；反之当Bias、RMSE和STD较大时，DTW时序检验结果也相对较大。在评价时序预测精度上，DTW与已有的评价指标(如Bias、RMSE和STD)具有较高的一致性。

图2 HJ-1B反演温度与地面站点、MODIS温度产品散点图Fig.2 Scatter plots of HJ-1B retrieval LST and ground temperature,and of HJ-1B retrieval LST and MODIS LST

续图2Fig.2 Continued

续图2Fig.2 Continued

图3 各站点反演温度与验证温度偏差值Fig.3 The bias between retrieved temperature and verified temperature at each site

续图3 Fig.3 Continued

表2 站点处反演温度评价指标

5 结论

本文使用Jiménez-Muňoz和Sobrino在2003年提出的普适性单通道算法对覆盖黑河流域研究区的2012-2018年的HJ-1B卫星数据进行地表温度反演，并基于地面实测站点数据和MODIS地表温度进行了精度分析和验证，利用DTW对时序不匹配的数据进行了时序相似性评价。通过对结果的分析和讨论，可以得出以下结论:

1)在本研究区内，对于地表覆盖类型为沙漠、荒漠的站点，由于地表辐射率变化不大，反演结果与地面站点相差最小(最小可达0.015 K)，精度最高；对于地表覆盖类型为植被或农作物的区域，地表辐射率随时间发生变化，反演结果精度较好，平均偏差均在2 K以内，相关性较高；对于地表覆盖类型复杂的区域，可能是由于地物比辐射率不能准确计算，反演结果与站点数据偏差较大。

2)HJ-1B反演结果与MODIS地表温度产品对比时，对于所有站点，可能是由于环境卫星和MODIS的过境时间不同，反演结果与验证数据的差值均大于其与地面观测温度的差值，但相差不足1 K，可以认为两种验证方法结果一致。

3)从反演结果来看，在本研究区内，下垫面为植被和农作物的站点在植被生长季时反演精度较低，在其他季节反演精度较高。因此在对植被覆盖区域进行地表温度反演预测时，应该选择植被非生长季时相的影像。

4)对于时序不匹配的验证数据和反演数据，利用DTW对其进行动态规整以评价两组数据的时序相似性，DTW的评价结果与现有的评价指标表现一致，为遥感影像的时序检验提供一种新的途径。