CN114089441A - 一种重力梯度仪测量系统数值仿真方法 - Google Patents

一种重力梯度仪测量系统数值仿真方法 Download PDF

Info

Publication number
CN114089441A
CN114089441A CN202111425819.9A CN202111425819A CN114089441A CN 114089441 A CN114089441 A CN 114089441A CN 202111425819 A CN202111425819 A CN 202111425819A CN 114089441 A CN114089441 A CN 114089441A
Authority
CN
China
Prior art keywords
accelerometer
coordinate system
gradiometer
time
specific force
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.)
Granted
Application number
CN202111425819.9A
Other languages
English (en)
Other versions
CN114089441B (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 CN202111425819.9A priority Critical patent/CN114089441B/zh
Publication of CN114089441A publication Critical patent/CN114089441A/zh
Application granted granted Critical
Publication of CN114089441B publication Critical patent/CN114089441B/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
    • G01V7/02—Details
    • G01V7/06—Analysis or interpretation of gravimetric records

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)
  • Length Measuring Devices With Unspecified Measuring Means (AREA)

Abstract

本发明公开一种重力梯度仪测量系统数值仿真方法，包括下述步骤：计算加速度计到检测质量的位置向量，并计算加速度计的自噪声；获得输入梯度仪的比力向量，角速度向量；计算加速度计的比力向量在梯度仪测量坐标系的坐标；计算四只加速度计在加速度计测量坐标系的输入轴的比力，输出轴的比力以及摆轴的比力；对重力梯度仪加速度计的比力进行滤波；计算加速度计的输出电压；计算重力梯度仪解调前的信号；判断更新后的时间是否小于仿真总时间，如是则返回，否则对信号Gout抗混叠滤波、降采样、带通滤波，然后进行正交幅值解调，恢复重力梯度信号。

Description

