CN103020478B - 一种海洋水色遥感产品真实性检验的方法 - Google Patents

Abstract

本发明公开了一种海洋水色遥感产品真实性检验的方法。本发明包括步骤如下：步骤(1).长时间序列遥感数据集的重构；步骤(2).基于遥感数据集的时间-空间变异特征计算与综合；步骤(3).用于检验遥感产品的实测数据集质量控制与筛选；步骤(4).基于时空匹配法则的实测-遥感数据检验数据集生成；步骤(5).符合条件的匹配数据集风险等级区划；步骤(6).遥感产品数据真实性检验的指标计算。本发明能够实现遥感产品数据与检验实测数据间的代表性与合理性进行检验，匹配数据集的筛选生成，得到不同风险等级区划的评价方法，使得评价结果具有多样性和合理性，以此降低检验过程的不确定性风险，提高检验可信度的一种实际操作方法。

Description

一种海洋水色遥感产品真实性检验的方法

技术领域

 本发明属于遥感技术领域，具体涉及是一种海洋水色遥感产品真实性检验的方法，尤其是对时间、空间变异存在的情况下造成检验结果不确定性的现实问题。

背景技术

     自1978年10月第一台水色卫星传感器CZCS（Coastal Zone Color Scanner）投入使用至今，水色卫星传感器的数量、光谱通道、辐射性能、空间分辨率等都得到了长足发展。众多的水色遥感产品的业务化应用，为海洋学相关研究提供了大范围、连续的观测数据。然而，水色遥感产品的精度验证问题一直是人们关注的热点。由于海洋表面存在复杂的时空变异特性，水色遥感产品的验证方法和评价体系仍然存在诸多争议问题。

  目前，国际相关组织（NASA Ocean Biology Processing Group's 简称：OBPG）及我国学者^[1]对水色遥感产品精度验证的主要方法是Sean W. Bailey等^[2]提出的验证算法。该方法通过时空窗口的约束法则和空间窗口内遥感数据的变异特性的判别来进行遥感产品精度验证。具体的方法是：首先通过时空窗口的约束来进行实测遥感数据集的匹配，同时，遥感空间窗口又必须满足空间均匀性的要求。时间窗口通常选择卫星国境时刻的±3h内，空间窗口可以选择3*3，5*5等。窗口内的数据同时满足有效性和均匀性的要求。

具体的匹配过程和判别法则如下：

(1)统计空间窗口内的有效像元个数(即非填充值的像元个数, NVP)和非陆地像元总数(NTP), 要求NVP>NTP/2+1 且NVP≥2, 以保证空间均值的代表性。

(2)计算有效像元的均值X 和标准差σ, 剔除X ±(1.5×σ)之外的像元, 以减少较异常数据对均值计算的影响。

(3)重新计算均值和标准差, 并计算方差系数CV=标准差/均值, 若CV≤0.15, 则通过均匀性判识。

     时间窗口的设置通常选择卫星成像时刻的                                                为时间窗口。在实际的匹配过程中，由于时间、空间变异性的存在，使得匹配得到的评价数据集具有很多不确定性。如：不同的点位会落在具有不同变异特征的海区内，得到的结果也会有所差异，导致遥感产品的评价结果就存在诸多不确定性。此外，评价指标的结果也相对单一，不能全面合理体现评价过程中存在的问题。

用于遥感数据评价的实测数据，是决定评价结果合理的重要参数之一。目前，在遥感数据精度评价和验证的实际操作过程中，往往并没有太多关注实测数据的质量控制。然而，实测数据的不确定性分析与质量控制在某种程度上已和数据本身具有同样的价值。如何量化和评价实测数据的不确定性，并对实测数据进行有效的质量控制，以及基于实测数据开展的遥感产品精度验证，是目前必须面对的一个重要问题。

    对于海洋环境观测数据中含有异常点的序列而言，目前常用的数据质量控制方法有极值检验、一致性判断、递增性判断、方法检验、格拉布斯检验以及狄克逊检验等^[3,4,5]。测量数据不确定性的来源及其随机特性决定了消除实测数据的不确定性是难以做到的，而随机现象的理想研究手段就是统计学规律分析，因此可以通过统计学的方法来降低测量数据的不确定性。对于异常值的处理上通常采用的做法是对序列数据进行各种滤波处理，以此来降低序列的随机误差。常见的滤波算法（如移动均值滤波、Savitzky-Golay滤波、Loess滤波等^[6,7]）总体上都是一些从数据间的关联性质进行考虑但物理机理不够明确的数据平滑技术。

引证文件

[1]. 孙凌 ,王晓梅,郭茂华,等. MODIS 水色产品在黄东海域的真实性检验[J]. 湖泊科学, 2009, 21(2): 298-306

[2]. Sean W. Bailey , P. Jeremy Werdell.  A multi-sensor approach for the on-orbit validation of ocean color satellite data products[J]. Remote Sensing of Environment 102 (2006) 12–23

