CN109164471A - 一种基于主成分分析的区域划分方法 - Google Patents

Abstract

本发明提供了一种基于主成分分析的区域划分方法，首先获取GNSS测站的坐标时间序列观测值，并根据坐标时间序列观测值获得残差坐标时间序列，然后计算所有测站的公共历元，并基于公共历元构建观测样本矩阵，接下来构建所有测站的协方差矩阵B_n×n，然后对协方差矩阵B_n×n正交分解得到正交矩阵P；接下来根据观测样本矩阵和正交矩阵获得载荷矩阵A，从正交矩阵P中获得特征矢量，根据特征矢量的特征值的累加结果，选择出对应的K个累积特征矢量；再对K个累积特征矢量进行再次处理，对正交矩阵P进行线性变换，获得线性变化后的载荷矩阵B’和主成分G。最后对线性变化后的载荷矩阵B’进行分析，获得区域划分结果。本发明可以大大提高区域划分的准确性。

Description

一种基于主成分分析的区域划分方法

技术领域

本发明涉及GNSS数据精密处理技术领域，尤其涉及一种基于主成分分析的区域划分方法。

背景技术

近年来，国内外进行了各种GPS/GNSS(Global Positioning System全球定位系统和Global Navigation Satellite System全球导航卫星系统)监测网络的建立，如中国地壳运动观测网络和中国大陆构造环境监测网络，美国的PBO网络，欧洲的EPN等。这些GPS/GNSS监测网络的投入运行，产生了大量的呈级数增长的观测数据。这些多年积累的观测数据构成了GNSS坐标时间序列。

目前，现有技术中通常采用叠加算法对于影响因素在空间分布上进行均匀性假设，然而GNSS坐标时间序列受到众多因素的影响，其在空间分布上并不均匀，从而使得利用其进行速度场估计的准确度大打折扣，也无法进行准确地区域划分。

由此可见，现有技术中的方法存在无法进行准确地区域划分的技术问题。

发明内容

本发明实施例提供了一种基于主成分分析的区域划分方法，用以解决或者至少部分解决现有技术中的方法存在无法进行准确地区域划分的技术问题。

本发明提供了一种基于主成分分析的区域划分方法，包括：

步骤S1：获取GNSS测站的坐标时间序列观测值，并根据所述坐标时间序列观测值获得所述GNSS测站的残差坐标时间序列；

步骤S2：计算GNSS网内所有测站的公共历元，并确定所述公共历元个数m，其中，m为大于0的整数；

步骤S3：基于所述公共历元构建所述GNSS网内所有测站在预设方向上的观测样本矩阵X(m,n)，其中，m表示公共历元个数，n表示GNSS网内测站的数目，n为大于0的整数，x_i,j为观测样本矩阵中的元素，表示第j个测站在第i个历元的误差观测值；

步骤S4：根据所有测站的残差坐标时间序列，构建所述所有测站在预设方向上的协方差矩阵B_n×n，其中，B_n×n中的对角元素用以表征测站的残差坐标时间序列的精度信息，非对角线元素用以表征相应测站之间的相关性；

步骤S5：对所述协方差矩阵B_n×n进行正交分解，获得正交矩阵P；

步骤S6：根据所述观测样本矩阵和所述正交矩阵，获得载荷矩阵A，其中，载荷矩阵中的A_i',j'表示第i'个主成分在第j'个观测量上的空间响应，所述主成分为所述正交矩阵P的主成分分量，所述观测量为协方差矩矩阵B_n×n中的协方差构成的观测量；

步骤S7：从所述正交矩阵P中获得特征矢量，并根据特征矢量的特征值的累加结果，选择出对应的K个累积特征矢量；

步骤S8：对所述K个累积特征矢量进行再次处理，对所述正交矩阵P中的m个变量的相关结构之间进行线性变换，并获得为线性变化后的载荷矩阵B’和线性变化后的主成分G，使得处理前后这K个空间型解释原场总方差对百分率保持不变；

步骤S9：对所述线性变化后的载荷矩阵B’进行分析，获得区域划分结果。

进一步地，步骤S1中，根据所述坐标时间序列观测值获得所述GNSS测站的残差坐标时间序列，包括：

通过对坐标时间序列观测值进行最小二乘拟合，获得所述GNSS测站的残差坐标时间序列；或

通过预设产品获得所述GNSS测站的残差坐标时间序列。

进一步地，步骤S2具体包括：

将每个测站的时间序列坐标时间序列观测值表示为位置随历元变化的函数，其中，一个历元对应一个位置；

根据上述函数，计算GNSS网内所有测站的公共历元，并确定所述公共历元的个数。

