一种极地反照率遥感数据的重构方法、装置及计算机设备

摘要

本发明公开一种极地反照率遥感数据的重构方法，方法包括以下步骤：获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，对空缺像元的反照率数值进行填补，生成拟合数据序列，计算重构的反照率数值，修正的反照率数值。本发明充分利用遥感冰雪反照率数据中的时空信息填补数据的缺失，并考虑冰雪反照率数据的时序依赖性采用Whittaker滤波法重建晴空反照率时序数据，最后基于冰雪反照率受云层的影响的物理特性，根据经验公式进一步修正云覆盖的对应冰雪覆盖的像元的反照率，该方法重建极地冰雪反照率简单有效，重建结果精度高。

说明书

技术领域

本发明涉及遥感技术领域，尤其涉及一种极地反照率遥感数据的重构方法、装置及计算机设备。

背景技术

极地地区常年覆盖有全球最大面积的冰雪地表，由于冰雪高反照率的特性，其微小的变化就可能对地气系统的能量收支形成正反馈作用，进一步加强地球变暖。因此，研究极地反照率的时空变化具有重要意义。由于极地环境的特殊性，开展大范围、长时序冰雪反照率的实地观测是极其困难的，光学卫星遥感提供了在大时空尺度上连续观测反照率的途径。

南北极上空常有层状云覆盖，北纬85°处6～8月总云量为80～90％，多为云底高度低于1公里、厚度350～500米的低云。云的遮挡导致光学遥感观测的地表信息出现缺失，严重影响遥感卫星图像的利用率，给后期的处理和应用带来困难。因此需要重建光学遥感卫星影像云覆盖区域的冰雪反照率信息。

发明内容

为了解决以上问题，本发明实施例提供一种极地反照率遥感数据的重构方法，具体包括以下步骤：

步骤S1：获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列；

步骤S2：根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值；

步骤S3：将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值；

步骤S4：分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值；

步骤S5：对于每个像元，计算权重W，权重W的计算公式为：

利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：

对所述拟合指数F进行迭代计算，当拟合指数F符合条件F

步骤S6：判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。

进一步地，所述步骤S2中，空缺像元的填补函数为：

其中，(m,n)代表空间上第(m,n)个像元的二维序号，所述像元具有固定的行列号，i代表时序为i的像元的序号，num表示相邻时空有效像元的数量，y

进一步地，所述步骤S3中，衡量所述第一像元组与所述第二像元组的反照率数值时序数据的相似性使用皮尔逊相关系数，其计算公式为：

其中，X表示第一像元组的反照率数值时序数据，Y表示第二像元组的反照率数值时序数据，ρ(X,Y)表示第一像元组与第二像元组的反照率数值时序数据的皮尔逊相关系数，E表示数学期望，D表示方差，con(X,Y)表示协方差。

优选地，ρ(X,Y)大于0.9时，所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件。

优选地，所述步骤S4中，拟合数据序列[Zi]与填补数据序列[Yi]的拟合关系式如下：

Z

其中，I为单位矩阵，λ为参数，Dz表示z

优选地，λ设置为2。

优选地，所述步骤S6中，对该空缺像元的所述重构的反照率数值进行修正的计算公式如下：

其中，

优选地，a、b、c、d分别设置为-0.0491243、1.06756、0.0217075、0.017950。

本发明实施例还提供一种极地反照率遥感数据的重构装置，包括：

数据获取模块，用于获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列；

数值填补模块，用于根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值；

数值替换模块，用于将数值填补模块中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值；

数值拟合模块，用于分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值；

迭代计算模块，用于对于每个像元，计算权重W，权重W的计算公式为：

利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：

对所述拟合指数F进行迭代计算，当拟合指数F符合条件F

修正模块，用于判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。

本发明实施例还提供一种计算机设备，包括处理器和存储器，存储在存储器上并可在处理器上运行的极地反照率遥感数据的重构程序，该重构程序被处理器执行时实现如上所述的重构方法的步骤。

与现有技术相比，本发明的有益效果包括：本发明充分利用遥感冰雪反照率数据中的时空信息填补数据的缺失，并考虑冰雪反照率数据的时序依赖性采用Whittaker滤波法重建晴空反照率时序数据，最后基于冰雪反照率受云层的影响的物理特性，根据经验公式进一步修正云覆盖的对应冰雪覆盖的像元的反照率，该方法重建极地冰雪反照率简单有效，重建结果精度高。

附图说明

图1是本发明重构方法的流程图；

图2本发明的重构装置的结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例1

本发明提供了一种极地反照率遥感数据的重构方法。请参考图1，图1是本发明重构方法的流程图；本发明方法包括以下步骤：

步骤S1：获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列。

