CN110472525B - 一种时间序列遥感植被指数的噪声检测方法 - Google Patents

Abstract

一种时间序列遥感植被指数(Time‑Series Vegetation Index，TSVI)的噪声检测方法。首先使用单位根检验将各个像素的观测值分为平稳序列或者非平稳序列；对于非平稳序列，利用一定的数学模型对离散的TSVI进行建模，再计算实际观测值与模型预测值之间的差异，记为偏差。由于偏差消除了季节性成分，从而将非平稳序列转换为平稳序列。对于平稳序列或者偏差数据采用观测值分布在均值附近一定范围内的假设，进行噪声检测；再对去除噪声保留之后的观测值，迭代进行拟合和噪声检测，直到达到最大迭代次数或在某次迭代不再有噪声检出。然后将时间序列转换回图像空间获得噪声掩膜并优化。本发明可以获得准确的噪声掩膜，并提高地表相关应用的可靠性。

Description

一种时间序列遥感植被指数的噪声检测方法

技术领域

本发明公开的技术涉及遥感影像处理，更具体地说，是一种从时间序列遥感植被指数中进行噪声检测的技术。

背景技术

与单期或多期影像相比，时间序列遥感影像能够提供地表特征随时间变化的信息，为精确地追踪地表演变过程提供了可能。近年来，随着中等分辨率影像(如Landsat系列卫星获得的TM、ETM+、OLI影像，Sentinel卫星获取的MIS影像)的开放发布和计算技术的快速发展，利用所有可用影像追踪地表时间演变得到了更加广泛地应用。多期的、多光谱遥感影像使得时间序列呈现出一个高维的空间，给数据分析和信息提取带来较大的困难。植被指数能够在一定程度上抑制辐射噪声并突出植被等信息，而时间序列植被指数(Time-Series Vegetation Index,TSVI)能够描述植被等特征随时间的变化，并已经广泛地应用于森林监测、农作物识别等应用中。

遥感成像过程中，云、云阴影、传感器坏道等多种因素都可能导致传感器的观测值偏离地表的真实值，从而导致数据丢失或失真，即噪声。噪声会使得信息提取和影像分类过程出现错误，因而，识别噪声是正确使用数据的前提。具体到TSVI中，噪声导致观测值在较短时间间隔内出现“陡升陡降”的现象，从而掩盖地表的真实变化。同样地，标记受污染的观测值的过程是正确使用TSVI数据提取地表信息的前提。因而，准确地识别其中的噪声具有重要意义。

云和云阴影(以下简称云影)是一种最为常见的噪声，当前大多数的噪声检测方法也是面向云影的噪声检测，本发明主要针对云影引起的噪声，但也可以用于其他类型的噪声。在单期遥感影像上，根据先验知识集成方式，云影检测包括阈值分割类方法和监督分类方法两种主要类型。

阈值分割及其改进方法尝试使用一个或多个阈值将影像分为云影及其干净观测。典型的方法如在多光谱空间中通过一系列规则构建决策树，对影像进行层次分割以获得云影掩膜。然而，在遥感影像上，云影像素和干净像素的光谱特征存在严重的重叠现象，使其难以在原始空间中通过一个或多个阈值来分离。对此，将原始影像投影到一些变换空间中可以提高云/干净像素的可区分性并简化阈值选择过程。但由于云和干净像素之间缺乏内在的可分性，转换空间的阈值分割依然是误检和漏检之间的一种平衡，基于阈值的方法始终无法彻底解决云像素和干净像素的混淆问题。

监督分类方法通过云影样本训练分类器，然后通过训练好的分类器对待处理影像进行处理，获得云影掩膜，该方法不仅可以获取样本呈现的云特征，还可以使用支持向量机(Support Vector Machine,SVM)等分类器获得非线性决策空间。在监督分类方法中，特征设计和分类器选择是决定云检测性能的两个重要因素，颜色、纹理等多种特征及其SVM，浅层神经网络等机器学习方法或者其组合已经用于云影掩膜。近年来，新兴的基于深度学习的方法将特征提取与分类集成在一个结构中，利用标记的样本进行协同优化，获得了优于传统方法的性能。

