一种遥感地表温度数据的时空一体化融合方法

摘要

本发明涉及定量遥感数据处理技术领域，公开了一种遥感地表温度数据的时空一体化融合方法，首先对遥感地表温度数据进行预处理，其次，利用多时相温度数据和对应的时间求解地表温度昼夜变化模型参数，实现空间无缝，最后在极轨卫星地表温度影像数据上逐窗口筛选相似像元，结合由多时相静止卫星与上述模型构成的时间差异、尺度差异、相似度差异以及欧氏距离等权重因子，共同构建时空一体化融合模型并求解生成空间无缝的高时空分辨率地表温度数据集。本发明充分利用极轨卫星与多时相静止卫星地表温度数据的时空互补信息和静止卫星地表温度的昼夜变化规律，不需要任何辅助数据，易于实现，实用价值高，可用于时空无缝的遥感地表温度数据业务化运行。

说明书

技术领域

本发明涉及定量遥感数据处理技术领域，尤其是涉及一种遥感地表温度数据的时空一体化融合方法。

背景技术

地表温度(LST，LandSurfaceTemperature)是研究地-气间物质与能量交换的重要参数，广泛应用于城市热岛效应、土壤湿度估算、森林火灾监测、气候变化等领域。卫星热红外遥感是获取地表温度的重要手段，然而，遥感观测系统需要在空间分辨率和时间分辨率之间进行权衡与取舍，难以同时获得高时间和高空间分辨率的温度数据，极大限制了其应用的广度与深度。例如：静止卫星地表温度数据的时间分辨率最高达到15分钟，却具有至少3千米的低空间分辨率；而极地卫星提供的地表温度数据虽然具有较高的空间分辨率(约100米-1000米)，但其时间分辨率一般都较低。另一方面，现有发布的地表温度数据由于云、气溶胶等不利天气的影响往往出现空值现象，导致空间不连续。因此，如何在不改变现有观测系统的条件下，实现遥感地表温度数据空间无缝的同时提高遥感LST的空间分辨率和时间分辨率，成为热红外遥感应用中的关键问题之一。

为了提高遥感LST数据的时空分辨率，国内外学者已经开展了一些研究。早期的研究主要有遥感LST空间降尺度/时间扩展：在空间降尺度方面，李小文院士提出将传统地学趋势面用于遥感尺度转换研究，其主要原理都是引入更高空间分辨率的地理要素构造趋势面，通过建立趋势面与LST间的关系完成降尺度转换。而对于遥感LST时间尺度扩展，一般是对多时相数据进行时相滤波、谐波分析等方法重建出中间时间的数据。以上针对单一遥感参量来提高时空分辨率的方法，对辅助数据质量或长时序性要求较高，且提高的分辨率程度有限，在某些应用也显得无能为力。因此，需发展一种新的LST时空降尺度方法。

多传感器遥感LST时空融合技术正是为了满足这种需求而发展起来的一种新技术。与上述针对单一遥感参量的空间降尺度、时间扩展方法不同的是，遥感LST时空融合是通过对具有不同空间尺度和时间尺度的LST进行融合处理，生成时间相对密集的高空间分辨率LST数据；其优势在于无需其它辅助数据条件，同时生成的数据具有输入数据中最高空间分辨率和最高时间分辨率。

然而，无论是传统时空降尺度还是时空融合，在算法实施过程中都假设所输入数据空间无缝，这个假设对于地表温度数据(由于受云等影响)经常不能满足；此外，在基于多时相、多尺度温度数据时空融合时，很少顾及参量昼夜变化规律、遥感观测条件等的差别，直接影响了融合影像的求解精度。因此如何利用多时相、多尺度时空互补信息实现遥感地表温度数据空间无缝同时提高其时空分辨率，需要研究新的时空一体化融合方法。

发明内容

本发明所要解决的问题是：针对现有观测系统的技术限制和遥感地表温度数据时空不连续现象，提供一种遥感地表温度数据的时空一体化融合方法。

为解决上述技术问题，本发明提出的技术方案为：一种遥感地表温度数据的时空一体化融合方法，设有一个空间无缝的极轨卫星地表温度数据M和同一地区的多时相静止卫星地表温度数据Gt_p(p＝0,1,…n且n>10)；进一步假设在时间序列t₀、t₁…t_n-1、t_n上有一景地表温度数据G与极轨卫星地表温度数据M获取时间上对应，则实现时空一体化融合包括以下步骤:

步骤一、对地表温度数据M和Gt_p(p＝0,1,…n)进行预处理，包括几何配准、裁剪、温度单位转换、空值标记、重采样；

