CN104776827A - Gps高程异常数据的粗差探测方法 - Google Patents
Gps高程异常数据的粗差探测方法 Download PDFInfo
- Publication number
- CN104776827A CN104776827A CN201510159639.9A CN201510159639A CN104776827A CN 104776827 A CN104776827 A CN 104776827A CN 201510159639 A CN201510159639 A CN 201510159639A CN 104776827 A CN104776827 A CN 104776827A
- Authority
- CN
- China
- Prior art keywords
- gross
- error
- data
- distance
- point
- 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.)
- Granted
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 21
- 238000000034 method Methods 0.000 claims abstract description 19
- 239000011159 matrix material Substances 0.000 claims description 9
- 230000002159 abnormal effect Effects 0.000 claims description 5
- 238000012937 correction Methods 0.000 claims description 2
- 238000005259 measurement Methods 0.000 abstract description 13
- 238000004458 analytical method Methods 0.000 abstract description 2
- 230000002547 anomalous effect Effects 0.000 abstract 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009897 systematic effect Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开一种GPS高程异常数据的粗差探测方法是基于异常点集多维粗差定位方法,具体为:1)数据采集;2)利用数据建立模型;3)计算Cook距离,根据各观测值的Cook距离D(i)和平均值D平均之差,确定强影响点号4)计算wk距离并依据分位值确定强影响点号;5)将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集,得异常点集。令不重复的强影响点号个数m为粗差的维数;6)做粗差的定位。经过大量实例应用结果分析,在粗差维数达到总样本的19%时,本发明仍能将粗差准确的定位。本发明对提升测量数据的应用质量具有重要意义。具有明显的社会和经济价值。
Description
技术领域
本发明涉及一种GPS高程异常数据的粗差探测方法,属于“测绘科学与技术”学科中的“大地测量学”技术领域。
背景技术
数据是进行科学研究的基础,为此人们运用各种手段获得测量或实验数据。由于种种原因,测量数据与被测原型的真值不完全一致,其中包括粗差、系统误差和偶然误差等测量误差。在应用测量数据时,要求数据仅含有偶然误差。系统误差可以采用一定的方法发现和修正。所谓检测粗差就是检测数据是否含有粗差,以及定位哪一个数据是粗差。如果不能有效地对粗差进行检验、识别,将会对参数的估计结果产生严重影响。检测粗差是有效的应用测量数据的前提。
对于粗差检测,先后进行了一维和多维的粗差检测研究。目前有以服从正态分布的标准残差为统计量的数据检测法;粗差判断的Bayes法;基于岭估计的粗差检测方法;粗差的拟准检定法;基于方差膨胀模型的多个粗差的检测法;应用偏相关系数区分多维粗差和用复相关系数对多维粗差进行总体显著性检测并定位等。综合来看,多维粗差检测和定位比较困难,主要是难以确定粗差的维数。如果用试探法同时确定多个粗差,试探的组合数可能相当大,也难以实现;一维粗差检测法虽简单可行,但未能考虑粗差对残差影响的相关性,用它逐次判断多个粗差的存在,也可能判断失误。
针对以上问题,本发明提出利用“奇异点”集,采用一定的工作流程来准确进行多维粗差检测及定位的方法。本发明的实质是将不同方法得到的强影响点的并集作为“奇异点”集,可以大大提高检测测量数据粗差的准确性。本发明对提升测量数据的应用质量具有重要意义。
发明内容
发明目的:本发明的目的在于本发明的目的是提供一种检测多维粗差的方法,使用该方法,能准确检测出测量数据中的粗差,使用方便。
技术方案:本发明所述的检测多维粗差的方法为:
1)利用观测值建立如下模型:
函数模型:L=AX+Δ (1)
随机模型:E(Δ)=0,
式中:为观测值向量;为未知参数向量;为未知参数的系数矩阵,rank(A)=p;为随机误差向量;E(·)为数学期望;D(·)为方差—协方差矩阵;为单位权方差;QLL为观测值协因数阵;P为权阵。
a)参数估计
X的估值用表示,并令N=ATPA。当rank(A)=p为列满秩阵时,得未知参数X最小二乘估计
误差向量的最小二乘估计
H=A(ATPA)-1ATP (5)
称为帽子矩阵。
残差的协因数矩阵
QVV=(I-H)QLL (6)
b)单位权方差的估值
2)计算Cook距离,根据各观测值的Cook距离D(i)和平均值D平均之差,确定强影响点号
式中:为从完全数据回归未知X的最小二乘估计;为删掉第i组数据后;从剩余的n-1组数据回归X的最小二乘估计。
当D(i)-D平均>0时确定第i号点为强影响点。
3)计算wk距离并依据分位值确定强影响点号
用式(9)计算Welsch—Kuh距离
当
时,第i个点为强影响点。
时,第i个点为强影响点。式中vi为(4)式中的第i(i=1,2,…n)个观测值的改正数,为去掉第i个观测值的σ0的估值,为(6)式中第i个主对角元,hii为(5)式中的第i个主对角元。为服从自由度为(n-t-1)的t分布值。
4)将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集,得异常点集。令不重复的强影响点号个数m为粗差的维数;
5)令
Qee=(BTPQVVPB)-1 (10)
则
Xe=-QeeBTPV (11)
式中B=(E1,E2,…,Em)为n×m矩阵,Ei为n×1向量,其中第i行元素为1,其余元素为零。
6)做粗差的定位。
令
当Fi>Fα,即第i(i=1,2,…n)号点为粗差点,粗差点个数为n1(n1=0或n1≤m);Fα为服从自由度为(1,n-t-m)的F分布值。
8)若第6)步中n1≠0,剔除n1个粗差数据,得n2=n-n1组观测数据,循环步骤1)—6)直至没有粗差点。若第6)步中n1=0,结束。
本发明与现有技术相比,其有益效果是:1)粗差定位准确高效。经过大量实例应用结果分析,在粗差维数达到总样本的19%时,本发明仍能将粗差准确的定位。
2)使用方便。整个多维粗差定位过程,无需编制专门的程序,仅需利用国际公认的标准计算软件MATLAB,以及t分布和F分布表即可完成粗差定位,快速、方便。
具体实施方式
下面对本发明技术方案进行详细说明,但是本发明的保护范围不局限于所述实施例。
实施例1:
在描述具体实施过程时,结合某具体工程实例,对本发明方法作详细说明。
(1)工程实例背景介绍
建立区域似大地水准面,需布设测量控制点。为获取有关数据,要求对测量控制点进行两项测量工作,GPS测量和水准测量。
①GPS测量:按照国家测量规范要求,对所有控制点进行相应等级的GPS测量,目的是获取各控制点的高斯平面坐标和大地高。
②水准测量:按照国家测量规范要求,对所有控制点进行相应等级的水准测量,目的是获取各控制点的正常高高程,并计算各控制点的高程异常。
具体实施例,建立区域面积约为220km2的似大地水准面,共布设36个测量控制点。按照国家测量规范要求,对所有控制点进行了C级GPS测量和三等水准测量,获得不含粗差的GPS高程异常数据如表1所示:
表1 GPS高程异常数据表
现对5、10、15、20、25、30及35分别施加3.5~10倍中误差的粗差。得含有粗差的数据如表2:
表2含有粗差的GPS高程异常数据表
(加*者为粗差点)
(2)利用观测值建立模型
建立二次曲面模型
ξ=b0+b1x+b2y+b3x2+b4xy+b5 2 (15)
式中,ξ为GPS点的高程异常;bi为二次曲面模型的待定参数(共6个);x、y为GPS点的高斯平面坐标。
具体实施例,将表2的36组GPS高程异常数据逐个代入(15)式,可以得到形如(1)的式矩阵形式。然后根据式(3)、(4)、(5)、(6)、(7)计算V、H、QVV及
(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,
式中:为观测值向量;为未知参数向量;为未知参数的系数矩阵,rank(A)=p;为随机误差向量;E(·)为数学期望;D(·)为方差—协方差矩阵;为单位权方差;QLL为观测值协因数阵;P为权阵。
2)参数估计
X的估值用表示,并令N=ATPA;当rank(A)=p为列满秩阵时,计算:
未知参数X最小二乘估计
误差向量的最小二乘估计
帽子矩阵
H=A(ATPA)-1ATP (5)
残差的协因数矩阵
QVV=(I-H)QLL (6)
单位权方差的估值
3)根据(8)式计算每个观测值的Cook距离,根据各观测值的Cook距离D(i)和平均值D平均之差,确定强影响点号
式中:为从完全数据回归未知X的最小二乘估计;为删掉第i组数据后;从剩余的n-1组数据回归X的最小二乘估计;
当D(i)-D平均>0时确定第i号点为强影响点。
4)计算wk距离并依据分位值确定强影响点号
用式(9)计算Welsch—Kuh距离
5)将Cook距离及Welsch—Kuh距离确定的强影响点号组成并集,得异常点集;令不重复的强影响点号个数m为粗差的维数;
6)令
Qee=(BTPQVVPB)-1 (10)
则
Xe=-QeeBTPV (11)
式中B=(E1,E2,…,Em)为n×m矩阵,Ei为n×1向量,其中第i行元素为1,其余元素为零;
7)粗差定位
令
当Fi>Fα,即第i(i=1,2,…n)号点为粗差点,粗差点个数为n1(n1=0或n1≤m);Fα为服从自由度为(1,n-t-m)的F分布值。
8)若第7)步中n1≠0,剔除n1个粗差数据,得n2=n-n1组观测数据,循环步骤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高程异常数据。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510159639.9A CN104776827B (zh) | 2015-04-03 | 2015-04-03 | Gps高程异常数据的粗差探测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510159639.9A CN104776827B (zh) | 2015-04-03 | 2015-04-03 | Gps高程异常数据的粗差探测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104776827A true CN104776827A (zh) | 2015-07-15 |
CN104776827B CN104776827B (zh) | 2017-04-05 |
Family
ID=53618442
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510159639.9A Expired - Fee Related CN104776827B (zh) | 2015-04-03 | 2015-04-03 | Gps高程异常数据的粗差探测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104776827B (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105953772A (zh) * | 2016-04-22 | 2016-09-21 | 深圳市华信天线技术有限公司 | 一种高程异常的测定方法及测定装置 |
CN108564298A (zh) * | 2018-04-26 | 2018-09-21 | 三重能有限公司 | 数据处理方法、装置及电子设备 |
CN109270560A (zh) * | 2018-10-12 | 2019-01-25 | 东南大学 | 区域高程异常数据的多维粗差定位及定值方法 |
CN110596736A (zh) * | 2019-10-15 | 2019-12-20 | 中国电子科技集团公司第五十四研究所 | 一种gnss观测异常值探测与隔离方法 |
CN114418886A (zh) * | 2022-01-19 | 2022-04-29 | 电子科技大学 | 一种基于深度卷积自编码器的鲁棒性去噪方法 |
CN114964152A (zh) * | 2022-05-25 | 2022-08-30 | 上海勘察设计研究院(集团)有限公司 | 一种基于多维特征的实时水准测量点号匹配方法 |
CN117473717A (zh) * | 2023-10-19 | 2024-01-30 | 同济大学 | 一种基于Bernoulli-Gaussian模型和EM算法的数据质量分析方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030171856A1 (en) * | 2002-03-07 | 2003-09-11 | Wilf Herbert S. | Global positioning system readout of recommended altitude in aircraft landing pattern |
CN101839710A (zh) * | 2010-06-12 | 2010-09-22 | 中国测绘科学研究院 | 一种似大地水准面计算的优化方法 |
CN103698790A (zh) * | 2013-12-30 | 2014-04-02 | 辽宁工程技术大学 | 北斗与gps双系统宽巷载波相位混频星间差分组合方法 |
CN104111061A (zh) * | 2014-06-30 | 2014-10-22 | 中国电力工程顾问集团中南电力设计院 | 一种基础资料缺乏地区测量控制点成果获取方法 |
CN104156968A (zh) * | 2014-08-19 | 2014-11-19 | 山东临沂烟草有限公司 | 大面积地形复杂区域无人机序列影像快速无缝拼接方法 |
-
2015
- 2015-04-03 CN CN201510159639.9A patent/CN104776827B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030171856A1 (en) * | 2002-03-07 | 2003-09-11 | Wilf Herbert S. | Global positioning system readout of recommended altitude in aircraft landing pattern |
CN101839710A (zh) * | 2010-06-12 | 2010-09-22 | 中国测绘科学研究院 | 一种似大地水准面计算的优化方法 |
CN103698790A (zh) * | 2013-12-30 | 2014-04-02 | 辽宁工程技术大学 | 北斗与gps双系统宽巷载波相位混频星间差分组合方法 |
CN104111061A (zh) * | 2014-06-30 | 2014-10-22 | 中国电力工程顾问集团中南电力设计院 | 一种基础资料缺乏地区测量控制点成果获取方法 |
CN104156968A (zh) * | 2014-08-19 | 2014-11-19 | 山东临沂烟草有限公司 | 大面积地形复杂区域无人机序列影像快速无缝拼接方法 |
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105953772A (zh) * | 2016-04-22 | 2016-09-21 | 深圳市华信天线技术有限公司 | 一种高程异常的测定方法及测定装置 |
CN108564298A (zh) * | 2018-04-26 | 2018-09-21 | 三重能有限公司 | 数据处理方法、装置及电子设备 |
CN108564298B (zh) * | 2018-04-26 | 2021-06-22 | 三一重能股份有限公司 | 数据处理方法、装置及电子设备 |
CN109270560A (zh) * | 2018-10-12 | 2019-01-25 | 东南大学 | 区域高程异常数据的多维粗差定位及定值方法 |
CN110596736A (zh) * | 2019-10-15 | 2019-12-20 | 中国电子科技集团公司第五十四研究所 | 一种gnss观测异常值探测与隔离方法 |
CN110596736B (zh) * | 2019-10-15 | 2021-04-02 | 中国电子科技集团公司第五十四研究所 | 一种gnss观测异常值探测与隔离方法 |
CN114418886A (zh) * | 2022-01-19 | 2022-04-29 | 电子科技大学 | 一种基于深度卷积自编码器的鲁棒性去噪方法 |
CN114964152A (zh) * | 2022-05-25 | 2022-08-30 | 上海勘察设计研究院(集团)有限公司 | 一种基于多维特征的实时水准测量点号匹配方法 |
CN114964152B (zh) * | 2022-05-25 | 2024-03-15 | 上海勘察设计研究院(集团)股份有限公司 | 一种基于多维特征的实时水准测量点号匹配方法 |
CN117473717A (zh) * | 2023-10-19 | 2024-01-30 | 同济大学 | 一种基于Bernoulli-Gaussian模型和EM算法的数据质量分析方法 |
Also Published As
Publication number | Publication date |
---|---|
CN104776827B (zh) | 2017-04-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104776827A (zh) | Gps高程异常数据的粗差探测方法 | |
CN102495413B (zh) | 输电线路杆塔坐标的获得方法 | |
CN105046046B (zh) | 一种集合卡尔曼滤波局地化方法 | |
CN109188481A (zh) | 一种拟合gps高程异常的新方法 | |
CN103278126A (zh) | 一种基于最小区域的零件球度误差评定方法 | |
CN103526783A (zh) | 一种测量建筑基坑水平位移的方法 | |
CN110333543A (zh) | 基于反射系数分析的低阻体解释及成像方法与系统 | |
CN103292773A (zh) | 一种基于最小区域的对称度误差评定方法 | |
CN111812730B (zh) | 一种用于滑坡探测的电阻率数据融合三维成像方法及系统 | |
CN111239838B (zh) | 一种磁探测精度的检测方法 | |
CN105651311A (zh) | 农机作业卫星导航自动驾驶精度的测试方法 | |
CN104180822B (zh) | 一种变形监测基准点稳定性检验方法 | |
CN101266153B (zh) | 测绘工程类陀螺全站仪精度评定方法 | |
CN105043390A (zh) | 基于泛克里金法的重力场插值方法 | |
CN107607904A (zh) | 一种基于测角信息的多目标同一性识别方法 | |
CN109683203A (zh) | 一种地震地电观测中偶极干扰源定位方法及系统 | |
CN104864906A (zh) | 海洋油气水下设备重量测量和重心检测方法 | |
CN104751005B (zh) | 一种基于正交实验的平面度误差评定方法 | |
CN114236624B (zh) | 基于电磁法估算压裂改造空间体积的方法和系统 | |
CN105260610A (zh) | 一种多探测器坐标系转化及误差纠正方法 | |
CN109270560B (zh) | 区域高程异常数据的多维粗差定位及定值方法 | |
CN106595472B (zh) | 摄影测量系统的精度确定方法 | |
CN111681124B (zh) | 一种深部砂岩型铀矿化信息三维氡异常识别方法及系统 | |
CN111597752B (zh) | 平衡孔间敏感性的跨孔电阻率ct深度学习反演方法及系统 | |
Hao et al. | Kriging-based study on the visualization of magnetic method data |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
EXSB | Decision made by sipo to initiate substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170405 |