一种遥感地表温度数据的时空一体化融合方法
【技术领域】
[0001] 本发明涉及定量遥感数据处理技术领域，尤其是涉及一种遥感地表温度数据的时 空一体化融合方法。
【背景技术】
[0002] 地表温度（LST，Land Surface Temperature)是研究地-气间物质与能量交换的 重要参数，广泛应用于城市热岛效应、土壤湿度估算、森林火灾监测、气候变化等领域。卫星 热红外遥感是获取地表温度的重要手段，然而，遥感观测系统需要在空间分辨率和时间分 辨率之间进行权衡与取舍，难以同时获得高时间和高空间分辨率的温度数据，极大限制了 其应用的广度与深度。例如：静止卫星地表温度数据的时间分辨率最高达到15分钟，却具 有至少3千米的低空间分辨率；而极地卫星提供的地表温度数据虽然具有较高的空间分辨 率（约100米-1000米），但其时间分辨率一般都较低。另一方面，现有发布的地表温度数 据由于云、气溶胶等不利天气的影响往往出现空值现象，导致空间不连续。因此，如何在不 改变现有观测系统的条件下，实现遥感地表温度数据空间无缝的同时提高遥感LST的空间 分辨率和时间分辨率，成为热红外遥感应用中的关键问题之一。
[0003] 为了提高遥感LST数据的时空分辨率，国内外学者已经开展了一些研究。早期的 研究主要有遥感LST空间降尺度/时间扩展：在空间降尺度方面，李小文院士提出将传统地 学趋势面用于遥感尺度转换研究，其主要原理都是引入更高空间分辨率的地理要素构造趋 势面，通过建立趋势面与LST间的关系完成降尺度转换。而对于遥感LST时间尺度扩展，一 般是对多时相数据进行时相滤波、谐波分析等方法重建出中间时间的数据。以上针对单一 遥感参量来提高时空分辨率的方法，对辅助数据质量或长时序性要求较高，且提高的分辨 率程度有限，在某些应用也显得无能为力。因此，需发展一种新的LST时空降尺度方法。
[0004] 多传感器遥感LST时空融合技术正是为了满足这种需求而发展起来的一种新技 术。与上述针对单一遥感参量的空间降尺度、时间扩展方法不同的是，遥感LST时空融合是 通过对具有不同空间尺度和时间尺度的LST进行融合处理，生成时间相对密集的高空间分 辨率LST数据；其优势在于无需其它辅助数据条件，同时生成的数据具有输入数据中最高 空间分辨率和最高时间分辨率。
[0005] 然而，无论是传统时空降尺度还是时空融合，在算法实施过程中都假设所输入数 据空间无缝，这个假设对于地表温度数据（由于受云等影响）经常不能满足；此外，在基 于多时相、多尺度温度数据时空融合时，很少顾及参量昼夜变化规律、遥感观测条件等的差 另IJ，直接影响了融合影像的求解精度。因此如何利用多时相、多尺度时空互补信息实现遥感 地表温度数据空间无缝同时提高其时空分辨率，需要研究新的时空一体化融合方法。
【发明内容】
[0006] 本发明所要解决的问题是：针对现有观测系统的技术限制和遥感地表温度数据时 空不连续现象，提供一种遥感地表温度数据的时空一体化融合方法。
[0007] 为解决上述技术问题，本发明提出的技术方案为：一种遥感地表温度数据的时空 一体化融合方法，设有一个空间无缝的极轨卫星地表温度数据M和同一地区的多时相静止 卫星地表温度数据Gtp(p = 0, 1，"·η且n>10);进一步假设在时间序列tQ、1v"tn Mn上有 一景地表温度数据G与极轨卫星地表温度数据M获取时间上对应，则实现时空一体化融合 包括以下步骤：
[0008] 步骤一、对地表温度数据M和Gtp(p = 0, 1，…η)进行预处理，包括几何配准、裁剪、 温度单位转换、空值标记、重采样；
[0009] 步骤二、基于昼夜变化模型的空间修复：根据当前时相待修复像元在其它时相对 应有效像元值的数量分如下两种情况求解所述昼夜变化模型的参数，并进行空值修复； [0010] 情况1)若当前时相的待修复像元在其它时相上有足够非空值像元，直接利用这 些有效值与对应时间构成的数据对按需求分白天和晚上进行所述昼夜变化模型参数拟 合；
[0011] 情况2)若待修复像元的其它时相有效像元及其对应时间构成的数据对不足以 进行温度昼夜变化模型参数拟合时，则以相邻时相1对应的非空值像元为中心像元开一 个窗口，设窗口大小为waX wa，wa根据G的特征而定，在窗口内筛选相似像元，筛选条件 为：邻域像元与中心像元差的绝对值小于或等于该中心像元相关的函数，用公式表达为：
，其中G(i，tx)为时间^时静止卫星地表温度数 据G上窗口内邻域像元，G(wa/2, tx)为静止卫星地表温度数据G上窗口内中心像元，（13为 常量，使得筛选的相似像元满足其它时相上有足够多的有效数据对；利用这些有效数据对 进行所述昼夜变化模型参数的拟合数据数据。
[0012] 步骤三、构建遥感地表温度数据时空一体化融合模型：
[0013] 在地表温度数据M上逐个以像元为中心像元开一个大小为wbXwb的窗口，窗口大 小视M特征而定，按以下公式在窗口内筛选相似像元：
[0015] 其中M(i, t)为时间t时M上窗口内邻域像元，M(wb/2, t)为M上窗口内中心像元， db为常量，
[0016] 计算相似像元地表温度数据在当前时间和预测时间间的温度变化值与中心像元 所在昼夜变化模型的温度变化值的差异TD1如下，
[0017] TDi= IG (i, t) -G (i, t x) - (Temp (wb/2, t) -Temp (wb/2, tx))
[0018] 其中，G (i，t)为当前时间t的静止卫星地表温度数据G上窗口内第i个邻域像元， G (i，tx)为时间tx时静止卫星地表温度数据G上窗口内第i个邻域像元，i为正整数，其最 大值为窗口大小；
[0019] 结合TD1及其它参量构建以下相似像元的权重函数表达式：
[0020] W1= (!/(V1 · SD1)/ Σ !/(V1 · SD1))
[0021] 其中SD1是窗口内相似像元与中心像元的相似度、V1是由M上窗口内相似像元与 中心像元温度变化差异、欧式距离以及不同分辨率地表温度数据间尺度差异共同构成；假 设地表温度数据M和Gtn在时间t n时对应，作为基础时间对，其它时间序列t xe {t。、心、… tn J作为预测时间，得到时空一体化融合模型： CN 105184076 A 说明书 3/6 页
[0023] 其中，M(wb/2, tx)为融合的空间无缝且高时空分辨率的地表温度数据，具有静止 卫星G的时间分辨率和极轨卫星M的空间分辨率。
[0024] 进一步的，在所述步骤二中，由于考虑到观测条件、传感器误差、反演误差的原因 使得现有地表温度数据存在包含一些噪声、误差较大的离群点，采用M估计理论提高模型 参数拟合精度。
[0025] 所述步骤二中，当待修复像元的其它时相像元都为空值时，选择当前时相邻域有 效像元值，包括原有的和建模求解的；利用空间内插技术完成修复。
[0026] 本发明充分结合多尺度、多时相传感器间的时空互补信息和静止卫星地表温度的 昼夜变化规律，利用静止卫星高时间分辨率特性建立地表温度昼夜变化模型，进而实现空 值像元的修复，同时与较高分辨率极轨卫星温度数据构建时空融合模型，最终生成空间无 缝的高时空分辨率地表温度数据集。而且不需要任何辅助数据，易于实现，实用价值高，可 用于业务化运行。本发明与现有技术相比具有空间修复精确、时空权重函数合理、时空一体 化融合等优点。主要体现在：首先，根据待修复空值像元在时间维上的有效像元数分情况求 解地表温度昼夜变化模型参数，并利用M估计等稳健回归方法使模型参数最优化，从而提 高修复精度；其次，通过引入相似像元本身的相似度差异因子、基于温度昼夜变化模型的时 间差异因子等构建更加合理的时空权重函数；最后，利用多时相、多尺度地表温度数据进行 时空一体化融合可实现空间无缝的同时提高其时空分辨率。总之，新方法可以批量生成空 间上连续且时间上相对密集的地表温度数据集，特别适用于具有高时间分辨率的静止卫星 温度数据。
【附图说明】
[0027] 图1是本发明的实施例流程图。
【具体实施方式】
[0028] 以下结合附图和实施例详细说明本发明技术方案。
[0029] 假设有一个空间无缝的高空间分辨率极轨卫星地表温度数据M和同一地区的多 时相静止卫星地表温度数据Gtp(p = 0, 1，…n，且n>10)。进一步假设在时间序列t。、tr·· tn i、tn上有一景地表温度数据G与极轨卫星地表温度数据M获取时间上对应，通过时空一 体化融合就能获取空间无缝且具有极轨卫星地表温度数据空间分辨率的Mtp(p = 0, 1，… η) 〇
[0030] 如图1所示，本发明的实施例流程可以分为三个阶段：
[0031] 第一阶段：地表温度数据的预处理
[0032] 对地表温度数据M和Gtp (ρ = 0, 1，…η)进行预处理，一般包括几何配准、裁剪、温 度单位转换、空值标记、重采样。几何配准是通过地面控制点及相应的经炜度信息完成投影 转换与定位，根据所感兴趣的研究区域范围分别对所有温度数据进行裁剪，将所有温度数 据转换成相同的单位，对地表温度数据Gtp空值像元进行标记，在完成空间修复后将Gt ρ (ρ =0, 1，…η)重采样与地表温度数据M相匹配的空间分辨率。
[0033] 第二阶段、空间无缝处理
[0034] 步骤2-1、根据预处理中空值像元标记结果，若所有时相Gtp均无空值，直接进入第 三阶段；否则，进入步骤2-2;
[0035] 步骤2-2、确定地表温度昼夜变化模型参数求解的数据对（温度与时间）。根据当 前时相待修复像元在其它时相对应有效像元值的数量分如下三种情况确定数据对：
[0036] 情况一、若当前时相的待修复空值像元在其它时相上有足够数量（大于6)非空值 像元，可确定这些有效温度Temp (t)及其观测时间t作为如下地表温度昼夜变化模型参数 求解的数据对：
[0038] 其中，bJP b 2表示如下：
[0039] bi = T 〇+Tacos ( β (ts-tm)) -b2 (2)
[0040] b2= [_β Tasin(β (ts_tm))]/a (3)
[0041] 式中个参数定义如表1所示
[0042] 表1 :各参数的意义
[0044] 情况二、当待修复像元的其它时相有效像元数及其观测时间所组成的数据对不足 以求解温度昼夜变化模型参数时，则以相邻时相1对应的非空值像元为中心像元开一个窗 口，设窗口大小为wa X wa，wa根据G的特征而定，在窗口内筛选相似像元，筛选条