一种重力梯度仪测量系统数值仿真方法
技术领域
本发明属于重力梯度仪数值仿真技术领域，更具体地，涉及一种旋转加速度计重力梯度仪测量系统的数值仿真方法。
背景技术
目前的重力梯度仪数值仿真方法(专利公开号为CN109085654B及论文NumericalModel of Moving-Base Rotating Accelerometer Gravity Gradiometer)，使用的加速度计输入输出模型，无法完全表征加速度计的输入输出特性差异：比如带宽差异，群延时差异，通带内纹波差异等。旋转加速度计重力梯度仪的运动激励主要包括移动载体的运动激励及梯度仪内部安装加速度计的旋转圆盘的运动激励。上述重力梯度仪数值仿真方法中，没有引入梯度仪内部的运动激励。旋转加速度计重力梯度仪测量系统主要包括移动载体、隔振器、惯性稳定平台、旋转加速度计重力梯度仪，其中移动载体，用于改变梯度仪的测量位置，隔振器用于衰减载体的平动、角运动激励，惯性稳定平台使梯度仪测量坐标系跟踪地理坐标系，并衰减载体的角运动激励。
且上述重力梯度仪数值仿真方法中，仅仿真旋转加速度计重力梯度仪的工作过程，没有对运动激励的传递过程进行仿真，因此不是系统级的仿真，不能满足梯度仪测量系统各组件指标的计算的需求。
发明内容
针对现有技术的缺陷，本发明的目的在于提供一种旋转加速度计重力梯度仪测量系统的数值仿真方法，旨在解决现有的仿真方法中由于加速度计模型无法完全表征加速度计输入输出特征差异，没有引入梯度仪内部运动激励，没有模拟梯度仪运动激励传递过程，进而导致仿真精度低，不能满足梯度仪测量系统各组件指标的计算的需求的问题。
本发明提供了一种重力梯度仪测量系统数值仿真方法，包括下述步骤：
(1)设置或导入如下参数：加速度计安装位置参数，加速度计输入轴失准参数，加速度计输入输出模型参数，加速度计自噪声参数，梯度仪测量系统的平动隔离度函数、角运动隔离度函数，梯度仪测量系统的平动输入、角运动输入，梯度仪的内部运动参数，引力源参数，加速度计温度、磁场参数，仿真采样率和仿真时长参数；
(2)根据步骤(1)中加速度计的安装参数及梯度仪运动参数，计算t时刻稳定平台坐标系到梯度仪测量坐标系的变换矩阵
Figure BDA0003378375470000021
t时刻载体坐标系到梯度仪测量坐标系的变换矩阵
Figure BDA0003378375470000022
t时刻加速度计A1在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000023
t时刻梯度仪测量坐标系到加速度计A1的测量坐标系的变换矩阵
Figure BDA0003378375470000024
t时刻加速度计A2在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000025
t时刻梯度仪测量坐标系到加速度计A2的测量坐标系的变换矩阵
Figure BDA0003378375470000026
t时刻加速度计A3在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000027
t时刻梯度仪测量坐标系到加速度计A3的测量坐标系的变换矩阵
Figure BDA0003378375470000028
t时刻加速度计A4在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000029
以及t时刻梯度仪测量坐标系到加速度计A4的测量坐标系的变换矩阵
Figure BDA00033783754700000210
(3)计算t时刻加速度计A1到第p个检测质量的位置向量
Figure BDA00033783754700000211
t时刻加速度计A2到第p个检测质量的位置向量
Figure BDA00033783754700000212
t时刻加速度计A3到第p个检测质量的位置向量
Figure BDA00033783754700000213
以及t时刻加速度计A4到第p个检测质量的位置向量
Figure BDA00033783754700000214
并根据功率谱密度模型计算t时刻加速度计A1的噪声f1noise、加速度计A2的噪声f2noise、加速度计A3的噪声f3noise、加速度计A4的噪声f4noise；
(4)衰减输入梯度仪测量系统的平动，角运动激励，获得输入梯度仪的比力向量
Figure BDA0003378375470000031
角速度向量
Figure BDA0003378375470000032
(5)分别计算t时刻加速度计A1的比力向量f1在梯度仪测量坐标系的坐标，t时刻加速度计A2的比力向量f2在梯度仪测量坐标系的坐标，t时刻加速度计A3的比力向量f3在梯度仪测量坐标系的坐标以及t时刻加速度计A4的比力向量f4在梯度仪测量坐标系的坐标；
(6)分别计算t时刻四只加速度计在加速度计测量坐标系的输入轴的比力，输出轴的比力以及摆轴的比力；
(7)对重力梯度仪加速度计的比力进行滤波；
(8)分别计算t时刻加速度计A1的输出电压V1、加速度计A2的输出电压V2、加速度计A3的输出电压V3、加速度计A4的输出电压V4：
(9)计算t时刻重力梯度仪解调前的信号Gout；
(10)令t＝t+1/fs，fs是采样频率，判断更新后的时间t是否小于仿真总时间Ttotal，如是，则返回步骤(2)，否则进入步骤(11)；
(11)对信号Gout抗混叠滤波、降采样、带通滤波，然后进行正交幅值解调，恢复重力梯度信号。
更进一步地，在步骤(3)中，根据公式
Figure BDA0003378375470000033
计算加速度计A1到第p个检测质量的位置向量
Figure BDA0003378375470000034
加速度计A2到第p个检测质量的位置向量
Figure BDA0003378375470000035
加速度计A3到第p个检测质量的位置向量
Figure BDA0003378375470000036
以及加速度计A4到第p个检测质量的位置向量
Figure BDA0003378375470000037
式中
Figure BDA0003378375470000038
表示第p个检测质量在载体坐标系中的坐标，
Figure BDA0003378375470000039
表示载体坐标系原点到稳定平坐标系原点的位置向量，在载体坐标系中的坐标，
Figure BDA00033783754700000310
表示稳定平台坐标系原点到梯度仪测量坐标系原点的位置向量，在稳定平台坐标系中的坐标。
其中，功率谱密度模型为Φ(f)noise＝αf-b+ωT，式中α，b分别表示加速度计的红噪声参数，ωT分别表示加速度计的白噪声参数，f表示噪声的频率。
更进一步地，在步骤(4)中，根据公式
Figure BDA0003378375470000041
计算平动隔离后输入梯度仪的比力向量，式中，F(·)表示滤波，
Figure BDA0003378375470000042
是重力梯度仪测量系统的隔振器的比力输入；
Figure BDA0003378375470000043
是隔振器隔振后，梯度仪的比力输入。
其中，在步骤(4)中，可以根据公式
Figure BDA0003378375470000044
计算隔离角运动后输入梯度仪的角速度向量，式中，F(·)表示滤波，
Figure BDA0003378375470000045
是重力梯度仪测量系统的角速度输入，
Figure BDA0003378375470000046
是隔振器及惯性稳定平台隔离后，梯度仪的角速度输入。
更进一步地，在步骤(5)中，根据公式
Figure BDA0003378375470000047
Figure BDA0003378375470000048
Figure BDA0003378375470000049
Figure BDA00033783754700000410
分别计算t时刻加速度计A1的比力向量f1在梯度仪测量坐标系的坐标，t时刻加速度计A2的比力向量f2在梯度仪测量坐标系的坐标，t时刻加速度计A3的比力向量f3在梯度仪测量坐标系的坐标，t时刻加速度计A4的比力向量f4在梯度仪测量坐标系的坐标，式中，
Figure BDA0003378375470000051
分别是
Figure BDA0003378375470000052
的导数，
Figure BDA0003378375470000053
表示稳定平台坐标系(p坐标系)相对惯性坐标系(i坐标系)的角速度向量，在梯度仪测量坐标系中的坐标，
Figure BDA0003378375470000054
表示梯度仪测量坐标系(m坐标系)相对稳定平台坐标系(p坐标系)的角速度向量，在梯度仪测量坐标系中的坐标，
Figure BDA0003378375470000055
表示转台固连坐标系(m1坐标系)相对梯度仪测量坐标系(m坐标系)的角速度向量，在梯度仪测量坐标系中的坐标，×为向量叉乘，G为万有引力常数；∑表示求和。
更进一步地，在步骤(6)中，根据公式
Figure BDA0003378375470000056
Figure BDA0003378375470000057
分别计算t时刻四只加速度计，在加速度计测量坐标系的输入轴的比力，输出轴的比力以及摆轴的比力，式中f1i,f1o,f1p分别是加速度计A1在其测量坐标系的输入轴、输出轴、摆轴的比力；f2i,f2o,f2p分别是加速度计A2在其测量坐标系的输入轴、输出轴、摆轴的比力，f3i,f3o,f3p分别是加速度计A3在其测量坐标系的输入轴、输出轴、摆轴的比力，f4i,f4o,f4p分别是加速度计A4在其测量坐标系的输入轴、输出轴、摆轴的比力。
更进一步地，在步骤(7)中，根据
Figure BDA0003378375470000058
Figure BDA0003378375470000059
对重力梯度仪加速度计的比力进行滤波，式中，F(·)表示滤波，
Figure BDA00033783754700000510
是加速度计A1测量坐标系中的比力f1i,f1o,f1p，经过传递函数tf1i(s),tf1o(s),tf1p(s)滤波后的值；
Figure BDA0003378375470000061
是加速度计A2测量坐标系中的比力f2i,f2o,f2p，经过传递函数tf2i(s),tf2o(s),tf2p(s)滤波后的值；
Figure BDA0003378375470000062
是加速度计A3测量坐标系中的比力f3i,f3o,f3p，经过传递函数tf3i(s),tf3o(s),tf3p(s)滤波后的值；
Figure BDA0003378375470000063
是加速度计A4测量坐标系中的比力f4i,f4o,f4p，经过传递函数tf4i(s),tf4o(s),tf4p(s)滤波后的值。
更进一步地，在步骤(8)中，根据公式
Figure BDA0003378375470000064
分别计算t时刻加速度计A1的输出电压V1、加速度计A2的输出电压V2、加速度计A3的输出电压V3以及加速度计A4的输出电压V4，式中，gBt(B,T)是加速度计对磁场和温度的响应函数。
更进一步地，在步骤(9)中，根据公式Gout＝V1+V2-V3-V4计算t时刻重力梯度仪解调前的信号Gout。
本发明在加速度计输入输出模型中，引入传递函数，进而可以表征真实梯度仪加速度计的带宽差异，通带纹波差异，群延时差异；通过动力学分析，在梯度仪加速度计比力模型中，增加梯度仪内部旋转圆盘运动激励的输入；根据梯度仪测量系统的结构，引入隔振器及惯性稳定平台隔离度函数，对载体运动激励进行衰减。改进后的梯度仪测量系统数值模型，考虑了加速度计的安装位置误差，敏感轴失准，输入输出特性差异，梯度仪内部运动，梯度仪外部运动及运动传递特性，和真实的旋转加速度计重力梯度仪系统的特性一致，可解决梯度仪研制中，加速度计的安装，梯度仪测量系统各组件指标的计算，误差补偿、故障诊断等各类技术方案的仿真验证问题。
附图说明
图1是重力梯度仪测量系统示意图；
图2是本发明实施例提供的重力梯度仪测量系统输入输出原理框图；
图3是本发明实施例提供的重力梯度仪测量系统中加速度计输入输出模型结构示意图；
图4是四只加速度计的传递函数相对应的幅频特性曲线示意图；
图5是四只加速度计的传递函数相对应的相频特性曲线示意图；
图6是与隔振器的传递函数，对应的隔离度曲线；
图7是输入重力梯度仪测量系统的比力数据；
图8是经过隔振器隔振后，输入梯度仪的比力数据；
图9是角运动隔离传递函数对应的角运动隔离度曲线；
图10是输入重力梯度仪测量系统的角速度数据；
图11是经过隔振器及惯性稳定平台隔离后的梯度仪的角速度输入；
图12是加速度计安装转台中心点的位移变化，它由转台的跳动引起；
图13是加速度计安装转台的进动或章动引起的角速度；
图14是加速度计安装转台的旋转相位；
图15是旋转加速度计重力梯度仪解调前的输出，它是梯度仪内部运动，外部运动，引力梯度共同激励的输出；
图16是运动误差补偿后，梯度仪的正弦通道的输出与理论测量值的对比；
图17是运动误差补偿后，梯度仪的余弦通道的输出与理论测量值的对比。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明在加速度计输入输出模型中，引入传递函数，进而可以表征真实梯度仪加速度计的带宽差异，通带纹波差异，群延时差异；通过动力学分析，在梯度仪加速度计比力模型中，增加梯度仪内部旋转圆盘运动激励的输入；根据梯度仪测量系统的结构，引入隔振器及惯性稳定平台隔离度函数，对载体运动激励进行衰减。改进后的梯度仪测量系统数值模型，考虑了加速度计的安装位置误差，敏感轴失准，输入输出特性差异，梯度仪内部运动，梯度仪外部运动及运动传递特性，和真实的旋转加速度计重力梯度仪系统的特性一致，可解决梯度仪研制中，加速度计的安装，梯度仪测量系统各组件指标的计算，误差补偿、故障诊断等各类技术方案的仿真验证问题。
本发明提供了一种旋转加速度计重力梯度仪测量系统的数值仿真方法，包括以下步骤：
(1)设置或导入如下参数：加速度计安装位置参数，加速度计输入轴失准参数，加速度计输入输出模型参数，加速度计自噪声参数，梯度仪测量系统的平动隔离度函数、角运动隔离度函数，梯度仪测量系统的平动输入、角运动输入，梯度仪的内部运动参数，引力源参数，加速度计温度、磁场参数，仿真采样率、仿真时长参数，具体为：
设置加速度计A1的安装位置参数：径向安装距离R1，高度角β1x，初相角β1z；设置加速度计A1的输入轴失准参数：θ1x、θ1y、θ1z；设置加速度计A1的输入输出模型参数：零偏k10，标度系数k11，二阶误差系数k12，k14，k15，k16，k17，k18，电流转电压增益k1V/I，传递函数：tf1i(s),tf1o(s),tf1p(s)；设置加速度计A1的温度参数Ta1，磁场向量Ba1＝[Ba1x,Ba1y,Ba1z]；
设置加速度计A2的安装位置参数：径向安装距离R2，高度角β2x，初相角β2z；设置加速度计A2的输入轴失准参数：θ2x、θ2y、θ2z；设置加速度计A2的输入输出模型参数：零偏k20，标度系数k21，二阶误差系数k22，k24，k25，k26，k27，k28，电流转电压增益k2V/I，传递函数：tf2i(s),tf2o(s),tf2p(s)；设置加速度计A2的温度参数Ta2，磁场向量Ba2＝[Ba2x,Ba2y,Ba2z]；
设置加速度计A3的安装位置参数：径向安装距离R3，高度角β3x，初相角β3z；设置加速度计A3的输入轴失准参数：θ3x、θ3y、θ3z；设置加速度计A3的输入输出模型参数：零偏k30，标度系数k31，二阶误差系数k32，k34，k35，k36，k37，k38，电流转电压增益k3V/I，传递函数：tf3i(s),tf3o(s),tf3p(s)；设置加速度计A3的温度参数Ta3，磁场向量Ba3＝[Ba3x,Ba3y,Ba3z]；
设置加速度计A4的安装位置参数：径向安装距离R4，高度角β4x，初相角β4z；设置加速度计A4的输入轴失准参数：θ4x、θ4y、θ4z；设置加速度计A4的输入输出模型参数：零偏k40，标度系数k41，二阶误差系数k42，k44，k45，k46，k47，k48，电流转电压增益k4V/I，传递函数：tf4i(s),tf4o(s),tf4p(s)；设置加速度计A4的温度参数Ta4，磁场向量Ba4＝[Ba4x,Ba4y,Ba4z]；
设置加速度计的自噪声参数：α，b，ωT；
设置梯度仪测量系统的平动隔离度函数：gax(s)，gay(s)，gaz(s)，角运动隔离度函数：gωx(s)，gωy(s)，gωz(s)；
设置梯度仪载体的平动输入(比力向量),在载体坐标系(b坐标系)中的坐标
Figure BDA0003378375470000091
梯度仪载体的角运动输入(角速度向量)，在载体坐标系(b坐标系)中的坐标
Figure BDA0003378375470000092
设置梯度仪的内部运动参数：梯度仪内部加速度计安装圆盘进动及章动引起的角速度向量，在梯度仪测量坐标系(m坐标系)中的坐标
Figure BDA0003378375470000093
梯度仪内部加速度计安装圆盘角速度向量,在转盘固连坐标系(m1坐标系)中的坐标
Figure BDA0003378375470000094
转盘中心点，在稳定平台坐标系(p坐标系)中的位置向量
Figure BDA0003378375470000095
稳定平台坐标系原点，在载体坐标系(b坐标系)中的位置向量
Figure BDA0003378375470000096
设置N个检测质量的参数：设第p个检测质量的质量为mp，第p个检测质量在载体坐标系(b坐标系)中的位置向量
Figure BDA0003378375470000101
设置采样频率fs，仿真总时间Ttotal。
(2)根据步骤(1)中加速度计的安装参数及梯度仪运动参数，计算t时刻稳定平台坐标系(p坐标系)到梯度仪测量坐标系(m坐标系)的变换矩阵
Figure BDA0003378375470000102
计算t时刻载体坐标系(b坐标系)到梯度仪测量坐标系(m坐标系)的变换矩阵
Figure BDA0003378375470000103
计算t时刻加速度计A1在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000104
计算t时刻梯度仪测量坐标系到加速度计A1的测量坐标系的变换矩阵
Figure BDA0003378375470000105
计算t时刻加速度计A2在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000106
计算t时刻梯度仪测量坐标系到加速度计A2的测量坐标系的变换矩阵
Figure BDA0003378375470000107
计算t时刻加速度计A3在重力梯度仪测量坐标系中的位置向量
Figure BDA0003378375470000108
计算t时刻梯度仪测量坐标系到加速度计A3的测量坐标系的变换矩阵
Figure BDA0003378375470000109
计算t时刻加速度计A4在重力梯度仪测量坐标系中的位置向量
Figure BDA00033783754700001010
计算t时刻梯度仪测量坐标系到加速度计A4的测量坐标系的变换矩阵
Figure BDA00033783754700001011
(3)根据下式，计算t时刻加速度计A1到第p个检测质量的位置向量
Figure BDA00033783754700001012
计算t时刻加速度计A2到第p个检测质量的位置向量
Figure BDA00033783754700001013
计算t时刻加速度计A3到第p个检测质量的位置向量
Figure BDA00033783754700001014
计算t时刻加速度计A4到第p个检测质量的位置向量
Figure BDA00033783754700001015
Figure BDA00033783754700001016
根据以下功率谱密度模型计算t时刻加速度计A1的噪声f1noise、加速度计A2的噪声f2noise、加速度计A3的噪声f3noise、加速度计A4的噪声f4noise：Φ(f)noise＝αf-b+ωT；
(4)衰减输入梯度仪测量系统的平动，角运动激励，获得输入梯度仪测量坐标系的平动、角运动激励：根据下式计算平动隔离后，输入梯度仪的比力向量
Figure BDA0003378375470000111
式中，F(·)表示滤波，
Figure BDA0003378375470000112
是重力梯度仪测量系统的隔振器的比力输入；
Figure BDA0003378375470000113
是隔振器隔振后，梯度仪的比力输入；根据下式计算隔振器及惯性稳定平台，隔离角运动后，输入梯度仪的角速度向量
Figure BDA0003378375470000114
Figure BDA0003378375470000115
式中，F(·)表示滤波，
Figure BDA0003378375470000116
是重力梯度仪测量系统的角速度输入，
Figure BDA0003378375470000117
是隔振器及惯性稳定平台隔离后，梯度仪的角速度输入；
(5)根据下式分别计算t时刻加速度计A1的比力向量f1在梯度仪测量坐标系的坐标，t时刻加速度计A2的比力向量f2在梯度仪测量坐标系的坐标，t时刻加速度计A3的比力向量f3在梯度仪测量坐标系的坐标，t时刻加速度计A4的比力向量f4在梯度仪测量坐标系的坐标：
Figure BDA0003378375470000121
Figure BDA0003378375470000122
Figure BDA0003378375470000123
Figure BDA0003378375470000124
式中，
Figure BDA0003378375470000125
分别是
Figure BDA0003378375470000126
的导数，×为向量叉乘，G为万有引力常数；∑表示求和；
(6)根据下式分别计算t时刻四只加速度计，在加速度计测量坐标系的输入轴的比力，输出轴的比力，摆轴的比力：
Figure BDA0003378375470000127
式中f1i,f1o,f1p分别是加速度计A1在其测量坐标系的输入轴、输出轴、摆轴的比力；f2i,f2o,f2p分别是加速度计A2在其测量坐标系的输入轴、输出轴、摆轴的比力，f3i,f3o,f3p分别是加速度计A3在其测量坐标系的输入轴、输出轴、摆轴的比力，f4i,f4o,f4p分别是加速度计A4在其测量坐标系的输入轴、输出轴、摆轴的比力；
(7)根据下式对重力梯度仪加速度计的比力进行滤波：
Figure BDA0003378375470000128
式中，F(·)表示滤波，
Figure BDA0003378375470000129
是加速度计A1测量坐标系中的比力f1i,f1o,f1p，经过传递函数tf1i(s),tf1o(s),tf1p(s)滤波后的值；
Figure BDA00033783754700001210
是加速度计A2测量坐标系中的比力f2i,f2o,f2p，经过传递函数tf2i(s),tf2o(s),tf2p(s)滤波后的值；
Figure BDA0003378375470000131
是加速度计A3测量坐标系中的比力f3i,f3o,f3p，经过传递函数tf3i(s),tf3o(s),tf3p(s)滤波后的值；
Figure BDA0003378375470000132
是加速度计A4测量坐标系中的比力f4i,f4o,f4p，经过传递函数tf4i(s),tf4o(s),tf4p(s)滤波后的值；
(8)根据下式分别计算t时刻加速度计A1的输出电压V1、加速度计A2的输出电压V2、加速度计A3的输出电压V3、加速度计A4的输出电压V4：
Figure BDA0003378375470000133
式中，gBt(B,T)是加速度计对磁场和温度的响应函数；
(9)根据下式计算t时刻重力梯度仪解调前的信号Gout：Gout＝V1+V2-V3-V4；
(10)令t＝t+1/fs，然后判断更新后的时间t是否小于仿真总时间Ttotal，如是，则返回步骤2)，否则进入步骤11)；
(11)对信号Gout抗混叠滤波、降采样、带通滤波，然后进行正交幅值解调，恢复重力梯度信号。
在本发明实施例中，惯性坐标系，载体坐标系，平台坐标系，梯度仪测量坐标系，圆盘固连坐标系，加速度计测量坐标系的定义如下：惯性坐标系，定义为地球中心惯性坐标系；载体坐标系的原点位于载体中心，x轴沿载体横轴向右，y轴沿载体纵轴向前，z轴垂直xy平面向上；稳定平台坐标系的原点位于稳定平台内框，梯度仪的底座中心，x轴为内框轴方向，y轴为中框轴方向，z轴垂直xy平面，它跟踪当地地理坐标系；梯度仪测量坐标系的原点位于圆盘中心，它的x轴指向加速度计A1的初始安装位置，y轴指向加速度计A3的初始安装位置，z轴垂直xy平面；圆盘固连坐标系，它的原点位于圆盘中心，该坐标系和圆盘固连，跟随圆盘旋转，初始状态下和梯度仪测量坐标系重合；加速度计测量坐标系，它的原点位于检测质量的中心，x轴和输入轴方向一致，y轴和输出轴方向一致，z轴和摆轴方向一致。
为了对本发明实施例中的方案进行验证，可以进行如下仿真分析：图1是旋转加速度计重力梯度仪测量系统示意图，它主要包括，载体，引力源，隔振器，稳定平台，梯度仪；在数值模型中，定义了六个坐标系：惯性坐标系，载体坐标系，平台坐标系，梯度仪测量坐标系，圆盘固连坐标系，加速度计测量坐标系。惯性坐标系，定义为地球中心惯性坐标系；载体坐标系的原点位于载体中心，x轴沿载体横轴向右，y轴沿载体纵轴向前，z轴垂直xy平面向上；稳定平台坐标系的原点位于稳定平台内框，梯度仪的底座中心，x轴为内框轴方向，y轴为中框轴方向，z轴垂直xy平面，它跟踪当地地理坐标系；梯度仪测量坐标系的原点位于圆盘中心，它的x轴指向加速度计A1的初始安装位置，y轴指向加速度计A3的初始安装位置，z轴垂直xy平面；圆盘固连坐标系，它的原点位于圆盘中心，该坐标系和圆盘固连，跟随圆盘旋转，初始状态下和梯度仪测量坐标系重合；加速度计测量坐标系，它的原点位于检测质量的中心，x轴和输入轴方向一致，y轴和输出轴方向一致，z轴和摆轴方向一致。
图2是梯度仪测量系统的输入输出，它描述了梯度仪的输入输出过程；仿真的重力梯度仪的圆盘半径R＝0.1m，重力梯度仪加速度计安装参数、加速度计的线性标度系数及二阶非线性误差系数给出在下表：
Figure BDA0003378375470000151
图3是梯度仪加速度计的输入输出模型，四只加速度计的传递函数给出如下：
Figure BDA0003378375470000161
Figure BDA0003378375470000162
在实际中，由于制造加速度计的材料等存在差异，梯度仪的加速度计的输入输出特性存在差异，表现为不同频率上的增益及延时存在差异。图4是四只加速度计的传递函数相对应的幅频特性曲线，它的横轴表示频率，纵轴表示增益，从图中可以看出四只加速度计，在不同频率上，它的增益存在差异。图5是上述四只加速度计的传递函数相对应的相频特性曲线，它的横轴表示频率，纵轴表示相位，从图中可以看出四只加速度计，在不同频率上，它的相位存在差异。综合图4和图5，四只加速度计，对某频率的激励信号，它的增益、延时均存在差异，这与梯度仪的实际情况相符合，这些差异，是梯度仪运动误差不能完全抑制的重要因素。
仿真的重力梯度仪测量系统的隔振器，它的传递函数如下：
Figure BDA0003378375470000163
图6是与隔振器的传递函数，对应的隔离度曲线(幅频特性曲线)；它的横轴表示频率，纵轴表示增益，从图中可以看出，它对大于0.7Hz的频率成分有隔离作用。图7是输入重力梯度仪测量系统的比力数据，它的横轴表示样本数，纵轴表示比力，单位为g(1g＝9.8m/s2)。图8是经过隔振器隔振后，输入梯度仪的比力数据，也就是将图7中的数据，通过隔振器的传递函数进行处理后的输出，由于隔振器对大于0.7Hz的信号有隔离作用(衰减)，和图7进行比较，图8中的信号幅度窄了很多。惯性稳定平台及隔振器对角运动的隔离度函数如下：
Figure BDA0003378375470000171
图9是上述传递函数对应的角运动隔离度曲线(幅频特性曲线)，它的横轴表示频率，纵轴表示增益。图10是输入重力梯度仪测量系统的角速度数据，它的横轴表示样本数，纵轴表示角速度大小，单位是rad/s。图11是经过隔振器及惯性稳定平台隔离后梯度仪的角速度输入，和图10相比，可以看出经过隔离后，角速度衰减了约3个量级。测试质量680Kg，在重力梯度仪测量坐标系的初始位置为(0.9,0,0)，绕重力梯度仪的旋转角速度为ω(t)＝2600+60sin(0.0628t)deg/h，产生引力梯度激励。图12是加速度计安装转台中心点位移的变化量，它由转台的跳动引起，它的横轴表示样本数，纵轴表示跳动距离，单位为m，本仿真中，施加的跳动为高斯白噪声，跳动标准差为9nm。图13是加速度计安装转台的进动或章动引起的角速度，它的横轴表示时间，纵轴表示角速度，单位为rad/s，实际中由重力矩、摩擦力矩、机械限位等因素共同引起。图14是加速度计安装转台的旋转相位，它的横轴表示时间，纵轴表示旋转的角度，圆盘旋转一周，圆盘旋转相位从0增大到2π，呈周期变化。图15是旋转加速度计重力梯度仪解调前的输出，它是梯度仪内部运动，外部运动，引力梯度共同激励的输出。记录梯度仪圆盘中心点的比力及角速度向量，根据梯度仪解析模型，可以扣除梯度仪输出中的运动误差，图16和图17分别是运动误差补偿后，梯度仪的正弦通道、余弦通道的输出与理论测量值的对比，横轴表示时间，纵轴表示梯度，单位为E，从图中可以看出，扣除运动误差后，梯度仪的测量和理论值相一致。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种重力梯度仪测量系统数值仿真方法，其特征在于，包括下述步骤：
(1)设置或导入如下参数：加速度计安装位置参数，加速度计输入轴失准参数，加速度计输入输出模型参数，加速度计自噪声参数，梯度仪测量系统的平动隔离度函数、角运动隔离度函数，梯度仪测量系统的平动输入、角运动输入，梯度仪的内部运动参数，引力源参数，加速度计温度、磁场参数，仿真采样率和仿真时长参数；
(2)根据步骤(1)中加速度计的安装参数及梯度仪运动参数，计算t时刻稳定平台坐标系到梯度仪测量坐标系的变换矩阵
Figure FDA0003378375460000011
t时刻载体坐标系到梯度仪测量坐标系的变换矩阵
Figure FDA0003378375460000012
t时刻加速度计A1在重力梯度仪测量坐标系中的位置向量
Figure FDA0003378375460000013
t时刻梯度仪测量坐标系到加速度计A1的测量坐标系的变换矩阵
Figure FDA0003378375460000014
t时刻加速度计A2在重力梯度仪测量坐标系中的位置向量
Figure FDA0003378375460000015
t时刻梯度仪测量坐标系到加速度计A2的测量坐标系的变换矩阵
Figure FDA0003378375460000016
t时刻加速度计A3在重力梯度仪测量坐标系中的位置向量
Figure FDA0003378375460000017
t时刻梯度仪测量坐标系到加速度计A3的测量坐标系的变换矩阵
Figure FDA0003378375460000018
t时刻加速度计A4在重力梯度仪测量坐标系中的位置向量
Figure FDA0003378375460000019
以及t时刻梯度仪测量坐标系到加速度计A4的测量坐标系的变换矩阵
Figure FDA00033783754600000110
(3)计算t时刻加速度计A1到第p个检测质量的位置向量
Figure FDA00033783754600000111
t时刻加速度计A2到第p个检测质量的位置向量
Figure FDA00033783754600000112
t时刻加速度计A3到第p个检测质量的位置向量
Figure FDA00033783754600000113
以及t时刻加速度计A4到第p个检测质量的位置向量
Figure FDA00033783754600000114
并根据功率谱密度模型计算t时刻加速度计A1的噪声f1noise、加速度计A2的噪声f2noise、加速度计A3的噪声f3noise、加速度计A4的噪声f4noise；
(4)衰减输入梯度仪测量系统的平动，角运动激励，获得输入梯度仪的比力向量
Figure FDA0003378375460000021
角速度向量
Figure FDA0003378375460000022
(5)分别计算t时刻加速度计A1的比力向量f1在梯度仪测量坐标系的坐标，t时刻加速度计A2的比力向量f2在梯度仪测量坐标系的坐标，t时刻加速度计A3的比力向量f3在梯度仪测量坐标系的坐标以及t时刻加速度计A4的比力向量f4在梯度仪测量坐标系的坐标；
(6)分别计算t时刻四只加速度计在加速度计测量坐标系的输入轴的比力，输出轴的比力以及摆轴的比力；
(7)对重力梯度仪加速度计的比力进行滤波；
(8)分别计算t时刻加速度计A1的输出电压V1、加速度计A2的输出电压V2、加速度计A3的输出电压V3、加速度计A4的输出电压V4：
(9)计算t时刻重力梯度仪解调前的信号Gout；
(10)令t＝t+1/fs，fs是采样频率，判断更新后的时间t是否小于仿真总时间Ttotal，如是，则返回步骤(2)，否则进入步骤(11)；
(11)对信号Gout抗混叠滤波、降采样、带通滤波，然后进行正交幅值解调，恢复重力梯度信号。
2.如权利要求1所述的数值仿真方法，其特征在于，在步骤(3)中，根据公式
Figure FDA0003378375460000023
计算加速度计A1到第p个检测质量的位置向量
Figure FDA0003378375460000024
加速度计A2到第p个检测质量的位置向量
Figure FDA0003378375460000025
加速度计A3到第p个检测质量的位置向量
Figure FDA0003378375460000026
以及加速度计A4到第p个检测质量的位置向量
Figure FDA0003378375460000027
式中
Figure FDA0003378375460000028
表示第p个检测质量在载体坐标系中的坐标，
Figure FDA0003378375460000029
表示载体坐标系原点到稳定平坐标系原点的位置向量，在载体坐标系中的坐标，
Figure FDA0003378375460000031
表示稳定平台坐标系原点到梯度仪测量坐标系原点的位置向量，在稳定平台坐标系中的坐标。
3.如权利要求1或2所述的数值仿真方法，其特征在于，在步骤(3)中，所述功率谱密度模型为Φ(f)noise＝αf-b+ωT，式中α，b分别表示加速度计的红噪声参数，ωT分别表示加速度计的白噪声参数，f表示噪声的频率。
4.如权利要求1-3任一项所述的数值仿真方法，其特征在于，在步骤(4)中，根据公式
Figure FDA0003378375460000032
计算平动隔离后输入梯度仪的比力向量，式中，F(·)表示滤波，
Figure FDA0003378375460000033
是重力梯度仪测量系统的隔振器的比力输入；
Figure FDA0003378375460000034
是隔振器隔振后，梯度仪的比力输入。
5.如权利要求1-4任一项所述的数值仿真方法，其特征在于，在步骤(4)中，根据公式
Figure FDA0003378375460000035
计算隔离角运动后输入梯度仪的角速度向量，式中，F(·)表示滤波，
Figure FDA0003378375460000036
是重力梯度仪测量系统的角速度输入，
Figure FDA0003378375460000037
是隔振器及惯性稳定平台隔离后，梯度仪的角速度输入。
6.如权利要求1-5任一项所述的数值仿真方法，其特征在于，在步骤(5)中，根据公式
Figure FDA0003378375460000041
Figure FDA0003378375460000042
Figure FDA0003378375460000043
Figure FDA0003378375460000044
分别计算t时刻加速度计A1的比力向量f1在梯度仪测量坐标系的坐标，t时刻加速度计A2的比力向量f2在梯度仪测量坐标系的坐标，t时刻加速度计A3的比力向量f3在梯度仪测量坐标系的坐标，t时刻加速度计A4的比力向量f4在梯度仪测量坐标系的坐标，式中，
Figure FDA0003378375460000045
分别是
Figure FDA0003378375460000046
的导数，
Figure FDA0003378375460000047
表示稳定平台坐标系相对惯性坐标系的角速度向量，在梯度仪测量坐标系中的坐标，
Figure FDA0003378375460000048
表示梯度仪测量坐标系相对稳定平台坐标系的角速度向量，在梯度仪测量坐标系中的坐标，
Figure FDA0003378375460000049
表示转台固连坐标系相对梯度仪测量坐标系的角速度向量，在梯度仪测量坐标系中的坐标，×为向量叉乘，G为万有引力常数；∑表示求和。
7.如权利要求1-6任一项所述的数值仿真方法，其特征在于，在步骤(6)中，根据公式
Figure FDA00033783754600000410
分别计算t时刻四只加速度计，在加速度计测量坐标系的输入轴的比力，输出轴的比力以及摆轴的比力，式中f1i,f1o,f1p分别是加速度计A1在其测量坐标系的输入轴、输出轴、摆轴的比力；f2i,f2o,f2p分别是加速度计A2在其测量坐标系的输入轴、输出轴、摆轴的比力，f3i,f3o,f3p分别是加速度计A3在其测量坐标系的输入轴、输出轴、摆轴的比力，f4i,f4o,f4p分别是加速度计A4在其测量坐标系的输入轴、输出轴、摆轴的比力。
8.如权利要求1-7任一项所述的数值仿真方法，其特征在于，在步骤(7)中，根据
Figure FDA0003378375460000051
对重力梯度仪加速度计的比力进行滤波，式中，F(·)表示滤波，
Figure FDA0003378375460000052
是加速度计A1测量坐标系中的比力f1i,f1o,f1p，经过传递函数tf1i(s),tf1o(s),tf1p(s)滤波后的值；
Figure FDA0003378375460000053
是加速度计A2测量坐标系中的比力f2i,f2o,f2p，经过传递函数tf2i(s),tf2o(s),tf2p(s)滤波后的值；
Figure FDA0003378375460000054
是加速度计A3测量坐标系中的比力f3i,f3o,f3p，经过传递函数tf3i(s),tf3o(s),tf3p(s)滤波后的值；
Figure FDA0003378375460000055
是加速度计A4测量坐标系中的比力f4i,f4o,f4p，经过传递函数tf4i(s),tf4o(s),tf4p(s)滤波后的值。
9.如权利要求1-8任一项所述的数值仿真方法，其特征在于，在步骤(8)中，根据公式
Figure FDA0003378375460000056
分别计算t时刻加速度计A1的输出电压V1、加速度计A2的输出电压V2、加速度计A3的输出电压V3以及加速度计A4的输出电压V4，式中，gBt(B,T)是加速度计对磁场和温度的响应函数。
10.如权利要求1-9任一项所述的数值仿真方法，其特征在于，在步骤(9)中，根据公式Gout＝V1+V2-V3-V4计算t时刻重力梯度仪解调前的信号Gout。
CN202111425819.9A 2021-11-26 2021-11-26 一种重力梯度仪测量系统数值仿真方法 Active CN114089441B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111425819.9A CN114089441B (zh) 2021-11-26 2021-11-26 一种重力梯度仪测量系统数值仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111425819.9A CN114089441B (zh) 2021-11-26 2021-11-26 一种重力梯度仪测量系统数值仿真方法