[3]. 国家质量技术监督局2712海洋监测规范(第2部分)：数据处理与分析质量控制[S] GB17378.2,1998年.

[4]. 陈上及，海洋数据处理分析方法及其应用[M]，海洋出版社，1991.

[5]. 黄谟涛、翟国君、王瑞等，海洋异常数据的检测[J]，测绘学报，1999，28（3），269‐277。

[6]. Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing[M], Prentice-Hall, 1989.

[7]. Mitra, S.K., Digital Signal Processing, 2nd ed[M]., McGraw-Hill, 2001。

发明内容

    本发明的目的是针对现有技术的不足，提出了一种海洋水色遥感产品真实性检验的方法。

    本发明解决其技术问题所采用的技术方案具体包括如才步骤：

    步骤(1).长时间序列遥感数据集的重构；

    步骤(2).基于遥感数据集的时间-空间变异特征计算与综合；

    步骤(3).用于检验遥感产品的实测数据集质量控制与筛选；

    步骤(4).基于时空匹配法则的实测-遥感数据检验数据集生成；

    步骤(5).符合条件的匹配数据集风险等级区划；

    步骤(6).遥感产品数据真实性检验的指标计算。

所述的步骤(1)中长时间序列遥感数据集的重构通过DINEOF数据重构方法实现，该方法是基于正交函数分解方法（EOF）发展起来的一种数据重构方法，假定原始序列数据集为维的二维数据矩阵，每一行代表某一空间位置点的时间序列值，每一列代表某一时刻所有空间序列值，原始序列数据集中的一些缺测点（无数据，如云覆盖、卫星轨道未覆盖区，数据不可靠点等等），用表示，具体步骤包括如下：

1-1.数据去均值处理，计算原始序列数据集中有效数据的均值，然后令；在二维数据矩阵中随机挑选1%的有效数据组成原始交叉验证点集，用于判断最佳重构模态数，对二维数据矩阵中处于原始交叉验证点集位置的有效数据同样赋值为；

1-2.用0替换二维数据矩阵中赋值为的点，使二维数据矩阵中缺测点的初始值为数据集的无偏估计值，并令特征模态参数；

1-3.对二维数据矩阵用式（1.1）进行EOF分解，得到个特征模态，使用式（1.2）计算二维数据矩阵中缺测点的重构值，大于1；

         （1.1）

其中分别为二维数据矩阵奇异值分解（SVD分解）后对应空间特征模态、奇异值矩阵和时间特征模态，T表示矩阵转置。

     （1.2）

其中为二维数据矩阵的空间与时间下标，和分别是空间特征模态和时间特征模态的第列，为相应的奇异值，；对步骤1-3迭代NITEMAX次（NITEMAX为程序预先设定的最大迭代次数），当时，计算原始交叉验证点集重构值与原始值的均方根误差。

1-4.令，重复步骤1-3计算出对应的均方根误差，然后比较得出值最小时对应的特征模态数，其中根据观测数据的时间维数来确定。

1-5.在特征模态数为时，计算二维数据矩阵中缺测点的重构值，并用计算出的重构值替换二维数据矩阵中原有缺测点的值；二维数据矩阵中处于原始交叉验证点集位置的有效数据用原始交叉验证点集中的原始值还原，得到新的数据集；令，得到重构序列数据集。

最终的重构误差计算如式（1.3）下：

     （1.3）

其中为新的数据集中有效点的总数，分别为原始序列数据集和重构序列数据集。

所述的步骤(2)中基于遥感数据集的时间-空间变异特征计算与综合，具体如下；

2-1.时间序列变异系数直接利用时间尺度内的像元属性值进行统计计算，统计标准差计算公式如下：

          （1.4）

式中，为时间统计序列标准差，为统计序列的均值，为统计序列的样本个数（与时间尺度大小以及历史样本的统计年数相关），代表时序遥感数据属性；

时间序列变异系数的计算公式为：

                  （1.5）

2-2.空间变异即空间窗口内相关属性数据观测值的变异程度，空间窗口的大小可取，其中a为大于3的奇数；窗口内的变异系数计算方法是窗口内的标准差除以窗口内的均值，形成空间变异的时间序列数据，通过采用步骤2-1所述方法计算空间变异时间序列的变异特性；时间、空间变异系数计算完成后，得到相应属性的时间空间变异数据的网格数据，即时间变异系数评价网格数据，空间变异变异系数评价网格数据。

2-3.对时间变异变异系数评价网格数据和空间变异系数评价网格数据，采用归一化网格图层叠加的形式进行综合，得到综合变异系数评价网格数据。

所述的步骤（3）中用于检验遥感产品的实测数据集质量控制与筛选，具体如下：

