一种地壳运动GPS水平速度场自适应最小二乘拟合推估算法的制作方法

本发明涉及一种针对地壳运动gps水平速度场有效拟合推估的新算法，具体涉及一种顾及距离尺度因子与自适应调整融合的最小二乘拟合推估算法。

背景技术：

地壳运动是由地壳内部应力变化引起的地壳结构改变、地壳内部物质变位，从而导致地表产生位移的现象，根据其运动方向可分为水平运动和垂直运动。对地壳运动进行高精度的监测可获取到地壳运动速度场，该速度场可为深入研究地壳深部构造运移与变形特性提供有价值的参考。随着现代空间大地测量技术的迅猛发展，特别是以全球导航卫星系统(gps，globalpositioningsystem；beidou，navigationsatellitesystem等)为代表的空间监测技术的快速现代化建设，使得利用这些现代空间大地测量技术在上百公里甚至上千公里距离范围内，可以厘米级甚至毫米级的精度实现对地壳水平运动的高精度监测。然而，由于受自然环境、构造位置、经济费用等因素的影响，使得在一些地壳运动活跃区域难以开展大面积且覆盖均匀的gps监测站点建设，因而难以获取到区域多期高精度重复gps观测站点资料，导致无法准确掌握研究区域地壳实际运动与变形状态。因此，需基于现有的gps监测数据，结合数理知识构建合理的拟合推估模型模型(函数模型与随机模型)，获取区域连续地壳运动形变场为进一步研究区域地壳运动与变形特征服务。

常见的函数模型有二次曲面模型、移动曲面模型、多面函数模型等；常见的随机模型有：权中数法、kringring法等。

上述函数模型虽能较准确反映出研究域趋势性变化，但均未顾及逼近场的随机性；而随机模型多采用统计法建模，只顾及到逼近场的随机性，却忽略了其整体趋势性。

技术实现要素：

本发明的目的是提供一种地壳运动gps水平速度场自适应最小二乘拟合推估算法，用以更加精确地拟合推估反映出区域地壳构造运动特性，为进一步研究区域构造动力学特征提供有价值的基础参考。

为了实现上述任务，本发明采用以下技术方案：

一种地壳运动gps水平速度场自适应最小二乘拟合推估算法，包括以下步骤：

步骤1，基于观测得到的研究域内所有gps监测站点水平运动速度场，求解出研究域的先验整体运动趋势；计算研究域整体运动速度场，进而得到研究域地壳运动观测信号，即验先信号值；

步骤2，根据所述验先信号值，构建拟合经验高斯函数所需的待定参数；

步骤3，引入顾及研究域范围的距离尺度因子，通过选取不同的距离尺度因子，确定经验高斯函数中另一个待定参数，然后构建经验高斯函数；

步骤4，根据经验高斯函数，进行观测点信号协方差矩阵的构建，并通过平差计算求解后验欧拉旋转参数和后验信号值；

步骤5，引入自适应因子，进行自适应迭代确定出欧拉旋转参数的最佳估值和后验信号最佳估值；

步骤6，依据引入距离尺度因子和自适应因子后确定的最佳信号协方差矩阵和欧拉旋转参数最佳估值，可确定出推估点与观测点间的信号协方差矩阵，从而推估出研究域任意点的信号值。

进一步地，所述的先验整体运动趋势的求解方程为：

上式中，vo为研究域内所有gps监测站点水平运动速度场构成的速度矩阵，即观测值；n为观测误差噪声，ω为欧拉旋转参数，a为ω系数阵；

所述的验先信号值的获得方法为：

根据所述的先验整体运动趋势，利用欧拉方程可计算获得研究域整体运动速度场，将vo扣除所述的研究域整体运动速度场，即可获得研究域地壳运动观测信号，即验先信号值。

进一步地，所述的根据所述验先信号值，构建拟合经验高斯函数所需的待定参数，公式为：

其中，c(0)为待定参数，q为参与建模的监测站点个数，为gps监测站点i的观测速度剔除趋势项后的验先信号值。

进一步地，所述的距离尺度因子的计算公式如下：

上式中，r为距离尺度因子，k为经验高斯参数中另一个待定参数，dmax表示研究域内相距最远的两个gps监测站点间的距离；

构建的经验高斯函数如下：

f(d)＝c(0)exp(-k²d²)

其中d为gps监测站点间的距离。

进一步地，所述的观测点信号协方差矩阵的构建公式为：

cij＝f(dij)

