CN112052589A - 基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置 - Google Patents
基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置 Download PDFInfo
- Publication number
- CN112052589A CN112052589A CN202010914015.4A CN202010914015A CN112052589A CN 112052589 A CN112052589 A CN 112052589A CN 202010914015 A CN202010914015 A CN 202010914015A CN 112052589 A CN112052589 A CN 112052589A
- Authority
- CN
- China
- Prior art keywords
- pixel
- day
- data
- snow
- processed
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 43
- 238000012545 processing Methods 0.000 claims abstract description 21
- 238000004458 analytical method Methods 0.000 claims abstract description 17
- 238000004364 calculation method Methods 0.000 claims abstract description 11
- 238000012360 testing method Methods 0.000 claims abstract description 10
- 230000015572 biosynthetic process Effects 0.000 claims description 8
- 238000003786 synthesis reaction Methods 0.000 claims description 8
- 238000001914 filtration Methods 0.000 claims description 7
- 238000012937 correction Methods 0.000 claims description 5
- 238000007689 inspection Methods 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 4
- 238000000611 regression analysis Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000001131 transforming effect Effects 0.000 claims description 4
- 238000012952 Resampling Methods 0.000 claims description 3
- 230000001419 dependent effect Effects 0.000 claims description 3
- 230000002194 synthesizing effect Effects 0.000 claims description 3
- 230000002123 temporal effect Effects 0.000 claims 1
- 238000012544 monitoring process Methods 0.000 abstract description 6
- 230000000875 corresponding effect Effects 0.000 description 14
- 230000002354 daily effect Effects 0.000 description 6
- 230000007547 defect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000007635 classification algorithm Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 230000003203 everyday effect Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Astronomy & Astrophysics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Multimedia (AREA)
- Remote Sensing (AREA)
- Evolutionary Biology (AREA)
- Geometry (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Evolutionary Computation (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Image Processing (AREA)
Abstract
本申请涉及基于MODIS逐日积雪产品的积雪覆盖比例的估算方法和装置，属于积雪遥感产品处理技术领域，本申请的方法包括，对获取的积雪覆盖比例产品数据进行像元值重编码处理，确定待处理像元；之后依次基于时间关系的填补分析、局地雪线分析对待处理像元进行积雪覆盖度填补；再之后针对剩余的待处理像元分别构建时间‑空间局部加权回归模型并对构建的模型进行显著性检验，检验通过时，应用该模型对相应的待处理像元进行积雪覆盖度填补，否则，保持该待处理像元的编码不变；之后再采用临近日积雪覆盖度迭代计算进行积雪覆盖度填补，以得到逐日积雪覆盖比例的估算结果。本申请可显著提高积雪监测的精度。
Description
技术领域
本申请属于积雪遥感产品处理技术领域，具体涉及一种基于MODIS逐日积雪产品的积雪覆盖比例的估算方法和装置。
背景技术
积雪是地球表面变化最活跃的要素之一，对地气系统的辐射平衡、能量交换、水分循环具有重要的影响。准确、快速的监测积雪的变化过程尤为重要。卫星遥感技术具有宏观、多波段、多尺度、多时相等优势，自20世纪60年代以来发展迅速，在积雪快速周期变化监测中发挥着重要的作用。Terra和Aqua卫星携带的MODIS传感器兼顾时间分辨率、空间分辨率和光谱分辨率，实现了对地表主要地物的逐日周期性数据获取和监测。基于MODIS数据，已经开发出了具有较高精度的积雪覆盖比例(Fractional Snow Cover)产品。但是，受云层覆盖的影响，该产品中存在大量的云像元，极大地限制了积雪产品的应用。因此，如何结合时间、空间上下文的特征获取云覆盖下的积雪覆盖比例是开展积雪监测的基础。
为了消除云覆盖对积雪产品的影响，近年来国内外学者开展了相关研究，并发表了部分论文。这些研究可大致分为三类：第一类是通过多日合成的方式消除或降低云的影响。第二类是结合微波遥感的数据，利用其能够穿云透雾的特性来填补云覆盖区积雪信息。第三类是利用积雪覆盖在时间和空间上相关的特性，通过制定一定的判别规则来区分“雪像元”和“非雪像元”。上述几种方法得到的去云积雪产品，多为二值积雪产品，由于二值分类算法的局限性，在积雪混合像元区域，会明显高估或低估地表实际积雪覆盖情况。多日合成的去云算法会降低积雪合成产品的时间分辨率，而结合微波遥感数据的方法得到的结果会存在明显的分块效应，这是由于现有的逐日覆盖的微波遥感产品分辨率显著低于MODIS光学遥感产品而造成的。
已有关于逐日积雪覆盖比例估算的专利有“一种基于MODIS逐日积雪覆盖图像的去云方法及装置(CN106127717B)”和“一种对逐日积雪覆盖率图像处理的方法(CN108053440B)”两项。专利CN106127717B实现了对MODIS逐日积雪影像的去云处理，但仅是通过二值判断的方式将像元分为“雪”与“非雪”两类，不能获得逐日积雪覆盖比例数据。专利CN106127717B实现了对逐日积雪覆盖度的计算，在通过一定的规则判定为积雪的像元，其积雪覆盖率表达为“时间相邻像元所对应的积雪覆盖率的平均值(时间平均)”或“空间相邻类型为积雪的像元所对应的积雪覆盖率的平均值(空间平均)”。在积雪覆盖率计算中，未将云像元一定时间和空间范围内的积雪覆盖像元纳入到同一个积雪覆盖比例估算模型，未引入地形因素对积雪覆盖比例进行估算，所得结果没有实现对全部云覆盖像元的去除。
上述内容仅用于辅助理解本发明的技术方案，并不代表承认上述内容是现有技术。
发明内容
为至少在一定程度上克服相关技术中存在的问题，本申请提供一种基于MODIS逐日积雪产品的积雪覆盖比例的估算方法和装置，有助于提高积雪覆盖度的估算精度，并实现了对估算结果精度的检验、估算方法不确定性的评估。
为实现以上目的，本申请采用如下技术方案：
第一方面，
本申请提供一种基于MODIS逐日积雪产品的积雪覆盖比例的估算方法，该方法包括：
步骤1、对获取的MODIS积雪覆盖比例产品数据进行像元值重编码处理，确定待处理像元，得到第一中间数据；
步骤2、针对所述第一中间数据，基于时间关系的填补分析对待处理像元进行积雪覆盖度填补，得到第二中间数据；
步骤3、针对所述第二中间数据，基于局地雪线分析对待处理像元进行积雪覆盖度填补，得到第三中间数据；
步骤4、针对所述第三中间数据中的待处理像元，分别构建时间-空间局部加权回归模型并对构建的模型进行显著性检验，
当检验通过时，应用该模型对相应的待处理像元进行积雪覆盖度填补，否则，保持该待处理像元的编码不变，
其中，所述时间-空间局部加权回归模型的因变量为像元的积雪覆盖度，自变量为对应像元的高程、坡度、坡向和时间变量；
步骤5、针对步骤4完成后得到数据，基于临近日积雪覆盖度迭代计算进行积雪覆盖度填补，以完全消除云覆盖像元得到逐日积雪覆盖比例的估算结果；
步骤6、对步骤5中得到的估算结果进行精度检验，以确定估算结果的准确性和方法的可靠性。
可选地，所述步骤1具体为：
下载MODIS积雪覆盖比例产品数据，对所述MODIS积雪覆盖比例数据进行镶嵌、裁切预处理，并对预处理后的数据按时间顺序进行波段合成；
根据第一预定规则，将波段合成后的数据中各像元的数据编码值进行编码变换，以得到第一中间数据；
其中，所述第一预设规则包括，
将指示内陆、海洋的数据编码值变换为指示无积雪覆盖的数据编码值，
将指示夜晚、云、数据缺失、无意义数据的数据编码值变换为指示待处理的数据编码值。
可选地，所述步骤2中基于时间关系的填补分析对待处理像元进行积雪覆盖度填补，包括，
根据第二预定规则，基于对所述第一中间数据中同一位置像元前后多日的数据进行合成，以实现对第一中间数据中的待处理像元的积雪覆盖度填补；
其中，所述第二预设规则包括，
若某日某像元为待处理像元，当该像元前一日的积雪覆盖度、后一日的积雪覆盖度均为0时，则判定该日该像元的积雪覆盖度为0并进行填补，
若某日某像元为待处理像元，当该像元前一日的积雪覆盖度为0，该像元后一日的数据编码值指示为待处理，且该像元后二日的积雪覆盖度为0，则判定该像元当日和后一日的积雪覆盖度为0并进行填补，
若某日某像元为待处理像元，当该像元前一日、后一日的数据编码均指示为待处理，且该像元之前及之后第二日的积雪覆盖度均为0，则判定该像元该日及前后两日的积雪覆盖度为0并进行填补，
若某日某像元为待处理像元，当该像元前一日、后一日均存在积雪覆盖度值，且两日的积雪覆盖度值相差不超过20％，则判定该像元该日的积雪覆盖度为前后两日积雪覆盖度的平均值并进行填补。
可选地，所述步骤3中基于局地雪线分析对待处理像元进行积雪覆盖度填补包括，针对所述第二中间数据中的每一天的数据，分别进行如下处理：
建立指定大小的空间滤波窗口，遍历该日数据中所有空间区域并比较空间滤波窗口内每个像元的高程值大小，若某个待处理像元的高程值小于该窗口内非积雪像元高程值的平均值，判定该日该像元的积雪覆盖度为0并进行填补。
可选地，所述步骤4中针对所述第三中间数据中的待处理像元，分别构建时间-空间局部多元加权回归模型，包括针对所述第三中间数据中每一天的数据，分别进行如下处理：
对该日数据中每一待处理像元分别进行如下处理，
以该待处理像元为中心，选取指定的时间范围和空间范围，形成时间-空间立方体；
基于该时间-空间立方体内像元的积雪覆盖度数据、像元对应地理位置的高程数据、坡度数据和坡向数据，以及时间变量进行回归分析，构建该待处理像元所对应的时间-空间局部加权回归模型。
可选地，所述坡向数据为基于SRTM DEM数据进行修正后得到的修正坡向数据；所述修正过程，具体为：
基于下载的SRTM DEM数据完成数字地形高程模型的镶嵌与裁切，将数字地形高程模型分辨率重采样为与所述MODIS数据一致；
从重采样后的数字地形高程模型中提取坡向参数，使其表示为与北方向的夹角并进行更新。
可选地，所述步骤4中采用F检验对构建的模型进行显著性检验。
可选地，所述步骤5中基于临近日积雪覆盖度迭代计算进行积雪覆盖度填补，具体为，
尝试将待处理像元的数据编码值赋值为其前一日和后一日的积雪覆盖度的平均值，若其前一日或后一日为待处理像元，
则扩大时间窗口直至其前后日不为云，确定该日相对于时间窗口两端的时间距离，并根据所述时间距离将该日待处理像元的数据编码值赋值为时间窗口两端日的积雪覆盖度的加权平均值。
可选地，所述MODIS积雪覆盖比例产品数据来源配置为上午星，所述步骤6中，采用来源为下午星的相应MODIS积雪覆盖比例产品数据，对估算结果进行精度检验。
第二方面，
本申请提供一种估算装置，包括：
存储器，其上存储有可执行程序；
处理器，用于执行所述存储器中的所述可执行程序，以实现上述所述方法的步骤。
本申请采用以上技术方案，至少具备以下有益效果：
本申请综合多个去云步骤，在保持积雪产品时间和空间分辨率不变的同时，充分考虑积雪时空集中分布的特性，将云像元一定时间和空间范围内的积雪覆盖像元纳入到同一个积雪覆盖比例估算模型，并结合数字高程模型、坡度、坡向等地形信息，实现了对云覆盖像元下的积雪覆盖比例的估算，可有效地弥补二值积雪产品的不足，显著提高积雪监测的精度。本申请通过临近日迭代计算能够消除全部云覆盖像元，获得逐日无云积雪覆盖比例产品。在此基础上，利用空间尺度一致、精度相对可靠的MODIS下午星数据对去云效果进行检验，实现对去云结果精度、去云方法可靠性的评价。
本发明的其他优点、目标，和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书，权利要求书，以及附图中所特别指出的结构来实现和获得。
附图说明
附图用来提供对本申请的技术方案或现有技术的进一步理解，并且构成说明书的一部分。其中，表达本申请实施例的附图与本申请的实施例一起用于解释本申请的技术方案，但并不构成对本申请技术方案的限制。
图1为本申请一个实施例提供的基于MODIS逐日积雪产品的积雪覆盖比例的估算方法的流程示意图；
图2为本申请一个实施例提供的估算装置的结构示意图。
具体实施方式
为使本申请的目的、技术方案和优点更加清楚，下面将对本申请的技术方案进行详细的描述。显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本申请所保护的范围。
针对现有技术中缺陷，本申请的技术方案构建局部时间空间多元统计回归模型，充分考虑积雪分布时空相关特性，来实现综合时间、空间、地形等多要素的逐日积雪覆盖度估算。
图1所示，为本申请一实施例中基于MODIS逐日积雪产品的积雪覆盖比例的估算方法的流程示意图。该实施例中，估算方法包括如下步骤：
步骤1，对获取的MODIS积雪覆盖比例产品数据进行像元值重编码处理，确定待处理像元，得到第一中间数据。
具体的，步骤1中，首先从美国冰雪数据中心(NSIDC)下载MODIS积雪覆盖比例产品数据，之后对MODIS积雪覆盖比例产品数据进行镶嵌、裁切预处理(该实施例中，利用R语言完成数据的下载，以及之后的对其的镶嵌、裁切等预处理)，并对预处理后的数据按时间顺序进行波段合成；
再根据第一预设规则，将波段合成后的数据中各像元的数据编码值进行编码变换(或者说重分类)，以得到第一中间数据；其中，第一预设规则包括，将指示内陆、海洋的数据编码值变换为指示无积雪覆盖的数据编码值，将指示夜晚、云、数据缺失、无意义数据的数据编码值变换为指示待处理的数据编码值，需要说明的是，本申请中像元的编码值为待处理的值时，认为其对应的像元为被云覆盖的像元，该类像元是后续需要进行处理消除的目标。
具体的，在该实施例中，MODIS积雪覆盖比例产品数据的空间分辨率为500米，时间跨度为4年，数据编码值规则为：0-100为积雪覆盖比例，211为夜晚，237为内陆水体，239为海洋，250为云，200、201、254、255为数据缺失或无意义数据。
基于此，该实施例中，第一预设规则所体现的数据编码值变换规则具体如下：
数据编码值为237和239的像元、其数据编码值重编码为0；数据编码值编为211、250、200、201、254和255的像元、其数据编码值重编码为150，换言之，在该实施例中，数据编码值为150的像元为待处理像元(或者称为云覆盖像元)。
步骤1之后，进行步骤2，针对第一中间数据，基于时间关系的填补分析对待处理像元进行积雪覆盖度填补，得到第二中间数据；
基于时间关系的填补分析的原理依据为，降雪具有时间相关的特点，积雪会在地面上存留一段时间，基于此可对待处理像元进行积雪覆盖度填补，该步骤中，基于时间关系的填补分析处理的过程，包括：
根据第二预定规则，基于对第一中间数据中同一位置像元前后多日的数据进行合成，以实现对第一中间数据中的待处理像元的积雪覆盖度填补。其中，第二预设规则包括：
若某日某像元为待处理像元，当该像元前一日的积雪覆盖度、后一日的积雪覆盖度均为0时，则判定该日该像元的积雪覆盖度为0并进行填补，即被云覆盖的像元，前后两日均无积雪，则判定该日该像元无积雪：
若某日某像元为待处理像元，当该像元前一日的积雪覆盖度为0，该像元后一日的数据编码值指示为待处理，且该像元后二日的积雪覆盖度为0，则判定该像元当日和后一日的积雪覆盖度为0并进行填补，即连续两日被云覆盖的像元，若前一日和后第二日积雪覆盖度为0，则判定上述连续两日被云覆盖的像元积雪覆盖度为0；
若某日某像元为待处理像元，当该像元前一日、后一日的数据编码均指示为待处理，且该像元之前及之后第二日的积雪覆盖度均为0，则判定该像元该日及前后两日的积雪覆盖度为0并进行填补，即连续三日被云覆盖的像元，若前第二日和后第二日积雪覆盖度为0，则认为上述连续三日被云覆盖的像元积雪覆盖度为0；
若某日某像元为待处理像元，当该像元前一日、后一日均存在积雪覆盖度值，且两日的积雪覆盖度值相差不超过20％，则判定该像元该日的积雪覆盖度为前后两日积雪覆盖度的平均值并进行填补。
此外针对基于时间关系的填补分析处理，还需要说明的是，步骤2中完整的分析处理过程需依时间顺序针对每天的数据依次进行处理更新，直至全部数据分析处理完，且容易理解的是，为实现第二预定规则，需要在处理的时间段两端外额外获得前后几天的数据，如本实施例中是对四年的数据进行处理，则需额外获取该四年时间段前后各两天的数据。
步骤2之后，为进一步消除云覆盖像元，继续进行步骤3，针对第二中间数据，基于局地雪线分析对待处理像元进行积雪覆盖度填补，得到第三中间数据。
局部雪线分析处理的原理依据为，积雪分布具有空间相关规律，在局部区域相近地理单元内，低海拔地区有积雪覆盖，则高海拔地区有积雪覆盖的概率极大，反之亦然，小区域内高海拔地区无积雪覆盖则低海拔地区无积雪覆盖的概率极大，进而可以根据积雪在空间分布上与高程相关的特征来将标记为云的像元进行判定处理。
具体的，步骤3中基于局地雪线分析进行积雪覆盖度填补，包括针对第二中间数据中的每一天的数据，分别进行如下处理：
建立指定大小的空间滤波窗口，遍历该日数据中所有空间区域并比较空间滤波窗口内每个像元的高程值大小，若某个待处理像元的高程值小于该窗口内非积雪像元高程值的平均值时，判定该日该像元的积雪覆盖度为0并进行填补。
该实施例中，空间滤波窗口的大小为11像元*11像元，对应于实际区域为5.5km*5.5km。
通过上述步骤2的基于时间关系的填补分析处理、以及步骤3的局地雪线分析处理可以消除一部分云像元，但上述步骤对时间窗口长度内，位置没有发生显著变化的大面积云无法消除。
本申请针对于此，如图1所示，在步骤3之后进行步骤4，针对第三中间数据中的待处理像元(步骤3后仍存在的)，分别构建时间-空间局部加权回归模型并对构建的模型进行显著性检验，
当检验通过时，应用该模型对相应的待处理像元进行积雪覆盖度填补，否则，保持该待处理像元的编码不变(即其数据编码值仍指示为待处理)，这里时间-空间局部加权回归模型的因变量为像元的积雪覆盖度，自变量为对应像元的高程、坡度、坡向和时间变量。
具体的，步骤4中，针对第三中间数据中的待处理像元，分别构建时间-空间局部多元加权回归模型，包括针对第三中间数据中每一天的数据，分别进行如下处理：
对该日数据中每一待处理像元分别进行如下处理，
以该待处理像元为中心，选取指定的时间范围和空间范围，形成时间-空间立方体；
基于该时间-空间立方体内像元的积雪覆盖度数据、像元对应地理位置的高程数据、坡度数据和坡向数据，以及时间变量进行回归分析，来构建该待处理像元所对应的时间-空间局部加权回归模型。
举例而言，该实施例中，以该待处理像元(云覆盖像元)为中心，选取11像元*11像元*5天(该日前后2天)的时间-空间立方体数据来进行回归分析、并构建时间-空间局部加权回归模型，时间-空间局部加权回归具有如下表达形式：
表达式(1)中，
表示第j个时空立方体中第i点的积雪覆盖度，ELEVi j、
M_ASPECTi j、Ti j分别表示第j个时空立方体中第i点的高程、坡度、坡向和时间次序，β0、β1、β2、β3、β4为局域回归系数，
为误差项。
需要说明的是，上述构建过程中，坡度数据、坡向数据均基于SRTM DEM数据获得，其中坡向数据为基于SRTM DEM数据进行修正后得到的修正坡向数据，该修正过程，具体如下：
基于下载的SRTM DEM数据完成数字地形高程模型的镶嵌与裁切，将数字地形高程模型分辨率重采样为与所述MODIS数据一致；
从重采样后的数字地形高程模型中提取坡向参数，使其表示为与北方向的夹角并进行更新，即保持原始坡向参数中与北方向夹角小于等于180°的参数值不变，将大于180°的参数值，以如下表示式进行重新定义：
M_Aspectx,y＝360°-Aspectx,y (2)
表达式(2)中M_Aspectx,y为空间位置(x,y)处修正后的坡向值，Aspectx,y为修正前的坡向值。
采用如此方式得到的修正坡向数据，综合考虑了太阳光线角度因素，以及坡向对积雪覆盖分布的影响要求，修正后的坡向能够更加准确的反映太阳照射对积雪覆盖分布的影响，物理意义更加明确，据此进行模型构建可提高后续估算结果的精度。
在该实施例中，步骤4中采用F检验对构建的模型进行显著性检验，其过程为，构建如下形式的统计量F：
上述表达式(3)式中ESS为回归平方和，RSS为剩余平方和，n为参与回归的积雪覆盖像元个数，k为模型中自变量的个数，这里取值为4；
基于表达式(3)，给定显著性水平α，查F分布表可得到临界值Fα(k,n-k-1)k，以此来判定模型的线性关系是否显著成立。
如图1所示，步骤4后进行步骤5，针对步骤4完成后得到数据，基于临近日积雪覆盖度迭代计算进行积雪覆盖度填补，以完全消除云覆盖像元得到逐日积雪覆盖比例的估算结果。
步骤4完成后得到数据中仍存在未能去云的待处理像元，在步骤5中，对这些待处理像元采用临近日积雪覆盖度迭代计算进行积雪覆盖度填补，临近日积雪覆盖度迭代计算的具体方式为，
尝试将待处理像元的数据编码值赋值为其前一日和后一日的积雪覆盖度的平均值，若其前一日或后一日为待处理像元，则扩大时间窗口直至其前后日不为云，确定该日相对于时间窗口两端的时间距离，并根据时间距离将该日待处理像元的数据编码值赋值为时间窗口两端日的积雪覆盖度的加权平均值。
举例而言，若某日某待处理像元的前一日和后一日的数据编码值分别为40、50，则该待处理像元的数据编码值赋值为45，
而若该日该待处理像元的前一日为待处理像元，后一日的数据编码值为50，则向前扩大时间窗口为前两日，并再次判断前第二日是否均指示为积雪覆盖度……如此直至时间窗口前后两端日均指示为积雪覆盖度，假设扩大后前第二日的积雪覆盖度为20，则将该待处理像元的数据编码值赋值为40。
如此，经过上述步骤1至5后，待处理像元(云覆盖像元)被完全消除，得出了估算结果，该估算结果是与原始积雪覆盖比例产品同分辨率的逐日无云积雪覆盖比例产品，如该实施例中，得到的是500米分辨率的逐日无云积雪覆盖比例产品。
如图1所示，步骤5后还包括步骤6，对得到的估算结果进行精度检验，以确定估算结果的准确性和方法的可靠性。
由于云具有移动的特点，MODIS积雪覆盖比例产品中，上、下午星获取的积雪产品MOD10A1及MYD10A1云覆盖情况略有差异，因此可以仅利用MOD10A1或MYD10A1产品进行处理，待上述去云处理(消除待处理像元)完成后，采用另一产品进行精度的验证。采用这种方式可保证每日积雪覆盖比例影像中都存在大量精度检测点，且因上、下午星使用的遥感传感器相同，可避免由于空间分辨率、波段设置等方面的差异对评价结果造成影响。
已有研究表明，下午星因星上仪器设备故障与上午星相比精度略低，因此在该实施例中，利用上午星获取的MOD10A1产品进行去云处理，而将下午星数据获取的MYD10A1积雪覆盖比例作为“真值”，对生成的无云积雪覆盖比例产品的精度进行验证。具体在该实施例中，MODIS积雪覆盖比例产品数据来源配置为上午星，步骤6中采用来源为下午星的相应MODIS积雪覆盖比例产品数据，对估算结果进行精度检验。
本申请综合多个去云步骤，在保持积雪产品时间和空间分辨率不变的同时，充分考虑积雪时空集中分布的特性，将云像元一定时间和空间范围内的积雪覆盖像元纳入到同一个积雪覆盖比例估算模型，并结合数字高程模型、坡度、坡向等地形信息，实现了对云覆盖像元下的积雪覆盖比例的估算，可有效地弥补二值积雪产品的不足，显著提高积雪监测的精度。本申请通过临近日迭代计算能够消除全部云覆盖像元，获得逐日无云积雪覆盖比例产品。在此基础上，利用空间尺度一致、精度相对可靠的MODIS下午星数据对去云效果进行检验，实现对去云结果精度、去云方法可靠性的评价。
图2为本申请一个实施例提供的估算装置的结构示意图，如图2所示，该估算设备200包括：
存储器201，其上存储有可执行程序；
处理器202，用于执行存储器201中的可执行程序，以实现上述方法的步骤。
关于上述实施例中的估算装置200，其处理器202执行存储器201中的程序的具体方式已经在有关该方法的实施例中进行了详细描述，此处将不做详细阐述说明。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人员在本发明所揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。
Claims (10)
1.一种基于MODIS逐日积雪产品的积雪覆盖比例的估算方法，其特征在于，包括：
步骤1、对获取的MODIS积雪覆盖比例产品数据进行像元值重编码处理，确定待处理像元，得到第一中间数据；
步骤2、针对所述第一中间数据，基于时间关系的填补分析对待处理像元进行积雪覆盖度填补，得到第二中间数据；
步骤3、针对所述第二中间数据，基于局地雪线分析对待处理像元进行积雪覆盖度填补，得到第三中间数据；
步骤4、针对所述第三中间数据中的待处理像元，分别构建时间-空间局部加权回归模型并对构建的模型进行显著性检验，
当检验通过时，应用该模型对相应的待处理像元进行积雪覆盖度填补，否则，保持该待处理像元的编码不变，
其中，所述时间-空间局部加权回归模型的因变量为像元的积雪覆盖度，自变量为对应像元的高程、坡度、坡向和时间变量；
步骤5、针对步骤4完成后得到数据，基于临近日积雪覆盖度迭代计算进行积雪覆盖度填补，以完全消除云覆盖像元，得到逐日积雪覆盖比例的估算结果；
步骤6、对步骤5中得到的估算结果进行精度检验，以确定估算结果的准确性和方法的可靠性。
2.根据权利要求1所述的估算方法，其特征在于，所述步骤1具体为：
下载MODIS积雪覆盖比例产品数据，对所述MODIS积雪覆盖比例数据进行镶嵌、裁切预处理，并对预处理后的数据按时间顺序进行波段合成；
根据第一预定规则，将波段合成后的数据中各像元的数据编码值进行编码变换，以得到第一中间数据；
其中，所述第一预设规则包括，
将指示内陆、海洋的数据编码值变换为指示无积雪覆盖的数据编码值，
将指示夜晚、云、数据缺失、无意义数据的数据编码值变换为指示待处理的数据编码值。
3.根据权利要求1所述的估算方法，其特征在于，所述步骤2中基于时间关系的填补分析对待处理像元进行积雪覆盖度填补，包括，
根据第二预定规则，基于对所述第一中间数据中同一位置像元前后多日的数据进行合成，以实现对第一中间数据中的待处理像元的积雪覆盖度填补；
其中，所述第二预设规则包括，
若某日某像元为待处理像元，当该像元前一日的积雪覆盖度、后一日的积雪覆盖度均为0时，则判定该日该像元的积雪覆盖度为0并进行填补，
若某日某像元为待处理像元，当该像元前一日的积雪覆盖度为0，该像元后一日的数据编码值指示为待处理，且该像元后二日的积雪覆盖度为0，则判定该像元当日和后一日的积雪覆盖度为0并进行填补，
若某日某像元为待处理像元，当该像元前一日、后一日的数据编码均指示为待处理，且该像元之前及之后第二日的积雪覆盖度均为0，则判定该像元该日及前后两日的积雪覆盖度为0并进行填补，
若某日某像元为待处理像元，当该像元前一日、后一日均存在积雪覆盖度值，且两日的积雪覆盖度值相差不超过20％，则判定该像元该日的积雪覆盖度为前后两日积雪覆盖度的平均值并进行填补。
4.根据权利要求1所述的估算方法，其特征在于，所述步骤3中基于局地雪线分析对待处理像元进行积雪覆盖度填补包括，针对所述第二中间数据中的每一天的数据，分别进行如下处理：
建立指定大小的空间滤波窗口，遍历该日数据中所有空间区域并比较空间滤波窗口内每个像元的高程值大小，若某个待处理像元的高程值小于该窗口内非积雪像元高程值的平均值，判定该日该像元的积雪覆盖度为0并进行填补。
5.根据权利要求1所述的估算方法，其特征在于，所述步骤4中针对所述第三中间数据中的待处理像元，分别构建时间-空间局部多元加权回归模型，包括针对所述第三中间数据中每一天的数据，分别进行如下处理：
对该日数据中每一待处理像元分别进行如下处理，
以该待处理像元为中心，选取指定的时间范围和空间范围，形成时间-空间立方体；
基于该时间-空间立方体内像元的积雪覆盖度数据、像元对应地理位置的高程数据、坡度数据和坡向数据，以及时间变量进行回归分析，构建该待处理像元所对应的时间-空间局部加权回归模型。
6.根据权利要求5所述的估算方法，其特征在于，
所述坡向数据为基于SRTM DEM数据进行修正后得到的修正坡向数据；所述修正过程，具体为：
基于下载的SRTM DEM数据完成数字地形高程模型的镶嵌与裁切，将数字地形高程模型分辨率重采样为与所述MODIS数据一致；
从重采样后的数字地形高程模型中提取坡向参数，使其表示为与北方向的夹角并进行更新。
7.根据权利要求1所述的估算方法，其特征在于，所述步骤4中采用F检验对构建的模型进行显著性检验。
8.根据权利要求1所述的估算方法，其特征在于，所述步骤5中基于临近日积雪覆盖度迭代计算进行积雪覆盖度填补，具体为，
尝试将待处理像元的数据编码值赋值为其前一日和后一日的积雪覆盖度的平均值，若其前一日或后一日为待处理像元，
则扩大时间窗口直至其前后日不为云，确定该日相对于时间窗口两端的时间距离，并根据所述时间距离将该日待处理像元的数据编码值赋值为时间窗口两端日的积雪覆盖度的加权平均值。
9.根据权利要求1所述的估算方法，其特征在于，所述MODIS积雪覆盖比例产品数据来源配置为上午星，所述步骤6中，采用来源为下午星的相应MODIS积雪覆盖比例产品数据，对估算结果进行精度检验。
10.一种估算装置，其特征在于，包括：
存储器，其上存储有可执行程序；
处理器，用于执行所述存储器中的所述可执行程序，以实现权利要求1-9中任一项所述方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010914015.4A CN112052589A (zh) | 2020-09-03 | 2020-09-03 | 基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010914015.4A CN112052589A (zh) | 2020-09-03 | 2020-09-03 | 基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112052589A true CN112052589A (zh) | 2020-12-08 |
Family
ID=73608407
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010914015.4A Pending CN112052589A (zh) | 2020-09-03 | 2020-09-03 | 基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112052589A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114494864A (zh) * | 2022-01-17 | 2022-05-13 | 中国科学院地理科学与资源研究所 | 基于遥感数据的积雪物候信息的提取方法 |
| CN117635820A (zh) * | 2023-10-24 | 2024-03-01 | 中国科学院西北生态环境资源研究院 | 一种modis积雪面积产品缺失信息时空重建方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005003563A (ja) * | 2003-06-13 | 2005-01-06 | Niigata Denki Kk | 積雪情報提供システム |
| CN103984862A (zh) * | 2014-05-15 | 2014-08-13 | 中国科学院遥感与数字地球研究所 | 一种多元遥感信息协同的积雪参数反演方法 |
| CN106127717A (zh) * | 2016-06-15 | 2016-11-16 | 中国科学院遥感与数字地球研究所 | 一种基于modis逐日积雪覆盖图像的去云方法及装置 |
| CN108053440A (zh) * | 2017-12-27 | 2018-05-18 | 中国科学院遥感与数字地球研究所 | 一种对逐日积雪覆盖率图像处理的方法 |
-
2020
- 2020-09-03 CN CN202010914015.4A patent/CN112052589A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005003563A (ja) * | 2003-06-13 | 2005-01-06 | Niigata Denki Kk | 積雪情報提供システム |
| CN103984862A (zh) * | 2014-05-15 | 2014-08-13 | 中国科学院遥感与数字地球研究所 | 一种多元遥感信息协同的积雪参数反演方法 |
| CN106127717A (zh) * | 2016-06-15 | 2016-11-16 | 中国科学院遥感与数字地球研究所 | 一种基于modis逐日积雪覆盖图像的去云方法及装置 |
| CN108053440A (zh) * | 2017-12-27 | 2018-05-18 | 中国科学院遥感与数字地球研究所 | 一种对逐日积雪覆盖率图像处理的方法 |
Non-Patent Citations (4)
| Title |
|---|
| 侯小刚;郑照军;李帅;陈雪华;崔宇;: "近15年新疆逐日无云积雪覆盖产品生成及精度验证", 国土资源遥感, no. 02, pages 214 - 222 * |
| 唐志光;王建;李弘毅;彦立利;梁继;: "青藏高原MODIS积雪面积比例产品的精度验证与去云研究", 遥感技术与应用, no. 03 * |
| 王子龙;胡石涛;付强;姜秋香;印玉明;: "基于多源数据的松嫩平原黑土区亚像元雪盖率算法研究", 农业机械学报, no. 02 * |
| 邱玉宝;张欢;除多;张雪成;于小淇;郑照军;: "基于MODIS的青藏高原逐日无云积雪产品算法", 冰川冻土, no. 03 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114494864A (zh) * | 2022-01-17 | 2022-05-13 | 中国科学院地理科学与资源研究所 | 基于遥感数据的积雪物候信息的提取方法 |
| CN117635820A (zh) * | 2023-10-24 | 2024-03-01 | 中国科学院西北生态环境资源研究院 | 一种modis积雪面积产品缺失信息时空重建方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111078678B (zh) | 一种基于多源信息融合与降尺度的卫星降水数据校正方法 | |
| Verger et al. | The CACAO method for smoothing, gap filling, and characterizing seasonal anomalies in satellite time series | |
| Hame et al. | Improved mapping of tropical forests with optical and SAR imagery, Part II: Above ground biomass estimation | |
| Davaze et al. | Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data | |
| CN109308688B (zh) | 一种可见光和近红外波段厚云及阴影去除方法 | |
| CN114627087B (zh) | 一种多时相卫星遥感图像的地物变化自动检测方法及系统 | |
| CN112052589A (zh) | 基于modis逐日积雪产品的积雪覆盖比例的估算方法和装置 | |
| CN114092835B (zh) | 基于不同时空分辨率的归一化植被指数数据时空融合方法 | |
| CN114819737B (zh) | 公路路域植被的碳储量估算方法、系统及存储介质 | |
| Hammad et al. | Land cover change investigation in the southern Syrian coastal basins during the past 30-years using Landsat remote sensing data | |
| CN111798132B (zh) | 基于多源时序遥感深度协同下的耕地动态监测方法及系统 | |
| CN111982822A (zh) | 一种长时间序列高精度植被指数改进算法 | |
| KR101089220B1 (ko) | 시공간 연속성을 이용한 식생 지수 보정방법 | |
| Volpe et al. | Operational interpolated ocean colour product in the mediterranean sea | |
| Wei et al. | A comparative assessment of multisensor data merging and fusion algorithms for high-resolution surface reflectance data | |
| Verger et al. | GEOV2/VGT: Near real time estimation of global biophysical variables from VEGETATION-P data | |
| CN115079172A (zh) | 一种MTInSAR滑坡监测方法、设备及存储介质 | |
| Sobrino et al. | Near real-time estimation of Sea and Land surface temperature for MSG SEVIRI sensors | |
| CN107576399B (zh) | 面向modis林火探测的亮温预测方法和系统 | |
| CN115222837A (zh) | 真彩云图生成方法、装置、电子设备及存储介质 | |
| CN115452167A (zh) | 基于不变像元的卫星遥感器交叉定标方法和装置 | |
| CN113936009A (zh) | 一种气象卫星洪涝监测的云阴影去除方法、装置及设备 | |
| Bannari et al. | Potential of WorldView-3 for soil salinity modeling and mapping in an arid environment | |
| Amos et al. | A continuous vertically resolved ozone dataset from the fusion of chemistry climate models with observations using a Bayesian neural network | |
| Moisture | Algorithm theoretical baseline document (ATBD) D2. 1 Version 04.2 merging active and passive soil moisture retrievals |
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 |