3-1.用ADF方法检验时间序列经过二阶差分处理后的平稳性，通过大量实测数据试验得出，80%的实测数据满足一阶差分平稳性检验，95%以上的实测数据都满足二阶差分平稳性检验，由此可假设测量时间序列的二阶导数是一个常数，模型可表示如下：

                 （1.6）

式中，为常数，为系统误差，为某一时间范围；但实际测量数据中，由于不确定性的存在，不可能表现为常数，为了便于利用统计学知识进行分析，将模型表示为：

                （1.7）

所述的实测数据的不确定性可分为系统误差和随机误差，系统误差通常表现为整体偏差，较难发现，随机误差则易于统计量化。为了便于分析，将式（1.7）中的两项表现误差参数合并为下式：

      （1.8）

对实测数据而言，其测量项目数据和时间数据都是以离散的数据集合来表示，如某一测量项目的数据集可表示为，时间数据集为，则式（1.8）可用二阶差分形式表示如下：

                    （1.9）

       （1.10）

3-2.对所有观测项目的二阶差分值进行归一化处理，以消除不同量纲带来的影响，对极差归一化方法为：

      （1.11）

3-3.对多参数的实测数据进行联合质量控制，首先判断多参数间的相关性，具体如下：

   （1.12）

 其次，根据相关性进行权值确定，即相关性越强，权值越大，具体如下：

                （1.13）

式中，为第个观测数据加权融合后的结果，为权值，通过对序列进行统计分析，将序列的以外的数据作为离异数据来处理，对离异后的数据进行识别并通过移动均值滤波、Savitzky-Golay滤波、Loess滤波对离异的数据进行重新赋值。

所述的步骤(4)中基于时空匹配法则的实测-遥感数据检验数据集生成，具体通过空间均匀性3个判识准则和时间窗口确定。

卫星数据的空间均匀性3个判识准则如下:

(1)统计空间窗口内的有效像元个数 (即非填充值的像元个数)和非陆地像元总数, 要求,且,以保证空间均值的代表性。

(2)计算有效像元的均值和标准差, 剔除之外的像元, 减少较异常数据对均值计算的影响.

(3)重新计算均值和标准差, 并计算方差系数，且=标准差/均值, 若≤0.15, 则通过均匀性判识。

     时间窗口的设置通常选择卫星成像时刻的为时间窗口。

所述的步骤(5)中符合条件的匹配数据集风险等级区划，具体如下：

根据步骤(2)中计算得出的时间变异变异系数评价网格数据、空间变异变异系数评价网格数据、综合变异系数评价网格数据，同时结合变异系数进行空间风险区划分，将上述3个网格数据划分为5个风险区划等级区间：

Level_1(等级一)：,

Level_2(等级二)：, 

Level_3(等级三)：, 

Level_4(等级四)：, 

Level_5(等级五)：。

    上述等级区划中，为统计区域的变异系数的均值，为统计区域变异系数的标准差。

时间、空间、以及时间-空间综合3种风险区划的生成后，根据实测数据的经纬度位置，匹配风险区域的等级，对实测、遥感数据集进行合理的精度评价。

 所述的步骤(6)中遥感产品数据真实性检验的指标计算，具体如下：

指标计算是遥感产品精度的直观结果，不同风险等级匹配数据下遥感产品数据真实性检验的结果包括如下5个指标：

（1）平均绝对相对误差RE：

（2）均方根误差：

（3）平均相对偏差：

上述3个计算公式中：为检验样本集数据对的个数，为遥感数据产品结果值，为实际测量数据结果值；

（4）相关系数（）：

            

其中，和分别为遥感数据产品和实际测量数据的算术平均值；值范围介于-1与+1之间，即-1≤r≤1；

当时直线的斜率为正，称正相关；

当时直线的斜率为负，称负相关；

当|r|＝1时全部数据点都落在拟合直线上；

若r＝0则与之间完全不相关；

值愈接近±1则它们之间的线性关系愈密切；

（5）误差标准差（SD）：

其中，为数据集的算术平均值，为检验样本集数据对的个数，为遥感数据产品结果值，为实际测量数据结果值。

本发明有益效果如下：

本发明能够实现遥感产品数据与检验实测数据间的代表性与合理性进行检验，匹配数据集的筛选生成，得到不同风险等级区划的评价方法，使得评价结果具有多样性和合理性，以此降低检验过程的不确定性风险，提高检验可信度的一种实际操作方法。

针对海洋水色遥感产品数据验证过程中，水色属性时空变异情况的存在而造成的精度检验的不确定性，该发明考虑了时空变异的影响，提供了不同风险等级区划的评价方法，使得评价结果具有多样性和合理性，是海洋水色遥感产品精度验证方法的一项创新。该方法具有实施过程可行强、评价计算过程和指标便于程序集成、评价结果严密合理，可用于海洋遥感领域相关产品验证，具有较高的行业推广价值和应用前景。该方法已经在海洋公益行业项目“中国近海海-气二氧化碳通量遥感监测评估系统研究与示范”的相关示范系统中进行集成，取得了良好的效果。

