CN104776827A - Gps高程异常数据的粗差探测方法 - Google Patents

Abstract

本发明公开一种GPS高程异常数据的粗差探测方法是基于异常点集多维粗差定位方法，具体为：1)数据采集；2)利用数据建立模型；3)计算Cook距离，根据各观测值的Cook距离D(i)和平均值D_平均之差，确定强影响点号4)计算wk距离并依据分位值确定强影响点号；5)将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集，得异常点集。令不重复的强影响点号个数m为粗差的维数；6)做粗差的定位。经过大量实例应用结果分析，在粗差维数达到总样本的19％时，本发明仍能将粗差准确的定位。本发明对提升测量数据的应用质量具有重要意义。具有明显的社会和经济价值。

Description

GPS高程异常数据的粗差探测方法

技术领域

本发明涉及一种GPS高程异常数据的粗差探测方法，属于“测绘科学与技术”学科中的“大地测量学”技术领域。

背景技术

数据是进行科学研究的基础，为此人们运用各种手段获得测量或实验数据。由于种种原因，测量数据与被测原型的真值不完全一致，其中包括粗差、系统误差和偶然误差等测量误差。在应用测量数据时，要求数据仅含有偶然误差。系统误差可以采用一定的方法发现和修正。所谓检测粗差就是检测数据是否含有粗差，以及定位哪一个数据是粗差。如果不能有效地对粗差进行检验、识别，将会对参数的估计结果产生严重影响。检测粗差是有效的应用测量数据的前提。

对于粗差检测,先后进行了一维和多维的粗差检测研究。目前有以服从正态分布的标准残差为统计量的数据检测法；粗差判断的Bayes法；基于岭估计的粗差检测方法；粗差的拟准检定法；基于方差膨胀模型的多个粗差的检测法；应用偏相关系数区分多维粗差和用复相关系数对多维粗差进行总体显著性检测并定位等。综合来看，多维粗差检测和定位比较困难，主要是难以确定粗差的维数。如果用试探法同时确定多个粗差，试探的组合数可能相当大，也难以实现；一维粗差检测法虽简单可行，但未能考虑粗差对残差影响的相关性，用它逐次判断多个粗差的存在，也可能判断失误。

针对以上问题，本发明提出利用“奇异点”集，采用一定的工作流程来准确进行多维粗差检测及定位的方法。本发明的实质是将不同方法得到的强影响点的并集作为“奇异点”集，可以大大提高检测测量数据粗差的准确性。本发明对提升测量数据的应用质量具有重要意义。

发明内容

发明目的：本发明的目的在于本发明的目的是提供一种检测多维粗差的方法，使用该方法，能准确检测出测量数据中的粗差，使用方便。

技术方案：本发明所述的检测多维粗差的方法为：

1)利用观测值建立如下模型：

函数模型：L＝AX+Δ    (1)

随机模型：E(Δ)＝0, $D\left(L\right)=D\left(\mathrm{\Δ}\right)={\mathrm{\σ}}_0^2{Q}_{\mathrm{LL}}={\mathrm{\σ}}_0^2{P}^{-1}---\left(2\right)$

式中:为观测值向量；为未知参数向量；为未知参数的系数矩阵，rank(A)＝p；为随机误差向量；E(·)为数学期望；D(·)为方差—协方差矩阵；为单位权方差；Q_LL为观测值协因数阵；P为权阵。

a)参数估计

X的估值用表示，并令N＝A^TPA。当rank(A)＝p为列满秩阵时，得未知参数X最小二乘估计

$\hat{X}=N^{-1}{A}^T\mathrm{PL}={\left(A^T\mathrm{PA}\right)}^{-1}{A}^T\mathrm{PL}---\left(3\right)$

误差向量的最小二乘估计

$V=A\hat{X}-L=(A{\left(A^T\mathrm{PA}\right)}^{-1}{A}^{T}P-I)L=(H-I)L---\left(4\right)$

H＝A(A^TPA)^-1A^TP    (5)

称为帽子矩阵。

残差的协因数矩阵

Q_VV＝(I-H)Q_LL    (6)

b)单位权方差的估值

${\widehat{\mathrm{\σ}}}_0^2=V^T\mathrm{PV}/r=V^T\mathrm{PV}/n-p---\left(7\right)$

2)计算Cook距离，根据各观测值的Cook距离D(i)和平均值D_平均之差，确定强影响点号

$D\left(i\right)=\frac{{(\hat{X}\left(i\right)-\hat{X})}^T{A}^T\mathrm{PA}(\hat{X}\left(i\right)-\hat{X})}{p\widehat{{\mathrm{\σ}}_0^2}}---\left(8\right)$

式中：为从完全数据回归未知X的最小二乘估计；为删掉第i组数据后；从剩余的n－1组数据回归X的最小二乘估计。

当D(i)－D_平均＞0时确定第i号点为强影响点。

3)计算wk距离并依据分位值确定强影响点号

用式(9)计算Welsch—Kuh距离

${\mathrm{wk}}_i=\left|\frac{v_i}{\widehat{{\mathrm{\σ}}_0}\left(i\right)\sqrt{Q_{\mathrm{vivi}}}}\right|{\left(\frac{h_{\mathrm{ii}}}{1-h_{\mathrm{ii}}}\right)}^{1/2}---\left(9\right)$

当

时，第i个点为强影响点。

时，第i个点为强影响点。式中v_i为(4)式中的第i(i＝1,2,…n)个观测值的改正数，为去掉第i个观测值的σ₀的估值，为(6)式中第i个主对角元，h_ii为(5)式中的第i个主对角元。为服从自由度为(n－t－1)的t分布值。

4)将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集，得异常点集。令不重复的强影响点号个数m为粗差的维数；

5)令

Q_ee＝(B^TPQ_VVPB)^-1    (10)

则

X_e＝-Q_eeB^TPV    (11)

${\hat{X}}^{*}={\left(A^T\mathrm{PA}\right)}^{-1}{A}^T\mathrm{PL}+Q_{\hat{x}\hat{x}}{A}^T\mathrm{PB}{Q}_{\mathrm{ee}}{B}^T\mathrm{PV}---\left(12\right)$

${\mathrm{\Ω}}_e={(A{\hat{X}}^{*}+B{X}_e-L)}^{T}P(A{\hat{X}}^{*}+B{X}_e-L)---\left(13\right)$

式中B＝(E₁，E₂，…，E_m)为n×m矩阵，E_i为n×1向量，其中第i行元素为1，其余元素为零。

6)做粗差的定位。

令

$F_i=\frac{{{\left(X_e\right)}_i}^2{\left(Q_{\mathrm{ee}}\right)}_i^{-1}}{{\mathrm{\Ω}}_e/(n-t-m)}~F_{(\mathrm{\α},1,n-t-m)}---\left(14\right)$

当F_i＞F_α，即第i(i＝1,2,…n)号点为粗差点，粗差点个数为n₁(n1＝0或n1≤m)；F_α为服从自由度为(1，n－t－m)的F分布值。

8)若第6)步中n1≠0，剔除n₁个粗差数据，得n₂＝n-n₁组观测数据,循环步骤1)—6)直至没有粗差点。若第6)步中n1＝0，结束。

本发明与现有技术相比，其有益效果是：1)粗差定位准确高效。经过大量实例应用结果分析，在粗差维数达到总样本的19％时，本发明仍能将粗差准确的定位。

2)使用方便。整个多维粗差定位过程，无需编制专门的程序，仅需利用国际公认的标准计算软件MATLAB，以及t分布和F分布表即可完成粗差定位，快速、方便。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

实施例1：

在描述具体实施过程时，结合某具体工程实例，对本发明方法作详细说明。

(1)工程实例背景介绍

建立区域似大地水准面，需布设测量控制点。为获取有关数据，要求对测量控制点进行两项测量工作，GPS测量和水准测量。

①GPS测量：按照国家测量规范要求，对所有控制点进行相应等级的GPS测量，目的是获取各控制点的高斯平面坐标和大地高。

②水准测量：按照国家测量规范要求，对所有控制点进行相应等级的水准测量，目的是获取各控制点的正常高高程，并计算各控制点的高程异常。

具体实施例，建立区域面积约为220km²的似大地水准面，共布设36个测量控制点。按照国家测量规范要求，对所有控制点进行了C级GPS测量和三等水准测量，获得不含粗差的GPS高程异常数据如表1所示：

表1 GPS高程异常数据表

现对5、10、15、20、25、30及35分别施加3.5～10倍中误差的粗差。得含有粗差的数据如表2：

表2含有粗差的GPS高程异常数据表

(加*者为粗差点)

(2)利用观测值建立模型

建立二次曲面模型

ξ＝b₀+b₁x+b₂y+b₃x²+b₄xy+b₅ ²    (15)

式中，ξ为GPS点的高程异常；b_i为二次曲面模型的待定参数(共6个)；x、y为GPS点的高斯平面坐标。

具体实施例，将表2的36组GPS高程异常数据逐个代入(15)式，可以得到形如(1)的式矩阵形式。然后根据式(3)、(4)、(5)、(6)、(7)计算V、H、Q_VV及

(3)根据(8)式计算Cook距离见表3，且D_平均＝0.0319。经过计算1、4、5、8、10、15、16、20、25、30及35号点的D(i)－D_平均＞0，故其为强影响点。

(4)用式(9)计算Welsch—Kuh距离见表3，并依据分位值确定强影响点号。由表3知25及35号点的Welsch—Kuh距离(加*号者)大于其分位置(α＝0.05)，故25及35号点为强影响点。

表3 Cook距离及Welsch—Kuh距离(wk)表

(5)由步骤(3)及(4)的并集得异常点集为1、4、5、8、10、15、16、20、25、30及35号点。m＝11。

(6)根据式(10)—式(14)进行粗差定位。其中t＝6,α＝0.05及

定位结果见表4，7个粗差点全部定位正确。

表4粗差定位结果表

如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。

Claims (4)

1.一种GPS高程异常数据的粗差探测方法，其特征在于，包括以下步骤：

1)获取GPS高程异常数据作为观测值建立如下模型：

函数模型：L＝AX+Δ           (1)

