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

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

Info

Publication number
CN114936202B
CN114936202B CN202210527353.1A CN202210527353A CN114936202B CN 114936202 B CN114936202 B CN 114936202B CN 202210527353 A CN202210527353 A CN 202210527353A CN 114936202 B CN114936202 B CN 114936202B
Authority
CN
China
Prior art keywords
albedo
pixels
value
data
pixel
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210527353.1A
Other languages
English (en)
Other versions
CN114936202A (zh
Inventor
程青
叶凡
张泽君
张景航
王力哲
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Geosciences
Original Assignee
China University of Geosciences
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Geosciences filed Critical China University of Geosciences
Priority to CN202210527353.1A priority Critical patent/CN114936202B/zh
Publication of CN114936202A publication Critical patent/CN114936202A/zh
Application granted granted Critical
Publication of CN114936202B publication Critical patent/CN114936202B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/21—Design, administration or maintenance of databases
    • G06F16/215—Improving data quality; Data cleansing, e.g. de-duplication, removing invalid entries or correcting typographical errors
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2462—Approximate or statistical queries
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2474—Sequence data queries, e.g. querying versioned data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Mathematical Physics (AREA)
  • Fuzzy Systems (AREA)
  • Computational Linguistics (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Quality & Reliability (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

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

Description

一种极地反照率遥感数据的重构方法、装置及计算机设备
技术领域
本发明涉及遥感技术领域，尤其涉及一种极地反照率遥感数据的重构方法、装置及计算机设备。
背景技术
极地地区常年覆盖有全球最大面积的冰雪地表，由于冰雪高反照率的特性，其微小的变化就可能对地气系统的能量收支形成正反馈作用，进一步加强地球变暖。因此，研究极地反照率的时空变化具有重要意义。由于极地环境的特殊性，开展大范围、长时序冰雪反照率的实地观测是极其困难的，光学卫星遥感提供了在大时空尺度上连续观测反照率的途径。
南北极上空常有层状云覆盖，北纬85°处6～8月总云量为80～90％，多为云底高度低于1公里、厚度350～500米的低云。云的遮挡导致光学遥感观测的地表信息出现缺失，严重影响遥感卫星图像的利用率，给后期的处理和应用带来困难。因此需要重建光学遥感卫星影像云覆盖区域的冰雪反照率信息。
发明内容
为了解决以上问题，本发明实施例提供一种极地反照率遥感数据的重构方法，具体包括以下步骤：
步骤S1：获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列；
步骤S2：根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值；
步骤S3：将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值；
步骤S4：分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值；
步骤S5：对于每个像元，计算权重W，权重W的计算公式为：
Figure BDA0003645088310000021
利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：
Figure BDA0003645088310000022
对所述拟合指数F进行迭代计算，当拟合指数F符合条件Fk-1≥Fk≤Fk+1时退出迭代，计算得到该像元的重构的反照率数值
Figure BDA0003645088310000023
其中k表示迭代次数，k≥2，/>
Figure BDA0003645088310000024
步骤S6：判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。
进一步地，所述步骤S2中，空缺像元的填补函数为：
Figure BDA0003645088310000031
其中，(m,n)代表空间上第(m,n)个像元的二维序号，所述像元具有固定的行列号，i代表时序为i的像元的序号，num表示相邻时空有效像元的数量，yM,N,I表示空缺像元相邻时空的某个有效像元的反照率数值，ym,n,i表示被填补的空缺像元的反照率数值。
进一步地，所述步骤S3中，衡量所述第一像元组与所述第二像元组的反照率数值时序数据的相似性使用皮尔逊相关系数，其计算公式为：
Figure BDA0003645088310000032
其中，X表示第一像元组的反照率数值时序数据，Y表示第二像元组的反照率数值时序数据，ρ(X,Y)表示第一像元组与第二像元组的反照率数值时序数据的皮尔逊相关系数，E表示数学期望，D表示方差，con(X,Y)表示协方差。
优选地，ρ(X,Y)大于0.9时，所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件。
优选地，所述步骤S4中，拟合数据序列[Zi]与填补数据序列[Yi]的拟合关系式如下：
Zi＝(I+λDz′Dz)-1Yi
其中，I为单位矩阵，λ为参数，Dz表示zi的二阶差分，Dz＝zi-3zi-1+3zi-2-zi-3，Dz′表示Dz对zi的导数。
优选地，λ设置为2。
优选地，所述步骤S6中，对该空缺像元的所述重构的反照率数值进行修正的计算公式如下：
Figure BDA0003645088310000041
其中，
Figure BDA0003645088310000042
τ、S分别表示修正的反照率数值、重构的反照率数值、云光学厚度数值、太阳天顶角数值，a、b、c、d为回归系数。
优选地，a、b、c、d分别设置为-0.0491243、1.06756、0.0217075、0.017950。
本发明实施例还提供一种极地反照率遥感数据的重构装置，包括：
数据获取模块，用于获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列；
数值填补模块，用于根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值；
数值替换模块，用于将数值填补模块中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值；
数值拟合模块，用于分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值；
迭代计算模块，用于对于每个像元，计算权重W，权重W的计算公式为：
Figure BDA0003645088310000051
利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：
Figure BDA0003645088310000052
对所述拟合指数F进行迭代计算，当拟合指数F符合条件Fk-1≥Fk≤Fk+1时退出迭代，计算得到该像元的重构的反照率数值
Figure BDA0003645088310000053
其中k表示迭代次数，k≥2，/>
Figure BDA0003645088310000054
修正模块，用于判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。
本发明实施例还提供一种计算机设备，包括处理器和存储器，存储在存储器上并可在处理器上运行的极地反照率遥感数据的重构程序，该重构程序被处理器执行时实现如上所述的重构方法的步骤。
与现有技术相比，本发明的有益效果包括：本发明充分利用遥感冰雪反照率数据中的时空信息填补数据的缺失，并考虑冰雪反照率数据的时序依赖性采用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：根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值。
具体的，某一个像元的填补函数为：
Figure BDA0003645088310000081
其中，(m,n)代表空间上第(m,n)个像元的二维序号，所述像元具有固定的行列号，i代表时序为i的像元的序号，num表示相邻时空有效像元的数量，yM,N,I表示空缺像元相邻时空的某个有效像元的反照率数值，ym,n,i表示被填补的空缺像元的反照率数值。
通过上述填补函数按照时空顺序进行循环填补，便可以将所有空缺像元缺省的反照率数值填补完全。
步骤S3：将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值。
空缺像元缺省的反照率数值通过时空信息填补之后，考虑到冰雪反照率数据具有很强的时序依赖性，本发明通过计算空缺像元与相同时间上的对比像元分别组成的连续像元组是否具有时序相似性，如果有则用所述对比像元的反照率数值替换所述目标像元的反照率数值，能够进一步提高空缺像元的反照率数值重构的准确性。
第一像元组的多个像元在空间位置上相同，第二像元组的多个像元在空间位置上相同，但第一像元组的多个像元在空间位置上与第二像元组的多个像元在空间位置上不同，目标像元在时间上与对比像元为同一天，第一像元组的多个像元在时序上与第二像元组的多个像元在时序上分别一一对应。
优选的，衡量所述第一像元组与所述第二像元组的反照率数值时序数据的相似性使用皮尔逊相关系数，其计算公式为：
Figure BDA0003645088310000091
其中，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＝∑i(Yi-Zi)2，R＝∑i(Zi-3Zi-1+3Zi-2-Zi-3)2
其中，S、R分别表示保真度和粗糙度。
记：Dz＝zi-3zi-1+3zi-2-zi-3
其中，Dz表示zi的二阶差分。
记：Q＝min(S+λR)＝min(∑i(yi-zi)2+λ∑i(DZi)2)
其中，Q表示拟合效果。
将Q对zi求导得到：
Figure BDA0003645088310000101
其中，Dz′表示Dz对zi的导数。
根据上式求得拟合数据序列[Zi]与填补数据序列[Yi]的拟合关系式如下：
Zi＝(I+λDz′Dz)-1Yi
其中，I为单位矩阵，λ为参数，可根据拟合需要进行调整，优选的设置为2。
步骤S5：对于每个像元，计算权重W，权重W的计算公式为：
Figure BDA0003645088310000111
利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：
Figure BDA0003645088310000112
对所述拟合指数F进行迭代计算，当拟合指数F符合条件Fk-1≥Fk≤Fk+1时退出迭代，计算得到该像元的重构的反照率数值
Figure BDA0003645088310000113
其中k表示迭代次数，k≥2，/>
Figure BDA0003645088310000114
步骤S2、S3根据时空信息填补的反照率时序数据和Whittaker算法拟合的反照率时序数据在一定程度上都反映了反照率的长期变化，在重构反照率数据时同时考虑填补数据和拟合数据，并根据拟合数据和填补数据的差值计算反照率时序数据中每个数据的权重；其中，差值越小的数据设置越高的权重，说明拟合数据中对应的数据质量越高，重建反照率数值需要最大程度接近填补数据和拟合数据，即最大程度靠近权重高的数据，因此使用拟合指数进行迭代进一步提高重建的反照率数值的准确性。
上述步骤会对每个像元计算相应的重构的反照率数值
Figure BDA0003645088310000115
包括有效像元和空缺像元，其目的在于计算出准确的空缺像元的反照率数值，对于有效像元的重构的反照率数值可在本步骤之后恢复为遥感数据集MOD10A1初始的反照率数值。
步骤S6：判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。
以上步骤S2-S5重建的是像元对应的区域假设未被云层遮盖时的晴空反照率，极地区域包括冰雪覆盖区域和陆地覆盖区域，考虑云层会增加冰雪的反照率，但不影响陆地的反照率，因此对应被冰雪覆盖区域的空缺像元的反照率数值需要进行修正。
计算修正的反照率数值需要先判断空缺像元对应的区域是否被冰雪覆盖，本发明中判断的方式为若该像元相同空间位置上在前一时序上的像元为冰雪，则判断该空缺像元对应的区域被冰雪覆盖，具体的，若该像元相同空间上前一天的像元为冰雪，则判断该空缺像元对应的区域被冰雪覆盖；反之，判断该空缺像元对应的区域为陆地。现有遥感数据集MOD10A1的有效像元的反照率数据值均有对应的属性，如陆地和冰雪覆盖，该属性是遥感卫星根据遥感数据光谱分析得到的，在数据集MOD10A1已有标定可以直接使用。
其中，利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值指步骤S1中得到的云光学厚度数值序列及太阳天顶角数值序列中，与空缺像元对应相同空间和时间的云光学厚度数值和太阳天顶角数值，如同一天同一位置的遥感云光学厚度数值和太阳天顶角数值。
对该空缺像元的所述重构的反照率数值进行修正的计算公式如下：
Figure BDA0003645088310000121
其中，
Figure BDA0003645088310000122
τ、S分别表示修正的反照率数值、重构的反照率数值、云光学厚度数值、太阳天顶角数值，a、b、c、d为回归系数，根据经验分别设置为-0.0491243、1.06756、0.0217075、0.017950。
本发明充分利用遥感冰雪反照率数据中的时空信息填补数据的缺失，并考虑冰雪反照率数据的时序依赖性采用Whittaker滤波法重建晴空反照率时序数据，最后基于冰雪反照率受云层的影响的物理特性，根据经验公式进一步修正云覆盖的对应冰雪覆盖的像元的反照率，该方法重建极地冰雪反照率简单有效，重建结果精度高。
实施例2
本发明提供了一种极地反照率遥感数据的重构装置，请参考图2，图2是本发明的重构装置100的结构示意图；本发明的重构装置100包括：
数据获取模块110，用于获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列。
数值填补模块120，用于根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值。
数值替换模块130，用于将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值。
数值拟合模块140，用于分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值。
迭代计算模块150，用于对于每个像元，计算权重W，权重W的计算公式为：
Figure BDA0003645088310000141
利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：
Figure BDA0003645088310000142
对所述拟合指数F进行迭代计算，当拟合指数F符合条件Fk-1≥Fk≤Fk+1时退出迭代，计算得到该像元的重构的反照率数值
Figure BDA0003645088310000143
其中k表示迭代次数，k≥2，/>
Figure BDA0003645088310000144
修正模块160，用于判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值。
实施例3
本发明还提供一种计算机设备，包括处理器和存储器，存储在存储器上并可在处理器上运行的极地反照率遥感数据的重构程序，该重构程序被处理器执行时实现上述极地反照率遥感数据的重构方法实施例的各个过程，且能达到相同的技术效果，为避免重复，这里不再赘述。
以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所做出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。

Claims (7)

1.一种极地反照率遥感数据的重构方法，其特征在于，包括以下的步骤：
步骤S1：获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列；
步骤S2：根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值，空缺像元的填补函数为：
Figure FDA0004224919900000011
其中，(m,n)代表空间上第(m,n)个像元的二维序号，所述像元具有固定的行列号，i代表时序为i的像元的序号，num表示相邻时空有效像元的数量，yM,N,I表示空缺像元相邻时空的某个有效像元的反照率数值，ym,n,i表示被填补的空缺像元的反照率数值；
步骤S3：将步骤S2中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值；
步骤S4：分别选取所述空缺像元所在空间位置上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值；拟合数据序列[Zi]与填补数据序列[Yi]的拟合关系式如下：
Zi＝(I+λDz'Dz)-1Yi
其中，I为单位矩阵，λ为参数，Dz表示zi的二阶差分，Dz＝zi-3zi-1+3zi-2-zi-3，Dz′表示Dz对zi的导数；
步骤S5：对于每个像元，计算权重W，权重W的计算公式为：
Figure FDA0004224919900000021
利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：
Figure FDA0004224919900000022
对所述拟合指数F进行迭代计算，当拟合指数F符合条件Fk-1≥fk≤fk+1时退出迭代，计算得到该像元的重构的反照率数值
Figure FDA0004224919900000023
其中k表示迭代次数，k≥2，/>
Figure FDA0004224919900000024
步骤S6：判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值；对该空缺像元的所述重构的反照率数值进行修正的计算公式如下：
Figure FDA0004224919900000025
其中，
Figure FDA0004224919900000026
τ、S分别表示修正的反照率数值、重构的反照率数值、云光学厚度数值、太阳天顶角数值，a、b、c、d为回归系数。
2.如权利要求1所述的重构方法，其特征在于：所述步骤S3中，衡量所述第一像元组与所述第二像元组的反照率数值时序数据的相似性使用皮尔逊相关系数，其计算公式为：
Figure FDA0004224919900000031
其中，X表示第一像元组的反照率数值时序数据，Y表示第二像元组的反照率数值时序数据，ρ(X,Y)表示第一像元组与第二像元组的反照率数值时序数据的皮尔逊相关系数，E表示数学期望，D表示方差，con(X,Y)表示协方差。
3.如权利要求2所述的重构方法，其特征在于：ρ(X,Y)大于0.9时，所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件。
4.如权利要求1所述的重构方法，其特征在于：λ设置为2。
5.如权利要求1所述的重构方法，其特征在于：a、b、c、d分别设置为-0.0491243、1.06756、0.0217075、0.017950。
6.一种极地反照率遥感数据的重构装置，其特征在于，包括：
数据获取模块，用于获取南北极区域反照率、云光学厚度和太阳天顶角的遥感数据，统一所述遥感数据的空间分辨率，得到反照率数值序列、云光学厚度数值序列及太阳天顶角数值序列；
数值填补模块，用于根据空缺像元相邻时间和空间的有效像元的反照率数值的均值填补空缺像元的反照率数值，空缺像元的填补函数为：
Figure FDA0004224919900000032
其中，(m,n)代表空间上第(m,n)个像元的二维序号，所述像元具有固定的行列号，i代表时序为i的像元的序号，num表示相邻时空有效像元的数量，yM,N,I表示空缺像元相邻时空的某个有效像元的反照率数值，ym,n,i表示被填补的空缺像元的反照率数值；
数值替换模块，用于将数值填补模块中填补的空缺像元分别作为目标像元，选取包括该目标像元的时序上连续的多个像元构成第一像元组，选取与该目标像元时间上相同的有效像元作为对比像元，选取包括该对比像元的时序上连续的多个像元构成第二像元组，所述第一像元组与所述第二像元组的时序相同，计算所述第一像元组与所述第二像元组的反照率数值时序数据的相似性，若所述第一像元组与所述第二像元组的反照率数值时序数据满足相似性条件，则使用所述对比像元的反照率数值替换所述目标像元的反照率数值；
数值拟合模块，用于分别选取所述空缺像元所在空间上的所有像元的反照率数值时序数据构成填补数据序列[Yi]，利用Whittaker算法对所述填补数据序列Yi进行滤波，生成拟合数据序列[Zi]；Yi表示时序为i的像元的反照率数值，Zi表示对应Yi的拟合数值；拟合数据序列[Zi]与填补数据序列[Yi]的拟合关系式如下：
Zi＝(I+hDz'Dz)-1Yi
其中，I为单位矩阵，λ为参数，Dz表示zi的二阶差分，Dz＝zi-3zi-1+3zi-2-zi-3，Dz′表示Dz对zi的导数；
迭代计算模块，用于对于每个像元，计算权重W，权重W的计算公式为：
Figure FDA0004224919900000041
利用Whittaker算法构建拟合指数F，拟合指数F的计算公式为：
Figure FDA0004224919900000051
对所述拟合指数F进行迭代计算，当拟合指数F符合条件Fk-1≥fk≤fk+1时退出迭代，计算得到该像元的重构的反照率数值
Figure FDA0004224919900000055
其中k表示迭代次数，k≥2，/>
Figure FDA0004224919900000052
修正模块，用于判断空缺像元是否对应被冰雪覆盖的区域，若是则利用所述空缺像元对应的云光学厚度数值和太阳天顶角数值对该空缺像元的所述重构的反照率数值进行修正，得到修正的反照率数值；对该空缺像元的所述重构的反照率数值进行修正的计算公式如下：
Figure FDA0004224919900000053
其中，
Figure FDA0004224919900000054
τ、S分别表示修正的反照率数值、重构的反照率数值、云光学厚度数值、太阳天顶角数值，a、b、c、d为回归系数。
7.一种计算机设备，其特征在于，包括处理器和存储器，存储在存储器上并可在处理器上运行的极地反照率遥感数据的重构程序，该重构程序被处理器执行时执行如权利要求1-5中任一项所述的重构方法的步骤。
CN202210527353.1A 2022-05-16 2022-05-16 一种极地反照率遥感数据的重构方法、装置及计算机设备 Active CN114936202B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210527353.1A CN114936202B (zh) 2022-05-16 2022-05-16 一种极地反照率遥感数据的重构方法、装置及计算机设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210527353.1A CN114936202B (zh) 2022-05-16 2022-05-16 一种极地反照率遥感数据的重构方法、装置及计算机设备

Publications (2)

Publication Number Publication Date
CN114936202A CN114936202A (zh) 2022-08-23
CN114936202B true CN114936202B (zh) 2023-07-11

Family

ID=82864150

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210527353.1A Active CN114936202B (zh) 2022-05-16 2022-05-16 一种极地反照率遥感数据的重构方法、装置及计算机设备

Country Status (1)

Country Link
CN (1) CN114936202B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2590135A1 (en) * 2010-06-30 2013-05-08 Nec Corporation Color image processing method, color image processing device, and color image processing program
CN107782700A (zh) * 2017-09-07 2018-03-09 北京师范大学 一种avhrr地表反射率重建方法、系统与装置
CN111798394A (zh) * 2020-06-30 2020-10-20 电子科技大学 一种基于多年时间序列数据的遥感图像云污染去除方法
CN114463202A (zh) * 2022-01-13 2022-05-10 武汉大学 结合矩阵补全和趋势滤波的植被指数时间序列重建方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2590135A1 (en) * 2010-06-30 2013-05-08 Nec Corporation Color image processing method, color image processing device, and color image processing program
CN107782700A (zh) * 2017-09-07 2018-03-09 北京师范大学 一种avhrr地表反射率重建方法、系统与装置
CN111798394A (zh) * 2020-06-30 2020-10-20 电子科技大学 一种基于多年时间序列数据的遥感图像云污染去除方法
CN114463202A (zh) * 2022-01-13 2022-05-10 武汉大学 结合矩阵补全和趋势滤波的植被指数时间序列重建方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
On the performance of remote sensing time series reconstruction methods – A spatial comparison;Jie Zhou;《Remote Sensing of Environment》;368-384页 *

Also Published As

Publication number Publication date
CN114936202A (zh) 2022-08-23

Similar Documents

Publication Publication Date Title
Mémin et al. Correcting GPS measurements for non-tidal loading
CN106780321B (zh) 一种cbers-02卫星hr传感器影像整体严密定向与纠正拼接方法
Arowolo et al. Comparison of spatial interpolation techniques to generate high‐resolution climate surfaces for Nigeria
CN107040695B (zh) 基于rpc定位模型的星载视频稳像方法及系统
CN108921035B (zh) 基于空间引力和像元聚集度的亚像元定位方法和系统
CN116011342B (zh) 一种高分辨率热红外地表温度全天候重建方法
Li et al. Synergistic data fusion of multimodal AOD and air quality data for near real-time full coverage air pollution assessment
CN111650579B (zh) 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质
Grgić et al. Absolute sea level surface modeling for the Mediterranean from satellite altimeter and tide gauge measurements
CN115327662A (zh) 一种基于统计学习的区域电离层tec预测方法
CN111666896A (zh) 一种基于线性融合模型的遥感影像时空融合方法
CN115629384A (zh) 一种时序InSAR误差的改正方法及相关设备
CN109188483B (zh) 一种时序化高精度外方位元素自动定标方法
CN105571598B (zh) 一种卫星激光高度计足印相机姿态的测定方法
CN114936202B (zh) 一种极地反照率遥感数据的重构方法、装置及计算机设备
CN107516291B (zh) 夜景影像正射纠正处理方法
CN112711022A (zh) 一种GNSS层析技术辅助的InSAR大气延迟改正方法
CN116029162B (zh) 利用星载gnss-r数据的洪涝灾害淹没范围监测方法和系统
CN114970222B (zh) 基于hasm的区域气候模式日平均气温偏差订正方法和系统
TWI717796B (zh) 以人工智慧估測地理位置日照量之系統、方法及儲存媒體
CN115755103A (zh) 一种抗差自适应的gnss水汽层析方法
CN113255148A (zh) 基于modis产品数据估算全天气气温及其时空分布方法
KR20090088131A (ko) 표고편차를 고려한 국지기온의 추정방법 및 그 추정시스템
CN116776651B (zh) 地表蒸散发测算方法、装置、电子设备及存储介质
CN115877420B (zh) 一种基于静止卫星的定位方法、系统、电子设备及介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant