CN114564767A - 一种基于太阳-云-卫星观测几何的云下地表温度估算方法 - Google Patents

Abstract

本发明属于遥感应用技术领域，公开了一种考虑太阳‑云‑卫星观测几何的云下地表温度估算方法，具体涉及了一种基于集成学习模型的FY‑3D MERSI‑II云下地表温度估算算法，云下地表温度的估算主要是针对FY‑3D MERSI‑II图像的考虑由于太阳‑地表‑卫星的几何观测结构，将云下像元分为云遮挡与云覆盖像元，并针对两部分像元不同的温度成因机制利用集成学习模型进行地表温度估算。估算过程中，结合了MODTRAN辐射传输模型，综合了太阳辐射、大气和地表和降水气象条件等相关的多重地表温度驱动因子，建立了一套实现自动化运算云下地表温度的流程方法，实现了FY‑3D MERSI‑II图像自动全天候地表温度图像估算，提高FY‑3D MERSI‑II高分辨率地表温度数据适用性。

Description

一种基于太阳-云-卫星观测几何的云下地表温度估算方法

技术领域

本发明属于遥感应用技术领域，具体涉及一种基于太阳-云-卫星观测几何的云下地表温度估算方法。

背景技术

地表温度(LST)是地表-大气相互作用在区域和全球尺度上能量交换结果的一种综合表征，广泛地应用于地表能量平衡、旱灾监测、农田土壤水分管理、农作物估产、数值天气预报和气候变化等领域中。遥感观测具有空间连续、周期短且成本低的特点，因此是获取区域以至全球尺度LST时空分布信息的有效技术手段。在实际发明与应用中，由于测量的便捷性和成熟的晴空LST反演算法，地表温度数据仍主要通过热红外(Thermal Infrared,TIR)遥感观测数据获取。目前基于热红外遥感的MODIS、GLASS、FY3C/VIRR(可见光红外扫描辐射计)等地表温度产品，已广泛应用于变化检测，干旱监测和气候变化发明。然而，TIRLST在有云情况下不可用，TIR信号不能穿透云层到达地表，热红外传感器所观测的云像元热辐射主要来自云层顶部，而不是云层下地表的热辐射。热红外遥感影像中云像元的温度实际是云顶温度，而不是云下的地表温度。云覆盖像元地表温度信息的缺失在很大程度上破坏了区域地表温度时空分布的完整性和连续性，严重阻碍和限制了地表温度产品的应用和发展。许多遥感应用中，如气候变化调查、蒸散量估算、城市热岛等，经常遇到因无法估算云像元地表温度而不能对发明区进行全面有效监测的难题。云覆盖导致的数据缺失是当前地表温度数据应用的重要障碍。全面掌握全天候LST空间格局分布和时间演变规律，对于地表资源环境监测、地气***能量平衡和生态***等相关科学领域发明都具有非常重要的实际意义。

近年来，集成学习(Ensemble Learning)成为机器学习发展的新趋势。通过将结合一系列的弱学习器(weak learner)组合成精度更高和泛化能力更强的新模型，可以改善原算法功能单一，易于出现过拟合的缺陷。随机森林(RF)和梯度增强回归树(GBRT)和自回归极端梯度上升算法(XGBoost)是近年来发展的新兴集成学习模型，不仅对不同场景有较好的适用性，且具有出众的准确率。其中，GBRT和XGBoost使用boosting的集成决策树算法，而随机森林(RF)是使用bagging的集成决策树算法。这三种集成算法已经在生物学、医学、地学、交通、能源、农业工程和工业等各个领域得到广泛应用，并取得了比较理想的分类或回归效果。云下地表温度估算问题属于多元回归问题，因此应用这三种集成回归算法进行FY-3D MERSI-II云下地表温度的模拟运算发明。

风云三号(FY-3)系列卫星是我国第二代极轨气象卫星，FY-3D(Fengyun-3D)是第四颗星，2017年11月15日在太原卫星发射中心发射，其轨道类型为近极地太阳同步轨道，共携带中分辨率光谱成像仪(MERSI-II)、微波成像仪(MWRI)、微波温度计(MWHS-2)等10个观测仪器。中分辨率光谱成像仪II(Medium Resolution Spectral Imager II,MERSI-II)是该卫星的主要传感器之一，其观测能力有较大提升，可以媲美美国新发射的联合极轨卫星。为充分发挥我国气象卫星数据在国内外基础数据资源和技术的优势，提升国产气象卫星的综合应用价值和服务能力，国产气象卫星数据的基础算法发明还需要进一步加强。

遥感影像上云空间分布与太阳高度角和卫星观测角度相关(王斐,2017；WANG etal.,2017)。当太阳入射角与卫星观测角度不一致时，卫星观测的云和实际位置的云产生位置偏差(SHI et al.,2017)。太阳入射方向和传感器观测方向的相对角度越大，这种偏移越大。具体而言，在实际应用中，遥感数据通过云检测获取云位置，来判定晴空和云的分布情况，确定像元适合于晴空还是多云天气算法。但是遥感影像检测的云像元，只是真实云团在传感器观测方向上的投影，而并不是对应云团在地表的正射投影位置，也并不完全对应于云对地表辐射的影响区域。真实受云影响的地表区域对应于整个云团在太阳入射角度方向上的形成的地面投影区，其中可视部分即为遥感影像上的阴影区。太阳-云-卫星几何关系对地表参数的这种影响，称为SCSG(Solar cloud satellite geometry)几何效应(WANG etal.,2017；WANG et al.,2019)。目前云下地表温度估算发明中往往忽视这些影响，地表温度和辐射通量的太阳-云-卫星几何效应(SCSG)问题还没有得到足够重视。