进一步地，步骤S3具体包括：

基于步骤S2中的公共历元，获得GNSS网内的任意一个测站均有m个历元的观测值，将其计作x，代表一个m行1列的向量，即x_m×1；

将满足上述的n个测站，表示为m行n列的向量，进而得到观测样本矩阵X(m,n)。

进一步地，步骤S4具体包括：

通过下述运算计算获得所述协方差矩阵B_n×n，其中，m表示公共历元个数，X为为观测样本矩阵，X^T为矩阵X的转置。

进一步地，步骤S5具体包括：

通过如下运算B＝PΛP^T，获得正交矩阵P，其中，Λ为特征值构成的对角矩阵，矩阵P对应协方差矩阵B的主成分。

进一步地，步骤S6具体包括：

通过如下运算X_m×n＝A_m×nP_n×n，获得载荷矩阵A，其中，X_m×n为观测样本矩阵，P_n×n为正交矩阵。

进一步地，步骤S8的具体实现包括如下变换：

X_m×n＝A_m×nP_n×n

X_m×n＝A_m×nRR^TP_n×n

X_m×n＝B'G,B'＝A_m×nR,G＝R^TP_n×n

其中，B'为线性变化后的载荷矩阵，G为线性变化后的主成分，R为变换矩阵，R^T为矩阵R的转置。

进一步地，步骤S9具体包括：

根据线性变化后的载荷矩阵B’，寻找在矩阵B’上对应的前k个主成分，以所述前k个主成分对应的载荷系数进行区域分块的划分。

本申请实施例中的上述一个或多个技术方案，至少具有如下一种或多种技术效果：

在本发明提供的方法中，首先获取GNSS测站的坐标时间序列观测值，并根据所述坐标时间序列观测值获得所述GNSS测站的残差坐标时间序列；并计算GNSS网内所有测站的公共历元，然后基于所述公共历元构建所述GNSS网内所有测站在预设方向上的观测样本矩阵X(m,n)，接下来根据所有测站的残差坐标时间序列，构建所述所有测站在预设方向上的协方差矩阵B_n×n，跟着对所述协方差矩阵B_n×n进行正交分解，获得正交矩阵P；然后根据所述观测样本矩阵和所述正交矩阵，获得载荷矩阵A，再从所述正交矩阵P中获得特征矢量，并根据特征矢量的特征值的累加结果，选择出对应的K个累积特征矢量；接着对所述K个累积特征矢量进行再次处理，对所述正交矩阵P中的m个变量的相关结构之间进行线性变换，并获得为线性变化后的载荷矩阵B’和线性变化后的主成分G。使得处理前后这K个空间型解释原场总方差对百分率保持不变；最后对所述线性变化后的载荷矩阵B’进行分析，获得区域划分结果。本发明提供的方案，通过协方差进行降维处理，获取表征各个测站间关系的物理量以实现降维，得到载荷矩阵，并进行了进一步的约束：对特征矢量进行总方差不变情况下的线性变换，可以实现个空间型对原变量场区域结构对呈现。载荷矩阵B’的载荷向量在不同的约束标准下可以反映多个测站残差时间序列的区域性相关结构，从而可以进行准确地区域划分，适用于解决多个GNSS测站时间序列中出现共同误差及其区域划分等问题，并提高基于GNSS坐标时间序列的衍生品的质量。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例中一种基于主成分分析的区域划分方法的流程图。

具体实施方式

本发明实施例提供了一种基于主成分分析的区域划分方法，用以改善现有技术存在无法进行准确的区域划分的技术问题。

本申请实施例中的技术方案，总体思路如下：

针对GPS/GNSS网内多个测站受同一来源影响而呈现出区域性的这一特点，本发明提供了一种基于主成分分析的区域划分方法，对特征矢量进行总方差不变情况下的线性变换，可以实现个空间型对原变量场区域结构对呈现。载荷矩阵B’的载荷向量在不同的约束标准下可以反映多个测站残差时间序列的区域性相关结构，从而可以进行准确地区域划分，适用于解决多个GNSS测站时间序列中出现共同误差及其区域划分等问题，并提高基于GNSS坐标时间序列的衍生品的质量。突破了常用叠加算法对于影响因素在空间分布上的均匀性假设，也能够不受GNSS网尺寸大小限制，适用于任意尺寸的GNSS网测站坐标时间序列的区域划分。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例提供了一种基于主成分分析的区域划分方法，请参见图1，该方法包括：

步骤S1：获取GNSS测站的坐标时间序列观测值，并根据坐标时间序列观测值获得GNSS测站的残差坐标时间序列。