步骤二、基于昼夜变化模型的空间修复：根据当前时相待修复像元在其它时相对应有效像元值的数量分如下两种情况求解所述昼夜变化模型的参数，并进行空值修复；

情况1)若当前时相的待修复像元在其它时相上有足够非空值像元，直接利用这些有效值与对应时间构成的数据对按需求分白天和晚上进行所述昼夜变化模型参数拟合；

情况2)若待修复像元的其它时相有效像元及其对应时间构成的数据对不足以进行温度昼夜变化模型参数拟合时，则以相邻时相t_x对应的非空值像元为中心像元开一个窗口，设窗口大小为wa×wa，wa根据G的特征而定，在窗口内筛选相似像元，筛选条件为：邻域像元与中心像元差的绝对值小于或等于该中心像元相关的函数，用公式表达为：其中G(i,t_x)为时间t_x时静止卫星地表温度数据G上窗口内邻域像元，G(wa/2,t_x)为静止卫星地表温度数据G上窗口内中心像元，d_a为常量，使得筛选的相似像元满足其它时相上有足够多的有效数据对；利用这些有效数据对进行所述昼夜变化模型参数的拟合数据数据。

步骤三、构建遥感地表温度数据时空一体化融合模型：

在地表温度数据M上逐个以像元为中心像元开一个大小为wb×wb的窗口，窗口大小视M特征而定，按以下公式在窗口内筛选相似像元:

|M(i,t)-M(wb/2,t)|≤d_b×2^M(wb/2,t)

其中M(i,t)为时间t时M上窗口内邻域像元，M(wb/2,t)为M上窗口内中心像元，d_b为常量，

计算相似像元地表温度数据在当前时间和预测时间间的温度变化值与中心像元所在昼夜变化模型的温度变化值的差异TD_i如下，

TD_i＝|G(i,t)-G(i,t_x)-(Temp(wb/2,t)-Temp(wb/2,t_x))|

其中，G(i,t)为当前时间t的静止卫星地表温度数据G上窗口内第i个邻域像元，G(i,t_x)为时间t_x时静止卫星地表温度数据G上窗口内第i个邻域像元，i为正整数，其最大值为窗口大小；

结合TD_i及其它参量构建以下相似像元的权重函数表达式：

W_i＝(1/(V_i·SD_i)/∑1/(V_i·SD_i))

其中SD_i是窗口内相似像元与中心像元的相似度、V_i是由M上窗口内相似像元与中心像元温度变化差异、欧式距离以及不同分辨率地表温度数据间尺度差异共同构成；假设地表温度数据M和Gt_n在时间t_n时对应，作为基础时间对，其它时间序列t_x∈{t₀、t₁、…t_n-1}作为预测时间，得到时空一体化融合模型：

$>M(wb/2,t_x)=\underset{i=1}{\overset{w}{\Sigma }}{W}_i\times \left(M\right(i,t_n\left)+G\right(i,t_n)-G(i,t_x\left)\right)$>

其中，M(wb/2,t_x)为融合的空间无缝且高时空分辨率的地表温度数据，具有静止卫星G的时间分辨率和极轨卫星M的空间分辨率。

进一步的，在所述步骤二中，由于考虑到观测条件、传感器误差、反演误差的原因使得现有地表温度数据存在包含一些噪声、误差较大的离群点，采用M估计理论提高模型参数拟合精度。

所述步骤二中，当待修复像元的其它时相像元都为空值时，选择当前时相邻域有效像元值，包括原有的和建模求解的；利用空间内插技术完成修复。

本发明充分结合多尺度、多时相传感器间的时空互补信息和静止卫星地表温度的昼夜变化规律，利用静止卫星高时间分辨率特性建立地表温度昼夜变化模型，进而实现空值像元的修复，同时与较高分辨率极轨卫星温度数据构建时空融合模型，最终生成空间无缝的高时空分辨率地表温度数据集。而且不需要任何辅助数据，易于实现，实用价值高，可用于业务化运行。本发明与现有技术相比具有空间修复精确、时空权重函数合理、时空一体化融合等优点。主要体现在：首先，根据待修复空值像元在时间维上的有效像元数分情况求解地表温度昼夜变化模型参数，并利用M估计等稳健回归方法使模型参数最优化，从而提高修复精度；其次，通过引入相似像元本身的相似度差异因子、基于温度昼夜变化模型的时间差异因子等构建更加合理的时空权重函数；最后，利用多时相、多尺度地表温度数据进行时空一体化融合可实现空间无缝的同时提高其时空分辨率。总之，新方法可以批量生成空间上连续且时间上相对密集的地表温度数据集，特别适用于具有高时间分辨率的静止卫星温度数据。

