CN112131752B - 一种基于拟准检定的超强崩溃污染率抗差估计算法 - Google Patents

一种基于拟准检定的超强崩溃污染率抗差估计算法 Download PDF

Info

Publication number
CN112131752B
CN112131752B CN202011047884.8A CN202011047884A CN112131752B CN 112131752 B CN112131752 B CN 112131752B CN 202011047884 A CN202011047884 A CN 202011047884A CN 112131752 B CN112131752 B CN 112131752B
Authority
CN
China
Prior art keywords
quasi
value
observation
robust
error
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.)
Active
Application number
CN202011047884.8A
Other languages
English (en)
Other versions
CN112131752A (zh
Inventor
瞿伟
陈海禄
张勤
高源�
梁世川
韩亚茜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changan University
Original Assignee
Changan University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Changan University filed Critical Changan University
Priority to CN202011047884.8A priority Critical patent/CN112131752B/zh
Publication of CN112131752A publication Critical patent/CN112131752A/zh
Application granted granted Critical
Publication of CN112131752B publication Critical patent/CN112131752B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于拟准检定的超强崩溃污染率抗差估计算法，利用K均值聚类算法实现了任意粗差占比情况下对于拟准观测值的自动化选择，以实现粗差的粗识别，然后以拟准真误差作为抗差估计中等价权函数的初值进行迭代计算，以实现粗差的精识别和模型参数的超强崩溃污染抗差估计。本方法相较于常规抗差估计和基于残差中位数的抗差估计而言，可更加准确地实现对区域GNSS速度场中粗差数据的探测，实现区域地壳运动模型参数的超强崩溃污染率抗差估计，为后续进一步研究区域地壳形变特征提供更真实、更有价值的基础数据，为复杂场景中地壳形变监测数据的粗差探测和模型参数估计提供一种有效的处理方法。

Description

一种基于拟准检定的超强崩溃污染率抗差估计算法
技术领域
本发明属于高精度地壳形变监测数据处理领域，涉及到一种基于拟准检定的地壳运动GNSS水平速度场粗差探测与模型参数估计技术，该算法以高精度地壳形变监测为实际应用背景，可用于观测环境复杂、局部地壳活动活跃等场景中的地壳形变高精度监测应用方向。
背景技术
随着现代空间大地测量技术的迅猛发展，特别是以全球导航卫星系统(GNSS)为代表的空间监测技术的现代化建设，可利用其以厘米级甚至毫米级的精度实现对地壳水平运动的高精度监测。但受监测环境干扰、局部地壳活动或地震等因素影响，导致所观测到的GNSS速度场数据中通常含有粗差，若粗差无法被有效的剔除，会对地壳形变分析模型参数估计造成极大的干扰，给实际应用带来不利。因此，粗差探测是高精度地壳形变监测数据后处理的重要内容。
目前在GNSS地壳形变数据处理中，常用的粗差探测策略有：人为判定法、粗差探测法、抗差估计法和拟准检定法。对于人为判定法而言，国内外研究主要依据人为经验和先验信息对粗差进行探测与剔除，但当先验信息不足时只依赖人为经验会出现误判的弊端；对于粗差探测法而言，国内外研究主要是对观测值后验残差进行假设检验，以实现对粗差的探测与辩识，常用的检验方法有正态分布检验、均匀分布检验、F分布检验等，但当观测值中含有多个粗差时，假设检验过程易受到干扰，此外，观测值后验残差易受模型平差结构的影响，难以全面反映出实际的误差水平；对于抗差估计而言，国内外研究主要通过构建等价权函数，如IGG方案、Huber权函数和双因子等价权等对含有粗差数据的观测值进行降权处理，从而达到削弱含粗差观测值对于模型参数估计影响的目的，但该方法易受最小二乘残差初值的影响而陷入局部收敛的误区，从而导致模型参数估计的不准确。对于拟准检定而言，国内外研究主要借鉴拟稳平差的思想，建立真误差数学模型，提出了“拟准观测”的概念，通过拟准检定解算粗差，粗差检测效果良好。该方法的关键在于拟准观测值的选择，当粗差占比过高时，现有的技术无法有效自动化的选择拟准观测值，从而致使模型参数估计陷入崩溃。
发明内容
本发明的目的在于克服现有技术的缺陷，提出了一种基于拟准检定的超强崩溃污染率抗差估计算法，该方法首先利用K均值聚类算法实现了任意粗差占比情况下对于拟准观测值的自动化选择，以实现粗差的粗识别，然后以拟准真误差作为抗差估计中等价权函数的初值进行迭代计算，以实现粗差的精识别和模型参数的超强崩溃污染抗差估计。该方法可以有效避免粗差占比过高时传统粗差探测策略易崩溃失效的问题，为复杂场景中地壳形变监测数据的粗差探测和模型参数估计提供一种有效的处理方法。
本发明的技术方案如下：
一种基于拟准检定的超强崩溃污染率抗差估计算法，包括以下步骤
(1)基于研究区域内所有的GNSS水平速度场，建立整体旋转与均匀应变模型，在最小二乘准则下求解出模型参数与对应的最小二乘改正数V；
(2)求出全部最小二乘改正数V的绝对值的中位值，并将其作为单位权中误差的估计值，据此计算拟准观测值初选指标，并对其按从小到大的顺序进行排序；
(3)根据排序结果进行多维完整搜索，利用K均值聚类对拟准观测值和非拟准观测值进行智能聚类，实现拟准观测值的自动化选择；
(4)根据选定的拟准观测值，执行拟准检定解算，获得拟准真误差ΔN；
(5)选择抗差等价权函数，以真误差拟准检定解ΔN为初值，构建组合算法，迭代计算至模型参数收敛，获得区域地壳运动模型参数的超强崩溃污染率抗差解
Figure GDA0003130083790000021
并探测粗差大小及其位置。
进一步地，所述的整体旋转与均匀应变模型为：
Figure GDA0003130083790000022
上式中，L为观测向量，包含东向和北向速度；Ve、Vn分别表示东向和北向速度；R′表示地球半径，
Figure GDA0003130083790000023
表示测站经纬度位置；
Figure GDA0003130083790000024
表示研究区域的几何中心；εee、εen、εnn分别表示三个主应变参数；ωx、ωy、ωz分别表示三个欧拉参数。
通过最小二乘准则求解模型参数与对应的最小二乘改正数V的函数模型如下：
V＝Aβ-L
Figure GDA0003130083790000025
即整体旋转与均匀应变模型系数矩阵；L为观测向量；
Figure GDA0003130083790000031
进一步地，所述的拟准观测值初选指标的获得方法为：
根据解算的最小二乘改正数V，取其绝对值的中位数为相应的单位权中误差估计值，将最小二乘改正数的绝对值除以单位权中误差和正交补投影矩阵对应元素的乘积，即可获得拟准观测值初选指标。
进一步地，所述的K均值聚类算法，选取欧式距离为相似度指标，明确待分类的样本数据集和要分类形成簇(相似度差异较小的样本)的数量后，通过启发式迭代方法，使得聚类目标实施的各类簇每一个点到相应簇中心的聚类平方和最小，即最小化，直至迭代过程中各簇中心的位置不发生变动。
Figure GDA0003130083790000032
上式中，n为分类样本个数；k为类别数量；xi为样本值；dk代表与xi最近的簇中心点，即对应类别中各数据点的平均值。
进一步地，所述的拟准检定的计算公式如下：
RΔ＝-RL+X
GΔ＝w
上式中，X是拟合残差；w是常数向量；G是系数阵；L为观测向量；R为正交补投影矩阵；Δ为真误差。
通过构造拉格朗日函数求解条件极值，得到的真误差拟准解的模型为：
Figure GDA0003130083790000033
上式中，Δa为拟准观测值对应的真误差；Δb为非拟准观测值对应的真误差；GN＝(0Aa T)；Aa是选定的拟准观测值所对应的模型系数矩阵。
进一步地，所述的抗差等价权函数公式为：
Figure GDA0003130083790000034
上式中，
Figure GDA0003130083790000041
为最小二乘改正数，V中第i个元素的标准化结果；k0和k1为常量；
Figure GDA0003130083790000042
Vvi为第i个观测值的观测误差，可由观测文件读取；k0取值1.5，k1取值3.0。
区域地壳运动模型参数的迭代抗差解算式为：
Figure GDA0003130083790000043
上式中，k为迭代次数；A为整体旋转与均匀应变模型系数矩阵；
Figure GDA0003130083790000044
为抗差等价权函数计算得到的等价权重；L为观测向量。
通过迭代计算获得模型参数的超强崩溃污染率抗差估计值，并识别定位粗差大小及位置，包括：
执行抗差等价权函数，多次迭代计算，直至相邻两次迭代计算的模型参数抗差估计值之差小于给定的阈值，至此获得模型参数的超强崩溃污染率抗差解；识别观测值中等价权为0的位置作为粗差估计位置，并以相对应的最后一次迭代计算得到的最小二乘残差作为粗差大小的估计值。
本发明与现有技术相比具有以下优点：
本发明先结合K均值聚类算法，利用其智能分类的优势实现不受崩溃污染率限制的拟准观测值自动化选择和粗差数据的粗识别，有效避免了粗差占比超过50％时，现有算法易崩溃的缺陷；然后采取拟准检定算法与等价权函数相组合的策略，以拟准真误差和拟准观测值对应拟准真误差的单位权中误差作为等价权函数的初值，从而使迭代模型的初值更符合实际情况，有效避免了因初值不准确而引起的迭代计算陷入局部收敛的弊端。该方法实现过程简单，应用效果明显，可实现区域地壳运动模型参数的超强崩溃污染率抗差估计。
附图说明
图1是本发明流程图；
图2是中国地壳运动监测网络和中国陆态观测网络获取的青藏高原东北缘地壳运动高精度GNSS水平速度场；
图3分别为最小二乘估计、常规抗差估计、基于残差中位数的抗差估计、基于拟准检定的超强崩溃污染率抗差估计算法，对区域整体旋转与均匀应变模型的部分参数(欧拉参数)估计的残差分布图；
具体实施方式
下面结合附图对本发明做进一步的描述。
参照图1，本发明的具体实施步骤如下：
步骤1.基于研究区域内所有的GNSS水平速度场，建立整体旋转与均匀应变模型，根据最小二乘原理，求解模型参数与对应的最小二乘改正数V；
其中，所述的整体旋转与均匀应变的模型为：
Figure GDA0003130083790000051
上式中，L为观测向量，包含东向和北向速度；Ve、Vn分别表示东向和北向速度；R′表示地球半径；
Figure GDA0003130083790000052
表示测站经纬度位置；
Figure GDA0003130083790000053
表示研究区域的几何中心；εee、εen、εnn分别表示三个主应变参数，分别为东西向线应变、东西南北向剪应变和南北向线应变；ωx、ωy、ωz分别表示三个欧拉参数，分别为东西向、南北向、垂向欧拉参数。
通过最小二乘准则计算在最小二乘准则下求解出模型参数与对应的最小二乘改正数V的公式如下：
V＝Aβ-L
Figure GDA0003130083790000054
即整体旋转与均匀应变模型系数矩阵；
Figure GDA0003130083790000055
步骤2.根据求出的最小二乘改正数V，取其绝对值的中位数作为单位权中误差的估计值
Figure GDA0003130083790000056
据此计算拟准观测值初选指标u：
Figure GDA0003130083790000057
上式中，V为最小二乘改正数；
Figure GDA0003130083790000058
为单位权中误差；R为正交投影矩阵。
为便于后续多维完整搜索计算，对上述拟准初选指标按从小到大的顺序进行排序得到u′。
下面的步骤开始进行拟准观测值的选择，本发明提出的利用K均值聚类算法自动化选择拟准观测值也是从该步骤分步引入的。
步骤3.根据排序结果u′进行多维完整搜索，利用K均值聚类对拟准观测值和非拟准观测值进行智能聚类，实现拟准观测值的自动化选择；
其中，K均值聚类算法采用现有的常规算法，选取欧式距离为相似度指标，明确待分类的样本数据集和要分类形成簇(相似度差异较小的样本)的数量后，通过启发式迭代方法确定各个簇中心(聚类中心)，使得聚类目标实施的各类簇每一个点到相应簇中心的聚类平方和J最小，即最小化，直至迭代过程中各簇中心的位置不发生变动。
Figure GDA0003130083790000061
上式中，n为分类样本个数；k为类别数量；xi为样本值；dk代表与xi最近的簇中心点，即对应类别中各数据点的平均值。
而在经典拟准检定算法中选择拟准观测值时，首先对拟准观测值初选指标进行人为选择或者设置阈值，阈值取值范围通常为1.0-3.0，初步选择拟准观测值，由初选的拟准观测值计算对应的真误差估计值；然后计算拟准观测值复选指标W
Figure GDA0003130083790000062
上式中，Δ为所有观测值对应的拟准真误差估计值；Δ0为初选拟准观测值对应的真误差估计值；Δ0为初选拟准观测值对应的真误差估计值；r为初选拟准观测值个数。
在此基础上，执行复选过程，同初选过程类似，对拟准观测值复选指标进行人为选择或者设置阈值，取值范围通常为1.0-3.0，一般而言复选过程的阈值要小于初选过程的阈值，然后根据复选拟准观测值计算对应的真误差估计值，以此作为真误差拟准检定解，非拟准观测值则被判定为粗差数据。
然而，上述拟准检定算法选择拟准观测值时，通过对复选指标和初选指标人为判定或者设置阈值，会导致选择人为主观性强、无法自动化实现和粗差占比过高时，阈值失效导致算法崩溃，粗差数据的识别定位也随之失效。本发明为了在提高算法的崩溃污染率的同时来实现拟准观测值的自动化选择。首先对拟准初选指标从从小到大进行排序，取排序结果的前m+1个最小值作为拟准观测值，执行拟准鉴定解算，获取拟准观测值对应的拟准真误差方差
Figure GDA0003130083790000063
1表示添加一个多余观测后参与计算，然后逐个增加多余观测值；经过n-m次循环计算后，所有观测数据纳为拟准观测值，每纳入一个多余观测对应计算出的方差均被记录
Figure GDA0003130083790000064
这个过程称为n-m维完整搜索。对上述记录的方差序列结果求差，采用K均值聚类方法对方差之差分两类聚类处理，选择方差较小的一类作为拟准观测值，其中m为模型参数个数，n为观测值个数，这一选择策略不会因粗差占比过高而使算法崩溃。同时不受人为干扰，得以自动化实现。
下面通过步骤3选择出的拟准观测值，构建拟准检定模型以进行后续计算。
步骤4.在最小二乘准则下推导出观测值L与其真误差Δ的对应关系：
RΔ＝RL
上式中，R为正交投影矩阵；Δ为真误差；L为观测向量。
步骤3选定的拟准观测值对应的真误差为Δa，非拟准观测值对应的真误差为Δb，以拟准观测值对应的拟准真误差平方和极小为条件解算真误差以观测值的关系式。
拟准真误差平方和极小条件为：
Figure GDA0003130083790000071
在此基础上，借用拟稳平差的思路作为附加条件，构建真误差与观测值的解算模型：
Figure GDA0003130083790000072
上式中，X＝R(Δ+L)是拟合残差；w是常数向量；G是系数阵。
选定拟准观测值后，真误差与观测值的解算模型可进一步表示为：
Figure GDA0003130083790000073
上式中，Aa是选定的拟准观测值所对应的模型系数矩阵；GN＝(0Aa T)；Δa为拟准观测值对应的真误差；Δb为非拟准观测值对应的真误差；通过构造拉格朗日函数求解条件极值，便可求出真误差的拟准解：
Figure GDA0003130083790000074
至此，已通过拟准检定算法解算出真误差拟准解，实现粗差数据的初步识别，拟准检定计算在这一步已经结束。
但本发明仅通过具有自动化选择拟准观测值的拟准检定算法只能实现在任意粗差占比情况下的粗差数据的粗识别，并不能实现粗差数据的精确识别和模型参数的抗差估计，其中部分中小粗差数据虽然不至于被剔除模型，但其对模型参数估计的干扰同样不可忽略。因此，本发明在拟准检定算法解算真误差拟准解的基础上再进一步引入抗差等价权函数，构建组合算法，通过对粗差数据赋0权和中小粗差赋弱权，实现模型参数的超强崩溃污染率抗差估计和粗差数据的精细识别；下面进行组合算法处理。
步骤5.在拟准检定解算基础上，引入IGG3抗差等价权函数，对中小粗差数据进行降权处理，以获得更为可靠的区域地壳运动整体旋转与均匀应变模型参数估计值
Figure GDA0003130083790000075
具体步骤如下：
建立IGG3抗差等价权函数模型：
Figure GDA0003130083790000081
上式中，
Figure GDA0003130083790000082
为最小二乘改正数，V中第i个元素的标准化结果；k0和k1为常量；
Figure GDA0003130083790000083
Vvi为第i个观测值的观测误差，可由观测文件读取；k0＝1.0～1.5；k1＝2.5～3.0。
在拟准检定解算基础上，以拟准真误差代替最小二乘残差，拟准观测值对应的拟准真误差单位权中误差代替最小二乘残差单位权中误差作为上述抗差等价权函数的初值，继而建立组合算法。
通过上述解算的抗差等价权，执行区域地壳运动模型参数的迭代抗差解算式：
Figure GDA0003130083790000084
上式中，k为迭代次数；A为整体旋转与均匀应变模型系数矩阵；
Figure GDA0003130083790000085
为抗差等价权函数计算得到的等价权；L为观测向量。
通过迭代计算，直至相邻两次迭代计算的模型参数抗差估计值之差小于给定阈值而停止计算，最终获得模型参数的超强崩溃污染率抗差解。并通过判断等价权为0的观测值和其对应的抗差最小二乘残差解来识别定位粗差的位置及大小，具体为将权重为0观测值的位置(经纬度)作为粗差出现的位置，对应的抗差最小二乘残差值作为粗差大小的估计值。
利用解算出的区域地壳运动模型参数的超强崩溃污染率抗差估计值，并结合研究域已有地质勘查资料，可进行区域地壳构造运动特性的研究，并进一步剔除粗差数据，获得干净可靠的地壳形变监测数据为相关地壳形变特征研究提供基础资料。
本发明的效果可以通过下述算例加以说明：
1.算例设置
这里为证明本发明提出新算法的准确性及有效性，设置模拟算例。
选取我国典型地壳构造强活动区域青藏高原东北缘区域，利用该区域地壳运动高精度GNSS水平速度场模拟值来检验本发明提出的基于拟准检定的超强崩溃污染率抗差估计新算法。GNSS观测数据来源于中国地壳运动观测网络和中国陆态观测网络，时间跨度为1999-2018年共217个站点数据，建立区域整体旋转与均匀应变模型反演其模型参数，并以此为参数真值在0.5°×0.5°格网间隔正演该区域GNSS速度场模拟值，GNSS站点速度场如图2所示。通过随机添加粗差数据设置模拟试验。
2.算例结果
为了对比不同方法对参数估计和粗差探测的影响，分别利用四种算法对含有粗差数据的青藏高原东北缘地壳运动GNSS水平速度场模拟值进行相应的粗差探测和参数抗差估计，得到四种算法对应的参数估计残差统计结果和粗差探测成功率。限于篇幅，只给出粗差占比为60％时对应的结果。最小二乘估计、常规抗差估计、基于残差中位数的抗差估计和基于拟准检定的超强崩溃污染率抗差估计分别用模型1、2、3和4表示。
图2给出了实测GNSS速度场的位移空间分布图，从图中可以看出研究区域处于青藏高原东北缘以及各观测点位的噪声水平。
图3给出了四种处理模型对应的部分参数(ωx、ωy和ωz)估计残差，可以看出当粗差占比为60％时，模型1、2和3所估计的参数残差异常增大，模型处于崩溃状态，而模型4所估计的参数残差则明显处于较低水平，说明模型4能够更好地抵御粗差影响，可对参数进行超强崩溃污染率抗差估计。
表1给出了模型参数反演结果。建立区域整体旋转与均匀应变模型，在实测数据的约束下反演出模型参数，并估计其中误差信息，模型参数反演值视作真值用于模拟算例参数估计的精度评定。
表1模型参数反演结果(1×10-9)
Figure GDA0003130083790000091
表2统计了当粗差占比为60％时，四种模型对于参数估计残差的相关统计指标(最大残差、最小残差、残差平均值和均方根误差)。从表2可以看出粗差占比为60％时，模型1的参数估计效果最差，其次是模型2和模型3，模型4的参数估计效果最好。
表2模型参数估计残差统计结果(1×10-9)
Figure GDA0003130083790000092
Figure GDA0003130083790000101
表3统计了模型2、模型3和模型4的粗差探测成功率比较结果，从表中可以看出当粗差占比高于60％时，只有本发明提出的基于拟准检定的超强崩溃污染率抗差估计法(模型4)能够准确探测粗差，其他两种方法则完全失效，这表明本发明提出的方法更具有普适性。
表3粗差探测成功率比较
Figure GDA0003130083790000102

Claims (9)

1.一种基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于，包括以下步骤：
(1)基于研究区域内所有的GNSS水平速度场，建立整体旋转与均匀应变模型，在最小二乘准则下求解出模型参数与对应的最小二乘改正数V；
(2)求出全部最小二乘改正数V的绝对值的中位值，并将其作为单位权中误差的估计值，据此计算拟准观测值初选指标，并对其按从小到大的顺序进行排序；
具体实施步骤为：根据解算的最小二乘改正数V，取其绝对值的中位数作为单位权中误差的估计值
Figure FDA0003185694330000011
将最小二乘改正数的绝对值除以单位权中误差和正交补投影矩阵对应元素的乘积，据此计算拟准观测值初选指标u：
Figure FDA0003185694330000012
上式中，V为最小二乘改正数；
Figure FDA0003185694330000013
为单位权中误差；R为正交投影矩阵；
(3)根据排序结果进行多维完整搜索，利用K均值聚类对拟准观测值和非拟准观测值进行智能聚类，实现拟准观测值的自动化选择；
具体实施步骤为：首先对拟准初选指标从小到大进行排序，取排序结果的前m+1个最小值作为拟准观测值，执行拟准鉴定解算，获取拟准观测值对应的拟准真误差方差
Figure FDA0003185694330000014
然后逐个增加多余观测值；经过n-m次循环计算后，所有观测数据纳为拟准观测值，每纳入一个多余观测对应计算出的方差均被记录为
Figure FDA0003185694330000015
这个过程称为n-m维完整搜索，对上述记录的方差序列结果求差，采用K均值聚类方法对方差之差分两类聚类处理，选择方差较小的一类作为拟准观测值，其中m为模型参数个数，n为观测值个数；
(4)根据选定的拟准观测值，执行拟准检定解算，获得拟准真误差ΔN；
(5)选择抗差等价权函数，以真误差拟准检定解ΔN为初值，构建组合算法，迭代计算至模型参数收敛，获得区域地壳运动模型参数的超强崩溃污染率抗差解
Figure FDA0003185694330000016
并探测粗差大小及其位置；
具体实施步骤为：在拟准检定解算基础上，以拟准真误差代替最小二乘残差，拟准观测值对应的拟准真误差单位权中误差代替最小二乘残差单位权中误差作为上述抗差等价权函数的初值，继而建立组合算法；
通过上述解算的抗差等价权，执行区域地壳运动模型参数的迭代抗差解算式：
Figure FDA0003185694330000017
上式中，k为迭代次数；A为整体旋转与均匀应变模型系数矩阵；
Figure FDA0003185694330000021
为抗差等价权函数计算得到的等价权；L为观测向量。
2.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于，步骤(1)中所述的整体旋转与均匀应变模型为：
Figure FDA0003185694330000022
上式中，L为观测向量，包含东向和北向速度；Ve、Vn分别表示东向和北向速度；R′表示地球半径，
Figure FDA0003185694330000023
表示测站经纬度位置；
Figure FDA0003185694330000024
表示研究区域的几何中心；εee、εen、εnn分别表示三个主应变参数；ωx、ωy、ωz分别表示三个欧拉参数。
3.根据权利要求2所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于，步骤(1)中通过最小二乘准则求解模型参数与对应的最小二乘改正数V的函数模型如下：
V＝Aβ-L
Figure FDA0003185694330000025
上式中，A为整体旋转与均匀应变模型系数矩阵；
Figure FDA0003185694330000026
L为观测向量，包含东向和北向速度；R′表示地球半径，
Figure FDA0003185694330000027
表示测站经纬度位置；
Figure FDA0003185694330000028
表示研究区域的几何中心；εee、εen、εnn分别表示三个主应变参数；ωx、ωy、ωz分别表示三个欧拉参数。
4.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(2)中拟准观测值初选指标的获得方法为：根据解算的最小二乘改正数V，取其绝对值的中位数为相应的单位权中误差估计值，将最小二乘改正数的绝对值除以单位权中误差和正交补投影矩阵对应元素的乘积，即可获得拟准观测值初选指标。
5.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(3)中所述的K均值聚类算法是选取欧式距离为相似度指标，明确待分类的样本数据集和待分类形成簇的数量后，通过启发式迭代方法确定各个簇中心，使得聚类目标实施的各类簇每一个点到相应簇中心的聚类平方和J最小，直至迭代过程中各簇中心的位置不发生变动；
Figure FDA0003185694330000031
上式中，n为分类样本个数；k为类别数量；xi为样本值；dk代表与xi最近的簇中心点，即对应类别中各数据点的平均值。
6.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(4)中所述的拟准检定的计算公式如下：
RΔ＝-RL+X
GΔ＝w
上式中，X是拟合残差；w是常数向量；G是系数阵；L为观测向量；R为正交补投影矩阵；Δ为真误差。
7.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(5)中所述的抗差等价权函数为：
Figure FDA0003185694330000032
上式中，
Figure FDA0003185694330000033
为最小二乘改正数V中第i个元素的标准化结果；k0和k1为常量；
Figure FDA0003185694330000034
Vvi为第i个观测值的观测误差，可由观测文件读取；k0取值1.5，k1取值3.0。
8.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(5)中区域地壳运动模型参数的超强崩溃污染率抗差解算式为：
Figure FDA0003185694330000035
上式中，k为迭代次数；A为整体旋转与均匀应变模型系数矩阵；
Figure FDA0003185694330000036
为抗差等价权函数计算得到的等价权重，L为观测向量。
9.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于，步骤(5)中粗差大小及其位置的探测方法为：
执行抗差等价权函数，多次迭代计算，直至相邻两次迭代计算的模型参数抗差估计值之差小于给定的阈值，至此获得模型参数的超强崩溃污染率抗差解；
识别观测值中等价权为0的位置作为粗差估计位置，并以相对应的最后一次迭代计算得到的最小二乘残差作为粗差大小的估计值。
CN202011047884.8A 2020-09-29 2020-09-29 一种基于拟准检定的超强崩溃污染率抗差估计算法 Active CN112131752B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011047884.8A CN112131752B (zh) 2020-09-29 2020-09-29 一种基于拟准检定的超强崩溃污染率抗差估计算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011047884.8A CN112131752B (zh) 2020-09-29 2020-09-29 一种基于拟准检定的超强崩溃污染率抗差估计算法

Publications (2)

Publication Number Publication Date
CN112131752A CN112131752A (zh) 2020-12-25
CN112131752B true CN112131752B (zh) 2021-09-10

Family

ID=73844570

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011047884.8A Active CN112131752B (zh) 2020-09-29 2020-09-29 一种基于拟准检定的超强崩溃污染率抗差估计算法

Country Status (1)

Country Link
CN (1) CN112131752B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114117802B (zh) * 2021-11-29 2023-09-26 同济大学 一种基于极大后验估计的多粗差探测方法、装置及介质
CN114332389B (zh) * 2021-12-24 2022-11-08 中国测绘科学研究院 一种三维地壳形变模型的构建方法和装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101887128A (zh) * 2010-07-09 2010-11-17 中国科学院测量与地球物理研究所 确定全球卫星导航系统导航卫星频间偏差的方法
CN111077550A (zh) * 2019-12-26 2020-04-28 广东星舆科技有限公司 一种应用于智能终端rtd定位的粗差探测方法及系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101887128A (zh) * 2010-07-09 2010-11-17 中国科学院测量与地球物理研究所 确定全球卫星导航系统导航卫星频间偏差的方法
CN111077550A (zh) * 2019-12-26 2020-04-28 广东星舆科技有限公司 一种应用于智能终端rtd定位的粗差探测方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
利用L1范数和中位数选取拟准观测值;赵俊等;《武 汉 大 学 学 报 · 信 息 科 学 版》;20180831;全文 *
最小二乘配置与REHSM求解GPS应变场的方法;管守奎等;《大地测量与地球动力学》;20150831;全文 *

Also Published As

Publication number Publication date
CN112131752A (zh) 2020-12-25

Similar Documents

Publication Publication Date Title
CN109682382B (zh) 基于自适应蒙特卡洛和特征匹配的全局融合定位方法
CN112131752B (zh) 一种基于拟准检定的超强崩溃污染率抗差估计算法
CN109543356B (zh) 考虑空间非平稳性的海洋内部温盐结构遥感反演方法
CN110398753A (zh) Gnss测站坐标时间序列周期性探测方法及系统
CN106912105A (zh) 基于pso_bp神经网络的三维定位方法
CN110536257B (zh) 一种基于深度自适应网络的室内定位方法
CN112417573A (zh) 基于ga-lssvm与nsga-ⅱ盾构下穿既有隧道施工多目标优化的方法
CN111190211B (zh) 一种gps失效位置预测定位方法
CN108716904B (zh) 基于有限测斜仪测点测值的坝体挠度获取方法
CN110205909B (zh) 一种基于沥青层当量温度的路面结构弯沉系指标的温度修正方法
CN106096847A (zh) 一种模糊变权工程地质环境质量评价方法
CN109188481A (zh) 一种拟合gps高程异常的新方法
Qu et al. A robust estimation algorithm for the increasing breakdown point based on quasi-accurate detection and its application to parameter estimation of the GNSS crustal deformation model
CN102506753B (zh) 基于十四点球面小波变换的不规则零件形状差异检测方法
CN113486591B (zh) 一种卷积神经网络结果的重力多参量数据密度加权反演方法
CN106597415A (zh) 一种高斯噪声下提高稀疏孔径成像系统误差检测精度的方法
CN113449254A (zh) 任意网型变形监测稳定性分析方法及监控点位置确定方法
CN113139337A (zh) 一种用于湖泊地形模拟的分区插值处理方法与装置
CN106501815B (zh) 一种仅天基测角跟踪的空间目标轨道机动融合检测方法
CN105204047A (zh) 一种卫星导航系统中观测量单个粗差的探测与修复方法
CN108572361A (zh) 机载激光雷达系统设备集成安置角检校方法及装置
CN114189313B (zh) 一种电表数据重构方法及装置
CN112986948B (zh) 基于InSAR技术的建筑形变监测方法和装置
CN114488247A (zh) 一种基于高精度北斗差分定位分析装备机动能力的方法
CN106125149B (zh) 点质量模型中浅层高分辨率点质量最佳埋藏深度确定方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant