DZW重力仪非线性零漂的贝叶斯估计方法

陈善勇 杨 江,2,3

1 中国地震局地震研究所，武汉市洪山侧路40号，430071 2 武汉地震科学仪器研究院有限公司，武汉市洪山侧路40号，430071 3 湖北省重大工程地震监测与预警处置工程技术研究中心，湖北省咸宁市青龙路11号，437000

地球重力场的精准测量对计量科学、防震减灾、大地测量、地球物理等领域具有十分重要的科学意义。我国一直对重力相关领域投入巨大，也取得众多研究成果。尤其是高精密重力仪器的开发，如超导重力仪、小型绝对重力仪和海空重力仪等均取得突破性进展[1]。高精密重力仪器的发展对重力数据处理方法提出更高的精度需求。

DZW重力仪是我国自行研发的第一台高精度重力仪，分辨率可达到μGal级。由于其采用的是垂直悬挂弹性系统，弹簧张力的衰减以及未被补偿或屏蔽的外界作用会产生零漂[2]。DZW重力仪1 d的零漂有时可高达40 μGal，直接影响到重力仪的观测精度与准确性。

经典的线性零漂处理方法认为，零漂率在观测过程中为恒值，在1个月以上的长期实际外业测量中呈非线性[2]。因此，线性零漂处理方法在重力观测过程中的应用方案是以24 h作为时间节点进行分段线性拟合，得到各段的零漂率[3]。然而对于处理几个月以上的长期重力观测数据，该方法会耗费较多时间和物力。如果1 d数据中有温度、气压等其他非线性因素，也会导致经过线性零漂校正后观测值的精度难以保证。

文献[4-5]通过建立贝叶斯平差模型来解决非线性漂移观测误差，从而为提高数据精度提供新的解决方案。本文在其基础上进行优化尝试，贝叶斯平差模型中的研究对象为零漂率，本文贝叶斯模型将其变为零漂。这是由于零漂与重力观测值之间只涉及到加减运算，不涉及与时间的乘法运算，所需要的参量只有重力观测值以及理想的重力仪读数，因此可以简化运算并增加算法的计算效率。

为将本文算法运用到实际重力观测中，我们设计一套关于该算法的应用方案，使得该算法可以在每次获得新观测数据后自动进行计算，增加算法的实用性。该方案将算法运用到重力观测数据的实时处理中，在采样重力观测值的同时进行零漂计算，可以减少后期数据处理中消除零漂所消耗的人力和物力[6]。

1 模型构建

零漂的贝叶斯估计模型需要先定义与算法相关的矩阵以及基本公式，模型的核心公式为贝叶斯公式，需要分别确定先验分布形式、似然函数以及边缘密度函数。

1.1 矩阵定义

将经过潮汐校正的相对重力观测数据表示为：

x=[x1,x2,…,xn]T

将每个样本点对应的零漂表示为：

v=[v1,v2,…,vn]T

将每个样本点的重力仪读数(理想值)表示为：

g=[g1,g2,…,gn]T

将用于得到一组数据的梯度变化量的矩阵A(A∈Rn×n)表示为：

式中，ζ为方阵的补充值，取值趋近于0。

1.2 基本公式

对于经过潮汐校正后的相对重力观测数据，其与零漂之间的关系式可表示为[7]：

x-v-g=δ

(1)

式中，δ为公式计算产生的误差，满足以下分布：

δ～N(0,Σ1)

(2)

式中，Σ1为高斯分布的协方差矩阵。

由式(1)可推出以下分布：

x～N(v+g,Σ1)

(3)

1.3 先验分布确定

对零漂的先验分布形式进行确认，每个时刻的零漂是未知的，因此存在随机性[8-9]，可用随机变量表示：

(4)

假设零漂在连续时间中具有连续性[11]，则需要满足：

Av=ξ

(5)

式中，ξ为趋向于0的数组，满足以下分布：

ξ～N(0，Σ2)

(6)

式中，Σ2为高斯分布的协方差矩阵。

由式(5)可推出以下分布：

Av～N(0，Σ2)

(7)

已知v为高斯过程，根据高斯过程的线性性质，由式(7)可推算出v的具体分布为：

v～N(A-1·0,A-1Σ2(AT)-1)⟹

v～N(0,A-1Σ2(AT)-1)

(8)

先验分布的概率密度函数为：

(9)

1.4 似然函数确定

