CN112363247B - 一种重力梯度仪运动误差事后补偿方法 - Google Patents

一种重力梯度仪运动误差事后补偿方法 Download PDF

Info

Publication number
CN112363247B
CN112363247B CN202011163660.3A CN202011163660A CN112363247B CN 112363247 B CN112363247 B CN 112363247B CN 202011163660 A CN202011163660 A CN 202011163660A CN 112363247 B CN112363247 B CN 112363247B
Authority
CN
China
Prior art keywords
motion
gravity gradiometer
data
motion error
sin
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011163660.3A
Other languages
English (en)
Other versions
CN112363247A (zh
Inventor
喻名彪
涂良成
胡宸源
喻立
范继
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN202011163660.3A priority Critical patent/CN112363247B/zh
Publication of CN112363247A publication Critical patent/CN112363247A/zh
Application granted granted Critical
Publication of CN112363247B publication Critical patent/CN112363247B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V13/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices covered by groups G01V1/00 – G01V11/00
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Engineering & Computer Science (AREA)
  • Manufacturing & Machinery (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
  • Force Measurement Appropriate To Specific Purposes (AREA)

Abstract

本发明公开了一种重力梯度仪运动误差事后补偿方法，包括下述步骤：(1)获得重力梯度仪的输出数据、比力数据、角速度数据和角加速度数据，并对数据进行预处理，将数据分成N个数据块；(2)利用重力梯度仪运动误差模型对N个数据块进行运动误差补偿；(3)将补偿后的N个数据块合并，并进行解调或滤波获得万有引力梯度。本发明基于重力梯度仪运动误差模型，利用匀速勘探的情况下，重力梯度仪运动状态仅与运动误差相关和万有引力梯度响应无关的特性，构建损失函数，利用损失函数最小，估计运动误差传递系数，进而计算出运动误差，剔除重力梯度仪输出数据中的运动误差。本发明同时适用于调制类和非调制类重力梯度仪的运动误差补偿，且操作简单、易实施。

Description

一种重力梯度仪运动误差事后补偿方法
技术领域
本发明属于精密测量技术领域，更具体地，涉及一种重力梯度仪运动误差事后补偿方法。
背景技术
重力梯度仪通过差分万有引力加速度或万有引力力矩的方式，测量万有引力梯度。受制造工艺的限制，重力梯度仪存在许多非理想因素，比如用于差分的万有引力加速度传感器组，其输入输出特性存在差异、安装位置非对称、敏感轴方向存在失准等，使得差分过程中，非引力加速度不能完全扣除，引起测量误差。
重力梯度仪事后运动误差补偿，利用运动传感器监测动基座勘探中的梯度仪的运动状态，勘探结束后，根据运动状态数据及重力梯度仪运动误差模型，计算出运动误差并进行补偿。事后运动误差补偿，其补偿精度由运动误差模型精度及运动监测精度共同决定。当运动监测精度足够高时，事后运动误差补偿方法可工作的动态环境，由重力梯度仪运动误差模型的动态环境阈值决定。
目前文献所述的事后运动误差补偿方法，专利2018116511768及论文10.1109/TIM.2020.3010194，使用的重力梯度仪运动误差模型(18参数运动误差模型)，没有考虑线运动与角运动的耦合项及调制频率3倍频以上的项，忽略的误差项，使得运动误差模型精度随动态环境的增大变差，限制了运动误差模型的动态环境阈值，进而限制了事后运动误差补偿方法的环境适应性。
发明内容
针对现有技术的缺陷，本发明的目的在于提供一种重力梯度仪运动误差事后补偿方法，旨在解决现有技术中由于重力梯度仪运动误差模型没有考虑线运动与角运动的耦合项及调制频率3倍频以上的项，导致的重力梯度仪动态环境适应性差的问题。
本发明提供了一种重力梯度仪运动误差事后补偿方法，包括下述步骤：
(1)获得重力梯度仪的输出数据、比力数据、角速度数据和角加速度数据，对数据进行预处理后，将数据分成N个数据块；
(2)利用重力梯度仪运动误差模型对N个数据块的运动误差进行补偿；
(3)将补偿后的N个数据块合并，并进行解调或滤波后获得万有引力梯度。其中，对于调制类重力梯度仪，进行解调获取万有引力梯度，对于非调制类重力梯度仪，进行滤波，获得万有引力梯度。
更进一步地，采用的重力梯度仪运动误差模型，它的25参数模型为：
Figure BDA0002745115030000021
式中
Figure BDA0002745115030000022
是25参数运动误差模型t时刻的输出，m1×25(t)表示t时刻的运动矢量，C25×1表示运动误差传递系数矢量；m1×25(t)是包含25个元素的行向量：m1×25(t)＝[m1(t) m2(t) … m25(t)]；m1×25(t)的25个元素给出如下：
Figure BDA0002745115030000023
Figure BDA0002745115030000031
其中，Ω表示重力梯度仪旋转圆盘角频率，也就是重力梯度仪的调制频率，当重力梯度仪为非调制类重力梯度仪时，Ω的值为0；sin,cos分别表示正弦函数和余弦函数；ax(t)、ay(t)、az(t)分别表示t时刻重力梯度仪测量坐标系原点的比力矢量，ωx(t)、ωy(t)、ωz(t)分别表示t时刻万有引力传感器相对惯性系的角速度矢量，ωax(t)、ωay(t)、ωaz(t)分别表示t时刻万有引力传感器相对惯性系的角加速度矢量。
更进一步优选地，采用25参数运动误差模型对N个数据块的运动误差进行补偿步骤具体包括：
(2.1.1)根据t时刻重力梯度仪测量坐标系原点的比力矢量ax(t)、ay(t)、az(t)，t时刻万有引力传感器相对惯性系的角速度矢量ωx(t)、ωy(t)、ωz(t)，t时刻万有引力传感器相对惯性系的角加速度矢量ωax(t)、ωay(t)、ωaz(t)和所述25参数运动误差模型，计算t时刻的运动矢量m1×25(t)＝[m1(t) m2(t) … m25(t)]；
(2.1.2)构建第一损失函数并利用第一损失函数最小作为约束来标定重力梯度仪误差传递系数矢量；
(2.1.3)根据标定后的运动误差传递系数矢量来剔除重力梯度仪运动误差，实现运动误差的补偿。
更进一步地，第一损失函数
Figure BDA0002745115030000041
其中，t表示数据的时间，Gout(t)表示t时刻重力梯度仪的输出，C25×1表示运动误差传递系数矢量，是25行1列的矢量，m1×25(t)表示t时刻的运动矢量，是一个1行25列矩阵，t0表示数据起始时刻，tp表示数据结束时刻。
更进一步地，在步骤(2.1.3)中，根据公式
Figure BDA0002745115030000042
来剔除重力梯度仪运动误差，其中，Gout(t)为t时刻重力梯度仪的原始输出，
Figure BDA0002745115030000043
为剔除了运动误差后t时刻的重力梯度仪的输出。
更进一步地，采用的重力梯度仪运动误差模型，它的54参数模型为：
Figure BDA0002745115030000044
式中
Figure BDA0002745115030000045
表示54参数运动误差模型t时刻的输出，m1×54(t)表示t时刻的运动矢量，C54×1表示运动误差传递系数矢量；m1×54(t)是包含54个元素的行向量：m1×54(t)＝[m1(t)m2(t) … m54(t)]；m1×54(t)的54个元素给出如下：
Figure BDA0002745115030000051
Figure BDA0002745115030000061
Figure BDA0002745115030000071
其中，Ω表示重力梯度仪旋转圆盘角频率，也就是重力梯度仪的调制频率，sin,cos分别表示正弦函数和余弦函数；ax(t)、ay(t)、az(t)分别表示t时刻重力梯度仪测量坐标系原点的比力矢量，ωx(t)、ωy(t)、ωz(t)分别表示t时刻万有引力传感器相对惯性系的角速度矢量，ωax(t)、ωay(t)、ωaz(t)分别表示t时刻万有引力传感器相对惯性系的角加速度矢量。
进一步优选地，采用54参数运动误差模型对N个数据块的运动误差进行补偿步骤具体包括：
(2.2.1)根据t时刻重力梯度仪测量坐标系原点的比力矢量ax(t)、ay(t)、az(t)，t时刻万有引力传感器相对惯性系的角速度矢量ωx(t)、ωy(t)、ωz(t)，t时刻万有引力传感器相对惯性系的角加速度矢量ωax(t)、ωay(t)、ωaz(t)和所述54参数运动误差模型，计算t时刻的运动矢量m1×54(t)＝[m1(t) m2(t) … m54(t)]；
(2.2.2)构建第二损失函数并利用第二损失函数最小作为约束来标定重力梯度仪误差传递系数矢量；
(2.2.3)根据标定后的运动误差传递系数矢量来剔除重力梯度仪运动误差，实现运动误差的补偿。
更进一步地，第二损失函数为
Figure BDA0002745115030000081
其中，t为数据的时间，Gout(t)表示t时刻重力梯度仪的输出，C54×1表示运动误差传递系数矢量，是54行1列的矢量，m1×54(t)表示t时刻的运动矢量，是一个1行54列矩阵，t0表示数据起始时刻，tp表示数据结束时刻。
更进一步地，在步骤(2.2.3)中，根据公式
Figure BDA0002745115030000082
来剔除重力梯度仪运动误差，其中，Gout(t)为t时刻重力梯度仪的原始输出，
Figure BDA0002745115030000083
为剔除了运动误差后t时刻的重力梯度仪的输出。
通过本发明所构思的以上技术方案，与现有技术相比，由于重力梯度仪运动误差模型补充考虑了线运动与角运动的耦合项及调制频率3倍频以上的项，其运动误差模型的精度不会随着动态环境的增大，迅速降低，能够取得显著提高事后运动误差补偿方法动态环境适应性的有益效果。
附图说明
图1为重力梯度仪运动监测传感器安装示意图；
图2为本发明实施例提供的事后运动误差补偿方法的实现流程图；
图3为一种旋转加速度计重力梯度仪运动监测方案；
图4为仿真的重力梯度仪x方向比力监测数据；
图5为仿真的重力梯度仪z方向比力监测数据；
图6为仿真的重力梯度仪x方向角速度监测数据；
图7为仿真的重力梯度仪z方向角速度监测数据；
图8为仿真的重力梯度仪的原始输出，包含运动误差和万有引力梯度响应；
图9为仿真的54参数模型进行事后运动误差补偿后的输出与理论值的对比；
图10为仿真的25参数模型进行事后运动误差补偿后的输出与理论值的对比；
图11为仿真的18参数模型进行事后运动误差补偿后的输出与理论值的对比。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明提供了两个重力梯度仪运动误差模型，第一个重力梯度仪运动误差模型，包含了25个运动误差传递系数，考虑了重力梯度仪的线运动响应、角运动响应，以及线运动与角运动的主要耦合项；第二个重力梯度仪运动误差模型，包含54个运动误差传递系数，它相比第一个运动误差模型，考虑了更多的线运动与角运动的耦合项及调制频率3倍频以上的项，相比第一个模型，第二个模型的精度更高，能够应对更大的重力梯度勘探动态环境。
由于第一个重力梯度仪运动误差模型比第二个模型更加精简，在一般重力梯度勘探动态环境下，使用第一个重力梯度仪运动误差模型进行事后运动误差补偿，占用计算资源少，更加高效。当勘探动态环境恶劣时，使用第二个重力梯度仪运动误差模型进行运动误差补偿，保证事后运动误差补偿的精度。这两个重力梯度仪运动误差模型，相比现有技术中的重力梯度仪运动误差模型，考虑了角运动与线运动的耦合项，有更高的精度和更优的动态环境适用性。该事后运动误差补偿方法，同时适用调制和非调制类重力梯度仪。
如图1所示，在重力梯度仪测量坐标系安装陀螺仪和加速度计，重力梯度勘探时，用于记录重力梯度仪经历的线运动(比力)和角运动(角速度，角加速度)。
如图2所示，重力梯度仪运动误差事后补偿方法，包括以下步骤：
对重力梯度仪的输出数据、线运动数据、角运动数据作预处理：主要包括滤波处理，降低数据噪声，采样率转换，使重力梯度仪的输出数据、比力数据、角速度数据、角加速度数据的采样率相同；当梯度仪运动监测加速度计没有安装在重力梯度仪测量坐标系原点时，还需利用角运动数据，去除杆臂效应，得到重力梯度仪测量坐标系原点的比力数据。
在整个航空重力梯度勘探过程中，重力梯度仪运动误差系数向量会随着时间缓慢变化，但在较短的时间段内，可以视为不变；为了提高事后运动误差补偿的精度，可以根据重力梯度仪的特性，将重力梯度仪的输出数据、比力数据、角速度数据、角加速度数据，按时间分割成N个数据块。比如航空重力梯度勘探的总时间长度为6小时，所用的重力梯度仪的误差传递系数在1小时内，变化很小，可以视为不变，那么可以将重力梯度仪的输出数据、比力数据、角速度数据、角加速度数据分割成6个连续的数据块，每个数据块的时间长度为1小时。依次剔除6个数据块中的运动误差，每个数据块剔除运动误差的方法相同。
作为本发明的一个实施例，利用25参数运动误差模型对数据块进行事后运动误差补偿，具体步骤如下：
(a)将t时刻重力梯度仪测量坐标系原点的比力矢量ax(t)、ay(t)、az(t)，t时刻万有引力传感器相对惯性系的角速度矢量ωx(t)、ωy(t)、ωz(t)，t时刻万有引力传感器相对惯性系的角加速度矢量ωax(t)、ωay(t)、ωaz(t)代入下式中，计算t时刻的运动矢量m1×25(t)＝[m1(t) m2(t) … m25(t)]：
Figure BDA0002745115030000111
Figure BDA0002745115030000112
Figure BDA0002745115030000113
Figure BDA0002745115030000114
m5(t)＝-(ωay(t)+ωx(t)ωz(t))sin(Ωt)-(ωax(t)-ωy(t)ωz(t))cos(Ωt)
m6(t)＝-(ωax(t)-ωy(t)ωz(t))sin(Ωt)+(ωay(t)+ωx(t)ωz(t))cos(Ωt)
m7(t)＝(ωax(t)+ωy(t)ωz(t))sin(Ωt)-(ωay(t)-ωx(t)ωz(t))cos(Ωt)
m8(t)＝-(ωay(t)-ωx(t)ωz(t))sin(Ωt)-(ωax(t)+ωy(t)ωz(t))cos(Ωt)
m9(t)＝ax(t)sin(Ωt)-ay(t)cos(Ωt)
m10(t)＝ay(t)sin(Ωt)+ax(t)cos(Ωt)
m11(t)＝ax(t)az(t)sin(Ωt)-ay(t)az(t)cos(Ωt)
m12(t)＝ay(t)az(t)sin(Ωt)+ax(t)az(t)cos(Ωt)
m13(t)＝ax(t)ωaz(t)sin(Ωt)-ay(t)ωaz(t)cos(Ωt)
Figure BDA0002745115030000121
m15(t)＝ay(t)ωaz(t)sin(Ωt)+ax(t)ωaz(t)cos(Ωt)
Figure BDA0002745115030000122
Figure BDA0002745115030000123
m19＝az(t)；m20＝1
Figure BDA0002745115030000124
m22＝ωaz(t)；m23＝-az(t)ωaz(t)
Figure BDA0002745115030000125
Figure BDA0002745115030000126
上式中，Ω表示重力梯度仪旋转圆盘角频率，也就是重力梯度仪的调制频率，当重力梯度仪为非调制类重力梯度仪时，Ω的值为0；sin,cos分别表示正弦函数和余弦函数；将计算的运动矢量元素m1(t) m2(t) … m25(t)，构成t时刻的重力梯度仪运动矢量：m1×25(t)＝[m1(t) m2(t) … m25(t)]。
(b)构建第一损失函数J1(C)如下：
Figure BDA0002745115030000127
式中t表示数据块的时间，Gout(t)表示t时刻重力梯度仪的输出，C25×1表示运动误差传递系数矢量，它是25行1列的矢量，m1×25(t)表示t时刻的运动矢量，它是一个1行25列矩阵，t0表示数据起始时刻，tp表示数据结束时刻；
(c)利用第一损失函数J1(C)最小，标定重力梯度仪误差传递系数矢量C：
Figure BDA0002745115030000128
式中
Figure BDA0002745115030000129
是重力梯度仪误差传递系数矢量C的估计值，
Figure BDA00027451150300001210
表示使得函数J1(C)取最小时的C值；
(d)将标定的运动误差传递系数矢量
Figure BDA00027451150300001211
代入下式，剔除重力梯度仪运动误差：
Figure BDA00027451150300001212
式中，Gout(t)是t时刻重力梯度仪的原始输出，
Figure BDA0002745115030000131
是剔除了运动误差后t时刻的重力梯度仪的输出。最后将6个剔除了线运动误差、角运动误差的数据块合并，解调，得到6小时勘探的万有引力梯度输出。
相比专利2018116511768及论文10.1109/TIM.2020.301019418中的18参数运动误差模型，25参数运动误差模型考虑了线运动与角运动的主要耦合项，其模型动态环境阈值得到了提升，基于该模型的运动误差事后补偿方法的动态适应性，得到显著改善。
作为本发明的另一个实施例，利用54参数运动误差模型对数据块进行事后运动误差补偿，具体步骤如下：
(a)将t时刻重力梯度仪测量坐标系原点的比力矢量ax(t)、ay(t)、az(t)，t时刻万有引力传感器相对惯性系的角速度矢量ωx(t)、ωy(t)、ωz(t)，t时刻万有引力传感器相对惯性系的角加速度矢量ωax(t)、ωay(t)、ωaz(t)代入下式中，计算t时刻的重力梯度仪运动矢量元素m1(t) m2(t) … m54(t)：
Figure BDA0002745115030000141
Figure BDA0002745115030000142
Figure BDA0002745115030000143
Figure BDA0002745115030000144
Figure BDA0002745115030000145
Figure BDA0002745115030000146
Figure BDA0002745115030000147
Figure BDA0002745115030000148
Figure BDA0002745115030000149
Figure BDA00027451150300001410
Figure BDA00027451150300001411
Figure BDA0002745115030000151
Figure BDA0002745115030000152
Figure BDA0002745115030000153
Figure BDA0002745115030000154
Figure BDA0002745115030000155
m17(t)＝(-ωay(t)-ωx(t)ωz(t))sin(Ωt)+(ωy(t)ωz(t)-ωax(t))cos(Ωt)
m18(t)＝(ωax(t)-ωy(t)ωz(t))sin(Ωt)+(-ωay(t)-ωx(t)ωz(t))cos(Ωt)
m19(t)＝ax(t)sin(Ωt)-ay(t)cos(Ωt)；m20(t)＝ay(t)sin(Ωt)+ax(t)cos(Ωt)
m21(t)＝ax(t)az(t)sin(Ωt)-ay(t)az(t)cos(Ωt)；
m22(t)＝-ay(t)az(t)sin(Ωt)-ax(t)az(t)cos(Ωt)
m23(t)＝-2ωax(t)sin(Ωt)+2ωay(t)cos(Ωt)；m24(t)＝2ωay(t)sin(Ωt)+2ωax(t)cos(Ωt)
m25(t)＝(ωax(t)+ωy(t)ωz(t))sin(Ωt)+(ωx(t)ωz(t)-ωay(t))cos(Ωt)
m26(t)＝(ωay(t)-ωx(t)ωz(t))sin(Ωt)+(ωax(t)+ωy(t)ωz(t))cos(Ωt)
Figure BDA0002745115030000156
Figure BDA0002745115030000157
Figure BDA0002745115030000161
m30(t)＝ax(t)ωaz(t)sin(Ωt)-ay(t)ωaz(t)cos(Ωt)
Figure BDA0002745115030000162
Figure BDA0002745115030000163
Figure BDA0002745115030000164
m34(t)＝ay(t)ωaz(t)sin(Ωt)+ax(t)ωaz(t)cos(Ωt)
m35(t)＝(2az(t)ωax(t)+2az(t)ωy(t)ωz(t))sin(Ωt)+(2az(t)ωx(t)ωz(t)-2az(t)ωay(t))cos(Ωt)
m36(t)＝(2az(t)ωx(t)ωz(t)-2az(t)ωay(t))sin(Ωt)+(-2az(t)ωax(t)-2az(t)ωy(t)ωz(t))cos(Ωt)
m37(t)＝(ωay(t)ωaz(t)-ωx(t)ωz(t)ωaz(t))sin(Ωt)+(ωax(t)ωaz(t)+ωy(t)ωz(t)ωaz(t))cos(Ωt)
Figure BDA0002745115030000165
m39(t)＝(-ωax(t)ωaz(t)-ωy(t)ωz(t)ωaz(t))sin(Ωt)+(ωay(t)ωaz(t)-ωx(t)ωz(t)ωaz(t))cos(Ωt)
Figure BDA0002745115030000166
Figure BDA0002745115030000167
Figure BDA0002745115030000168
Figure BDA0002745115030000169
Figure BDA00027451150300001610
Figure BDA00027451150300001611
Figure BDA00027451150300001612
Figure BDA00027451150300001613
m52(t)＝ωaz(t)
Figure BDA00027451150300001614
m54(t)＝1
式子中，Ω表示重力梯度仪旋转圆盘角频率，也就是重力梯度仪的调制频率，sin,cos分别表示正弦函数和余弦函数；将计算的运动矢量元素m1(t) m2(t) … m54(t)，构成t时刻的重力梯度仪运动矢量：m1×54(t)＝[m1(t) m2(t) … m54(t)]。
(b)构建第二损失函数J2(C)如下：
Figure BDA0002745115030000171
式中t表示数据块的时间，Gout(t)表示t时刻重力梯度仪的输出，C54×1表示运动误差传递系数矢量，它是54行1列的矢量，m1×54(t)表示t时刻的运动矢量，它是一个1行54列矩阵，t0表示数据起始时刻，tp表示数据结束时刻；
(c)利用第二损失函数J2(C)最小，标定重力梯度仪误差传递系数矢量C：
Figure BDA0002745115030000172
式中
Figure BDA0002745115030000173
是重力梯度仪误差传递系数矢量C的估计值，
Figure BDA0002745115030000174
表示使得函数J2(C)取最小时的C值，推荐的算法有最小二乘算法、梯度下降法等。
(d)将标定的运动误差传递系数矢量
Figure BDA0002745115030000175
代入下式，剔除重力梯度仪运动误差：
Figure BDA0002745115030000176
式中，Gout(t)是t时刻重力梯度仪的原始输出，
Figure BDA0002745115030000177
是剔除了运动误差后t时刻的重力梯度仪的输出；最后将6个剔除了线运动误差、角运动误差的数据块合并，并解调，得到6小时勘探的万有引力梯度输出。
为了对本发明实施例中的方案进行验证，这里以调制类重力梯度仪(旋转加速度计重力梯度仪)为例，进行仿真分析：
利用测试质量变转速绕重力梯度仪旋转，产生调频万有引力梯度激励，同时对旋转加速度计重力梯度仪施加扫频线振动和角振动，模拟动基座重力梯度勘探的情况；利用重力梯度仪数值模型，模拟上述过程。仿真的重力梯度仪的圆盘半径R＝0.1m，圆盘旋转角频率Ω＝1.57rad/s，重力梯度仪敏感万有引力加速度的加速度计，它的模型参数、安装参数、放大电路参数列出在下表：
Figure BDA0002745115030000181
旋转加速度计重力梯度仪，运动误差补偿的方式有：在线运动误差补偿和事后运动误差补偿。对于分辨率1E，基线距离0.1m的重力梯度仪，加速度计组合输出中，10-11g量级的扰动非引力加速度，不可忽略。如果重力梯度仪的初始共模抑制比10-3，不进行在线运动误差补偿，那么事后运动误差补偿对梯度仪的运动监测精度需要优于10-8g，它要求运动监测加速度计大动态范围且高精度，难以实现，或代价很高。反之，如果只进行在线运动误差补偿，梯度仪动态环境阈值0.01g，那么要求在线误差补偿系统的共模抑制比要优于10-9，同样难以实现，或代价很高。通过将在线运动误差补偿和事后运动误差补偿相结合，可以使梯度仪对在线误差补偿的共模抑制比的需求及事后运动误差补偿的运动监测精度需求，同时降低，能显著降低仪器制造的难度、复杂性和成本。事后运动误差补偿的运动监测精度，受在线运动误差补偿的共模抑制比影响，在本仿真中，通过设置失准角、加速度计二阶非线性误差系数、线性标度系数等参数，使在线误差补偿系统的共模抑制比接近10-6。
仿真的引力源，其质量为480Kg，在重力梯度仪测量坐标系的初始位置为(0.8,0,0)，绕重力梯度仪的旋转角速度为ω(t)＝3600+560sin(0.0628t)deg/h。给重力梯度仪施加高斯分布的扫频线振动，垂直方向的线振动加速度均值为1g，标准差0.1g；水平方向的线振动均值1mg，标准差0.03g。同时，给重力梯度仪的三个方向施加相同强度的高斯角振动，角振动均值等于600deg/h，标准差等于700deg/h。重力梯度仪运动参数如下表：
梯度仪运动参数 均值 标准差
a<sub>x</sub>/a<sub>y</sub> 1mg 30mg
a<sub>z</sub> 1g 0.1g
ω<sub>x</sub>/ω<sub>y</sub>/ω<sub>z</sub> 600deg/h 700deg/h
图3是旋转加速度计重力梯度仪的一种运动监测方案。旋转加速度计重力梯度仪的运动误差源有：加速度计安装旋转圆盘的端跳、径跳、wobble等引起的线运动，加速度计安装旋转圆盘的角运动，重力梯度仪整体的线运动，重力梯度仪整体的角运动。
本仿真中，设置的在线误差补偿系统的共模抑制比接近10-6，能够有效抑制小于10ug的共模运动，因此事后运动误差补偿，只需补偿和检测大于10ug的共模运动。留有一定的余量，仿真中，重力梯度仪运动监测加速度计的精度为1ug。由于重力梯度仪加速度计安装旋转圆盘引起的线运动，能够控制在10ug以内，它能够被在线误差补偿系统补偿掉，因此不需要对该项进行监测。在重力梯度仪的壳体上，对准梯度仪测量坐标系，安装三只加速度计，监测重力梯度仪整体的线运动。
重力梯度仪加速度计安装旋转圆盘的角运动中，z向角速度约为1.57rad/s，该项量级较大，与线运动的耦合效应不可忽略，需要进行监测。因此在旋转圆盘中心安装陀螺仪，监测旋转圆盘的角运动与重力梯度仪整体的角运动；同时，在重力梯度仪壳体上安装两只陀螺仪，对准梯度仪测量坐标系的x,y轴，监测梯度仪整体的角运动的x,y分量。本仿真中，重力梯度仪的运动误差源的补偿和监测方式如下表：
Figure BDA0002745115030000201
本仿真中，用于监测重力梯度仪运动状态的加速度计的精度1ug，安装失准角0.01deg，用于监测角运动的陀螺仪精度0.1arcsec/s安装失准角0.01deg，运动监测传感器的参数如下表：
Figure BDA0002745115030000202
将梯度仪运动监测加速度计监测的梯度仪的运动，去除杆臂效应，得到重力梯度仪测量坐标系原点的比力在重力梯度仪测量坐标系中的分量。图4中数据ax表示仿真的重力梯度仪测量坐标系原点的比力在梯度仪测量坐标系x方向分量的理论值；数据ax检测，表示仿真的运动监测传感器监测的重力梯度仪测量坐标系原点的比力在梯度仪测量坐标系x方向分量；数据ax检测噪声是数据ax与数据ax检测的差值，也就是运动监测噪声，由传感器噪声和安装失准引起。梯度仪y方向的运动监测情况和x方向类似，这里没有给出，图5是仿真的z方向的运动监测情况，图中数据的含义和图4一致。图6是仿真的重力梯度仪x方向角运动监测情况，图中数据ωx表示重力梯度仪本体相对惯性系的角速度，在梯度仪测量坐标系的x方向的分量的理论值；数据ωx检测表示传感器监测的重力梯度仪本体相对惯性系的角速度，在梯度仪测量坐标系的x方向的分量；数据ωx检测噪声是数据ωx检测与数据ωx的差，它是检测噪声，由传感器噪声和安装失准引起。y方向的角运动监测和x方向类似，这里没有给出。图7是仿真的z方向角运动监测的情况，陀螺仪安装在旋转圆盘中心，监测万有引力加速度传感器相对惯性系的角速度，它包括重力梯度仪本体的角速度和圆盘旋转角速度，数据含义与图6相同。
图8是仿真的重力梯度仪的原始输出，它是重力梯度仪线运动、角运动、万有引力梯度共同激励的输出。图9是54参数重力梯度仪运动误差模型进行事后运动误差补偿的情况，图中数据正弦通道理论值，余弦通道理论值是根据测试质量参数计算的重力梯度仪正弦通道和余弦通道的理论测量值；利用重力梯度仪的运动传感器监测数据及54参数运动误差模型，将图8中重力梯度仪的原始输出，去除运动误差，并进行解调，得到54参数运动误差模型事后补偿后的正弦通道和余弦通道数据。从图中，可以看出去除运动误差后的数据和理论数据趋势一致，对补偿误差进行统计，事后运动误差补偿精度约4E。图9和图10分别是25参数和18参数重力梯度仪运动误差模型进行事后运动误差补偿后的测量值与理论值的对比，经统计分析，25参数模型的补偿精度约为56E，18参数模型的补偿精度约为60E。从仿真结果，可以看出，54参数运动误差模型的动态环境适应性最佳，25参数模型稍优于18参数模型；本发明提供的事后补偿方法能够显著提高事后运动误差补偿算法的环境阈值和环境适应性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种重力梯度仪运动误差事后补偿方法，其特征在于，包括下述步骤：
(1)获得重力梯度仪的输出数据、比力数据、角速度数据和角加速度数据，对数据进行预处理后，将数据分成N个数据块；
(2)利用重力梯度仪运动误差模型对N个数据块的运动误差进行补偿；
(3)将补偿后的N个数据块合并，并进行解调或滤波后获得万有引力梯度；所述步骤(2)中采用的重力梯度仪运动误差模型，它的54参数模型为：
Figure FDA0003156578760000011
式中
Figure FDA0003156578760000012
表示54参数运动误差模型t时刻的输出，m1×54(t)表示t时刻的运动矢量，C54×1表示运动误差传递系数矢量；m1×54(t)是包含54个元素的行向量：m1×54(t)＝[m1(t) m2(t) … m54(t)]；m1×54(t)的54个元素给出如下：
Figure FDA0003156578760000021
Figure FDA0003156578760000022
Figure FDA0003156578760000023
Figure FDA0003156578760000024
Figure FDA0003156578760000025
Figure FDA0003156578760000026
Figure FDA0003156578760000027
Figure FDA0003156578760000028
Figure FDA0003156578760000029
Figure FDA00031565787600000210
Figure FDA00031565787600000211
Figure FDA0003156578760000031
Figure FDA0003156578760000032
Figure FDA0003156578760000033
Figure FDA0003156578760000034
Figure FDA0003156578760000035
m17(t)＝(-ωay(t)-ωx(t)ωz(t))sin(Ωt)+(ωy(t)ωz(t)-ωax(t))cos(Ωt)
m18(t)＝(ωax(t)-ωy(t)ωz(t))sin(Ωt)+(-ωay(t)-ωx(t)ωz(t))cos(Ωt)
m19(t)＝ax(t)sin(Ωt)-ay(t)cos(Ωt)；m20(t)＝ay(t)sin(Ωt)+ax(t)cos(Ωt)
m21(t)＝ax(t)az(t)sin(Ωt)-ay(t)az(t)cos(Ωt)；
m22(t)＝-ay(t)az(t)sin(Ωt)-ax(t)az(t)cos(Ωt)
m23(t)＝-2ωax(t)sin(Ωt)+2ωay(t)cos(Ωt)；m24(t)＝2ωay(t)sin(Ωt)+2ωax(t)cos(Ωt)
m25(t)＝(ωax(t)+ωy(t)ωz(t))sin(Ωt)+(ωx(t)ωz(t)-ωay(t))cos(Ωt)
m26(t)＝(ωay(t)-ωx(t)ωz(t))sin(Ωt)+(ωax(t)+ωy(t)ωz(t))cos(Ωt)
Figure FDA0003156578760000036
Figure FDA0003156578760000037
Figure FDA0003156578760000041
m30(t)＝ax(t)ωaz(t)sin(Ωt)-ay(t)ωaz(t)cos(Ωt)
Figure FDA0003156578760000042
Figure FDA0003156578760000043
Figure FDA0003156578760000044
m34(t)＝ay(t)ωaz(t)sin(Ωt)+ax(t)ωaz(t)cos(Ωt)
Figure FDA0003156578760000045
Figure FDA0003156578760000046
Figure FDA0003156578760000047
Figure FDA0003156578760000048
m39(t)＝(-ωax(t)ωaz(t)-ωy(t)ωz(t)ωaz(t))sin(Ωt)+(ωay(t)ωaz(t)-ωx(t)ωz(t)ωaz(t))cos(Ωt)
Figure FDA0003156578760000049
Figure FDA00031565787600000410
Figure FDA00031565787600000411
Figure FDA00031565787600000412
Figure FDA00031565787600000413
Figure FDA00031565787600000414
Figure FDA00031565787600000415
其中，Ω表示重力梯度仪旋转圆盘角频率，也就是重力梯度仪的调制频率，sin,cos分别表示正弦函数和余弦函数；ax(t)、ay(t)、az(t)分别表示t时刻重力梯度仪测量坐标系原点的比力矢量，ωx(t)、ωy(t)、ωz(t)分别表示t时刻万有引力传感器相对惯性系的角速度矢量，ωax(t)、ωay(t)、ωaz(t)分别表示t时刻万有引力传感器相对惯性系的角加速度矢量。
2.如权利要求1所述的重力梯度仪运动误差事后补偿方法，其特征在于，采用54参数运动误差模型对N个数据块的运动误差进行补偿步骤具体包括：
(2.2.1)根据t时刻重力梯度仪测量坐标系原点的比力矢量ax(t)、ay(t)、az(t)，t时刻万有引力传感器相对惯性系的角速度矢量ωx(t)、ωy(t)、ωz(t)，t时刻万有引力传感器相对惯性系的角加速度矢量ωax(t)、ωay(t)、ωaz(t)和所述54参数运动误差模型，计算t时刻的运动矢量m1×54(t)＝[m1(t) m2(t) … m54(t)]；
(2.2.2)构建第二损失函数并利用第二损失函数最小作为约束来标定重力梯度仪误差传递系数矢量；
(2.2.3)根据标定后的运动误差传递系数矢量来剔除重力梯度仪运动误差，实现运动误差的补偿。
3.如权利要求2所述的重力梯度仪运动误差事后补偿方法，其特征在于，所述第二损失函数为
Figure FDA0003156578760000051
其中，t为数据的时间，Gout(t)表示t时刻重力梯度仪的输出，C54×1表示运动误差传递系数矢量，是54行1列的矢量，m1×54(t)表示t时刻的运动矢量，是一个1行54列矩阵，t0表示数据起始时刻，tp表示数据结束时刻。
4.如权利要求3所述的重力梯度仪运动误差事后补偿方法，其特征在于，在步骤(2.2.3)中，根据公式
Figure FDA0003156578760000052
来剔除重力梯度仪运动误差，
其中，Gout(t)为t时刻重力梯度仪的原始输出，
Figure FDA0003156578760000053
为剔除了运动误差后t时刻的重力梯度仪的输出。
CN202011163660.3A 2020-10-27 2020-10-27 一种重力梯度仪运动误差事后补偿方法 Active CN112363247B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011163660.3A CN112363247B (zh) 2020-10-27 2020-10-27 一种重力梯度仪运动误差事后补偿方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011163660.3A CN112363247B (zh) 2020-10-27 2020-10-27 一种重力梯度仪运动误差事后补偿方法

Publications (2)

Publication Number Publication Date
CN112363247A CN112363247A (zh) 2021-02-12
CN112363247B true CN112363247B (zh) 2021-09-07

Family

ID=74510754

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011163660.3A Active CN112363247B (zh) 2020-10-27 2020-10-27 一种重力梯度仪运动误差事后补偿方法

Country Status (1)

Country Link
CN (1) CN112363247B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113777671B (zh) * 2021-08-10 2022-07-12 华中科技大学 一种旋转式重力梯度仪的磁场效应评估方法、装置及系统
CN113885098B (zh) * 2021-09-07 2023-04-28 中国船舶重工集团公司第七0七研究所 一种重力敏感器低频频率响应误差在线建模及补偿方法
CN114137447B (zh) * 2022-02-07 2022-07-12 西南民族大学 磁梯度仪摆动噪声补偿方法、装置、电子设备及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011097169A3 (en) * 2010-02-03 2011-11-17 Baker Hughes Incorporated Measurement method for a component of the gravity vector
CN109709628A (zh) * 2019-02-15 2019-05-03 东南大学 一种旋转加速度计重力梯度仪标定方法
CN109766812A (zh) * 2018-12-31 2019-05-17 东南大学 一种旋转加速度计重力梯度仪运动误差事后补偿方法
CN109781075A (zh) * 2018-12-13 2019-05-21 中国航空工业集团公司上海航空测控技术研究所 一种海洋浪高测量系统及方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN210605030U (zh) * 2019-09-24 2020-05-22 天津七所高科技有限公司 一种质量体角位移实验测量装置

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011097169A3 (en) * 2010-02-03 2011-11-17 Baker Hughes Incorporated Measurement method for a component of the gravity vector
CN109781075A (zh) * 2018-12-13 2019-05-21 中国航空工业集团公司上海航空测控技术研究所 一种海洋浪高测量系统及方法
CN109766812A (zh) * 2018-12-31 2019-05-17 东南大学 一种旋转加速度计重力梯度仪运动误差事后补偿方法
CN109709628A (zh) * 2019-02-15 2019-05-03 东南大学 一种旋转加速度计重力梯度仪标定方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
旋转加速度计重力梯度仪误差补偿及信号处理;喻名彪;《中国博士学位论文全文数据库 工程科技Ⅱ辑》;20200515(第05期);30、94-95 *

Also Published As

Publication number Publication date
CN112363247A (zh) 2021-02-12

Similar Documents

Publication Publication Date Title
CN112363247B (zh) 一种重力梯度仪运动误差事后补偿方法
CN108592952B (zh) 基于杆臂补偿与正反倍速率同时标定多mimu误差的方法
CN105043412B (zh) 一种惯性测量单元误差补偿方法
CN109709628B (zh) 一种旋转加速度计重力梯度仪标定方法
CN100405014C (zh) 一种载体姿态测量方法
US11372129B2 (en) Post-compensation method for motion errors of rotating accelerometer gravity gradiometer
CN111076748A (zh) 基于mems加速度计的水平倾角仪误差补偿方法及系统
CN106840100A (zh) 一种数字式倾角传感器及测量方法
CN110542439A (zh) 基于三维气浮的惯性器件残余力矩测量装置及方法
CN107870371A (zh) 一种动基座重力梯度仪自梯度补偿方法
CN107490809A (zh) 一种重力敏感器的惯组级标定方法
KR101314151B1 (ko) 주기적 회전 진동을 이용한 6축 진동 센서의 교정 방법 및 장치
CN109855653A (zh) 一种冗余式mems-imu的降噪处理后的标定方法
KR101658473B1 (ko) Mems자이로스코프의 가속도 민감도 보정 방법
CN108931824B (zh) 一种动基座旋转加速度计重力梯度仪误差增益系数标定方法
CN2901249Y (zh) 汽车运动状态测量系统
CN116972875B (zh) 一种基于陀螺仪的物体运动轨迹监测方法
Mayer et al. Modeling and experimental analysis of low-cost MEMS gyroscopes under PCB bending stress
Shen et al. A new calibration method for MEMS inertial sensor module
CN109085654B (zh) 一种旋转加速度计重力梯度仪数字建模仿真方法
Shen et al. A new calibration method for low cost MEMS inertial sensor module
Shen et al. Evaluation of MEMS inertial sensor module for underwater vehicle navigation application
CN113945230B (zh) 一种惯性器件的高阶误差系数的辨识方法
CN113255193B (zh) 一种基于有限元分析的旋转设备虚拟测点构建方法及系统
CN115876225A (zh) 基于二自由度转台的mems imu标定方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant