CN112131752A

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

Info

Publication number: CN112131752A
Application number: CN202011047884.8A
Authority: CN
Inventors: 瞿伟; 陈海禄; 张勤; 高源�; 梁世川; 韩亚茜
Original assignee: Changan University
Current assignee: Changan University
Priority date: 2020-09-29
Filing date: 2020-09-29
Publication date: 2020-12-25
Anticipated expiration: 2040-09-29
Also published as: CN112131752B

Abstract

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

Description

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

技术领域

本发明属于高精度地壳形变监测数据处理领域，涉及到一种基于拟准检定的地壳运动GNSS水平速度场粗差探测与模型参数估计技术，该算法以高精度地壳形变监测为实际应用背景，可用于观测环境复杂、局部地壳活动活跃等场景中的地壳形变高精度监测应用方向。

背景技术

随着现代空间大地测量技术的迅猛发展，特别是以全球导航卫星系统(GNSS)为代表的空间监测技术的现代化建设，可利用其以厘米级甚至毫米级的精度实现对地壳水平运动的高精度监测。但受监测环境干扰、局部地壳活动或地震等因素影响，导致所观测到的GNSS速度场数据中通常含有粗差，若粗差无法被有效的剔除，会对地壳形变分析模型参数估计造成极大的干扰，给实际应用带来不利。因此，粗差探测是高精度地壳形变监测数据后处理的重要内容。

目前在GNSS地壳形变数据处理中，常用的粗差探测策略有：人为判定法、粗差探测法、抗差估计法和拟准检定法。对于人为判定法而言，国内外研究主要依据人为经验和先验信息对粗差进行探测与剔除，但当先验信息不足时只依赖人为经验会出现误判的弊端；对于粗差探测法而言，国内外研究主要是对观测值后验残差进行假设检验，以实现对粗差的探测与辩识，常用的检验方法有正态分布检验、均匀分布检验、F分布检验等，但当观测值中含有多个粗差时，假设检验过程易受到干扰，此外，观测值后验残差易受模型平差结构的影响，难以全面反映出实际的误差水平；对于抗差估计而言，国内外研究主要通过构建等价权函数，如IGG方案、Huber权函数和双因子等价权等对含有粗差数据的观测值进行降权处理，从而达到削弱含粗差观测值对于模型参数估计影响的目的，但该方法易受最小二乘残差初值的影响而陷入局部收敛的误区，从而导致模型参数估计的不准确。对于拟准检定而言，国内外研究主要借鉴拟稳平差的思想，建立真误差数学模型，提出了“拟准观测”的概念，通过拟准检定解算粗差，粗差检测效果良好。该方法的关键在于拟准观测值的选择，当粗差占比过高时，现有的技术无法有效自动化的选择拟准观测值，从而致使模型参数估计陷入崩溃。

发明内容

本发明的目的在于克服现有技术的缺陷，提出了一种基于拟准检定的超强崩溃污染率抗差估计算法，该方法首先利用K均值聚类算法实现了任意粗差占比情况下对于拟准观测值的自动化选择，以实现粗差的粗识别，然后以拟准真误差作为抗差估计中等价权函数的初值进行迭代计算，以实现粗差的精识别和模型参数的超强崩溃污染抗差估计。该方法可以有效避免粗差占比过高时传统粗差探测策略易崩溃失效的问题，为复杂场景中地壳形变监测数据的粗差探测和模型参数估计提供一种有效的处理方法。

本发明的技术方案如下：

一种基于拟准检定的超强崩溃污染率抗差估计算法，包括以下步骤

(1)基于研究区域内所有的GNSS水平速度场，建立整体旋转与均匀应变模型，在最小二乘准则下求解出模型参数与对应的最小二乘改正数V；

(2)求出全部最小二乘改正数V的绝对值的中位值，并将其作为单位权中误差的估计值，据此计算拟准观测值初选指标，并对其按从小到大的顺序进行排序；

(3)根据排序结果进行多维完整搜索，利用K均值聚类对拟准观测值和非拟准观测值进行智能聚类，实现拟准观测值的自动化选择；

(4)根据选定的拟准观测值，执行拟准检定解算，获得拟准真误差Δ_N；

(5)选择抗差等价权函数，以真误差拟准检定解Δ_N为初值，构建组合算法，迭代计算至模型参数收敛，获得区域地壳运动模型参数的超强崩溃污染率抗差解