附图说明

    图1为本发明中实测数据集质量控制与筛选流程图。

    图2为本发明本发明评价基本流程图。

具体实施方式

下面结合附图对本发明作进一步说明。

如图1、图2所示，一种海洋水色遥感产品真实性检验的方法，具体包括如下步骤：

    步骤(1).长时间序列遥感数据集的重构；

    步骤(2).基于遥感数据集的时间-空间变异特征计算与综合；

    步骤(3).用于检验遥感产品的实测数据集质量控制与筛选；

    步骤(4).基于时空匹配法则的实测-遥感数据检验数据集生成；

    步骤(5).符合条件的匹配数据集风险等级区划；

    步骤(6).遥感产品数据真实性检验的指标计算。

所述的步骤(1)中长时间序列遥感数据集的重构通过DINEOF数据重构方法实现，该方法是基于正交函数分解方法（EOF）发展起来的一种数据重构方法，假定原始序列数据集为维的二维数据矩阵，每一行代表某一空间位置点的时间序列值，每一列代表某一时刻所有空间序列值，原始序列数据集中的一些缺测点（无数据，如云覆盖、卫星轨道未覆盖区，数据不可靠点等等），用表示，具体步骤包括如下：

1-1.数据去均值处理，计算原始序列数据集中有效数据的均值，然后令；在二维数据矩阵中随机挑选1%的有效数据组成原始交叉验证点集，用于判断最佳重构模态数，对二维数据矩阵中处于原始交叉验证点集位置的有效数据同样赋值为；

1-2.用0替换二维数据矩阵中赋值为的点，使二维数据矩阵中缺测点的初始值为数据集的无偏估计值，并令特征模态参数；

1-3.对二维数据矩阵用式（1.1）进行EOF分解，得到个特征模态，使用式（1.2）计算二维数据矩阵中缺测点的重构值，大于1；

         （1.1）

其中分别为二维数据矩阵奇异值分解（SVD分解）后对应空间特征模态、奇异值矩阵和时间特征模态，T表示矩阵转置。

     （1.2）

其中为二维数据矩阵的空间与时间下标，和分别是空间特征模态和时间特征模态的第列，为相应的奇异值，；对步骤1-3迭代NITEMAX次（NITEMAX为程序预先设定的最大迭代次数），当时，计算原始交叉验证点集重构值与原始值的均方根误差。

1-4.令，重复步骤1-3计算出对应的均方根误差，然后比较得出值最小时对应的特征模态数，其中根据观测数据的时间维数来确定。

1-5.在特征模态数为时，计算二维数据矩阵中缺测点的重构值，并用计算出的重构值替换二维数据矩阵中原有缺测点的值；二维数据矩阵中处于原始交叉验证点集位置的有效数据用原始交叉验证点集中的原始值还原，得到新的数据集；令，得到重构序列数据集。

最终的重构误差计算如式（1.3）下：

     （1.3）

其中为新的数据集中有效点的总数，分别为原始序列数据集和重构序列数据集。

所述的步骤(2)中基于遥感数据集的时间-空间变异特征计算与综合，具体如下；

2-1.时间序列变异系数直接利用时间尺度内的像元属性值进行统计计算，统计标准差计算公式如下：

          （1.4）

式中，为时间统计序列标准差，为统计序列的均值，为统计序列的样本个数（与时间尺度大小以及历史样本的统计年数相关），代表时序遥感数据属性；

时间序列变异系数的计算公式为：

                  （1.5）

2-2.空间变异即空间窗口内相关属性数据观测值的变异程度，空间窗口的大小可取，其中a为大于3的奇数；窗口内的变异系数计算方法是窗口内的标准差除以窗口内的均值，形成空间变异的时间序列数据，通过采用步骤2-1所述方法计算空间变异时间序列的变异特性；时间、空间变异系数计算完成后，得到相应属性的时间空间变异数据的网格数据，即时间变异系数评价网格数据，空间变异变异系数评价网格数据。

2-3.对时间变异变异系数评价网格数据和空间变异系数评价网格数据，采用归一化网格图层叠加的形式进行综合，得到综合变异系数评价网格数据。

所述的步骤（3）中用于检验遥感产品的实测数据集质量控制与筛选过程中，实测数据的可靠程度，在一定程度上也直接影响遥感精度验证的结果。因此，对实测数据进行质量控制也显得尤为重要，这里提出一种针对海量、高频连续观测海洋属性相关实测数据（如走航数据、浮标数据等）差分统计的质量控制方法，具体如下：

3-1.用ADF方法检验时间序列经过二阶差分处理后的平稳性，通过大量实测数据试验得出，80%的实测数据满足一阶差分平稳性检验，95%以上的实测数据都满足二阶差分平稳性检验，由此可假设测量时间序列的二阶导数是一个常数，模型可表示如下：

                 （1.6）

