水平张量重力梯度仪垂向运动误差实时补偿方法

李 达，赵 明，范士锋，李 中，李城锁，赵 琳

（1.哈尔滨工程大学 智能科学与工程学院，哈尔滨 150001；2.天津航海仪器研究所，天津 300131；3.中国人民解放军海军装备部，北京 100071；4.中国船舶航海保障技术实验室，天津 300131）

重力梯度信号是重力矢量的空间微分信号，与传统的重力异常相比，该信号拥有更高分辨率的地球重力场信息，在油气资源勘探、地球内部构造研究以及惯性重力组合导航等领域具有重要意义[1-3]。目前存在静电悬浮、冷原子干涉的超导等多种不同原理的重力梯度测量方案，但迄今仅有旋转加速度计原理的重力梯度测量方案实现了近地表动态测量工程应用[4]。按照测量分量的不同，该原理的重力梯度仪分为水平张量式和全张量式两种[5,6]。本文研究的重力梯度仪为水平张量式，能够同时测量Γuv和Γxy两种水平重力梯度张量分量。该原理的重力梯度仪将高分辨力加速度计敏感轴水平布放，用于抑制动态条件下载体垂向运动加速度产生的动态测量误差。但由于加速度计安装误差的存在，载体垂向运动会以一定比例进入重力梯度测量通道，而真实的重力梯度信号极其微小。为此，一方面通过精密加工装配尽可能减小安装误差[7]，另一方面，需要通过数据处理进一步抑制该项误差对重力梯度测量的影响[8]。

在重力梯度垂向运动测量误差数据处理方面，国外成熟商用重力梯度仪FALCONTM和Air-FTGTM的研究报告中，将垂向运动补偿作为重力梯度数据处理的关键技术列出[9]，但详细的数据处理方法并不对外公布。国内通常采用递推最小二乘法估计加速度计敏感轴垂向角度综合安装误差角，以此结合载体实测垂向加速度完成垂向运动误差补偿。该方法在短时间、海况平稳的情况下能够较好地完成垂向运动误差补偿，但随着测量时间的增加或者海况变化时，该方法对综合安装误差角的估计精度变差，从而影响补偿效果。本文根据船载实测重力梯度数据特征提出含有遗忘因子且具备自适应能力的垂向运动误差实时补偿方法，通过调整补偿参数的方法能够在不同测量工况下完成重力梯度垂向运动误差的实时补偿，为实现高精度动态重力梯度测量提供了一种数据实时处理方法。

1 垂向运动误差机理分析

1.1 测量原理

重力梯度仪样机由重力梯度敏感器和稳定平台组成，并配套相应的减震器、环境控制装置和显控机箱等辅助测量设备，其中敏感器与惯性稳定平台如图1所示，敏感器用于测量水平分量重力梯度信号，稳定平台则为其隔离载体角运动，并在动态环境下提供可靠的指向基准。重力梯度敏感器将四只高分辨力加速度计均匀布设在一个旋转圆盘上，对称安装的加速度计敏感轴指向相反，用于减弱载体线运动对重力梯度测量的影响，正交安装的加速度计则可以抑制旋转误差对重力梯度测量的影响，其原理图如图2 所示，测量原理公式[10]为：

图1 重力梯度仪样机组成Fig.1 Composition of gravity gradiometer

图2 重力梯度敏感器原理图Fig.2 Schematic diagram of gravity gradient sensor

式(1)中，ai(i=1,2,3,4)是第i只加速度计输出的加速度测量值，l是加速度计旋转半径，Γuv、Γxy是重力梯度张量水平分量，其单位为E（1E=10-9s-2），ω是旋转圆盘的旋转角速度，一般为

1.2 误差机理分析

重力梯度敏感器中的各加速度计敏感轴垂向角度安装误差如图3 所示。

图3 加速度计敏感轴垂向角度安装误差示意图Fig.3 Schematic diagram of vertical angle installation error of accelerometer sensitive axis

该安装误差会使动态条件下的载体垂向运动干扰重力梯度敏感器测量信号，记asum=(a1+a3) -(a2+a4)，则有：

式(2)中δasum,az是由载体垂向运动引起的重力梯度敏感器测量误差，az是载体垂向线加速度，βi(i=1,2,3,4)是第i只加速度计敏感轴的垂向角度安装误差。结合式(1)和式(2)可知，载体垂向运动中的2ω频率分量，与综合安装误差角(β1+β3) -(β2+β4)相耦合，形成重力梯度测量误差。因此在重力梯度数据处理中，需要利用重力梯度仪惯性稳定平台中垂向加速度计测量得到的载体垂向加速度，结合综合安装误差角，进行载体垂向运动误差补偿。对式(2)进行全微分，同时记βsum=(β1+β3) -(β2+β4)，得到：