Publications (2)

Publication Number Publication Date
CN114089441A true CN114089441A (zh) 2022-02-25
CN114089441B CN114089441B (zh) 2022-08-30

Family

ID=80305122

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111425819.9A Active CN114089441B (zh) 2021-11-26 2021-11-26 一种重力梯度仪测量系统数值仿真方法

Country Status (1)

Country Link
CN (1) CN114089441B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002044757A2 (en) * 2000-11-28 2002-06-06 Business Arts Inc. Gravity gradiometry
US20130055808A1 (en) * 2010-03-29 2013-03-07 Frank Joachim Van Kann Gravity gradiometer with correction of external disturbance
CN104459826A (zh) * 2014-11-03 2015-03-25 东南大学 旋转加速度计重力梯度仪重力梯度信号仿真方法
WO2018161474A1 (zh) * 2017-03-09 2018-09-13 中国科学院电工研究所 重力梯度测量方法及装置
CN109084756A (zh) * 2018-06-20 2018-12-25 东南大学 一种重力视运动参数辨识与加速度计零偏分离方法
CN109085654A (zh) * 2018-06-11 2018-12-25 东南大学 一种旋转加速度计重力梯度仪数字建模仿真方法
CN109709628A (zh) * 2019-02-15 2019-05-03 东南大学 一种旋转加速度计重力梯度仪标定方法
CN112325886A (zh) * 2020-11-02 2021-02-05 北京航空航天大学 一种基于重力梯度仪和陀螺仪组合的航天器自主定姿系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002044757A2 (en) * 2000-11-28 2002-06-06 Business Arts Inc. Gravity gradiometry
US20130055808A1 (en) * 2010-03-29 2013-03-07 Frank Joachim Van Kann Gravity gradiometer with correction of external disturbance
CN104459826A (zh) * 2014-11-03 2015-03-25 东南大学 旋转加速度计重力梯度仪重力梯度信号仿真方法
WO2018161474A1 (zh) * 2017-03-09 2018-09-13 中国科学院电工研究所 重力梯度测量方法及装置
CN109085654A (zh) * 2018-06-11 2018-12-25 东南大学 一种旋转加速度计重力梯度仪数字建模仿真方法
CN109084756A (zh) * 2018-06-20 2018-12-25 东南大学 一种重力视运动参数辨识与加速度计零偏分离方法
CN109709628A (zh) * 2019-02-15 2019-05-03 东南大学 一种旋转加速度计重力梯度仪标定方法
CN112325886A (zh) * 2020-11-02 2021-02-05 北京航空航天大学 一种基于重力梯度仪和陀螺仪组合的航天器自主定姿系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
MINGBIAO YU 等: "Post Error Compensation of Moving-base Rotating", 《IEEET RANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》 *
刘杰 等: "重力梯度仪仿真系统设计", 《吉林大学学报》 *
蔡体菁 等: "旋转加速度计重力梯度仪重力梯度信号仿真", 《物探与化探》 *

