一种异质下垫面地表温度全天时采样优化布设方法与流程

本发明涉及一种采样优化布设方法，尤其涉及一种异质下垫面地表温度全天时采样优化布设方法，属于遥感定标与真实性检验技术领域。

背景技术：

地表温度在区域及全球尺度的地表能量平衡、水循环过程、气候变化中扮演着重要的角色。近几十年来，基于卫星热红外数据反演地表温度的研究一直是定量遥感研究领域的重点内容之一，迫使地表温度产品真实性检验的应用需求与日俱增，引发了对地表温度采样布设方法研究的强烈关注。目前，遥感反演的地表温度产品空间分辨率大多介于百米至公里级范围内不等，如何有效、准确的布设地表温度的采样点，进而准确获取像元尺度的地表温度真值，实现对遥感反演的地表温度产品的精度验证已成为研究的热点。

传统的地表温度采样点的布设方法通常即未考虑下垫面的空间异质性，也未考虑地表温度的时间变化特征，只能适用于内陆湖泊、草地和雪地等均匀下垫面，而对于由多种不同下垫面构成的异质性地表，传统方法主要存在三个缺点：第一，采样点的个数选择上具有较大的随意性，而且采样点的位置也简单通过均匀分布或者随机分布的方式来确定，导致采样点的布设方案缺乏一定的科学性；第二，采样点的空间代表性有限，对于均匀下垫面，布设的采样点尚能代表像元尺度的地表温度真值，而对于异质下垫面，则无法准确表征像元尺度的地表温度真值；第三，采样点的时间代表性也有限，不能保证布设的采样点能准确代表任意时刻的像元尺度的地表温度真值。

因此，受制于地表温度时空多变的特点，当前地表温度的采样优化布设方案仅能够用于均匀下垫面的验证场，而对于由多种不同下垫面构成的异质性地表，无法准确表征像元尺度的地表温度真值，不利于实现像元尺度地表温度的精确测量。

技术实现要素：

为了解决上述技术所存在的不足之处，本发明提供了一种异质下垫面地表温度全天时采样优化布设方法。

为了解决以上技术问题，本发明采用的技术方案是：一种异质下垫面地表温度全天时采样优化布设方法，包括以下步骤：

一、相关产品的获取与预处理；

分别获取研究区的地表覆盖产品以及对应的地表温度产品；所述地表覆盖产品需选择高分辨率且下垫面地表类型没有发生明显变化的时期；所述地表温度产品需同时获取瞬时高空间分辨率极轨卫星的地表温度产品和高时间分辨率静止气象卫星的地表温度产品；

对获取的地表覆盖产品以及对应的地表温度产品进行预处理，使所有的产品具有相同的空间分辨率和覆盖范围；

二、静止气象卫星逐像元地表温度日周期变化曲线的拟合；

针对预处理完的高时间分辨率的地表温度产品，逐像元的进行温度日周期变化曲线的拟合，用以获取描述地表温度日周期变化的关键参数；

三、高空间分辨率极轨卫星逐小时地表温度产品的重建；

融合静止气象卫星和极轨卫星各自优点，重建出高空间分辨率极轨卫星逐小时地表温度产品，实现高空间和高时间分辨率的地表温度产品的获取；

四、地表温度采样代价函数的构建；

通过比较采样点和整个研究区的地表覆盖类型直方图的差异以及不同时刻地表温度直方图的差异，同时考虑采样点的空间聚集分布特征来构建采样代价函数cf；当采样代价函数取最小值的时候，采样点的布设方式在时空范围的代表性达到最佳；

五、采样样本数量的设置；

采用迭代遍历的方式确定最佳的采样样本数量；

首先，判断是否是初次设置采样样本的数量；如果是，则将采样样本数量设置为研究区地表覆盖类型的总个数，即sn＝c，其中sn为当前设置的采样样本总数；c为研究区地表覆盖类型的总个数；如果否，则将采样样本数量设置为前一次采样样本数量加c，即sn＝sn′+c，其中sn′为前一次设置的采样样本总数；

六、基于模拟退火算法的最佳样本布设位置确定；

