CN113885099B - 重力梯度仪加速度计标度因数不一致性动态实时估计方法 - Google Patents

Abstract

本发明涉及一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，包括以下步骤：步骤1、在动态条件下，对重力梯度仪输出信号在旋转频率ω上进行解调，得到重力梯度仪输出信号的一倍频分量D₁和D₂；步骤2、基于步骤1获得的重力梯度仪输出信号的一倍频分量，在不一致信息提取中引入Kalman滤波环节，对两路信息实现最优估计；步骤3、在步骤2中引入自适应滤波的方法，得到两路对称安装的加速度计标度因数不一致信号的实时最优估计。本发明能够提高加速度计标度因数动态调整精度，进而提高旋转加速度计式重力梯度仪动态测量精度。

Description

重力梯度仪加速度计标度因数不一致性动态实时估计方法

技术领域

本发明属于重力梯度仪技术领域，尤其是一种重力梯度仪加速度计标度因数不一致性动态实时估计方法。

背景技术

重力梯度是重力矢量的空间变化率，在矿产资源勘探、地球科学研究以及惯性导航等方面具有重要意义。重力梯度仪是用于测量重力梯度的精密设备，基于BellAerospace公司提出的旋转加速度计测量原理的重力梯度仪是迄今唯一实用的近地表动态重力梯度仪。如图1所示，作为核心敏感器的重力梯度测量组件基于加速度计位置差分测量原理，通过机械旋转的方式将重力梯度张量分量调制到系统旋转频率的二倍频处，重力梯度敏感器输出与重力梯度张量分量之间的关系可表示为：

