一种融合多源卫星遥感的全天候地表温度近实时反演方法

技术领域

本发明属于卫星遥感地表温度技术领域，具体涉及一种融合多源卫星遥感的全天候地表温度近实时反演方法。

背景技术

地表温度(LST)是地表与大气相互作用过程中的关键物理量。掌握区域和全球尺度LST的时空分布，特别是全面、完整、连续的信息具有重要意义。随着遥感技术的不断进步和发展，各种平台的传感器的发射升空使快速、实时获取区域和全球LST成为可能。LST对于地-气辐射与平衡之间的地表发射率(LSE)研究具有重要意义，该领域已成为定量遥感领域的重要组成部分。使用光学遥感反演LST的算法较为成熟，产品精度高，分辨率基本可以满足实际应用的需求。但是，由于云覆盖的影响，极大地破坏了LST时空分布的完整性和连续性，严重阻碍和限制了LST产品的应用和发展。在区域和全球范围内对资源和环境进行动态监测(如气候变化监测、农业干旱预测等)需要高分辨率的全天候LST。在过去的5年中，集成卫星热红外遥感(TIR)和被动微波(MW)已被证明是估算全天候LST的一种实用可行的方法。但目前，所有方法都集中与历史存档数据的重建，无法满足对地表进行实时或近实时(NRT)监测的要求。同时，被动微波数据(如AMSR2)因卫星极轨运行方式导致轨道间隙区域数据缺失的问题也没有得到很好的解决，换言之，被动微波与热红外遥感相融合也不能实现真正意义上的全天候。

发明内容

针对上述问题，本发明提供了一种融合多源卫星遥感的全天候地表温度近实时反演方法。

本发明的技术方案包括下列步骤：

步骤1：获取目标区域的亮温数据，被动微波的亮温数据，以及目标区域中分辨率成像光谱仪的地表温度数据；

步骤2：构建微波辐射计的亮温数据与被动微波亮温数据之间的回归映射模型：

BT

其中，n表示目标时段的第任意天，d

步骤3：构建轨道间隙区域内的被动微波亮度温度的估算模型：

BT

其中，BT

以轨道间隙区域外的原始被动微波的亮温数据作为验证集，根据对被动微波亮度温度的估算精度对随机森林映射RF

基于筛选结果，根据BT

步骤4：基于目标区域的被动微波数据的微波亮温数据和中分辨率成像光谱仪的地表温度构建随机森林回归关系：

T

其中，t

步骤5：对目标时段的第n天，取指定周期的对应时刻与该第n天组成重构周期，基于所述重构周期内的无缝被动微波亮度温度进行全天候地表温度估计，得到所述重构周期内的全天候地表温度估计值；

所述全天候地表温度估计为：

对于所述重构周期内的任意天的任意时刻t

其中，T'

T'

T'

其中，BT

步骤6：获取所述重构周期内的已发布的中分辨率成像光谱仪的地表温度组成时间序列，并基于步骤5所获得的全天候地表温度估计值，对目标时段的第n天的全天候地表温度进行校正处理(即对系统误差进行校正)：

其中，T'

T'

本发明提供的技术方案至少带来如下有益效果：

首先实现了被动微波数据轨道间隙区域的缺失值填补，耗时短，结果精度高。集成热红外卫星遥感数据和被动微波数据生产全天候地表温度的方法，长期以来不能很好地解决轨道间隙区域测量值的获取问题，往往简单地使用三次样条插值来进行简单填补，结果的实际物理意义不明确。因此该方法的结果并不能实现真正意义上的“全天候”。其次，实现了指定距离(例如1km)全天候地表温度数据的近实时获取。一直以来，全天候地表温度生成(估算)方法仅在具有完整时间周期的历史数据集上进行了应用，这与其方法大多建立在地表温度在一个完整时段(例如年尺度/日尺度)是周期变化且具有规律的前提上有关。国际上主流的卫星热红外遥感地表温度产品(例如中分辨率成像光谱仪)，已经能够实现近实时获取(普遍滞后1.5-2天)，全天候地表温度估算方法在时效性上的短板极大地限制了其进一步的应用。本发明通过灵活重构年尺度周期，使用机器学习方法建立被动微波的亮温值与中分辨率成像光谱仪探测的亮温值之间的回归模型，然后将其延拓到相邻的包含最新时刻(目标点)上，得到与被动微波的亮温值和中分辨率成像光谱仪的亮温值同步时间的近实时1km全天候地表温度产品。本发明摆脱了传统集成方法对完整周期数据的依赖，能够满足不同应用需求、不同长度监测时段的数据生产需求，逐日产品的质量稳定，同传统的热红外卫星遥感地表温度产品具有良好的一致性，且经过实测地表温度的检验，精度在可接受范围，图像质量较高，能展现更多地表温度的空间细节信息。