式中，为常数，为系统误差，为某一时间范围；但实际测量数据中，由于不确定性的存在，不可能表现为常数，为了便于利用统计学知识进行分析，将模型表示为：

                （1.7）

所述的实测数据的不确定性可分为系统误差和随机误差，系统误差通常表现为整体偏差，较难发现，随机误差则易于统计量化。为了便于分析，将式（1.7）中的两项表现误差参数合并为下式：

      （1.8）

对实测数据而言，其测量项目数据和时间数据都是以离散的数据集合来表示，如某一测量项目的数据集可表示为，时间数据集为，则式（1.8）可用二阶差分形式表示如下：

                    （1.9）

    （1.10）

3-2. 3-2.设某仪器测量平台上有个传感器，能够同步观测个参数，由于观测参数的量纲不同，所以得到的二阶差分的值必然会不同，因此，在进行分析之前需要对所有观测项目的二阶差分值进行归一化处理，以消除不同量纲带来的影响，这里采用极差归一化方法为：

      （1.11）

3-3.对多参数的实测数据进行联合质量控制，首先判断多参数间的相关性，具体如下，参看图1：

   （1.12）

 其次，根据相关性进行权值确定，即相关性越强，权值越大，具体如下：

                （1.13）

式中，为第个观测数据加权融合后的结果，为权值，通过对序列进行统计分析，将序列的以外的数据作为离异数据来处理，对离异后的数据进行识别并通过移动均值滤波、Savitzky-Golay滤波、Loess滤波对离异的数据进行重新赋值；

对测量数据的不确定分析与质量控制，通常都是通过整体的统计特征来表现，而对每一个测量值的个体而言并不能给出一个确切的不确定性表示方法，最直观的描述就是偏移测量序列统计平均值的大小。所以，对每一个测量值的个体，处理原则是尽可能的保护原始数据，而对确实存在不符合整体统计规律的个体，才进行识别和处理，以降低其不确定性。在统计学中，样本序列的标准差表示样本内各个观测值的变异程度，通常将测量序列的以外的数据作为离异数据来处理，对于识别出来的离异数据，其估测值可以利用前后观测值的相互关系进行内插估算，具体的内插方法可以参考移动均值滤波、Savitzky-Golay滤波、Loess滤波等成熟的滤波技术。

所述的步骤(4)中基于时空匹配法则的实测-遥感数据检验数据集生成，具体通过空间均匀性3个判识准则和时间窗口确定，参看图2。

卫星数据的空间均匀性3个判识准则如下:

(1)统计空间窗口内的有效像元个数 (即非填充值的像元个数)和非陆地像元总数, 要求,且,以保证空间均值的代表性。

(2)计算有效像元的均值和标准差, 剔除之外的像元, 减少较异常数据对均值计算的影响.

(3)重新计算均值和标准差, 并计算方差系数，且=标准差/均值, 若≤0.15, 则通过均匀性判识。

     时间窗口的设置通常选择卫星成像时刻的为时间窗口。

所述的步骤(5)中符合条件的匹配数据集风险等级区划，具体如下：

如图2所示，根据步骤(2)中计算得出的时间变异变异系数评价网格数据、空间变异变异系数评价网格数据、综合变异系数评价网格数据，同时结合变异系数进行空间风险区划分，将上述3个网格数据划分为5个风险区划等级区间：

Level_1(等级一)：,

Level_2(等级二)：, 

Level_3(等级三)：, 

Level_4(等级四)：, 

Level_5(等级五)：。

    上述等级区划中，为统计区域的变异系数的均值，为统计区域变异系数的标准差。

时间、空间、以及时间-空间综合3种风险区划的生成后，根据实测数据的经纬度位置，匹配风险区域的等级，对实测、遥感数据集进行合理的精度评价。

 所述的步骤(6)中遥感产品数据真实性检验的指标计算，具体如下：

指标计算是遥感产品精度的直观结果，不同风险等级匹配数据下遥感产品数据真实性检验的结果包括如下5个指标：

（1）平均绝对相对误差RE：

（2）均方根误差：

（3）平均相对偏差：

上述3个计算公式中：为检验样本集数据对的个数，为遥感数据产品结果值，为实际测量数据结果值；

（4）相关系数（）：

            

其中，和分别为遥感数据产品和实际测量数据的算术平均值；值范围介于-1与+1之间，即-1≤r≤1；

当时直线的斜率为正，称正相关；

当时直线的斜率为负，称负相关；

当|r|＝1时全部数据点都落在拟合直线上；

若r＝0则与之间完全不相关；

值愈接近±1则它们之间的线性关系愈密切；

（5）误差标准差（SD）：

其中，为数据集的算术平均值，为检验样本集数据对的个数，为遥感数据产品结果值，为实际测量数据结果值。

Claims (1)