尽管上述面向单期影像的云影掩膜方法取得了巨大的成功，促进了地表相关应用的发展，但其检测精度远非完美，尤其是上述方法难以实现对薄云和高反射率干净像素的区分。而在TSVI中，残留的薄云导致了时间相邻植被指数出现“陡升陡降”现象，这将掩藏实际的土地覆盖变化并降低时间序列应用的可信度。由于时间序列中影像数目的增多，通过人工解译方式获取薄云区域的劳动和财力消耗将更加巨大，因此，发展面向TSVI的噪声检测方法具有非常重要的意义。

传统噪声检测方法大多假设不同时间的观测值分布在平均值附近的一定范围内(如均值加减2倍标准差)变化，而将分布在此范围之外的观测值作为噪声。对于建筑物、道路和干净水体等伪不变特征，其植被指数不会随时间变化，而实际中获取的影像主要是由于云影等噪声因素引起，上述方法能够较好地适用于伪不变特征的噪声检测。然而，植被等多种土地覆盖类型具有随季节性变化的特征，如森林在春季生长，秋季落叶，其植被指数不仅随云影等噪声的变化，也存在本身特征随时间的变化，即特征是不平稳的，因而其难以利用上述方法进行噪声检测。

发明内容

为克服现有技术的上述缺点，本发明提供一种对TSVI进行噪声检测的新方法，其核心思想是从非平稳TSVI中分离出周期性变化成分获得偏差，从而将其转化为平稳序列。根据以上思想，本发明的实施方法为：

一种时间序列遥感植被指数的噪声检测方法，包括以下步骤：

步骤1，计算时序特征；

假设研究区有n期不同时间获取的影像I构成：

I＝<I₁,I₂,...,I_n> (1)

各期影像的获取日期T表示为：

T＝＜t₁,t₂,...,t_n＞ (2)

获取日期之间存在关系t_i-1≤t_i，i＝2,3,...,n，即t_i-1在t_i之前。

研究区有m个像素P表示为：

P＝{P¹,P²，...,P^m}^T (3)

像素

提取的TSVIX^j表示为：

植被指数可以使用归一化植被指数NDVI，也可以使用增强型植被指数EVI等其他植被指数，其特征在于该植被指数能够消除地表的辐射背景，并突出植被特征。

步骤2，对序列进行分类；

判断各个像素的TSVI的变化情况，将其划分到平稳序列集合S或者非平稳序列集合NS中，S和NS满足：

S∪NS＝P (6)

其中，平稳序列是指均值和标准差不随时间变化的序列，而非平稳序列是指均值和标准差会随时间变化的序列。平稳序列检验可以采用单位根检测方法，如AugmentedDickey-Fuller检验(ADFTest)等经典方法，其特征在于能够将各个像素划分到平稳序列集合或非平稳序列集合中。

步骤3，非平稳序列拟合；

对于非平稳序列集合NS中的像素，利用数学模型对离散的TSVI进行拟合。这里以高斯模型为例，说明其原理。

拟合过程中，自变量为T，应变量为X，拟合结果为

复合高斯模型表示为：

其中，a_k，μ_k，σ_k是待估算的参数，描述了高斯模型；α_K表示高斯成分的比例，即该峰值所占的比例；K表示高斯成分的数目，即曲线中包含几个峰谷组合；μ_k表示均值，即达到最大值时候的t值；σ_K决定了峰的陡峭程度。

本步骤中可以采用高斯模型及其改进的模型，如非对称高斯模型，也可以采用double logistic等其他模型，其特点在于该模型能够描述植被的周期性变化。

步骤4，计算偏差；

对于像素P^j∈NS，在t_i∈T时存在模型计算值

和实际观测值

其偏差

表示为：

其上标1表示第一次迭代结果。

表示了模型计算值

与实际观测值

之间的接近程度。

越大，表示二者差异越大，拟合效果越差；

越小，表示真实值和拟合值越接近，拟合效果越好。

像素P^j∈P在t_i∈T时的第一次偏差计算结果表示为：

由于偏差D^j,1消除了季节性成分，从而将非平稳序列转换为平稳序列；

步骤5，检测噪声；

对于转换后的偏差序列D^j，1和平稳序列S，分别采用观测值分布在均值附近一定范围内的假设，进行噪声检测。

对于P^j∈S，采用原始的序列X^j实施噪声检测；而对于P^j∈NS，采用偏差序列D^j,1，为了论述的方便，统一记为s^j：