附图说明

图1是本发明实施中，训练时段与预测时段的示意图。

图2是本发明实施中，包含实验区域的区位图。

图3是本发明实施例提供的融合多源卫星遥感的全天候地表温度近实时反演方法的处理过程示意图。

图4是本发明实施中，轨道间隙区域缺失值填补效果对比图。

图5是本发明实施填补后BT与原始AMSR2 BT的定量对比散点图。

图6是本发明实施处理得到的近实时全天候LST与MODIS LST的定量对比结果。

图7是本发明实施例中针对第一个目标时段的处理得到的全天候LST与MODIS LST的空间分布图。

图8是本发明实施例中针对第二个目标时段的处理得到的全天候LST与MODIS LST的空间分布图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

本发明实施例提供的融合多源卫星遥感的全天候地表温度近实时反演方法，首先基于中国风云气象卫星搭载的微波辐射计(MWRI)提供的微波亮温数据(BT)填补了AMSR2(Advanced Microwave Scanning Radiometer 2)微波亮温数据轨道间隙内的缺失值，然后建立空间无缝的AMSR2 BT和TIR LST(如MODIS)的年尺度模型，并将其应用于相邻年尺度时间段以获得最新目标时刻的LST。

本发明实施例提供的的融合多源卫星遥感的全天候地表温度近实时反演方法可以分为两个部分。在第一部分，在卫星搭载的微波辐射计(MWRI)和高级微波扫描辐射计(AMSR2)的微波亮温数据之间建立随机森林回归映射关系以获得近实时的空间无缝移动微波的微波亮温数据(MW BT)。第二步，通过建立空间无缝的高级微波扫描辐射计的微波亮温数据(AMSR2BT)和TIR LST的年度尺度模型并将其应用于目标时间，最终得到近实时的全天候LST。其中，卫星搭载的微波辐射计可以是中国风云气象卫星搭载的微波辐射计。

虽然卫星搭载的微波辐射计和中分辨率成像光谱仪(Moderate-resolutionImaging Spectroradiometer，MODIS)具有不同的扫描路径和近似的观测时间。然而，本发明的发明人基于统计分析发现MWRI与AMSR2之间仍存在36-120分钟的观测时间差。因此，MWRI的微波亮温数据不能直接用于填补AMSR2的轨道间隙区域。但是，本发明的发明人基于统计分析还发现，实验区域的MWRI BT与AMSR2 BT之间具有高度的相关性，

中国风云气象卫星搭载的微波辐射计MWRI和Aqua MODIS具有不同的扫描路径和近似的观测时间。然而，根据我们前期的统计结果，发现MWRI与AMSR2之间仍存在36-120分钟的观测时间差。因此，MWRI BT不能直接用于填补AMSR2的轨道间隙区域。根据我们的统计分析，实验区域的MWRI BT与AMSR2 BT之间具有高度的相关性，相关度R>0.85。

式(1)中，i为AMSR2轨道间隙区域外有效值的序号，BT

因此，可以通过建立MWRI BT与AMSR2 BT之间的回归映射关系来估算AMSR2轨道间隙区域内的亮温值。同其他机器学习方法相比，随机森林(RF)可以更充分地表征回归参数与目标特征之间的关系，并以更高的准确性和泛化能力构建回归映射模型。因此，本发明实施例使用RF建立MWRI每个通道的亮温值与AMSR2 BT之间的回归关系，然后在AMSR2的轨道间隙区域内进行应用。

在实验区域(目标区域)的目标时段的第n天，在AMSR2的轨道间隙区域外建立AMSR2BT和MWRI BT的RF回归映射模型：

BT

式(2)中，f和p分别为观测通道的频率和极化方式，BT

然后，将式(2)构建的RF关系应用到AMSR2的轨道间隙区域内，实现轨道间隙区域的缺失值填补：

BT

式(3)中，BT

以上为第一部分，最终得到空间无缝的AMSR2 BT。

然后建立空间无缝AMSR2 BT和MODIS LST的年度尺度模型并将其应用于目标时间，最终得到近实时的全天候LST。

首先，将研究时段分为训练时段与预测时段。训练时段为需要近实时估算的日期的前一个年尺度周期(图1所示，reconstruction cycle-1)，预测时段即为需要近实时估算的目标时段，具体时间长度不定。