具体的，南北极区域反照率的遥感数据来源于现有遥感数据集MOD10A1，其中的反照率的遥感数据为三维数据，即空间上的栅格数据和时间上的时序数据。该反照率的遥感数据在空间上的分辨率为500m，即一个栅格的大小为500m*500m，一个栅格对应一个反照率时序数据，该时序数据的时间间隔是1天，都是一天内相同时刻的遥感测量数据。由于云层是移动的，一个栅格对应的空间位置在某天可能被云层遮盖，在另一天可能未被云层遮盖。某一时间某个栅格对应的空间位置未被云层遮盖时，该时间的栅格为有效像元，该有效像元具有一个反照率数据值；某一时间某个栅格对应的空间被云层遮盖时，该时间的栅格为空缺像元，该空缺像元的反照率数据值为缺省的空值，需要进行填补重构。现有遥感数据集MOD10A1目前具有20多年的观测值，本发明的方法可以对该数据集所有的空缺像元的反照率数据进行重构。可以理解的是，本发明的方法也可以对该数据集的某一时间段的空缺像元的反照率数据进行重构，如近10年、5年、3年等。

云光学厚度的遥感数据来源于现有遥感数据集AVHRR Polar Pathfinder，其中的云光学厚度的遥感数据也为三维数据，即空间上的栅格数据和时间上的时序数据。该云光学厚度的遥感数据在空间上的分辨率为25km，即一个栅格的大小为25km*25km，一个栅格对应一个云光学厚度时序数据，该时序数据的时间间隔是1天，都是一天内相同时刻的遥感测量数据，一天中云光学厚度的遥感数据对应的时刻与反照率的遥感数据对应的时刻相差3个小时左右。

太阳天顶角的遥感数据来源于现有遥感数据集MOD09GA，其中的太阳天顶角的遥感数据也为三维数据，即空间上的栅格数据和时间上的时序数据。该太阳天顶角的遥感数据在空间上的分辨率为500m，即一个栅格的大小为500m*500m，一个栅格对应一个太阳天顶角时序数据，该时序数据的时间间隔是1天，都是一天内相同时刻的遥感测量数据，一天中太阳天顶角的遥感数据对应的时刻与反照率的遥感数据对应的时刻相同。

统一所述遥感数据的空间分辨率时，将云光学厚度的遥感数据重采样至与反照率的遥感数据具有相同的空间分辨率，得到云光学厚度时序数据，该云光学厚度时序数据在空间上的分辨率为500m，即一个栅格的大小为500m*500m，可理解为将大栅格划分为小栅格，小栅格的云光学厚度时序数据为原来大栅格的云光学厚度时序数据。反照率时序数据和太阳天顶角时序数据为初始获取的南北极区域反照率和太阳天顶角的遥感数据。

应当理解的是，也可以根据需要对反照率、云光学厚度和太阳天顶角的遥感数据都进行重采样，将它们统一到一个更小的空间分辨率上。

步骤S2：根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值。

具体的，某一个像元的填补函数为：

其中，(m,n)代表空间上第(m,n)个像元的二维序号，所述像元具有固定的行列号，i代表时序为i的像元的序号，num表示相邻时空有效像元的数量，y

通过上述填补函数按照时空顺序进行循环填补，便可以将所有空缺像元缺省的反照率数值填补完全。

步骤S3：将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值。

空缺像元缺省的反照率数值通过时空信息填补之后，考虑到冰雪反照率数据具有很强的时序依赖性，本发明通过计算空缺像元与相同时间上的对比像元分别组成的连续像元组是否具有时序相似性，如果有则用所述对比像元的反照率数值替换所述目标像元的反照率数值，能够进一步提高空缺像元的反照率数值重构的准确性。

第一像元组的多个像元在空间位置上相同，第二像元组的多个像元在空间位置上相同，但第一像元组的多个像元在空间位置上与第二像元组的多个像元在空间位置上不同，目标像元在时间上与对比像元为同一天，第一像元组的多个像元在时序上与第二像元组的多个像元在时序上分别一一对应。

优选的，衡量所述第一像元组与所述第二像元组的反照率数值时序数据的相似性使用皮尔逊相关系数，其计算公式为：

其中，X表示第一像元组的反照率数值时序数据，Y表示第二像元组的反照率数值时序数据，ρ(X,Y)表示第一像元组与第二像元组的反照率数值时序数据的皮尔逊相关系数，E表示数学期望，D表示方差，con(X,Y)表示协方差。

例如，X为目标像元前5天至后4天共10天的反照率数值时序数据，Y为对比像元前5天至后4天共10天的反照率数值时序数据。

优选的，ρ(X,Y)大于0.9时，所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件。