附图说明

图1是本发明的实施例流程图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

假设有一个空间无缝的高空间分辨率极轨卫星地表温度数据M和同一地区的多时相静止卫星地表温度数据Gt_p(p＝0,1,…n,且n>10)。进一步假设在时间序列t₀、t₁…t_n-1、t_n上有一景地表温度数据G与极轨卫星地表温度数据M获取时间上对应，通过时空一体化融合就能获取空间无缝且具有极轨卫星地表温度数据空间分辨率的Mt_p(p＝0,1,…n)。

如图1所示，本发明的实施例流程可以分为三个阶段：

第一阶段：地表温度数据的预处理

对地表温度数据M和Gt_p(p＝0,1,…n)进行预处理，一般包括几何配准、裁剪、温度单位转换、空值标记、重采样。几何配准是通过地面控制点及相应的经纬度信息完成投影转换与定位，根据所感兴趣的研究区域范围分别对所有温度数据进行裁剪，将所有温度数据转换成相同的单位，对地表温度数据Gt_p空值像元进行标记，在完成空间修复后将Gt_p(p＝0,1,…n)重采样与地表温度数据M相匹配的空间分辨率。

第二阶段、空间无缝处理

步骤2-1、根据预处理中空值像元标记结果，若所有时相Gt_p均无空值，直接进入第三阶段；否则，进入步骤2-2；

步骤2-2、确定地表温度昼夜变化模型参数求解的数据对(温度与时间)。根据当前时相待修复像元在其它时相对应有效像元值的数量分如下三种情况确定数据对：

情况一、若当前时相的待修复空值像元在其它时相上有足够数量(大于6)非空值像元，可确定这些有效温度Temp(t)及其观测时间t作为如下地表温度昼夜变化模型参数求解的数据对：

$>Temp\left(t\right)=\left(\begin{array}{c}{T}_0+T_a\cos \left(\beta \right(t-t_m\left)\right),(t<t_{s}day)\\ b_1+b_2{e}^{\alpha (t-t_s)},(t<t_s,night)\end{array}\right)---\left(1\right)$>

其中，b₁和b₂表示如下：

b₁＝T₀+T_acos(β(t_s-t_m))-b₂(2)

b₂＝[-βT_asin(β(t_s-t_m))]/α(3)

式中个参数定义如表1所示

表1：各参数的意义

情况二、当待修复像元的其它时相有效像元数及其观测时间所组成的数据对不足以求解温度昼夜变化模型参数时，则以相邻时相t_x对应的非空值像元为中心像元开一个窗口，设窗口大小为wa×wa，wa根据G的特征而定，在窗口内筛选相似像元，筛选条件为：邻域像元与中心像元差的绝对值小于或等于该中心像元相关的函数，用公式表达为：

$>|G(i,t_x)-G(wa/2,t_x)|\le d_a\times 2^{G(wa/2,t_x)}$>

其中G(i,t_x)为时间t_x时静止卫星地表温度数据G上窗口内第i个邻域像元，i为正整数，其最大值为窗口大小，G(wa/2,t_x)为静止卫星地表温度数据G上窗口内中心像元，d_a为常量，使得筛选的相似像元满足其它时相上有足够多的有效像元。利用这些有效像元值及其观测时间来求温度昼夜变化模型参数。

情况三、当待修复像元的其它时相像元都为空值时，暂时不予确定数据对。

步骤2-3、温度昼夜变化模型参数求解。依次针对每一个能求其模型参数的空值像元，将静止卫星观测时间与对应的地表温度值带入公式(1)进行参数拟合求解。由于观测条件、传感器误差、反演误差等原因使得现有地表温度数据可能包含一些噪声、误差较大的离群点，为了提高参数求解的精度，必须加以剔除。为此，引入了基于M-估计理论的稳健拟合方法，计算出最优的参数。

令ε_i,t为(i,t)像素点的地表温度差异，即：

$>{ϵ}_{i,t}=\left(\begin{array}{c}{T}_0+T_a\cos \left(\beta \right(t-t_m\left)\right)-Temp\left(t\right),t<t_s\\ b_1+b_2{e}^{\alpha (t-t_s)}-Temp\left(t\right),t\ge t_s\end{array}\right).$>

使其满足所有像元的地表温度差异和最小：其中ρ(·)有多种代价函数，本发明中使用经典的Huber代价函数：