基于以上发明背景，本次发明将针对FY-3D MERSI-II数据进行云下地表温度估算发明。首先基于两因素***窗算法(Two Factor SplitWindowAlgorithm,TFSWA)进行晴空地表温度反演，考虑地表温度和辐射通量的太阳-云-卫星几何效应，基于集成学习模型发展适用于FY-3D MERSI-II数据的云下地表温度估算算法。最终，构建FY-3DMERSI-II全天候地表温度数据，推动热红外遥感的进一步发展，提升国产卫星的综合效益与服务能力。

总体而言，本发明所发展的云下地表温度估算方法，对不同地表性质和气候区具有较好的适用性，可以为极轨卫星数据的云下地表温度估算提供方法参考。

发明内容

本发明属于遥感应用技术领域，公开了一种基于太阳-云-卫星几何机构用于云下地表温度估算新方法。为了地表温度全天候遥感估算，把像元分为晴空像元和有云像元，并进一步把有云像元分为云遮挡像元和云覆盖像元。考虑地表温度的太阳-云-卫星几何效应，将遥感影像上的云像元分为云覆盖(地表入射太阳辐射受到云影响，B区)和云遮挡(地表入射太阳辐射未受云影响，C区)两种类型，分别进行云下LST估算。其中C区为太阳直射区域，理论晴空地表温度代表该区域的云像元地表温度；云覆盖区(B区)是云在太阳观测方向的真实投影，太阳辐射受到云遮挡影响。图2显示的是太阳-云-卫星观测几何关系示意图。云下地表温度具有多因素耦合复杂性，集成回归模型对于非线性与非平稳过程特征挖掘具有优越性能。本发明给出了基于集成回归学习模型的云下地表温度估算的新方法，主要包括以下步骤：

步骤一、遥感影像辐射定标和几何校正

辐射定标是定量遥感应用的必要步骤，将传感器原始观测记录的无量纲DN值，转换为有实际物理意义的光谱反射率或辐亮度数据。采集L1级卫星图像数据,计算反射通道星上反射率、热发射通道亮度温度，并进行几何校正工作。

步骤二、晴空地表温度反演

将普朗克函数线性展开，经过一系列推导提出两因素***窗算法(TFSWA)，使用两个亮度温度通道、相应的波段平均地表发射率和大气透过率的一系列线性组合来估算LST，如下所示：

Ts＝A0+A1T₂₄-A2T₂₅ (1-1)

A0＝a₂₄[D₂₅(1-C₂₄-D₂₄)/(C₂₄D₂₅-C₂₅D₂₄)]-a₂₅[D₂₄(1-C₂₅-D₂₅)/(C₂₄D₂₅-C₂₅D₂₄)] (1-2)

A1＝1+D₂₄/(C₂₄D₂₅-C₂₅D₂₄)+b₂₄[D₂₅(1-C₂₄-D₂₄)/(C₂₄D₂₅-C₂₅D₂₄)] (1-3)

A2＝D₂₄/(C₂₄D₂₅-C₂₅D₂₄)+b₂₅[D₂₄(1-C₂₅-D₂₅)/(C₂₄D₂₅-C₂₅D₂₄)] (1-4)

C₂₄＝ε₂₄τ_24n (1-5)

D₂₄＝[1-τ_24n][1+(1-ε₂₄)τ_24n] (1-6)

C₂₅＝ε₂₅τ_25n (1-7)

D₂₅＝[1-τ_25n][1+(1-ε₂₅)τ_25n] (1-8)

式中，Ts为地表温度(K)；A0、A1、A2为中间参数，T_i为亮温，τin为大气透过率，i为热红外波段，εi为热红外波段相应的地表发射率，a₂₄、a₂₅、b₂₄、b₂₅为常数。

步骤三、晴空、云阴影、云遮挡与云覆盖区检测

针对云像元进行地表示温度估算，需要根据云检测产品区分云和晴空像元。针对云下像元对应的地表太阳辐射是否受到云削弱作用影响采用不同的方法估算云下地表温度，需要区分云遮挡和云覆盖区。因此需要检测云阴影、云遮挡与云覆盖区对应的具***置。在此过程中，云的正射校正是正确表征太阳-云-卫星几何关系的关键步骤(SHI etal.,2017)。根据太阳-云-卫星几何关系，通过云层高度、太阳和卫星的天顶角和方位角以及云层在遥感影像上的位置作为输入数据，就可以直接确定云的实际位置，依据正射位置可以推导出地表阴影空间位置(WANG et al.,2017)。角度和云高由几何定位文件和云高产品导出。使用公式(1-9)确定云的正射位置，然后使用公式(1-10)估计对应云阴影的正射投影位置(WANG et al.,2017；WANG et al.,2019)。