并探测粗差大小及其位置。

进一步地，所述的整体旋转与均匀应变模型为：

上式中，L为观测速度向量，包含东向和北向速度；V_e、V_n分别表示东向和北向速度；R_earth表示地球半径，

表示测站经纬度位置；

表示研究区域的几何中心；ε_ee、ε_en、ε_nn分别表示三个主应变参数；ω_x、ω_y、ω_z分别表示三个欧拉参数。

通过最小二乘准则求解模型参数与对应的最小二乘改正数V的函数模型如下：

V＝Aβ-L

上式中

即整体旋转与均匀应变模型系数矩阵；L为观测向量；

进一步地，所述的拟准观测值初选指标的获得方法为：

根据解算的最小二乘改正数V，取其绝对值取中位数为相应的单位权中误差估计值，将最小二乘改正数的绝对值除以单位权中误差和正交补投影矩阵对应元素的乘积，即可获得拟准观测值初选指标。

进一步地，所述的K均值聚类算法，选取欧式距离为相似度指标，明确待分类的样本数据集和要分类形成簇(相似度差异较小的样本)的数量后，通过启发式迭代方法，使得聚类目标实施的各类簇每一个点到相应簇中心的聚类平方和最小，即最小化，直至迭代过程中各簇中心的位置不发生变动。

上式中，n为分类样本个数；k为类别数量；x_i为样本值；d_k代表与x_i最近的簇中心点，即对应类别中各数据点的平均值。

进一步地，所述的拟准检定的计算公式如下：

RΔ＝-RL+X

GΔ＝w

上式中，X是拟合残差；w是常数向量；G是系数阵；L为观测向量；R为正交补投影矩阵；Δ为真误差。

通过构造拉格朗日函数求解条件极值，得到的真误差拟准解的模型为：

上式中，Δ_a为拟准观测值对应的真误差；Δ_b为非拟准观测值对应的真误差；

G_N(0A_a ^T)；A_a是选定的拟准观测值所对应的模型系数矩阵。

进一步地，所述的抗差等价权函数公式为：

上式中，

为最小二乘改正数向量V中第i个元素的标准化结果；k₀和k₁为常量；通常k₀＝1.0～1.5；k₁＝2.5～3.0；

Vv_i为第i个观测值的观测误差，可由观测文件读取；本发明中k₀取值1.5，k₁取值3.0。

区域地壳运动模型参数的迭代抗差解算式为：

上式中，k为迭代次数；A为整体旋转与均匀应变模型系数矩阵；

为抗差等价权函数计算得到的等价权重。

通过迭代计算获得模型参数的超强崩溃污染率抗差估计值，并识别定位粗差大小及位置，包括：

执行抗差等价权函数，多次迭代计算，直至相邻两次迭代计算的模型参数抗差估计值之差小于给定的阈值，至此获得模型参数的超强崩溃污染率抗差解；识别观测值中等价权为0的位置作为粗差估计位置，并以相对应的最后一次迭代计算得到的最小二乘残差作为粗差大小的估计值。

本发明与现有技术相比具有以下优点：

本发明先结合K均值聚类算法，利用其智能分类的优势实现不受崩溃污染率限制的拟准观测值自动化选择和粗差数据的粗识别，有效避免了粗差占比超过50％时，现有算法易崩溃的缺陷；然后采取拟准检定算法与等价权函数相组合的策略，以拟准真误差和拟准观测值对应拟准真误差的单位权中误差作为等价权函数的初值，从而使迭代模型的初值更符合实际情况，有效避免了因初值不准确而引起的迭代计算陷入局部收敛的弊端。该方法实现过程简单，应用效果明显，可实现区域地壳运动模型参数的超强崩溃污染率抗差估计。

附图说明

图1是本发明流程图；

图2是中国地壳运动监测网络和中国陆态观测网络获取的青藏高原东北缘地壳运动高精度GNSS水平速度场；

图3分别为最小二乘估计、常规抗差估计、基于残差中位数的抗差估计、基于拟准检定的超强崩溃污染率抗差估计算法，对区域整体旋转与均匀应变模型的部分参数(欧拉参数)估计的残差分布图；

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1.基于研究区域内所有的GNSS水平速度场，建立整体旋转与均匀应变模型，根据最小二乘原理，求解模型参数与对应的最小二乘改正数向量V；

其中，所述的整体旋转与均匀应变的模型为：

上式中，L为观测速度向量，包含东向和北向速度；V_e、V_n分别表示东向和北向速度；R_earth表示地球半径；

表示测站经纬度位置；

表示研究区域的几何中心；ε_ee、ε_en、ε_nn分别表示三个主应变参数，分别为东西向线应变、东西南北向剪应变和南北向线应变；ω_x、ω_y、ω_z分别表示三个欧拉参数，分别为东西向、南北向、垂向欧拉参数。

通过最小二乘准则计算在最小二乘准则下求解出模型参数与对应的最小二乘改正数V的公式如下：

V＝Aβ-L

上式中

即整

体旋转与均匀应变模型系数矩阵；

步骤2.根据求出的最小二乘改正数V，取其绝对值的中位数作为单位权中误差的估计值

据此计算拟准观测值初选指标u：

上式中，V为最小二乘改正数向量；

为单位权中误差；R为正交投影矩阵。

为便于后续多维完整搜索计算，对上述拟准初选指标按从小到大的顺序进行排序得到u′。

下面的步骤开始进行拟准观测值的选择，本发明提出的利用K均值聚类算法自动化选择拟准观测值也是从该步骤分步引入的。

步骤3.根据排序结果u′进行多维完整搜索，利用K均值聚类对拟准观测值和非拟准观测值进行智能聚类，实现拟准观测值的自动化选择；

其中，K均值聚类算法采用现有的常规算法，选取欧式距离为相似度指标，明确待分类的样本数据集和要分类形成簇(相似度差异较小的样本)的数量后，通过启发式迭代方法确定各个簇中心(聚类中心)，使得聚类目标实施的各类簇每一个点到相应簇中心的聚类平方和J最小，即最小化，直至迭代过程中各簇中心的位置不发生变动。

而在经典拟准检定算法中选择拟准观测值时，首先对拟准观测值初选指标进行人为选择或者设置阈值，阈值取值范围通常为1.0-3.0，初步选择拟准观测值，由初选的拟准观测值计算对应的真误差估计值；然后计算拟准观测值复选指标W

上式中，Δ为所有观测值对应的拟准真误差估计值；Δ₀为初选拟准观测值对应的真误差估计值；Δ₀为初选拟准观测值对应的真误差估计值；r为初选拟准观测值个数。

在此基础上，执行复选过程，同初选过程类似，对拟准观测值复选指标进行人为选择或者设置阈值，取值范围通常为1.0-3.0，一般而言复选过程的阈值要小于初选过程的阈值，然后根据复选拟准观测值计算对应的真误差估计值，以此作为真误差拟准检定解，非拟准观测值则被判定为粗差数据。

然而，上述拟准检定算法选择拟准观测值时，通过对复选指标和初选指标人为判定或者设置阈值，会导致选择人为主观性强、无法自动化实现和粗差占比过高时，阈值失效导致算法崩溃，粗差数据的识别定位也随之失效。本发明为了在提高算法的崩溃污染率的同时来实现拟准观测值的自动化选择。首先对拟准初选指标从从小到大进行排序，取排序结果的前m+1个(m为模型参数个数)最小值作为拟准观测值，执行拟准鉴定解算，获取拟准观测值对应的拟准真误差方差

(1表示添加一个多余观测后参与计算)，然后逐个增加多余观测值。经过n-m次(n为观测值个数)循环计算后，所有观测数据纳为拟准观测值，每纳入一个多余观测对应计算出的方差均被记录

这个过程称为n-m维完整搜索。对上述记录的方差序列结果求差，采用K均值聚类方法对方差之差分两类聚类处理，选择方差较小的一类作为拟准观测值，这一选择策略不会因粗差占比过高而使算法崩溃。同时不受人为干扰，得以自动化实现。

下面通过步骤3选择出的拟准观测值，构建拟准检定模型以进行后续计算。

步骤4.在最小二乘准则下推导出观测值L与其真误差Δ的对应关系：

RΔ＝RL

上式中，R为正交投影矩阵；Δ为真误差；L为观测值。

步骤3选定的拟准观测值对应的真误差为Δ_a，非拟准观测值对应的真误差为Δ_b，以拟准观测值对应的拟准真误差平方和极小为条件解算真误差以观测值的关系式。