$>\rho \left({ϵ}_{i,t}\right)=\{\begin{array}{cc}{{ϵ}_{i,t}}^2/2& \left|{ϵ}_{i,t}\right|\le k\\ k\left|{ϵ}_{i,t}\right|-k^2/2& \left|{ϵ}_{i,t}\right|>k\end{array},$>

其中，k为Huber参数。

为了方便求解，将尺度差异和最小化问题转化为迭代重加权最小二乘问题，即

$>\min \underset{i,t}{\Sigma }{\left(\omega \right({ϵ}_{i,t})*{ϵ}_{i,t})}^2,$>

这里

$>\omega \left({ϵ}_{i,t}\right)=\frac{{\rho }^{\prime }\left({ϵ}_{i,t}\right)}{{ϵ}_{i,t}}$>

表示权重,ρ′(ε_i,t)是Huber函数ρ(ε_i,t)的导数，通过这种方式依次可以求出所有能建模空值像元的模型参数值，用于后续像元填充。

步骤2-4、空值像元填充。当空值像元所在的模型参数获得后，将上述空值像元观测时间带入模型，即可求得改时间所对应的温度值，完成像元填充。

在情况三下，由于没有建模，部分像元依然会存在空值，此时可选择当前时相邻域有效像元值(包括原有的和建模求解的)，利用空间内插技术完成修复。由于静止卫星地表温度数据数量能够保证(目前发布的静止卫星一天最少24景，最多96景)，这种特殊情况一般很少见。

第三阶段、时空一体化融合处理

步骤3-1、在M上以逐个以像元为中心像元开一个大小为wb×wb的窗口，wb根据M的特征而定，在窗口内筛选相似像元，筛选条件为：邻域像元与中心像元差的绝对值小于或等于该中心像元相关的函数，用公式表达为：

|M(i,t)-M(wb/2,t)|≤d_b×2^M(wb/2,t)

其中M(i,t)为时间t时地表温度数据M上窗口内邻域像元，M(wb/2,t)为相应窗口内中心像元，d_b为常量，对不同的传感器会有稍微的变化；窗口内满足筛选条件的邻域像元即为相似像元。

步骤3-2、实施例根据M和Gt_p之间的尺度差异、相似像元相似度、不同时相Gt_p之间的时间差异以及欧氏距离，构建相似像元的权重函数表达式。对当前窗口内相似像元赋权，首先计算窗口内相似像元到中心像元的距离D_i，表达式为：

$>D_i=1+\sqrt{{(x_{w_b/2}-x_i)}^2+{(y_{w_b/2}-y_i)}^2}/(w_b/2)$>

这样的表达式是为了控制距离因子在1到1+2^0.5之间。

其次，计算窗口内各相似像元M(i,t)与中心像元M(wb/2,t)的相似度SD_i，在计算之前须先求出它们的相似性差异S_i如下，

S_i＝|M(i,t)-M(wb/2,t)|

由于相似性差异越小，表明该相似像元具有与中心像元更高的相似度，计算公式如下，

SD_i＝exp(-S_i)/∑exp(-S_i)

再次，计算当前时间不同分辨率地表温度数据M(i,t)与G(i,t)间的尺度差异R_i如下，

R_i＝|M(i,t)-G(i,t)|

最后，计算相似像元地表温度数据在当前时间和预测时间间的温度变化值与中心像元所在模型的温度变化值的差异TD_i如下，

TD_i＝|G(i,t)-G(i,t_x)-(Temp(wb/2,t)-Temp(wb/2,t_x))|

其中，G(i,t)为当前时间t的静止卫星地表温度数据G上窗口内第i个邻域像元；

合并上述尺度因子、时间因子和距离因子，得到E_i如下，

E_i＝ln(R_i*100+1)*ln(TD_i*100+1)*D_i

为了与相似度因子匹配，需要对E_i进行归一化，得到

V_i＝E_i/∑E_i(14)

为保证权重总和为1,建立如下归一化的权重函数W_it

W_i＝(1/(V_i·SD_i)/∑1/(V_i·SD_i))

该式表示尺度差异越小、时间差异越小、相似度越高且离中心像元距离越近的邻域像元获得的权重值越大。

步骤3-3、计算融合后的地表温度值。设地表温度数据M和Gt_n在时间t_n时对应，可作为基础时间对，其它时间序列t_x∈{t₀、t₁、…t_n-1}作为预测时间，由地表温度数据M和Gt_n和窗口内相似像元及其权重函数加权计算窗口内中心像元的温度值，

$>M(wb/2,t_x)=\underset{i=1}{\overset{w}{\Sigma }}{W}_i\times \left(M\right(i,t_n\left)+G\right(i,t_n)-G(i,t_x\left)\right)$>

其中，M(wb/2,t_x)为融合的空间无缝且高时空分辨率的地表温度数据，具有静止卫星G的时间分辨率和极轨卫星M的空间分辨率；

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。