似然函数是样本信息的集中体现，由式(3)可推出其表达式为：

(10)

1.5 边缘密度函数确定

边缘密度函数本身并不会影响后验分布的分布形式，但会影响后验分布的幅值变化[12]，可通过式(9)、(10)推出边缘密度函数的表达式为：

(11)

1.6 后验分布计算

将先验分布的概率密度函数、似然函数和边缘密度函数代入贝叶斯公式，推出后验分布的概率密度函数：

(12)

正态均值(方差已知)的共轭先验分布为正态分布。后验分布满足以下分布形式：

(13)

1.7 超参数确定

由式(13)可知，要得到确定的后验分布，需要确定两组超参数值，分别是Σ1和Σ2。其中，Σ1为样本所产生的方差，代表重力仪的属性，视为恒定值[13]，可以取其中1 d的数据进行线性拟合得到零漂率，从而得到零漂校正后的数据方差；Σ2为零漂先验分布的超参数，可通过增加计算次数得到收敛极限进行确定，初始值可以取与Σ1同一量级的数值。

对于零漂后验分布均值和方差的确定，方法如下：

(14)

经过贝叶斯公式计算得到的零漂后验分布方差是两者的综合体现。将该后验分布作为下一次贝叶斯公式的先验分布，则得到的第二次后验分布的方差可表示为：

(15)

第n次后验分布的方差为：

(16)

由式(16)可知，后验分布方差会随着计算层数的增加而不断减小并趋向于零，表明后验分布具有收敛性，并且收敛极限为0。

后验分布的均值经过化简可表示为：

(17)

2 算法数据验证

本文算法在经过前文数学逻辑推导后，需要通过实际重力数据来进一步验证。首先，对算法中超参数进行初步估计，通过数据进一步验证零漂超参数的收敛性；然后，通过数据来确认均值的收敛极限性质。在算法本身逻辑性得到证明的基础上，需要确定应用方案来模拟测量过程，对所得结果的可用性和准确性进行分析；最后，通过与传统线性拟合去零漂方案进行对比分析，探讨算法是否实现去非线性零漂的目标。

2.1 超参数初步估计

将九峰台站DZW重力仪在2019年的连续观测数据作为算法的验证数据源。该数据已经过精度为10-5μGal的固体潮校正，可不考虑温度、气压等因素的影响。通过线性拟合，得到图1中表示线性零漂的直线，与潮汐校正后的数据点进行对比。拟合直线斜率即为估计的零漂率，为0.014 25 μGal/min=0.855 μGal/h，由给出的均方根误差RMSE(root mean squared error)推出样本的数据方差Σ1=RMSE2=0.184 92μGal2=0.034 19 μGal2，Σ1即为0.034 19的对角矩阵。由于零漂的先验分布的协方差矩阵Σ2为未知参数，但在多层运算后将会趋近于0，因此不妨设其为与Σ1同数量级的矩阵，令其为0.01的对角矩阵。

图1 对1 d重力观测数据进行线性拟合所得结果Fig.1 The results obtained by linear fitting of one-day gravity observation data

2.2 超参数收敛性验证

为确保计算次数足够多，将计算次数设置为1 000。1 000次计算中后验分布的方差变化如图2所示。

图2 后验分布方差变化Fig.2 Variance variation of posteriori distribution

在前100次计算中，方差下降速度较快；在后900次计算中，下降速度趋于平缓，总体呈收敛态势。为进一步验证均值的收敛性，需要研究均值的收敛方向和收敛极限。为遍历所有点的均值收敛情况，将计算次数减少至100。经过研究发现，每个样本点的收敛方向不同，收敛极限也存在差异。以第11个和第1 003个样本点的均值变化为例，结果见图3(a)和3(b)。

图3 样本点均值变化Fig.3 Meanvalues change of sample points

通过上述分析可知，随着计算次数不断增加，后验分布的方差趋近于0，各样本点的均值沿着各自方向趋近于不同的收敛极限。但收敛极限是否为零漂的真实值需进一步确认。

2.3 均值收敛极限性质确认

将重力观测值与后验分布均值作差，理论结果为真实重力观测值，是一个恒值，结果如图4所示。

图4 观测数据在去除不同计算次数所得零漂结果Fig.4 Drift results of observed data after removing different calculation times

图4(a)中方差为0.000 343 μGal2，图4(b)中方差为1.803 78×10-8μGal2，后者方差明显小于前者，说明随着计算次数的增加，均值不断逼近，使去零漂结果为常数数值，均值的收敛极限即为真实的零漂值。