Also Published As

Publication number Publication date
CN114089441B (zh) 2022-08-30

Similar Documents

Publication Publication Date Title
CN107015287B (zh) 一种重力梯度测量装置及测量方法
CN105043412B (zh) 一种惯性测量单元误差补偿方法
US5357802A (en) Rotating accelerometer gradiometer
US11372129B2 (en) Post-compensation method for motion errors of rotating accelerometer gravity gradiometer
CN106052595B (zh) 基于激光陀螺捷联惯导的三轴转台轴线垂直度检测方法
CN109782023B (zh) 一种通过旋转调制法测量加速度计高阶项数系数的方法
CN109709628B (zh) 一种旋转加速度计重力梯度仪标定方法
US8359920B2 (en) Gravity sensing instrument
CN105137804A (zh) 一种针对飞行姿态扰动的实验室模拟方法
CN107064559A (zh) 一种基于角摇摆运动的sins加速度计频率特性测试方法
CN112363247A (zh) 一种重力梯度仪运动误差事后补偿方法
CN113945230B (zh) 一种惯性器件的高阶误差系数的辨识方法
CN107289921A (zh) 一种基于椭圆拟合的对抛式冷原子干涉陀螺仪的转动角速度测量方法
CN108931824B (zh) 一种动基座旋转加速度计重力梯度仪误差增益系数标定方法
CN111650664A (zh) 一种航空重力梯度仪实时重力梯度解调方法及装置
CN105478245A (zh) 基于主轴振动检测的双自由度精密离心机副轴动不平衡量辨识方法
CN114089441B (zh) 一种重力梯度仪测量系统数值仿真方法
CN113701747B (zh) 一种基于离心机激励的惯性测量系统姿态角误差分离方法
Yingbo et al. Calibration method of quartz accelerometer on dynamic centrifuge
CN112683303B (zh) 一种惯性测量单元陀螺位置补偿方法
CN113916219A (zh) 一种基于离心机激励的惯性测量系统误差分离方法
CN109085654B (zh) 一种旋转加速度计重力梯度仪数字建模仿真方法
Flack et al. Comparison of the unbalance responses of Jeffcott rotors with shaft bow and shaft runout
CN105716626A (zh) 一种悬浮类陀螺仪的定子旋转调制误差补偿方法
Yu et al. Angular rate measurement method of magnetically suspended control and sensing gyroscope using component-level rotation modulation

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