1.一种海洋水色遥感产品真实性检验的方法,其特征在于包括如下步骤：

步骤(1).长时间序列遥感数据集的重构；

步骤(2).基于遥感数据集的时间-空间变异特征计算与综合；

步骤(3).用于检验遥感产品的实测数据集质量控制与筛选；

步骤(4).基于时空匹配法则的实测-遥感数据检验数据集生成；

步骤(5).符合条件的匹配数据集风险等级区划；

步骤(6).遥感产品数据真实性检验的指标计算；

步骤(1)中长时间序列遥感数据集的重构通过DINEOF数据重构方法实现，该方法是基于正交函数分解方法发展起来的一种数据重构方法，假定原始序列数据集X⁰为m0×n0维的二维数据矩阵，每一行代表某一空间位置点的时间序列值，每一列代表某一时刻所有空间序列值，原始序列数据集X⁰中的一些缺测点，用NaN表示，具体步骤包括如下：

1-1.数据去均值处理，计算原始序列数据集X⁰中有效数据的均值然后令在二维数据矩阵X中随机挑选1％的有效数据组成原始交叉验证点集X_{cross-validation}，用于判断最佳重构模态数，对二维数据矩阵X中处于原始交叉验证点集X_{cross-validation}位置的有效数据同样赋值为NaN；

1-2.用0替换二维数据矩阵X中赋值为NaN的点，使二维数据矩阵X中缺测点的初始值为数据集的无偏估计值，并令特征模态参数P＝1；

1-3.对二维数据矩阵X用式(1.1)进行EOF分解，得到P个特征模态，使用式(1.2)计算二维数据矩阵X中缺测点的重构值X_i,j，P大于1；

X_i,j＝USV^T   (1.1)

其中U,S,V分别为二维数据矩阵X奇异值分解后对应空间特征模态、奇异值矩阵和时间特征模态，T表示矩阵转置；

$X_{i,j}=\underset{i=N}{\overset{P}{\Σ}}{\mathrm{\ρ}}_{t}u{\left(t\right)}_{i}v{\left(t\right)}_j---\left(1.2\right)$

其中i,j为二维数据矩阵X的空间与时间下标，u(t)和v(t)分别是空间特征模态U和时间特征模态V的第t^th列，ρ_t为相应的奇异值，t＝1,...,P；对步骤1-3迭代NITEMAX次，所述的NITEMAX为程序预先设定的最大迭代次数；当P＝1时，计算原始交叉验证点集X_{cross-validation}重构值与原始值的均方根误差Rp；

1-4.令P＝2,...,k_max，重复步骤1-3计算出对应的均方根误差R^P，然后比较得出R^P值最小时对应的特征模态数其中k_max根据观测数据的时间维数来确定；

1-5.在特征模态数为时，计算二维数据矩阵X中缺测点的重构值并用计算出的重构值替换二维数据矩阵X中原有缺测点的值；二维数据矩阵X中处于原始交叉验证点集X_{cross-validation}位置的有效数据用原始交叉验证点集X_{cross-validation}中的原始值还原，得到新的据集X1；令X^re＝X1+X_mean，得到重构序列数据集X^re，其中，X_mean为二维数据矩阵X的平均值；

最终的重构误差计算如式(1.3)下：

$\stackrel{\‾}{R}=\sqrt{\frac{1}{N}\underset{t=1}{\overset{N}{\Σ}}{(X_t^{re}-X_t^0)}^2}---\left(1.3\right)$

其中N为新的数据集X1中有效点的总数，分别为原始序列数据集和重构序列数据集；

步骤(2)中基于遥感数据集的时间-空间变异特征计算与综合，具体如下；

2-1.时间序列变异系数直接利用时间尺度内的像元属性值进行统计计算，统计标准差计算公式如下：

$S_l\left(X_{i,j}\right)=\sqrt{\frac{\Σ{(x_{i,j}-\stackrel{\‾}{x})}^2}{num-1}}---\left(1.4\right)$

式中，S_l(X_i,j)为时间统计序列标准差，为统计序列的均值，num为统计序列的样本个数，l代表时序遥感数据属性；

时间序列变异系数的计算公式为：

${\mathrm{CV}}_l=\frac{S_l\left(X_{i,j}\right)}{\stackrel{\‾}{x}}\×100\%---\left(1.5\right)$

2-2.空间变异即空间窗口内相关属性数据观测值的变异程度，空间窗口的大小为a×a，其中a为大于3的奇数；窗口内的变异系数计算方法是窗口内的标准差除以窗口内的均值，形成空间变异的时间序列数据，通过采用步骤2-1所述方法计算空间变异时间序列的变异特性；时间、空间变异系数计算完成后，得到相应属性的时间空间变异数据的网格数据，即时间序列变异系数评价网格数据，空间变异系数评价网格数据；