根据之前基于历史数据重建的全天候LST估算研究表明，在两个相邻的年尺度周期内，TIR LST和MW BT之间的映射关系是相似的。因此，根据从前一个年尺度周期(训练时段)获得的映射关系，可以应用到预测时段内以得到近实时的全天候LST。通过建立线性回归模型能够较好地刻画TIR LST与MW BT之间的映射关系，并且该方法估算的LST精度较高。本发明方法考虑到映射关系的稳定性和描述非线性关系的能力，使用RF回归方法代替线性或非线性回归方法。建立晴空指定距离(例如1km)的MODIS LST像元与多通道MW BT之间的RF模型：

式(4)中，t

然后，将公式(4)构建的RF映射(RF

式中，T'

MODIS LST作为应用较为广泛的热红外地表温度产品，本发明实施例提供的实时反演方法需将估算的近实时全天候LST“系统性”地纠正到MODIS LST水平上(即尽可能减少与MODIS LST的系统偏差)，校正后的近实时全天候LST如公式(5)所示：

公式(6)中，T'

示例性的，本发明实施例实际实施区域选择黑河流域及其周边区域，其区位图如图2所示。黑河是中国第二大内陆河，总流域面积为14.3万km

本实施例中具体采用的数据集包括：

1)MODIS产品：MODIS逐日地表温度产品(MYD11A1,1km)，逐日MODIS云量产品(MOD06_L2,5km)，逐日地表反照率产品(MCD43A3，1km)，逐日归一化雪指数(NDSI)产品(MOD10A1，500m)，16天合成归一化植被指数(NDVI)(MOD12A2,1km)；

2)被动微波数据集：逐日近实时AMSR2 BT三级产品，分辨率0.1°，中国风云气象卫星微波辐射计(MWRI)提供的近实时三级亮温产品(FY-3B/3D MWRI BT,0.1°)；

3)地面实测数据集：黑河流域生态—水文过程综合遥感观测联合试验提供的2014年6个自动气象站的长波辐射数据集。

参见图3，具体实现对全天候地表温度近实时反演处理的过程包括：

1.数据预处理与时空匹配。

针对的数据对象包括近实时数据和长时序数据。

首先，使用IDL(Interface Description Language)语言和MRT(MultispectralRefraction Topography)工具对MYD11A1进行批量的格式转换、投影变换、拼接和裁剪工作，得到包含黑河流域及其周边地区的MODIS LST数据集)；对于MYD11A1数据，使用质量控制(QC)图层，剔除低质量的MODIS LST像元。

其次，使用IDL语言对AMSR2 BT数据进行批量的格式转换、投影以及裁剪(剔除与MODIS LST观测时间差大于指定值(例如5min)的像元)，然后分别读取AMSR2 BT和MODISLST的经纬度信息并保存为数值矩阵。在MATLAB中根据读取经纬度信息，进行空间匹配，生成对应位置关系的查找表。对于MWRI BT数据，由于存在轨道漂移，需要在MATLAB平台对MWRI BT与AMSR2 BT逐日进行空间位置的匹配，同时剔除观测时间差超过指定值(例如5min)的像元；

再其次，分别将5km的MODIS云量产品(MOD06_L2)和500m的逐日归一化雪指数(NDSI)产品MOD10A1重采样到1km分辨率，以匹配MODIS LST数据的空间分辨率。16天1km分辨率的NDVI数据则进行时间插值(三次样条插值)，以匹配MODIS LST逐日的时间分辨率。

最后，在MATLAB上使用时间滤波器方法填补NDSI因云覆盖导致的数据缺失部分。

即数据预处理的目的是为实现不同数据源之间的空间分辨率匹配和时间分辨率，其中数据预处理包括：格式转换、投影变换、拼接和裁剪。

2.逐日的AMSR2 BT轨道间隙区域缺失数据填补。

2-1)根据公式(2)，在MATLAB上对目标时段某一天的MWRI和AMSR2 BT建立轨道间隙区域外的基于RF的回归关系。

2-2)根据公式(2)得到的回归结果，利用公式(3)估算轨道间隙区域内的AMSR2 BT值。

2-3)根据提前设定的回归树数量，将轨道间隙区域外的原始AMSR2 BT作为验证集，利用公式(3)估算AMSR2 BT值，与验证集相比较，选择精度最好的回归树的结果作为最终的轨道间隙区域填补结果。

3.近实时全天候地表温度估算。

3-1)对于目标时段的第n天，在其前一年或多年具备完整历史周期的AMSR2 BT与MODIS LST之间建立随机森林回归关系(公式(4))，目的是获取足够的当地LST周期变化的时间规律信息。

3-2)对于目标时段的第n天，取前一个年尺度周期的对应时刻与其组成重构周期，该周期内的空间无缝AMSR2 BT带入公式(5)中，得到该周期内的全天候地表温度估计值。

3-3)在步骤3-2)中的重构周期内，获取已发布的MODIS LST组成时间序列(不一定完整)，结合公式(5)得到的初始估计值，带入到公式(6)中，得到与MODIS LST一致的、目标时段第n天的全天候LST。