通过上述将步骤S2中所有填补的空缺像元分别计算是否存在时序相似性的有效像元，若存在，对相应的空缺像元的反照率数值用相应的对比像元的反照率数值进行覆盖。

步骤S4：分别选取所述空缺像元所在空间上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值。

具体的，Yi表示某个栅格上第i天的反照率数值，Zi表示该栅格对应的第i天的反照率数值的拟合数值。

空缺像元所在空间位置上的所有像元指与空缺像元位于同一栅格上的不同时间上的像元，包括有效像元和空缺像元，当一个空间位置上存在多个空缺像元时，只生成一个填补数据序列[Yi]。对不同空间位置上的栅格分别采用步骤S4进行滤波，可以得到所有空间位置上反照率数值的拟合数据序列[Zi]。

根据Whittaker算法，S＝∑

其中，S、R分别表示保真度和粗糙度。

记：Dz＝z

其中，Dz表示z

记：Q＝min(S+λR)＝min(∑

其中，Q表示拟合效果。

将Q对z

其中，Dz′表示Dz对z

根据上式求得拟合数据序列[Zi]与填补数据序列[Yi]的拟合关系式如下：

Z

其中，I为单位矩阵，λ为参数，可根据拟合需要进行调整，优选的设置为2。

步骤S5：对于每个像元，计算权重W，权重W的计算公式为：

利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：

对所述拟合指数F进行迭代计算，当拟合指数F符合条件F

步骤S2、S3根据时空信息填补的反照率时序数据和Whittaker算法拟合的反照率时序数据在一定程度上都反映了反照率的长期变化，在重构反照率数据时同时考虑填补数据和拟合数据，并根据拟合数据和填补数据的差值计算反照率时序数据中每个数据的权重；其中，差值越小的数据设置越高的权重，说明拟合数据中对应的数据质量越高，重建反照率数值需要最大程度接近填补数据和拟合数据，即最大程度靠近权重高的数据，因此使用拟合指数进行迭代进一步提高重建的反照率数值的准确性。

上述步骤会对每个像元计算相应的重构的反照率数值

步骤S6：判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。

以上步骤S2-S5重建的是像元对应的区域假设未被云层遮盖时的晴空反照率，极地区域包括冰雪覆盖区域和陆地覆盖区域，考虑云层会增加冰雪的反照率，但不影响陆地的反照率，因此对应被冰雪覆盖区域的空缺像元的反照率数值需要进行修正。

计算修正的反照率数值需要先判断空缺像元对应的区域是否被冰雪覆盖，本发明中判断的方式为若该像元相同空间位置上在前一时序上的像元为冰雪，则判断该空缺像元对应的区域被冰雪覆盖，具体的，若该像元相同空间上前一天的像元为冰雪，则判断该空缺像元对应的区域被冰雪覆盖；反之，判断该空缺像元对应的区域为陆地。现有遥感数据集MOD10A1的有效像元的反照率数据值均有对应的属性，如陆地和冰雪覆盖，该属性是遥感卫星根据遥感数据光谱分析得到的，在数据集MOD10A1已有标定可以直接使用。

其中，利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值指步骤S1中得到的云光学厚度数值序列及太阳天顶角数值序列中，与空缺像元对应相同空间和时间的云光学厚度数值和太阳天顶角数值，如同一天同一位置的遥感云光学厚度数值和太阳天顶角数值。

对该空缺像元的所述重构的反照率数值进行修正的计算公式如下：

其中，

本发明充分利用遥感冰雪反照率数据中的时空信息填补数据的缺失，并考虑冰雪反照率数据的时序依赖性采用Whittaker滤波法重建晴空反照率时序数据，最后基于冰雪反照率受云层的影响的物理特性，根据经验公式进一步修正云覆盖的对应冰雪覆盖的像元的反照率，该方法重建极地冰雪反照率简单有效，重建结果精度高。

实施例2

本发明提供了一种极地反照率遥感数据的重构装置，请参考图2，图2是本发明的重构装置100的结构示意图；本发明的重构装置100包括：

数据获取模块110，用于获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列。

数值填补模块120，用于根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值。

数值替换模块130，用于将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值。

数值拟合模块140，用于分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值。

迭代计算模块150，用于对于每个像元，计算权重W，权重W的计算公式为：

利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：

对所述拟合指数F进行迭代计算，当拟合指数F符合条件F

修正模块160，用于判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。

实施例3

本发明还提供一种计算机设备，包括处理器和存储器，存储在存储器上并可在处理器上运行的极地反照率遥感数据的重构程序，该重构程序被处理器执行时实现上述极地反照率遥感数据的重构方法实施例的各个过程，且能达到相同的技术效果，为避免重复，这里不再赘述。

以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所做出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。