具体来说，PCA(Principal Component Analysis)是一种数据分析方法，它通过线性变换将原始数据变换为一组各维度线性无关的表示，用于提取数据的主要特征分量，可以用于高维数据的降维。一般数据的分布仅限于2维或者3维，低维数据可以在数值算法中明显减少运行时间，而且许多统计模型存在协变量的高度相关，PCA就可以产生不相关协变量的线性组合。

GNSS测站的坐标时间序列观测值可以通过原始观测数据处理获得，也可通过公开的产品获取；GNSS测站的残差坐标时间序列可以通过对坐标时间序列观测值进行最小二乘拟合后获取，也可通过预设的产品获取，其中，预设的产品为现有技术中公开的软件等，在此不做具体限定。在具体的实施过程中，可通过构建模型来建立坐标时间序列观测值和残差坐标时间序列之间的关系，具体参数包含长期趋势项、周年/半周年项等，然后采用最小二乘方法进行拟合。

步骤S2：计算GNSS网内所有测站的公共历元，并确定公共历元个数m，其中，m为大于0的整数。

具体来说，将每个测站的时间序列坐标时间序列观测值表示为位置随历元变化的函数，其中，一个历元对应一个位置；

根据上述函数，计算GNSS网内所有测站的公共历元，并确定公共历元的个数。

在具体的实施过程中，每个测站的时间序列均可表示为位置随时间(历元)的函数，其中，一个历元对应一个位置，如果对测站1有m个历元，测站2有k个历元，由于数据缺失、时间序列起止时间不同等原因，出现m与k不相等的情况，或者即使两者相等，代表的也不是同一跨度内具有共同观测时刻的时间序列。

步骤S3：基于公共历元构建GNSS网内所有测站在预设方向上的观测样本矩阵X(m,n)，其中，m表示公共历元个数，n表示GNSS网内测站的数目，n为大于0的整数，x_i,j为观测样本矩阵中的元素，表示第j个测站在第i个历元的误差观测值。

具体来说，预设方向可以是东方向，北方向和高程方向，其中方向用东西方向分量，简称为东方向；南北方向分量，简称为北方向，另一个为高程方向。

其中，步骤S3具体包括：

基于步骤S2中的公共历元，获得GNSS网内的任意一个测站均有m个历元的观测值，将其计作x，代表一个m行1列的向量，即x_m×1；

将满足上述的n个测站，表示为m行n列的向量，进而得到观测样本矩阵X(m,n)。

步骤S4：根据所有测站的残差坐标时间序列，构建所有测站在预设方向上的协方差矩阵B_n×n，其中，B_n×n中的对角元素用以表征测站的残差坐标时间序列的精度信息，非对角线元素用以表征相应测站之间的相关性。

具体地，通过下述运算计算获得所述协方差矩阵B_n×n，其中，m表示公共历元个数，X为观测样本矩阵，X^T为矩阵X的转置。

步骤S5：对协方差矩阵B_n×n进行正交分解，获得正交矩阵P。

其中，步骤S5具体包括：

通过如下运算B＝PΛP^T，获得正交矩阵P，其中，Λ为特征值构成的对角矩阵，矩阵P对应协方差矩阵B的主成分。

具体来说，本实施例中的正交分解是矩阵降维的重要方法，主要作用是为了获取步骤S6中的主成分。矩阵P是一个正交矩阵，对应的是协方差矩阵B的主成分，Λ是由特征值构成的对角矩阵。其中，上述计算所得的P矩阵是唯一的，其是中间计算的一个产物，与载荷矩阵一起构成对原始观测值的重构及解释。

步骤S6：根据观测样本矩阵和正交矩阵，获得载荷矩阵A，其中，载荷矩阵中的A_i',j'表示第i'个主成分在第j'个观测量上的空间响应，主成分为正交矩阵P的主成分分量，观测量为协方差矩阵B_n×n中的协方差构成的观测量。

步骤S6具体包括：

通过如下运算X_m×n＝A_m×nP_n×n，获得载荷矩阵A，其中，X_m×n为观测样本矩阵，P_n×n为正交矩阵。

具体来说，其中，矩阵A的元素A_i,j表示第i个主成分在第j个观测量上的空间响应，载荷系数A_i,j越大，响应的主成分对原观测量对解释程度越大，即空间响应越强。此步骤是通过对构建的协方差矩阵进行降维分解，因此观测量是协方差构成的观测量，主成分对应矩阵P矩阵中的分量。即矩阵P是主成分分量，矩阵A是载荷矩阵，反映的是对应主成分分量在空间上的响应。举例来说，1，2，3，4这四个测站通过上述的步骤，可以获得这4个测站降维以后的前四个主成分的载荷系数，对测站1可能为0.3，0.1，0.2，0.05，这种响应表明测站1中第一个主成分可以解释原观测量的30％，第二主成分10％，第三主成分20％，第四主成分5％等。