式(3)中dδasum,az为δasum,az的全微分函数，是垂向运动补偿的残余误差，δaz和δβsum分别表示载体垂向加速度az和综合安装误差角βsum的估计误差。由于重力梯度信号通过机械旋转的方式调制在敏感器测量信号的2ω频率分量上，加速度计在2ω频率分量上的测量精度约为水平，通过对加速度计敏感轴精密装调，可将加速度计敏感轴垂向角度综合安装误差控制在0.5′以内，即由此可计算载体垂向加速度测量误差产生的重力梯度垂向运动补偿残余误差在0.5E 水平，可忽略不计。

某次重力梯度船载测量实验中，载体垂向加速度的时域图和频域图如图4 所示，其在2ω频率分量上的幅值约为量级，以敏感器旋转半径100 mm 和滤波时间100 s 计算，为使垂向运动补偿的残余误差控制在10E 水平，需要加速度计敏感轴垂向角度综合安装误差的估计精度δβsum控制在0.5″以内，这对于高精度重力梯度动态测量造成了严峻挑战。

图4 船载实验中载体垂向加速度的时域图和频域图Fig.4 Time-domain and frequency-domain plots of carrier vertical acceleration in shipboard experiments

2 综合安装误差角实时估计方法

综合安装误差的估计精度δβsum要求控制在0.5″以内，只能利用载体垂向运动与敏感器输出之间的相关性，实时估计该综合安装误差角。

2.1 回归方程构建

在式(1)的测量方程中考虑载体垂向加速度的影响，得：

通过递推最小二乘的方法求取加速度计敏感轴垂向角度综合安装误差，其递推公式为：

式中，K是修正系数矩阵，P是状态估计的均方误差阵，是状态的最优估计，I 是单位矩阵，R是正定的量测均方误差阵，即递推中的初值为：

利用载体垂向运动产生敏感器测量误差的解析关系建立回归方程，并通过递推最小二乘的方法实时估计综合安装误差角，从而结合载体垂向加速度完成重力梯度垂向运动补偿，提高重力梯度实时测量精度。但递推最小二乘方法是综合利用了历史所有量测值Zk，在经过较长时间递推后，Pk阵逐渐收敛，修正系数矩阵Kk逐渐减小，使得新的量测值对状态估计的修正作用逐渐减弱，致使状态的最优估计不能很好地跟踪状态实际变化。此外，量测均方误差阵R则与载体动态密切相关，同一数值无法满足全测量周期内不同的载体运动扰动下滤波要求。而R阵的设置误差会引起滤波精度降低，严重时还可导致滤波器发散，在实际工程应用中必须予以考虑。

2.2 含有遗忘因子与渐消因子的估计方法

针对上述问题，首先在递推公式中引入遗忘因子，在状态估计中刻意改变了递推误差与量测误差之间的权重，从而提高了状态估计对实际状态的跟踪能力。可将式(6)改写为：

式(8)中s是略小于1 的实数比例因子，称为遗忘因子。工程实现中可通过加速度计敏感轴垂向角度综合安装误差的短期变化规律确定遗忘因子s的数值。

此外，为使得状态估计递推方程具备对不同工况的适应能力，在滤波中增加Sage-Husa 自适应环节[10]。量测矩阵的预测均方误差矩阵可定义为：

由此得到R阵的表达式为：

为提高系统实时运算效率，将式(11)写成递推估计的形式，即：

式中初值λ0为1，c（0＜c＜ 1）是渐消因子，当k→∞时，有λk→ 1-c，此时可使滤波器始终拥有自适应能力。此外，通过序贯滤波的方法限制R阵的对角线元素数值，以保障R阵的正定性，满足滤波器的实现要求。记则：

图5 含有遗忘因子的Sage-Husa 自适应滤波算法流程图Fig.5 Flowchart of Sage-Husa adaptive filtering algorithm with forgetting factor

3 遗忘因子和渐消因子的数值确定

3.1 遗忘因子的数值确定

在递推最小二乘算法中引入遗忘因子s，能够改善状态估计的跟踪性能，但也以s的比例提高了新量测信息的精度，等效降低了状态估计算法的稳定性。因此，在遗忘因子s的参数选择中，既要考虑状态估计的跟踪性能，又要考虑算法的鲁棒性。本文通过半物理仿真的方法确定该参数的数值，其中综合安装误差角βsum在动态下的测量规律由静态下获取，并结合某次船载重力梯度测量中的实测运动数据，构建仿真数据集。仿真中选取船载实验实测数据海况为3 级，舰船排水量为3000 t，重力梯度测量过程中航速为10 knot。无需考虑综合安装误差角实时估计中R阵的自适应问题，因此在本次仿真中将R阵设置为固定值。对该敏感器仿真信号使用递推最小二乘法（s=1）以及含有遗忘因子典型值的递推最小二乘方法对综合安装误差角βsum进行状态估计，多次仿真后某次典型结果如图6 所示，待状态估计收敛后，通过计算均方根（RMS）的方法评价状态估计精度，相关统计结果如图7 所示。结合图6 和图7 可知，当遗忘因子s取0.995 时，对综合安装误差角βsum的估计精度最高，优于0.07″，因此在本研究中遗忘因子s取值为0.995。

图6 不同遗忘因子情况下综合安装误差角估计结果Fig.6 Estimation results of integrated installation error angle for different forgetting factor cases

3.2 渐消因子的数值确定

与遗忘因子的取值原则类似，渐消因子c的取值越小，对最新量测噪声变化的适应能力越强，但若渐消因子c的取值过小，噪声均方误差阵的估计结果会剧烈变化，从而导致状态估计精度降低，严重时会引起系统发散。一般来讲，渐消因子c的取值范围在0.9和0.999 之间。本文同样采用半物理仿真的方法确定渐消因子的数值，仿真数据的构建方法与遗忘因子仿真数据的构建方法一致，但为有针对性提高量测噪声的自适应能力，仿真中载体运动选取船载实验实测数据海况为3-5 级，其中两端数据部分海况为3 级，中间数据部分海况为5 级，遗忘因子s设置为0.995。

对该敏感器仿真信号使用不同渐消因子的最小二乘方法对综合安装误差角βsum进行状态估计，多次仿真后某次典型结果如图8 所示，待状态估计收敛后，通过计算均方根（RMS）的方法评价状态估计精度，相关统计结果如图9 所示。结合图8 和图9 可知，当渐消因子c取0.9385 时，对综合安装误差角βsum的估计精度最高，优于0.23″，因此在本研究中渐消因子c的数值取0.9385。

图8 不同渐消因子情况下综合安装误差角估计结果Fig.8 Estimation results of the integrated installation error angle for different asymptotic factors

图9 不同渐消因子情况下综合安装误差角估计误差Fig.9 Combined installation error angle estimation error with different asymptotic factors

4 实验验证

4.1 综合安装误差估计结果

项目组在2022 年在南海某海域开展船载水平张量重力梯度测量实验，实验现场如图10 所示。实验中有效重力梯度测量时间为8 天，其中第3-4 天海况为5 级，其余时间海况为3-4 级。分别使用递推最小二乘方法和本文提出的含有遗忘因子与自适应滤波的估计方法处理仪器实测数据，两种方法得到的综合安装误差角βsum估计值如图11 所示。

图10 船载水平张量重力梯度测量实验现场图Fig.10 Site map of shipboard horizontal tensor gravity gradient measurement experiment

图11 两种估计方法得到的 βsum估计结果Fig.11 βsumestimates obtained by both estimation methods

由图11 可知，采用递推最小二乘方法得到的βsum估计值无法反映βsum的时变特性，而采用遗忘因子与自适应滤波的估计方法得到的βsum估计值存在24 h的周期特性，能够反映βsum随温度的变化特性，表明该方法能够提高重力梯度仪对βsum的实时估计精度。

4.2 重力梯度测量结果及精度统计

分别使用递推最小二乘方法和本文提出的含有遗忘因子与自适应滤波的估计方法处理仪器实测数据，得到两路重力梯度测量信号重复线测量结果如图12和图13 所示。

图12 递推最小二乘方法下水平张量重力梯度测量结果Fig.12 Measurement results of horizontal tensor gravity gradient using recursive least squares method

图13 遗忘因子与自适应滤波的估计方法下两路重力梯度信号测量结果Fig.13 Measurements of two gravity gradient signals under the estimation method of forgetting factor and adaptive filtering

重力梯度动态测量通过内符合中误差的数值来评价测量精度，其计算公式为[12]：

式(16)中，εj是重力梯度j分量的内符合中误差，Γj1(i)和Γj2是两次重力梯度动态测量j分量在第i点的测量值，n是测量总数。不同处理方法重力梯度测量重复线内符合精度统计结果如表1 所示，结果表明本文提出的估计方法能够将内符合精度由30E@1km提高至15E@1km。

表1 不同补偿方法重力梯度测量内符合精度统计Tab.1 Statistics on the internal compliance accuracy of gravity gradient measurements with different compensation methods

5 结论

本文针对基于旋转加速度计原理的水平张量重力梯度仪在动态测量中加速度计敏感轴垂向综合安装误差角的时变特性，以及状态估计中易受到载体运动干扰的问题，提出了一种含有遗忘因子与自适应滤波的综合安装误差实时估计方法。通过在状态估计中增加遗忘因子的方式抑制综合安装误差角的时变特性，并在滤波中引入渐消因子提高该方法对载体动态的适应性，实现了对综合安装误差角进行次优估计，以此完成重力梯度信号的实时垂向运动补偿。船载重力梯度测量实验结果表明，与递推最小二乘方法相比，本文所提方法实时得到的综合安装误差角能够反映其随温度的变化特性，可将重力梯度船载测量内符合精度由30E@1km 提高至15E@1km，具有工程指导意义。