2-3.对时间序列变异系数评价网格数据和空间变异系数评价网格数据，采用归一化网格图层叠加的形式进行综合，得到综合变异系数评价网格数据；

步骤(3)中用于检验遥感产品的实测数据集质量控制与筛选，具体如下：

3-1.用ADF方法检验时间序列经过二阶差分处理后的平稳性，通过大量实测数据试验得出，80％的实测数据满足一阶差分平稳性检验，95％以上的实测数据都满足二阶差分平稳性检验，由此可假设测量时间序列f(t)的二阶导数是一个常数，模型可表示如下：

$\{\begin{array}{c}\frac{d^{2}f\left(t\right)}{{\mathrm{dt}}^2}=a+\mathrm{\Δ}\mathrm{\ϵ}\left(t\right)\\ t\∈\[t-\mathrm{\Δ}T,t+\mathrm{\Δ}T\]\end{array}---\left(1.6\right)$

式中，a为常数，△ε(t)为系统误差，2△T为某一时间范围；由于实测数据的不确定性，同时为了便于统计学知识分析，将公式(1.6)模型表示为：

$\frac{d^{2}f\left(t\right)}{{\mathrm{dt}}^2}=a\left(t\right)+\mathrm{\Δ}\mathrm{\ϵ}\left(t\right)---\left(1.7\right)$

所述的实测数据的不确定性包括系统误差△ε(t)和随机误差a(t)，为了便于分析，将式(1.7)中的两项表现误差参数合并为下式：

$\frac{d^{2}f\left(t\right)}{{\mathrm{dt}}^2}=a\left(t\right)+\mathrm{\Δ}\mathrm{\ϵ}\left(t\right)=\mathrm{\Δ}E\left(t\right)---\left(1.8\right)$

对实测数据而言，其测量项目数据和时间数据都是以离散的数据集合来表示，即：测量项目i的数据集可表示为{X(i)₁,X(i)₂,…,X(i)_k,…,X(i)_n}，测量项目j的数据集可表示为{X(j)₁,X(j)₂,…,X(j)_k,…,X(j)_n}，时间数据集为{T₁,T₂,…,T_k,…,T_n}，则式(1.8)对于测量项目i和j可用二阶差分形式表示如下：

$F^{\′}(X{\left(i\right)}_k,T_k)=\frac{X{\left(i\right)}_{k+1}-X{\left(i\right)}_k}{T_{k+1}-T_k},$

$F^{\′}(X{\left(j\right)}_k,T_k)=\frac{X{\left(j\right)}_{k+1}-X{\left(j\right)}_k}{T_{k+1}-T_k}---\left(1.9\right)$

$F^{\′\′}(X{\left(i\right)}_k,T_k)=\frac{F^{\′}(X{\left(i\right)}_{k+1},T_{k+1})-F^{\′}(X{\left(i\right)}_k,T_k)}{T_{k+1}-T_k},$

$F^{\′\′}(X{\left(j\right)}_k,T_k)=\frac{F^{\′}(X{\left(j\right)}_{k+1},T_{k+1})-F^{\′}(X{\left(j\right)}_k,T_k)}{T_{k+1}-T_k}---\left(1.10\right)$

3-2.对所有观测项目的二阶差分值进行归一化处理，以消除不同量纲带来的影响，对测量项目i和j极差归一化方法为：

$f^{\′\′}(X{\left(i\right)}_k,T_k)=\frac{F^{\′\′}(X{\left(i\right)}_k,T_k)-\min \left(F^{\′\′}\right(X{\left(i\right)}_k,T_k)}{\max (F^{\′\′}\left(X{\left(i\right)}_k,T_k\right)-\min (F^{\′\′}\left(X{\left(i\right)}_k,T_k\right)},$

$f^{\′\′}(X{\left(j\right)}_k,T_k)=\frac{F^{\′\′}(X{\left(j\right)}_k,T_k)-\min \left(F^{\′\′}\right(X{\left(j\right)}_k{T}_k)}{\max (F^{\′\′}\left(X{\left(j\right)}_k,T_k\right)-\min (F^{\′\′}\left(X{\left(j\right)}_k,T_k\right)}---\left(1.11\right)$

3-3.对多参数的实测数据进行联合质量控制，首先判断多参数间的相关性，具体如下：

$\left\{\begin{array}{c}r\left(ij\right)=\frac{n\Σf^{\′\′}(X{\left(i\right)}_k,T_k)f^{\′\′}(X{\left(j\right)}_k,T_k)-\Σf^{\′\′}(X{\left(i\right)}_k,T_k)\Σf^{\′\′}(X{\left(j\right)}_k,T_k)}{\sqrt{n\Σf^{\′\′}{(X{\left(i\right)}_k,T_k)}^2-{(\Σf^{\′\′}(X{\left(i\right)}_k,T_k\left)\right)}^2}\sqrt{n\Σf^{\′\′}{(X{\left(j\right)}_k,T_k)}^2-{(\Σf^{\′\′}(X{\left(j\right)}_k,T_k\left)\right)}^2}}---\left(1.12\right)\\ i=1,2,\mathrm{...},m;j=1,2,\mathrm{...},m\end{array}\right.$

其次，根据相关性进行权值确定，即相关性越强，权值越大，具体如下：

$\left\{\begin{array}{c}I(X{\left(i\right)}_k,T_k)=\underset{j=1}{\overset{m}{\Σ}}\mathrm{\ω}\left(ij\right)f^{\′\′}(X{\left(j\right)}_k,T_k)\\ \mathrm{\ω}\left(ij\right)=\frac{\left|r\left(ij\right)\right|}{\underset{j=1}{\overset{m}{\Σ}}\left|r\left(ij\right)\right|}\\ i=1,2,\mathrm{...},m;j=1,2,\mathrm{...},m\end{array}\right.---\left(1.13\right)$

式中，I(X(i)_k,T_k)为第i个观测数据加权融合后的结果，ω(ij)为权值，通过对I(X(i)_k,T_k)序列进行统计分析，将序列[μ-3σ,μ+3σ]的以外的数据作为离异数据来处理，对离异后的数据进行识别并通过移动均值滤波、Savitzky-Golay滤波、Loess滤波对离异的数据进行重新赋值；其中，σ为有效像元的标准差，μ为统计区域的变异系数的均值；

步骤(4)中基于时空匹配法则的实测-遥感数据检验数据集生成，具体通过空间均匀性3个判识准则和时间窗口确定：

卫星数据的空间均匀性3个判识准则如下:

(1)统计空间窗口内的有效像元个数NVP和非陆地像元总数NTP,要求 $NVP>\frac{NTP}{2}+1,$ 且NVP≥2；

(2)计算有效像元的均值S和标准差σ,剔除S±(1.5×σ)之外的像元；

(3)重新计算均值和标准差,并计算方差系数CV，且CV＝标准差/均值,若CV≤0.15,则通过均匀性判识；

时间窗口的设置通常选择卫星成像时刻的±3h为时间窗口；

步骤(5)中符合条件的匹配数据集风险等级区划，具体如下：

根据步骤(2)中计算得出的时间变异系数评价网格数据、空间变异系数评价网格数据、综合变异系数评价网格数据，同时结合变异系数进行空间风险区划分，将上述3个网格数据划分为5个风险区划等级区间：

Level_1(等级一)：[0,μ-2std),

Level_2(等级二)：[μ-2std,μ-std),

Level_3(等级三)：[μ-std,μ+std],

Level_4(等级四)：(μ+std,μ+2std],

Level_5(等级五)：>μ+2std；

上述等级区划中，μ为统计区域的变异系数的均值，std为统计区域变异系数的标准差；

时间、空间、以及时间-空间综合3种风险区划的生成后，根据实测数据的经纬度位置，匹配风险区域的等级，对实测、遥感数据集进行合理的精度评价；

步骤(6)中遥感产品数据真实性检验的指标计算，具体如下：

指标计算是遥感产品精度的直观结果，不同风险等级匹配数据下遥感产品数据真实性检验的结果包括如下5个指标：

(1)平均绝对相对误差RE：

$RE=\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}(\frac{\left|r_i-s_i\right|}{\left|s_i\right|}\×100)$

(2)均方根误差：

$RMSE=\sqrt{\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}{(r_i-s_i)}^2}$

(3)平均相对偏差：

$BIAS=\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}(\frac{(r_i-s_i)}{\left|s_i\right|}\×100)$

上述3个计算公式中：M为检验样本集数据对的个数，r_i为遥感数据产品结果值，s_i为实际测量数据结果值；

(4)相关系数γ：

$\mathrm{\γ}=\frac{\Σ(r_i-\stackrel{\‾}{r})(s_i-\stackrel{\‾}{s})}{\sqrt{\Σ{(r_i-r)}^2\Σ{(r_i-\stackrel{\‾}{s})}^2}}$

其中，和分别为遥感数据产品r和实际测量数据s的算术平均值；γ值范围介于-1与+1之间，即-1≤r≤1；

当γ>0时直线的斜率为正，称正相关；

当γ<0时直线的斜率为负，称负相关；

当|r|＝1时全部数据点(r,s)都落在拟合直线上；

若r＝0则r与s之间完全不相关；

γ值愈接近±1则它们之间的线性关系愈密切；

(5)误差标准差SD：

$\left\{\begin{array}{c}SD=\sqrt{\frac{1}{M}\underset{i=1}{\overset{M}{\&Sigma;}}{(d_i-\stackrel{\&OverBar;}{d})}^2}\\ d_i=r_i-s_i\end{array}\right.$

其中，为数据集d的算术平均值，M为检验样本集数据对的个数，r_i为遥感数据产品结果值，s_i为实际测量数据结果值。