式中，(X_cld，Y_cld)为云在地表的正射投影；(X_shw，Y_shw)为云影在地表上的投影位置；(x_map,y_map)是遥感影像中的云位置即云团沿卫星观测方向的投影位置；而H是地面以上的云顶高度，由MERSI-II的云顶高度(云高产品)减去地面高度(DEM产品)计算。θ_v和φ_v是卫星过境时刻传感器观测天顶角和观测方位角，θ_s和φ_s分别是传感器观测时太阳的天顶角和方位角。偏移位置计算时，方位角需要转换为-180°-180°之间，东西南北方向分别是：90°，-90°，180°，0°。

步骤四、短波净辐射反演及其观测几何效应校正

在地球表面常被云层所覆盖，是全球或区域尺度天气或气候变化的主要因素之一。云对太阳辐射有透射、吸收和反射作用，自身发射长波辐射。云控制着太阳入射辐射和大气长波辐射。因此，云对地面辐射通量收支，地表温度变化和地表能量平衡等有决定性的影响。JIN et al.(2000)基于地表能量平衡模型(Surface Energy Balance,SEB)提出了邻近像元插值算法(NP)用于估算云下像元地表温度。该算法的基本假设是云覆盖像元与时空邻近的相似晴空像元，主要由太阳辐射差异引起地表温度差。通过建立地表长/短波辐射、潜/感热通量与地表温度的线性关系估算云下像元地表温度。本发明用短波净辐射作为云下地表温度估算的辐射特征变量。

(1)短波净辐射反演

本发明采用参数化算法反演短波净辐射，计算过程中只要输入太阳天顶角、大气水汽、相对方位角和大气水气含量以及大气层顶反照率作为输入数据反演短波净辐射，该方法具有计算简单、对云反演依赖性小的优点。

短波净辐射反演首先将窄带反射率ρ_i转换为TOA短波段宽带反照率γ，然后用参数化方法得到经验系数来计算NSSR。基于MODTRAN模型可以建立窄波段反射率与宽波段反照率的线性转换公式：

γ＝b₀+b₁ρ₁+b₂ρ₂+b₃ρ₃+b₄ρ₄+b₅ρ₅+b₆ρ₆+b₇ρ₇ (1-11)

式中，γ为星上短波宽波段反照率，ρ_i为MERSI-II波段i的星上窄波段反射率，b_i为MERSI-II波段i对应的拟合参数。反照率随方位角的变化小于天顶角的变化，因此，选择范围从0°～60°、60°～120°和120°～180°的三个粗略间隔，分别表示正向、侧向和反向的近似反射。

对于给定的观测天顶角(VZA)和太阳天顶角(SZA)对，使用公式(5-20)来建立获得的系数b₀～b₇和VZA之间的关系

式中，c_1i-c_4i分别是给定SZA＝0°、10°、20°、30°、40°、50°、60°和70°的常数。此处使用的不同SZAs、VZAs和相对方位角的转换参数。利用MODTRAN5.2以典型大气条件和地表边界参数为输入，通过辐射传输模型建立模拟数据库。根据模拟数据库，建立晴空和云下两种条件下，不同观测天顶角、太阳高度角和相对方位角对应的参数查找表，通过检索查找表线性插值得到星上宽波段反照率。

短波净辐射表示为如下公式：

NSSR＝a_sE₀μ/d² (1-13)

式中，E₀是星上太阳辐照度(1357W/m²)，d是地球-太阳距离。两个参数都是以天文单位表示，a_s为地表吸收系数。

地表吸收系数a_s可以表示为大气层顶宽波段反照率r的方程：

a_s(μ，ω，γ)＝α′-β′γ (1-14)

方程式(1-14)中的截距α′斜率β′分别表示为:

α′＝1-α₁μ^-1-a₂μ^-x-[1-exp(-μ)](a₃+a₄ω^y)μ^-1 (1-15)

式中，a₁-a₇是各种表面的常数，μ是SZA的余弦，x，y和z的值为0.5。

(2)短波净辐射几何效应校正

基于地表能量平衡模型估算云下地表温度，短波净辐射是关键输入参数。太阳-云-卫星的几何结构对地表辐射的影响具有普适性，短波净辐射数据受SCSG几何效应影响，有明显的高估或者低估现象。云的分布位置由于SCSG效应影响在遥感影像上出现一定的偏移。NSSR几何效应校正步骤为：首先主要利用方程1-9和1-10逐像元确定云及其阴影范围，并标定每个云像元与其阴影投影的位置关系。对于阴影区(A区)像元和云覆盖区(B区)云像元的NSSR被指定为其对应的云投影区NSSR。首先对于图像中的所有云覆盖区像元，逐像元迭代执行上述过程。对于云遮挡区(C区)，这些像元可以被太阳照射到，因此这些像元的NSSR应该与近距离内的晴空像元非常相似。在本文中，这些云遮挡区像元NSSR用邻近晴空像元插值计算或HANTS算法填补，最终生成经过几何效应校正的短波辐射数据NSSR_{corrected-cld}。

步骤五、特征变量选取与预处理

(1)特征变量选取

温度特征因子是训练数据及回归重建云像元LST的必要输入数据，用于实现整个重建过程。地表温度受到地理位置、植被、地形、微地貌、地表有效辐射、气象信息等的协同影响。