随机模型：E(Δ)＝0, $D\left(L\right)=D\left(\mathrm{\Δ}\right)={\mathrm{\σ}}_0^2{Q}_{\mathrm{LL}}={\mathrm{\σ}}_0^2{P}^{-1}---\left(2\right)$

式中:为观测值向量；为未知参数向量；为未知参数的系数矩阵，rank(A)＝p；为随机误差向量；E(·)为数学期望；D(·)为方差—协方差矩阵；为单位权方差；Q_LL为观测值协因数阵；P为权阵。

2)参数估计

X的估值用表示，并令N＝A^TPA；当rank(A)＝p为列满秩阵时，计算：

未知参数X最小二乘估计

$\hat{X}=N^{-1}{A}^T\mathrm{PL}={\left(A^T\mathrm{PA}\right)}^{-1}{A}^T\mathrm{PL}---\left(3\right)$

误差向量的最小二乘估计

$V=A\hat{X}-L=(A{\left(A^T\mathrm{PA}\right)}^{-1}{A}^{T}P-I)L=(H-I)L---\left(4\right)$

帽子矩阵

H＝A(A^TPA)^-1A^TP          (5)

残差的协因数矩阵

Q_VV＝(I-H)Q_LL              (6)

单位权方差的估值

${\widehat{\mathrm{\σ}}}_0^2=V^T\mathrm{PV}/r=V^T\mathrm{PV}/n-p---\left(7\right)$

3)根据(8)式计算每个观测值的Cook距离，根据各观测值的Cook距离D(i)和平均值D_平均之差，确定强影响点号

$D\left(i\right)=\frac{{(\hat{X}\left(i\right)-\hat{X})}^T{A}^T\mathrm{PA}(\hat{X}\left(i\right)-\hat{X})}{p\widehat{{\mathrm{\σ}}_0^2}}---\left(8\right)$

式中：为从完全数据回归未知X的最小二乘估计；为删掉第i组数据后；从剩余的n－1组数据回归X的最小二乘估计；

当D(i)－D_平均＞0时确定第i号点为强影响点。

4)计算wk距离并依据分位值确定强影响点号

用式(9)计算Welsch—Kuh距离

${\mathrm{wk}}_i=\left|\frac{v_i}{\widehat{{\mathrm{\σ}}_0}\left(i\right)\sqrt{Q_{\mathrm{vivi}}}}\right|{\left(\frac{h_{\mathrm{ii}}}{1-h_{\mathrm{ii}}}\right)}^{1/2}---\left(9\right)$ 当时，第i个点为强影响点。式中v_i为(4)式中的第i(i＝1,2,…n)个观测值的改正数，去掉第i个观测值的σ₀的估值，为(6)式中第i个主对角元，h_ii为(5)式中的第i个主对角元。为服从自由度为(n－t－1)的t分布值。

5)将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集，得异常点集；令不重复的强影响点号个数m为粗差的维数；

6)令

Q_ee＝(B^TPQ_VVPB)^-1         (10)

则

X_e＝-Q_eeB^TPV           (11)

${\hat{X}}^{*}={\left(A^T\mathrm{PA}\right)}^{-1}{A}^T\mathrm{PL}+Q_{\hat{X}\hat{X}}{A}^T{\mathrm{PBQ}}_{\mathrm{ee}}{B}^T\mathrm{PV}---\left(12\right)$

${\mathrm{\Ω}}_e={(A{\hat{X}}^{*}+{\mathrm{BX}}_e-L)}^{T}P(A{\hat{X}}^{*}+{\mathrm{BX}}_e-L)---\left(13\right)$

式中B＝(E1，E2，…，Em)为n×m矩阵，Ei为n×1向量，其中第i行元素为1，其余元素为零；

7)粗差定位

令 $F_i=\frac{{{\left(X_e\right)}_i}^2{\left(Q_{\mathrm{ee}}\right)}_i^{-1}}{{\mathrm{\Ω}}_e/(n-t-m)}~F_{(\mathrm{\α},1,n-t-m)}---\left(14\right)$

当F_i＞F_α，即第i(i＝1,2,…n)号点为粗差点，粗差点个数为n₁(n1＝0或n1≤m)；F_α为服从自由度为(1，n－t－m)的F分布值。

8)若第7)步中n1≠0，剔除n₁个粗差数据，得n₂＝n-n₁组观测数据,循环步骤1)—7)直至没有粗差点。若第7)步中n1＝0，结束。

2.根据权利要求1所述的GPS高程异常数据的粗差探测方法，其特征在于，所述的GPS高程数据内的粗差数据维数小于总样本数的19％。

3.根据权利要求1所述的GPS高程异常数据的粗差探测方法，其特征在于，计算Cook距离，当D(i)－D_平均＞0时确定第i号点为强影响点；计算wk距离并依据分位值确定强影响点号；将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集，得异常点集；令不重复的异常点号个数m为粗差的维数。

4.根据权利要求1所述的GPS高程异常数据的粗差探测方法，其特征在于，所述的GPS高程异常数据内的粗差数据均为统一为国际标准单位制基本单位的GPS高程异常数据。