(a₁+a₃)-(a₂+a₄)＝4R(Γ_uvsin2ωt+Γ_xycos2ωt) (1)式中a₁、a₂、a₃和a₄是四只加速度计的测量输出信号，R是加速度计检测质心到旋转中心的距离，Γ_uv和Γ_xy是对应方向上重力梯度张量分量(其中ω是旋转机构的旋转角速度。将敏感器输出的加速度计组合信号进行2ω频率上的同步解调得到最终的重力梯度张量信号Γ_uv和Γ_xy。为便于书写，记a_sum＝(a₁+a₃)-(a₂+a₄)，即用a_sum表示重力梯度敏感器测量输出信号。

根据旋转加速度计式重力梯度测量原理，该原理的重力梯度仪通过对称反向安装的加速度计抵御近地表水平线运动对重力梯度测量的影响，但近地表测量载体运动加速度通常比重力梯度引起的运动对重力梯度测量的影响，但近地表测量载体运动加速度通常比重力梯度引起的加速度差值信号大6个数量级，加速度计间标度因数的微小差异会使得载体残余水平加速度进入重力梯度敏感器输出信号，形成测量误差。为实现高精度重力梯度动态测量，需要加速度计间标度因数不一致性优于10^-7量级，对实时调整精度提出了严苛的要求。

在工程实现中，旋转加速度计式重力梯度仪通过设置加速度计标度因数一致性调整回路实现以实现对称安装的加速度计标度因数高度一致性。具体实现步骤为：(1)在加速度计表体磁路中设置调整机构，通过控制调整电流实现标度因数在线调整功能；(2)由技术设计和工艺控制保证配对的加速度计标度因数初始差异不大于1‰；(3)倾斜敏感器承受重力分量形成常值分量，激励两组加速度计标度因数不一致信息调制在旋转频率ω的正余弦分量上，再通过同步解调提取；(4)对标度因数不一致信息进行PID校正，控制加速度计表体调整电流，实现加速度计标度因数实时调整功能。加速度计标度因数一致性调整回路原理如图2所示，步骤3的原理公式为：

式中δa_sum是因四只加速度计标度因数不一致引起的重力梯度仪输出信号误差，K₁、K₂、K₃和K₄是四只加速度计的标度因数，g是当前测量位置的重力加速度，η重力梯度敏感器倾斜角度，a_x和a_y是重力梯度仪测量坐标下两个水平方向的运动加速度。为便于书写，记两组对称安装的加速度计标度因数不一致信号为S₁和S₂，即S₁＝K₁-K₃和S₂＝K₂-K₄。

鉴于加速度计标度因数在温控条件下变化缓慢，通过合理设置PID参数能够使实际加速度计标度因数不一致信号跟踪估计的不一致信号。因此，动态下影响加速度计标度因数一致性调整回路精度的主要因素是如何在动态运动干扰下估计对称安装的加速度计不一致信号的数值。为此，诸多研究人员提出了不同的加速度计标度因数调整方法，但绝对多数方案只是针对静态条件下的调整方法，在动态条件下如何应对动态干扰的问题缺乏深入研究。

因此，亟需提出一种动态加速度计标度因数不一致信号估计方法，提高加速度计标度因数动态调整精度，进而提高旋转加速度计式重力梯度仪动态测量精度。

经检索未发现和本发明相同或相似的现有技术的公开文献。

发明内容

本发明的目的在于克服现有技术的不足，提出一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，能够提高加速度计标度因数动态调整精度，进而提高旋转加速度计式重力梯度仪动态测量精度。

本发明解决其现实问题是采取以下技术方案实现的：

一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，包括以下步骤：

步骤1、在动态条件下，对重力梯度仪输出信号在旋转频率ω上进行解调，得到重力梯度仪输出信号的一倍频分量D₁和D₂；

步骤2、基于步骤1获得的重力梯度仪输出信号的一倍频分量，在不一致信息提取中引入Kalman滤波环节，对两路信息实现最优估计；

步骤3、在步骤2中引入自适应滤波的方法，得到两路对称安装的加速度计标度因数不一致信号的实时最优估计。

而且，所述步骤1的具体公式为：

式中，D₁和D₂是重力梯度仪输出信号的旋转频率ω正弦分量和余弦分量，F()表示对括号内的信号进行低通滤波。

而且，所述步骤2的具体步骤包括：

(1)系统的状态转移方程离散形式为：

X_k＝X_k-1+W_k-1 (4)

式中，X＝[S₁ S₂]^T，两路对称安装的加速度计标度因数不一致信号的向量形式，W是2维系统噪声列向量，下角标k和表示k时刻；

(2)以重力梯度仪旋转频率ω解调值作为观测量，系统的量测方程离散形式为：

Z_k＝H_kX_k+V_k (5)

式中Z＝[D₁ D₂]^T，V是2维量测噪声向量，量测矩阵H为：

(3)由此得到Kalman滤波过程为：

①状态一步预测：

②状态一步预测均方误差阵：

P_k/k-1＝P_k-1+Q_k-1 (8)

③滤波增益计算：

④状态估计：

⑤状态估计均方误差阵：

P_k＝(I-K_kH_k)P_k/k-1 (11)

式中，P是状态估计向量的均方误差阵，Q是系统噪声均方误差阵，是系统噪声向量W的方差值，K是滤波增益矩阵，R是量测噪声均方误差阵，是系统噪声向量V的方差值，是当前时刻最优估计结果，其他变量均为中间计算过程变量。

而且，所述步骤3的具体步骤包括：

(1)定义量测预测误差为：

对式(12)求方差，则有：

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

将式(14)写成递推估计的形式，即：

式中，初值可根据载体运动情况设为略大的方差阵。

(2)当k→∞时，有1/k→0，将等加权平均改为指数渐消记忆加权平均，即：

式中初值λ₀＝1，c是渐消因子，要求0＜c＜1,；当k→∞时有λ_k→1-c，使滤波器始终保持R阵的自适应能力；对R阵的对角线元素进行限制，简记则：

通过上述方法可将始终限制在[R_min,R_max]之内；式中R_min和R_max的具体数值可分别在最佳测量工况和最恶劣测量工况下的实测数据计算获得。

(3)由此得到两路对称安装的加速度计标度因数不一致信号的实时最优估计

本发明的优点和有益效果：

1、本发明提出一种旋转加速度计式重力梯度敏感器加速度计标度一致性动态实时估计方法，能够实时适应载体动态对加速度计标度因数不一致性估计的干扰，提高动态下加速度计标度因数不一致性估计精度，进而提高重力梯度动态测量精度。

2、本发明在Kalman滤波过程中增加了利用量测信息与一步预测信息的差异实时修正量测信息均方误差阵的环节，以适应不同载体的动态工况，能够提高重力梯度仪对加速度计标度因数不一致信号的动态估计精度，进而提高重力梯度动态测量精度。

附图说明

图1为旋转加速度计式重力梯度仪测量原理图；

图2为加速度计标度因数一致性调整回路原理图；

图3为本发明的加速度计标度一致性动态实时估计算法流程图。

具体实施方式

以下结合附图对本发明实施例作进一步详述：

一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，包括以下步骤：

步骤1、在动态条件下，对重力梯度仪输出信号在旋转频率ω上进行解调，得到重力梯度仪输出信号的一倍频分量D₁和D₂；

所述步骤1的具体公式为：

式中，D₁和D₂是重力梯度仪输出信号的旋转频率ω正弦分量和余弦分量，F()表示对括号内的信号进行低通滤波。

步骤2、为提高动态下两路加速度计标度因数不一致信息提取精度，基于步骤1获得的重力梯度仪输出信号的一倍频分量，在不一致信息提取中引入Kalman滤波环节，对两路信息实现最优估计；

所述步骤2的具体步骤包括：

(1)系统的状态转移方程离散形式为：

X_k＝X_k-1+W_k-1 (4)

式中，X＝[S₁ S₂]^T，两路对称安装的加速度计标度因数不一致信号的向量形式，W是2维系统噪声列向量，下角标k和表示k时刻；

(2)以重力梯度仪旋转频率ω解调值作为观测量，系统的量测方程离散形式为：

Z_k＝H_kX_k+V_k (5)

式中Z＝[D₁ D₂]^T，V是2维量测噪声向量，量测矩阵H为：

(3)由此得到Kalman滤波过程为：

①状态一步预测：

②状态一步预测均方误差阵：

P_k/k-1＝P_k-1+Q_k-1 (8)

③滤波增益计算：

④状态估计：

⑤状态估计均方误差阵：

P_k＝(I-K_kH_k)P_k/k-1 (11)

式中P是状态估计向量的均方误差阵，Q是系统噪声均方误差阵，是系统噪声向量W的方差值，K是滤波增益矩阵，R是量测噪声均方误差阵，是系统噪声向量V的方差值，是当前时刻最优估计结果，其他变量均为中间计算过程变量。

在Kalman滤波参数设置中，系统噪声均方误差阵Q的数值可通过重力梯度仪静态数据计算得到，但量测噪声均方误差阵R则与载体动态息息相关，很难通过事先预设的数值满足同一种载体不同工况下的滤波要求。而R阵的设置误差会导致Kalman滤波的精度降低，严重时还可能会引起滤波发散，在实际工程应用中必须予以考虑。

步骤3、为解决上述问题，在步骤2中引入自适应滤波的方法，得到两路对称安装的加速度计标度因数不一致信号的实时最优估计。

所述步骤3的具体步骤包括：

定义量测预测误差为：

对式(12)求方差，则有：

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

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

式中，初值可根据载体运动情况设为略大的方差阵。

当k→∞时有1/k→0，即长时间滤波后自适应能力将逐渐减弱，直至几乎失去自适应能力，为始终保持R阵的自适应能力，将等加权平均改为指数渐消记忆加权平均，即：

式中初值λ₀＝1，c是渐消因子，要求0＜c＜1。当k→∞时有λ_k→1-c，使滤波器始终保持R阵的自适应能力。此外，为保证R阵的正定性，对R阵的对角线元素进行限制，简记则：

通过上述方法可将始终限制在[R_min,R_max]之内，从而保证滤波器具有良好的自适应能力和滤波稳定性。式中R_min和R_max的具体数值可分别在最佳测量工况和最恶劣测量工况下的实测数据计算获得，具体滤波算法流程如图3所示。

由此得到两路对称安装的加速度计标度因数不一致信号的实时最优估计

本发明的工作原理是：

旋转加速度计式重力梯度仪通过对加速度计标度因数实时调整的方法，实现动态条件下对称安装的加速度计标度因数的高度一致性，以抵御水平线运动对重力梯度测量的影响，为提高动态下加速度计标度因数调整精度，需要提出一种加速度计标度因数一致性动态估计方法。

本发明通过对重力梯度仪输出信号解调，得到两组加速度计标度因数不一致信号的观测值，结合加速度计标度因数变化缓慢的特点建立Kalman滤波方程，并利用量测信息与一步预测信息的差异实时修正量测信息的均方误差阵，通过递推的方式逐步修改R阵以适应不同载体动态，提高系统对加速度计标度因数不一致信号的动态估计精度。

需要强调的是，本发明所述实施例是说明性的，而不是限定性的，因此本发明包括并不限于具体实施方式中所述实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。

Claims (3)

1.一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，其特征在于：包括以下步骤：

步骤1、在动态条件下，对重力梯度仪输出信号在旋转频率ω上进行解调，得到重力梯度仪输出信号的一倍频分量D₁和D₂；

步骤2、基于步骤1获得的重力梯度仪输出信号的一倍频分量，在不一致信息提取中引入Kalman滤波环节，对两路信息实现最优估计；

步骤3、在步骤2中引入自适应滤波的方法，得到两路对称安装的加速度计标度因数不一致信号的实时最优估计；

所述步骤3的具体步骤包括：

(1)定义量测预测误差为：

对式(10)求方差，则有：

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

将式(12)写成递推估计的形式，即：

式中，初值可根据载体运动情况设为略大的方差阵；

(2)当k→∞时，有1/k→0，将等加权平均改为指数渐消记忆加权平均，即：

式中初值λ₀＝1，c是渐消因子，要求0＜c＜1,；当k→∞时有λ_k→1-c，使滤波器始终保持R阵的自适应能力；对R阵的对角线元素进行限制，简记则：

通过上述方法可将始终限制在[R_min,R_max]之内；式中R_min和R_max的具体数值可分别在最佳测量工况和最恶劣测量工况下的实测数据计算获得；

(3)由此得到两路对称安装的加速度计标度因数不一致信号的实时最优估计

2.根据权利要求1所述的一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，其特征在于：所述步骤1的具体公式为：

式中，D₁和D₂是重力梯度仪输出信号的旋转频率ω正弦分量和余弦分量，F()表示对括号内的信号进行低通滤波。

3.根据权利要求1所述的一种重力梯度仪加速度计标度因数不一致性动态实时估计方法，其特征在于：所述步骤2的具体步骤包括：

(1)系统的状态转移方程离散形式为：

X_k＝X_k-1+W_k-1 (2)

式中，X＝[S₁ S₂]^T，两路对称安装的加速度计标度因数不一致信号的向量形式，W是2维系统噪声列向量，下角标k和表示k时刻；

(2)以重力梯度仪旋转频率ω解调值作为观测量，系统的量测方程离散形式为：

Z_k＝H_kX_k+V_k (3)

式中Z＝[D₁ D₂]^T，V是2维量测噪声向量，量测矩阵H为：

(3)由此得到Kalman滤波过程为：

①状态一步预测：

②状态一步预测均方误差阵：

P_k/k-1＝P_k-1+Q_k-1 (6)

③滤波增益计算：

④状态估计：

⑤状态估计均方误差阵：

P_k＝(I-K_kH_k)P_k/k-1 (9)

式中，P是状态估计向量的均方误差阵，Q是系统噪声均方误差阵，是系统噪声向量W的方差值，K是滤波增益矩阵，R是量测噪声均方误差阵，是系统噪声向量V的方差值，是当前时刻最优估计结果，其他变量均为中间计算过程变量。