本发明分别从地表性质、太阳辐射、大气条件和传感器几方面选取云下地表温度回归运算的特征变量。所有与地表温度相关的因素中，高程是首选重要影响因素，海拔与温度相关性是可量化的参数，高程对地表温度变化的影响在高原和山区尤为显著，海拔降低时，对辐射能量分配上的影响减弱。DEM、坡度和坡向选取为地表地形表征参数。植被选择性地吸收和反射太阳辐射，调节潜热和显热交换来影响地表热辐射信息，植被参数必不可少，选取NVI、EVI和LAI作为植被参数。此外，LST日变化曲线与地表吸收和产热的短波净辐射通量有关，日地表温度循环与地表太阳净辐射(短波净辐射)之间的密切关系也已被确定。因此，太阳辐射方面，考虑了地表短波净辐射和太阳高度角。根据短波净辐射反演过程可知，地表短波净辐射反演的输入数据包括太阳常数、太阳高度角和方位角、大气水汽(WVC)和传感器观测角度等，说明该特征变量是地理位置、太阳辐射、传感器观测条件和大气条件的综合表征。表1给出了各特征变量的时间分辨率、单位等详细信息。

表1云下地表温度估算特征变量

(2)NVI、EVI和LAI数据的谐波优化处理

NVI、EVI和LAI是指示地表植被覆盖特征的重要变量,因遥感数据观测时大气条件、气溶胶和传感器故障等因素的影响,在MERSI-II LAI数据产品中存在一定的数据缺失和低质量数据，在云下地表温度估算中，严重影响了模型训练应用。为了获取较高精度的模式模拟结果,利用HANTS算法对NVI、EVI和LAI数据的缺失像元进行填补，将受云污染或其他因素引起的低质量数据进行平滑和滤波处理，生成连续完整的长时间序列高质量的训练数据。

(3)特征变量共线性分析

为了选取高质量的特征变量，需要检查变量之间是否存在共线性问题。基于Seaborn统计分析可视化模块计算了特征变量之间的PEARSON特征相关矩阵热图，去掉相关性较大的变量。

步骤六、模型建立与参数优化

(1)训练数据集的确定

为了准确地表示LST与特征变量之间的复杂关系，构建各特征变量到地表温度的拟合回归模型，选取高质量的晴空和阴影数据为训练数据集,将云像元地表温度作为目标变量，建立回归模型。基于Python的第三方模块Scikit-learn和XGBoost，利用RF、GBRT和XGBoost三种回归模型对晴空和阴影像元的地表温度数据集进行学习训练，利用train_test_split函数，将样本数据随机抽取划分，80％训练分区和20％测试分区两个部分，分别应用于训练模型和测试模型。

(2)参数设置与优化

为了防止过拟合现象，调整三个模型中模型回归树的个数(max-depth)参数，通过10折交叉验证查看训练数据的得分。为了查看模型性能，针对不同的训练样本查看其训练时间。

步骤七、基于集成学习模型针对晴空和阴影区进行训练，建立特征变量与LST之间的回归预测模型

考虑云下地表温度具有多因素耦合的特性，集成回归算法对于非线性变换过程和非平稳特性的特征挖掘具有优良性能，在充分考虑地表温度和辐射的太阳-云-卫星几何效应的基础上，采用RF、GBRT和XGBoost三种集成回归算法，针对晴空区和阴影区，以目标变量(FY-3D MERSI-II地表温度)和预处理后的特征变量，进行特征训练模拟建立LST回归预测模型。同时基于Python第三方扩展包scikit-learn特征变量重要性计算方法(PermutationImportance)，计算各特征变量的重要性得分和排序。

步骤八、输入待计算图像，利用回归模型进行云遮挡和云覆盖区云下LST估算

采用上一步骤建立的云下LST回归模型对发明区内云遮挡和云覆盖区所有像元的地表温度进行预测，并利用CLDAS同化LST数据对估算的云下地表温度进行精度评价。

相较于针对其他的云下地表温度估算的方法，其优势在于：

1、集成回归算法进行云下地表温度估算，计算速度较快，提升了数据处理的效率。

2、该算法不会受到下垫面类型、云覆盖范围大小等因素的影响，提高了云像元LST估算的适用性。

3、采用三种决策树类型的集成回归模型进行云下地表温度估算，能够提高模型拟合的精确度，并且通过参数的优化，可以防止过拟合现象。

附图说明

图1考虑地表温度和辐射SCSG效应的云下地表温度集成回归算法流程图

图2太阳-云-卫星观测几何关系示意图

图3 FY-3D MERSI-II云阴影检测及云遮挡和云覆盖区像元统计结果

图4 FY-3D MERSI-II四个季节SCSG效应像元分类结果

图5 NSSR几何效应校正示意图

图6几何效应校正后短波净辐射数据

图7机器学习回归模型Pearson特征相关矩阵热图

图8 2020年5月2日晴空和阴影区各输入特征变量与LST密度分布图

图9最大树深度与模型验证集得分

图10三种模型对不同的训练样本数量所需要消耗的时间

图11训练精度与所用的时间情况

图12 RF，GBRT和Xgboost算法云下地表温度回归运算特征重要性

图13基于三种集成回归模型的FY-3D MERSI-II云下地表温度影像图

图14基于FY-3D MERSI-II数据云下地表温度与CLDAS温度散点图

具体实施方式