拟准真误差平方和极小条件为：

在此基础上，借用拟稳平差的思路作为附加条件，构建真误差与观测值的解算模型：

上式中，X＝R(Δ+L)是拟合残差；w是常数向量；G是系数阵。

选定拟准观测值后，真误差与观测值的解算模型可进一步表示为：

上式中，A_a是选定的拟准观测值所对应的模型系数矩阵；G_N＝(0A_a ^T)；Δ_a为拟准观测值对应的真误差；Δ_b为非拟准观测值对应的真误差；通过构造拉格朗日函数求解条件极值，便可求出真误差的拟准解：

至此，已通过拟准检定算法解算出真误差拟准解，实现粗差数据的初步识别，拟准检定计算在这一步已经结束。

但本发明仅通过具有自动化选择拟准观测值的拟准检定算法只能实现在任意粗差占比情况下的粗差数据的粗识别，并不能实现粗差数据的精确识别和模型参数的抗差估计，其中部分中小粗差数据虽然不至于被剔除模型，但其对模型参数估计的干扰同样不可忽略。因此，本发明在拟准检定算法解算真误差拟准解的基础上再进一步引入抗差等价权函数，构建组合算法，通过对粗差数据赋0权和中小粗差赋弱权，实现模型参数的超强崩溃污染率抗差估计和粗差数据的精细识别；下面进行组合算法处理。

步骤5.在拟准检定解算基础上，引入IGG3抗差等价权函数，对中小粗差数据进行降权处理，以获得更为可靠的区域地壳运动整体旋转与均匀应变模型参数估计值

具体步骤如下：

建立IGG3抗差等价权函数模型：

上式中，

为最小二乘改正数向量V中第i个元素的标准化结果；k₀和k₁为常量；

Vv_i为第i个观测值的观测误差，可由观测文件读取；通常k₀＝1.0～1.5；k₁＝2.5～3.0。

在拟准检定解算基础上，以拟准真误差代替最小二乘残差，拟准观测值对应的拟准真误差单位权中误差代替最小二乘残差单位权中误差作为上述抗差等价权函数的初值，继而建立组合算法。

通过上述解算的抗差等价权，执行区域地壳运动模型参数的迭代抗差解算式：

为抗差等价权函数计算得到的等价权。

通过迭代计算，直至相邻两次迭代计算的模型参数抗差估计值之差小于给定阈值而停止计算，最终获得模型参数的超强崩溃污染率抗差解。并通过判断等价权为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)

表2统计了当粗差占比为60％时，四种模型对于参数估计残差的相关统计指标(最大残差、最小残差、残差平均值和均方根误差)。从表2可以看出粗差占比为60％时，模型1的参数估计效果最差，其次是模型2和模型3，模型4的参数估计效果最好。

表2模型参数估计残差统计结果(1×10^-9)

表3统计了模型2、模型3和模型4的粗差探测成功率比较结果，从表中可以看出当粗差占比高于60％时，只有本发明提出的基于拟准检定的超强崩溃污染率抗差估计法(模型4)能够准确探测粗差，其他两种方法则完全失效，这表明本发明提出的方法更具有普适性。

表3粗差探测成功率比较

Claims

1.一种基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于，包括以下步骤：

并探测粗差大小及其位置。

2.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(1)中所述的整体旋转与均匀应变模型为：

表示测站经纬度位置；

3.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(1)中通过最小二乘准则求解模型参数与对应的最小二乘改正数V的函数模型如下：

V＝Aβ-L

上式中

A为整体旋转与均匀应变模型系数矩阵；L为观测向量；

4.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(2)中拟准观测值初选指标的获得方法为：根据解算的最小二乘改正数V，取其绝对值取中位数为相应的单位权中误差估计值，将最小二乘改正数的绝对值除以单位权中误差和正交补投影矩阵对应元素的乘积，即可获得拟准观测值初选指标。

5.根据权利要求1所述的基于拟准检定的超强崩溃污染率抗差估计算法，其特征在于：步骤(3)中所述的K均值聚类算法是选取欧式距离为相似度指标，明确待分类的样本数据集和待分类形成簇的数量后，通过启发式迭代方法确定各个簇中心，使得聚类目标实施的各类簇每一个点到相应簇中心的聚类平方和J最小，直至迭代过程中各簇中心的位置不发生变动；