然后按照以下方法判定噪声，

为掩膜数据，表示像素P^j∈P在t_i∈T是否为有效观测，1表示有效观测，而0表示噪声。

其中，μ(s^j)表示s^j的均值，σ(s^j)表示s^j的标准差，λ表示标准差的倍数。该算法的意义为，将距拟合曲线的距离分布在均值加减λ倍标准差的观测值作为有效观测，而分布在该范围之外的观测作为噪声。根据上述过程，可以获得像素P^j∈P的掩膜M^j。

步骤6，迭代处理；

对去除噪声保留之后的观测值，即不考虑

的像素，重新进行非平稳序列拟合(步骤3)、计算偏差(步骤4)和噪声检测(步骤5)，重复上述过程直到迭代达到最大次数或在某次迭代过程中不再有噪声检出。

步骤7，噪声掩膜；

将各个像素P^j∈P的掩膜数据转换回图像空间获得噪声掩膜M。

M＝＜M₁,M₂,...,M_n＞ (13)

其中，M_i∈M表示了影像I_i∈I上的各个观测值是否为有效观测值，在后续的应用中，需要将噪声进行剔除，以获得准确的应用效果。

步骤8，掩膜优化；

根据数学形态学方法对噪声进行优化，例如对孤立的椒盐噪声、大块噪声中的孤立有效数据点进行去除，获得平滑的掩膜结果。

本发明首先将序列分为平稳序列或非平稳序列，然后通过去除非平稳TSVI中的季节分量并获得偏差，偏差可以视为遵循高斯正态分布的平稳序列，并且大部分有效观测值分布在均值周围一定范围内，因此噪声可以被识别为远离平均值的观测值。

本发明的优点是：可以获得准确的噪声掩膜，并提高地表相关应用的可靠性。

附图说明

图1是研究区的时间序列影像示意图，研究区的土地覆盖类型包括建筑物、水体、森林和农田等，影像上的部分区域被厚云或薄云污染，遮挡了真实的地表状况。

图2A～图2D是建筑物、水体、森林和农田四种典型地物的归一化差异植被指数(NDVI)随时间变化的曲线，其中图2A为建筑物；图2B为水体；图2C为森林；图2D为农田。

图3是本发明的噪声检测流程图。

图4A～图4D是噪声检测迭代过程。图4A为原始时间序列特征的拟合结果；图4B为偏差及其第一次迭代的噪声检测结果；图4C为去除噪声之后的时间序列拟合结果；图4D为时序拟合过程及其噪声检测结果。

图5是图1中部分时间序列影像的噪声检测结果。

图6是实现本发明所叙述噪声检测的一个代表性系统。

具体实施方式

本发明的实施例提供了一种遥感TSVI进行噪声检测的方法，其具体实施方式如下：

图1展示了不同日期同一研究区的NDVI时间序列影像，颜色由灰变黑表示NDVI值逐渐减低。NDVI可以由相同或者不同的传感器，甚至在不同的平台上获取的影像进行计算得到，图1中所示的NDVI由Sentinel-2A和Sentinel-2B卫星上的多光谱仪成像仪(MSI)获得的影像计算得到。其他卫星或者传感器获取的影像也可用于构建时间序列。如果原始时间序列影像的空间分辨率不同，则需要将影像重新采样到相同的分辨率。

图1仅是研究区的部分时间序列影像，在数据分发网站一共有115景影像覆盖了部分或全部的研究区域，这里仅每个月选取1景影像进行示意。影像数目不受约束，更多的影像以及更短时间的间隔将提供更丰富的地表演化过程信息，从而提高拟合的准确度。此外，我们建议时间跨度设置为一年或者更长，这样能使时间序列影像包含一个或多个完整的植被和农作物生长周期，但这不是强制性的，只是更多的具有覆盖植被周期的影像将提升噪声检测性能。

为了进行逐像素的时间序列分析，影像必须经过严格的几何校正并达到亚像素精度。此外，影像经过辐射校正去除传感器和大气相关的噪声，并转换为大气顶部或底部反射率(TOA/BOA)，这将提高在不同日期获取的影像之间的一致性。作为一种替代方案，辐射归一化也可以用于影像辐射处理，以消除不同时间影像之间的辐射差异，使其具有类似的辐射采集条件。