上式中，dij是观测点i和观测点j的球面距离；

对构建出的协方差矩阵求逆可得到初始信号权阵ps⁽⁰⁾；约定观测噪声和信号的初始先验单位权方差都为1，即则有：

式中，coo为速度矩阵vo的协方差矩阵，cnn为vo观测误差自协方差矩阵，cuo为推估点信号与观测点信号之间的协方差矩阵。

进一步地，引入自适应因子，进行自适应迭代确定出欧拉旋转参数的最佳估值和后验信号最佳估值，包括：

建立基于helmert方差分量的自适应拟合推估的函数模型，然后依据参数平差原理，组建法方程；根据法方程构造自适应因子，然后用自适应因子对信号权阵进行调整，得到新的信号权阵；

重新进行平差计算，同时求解新的后验欧拉矢量参数和后验信号值，然后进行迭代计算，经过λ次迭代后，当拟合点残差中误差达到最小时，此时认为观测噪声与信号方差分量协调一致，由此最终得到欧拉旋转参数最佳估值后验信号最佳估值。

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

1.本发明针对建立距离与协方差出现的负定问题，引入距离尺度因子来建立高斯函数，提高数据利用率，有效避免了互协方差负值对构建经验协方差模型所造成的影响。

2.本发明针对观测信息和先验信息对模型参数的影响不平衡问题，引进自适应因子进行调整，顾及距离尺度因子后所建立的信号协方差阵与观测值信号协方差阵的协调性较常规算法有所提高，但在此基础上进一步加入自适应调整过程，在此过程中自适应因子的物理意义更加明确，可有效避免纯数学算法计算时造成的信号失真。

3.本发明提出的顾及距离尺度因子与自适应调整相融合的最小二乘拟合推估算法，可更准确地拟合推估反映出区域地壳构造运动特性，可为后续进一步研究区域构造动力学特征提供有价值的基础数据与参考。

附图说明

图1为本发明方法的流程图。

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

图3分别为常规最小二乘拟合推估法、单一顾及距离尺度因子最小二乘拟合推估法、单一自适应最小二乘拟合推估法、顾及距离尺度因子的自适应融合算法对青藏高原地壳运动高精度gps水平速度场拟合残差的统计直方图。

具体实施方式

在拟合推估gps水平地壳速度场研究时，可同时顾及趋势性与随机性的最小二乘拟合推估法近些年来得到了越来越广泛的应用。最小二乘拟合推估法可同时顾及研究域地壳运动整体趋势与信号随机性，保证了未知量具备随机性和解析性。该方法可顾及先验信号信息，在理论上最优且拟合推估值可靠，但该算法的核心是如何精确确定出用于构造协方差阵的协方差函数(武艳强，2009；江在森，2010)。针对最小二乘拟合推估法的这一关键核心问题，相关学者也提出了相应策略：如考虑研究域尺度范围与地壳活动的复杂性，拟合各向异性协方差函数以保证随机信号的各向异性等；考虑观测数据复杂性和各类数据的协调一致性问题，利用方差分量估计法调整观测信息和信号的先验权比等。然而，上述策略均只涉及或针对某一因素进行改进，特别是未同时顾及到研究域范围尺度和不同协方差矩阵关系不合理因素对构建合理经验协方差函数的综合影响，因此会导致速度场拟合推估结果不能达到最优。

因此，本发明兼顾考虑研究域范围尺度与不同协方差矩阵关系不合理因素，引入距离尺度因子构建更合理的协方差函数和协方差阵，在此基础上进一步利用方差分量估计以达到观测值和信号先验权比的协调一致。

本发明在对地壳运动gps水平速度场进行最小二乘拟合推估时，针对最小二乘拟合推估中存在的两个关键问题：

(1)在拟合经验高斯函数构建协方差矩阵时，通用方法对负值协方差多采取删减回避措施，降低了数据利用率导致拟合精度不高。

(2)观测信息和先验信息对模型参数解算影响不平衡，即观测值权阵与信号值的经验协方差矩阵方差因子没能协调一致，导致解算结果与实际不符。

针对上述问题，本发明首先引入距离尺度因子构建高斯经验协方差函数，解决数据利用率问题，有效避免了互协方差负值对构建经验协方差模型所造成的影响；在此基础上，再通过进一步引入自适应因子来调整观测信息和先验信息对模型参数的影响不平衡性，保证观测值权阵与信号值的经验协方差矩阵方差因子协调一致。新算法较常规以及仅顾及单一因素改进的最小二乘拟合推估法，显著提高了gps水平速度场的拟合精度，可更加准确地对区域gps速度场进行拟合推估，从而能更加准确地拟合推估反映出区域地壳构造运动特性，可为深入研究区域现今地壳构造运动状态及其内部机制特征提供十分重要且有价值的基础参考。本发明具体技术方案如下：