2.4 算法应用方案确定

第1天的1 440个样本点采用贝叶斯模型进行预处理，通过对其进行1 000次重复计算，得到精准的1 440个零漂值。对于新样本数据，通过左移操作实现样本数组的更新，为确保样本点与零漂值一一对应，零漂的均值矩阵与协方差矩阵中的对角元素数组同时进行左移操作(图5)。

图5 算法应用方案Fig.5 Scheme of algorithm application

图5中v0作为零漂数组的输出数据按序存放，形成的数据组即为精准的零漂数组。xn+1为新样本值，vn+1代表与新样本值对应的均值，是需要预测的参数。Sn+1代表与新样本值对应的方差，等于先验分布的方差。该方案的优点是对于新移入的样本值，其均值由移入至最后输出会经历1 440次贝叶斯算法计算，可保证结果的高精度。

2.5 贝叶斯算法结果分析

贝叶斯算法得到连续2 d的零漂以及零漂校正后的重力观测值如图6所示。

图6 贝叶斯算法所得结果Fig.6 Bayesian algorithm results

由图6(b)可知，样本矩阵、零漂矩阵、协方差矩阵的不断变化，会使第2天的重力观测值方差更大，第1天重力观测值方差为8.458 01×10-12μGal2，而第2天为1.595 53×10-10μGal2，方差数量级已符合高精度的要求。

2.6 干扰作用下的结果分析

人为添加白噪声作为干扰信号加入到重力观测值中，模拟在观测过程中遇到外部干扰的情况。通过算法得到连续2 d的零漂结果以及零漂校正后的重力观测值如图7所示。

图7 加入干扰后的算法结果Fig.7 Algorithm results after adding interference

由图7(a)可知，干扰信号会对零漂结果造成一定影响，影响在可控范围内的原因在于先验分布的协方差矩阵会约束零漂的变动幅度。由图7(b)可知，经过零漂校正后的重力观测值仍会将干扰信号特征体现出来，同时也说明算法具有一定的抗干扰能力，在干扰结束后仍然可以精准地计算零漂值。

2.7 线性掉格下的结果分析

人为在重力观测数据中添加线性参数，模拟重力仪由于人员操作或温度骤变等外界因素而使仪器出现线性掉格的情况，结果如图8(a)所示。

由图8(b)可知，重力仪出现掉格后，该算法的去零漂结果会出现微小跳变，经过300个数据点后零漂校正结果回归正常，说明该算法需要300个数据点进行恢复零漂计算，但由于跳动较小，因此可以视为在出现线性掉格时，该算法仍然可以正常工作并保持准确度。

图8 线性掉格情况下算法结果Fig.8 Algorithm results in case of linear drifting

2.8 经典算法结果对比

通过线性拟合对同样2 d的连续重力观测值进行零漂估计，拟合结果如图9所示。

图9 连续2 d重力观测值线性拟合结果Fig.9 Linear fitting results of gravity observations for two consecutive days

线性拟合的零漂率为0.014 02 μGal/min ，由此得到去除零漂的重力观测值如图10所示。

图10 线性拟合算法零漂校正Fig.10 Drift correction of linear fitting algorithm

对比图10与图6(b)发现，贝叶斯算法可将重力观测值中的非线性零漂去除，实现对非线性零漂的精准估计。

3 结 语

本文提出一种基于贝叶斯公式，以零漂的连续性作为先验约束条件，考虑重力观测数据浮动产生的偏差，将每个采样点零漂的均值和方差作为超参数进行求解，从而获得零漂估计值的算法。通过零漂与重力观测数据的关系式以及零漂的连续性得到零漂的分布规律，采用逻辑推导和数据验证证明，贝叶斯算法得到零漂后验分布的均值会随着计算次数的增加而向零漂真实值收敛。在模拟的数据采样过程中，通过贝叶斯算法得到零漂校正结果的方差为1.595 53×10-10μGal2，符合DZW重力仪的精度需求。在添加人为干扰后得到的零漂校正结果表明，该算法可以在保留重力信号特征的前提下精准去除零漂值。

本文提出的非线性零漂贝叶斯估计方法能有效地去除重力观测值中的非线性零漂，可以在重力分采样过程中实时对重力观测值进行预处理，从而为重力观测精度提供保障。