从时间序列影像中，我们可以观察到：(1)影像112上建筑物151在不同时间获取的NDVI值保持相对稳定，变化较小；(2)水体152的NDVI的值全年均降低，且变化较小；(3)农田153有两个峰值，一个峰值在三月(影像对应于112)和四月(影像对应于113)，另一个峰值在八月(影像对应于117)和九月(影像对应于118)；(4)森林154在五月份增长之后，保持旺盛直到9月，然后逐渐降低。

虽然，我们可以通过同时观察多期影像来发现土地覆盖的演化模式，然而，由于时间序列影像是高维度的，因此同时在视觉上解析多个影像并不直观，这是一种劳动密集的工作。从时间序列影像中提取明确定义的特征并构造特征向量将简化该过程。例如，归一化差值植被指数NDVI，增强植被指数EVI等植被指数不仅可以突出植被和背景之间的差异，还可以降低多光谱影像的维度，从而，TSVI为分析土地覆盖动态变化提供了方便有效的指标。

根据时间演化的模式，研究区域中存在着四种主要地物类型，即建筑物，水，森林和农田，图2A-D分别显示了上述土地覆盖类型的NDVI的时间变化情况。主要特点总结如下：(1)建筑物的NDVI理论上不随时间变化，通常称之为伪不变特征，NDVI的变化可归因于非地表相关因素导致的，该特征已经广泛地用于辐射归一化。(2)水面的NDVI变化不仅受水本身辐射特征影响，还受到叶绿素和沉积物含量的影响，该变化难以利用预定的模型进行表示和描述。(3)农田是一种受到人类活动控制下的土地覆盖类型，NDVI的峰值时刻可能会偏离自然植被的峰值时刻。此外，由于同一地块一年中种植多茬农作物的影响，在其时间轨迹中可以找到多对峰谷。(4)森林具有季节性变化，如春季生长，秋季衰落。与此规则一致，NDVI在春季增长，在夏季达到峰值，之后在秋季下降。

根据上述时间曲线可以看出：农田及其森林像素包含了植被周期性变化，难以利用噪声偏离均值较大的假设进行检测。幸运的是，研究表明由植被的生长和衰落导致的NDVI变化可以使用数学模型进行较好地模拟。受此机制的启发，我们可以从森林、农作物等的TSVI中去除季节性成分，将非平稳序列转化为平稳序列，再实现噪声检测。

图3展示了该实施例中公开的方法的流程图，图4展示了本发明实施所能达到的效果，本发明将实施方法划分为八个步骤，并结合图4详细介绍本发明实施方式：

步骤311：计算时序特征

在该步骤中，逐像素

地从时间序列影像I＝＜I₁,I₂,...,I_n＞中计算植被指数。其中，NDVI的计算方法为：

其中，NIR、R分别表示近红外波段和红波段的观测值，该值使用光谱反射率或者DigitalNumber(DN)值表示。据此，构造各个像素的特征向量

其中

表示影像在时刻t_j∈T时刻VI，X^j描述了像素P^j∈P随时间的变化的情况。但是由于噪声的影响，

的值包含了畸变，需要通过本发明的方法予以检测。

图4A展示了一个农作物像素NDVI随时间的变化。

步骤312：对序列进行分类；

由于本发明的目标不是将土地覆盖分类为具有明确地理意义的土地类型，而是测试时间序列序列的平稳性，然后将非平稳序列转换为平稳序列以进行噪声检测。在数学中，平稳序列是指序列的平均值和标准差不随时间变化。在TSVI中，建筑物的土地覆盖和清澈的水没有明显的季节变化，可以看作是平稳序列。而以植被为代表的像素具有明显的周期性变化，TSVI为非平稳序列。

单位根检验用于检测序列是否平稳，若序列具有单位根，则是非平稳的，反之序列是平稳的。现有的平稳序列测试技术，例如ADF测试，丹尼尔测试(Daniel Test)。在该实施例的公开内容中，我们使用ADF测试来检查序列是否平稳，需要置信度参数，表示结果具有多大的可信度，该参数可以根据需要进行调整。

该步骤中也可以采用其他现有的单位根检验技术，其目标是将序列精确地区分为平稳序列和非平稳序列。

如果序列是非平稳的，我们将通过313，314和315进行噪声检测。反之，如果序列是平稳的，我们可以跳过模型拟合313和偏差计算314这两个步骤，直接进行噪声检测315。