步骤S7：从正交矩阵P中获得特征矢量，并根据特征矢量的特征值的累加结果，选择出对应的K个累积特征矢量。

具体来说，可以将特征矢量按照特征值从大到小顺序排列，并根据特征值累计贡献率，得到累加结果，并根据累加结果选择出对应的K个累积特征矢量。

步骤S8：对K个累积特征矢量进行再次处理，对正交矩阵P中的m个变量的相关结构之间进行线性变换，并获得为线性变化后的载荷矩阵B’和线性变化后的主成分G。使得处理前后这K个空间型解释原场总方差对百分率保持不变。

具体来说，是对K个累积特征矢量进行进一步的处理，对正交矩阵P中变量进行线形变化。

步骤S8的具体实现包括如下变换：

X_m×n＝A_m×nP_n×n

X_m×n＝A_m×nRR^TP_n×n

X_m×n＝B'G,B'＝A_m×nR,G＝R^TP_n×n

其中，B'为线性变化后的载荷矩阵，G为线性变化后的主成分，R为变换矩阵，R^T为矩阵R的转置。

具体来说，K个累积特征矢量对应于对A矩阵进行了进一步的线性变换的矩阵中的分量。步骤S8的主要目的是使每个空间型(载荷向量)能反映原变量场的区域性相关结构，并且通过线性变换，使得载荷矩阵A的各列内部元素的平方之间的差异增大。首先根据观测样本矩阵X_m×n和正交矩阵P_n×n，得到载荷矩阵A_m×n，然后得到线形变化矩阵R和R^T。接下来，通过线形变化矩阵R与A矩阵相乘得到线形变化后的载荷矩阵B'，通过线形变化矩阵R^T与矩阵P相乘得到线性变化后的主成分G。

步骤S9：对线性变化后的载荷矩阵B’进行分析，获得区域划分结果。

具体来说，根据线性变化后的载荷矩阵B’，寻找在矩阵B’上对应的前k个主成分，以前k个主成分对应的载荷系数进行区域分块的划分。

划分的标准为，如果这k个主成分累计的载荷系数经过步骤S8中的线性变换后，对应的旋转主成分是不相关的，则将对应的测站分为一组。

相对于现有方法而言，本发明实施例提供的方法，通过构建协方差矩阵，并对协方差矩阵正交分解得到的主成分再次进行变换而得到载荷矩阵B’，即在现有方法的基础上进行了进一步的约束，然后通过载荷矩阵B’进行分析，获得区域划分结果的一种区域划分方法。通过对协方差矩阵进行降维处理，获取表征测站间关系的物理量，并实现降维，得到载荷矩阵，使得载荷矩阵中的载荷向量在不同的约束标准下能反映多个测站残差时间序列的区域性相关结构，从而达到准确地进行区域划分的技术效果。

本申请实施例中的上述一个或多个技术方案，至少具有如下一种或多种技术效果：

在本发明提供的方法中，首先获取GNSS测站的坐标时间序列观测值，并根据所述坐标时间序列观测值获得所述GNSS测站的残差坐标时间序列；并计算GNSS网内所有测站的公共历元，然后基于所述公共历元构建所述GNSS网内所有测站在预设方向上的观测样本矩阵X(m,n)，接下来根据所有测站的残差坐标时间序列，构建所述所有测站在预设方向上的协方差矩阵B_n×n，跟着对所述协方差矩阵B_n×n进行正交分解，获得正交矩阵P；然后根据所述观测样本矩阵和所述正交矩阵，获得载荷矩阵A，再从所述正交矩阵P中获得特征矢量，并根据特征矢量的特征值的累加结果，选择出对应的K个累积特征矢量；接着对所述K个累积特征矢量进行再次处理，对所述正交矩阵P中的m个变量的相关结构之间进行线性变换，并获得为线性变化后的载荷矩阵B’和线性变化后的主成分G。使得处理前后这K个空间型解释原场总方差对百分率保持不变；最后对所述线性变化后的载荷矩阵B’进行分析，获得区域划分结果。本发明提供的方案，通过协方差进行降维处理，获取表征各个测站间关系的物理量以实现降维，得到载荷矩阵，对特征矢量进行总方差不变情况下的线性变换，可以实现个空间型对原变量场区域结构对呈现。载荷矩阵B’的载荷向量在不同的约束标准下可以反映多个测站残差时间序列的区域性相关结构，从而可以进行准确地区域划分，适用于解决多个GNSS测站时间序列中出现共同误差及其区域划分等问题，并提高基于GNSS坐标时间序列的衍生品的质量。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明实施例进行各种改动和变型而不脱离本发明实施例的精神和范围。这样，倘若本发明实施例的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (9)