为了实现所提出的方法，我们以FY-3D中分辨率成像仪MERSI-II为例，进行云下地表温度估算的实验。我们需要在目标图像和参考图像之间进行几何配准，并确定像元是晴空还是云覆盖。其过程共有三部分：参考图像的辐射定标和几何校正，晴空地表温度反演，云检测识别云遮挡、云覆盖像元、阴影像元和晴空像元，特征变量的选择和预处理，模型建立和训练，云下地表温度估算。

FY-3D我国第二代极轨气象卫星FY-3系列卫星的第四颗星，中分辨率光谱成像仪II(Medium Resolution Spectral Imager,MERSI-II)是搭载在FY-3D上的核心仪器之一，其整合了原有风云三号卫星成像仪器MERSI和VIRR的功能，共设有25个通道，其中通道1-19为太阳反射通道(0.4～2.1μm)，20-25为红外发射通道(3.8～12.5μm)，通道1-4、24-25的空间分辨率为250米，通道5-23的空间分辨率为1000m。FY-3D MERSI-II目前主要提供L1、L2和L3三级数据产品，共计25个数据集，数据存储格式为HDF5科学数据集格式。

本发明以太阳-云-卫星几何效应为出发点，发展FY-3D MERSI-II云下地表温度集成回归算法。首先从太阳辐射、地表和大气条件三方面对地表温度的影响出发，选取了7个变量作为地表热辐射关联因子，重点对关键特征变量短波净辐射进行了太阳-云-卫星的几何效应的校正，并对各个特征变量进行了共线性分析等预处理。采用随机森林(RandomForest,RF)、梯度回归树(Gradient Boosting Tree,GBRT)和自回归极端梯度上升算法(Extreme Gradient Boosting,XGBoost)三种集成回归模型，以FY-3D MERSI-II晴空和云阴影区域数据作为训练数据集，云遮挡与云覆盖区地表温度作为目标变量进行模拟预测，建立了一套充分考虑地表温度和辐射通量太阳-云-卫星几何效应的云下LST集成回归算法。基于该算法采用2020年华北区域不同季节共四天数据，进行云下地表温度预测。基于FY-3D MERSI-II开展晴空和云下地表温度估算，需要MERSI-II几何定位数据和部分数据产品(云高、云检测、NVI和LAI旬产品)作为辅助数据。结合附图对本发明作进一步详细的说明。表2显示了全天候地表温度估算所用的FY-3D MERSI-II数据集。

表2全天候地表温度估算所用的FY-3D MERSI-II数据集

步骤一、FY-3D MERSI-II数据预处理

根据实现选定的华北发明区范围，采集FY-3D MERSI-II的L1级卫星图像数据。该部分工作主要包括反射通道星上反射率计算、热发射通道亮度温度计算和几何校正三部分工作。

1.星上反射率计算

MERSI-II L1级数据采用无符号整型数据存储原始观测像元以减少数据量大小，灰度值经过多探元归一化订正后的输出。辐射定标是定量遥感应用的必要步骤，将传感器原始观测记录的无量纲DN值，转换为有实际物理意义的光谱反射率或辐亮度数据。

太阳反射波段CH1～19星上反射率计算，首先需要将DN值转换为可见光定标反射率Ref。数据集EV_1KM_RefSB和EV_250_Aggr.1KM_RefSB分别为1km反射通道(CH5～19)和250m反射通道(CH1～4)的归一化的DN值，其定标反射率转换公式如下：

dn＝DN*Slope+Intercept (1-17)

Ref＝Cal₂*dn²+Cal₁*dn+Cal₀ (1-18)

式中，Cal₀、Cal₁和Cal₂分别为数据集VIS_Cal_Coeff中相应通道的定标系数(分别对应第1，2，3列)，Slope和Intercept为数据集EV_1KM_RefSB和EV_250_Aggr.1KM_RefSB属性。

之后将可见光定标反射率Ref转换为入瞳辐亮度L_toa：

再结合入瞳辐亮度L_toa即可转换为表观反射率ρ_toa：

式中，E₀为通道大气层外太阳辐照度，可以从L1数据集中属性Solar_Irradiance，DES取自L1数据集中属性EarthSunDistanceRatio。μ为cosθ，θ为太阳天顶角，可以从MERSI-IIL1科学数据集SolarZenith属性中获取，由于存储方式是缩放后的整型数据，需要根据以下公式用相应的slope和Intercept系数，计算真实的太阳天顶角θ。

θ＝SolarZenith×Slope+Intercept (1-21)

2.亮度温度计算

FY-3D MERSI-II数据EV_1KM_Emissive和EV_250_Aggr.1KM_Emissive分别为1km发射通道(CH20～23)和250m发射通道(CH24～25)的放大后辐亮度值，根据属性中获取每个波段i增益(scales)和偏置(offsets)，根据以下公式计算波段i辐亮度I_i(W·m^-2·sr^-1·μm^-1)(徐娜等，2019)：

I_i＝(DN-offset)·scales (1-22)

热红外通道(CH24～25)辐亮度转换为黑体亮温，通过等效中心波长effect_Center_WaveLength(单位μm)以及通道辐亮度，用Planck函数逆变换计算得到等效黑体亮温T_ei，计算公式如下：