图4A所示的一个农作物像素NDVI随时间变化曲线图具有明显的周期性变化，根据该步骤会将其判定为非平稳序列。

步骤313：模型拟合；

本步骤使用一定数学模型拟合离散的TSVI，数学模型可以采用非对称高斯，非线性傅里叶和double logistic等，也可以使用其他新的数学模型，其特征是尽量准确地描述植被的变化，以减小模型拟合的误差。本实施案例中使用高斯模型。

由于农作物的时间序列中可能存在多对峰谷，因而本实施案例中使用多峰模型，需要确定峰谷的对数K：我们分别设置K＝1,2,...3，并选择具有最小拟合误差的K作为参数。

计算模型参数可以使用解析法(如最小二乘法)、启发式方法(如随机梯度下降法)或者优化方法。该参数也可以根据需要进行调整。

图4A展示了时间序列拟合结果，可以看出本文方法能够有效拟合出时间序列结果。

步骤314：计算偏差；

我们通过模型拟合313，可以得到像素P^j∈P在t_i∈T时刻的预测值

差值

表示模型拟合值和曲线的拟合程度。对于P^j∈P的偏差D^j，1表示为：

通过该步骤，由于公式(7)表示的模型具有季节分量，

描述了像素实际观测与模型预测的近似程度。由云影引起的噪声将远离拟合曲线，具有较大的偏差，而干净观测值将分布在拟合曲线附近，其具有较小的偏差。

图4B所示为偏差计算结果，可以看出偏差序列去除了时间序列观察中的周期性成分。

步骤315：识别噪声；

由于时间序列偏差中去除了季节性成分，理想情况下偏差的平均值将不随时间而改变，即平稳序列。因而，干净观测值将分布在平均值附近的狭窄范围内，而噪声将远离拟合曲线。与传统的异常值检测方法类似，我们可以将噪声定义为分布均值μ加减λ倍标准差σ的范围之外的观测值。其中，λ是由用户设置的参数，我们可以通过设置λ以检测噪声。

这里可以融入先验知识进行噪声检测，如受到噪声污染的植被的NDVI被低估，可以仅将偏差小于μ-λσ的观测值视为噪声。

图4B标记了上述的方法的噪声检测结果，噪声使用十字标记。

步骤316：迭代处理；

由于我们第一次进行模型拟合313的时候，使用包含噪声的观测值，因此所获得的模型可能存在偏差，如图4A所示，其拟合曲线偏离原始观测值较大，不利于偏差的估算。去除了第一次迭代中的噪声后，可以使用如图4C所示的剩余观察来获得新函数。由于已经排除了用于模型拟合的噪声，所以预期新的估算将与实际的轮廓更加一致。通过这个新函数，我们可以重新更准确地计算偏差。并且我们可以重复步骤313，314和315以去除剩余的异常值，直到迭代次数达到最大值或者在一次迭代中没有检测到异常值。

最终的模型拟合和异常值结果见图4D。我们可以观察到剩余的观测结果与实际值一致，并且几乎所有的异常值都已被去除。

步骤317：掩膜计算；

我们将原始时间序列VI转换回影像空间并获得噪声掩膜。每个影像都有一个噪声掩膜，该掩膜的像素表示它是噪声，还是干净观测。

步骤318：噪声优化；

由于阈值设置不当或无法准确拟合曲线，因此我们可以将空间信息合并到噪声检测中。鉴于云的范围较大，我们可以删除小于阈值T_h的有效观测中的小块云或者云中的小块有效观测。

图5展示了使用本实施例得到的噪声检测结果。511中的灰色区域表示噪声掩膜，这些区域是缺失的或是被污染的，它们不再进入后续的应用或者通过一定的方法进行填补。白色区域表示清晰观测，它能够可靠地用于与土地相关的应用。

如果我们的目标是准确地获得精确的曲线，则进行步骤317和318是可选的。