1.一种基于主成分分析的区域划分方法，其特征在于，包括：

步骤S1：获取GNSS测站的坐标时间序列观测值，并根据所述坐标时间序列观测值获得所述GNSS测站的残差坐标时间序列；

步骤S2：计算GNSS网内所有测站的公共历元，并确定所述公共历元个数m，其中，m为大于0的整数；

步骤S3：基于所述公共历元构建所述GNSS网内所有测站在预设方向上的观测样本矩阵X(m,n)，其中，m表示公共历元个数，n表示GNSS网内测站的数目，n为大于0的整数，x_i,j为观测样本矩阵中的元素，表示第j个测站在第i个历元的误差观测值；

步骤S4：根据所有测站的残差坐标时间序列，构建所述所有测站在预设方向上的协方差矩阵B_n×n，其中，B_n×n中的对角元素用以表征测站的残差坐标时间序列的精度信息，非对角线元素用以表征相应测站之间的相关性；

步骤S5：对所述协方差矩阵B_n×n进行正交分解，获得正交矩阵P；

步骤S6：根据所述观测样本矩阵和所述正交矩阵，获得载荷矩阵A，其中，载荷矩阵中的A_i',j'表示第i'个主成分在第j'个观测量上的空间响应，所述主成分为所述正交矩阵P的主成分分量，所述观测量由协方差矩阵B_n×n中的协方差构成；

步骤S7：从所述正交矩阵P中获得特征矢量，并根据特征矢量的特征值的累加结果，选择出对应的K个累积特征矢量；

步骤S8：对所述K个累积特征矢量进行再次处理，对所述正交矩阵P中的m个变量的相关结构之间进行线性变换，并获得为线性变化后的载荷矩阵B’和线性变化后的主成分G，使得处理前后这K个空间型解释原场总方差对百分率保持不变；

步骤S9：对所述线性变化后的载荷矩阵B’进行分析，获得区域划分结果。

2.如权利要求1所述的方法，其特征在于，步骤S1中，根据所述坐标时间序列观测值获得所述GNSS测站的残差坐标时间序列，包括：

通过对坐标时间序列观测值进行最小二乘拟合，获得所述GNSS测站的残差坐标时间序列；或

通过预设产品获得所述GNSS测站的残差坐标时间序列。

3.如权利要求1所述的方法，其特征在于，步骤S2具体包括：

将每个测站的时间序列坐标时间序列观测值表示为位置随历元变化的函数，其中，一个历元对应一个位置；

根据上述函数，计算GNSS网内所有测站的公共历元，并确定所述公共历元的个数。

4.如权利要求1所述的方法，其特征在于，步骤S3具体包括：

基于步骤S2中的公共历元，获得GNSS网内的任意一个测站均有m个历元的观测值，将其计作x，代表一个m行1列的向量，即x_m×1；

将满足上述的n个测站，表示为m行n列的向量，进而得到观测样本矩阵X(m,n)。

5.如权利要求1所述的方法，其特征在于，步骤S4具体包括：

通过下述运算计算获得所述协方差矩阵B_n×n，其中，m表示公共历元个数，X为观测样本矩阵，X^T为矩阵X的转置。

6.如权利要求1所述的方法，其特征在于，步骤S5具体包括：

通过如下运算B＝PΛP^T，获得正交矩阵P，其中，Λ为特征值构成的对角矩阵，矩阵P对应协方差矩阵B的主成分。

7.如权利要求1所述的方法，其特征在于，步骤S6具体包括：

通过如下运算X_m×n＝A_m×nP_n×n，获得载荷矩阵A，其中，X_m×n为观测样本矩阵，P_n×n为正交矩阵。

8.如权利要求1所述的方法，其特征在于，步骤S8的具体实现包括如下变换：

X_m×n＝A_m×nP_n×n

X_m×n＝A_m×nRR^TP_n×n

X_m×n＝B'G,B'＝A_m×nR,G＝R^TP_n×n

其中，B'为线性变化后的载荷矩阵，G为线性变化后的主成分，R为变换矩阵，R^T为矩阵R的转置。

9.如权利要求8所述的方法，其特征在于，步骤S9具体包括：

根据线性变化后的载荷矩阵B’，寻找在矩阵B’上对应的前k个主成分，以所述前k个主成分对应的载荷系数进行区域分块的划分。