在给定的采样样本数量下，基于模拟退火算法寻找全局最佳的采样点布设方案，确定最佳的样本布设位置；

七、迭代停止条件判断；

进行迭代停止条件判断，如果当前设置的采样样本总数sn大于某一给定阈值，则采样数量停止增长，转至步骤八，否则返回步骤五；

八、研究区最佳采样样本数量和位置的确定；

绘制采样样本数量sn与对应最佳采样方案的地表温度采样代价函数值cf的变化曲线图，分析地表温度采样代价函数值cf随采样样本数量sn的变化情况，选择地表温度采样代价函数值cf减小速率缓慢时对应的采样样本数量sn为最佳的采样样本数量，其对应的基于模拟退火算法所确定的采样样本布设位置为最佳的采样样本的布设位置；两者共同形成异质性下垫面地表温度全天时采样优化布设方案。

进一步地，步骤二的具体过程为：

将地表温度日周期变化曲线用六参数进行描述；白天变化利用余弦形式进行表征，而晚上变化则利用指数形式进行表征，数学公式表示为：

其中，lstgeo(t)是静止气象卫星在t时刻的地表温度产品值；t0为日出时刻附近的初始温度；ta为温度的幅度；tm为最大温度出现时刻；ts为温度开始衰减时刻；α为衰减系数；β为角频率；

由此可见，晴空无云情况下日周期变化的地表温度产品可通过t0、ta、tm、ts、α、β这六个参数来描述；根据地表温度产品获取的时间t，结合预处理完的逐小时的地表温度产品，利用最小二乘非线性数学回归方法率定出公式1所涉及的地表温度日周期变化曲线的关键六参数。

进一步地，步骤三的具体过程为：

考虑到静止气象卫星与极轨卫星存在通道设置以及地表温度反演算法间的差异，它们存在一个系统的缩放和偏移；将日周期变化的极轨卫星地表温度产品表示为：

lstleo(t)＝dtca×lstgeo(t)+dtcb公式2

其中，lstleo(t)是高空间分辨率的极轨卫星在t时刻的地表温度产品值；dtca、dtcb分别为极轨卫星地表温度和静止气象卫星地表温度间的系统性缩放系数、偏移系数，这两系数通过结合同一天不同时刻多次观测的瞬时高空间分辨率极轨卫星的地表温度产品和对应时刻的静止气象卫星的地表温度产品，利用最小二乘线性数学回归方法率定获得；

至此，利用最小二乘数学方法率定的地表温度日周期变化的六个参数和相应的系统缩放系数dtca、偏移系数dtcb，利用公式1和2重建出高空间分辨率极轨卫星逐小时地表温度产品。

进一步地，步骤四的具体过程为：

采样代价函数的数学表达式为：

其中，cf为地表温度采样代价函数；elc为采样点和整个研究区的地表覆盖类型直方图的差异；elst为采样点和整个研究区的不同时刻地表温度直方图的差异；ai为采样点的空间集聚分布特征指标；

elc的计算公式如下：

其中，c为研究区地表覆盖类型的总个数；为全部采样点中地表覆盖类型属于第i类土地覆盖的比例；为研究区中地表覆盖类型属于第i类土地覆盖的比例；i为序列数；

elst计算公式如下：

其中，sn为采样样本总数；为第k时刻采样的地表温度落入第j个地表温度层内的比例，为第k时刻整个研究区的地表温度落入第j个地表温度层内的比例；k和j为序列数；

ai的计算公式如下：

其中，sn为采样样本总数；min(ds)为第s个样本距离周围样本集合的最小距离；a为研究区的面积；s为序列数。

进一步地，步骤六的具体过程为：

a、在研究区随机选取sn个采样样本，模拟退火算法迭代次数设置为0；

b、计算首次采样点布设方案下地表温度采样代价函数cf1的值；

c、更改首次采样点布设方案，随机替换一个已选样本，并计算当前采样点布设方案下地表温度采样代价函数cf2的值；

d、判断是否接受当前采样点布设方案；如果cfk≤cfk-1，则接受当前采样布设方案；

反之，如果cfk＞cfk-1，继续进行如下判断：

首先计算误差概率，计算公式表示为：