图6展示的是可以实现我们的方法的计算机系统600。601是被放入计算机系统610中的时间序列影像，这些影像的部分区域被云影污染而丢失，且这些部分的准确区域未知。计算机系统610执行一系列编程指令，这些指令对从遥感影像导出的TSVI执行噪声检测。程序、原始时间序列影像以及中间结果存储在诸如硬盘或远程存储服务器的存储介质611中。处理结果或处理中间状态可以由计算机显示器612进行可视化以及参数调整。计算机600从计算机存储介质，相机存储卡或通过通信网络接收时间序列影像，并如所描述的那样检测噪声观察。612表示视频监控器或其他显示器，在612上为用户显示所获得的影像方式掩膜602，也可以将影像方式的噪声掩膜存储或发送到远程位置以进行分析。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (1)

1.一种时间序列遥感植被指数的噪声检测方法，包括以下步骤：

步骤1，计算时序特征；

假设研究区有n期不同时间获取的影像I构成：

I＝<I₁,I₂,...,I_n> (1)

各期影像的获取日期T表示为：

T＝<t₁,t₂,...,t_n> (2)

获取日期之间存在关系t_i-1≤t_i，i＝2,3,...,n，即t_i-1在t_i之前；

研究区有m个像素P表示为：

P＝{P¹,P²，...,P^m}^T (3)

像素

提取的TSVIX^j表示为：

植被指数使用归一化植被指数NDVI，或者使用增强型植被指数EVI,或者其他植被指数，使用的植被指数特征在于能够消除地表的辐射背景，并突出植被特征；

步骤2，对序列进行分类；

判断各个像素

的TSVI的X^j变化情况，将其划分到平稳序列集合S或者非平稳序列集合NS中，S和NS满足：

S∪NS＝P (6)

其中，平稳序列是指均值和标准差不随时间的变化的序列，而非平稳序列是指均值和标准差会随时间而变化的序列；平稳序列检验采用单位根检测方法，单位根检测方法特征在于能够将影像划分为平稳序列和非平稳序列；

步骤3，非平稳序列拟合；

对于非平稳序列集合NS中的像素，利用复合高斯模型对离散的TSVI进行拟合；拟合过程中，自变量为T，应变量为X，拟合结果为

模型表示为：

其中，a_k，μ_k，σ_k是待估算的参数，描述了高斯模型；α_K表示高斯成分的比例；K表示高斯成分的数目，即曲线中包含几个峰谷组合；μ_k表示均值，即达到最大值时候的t值；σ_K决定了峰的陡峭程度；

步骤4，计算偏差；

对于像素P^j∈NS，在t_i∈T时存在模型计算值

和实际观测值

其偏差

表示为：

其上标1表示第一次迭代结果；

表示了模型计算值

与实际观测值

之间的接近程度；

越大，表示二者差异较大，拟合效果越差；

越小，表示真实值和拟合值越接近，拟合效果较好；

像素P^j∈P在t_i∈T时的第一次偏差计算结果表示为：

由于偏差D^j,1消除了季节性成分，从而将非平稳序列转换为平稳序列；

步骤5，检测噪声；

对于转换后的偏差序列D^j,1和平稳序列S，分别采用观测值分布在均值附近一定范围内的假设，进行噪声检测；

对于P^j∈S，采用原始的序列X^j实施噪声检测；而对于P^j∈NS，采用偏差序列D^j,1，为了论述的方便，统一记为s^j：

然后按照以下方法判定噪声，

为掩膜数据，表示像素P^j∈P在t_i∈T是否为有效观测，1表示有效观测，而0表示噪声；

其中，μ(s^j)表示s^j的均值，σ(s^j)表示s^j的标准差，λ表示标准差的倍数；公式(11)的意义为，将距拟合曲线的距离分布在均值加减λ倍标准差的观测值作为有效观测，而分布在之外的作为的噪声；可以获得像素P^j∈P的掩膜M^j；

步骤6，迭代处理；

对去除噪声保留之后的观测值，即不考虑

的像素，重新按照步骤2进行非平稳序列拟合、步骤3进行偏差计算和步骤4进行噪声检测，直到迭代达到最大次数或在某次迭代过程中不再有噪声检出；

步骤7，噪声掩膜；

将各个像素P^j∈P的掩膜数据转换回图像空间获得噪声掩膜M；

M＝<M₁,M₂,...,M_n> (13)

其中，M_i∈M表示了影像I_i∈I上的各个像素的观测值是否为有效观测值，在后续的应用中，需要将噪声进行剔除，以获得准确的应用效果；

步骤8，掩膜优化；

根据数学形态学方法对噪声进行，获得平滑的噪声掩膜。

Images