式中，λ_i是第i(i＝24，25)波段的有效中心波长，可以根据为L1文件里面的属性变量effect_Center_WaveLength(单位μm)获取；由表3-1可知，MERSI-II第24和25通道的中心波长λ_i分别取λ₂₄＝11.8μm和λ₂₅＝12.0μm；C₁和C₂分别为第1和2光谱常量，分别取C₁＝1.19104356×10^-16W/m²和C₂＝1.4387685×10⁴μm·K。最后需要利用通道亮温修正系数(TbbCorr_Coeff)，将T_ei转换为通道黑体亮温T_i：

T_i＝A*T_ei+B (1-24)

式中，T_i是MERSI-II第i(i＝24，25)波段的亮度温度，即热红外通道24和25对应的亮度温度T₂₄和T₂₅；A和B分别为L1文件里面的属性变量TBB_Trans_Coefficient_A和TBB_Trans_Coefficient_B.

3.几何校正

FY-3DMERSI-II数据采用多元探测器并扫技术，传感器成像几何特性、地形起伏、地表曲率及卫星抖动等因素导致图像几何畸变现象的产生。本发明采用地理查找表(GeographicLookupTable，GLT)方法进行***定位工作。具体而言，通过几何定位文件输入的对应数据的几何定位信息，生成地理位置查找表。通过查找表可以推算出原始数据中各像元实际地理位置。但是由于地形、地表和大气折射等诸多因素影响，这种初步的***校正结果仍存在2-10像元的几何定位误差，还不能满足应用需求。为此，本发明采用了基于Moravec角点检测算法的遥感影像自动匹配技术，对初步校正的数据进行几何精校正。该算法通过自动化控制点选取，有效减少了影像的局部几何畸变，提高了数据的几何定位精度，误差基本控制在3个像元内，几何校正精度基本能够满足需求。精确的几何定位为后续FY-3DMERSI-II数据的全天候地表温度估算工作的顺利开展，奠定了良好的几何基础。精确的几何定位为后续FY-3DMERSI-II数据的全天候地表温度估算工作的顺利开展，奠定了良好的几何基础。准备与数据相应区域和时间的NDVI、LAI、DEM、坡度、坡向、太阳高度角、卫星观测角度等数据，确保数据的投影和分辨率一致。

步骤二、FY-3D MERSI-II晴空地表温度反演

将普朗克函数线性展开，经过一系列推导提出两因素***窗算法(TFSWA)，使用两个亮度温度通道、相应的波段平均地表发射率和大气透过率的一系列线性组合来估算LST，在-50°～50°之间，MERSI-II TFSWA算法波段24常数项a₂₄＝-53.477，b₂₄＝0.3951；波段25常数项为a₂₅＝-57.087，b₂₅＝0.4292。

步骤三、云阴影、云遮挡与云覆盖区检测

MERSI-II云检测产品主要提供了MERSI-II卫星数据云和晴空分布信息，是用于云下地表温度估算研究的基础性辅助资料。本发明使用的是FY-3D MERSI-II云检测产品由多波段(3，4，6，7，19，20，21，24，25)特征阈值法综合分析得到的云检测结果。根据云下像元对应的地表太阳辐射是否受到云削弱作用影响，采用不同的方法估算云下地表温度，需要检测出云阴影、云遮挡与云覆盖区对应的具***置。因此结合FY-3DMERSI-II的角度和云高数据，利用步骤三中的方法针对华北区域开展2020年不同季节共4天数据进行了云遮挡、云覆盖和云阴影检测工作。图3为FY-3D MERSI-II云阴影、云遮挡、云覆盖和晴空像元统计结果，图4显示了2020年不同季节的四幅影像像元分类结果。

步骤四、短波净辐射反演和几何效应校正

短波净辐射是云下地表温度反演中重要的特征参量。本次采用参数化算法反演FY-3D MERSI-II短波净辐射，并对短波净辐射进行SCSG几何效应校正。MERSI-II有两组可见光波段，其中波段8～14的带宽较窄(20nm),主要用于海洋颜色、浮游生物等发明，1～7波段为陆地应用设计，MERSI-II星上反照率的反演主要使用了这七个反射太阳辐射(RSB)通道和1个水汽吸收通道(表3)。

表3用于NSSR反演的MERSI-II反射太阳辐射波段

短波净辐射反演首先将窄带反射率ρ_i转换为TOA短波段宽带反照率γ，然后用参数化方法得到经验系数来计算NSSR。基于MODTRAN模型可以建立窄波段反射率与宽波段反照率的线性转换公式：

γ＝b₀+b₁ρ₁+b₂ρ₂+b₃ρ₃+b₄ρ₄+b₅ρ₅+b₆ρ₆+b₇ρ₇+b₁₉ρ₁₉ (1-25)

式中，γ为星上短波宽波段反照率，ρ_i为MERSI-II波段i的星上窄波段反射率，b_i为MERSI-II波段i对应的拟合参数。反照率随方位角的变化小于天顶角的变化，因此，选择范围从0°～60°、60°～120°和120°～180°的三个粗略间隔，分别表示正向、侧向和反向的近似反射。

表4 MODTRAN5.2大气辐射传输模拟中主要参数设置

根据步骤四的参数化过程就可以估算FY-3D MERSI-II短波净辐射。

步骤五、特征变量选取与预处理

(1)特征变量选取

MERSI-II LAI叶面积指数旬产品和NVI植被指数旬产品数据主要用于相似像元查找和作为集成回归算法的特征变量，两个数据集分别是GLL等经纬度投影和Hammer投影，都是采用全球10°×10°分幅0.01°空间分辨率存储。

MERSI-II植被指数旬产品数据，包含NDVI、EVI、传感器和卫星角度信息以及质量控制信息等。MERSI-II叶面积指数旬产品主要包括叶面积指数和相应的数据质量控制信息。云顶高度指云层顶所处的高度，单位为hPa或者m，全部数据统一转换为几何高单位后用于计算云阴影位置及辐射数据的几何效应校正。将所有辅助数据统一为WGS84地理坐标系，0.01°分辨率等经纬度投影。

(2)变量预处理

a.短波净辐射几何效应校正

考虑到辐射数据SCSG效应的普适性，将用于云下地表温度估算的关键特征变量短波净辐射进行SCSG几何效应校正，具体方法参考步骤四。图5显示了短波辐射数据几何效应示意图，该图显示的是2020年5月2日观测的华北地区的MERSI-IINSSR，其中图5(a)是RGB真彩色影像图，m和n两直线用于发明NSSR校正前后插值剖面曲线。图5(b)为各类型像元分类结果：A区云阴影区，B区云覆盖区，C区云遮挡区，D区晴空区。为了更好地表征太阳-云-卫星几何结构对NSSR的影响，将SCSG几何效应校正后的NSSR减去原始NSSR数据，详见图5(c)。从差值图中，可以看到影像中的阴影区(A区)属于负值，云遮挡区(C区)为正值，说明如果忽略太阳-云-卫星的几何效应，NSSR在云遮档区(C区)会被低估，而在阴影区(A区)则会被高估。这种差异的具体大小随云的物理和光学特征而变化。在给出的例子中，这种差异可以达到约800W/m²，约占晴空NSSR的80％。图5(d)和(e)是m和n两直线处NSSR差值(NSSR_corrected-cld-NSSR_clear)剖面图，左侧虚线部分对应的是原始数据的云遮挡区，为低估区，右侧的虚线对应的是原始数据的阴影区，属于NSSR高估区。图6给出了四个季节不同时间(20200218、20200502、20200715和20201113)短波净辐射数据几何效应校正后的影像图，对应的云检测结果，即像元分类结果。

b.NVI、EVI和LAI数据的谐波优化处理

为了获取较高精度的模式模拟结果,利用HANTS算法对FY-3DMERSI-IINVI、EVI和LAI数据的缺失像元进行填补，将受云污染或其他因素引起的低质量数据进行平滑和滤波处理，生成连续完整的长时间序列高质量的训练数据。

c.特征变量共线性分析

为了选取高质量的特征变量，需要检查变量之间是否存在共线性问题。基于Seaborn统计分析可视化模块计算了以上9个特征变量之间的PEARSON特征相关矩阵热图(图7)，可以看出，SZA和NSSR相关性在0.7以上，存在共线性问题，因此去掉SZA变量。其他变量不存在明显的共线性问题。此外进行特征重要性分析，发现坡向得分最低，因此为了降低噪音影响，去掉坡向特征变量，因此确定了7个有效特征变量。2020年5月2日，有802468个晴空和阴影像元；这些像元广泛分布涵盖了不同的土地覆盖类型、植被覆盖程度、地形特征和太阳辐射条件等。基于2020年5月2日的观测数据，分析了晴空和云阴影区域LST即遥感可视区域各特征变量与温度之间的关系。图8显示了2020年5月2日，晴空和阴影区确定的7个特征变量与LST的密度分布图。由图可知，地表温度和各个参数的分布情况较为集中,云下地表温度与参数之间没有明显的线性关系。这种复杂关系证实了地表温度受许多地表变量的共同影响，各因素之间的耦合效应导致了这种复杂性。

步骤六、模型建立与参数优化

(1)训练数据集的确定

为了准确地表示LST与特征变量之间的复杂关系，构建各特征变量到地表温度的拟合回归模型，选取2020年5月2日高质量的晴空和阴影数据为训练数据集,将云像元地表温度作为目标变量，建立回归模型。基于Python的第三方模块scikit-learn和XGBoost库，利用RF、GBRT和XGBoost三种回归模型对晴空和阴影像元的地表温度数据集进行学习训练，在此过程中利用train_test_split函数，将样本数据随机抽取划分，80％训练分区和20％测试分区两个部分，分别应用于训练模型和测试模型。

(2)参数设置与优化

为了防止过拟合现象，调整三个模型中模型回归树的个数(max-depth)参数，通过10折交叉验证查看训练数据的得分，结果如图9，随着max-depth从0递增到20，模型的验证得分递增，RF、GBRT和XGBoost分别是max-depth＝8，5，4时，模型得分达到极值。选择max_depth＝8，5，4为RF、GBRT和XGBoost的最优参数。为了查看模型性能，针对不同的训练样本查看其训练时间(图10)，对相同训练样本数量GBRT所需要消耗的时间最短，XGBoost时间最长。图11显示了达到不同精度各模型所用时间，由图分析可知，三种集成树类方法随机森林和XGBoosting算法消耗的时间相对较长，但达到相对较高的精度。

步骤七、云下地表温度回归运算

针对亚洲中东部区域，在2020年四个季节选择若干影像进行云下地表温度重建实验，评估所提出的方法。基于Python第三方扩展包scikit-learn特征变量重要性计算方法(Permutation Importance)，计算各特征变量的重要性得分和排序。图11给出了三种回归模型特征重要性计算结果。由图可知，三种模型所使用的参数中，排名前三的特征变量一致，对云下地表温度影响最大的参数是DEM，特征重要性从XGBoost的0.87到RF的0.95，其次为NWDI(0.76～0.82)和NSSR(0.35～0.39)，最后是植被指数和叶面积指数。DEM不同，太阳辐射的光学路径长度的会发生变化，此外太阳光线到达地面时穿过的大气层厚度，海拔越高，大气层越稀薄，辐射越强。高程条件决定了能量分配，因此在局部范围内对地表热辐射影响较大。短波净辐射是太阳辐射、云和地表的综合表征，是地表温度变化的主导因素之一。总体来说在发明区内，对云下地表温度影响最大的条件是地形条件，其次是NDWI和地表短波净辐射，最后是植被覆盖条件。图12基于三种机器学习回归模型算法估算的FY-3DMERSI-II云下地表温度分布影像图。

利用中国气象局陆面数据同化***(CMA Land Data Assimilation System,CLDAS)地表温度数据对该算法进行精度验证。图14显示的是基于FY-3D MERSI-II数据云下地表温度与CLDAS温度散点图。与CLDAS地表温度数据对比分析，结果显示该方法RMSE精度约为2.84K-5.43K。结果表明，考虑太阳-云-卫星观测几何效应的云下地表温度集成回归算法估算的地表温度精度较好，可以较好地刻画地形、植被对地表温度的影响。

Claims (10)

1.一种基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，包括以下步骤：

S1：采集FY-3D MERSI-II的L1级卫星图像数据并对其进行预处理；

S2：晴空地表温度反演；

S3：云阴影、云遮挡与云覆盖区检测；

S4：采用参数化算法反演FY-3D MERSI-II短波净辐射，并对短波净辐射进行SCSG几何效应校正；

S5：选取云下地表温度回归运算的特征变量并对其进行预处理；

S6：针对三种集成回归模型进行参数设置与优化；

S7：针对晴空和阴影区的LST，建立各特征变量到云下地表温度的集成回归模型；

S8：输入待计算图像，采用建立的回归模型对云遮挡与云覆盖区进行LST估算。

2.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S1具体包括以下步骤：

S11：计算星上反射率；

S12：计算亮度温度；

S13：基于采用地理查找表(GLT)方法进行***定位工作，并用Moravec角点检测算法对初步校正的数据进行几何精校正。

3.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S2具体包括以下步骤：

使用两个亮度温度通道、相应的波段平均地表发射率和大气透过率的一系列方程来估算LST。

4.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S3具体包括以下步骤：

确定云的正射位置：

然后估计对应云阴影的正射投影位置：

式中，(X_cld，Y_cld)为云在地表的正射投影；(X_shw，Y_shw)为云影在地表上的投影位置；(x_map,y_map)是遥感影像中的云位置即云团沿卫星观测方向的投影位置；而H是地面以上的云顶高度，由MERSI-II的云顶高度减去地面高度计算；θ_v和

是卫星过境时刻传感器观测天顶角和观测方位角，θ_s和

分别是传感器观测时太阳的天顶角和方位角。

5.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S4具体包括以下步骤：

短波净辐射反演首先将窄带反射率ρ_i转换为TOA短波段宽带反照率γ，然后用参数化方法得到经验系数来计算NSSR；基于MODTRAN模型可以建立窄波段反射率与宽波段反照率的线性转换公式：

γ＝b₀+b₁ρ₁+b₂ρ₂+b₃ρ₃+b₄ρ₄+b₅ρ₅+b₆ρ₆+b₇ρ₇+b₁₉ρ₁₉ (1-11)

式中，γ为星上短波宽波段反照率，ρ_i为MERSI-II波段i的星上窄波段反射率，b_i为MERSI-II波段i对应的拟合参数；

然后，逐像元确定云及其阴影范围，标定每个云像元与其阴影投影的位置关系，采用Hants算法插值计算生成云遮挡区像元NSSR，云覆盖区的NSSR进行几何位置校正，最终生成经过几何效应校正的短波辐射数据。

6.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其中，云下地表温度回归运算的特征变量具体可以是以下一个或多个变量的组合：叶面积指数、归一化植被指数、增强型植被指数、数字地面高程、坡度、坡向、归一化水体指数、太阳高度角、短波净辐射。

7.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S6具体包括：

采用RF、GBRT和XGBoost三种回归模型对晴空和阴影像元的温度数据集进行学习训练，利用train_test_split函数，将样本数据随机抽取划分为训练分区和测试分区两个部分，分别应用于训练模型和测试模型。

8.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S7具体包括：

采用RF、GBRT和XGBoost三种集成回归算法，以FY-3D MERSI-II云下地表温度作为目标变量进行回归运算，并利用CLDAS同化LST数据对估算的云下地表温度进行精度评价。

9.根据权利要求1所述的基于太阳-云-卫星观测几何的云下地表温度估算方法，其特征在于，所述步骤S8具体包括：

采用构建的回归模型对发明区内云下地表温度进行预测，将预测的云像元LST与晴空LST合并，得到发明区的无缝全天候LST影像。

10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1-9中任一项所述的方法。

Images