按照当前采用的数据源来说，本发明实施例可以做到滞后两天发布，符合近实时监测的要求。

进一步地，具体的结果分析如下：

由于实施地面观测数据获取存很大困难，本发明实施例提供的实时反演方法使用黑河流域生态—水文过程综合遥感观测联合试验提供的2014年6个自动气象站的长波辐射数据集，以验证1km近实时全天候LST。本发明实施例提供的实时反演方法选择的评估指标包括平均偏差误差(MBE)，标准偏差(STD)，均方根误差(RMSE)和确定系数(R

图4以10GHz和89GHz水平通道为例，展示了实验区域内AMSR2 BT和填补后BT的空间分布。填补后的BT不仅补全了轨道间隙区域内AMSR2的缺失值，而且在空间分布和数值范围上与原始AMSR2 BT保持了高度一致性。此外，填补后的AMSR2 BT图像上，原轨道间隙区域的边缘过渡平滑，没有出现空间不连续现象。这表明本发明方法中的间隙填充算法具有很高的可靠性。

图5展示了填补后BT和AMSR2 BT之间的散点图。MBE为0.03-0.25K。RMSE为1.02-1.58K，R

图6展示了在晴空条件下，逐日1km的近实时全天候LST和MODIS LST产品的定量对比结果。近实时LST与MODIS LST相比有一定的系统误差，2014和2020年MODIS LST的MBE在白天为0.12-0.19K，夜间为-0.07-0.06K。STD小于1.8K，R

另外，MODIS LST和近实时全天候LST的图像质量也进行了比较。图7展示了2014年在四个不同季节(DOY 15、105、196和288)白天MODIS LST和近实时全天候LST的空间格局。无论哪个季节，近实时全天候LST在空间上都是无缝的，并且在空间格局和有效值范围与原始MODIS LST高度一致。同时，近实时全天候LST图像中没有马赛克或空间不连续性，这表明本发明方法的结果图像质量良好。更进一步，图8显示了2020年9月30日白天MODIS LST和近实时全天候LST的空间格局。近实时全天候LST是空间无缝的，并且在空间格局和有效值范围与原始MODIS LST高度一致。云覆盖范围内缺失的LST得到了很好的恢复，没有出现异常值。同时，在近实时全天候LST图像中没有马赛克或空间不连续性，这表明公式(4)中的RF映射在克服粗分辨率BT与其1km子像元(即，MODIS像元)LST之间的尺度不匹配方面确实有效。

表1和表2展示了近实时全天候LST和站点(Site)实测LST的验证结果。

表1

表2

在表1和表2中，“ARO”、“DAM”、“HZS”、“GB”、“SSW”、“HZZ”分别表示不同的站点名称，各站点的位置示意图可参考图2。

综合白天和夜间的验证结果，在晴朗的天空和非晴朗条件下，本发明实施例提供的实时反演方法得到的逐日近实时全天候LST通常低于站点实测LST。这种现象的主要原因是MYD11A1产品在干旱和半干旱地区存在冷偏差现象。在验证结果中，夜间的MBE介于-1.77K至-0.68K之间，而白天MBE介于-2.58-0.17K之间。除了HZZ站外，其他站点的RMSE和MBE波动均未超过1K，表明本发明实施例提供的实时反演方法在全天候条件下具有令人满意的性能且不确定性有限。

在所有六个站点中，GB，SSW，HZS和HZZ的验证结果均比其他站点差，RMSE在2.5K和4.94K之间。可能的原因有三个方面1)由被动微波热采样深度引起的地下热辐射；2)站点所在区域的空间异质性；3)MODIS LST本身存在的误差，主要是干旱和半干旱地区的温度值低估。

综上，可以表明本发明实施例提供的实时反演方法在不同的时间尺度上具有良好的适用性和鲁棒性。近实时全天候LST在晴空和非晴空条件下均具有稳定且令人满意的精度。尽管本发明实施例中使用AMSR2和MODIS作为数据源，但它有被移植到其他传感器平台的潜力。例如，本发明实施例提供的实时反演方法可以应用于集成地球静止气象卫星的TIRLST(例如ESA的地球温度产品)和MW BT，然后可以生成近实时的高分辨率的全天候LST，这对于区域或全球气象观测将具有很高的价值(天气预报，干旱和火灾监测以及径流预测等)。即本发明实施例提供的实时反演方法所获得的处理结果具有良好的精度，RMSE范围为2.46-4.62K。将有助于近实时全天候LST产品的生产，能够在时效性上对标国际上主流的热红外遥感地表温度产品。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

以上所述的仅是本发明的一些实施方式。对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。