一种地壳运动gps水平速度场自适应最小二乘拟合推估算法，包括以下步骤：

步骤1，基于观测得到的研究域内所有gps监测站点水平运动速度场，依据最小二乘原理，求解出研究域的先验整体运动趋势；

根据欧拉方程计算研究域整体运动速度场，进而得到研究域地壳运动观测信号，即验先信号值；

其中，所述的先验整体运动趋势的求解方程为：

上式中，vo为研究域内所有gps监测站点水平运动速度场构成的速度矩阵，即观测值；n为观测误差噪声，ω为欧拉旋转参数，a为ω系数阵，可由gps监测站点坐标确定，因此aω即表示研究域的先验整体运动趋势(这里不考虑随机信号值或将其归结为观测噪声)。在上式中，q×1、q×3、3×1、q×1分别表示矩阵vo、a、ω和n的维数。

其中，所述的验先信号值的获得方法为：

根据所述的先验整体运动趋势，利用欧拉方程可计算获得研究域整体运动速度场，将所述的研究域内所有gps监测站点水平运动速度场构成的速度矩阵vo扣除所述的研究域整体运动速度场，即可获得研究域地壳运动观测信号，即验先信号值，用于后续构造经验高斯函数。

步骤2，根据所述验先信号值，根据下式构建拟合经验高斯函数所需的待定参数：

其中，c(0)为gps监测站点自协方差，即经验高斯参数中的一个待定参数；q为参与建模的监测站点个数，含义同步骤1中的q；为gps监测站点i、j的观测速度剔除趋势项后的验先信号值，其计算表达式为：vs＝vo-aω，分别为vs中的第i、j行；qdi表示在di距离段范围内所包含的gps监测站点对的个数，c(di)为划分距离段后di距离段对应的互协方差。

下面的步骤开始进行经验高斯函数的拟合，本发明提出的距离尺度因子和自适应因子也是从该步骤分步引入的。

步骤3，引入顾及研究域范围的距离尺度因子，通过选取不同的距离尺度因子，确定经验高斯函数中另一个待定参数，然后构建经验高斯函数；

其中，r的计算公式如下：

上式中，k为经验高斯参数中另一个待定参数，dmax表示研究域内相距最远的两个gps监测站点间的距离。在gps监测站点间距离参数d取dmax时r趋近于无穷小的准则下通过对r取不同的值，确定相应的k，进而根据顾及距离尺度因子下的经验高斯函数表达式确定经验高斯函数，表达式如下：

f(d)＝c(0)exp(-k²d²)

其中d为gps监测站点间的球面距离。通过对拟合结果残差中误差进行分析，选取其中残差中误差较小k值，并顾及研究域范围大小，适当调整r，通过k控制高斯函数和研究域尺度更加匹配，以获得更好的拟合结果，确定引入距离尺度因子后的经验高斯函数具体数学表达式。

而在利用常规最小二乘拟合推估在此步中拟合高斯函数时，通常是根据步骤2中计算得到的c(0)，c(di)，将c(di)中的负值剔除，以保证用于拟合的数据同高斯函数正定性一致，然后再以各距离段为横轴，以每距离段站点间互协方差为纵轴(0距离段对应的站点间互协方差即为站点自协方差)，根据未顾及研究域距离尺度的经验高斯函数f(d)＝aexp(-k²d²)(常规最小二乘拟合推估的经验高斯函数系数a为最小二乘求得，而顾及距离尺度因子的经验高斯函数的系数c(0)为步骤2中直接算出)，用最小二乘分别拟合东、西向的高斯函数，分别得到常规最小二乘拟合推估法经验高斯函数的两个东向待定参数a，k和两个西向待定参数a，k，至此，在常规最小二乘拟合推估中，经验高斯函数的具体数学表达形式就认为是完全确定了，被用于后续各协方差阵的构建。

然而，上述常规最小二乘拟合推估算法中，对负值采取的删减与回避措施，会导致数据利用率大大降低，拟合模型可靠性也随之下降。本发明为了提高数据的利用率，引入顾及研究域范围的距离尺度因子r。在引入距离尺度因子r后，不仅顾及了研究域尺度范围，而且显著提高了数据的利用率，避免了互协方差负值对构建经验协方差模型所造成的影响，使经验高斯函数更加准确可靠。

下面通过步骤3构造的经验高斯函数，构造经验协方差矩阵以进行后续计算。

步骤4，根据经验高斯函数，进行观测点信号协方差矩阵的构建，并求解后验欧拉旋转参数和后验信号值；

协方差矩阵构建公式为：

cij＝f(dij)

上式中，dij是观测点i和观测点j的球面距离；

对构建出的协方差矩阵求逆可得到初始信号权阵ps⁽⁰⁾；约定观测噪声和信号的初始先验单位权方差都为1，即则有：

上式中，coo为速度矩阵vo的协方差矩阵，用来描述gps站间空间分布相关性与离散性；cnn为vo观测误差自协方差矩阵，用来描述gps站点速度自身的离散性，反映观测精度；cuo为推估点信号与观测点信号之间的协方差矩阵。上述三个协方差阵，cnn可根据gps速度场解算结果得出；coo、cuo根据引入距离尺度因子后构建的经验高斯函数确定。

通过上面三个公式平差计算可求解得到后验欧拉旋转参数和推估出的研究域地壳运动观测点后验信号值此时求解出的是自适应迭代前的，为便于与后续计算参数区别，故加上上缀0，将写为将写成

至此，已解出后验欧拉旋转参数和观测点后验信号值，常规最小二乘拟合推估法和仅顾及距离尺度因子的最小二乘拟合推估法在这一步已结束。但本发明考虑到仅顾及研究域范围尺度引入距离尺度因子后，并不能使构造的协方差阵直接达到最优，而不合理的协方差函数会造成观测值噪声和信号方差因子的不协调而产生系统误差，因此，本发明在估计距离尺度因子的常规最小二乘拟合推估的基础上再进一步融合引入自适应调整策略，使观测值噪声和信号方差因子达到协调一致，进一步提升对gps水平速度场拟合推估的精度；下面进行自适应调整。

步骤5，引入自适应因子，进行自适应迭代确定出欧拉旋转参数最佳估值和后验信号最佳估值，具体步骤如下：

建立基于helmert方差分量的自适应拟合推估的函数模型min：

min＝v^tpev+αs^tpss

上式中，α为自适应因子，v为观测值v0的改正数，s为信号值，pe为观测向量权阵，ps为信号权阵。该模型将信号视作虚拟观测值，与观测值信号一同形成两类观测量，可建立如下矩阵：

进一步地，依据参数平差原理，组建法方程：

上面的几个公式中，a的含义同步骤1，为欧拉刚性运动模型参数向量ω的系数阵，b为信号值的系数，其形如[i0]^t，假设推估点有n3个，则b中的i为n1×n1阶为单位阵，0为n3×n1阶的零矩阵；l的含义同步骤1中的vo，但因此处为误差方程形式，故用l表示；vs为信号值的改正数，此处vs即等于s；为观测信号单位权方差；σe²为观测噪声单位权方差；对于s和n矩阵而言，其下标i为数值1,2中的一个，n1＝tr(ps*ps^-1)，n2＝tr(pe*pe^-1)，j同i，即siisij对应上述法方程中的s11s12s21s22，根据上述各式可解出σe²和构造自适应因子α，然后利用自适应因子α对信号权阵ps进行升高或者降低调整，得到新的信号权阵ps⁽¹⁾：

ps⁽¹⁾＝αps

根据公式对新得到的信号值权阵ps⁽¹⁾进行求逆运算进而获得新的信号协方差阵后，再重新进行平差计算，同时求解新的后验欧拉矢量参数和后验信号值将步骤4中利用最小二乘拟合推估求解公式求解信号值和验后整体趋势项及步骤5中的一次方差分量估计重新确定权阵和σe²，作为一次迭代计算，经过λ次迭代后，当拟合点残差中误差达到最小时，此时认为观测噪声与信号方差分量协调一致，由此最终得到欧拉旋转参数最佳估值后验信号最佳估值

引进自适应因子进行调整后，顾及距离尺度因子后所建立的信号协方差阵与观测值信号协方差阵的协调性较常规算法进一步提高，在自适应调整过程中，有了针对于研究区域尺度范围的物理意义，可有效避免仅利用纯数学解算所造成的信号失真问题。融合距离尺度因子与自适应因子新算法，可更准确地拟合推估反映出区域地壳构造运动特性。

步骤6，依据步骤5顾及距离尺度因子和引入自适应因子后确定的最佳信号协方差矩阵和欧拉旋转参数最佳估值，可确定出推估点与观测点间的信号协方差矩阵，从而推估出研究域任意点的信号值：

式中，为的逆矩阵，coo为速度矩阵vo的协方差矩阵，用来描述gps站间空间分布相关性与离散性；cnn为vo观测误差自协方差矩阵，用来描述gps站点速度自身的离散性，反映观测精度；cuo为推估点信号与观测点信号之间的协方差矩阵。上述三个协方差阵，cnn可根据gps速度场解算结果得出；coo、cuo根据引入距离尺度因子后构建的经验高斯函数确定。表示推估出的待推估点的信号值，此处含义同步骤4中vo、a同步骤1中的vo、a，则表示欧拉旋转参数的最佳估值表示研究区域整体欧拉趋势项的最佳估值。

利用推估出的研究域任意点的信号值并结合研究域已有地质勘查资料，可进行区域地壳构造运动特性的研究。

实例：

为验证本发明提出算法的最优性及有效性，选取我国典型地壳构造强活动区域青藏高原区域，利用该区域地壳运动长期高精度gps水平速度场来检验本发明提出的最小二乘拟合推估新算法。gps观测数据来源于中国地壳运动观测网络和中国陆态观测网络，时间跨度为1999-2013年共672个站点数据，其中含有127个连续运行参考站和545个流动观测站，gps站点速度场信息如图2所示。分别利用常规、自适应、顾及距离尺度因子、以及融合顾及距离尺度因子与自适应的最小二乘拟合推估四种算法对青藏高原地壳运动gps水平速度场进行相应的拟合推估计算。本发明提出的融合距离尺度因子与自适应调整的最小二乘拟合推估法对青藏高原地壳运动gps水平速度场进行分析处理的方法如下：

(1)基于青藏高原地壳运动gps水平速度场，首先应用最小二乘法求解验前欧拉旋转参数。

(2)将上一步得到的青藏高原先验欧拉旋转参数，按公式vs＝vo-aω计算得到验前信号值，再利用具体实施步骤2中用验前信号值计算c(0)的公式，计算得到该研究域经验高斯函数中的参数c(0)。

(3)引入距离尺度因子r，挑选出研究域中相距最远的两站点的距离dmax，在确定dmax后实验选取不同的r确定该研究区域经验高斯函数中的参数k，进而确定经验高斯函数的具体数学表达式，用以后续构建经验协方差矩阵。

(4)由cij＝f(dij)构建研究域gps观测站点间观测速度向量vo的协方差coo，以及推估点信号与观测点信号之间的协方差阵cuo，进一步，根据步骤4中提到的常规最小二乘拟合推估解算公式进行后验整体趋势项和后验信号值的计算，即得到研究域在顾及距离尺度因子下的整体趋势项和测站点的位移信号值。

(5)引入自适应因子，进行研究域的观测值噪声和解出验后信号的自适应迭代计算。在每一次迭代中根据步骤5中helmert方差分量估计公式，解出σe²和构造自适应因子，然后利用自适应因子对信号权阵进行调整得到新的信号权阵ps⁽¹⁾，重新进行平差，同时求解新的后验欧拉旋转参数和后验信号值经过k次迭代，拟合点残差中误差最小时，认为观测噪声与信号方差分量协调一致，得到研究域欧拉旋转参数最佳估值和信号最佳估值

(6)将上述构造的cuo代入公式中，可计算得到研究域内任一推估点的信号值，再结合后验欧拉旋转参数解算出的块体整体运动趋势项，即可得到研究域任意区域内的地壳运动gps水平速度场。

分别利用四种算法对青藏高原地壳运动gps水平速度场进行相应的拟合推估计算，得到四种算法的拟合残差中误差，如表1所示。表1显示，改进算法的拟合精度显著高于其他三类算法。图3拟合残差统计直方图，也直观地显示出本文提出的改进算法具有较好拟合效果。

为进一步检验算法的有效性，在研究域内选取了点位精度较高、可靠性与稳定性较好且较均匀分布的20个gps监测站点作为检核点，检核点位置如图2所示，其余652个gps监测站点作为拟合点，检核点拟合残差精度结果如表1所示。表1所示结果也表明，本文提出的改进算法是四种方案中的最优算法，进一步验证了本文提出的该改进算法的有效性。

表1检核点拟合残差精度结果