其中，p为误差概率；f为衰减系数，f的初始值为1，每10次迭代后乘以0.95；随后，生成一个0到1之间的随机数r；最后通过比较误差概率p与随机数r之间的大小来判断是否接受当前采样布设方案；若p＞r，则接受当前采样点布设方案；反之，拒绝此次采样布设方案的更改，保留上一次迭代确定的采样点布设方案，模拟退火算法迭代次数加一；

e、继续迭代修改采样点布设方案；如果cfk≤cfk-1或者随机数r<0.5，则随机替换当前采样点中某个已选样本；否则，有选择性的替换已选样本，即对于每一时刻，统计采样的地表温度落入不同地表温度层内的比例，将比例最大值对应的地表温度分组中的所有样本挑选出来形成样本集合；将所有时刻挑选的样本集合求交集，若交集不为空，则随机替换交集中的某个已选样本；若交集为空，则对所有时刻挑选出来的样本集合中的样本分别统计其频数，随机替换样本频数最大的某个已选样本；

f、重复步骤d～e，直至地表温度采样代价函数值cf小于某一固定阈值或者模拟退火算法迭代次数超过某一固定阈值，则停止迭代；

g、保存地表温度采样代价函数值cf值最小的采样点布设结果作为某一采样样本数量设置下的最佳样本布设方案。

本发明利用模拟退火算法，遍历寻找最佳的采样点个数和位置，实现了采样点确定的科学性和合理性；通过对样本点与研究区地表温度空间变化的一致性约束，克服了采样点空间代表性有限的弊端；通过对样本点与研究区地表温度时间变化的一致性约束，解决了采样点时间代表性的问题；适用于异质下垫面任意时刻地表温度的最佳采样优化布设，能够全天时的准确获取像元尺度的地表温度真值。

附图说明

图1为本发明的整体流程示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1所示的一种异质下垫面地表温度全天时采样优化布设方法，包括以下步骤：

一、相关产品的获取与预处理；

分别获取研究区的地表覆盖产品以及对应的地表温度产品，其中地表覆盖产品需选择高分辨率且下垫面地表类型没有发生明显变化的时期，以确保地表覆盖数据类型的代表性；而地表温度产品需同时获取瞬时高空间分辨率极轨卫星的地表温度产品(例如，瞬时90m空间分辨率的aster地表温度产品)和高时间分辨率静止气象卫星的地表温度产品(例如，逐小时的5km空间分辨率的fy地表温度产品)。地表覆盖产品以及对应的地表温度产品获取完成后，需对其进行预处理，包括投影转换、重采样以及裁剪，使所有的产品具有相同的空间分辨率和覆盖范围，其中地表覆盖产品的重采样需采用最邻近采用方法以保证原始的地表覆盖类型值不发生改变。

二、静止气象卫星逐像元地表温度日周期变化曲线的拟合；

针对预处理完的高时间分辨率的地表温度产品，逐像元的进行温度日周期变化曲线的拟合，用以获取描述地表温度日周期变化的关键参数。

地表温度日周期变化曲线可以用六参数进行描述，即白天变化可以利用余弦形式进行表征，而晚上变化则可以利用指数形式进行表征，其数学公式可以表示为：

其中，lstgeo(t)是静止气象卫星在t时刻的地表温度产品值；t0为日出时刻附近的初始温度；ta为温度的幅度；tm为最大温度出现时刻；ts为温度开始衰减时刻；α为衰减系数；β为角频率；cos为数学上余弦运算符号；exp为数学上指数运算符号。由此可见，晴空无云情况下日周期变化的地表温度产品可以通过t0、ta、tm、ts、α、β这六个参数来描述。根据地表温度产品获取的时间t，结合预处理完的逐小时的地表温度产品，利用最小二乘非线性数学回归方法即可率定出公式1所涉及的地表温度日周期变化曲线的关键六参数。

三、高空间分辨率极轨卫星逐小时地表温度产品的重建；

由于静止气象卫星时间分辨率高，因此其能准确反映地表温度的时间变化特征，但是由于其空间分辨率低，损失了地表温度空间变化特征。然而，极轨卫星的地表温度温度产品更好相反，其空间分辨率高，但却缺乏逐小时地表温度时间变化信息的表征能力。融合静止气象卫星和极轨卫星各自优点，即可实现高空间和高时间分辨率的地表温度产品的获取。两种卫星的地表温度产品的日周期变化趋势(形状)近似一致，但考虑到静止气象卫星与极轨卫星存在通道设置以及地表温度反演算法间的差异，它们存在一个系统的缩放和偏移，即日周期变化的极轨卫星地表温度产品可以表示为：

lstleo(t)＝dtca×lstgeo(t)+dtcb公式2

其中，lstleo(t)是高空间分辨率的极轨卫星在t时刻的地表温度产品值；lstgeo(t)是静止气象卫星在t时刻的地表温度产品值；dtca和dtcb为极轨卫星地表温度和静止气象卫星地表温度间的系统性缩放和偏移系数。同样，结合同一天不同时刻多次观测的瞬时高空间分辨率极轨卫星的地表温度产品和对应时刻的静止气象卫星的地表温度产品，利用最小二乘线性数学回归方法即可率定出公式2所涉及的dtca和dtcb两系数。至此，利用最小二乘数学方法率定的地表温度日周期变化的六个参数(t0、ta、tm、ts、α、β)和相应的系统缩放系数(dtca)、偏移系数(dtcb)，利用公式1和2就能重建出高空间分辨率极轨卫星逐小时地表温度产品。

四、地表温度采样代价函数的构建；

通过比较采样点和整个研究区的地表覆盖类型直方图的差异以及不同时刻地表温度直方图的差异，同时考虑采样点的空间聚集分布特征来构建采样代价函数cf。当采样代价函数取最小值的时候，采样点的布设方式能够在时空范围的代表性达到最佳。采样代价函数的数学表达式表示为：

其中，cf为地表温度采样代价函数；elc为采样点和整个研究区的地表覆盖类型直方图的差异；elst为采样点和整个研究区的不同时刻地表温度直方图的差异；ai为采样点的空间集聚分布特征指标。

elc的计算公式如下：

其中，c为研究区地表覆盖类型的总个数；为全部采样点中地表覆盖类型属于第i类土地覆盖的比例；为研究区中地表覆盖类型属于第i类土地覆盖的比例；i为序列数。

elst计算公式如下：

其中，sn为采样样本总数，统计研究区不同时刻地表温度的变化范围，按变化范围将地表温度从小到大等温度间隔分为sn层。为第k时刻采样的地表温度落入第j个地表温度层内的比例，为第k时刻整个研究区的地表温度落入第j个地表温度层内的比例。k和j为序列数。

ai的计算公式如下：

其中，sn为采样样本总数；min(ds)为第s个样本距离周围样本集合的最小距离；a为研究区的面积；s为序列数。

五、采样样本数量的设置；

为了较准确的确定最佳的采样样本数量，这里采用迭代遍历的方式。首先，判断是否是初次设置采样样本的数量。如果是初次设置，则将采样样本数量设置为研究区地表覆盖类型的总个数，即sn＝c，其中sn为当前设置的采样样本总数；c为研究区地表覆盖类型的总个数。如果是非初次设置，则将采样样本数量设置为前一次采样样本数量加c，即sn＝sn′+c，其中sn′为前一次设置的采样样本总数。

六、基于模拟退火算法的最佳样本布设位置确定；

在给定的采样样本数量下，基于模拟退火算法寻找全局最佳的采样点布设方案，确定最佳的样本布设位置。具体的基于模拟退火方法确定最佳采样点布设的实施过程分为以下步骤：

a、在研究区随机选取sn个采样样本，模拟退火算法迭代次数设置为0；

b、计算首次采样点布设方案下地表温度采样代价函数cf1的值；

c、更改首次采样点布设方案，随机替换一个已选样本，并计算当前采样点布设方案下地表温度采样代价函数cf2的值；

d、判断是否接受当前采样点布设方案；如果cfk≤cfk-1(cfk表示当前的地表温度采样代价函数值；cfk-1表示前一次的地表温度采样代价函数值；k代表当前迭代，k-1代表前一次迭代)，即接受当前采样布设方案；

反之，如果cfk＞cfk-1，继续进行如下判断：首先计算误差概率，其计算公式可以表示为：

其中，p为误差概率；f为衰减系数，f的初始值为1，每10次迭代后乘以0.95。随后，生成一个0到1之间的随机数r。最后通过比较误差概率p与随机数r之间的大小来判断是否接受当前采样布设方案。若p＞r，则接受当前采样点布设方案；反之，拒绝此次采样布设方案的更改，保留上一次迭代确定的采样点布设方案，模拟退火算法迭代次数加一。

e、继续迭代修改采样点布设方案；如果cfk≤cfk-1或者随机数r<0.5，则随机替换当前采样点中某个已选样本；否则，有选择性的替换已选样本，即对于每一时刻，统计采样的地表温度落入不同地表温度层内的比例，将比例最大值对应的地表温度分组中的所有样本挑选出来形成样本集合；将所有时刻挑选的样本集合求交集，若交集不为空，则随机替换交集中的某个已选样本；若交集为空，则对所有时刻挑选出来的样本集合中的样本分别统计其频数，随机替换样本频数最大的某个已选样本；

f、重复步骤d～e直至地表温度采样代价函数值cf小于某一固定阈值(例如0.1)或者模拟退火算法迭代次数超过某一固定阈值(例如20000次)则停止迭代；

g、保存地表温度采样代价函数值cf值最小的采样点布设结果作为某一采样样本数量设置下的最佳样本布设方案。

七、迭代停止条件判断；

进行迭代停止条件判断，如果当前设置的采样样本总数sn大于某一给定阈值，例如研究区所有像元点总体数量的1％，则采样数量停止增长，转至步骤八，否则返回步骤五。

八、研究区最佳采样样本数量和位置的确定；

绘制采样样本数量sn与对应最佳采样方案的地表温度采样代价函数值cf的变化曲线图，分析地表温度采样代价函数值cf随采样样本数量sn的变化情况，选择地表温度采样代价函数值cf减小速率缓慢时(如减小速率为0.02)对应的采样样本数量sn为最佳的采样样本数量，其对应的基于模拟退火算法所确定的采样样本布设位置为最佳的采样样本的布设位置。两者共同形成了异质性下垫面地表温度全天时采样优化布设方案。

本发明的技术关键点以及优势在于：

1)根据下垫面地表覆盖类型确定采样点数量的初始值，并在此基础上采用遍历的方式按一定倍率逐步扩大采样点数量，并针对每一采样点数量，利用模拟退火算法确定最佳的采样位置，最终根据代价函数随采样点数量变化曲线给出适用于异质下垫面的地表温度全天时最佳的采样点个数和位置；

2)通过将采样点和整个研究区的地表覆盖类型直方图的一致性作为判读条件，挑选与研究区整体地表温度空间变化统计分布特征相一致的采样点，增加了采样点对地表温度空间变化的代表性，实现了采样点对地表温度空间变化的有效表征；

3)通过重构地表温度日周期变化曲线，将采样点和整个研究区的不同时刻地表温度直方图的一致性作为判读条件，挑选与研究区整体地表温度时间变化统计分布特征相一致的采样点，增加了采样点对地表温度时间变化的代表性，实现了采样点对地表温度时间变化的有效表征。

本发明明主要解决的技术问题为；

一、考虑了不同下垫面地表温度时空变化的差异，通过构建代价函数遍历寻找最佳的方式，解决了采样点个数选择的随意性问题；

二、通过挑选与研究区整体地表温度空间变化统计分布特征相一致的采样点，消除了地表温度空间异质性差异对采样点布设的影响，解决了传统采样点位置布设引发的空间代表性有限的问题；

三、通过温度日周期变化模型表征全天时地表温度，进一步挑选与研究区地表温度时间变化统计分布特征相一致的采样点，消除了地表温度时间变化差异对采样点布设的影响，解决了某一时刻最佳采样位置不适于另一时刻的弊端；

四、通过空间统计分布特征和时间统计分布特征双重标准的迭代筛选，同时满足时间和空间上的最佳，形成适用于异质下垫面的地表温度全天时的最佳采样优化布设方案，确定了最佳采样点个数以及采样点空间布设位置。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。