CN106405670B - 一种适用于捷联式海洋重力仪的重力异常数据处理方法 - Google Patents

一种适用于捷联式海洋重力仪的重力异常数据处理方法 Download PDF

Info

Publication number
CN106405670B
CN106405670B CN201610885518.7A CN201610885518A CN106405670B CN 106405670 B CN106405670 B CN 106405670B CN 201610885518 A CN201610885518 A CN 201610885518A CN 106405670 B CN106405670 B CN 106405670B
Authority
CN
China
Prior art keywords
speed
moment
coordinate system
strapdown
data processing
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
CN201610885518.7A
Other languages
English (en)
Other versions
CN106405670A (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.)
China Aerospace Times Electronics Corp
Original Assignee
China Aerospace Times Electronics Corp
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 China Aerospace Times Electronics Corp filed Critical China Aerospace Times Electronics Corp
Priority to CN201610885518.7A priority Critical patent/CN106405670B/zh
Publication of CN106405670A publication Critical patent/CN106405670A/zh
Application granted granted Critical
Publication of CN106405670B publication Critical patent/CN106405670B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
    • G01C25/005—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass initial alignment, calibration or starting-up of inertial devices

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Manufacturing & Machinery (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Navigation (AREA)

Abstract

本发明涉及一种适用于捷联式海洋重力仪的重力异常数据处理方法，该方法采用kalman滤波器，估计惯性导航解算中的姿态误差，对姿态矩阵和垂向比力分量进行修正，得到更为精确的比力信息，再由PPP技术得提供的位置、速度和高度信息计算重力各改正项，最后经低通滤波得到沿航线的重力异常信息，本发明可用于海洋重力仪数据的处理，抗干扰能力强，具有较高的数据处理精度，可以作为高精度海洋重力仪重力异常提取的数据处理方法，特别适用于高精度捷联式动基座海洋重力仪的重力异常提取。

Description

一种适用于捷联式海洋重力仪的重力异常数据处理方法
技术领域
本发明涉及一种适用于捷联式海洋重力仪的重力异常数据处理方法，属于航空航天高精度惯性元件的测试技术领域。
背景技术
目前，重力测量在资源勘探、固体潮监测等民用和科学研究方面应用广泛。动基座重力测量系统的出现使得高效、大范围的重力测量成为可能。对比于平台式重力仪，捷联式重力仪具有体积小、功耗低、结构简单、可靠性高等优点。捷联式动基座海洋重力仪重力异常提取的数据处理方法是捷联式重力仪应用于海洋重力测量的关键技术。
要求重力测量系统的精度优于1mGal(1mGal＝10-5m/s2＝1μg)。这对于重力异常提取的数据处理方法而言，提出了很高的要求。由于重力数据的敏感性，国外成熟产品的软件已经封存为程序包，购买国外的仪器并不能知晓其详细的处理过程和方法。目前，国内对复杂海况下的捷联式重力仪重力异常的提取方法并未形成好的解决方案。
发明内容
本发明的目的在于克服现有技术的上述不足，提供一种适用于捷联式海洋重力仪的重力异常数据处理方法，该方法可用于海洋重力仪数据的处理，抗干扰能力强，具有较高的数据处理精度，可以作为高精度海洋重力仪重力异常提取的数据处理方法，特别适用于高精度捷联式动基座海洋重力仪的重力异常提取。
本发明的上述目的主要是通过如下技术方案予以实现的：
一种适用于捷联式海洋重力仪的重力异常数据处理方法，包括如下步骤：
(1)、将捷联式海洋重力仪安装于船舱内，捷联式海洋重力仪实时记录船舶的速度增量和角速度增量，并投影至载体坐标系下得到速度增量fb和角速度增量
(2)、进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵所述姿态转移矩阵为对准时间段内最后时刻k时刻的姿态转移矩阵；
(3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵获得船舶导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，其中ti为当前时刻；
(4)、将GPS测量获得的船舶运动信息进行差分处理，获得差分GPS数据；
(5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，以及步骤(4)中差分GPS数据中的当前时刻的东速VE，北速VN、经度λ和纬度L，得到相应东速的差值ve(ti)-VE、北速的差值vn(ti)-VN、经度的差值lon(ti)-λ和纬度的差值lat(ti)-L，作为卡尔曼滤波器的观测量进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计；
(6)、根据步骤(5)中得到的当前时刻的姿态误差修正东北天向比力值，得到修正后的东北天向比力值fn'；
(7)、计算重力异常粗值δg，公式如下：
其中：gb为码头处的重力基准值；
fu为fn'中的天向比力值；
为码头处的天向比力初值；
au为天向运动加速度；
δaE为厄特弗斯改正；
δaF为自由空间改正；
γ0为正常重力改正；
δgdrift为零点漂移改正。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，获得船舶导航过程中当前时刻的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)的具体方法如下：
根据导航过程中起始时刻，即第k+1时刻的速度增量fb(tk+1)、角速度增量和所述姿态转移矩阵获得导航过程中船舶在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，根据船舶在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，以及第k+2时刻的速度增量fb(tk+2)、角速度增量获得的船舶在第k+2时刻的姿态转移矩阵东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)，依次类推，获得船舶在导航过程中当前时刻的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，姿态转移矩阵通过如下方法获得：
其中：
其中：为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ωie为地球自转角速度，t0为对准时间段内的初始对准时刻，tk为对准时间段内的任意时刻；
为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：
式中：q0 q1 q2 q3为对准数据段最后时刻k的四元素；
为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：
其中：g为地球重力值，分别计算Vi(tk1)和Vi(tk2)的值，tk1和tk2分别是对准时段内的两个时刻；
Δtk＝tk-t0，
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
fb(ti)为当前时刻的速度增量。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，载体坐标系到载体惯性坐标系的姿态转移矩通过如下方法获得：
(3.1)、初始对准时刻t0，载体坐标系到载体惯性坐标系的姿态转移矩阵表示如下：
其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t0)＝[1000]；
(3.2)、根据t0时刻的四元素Q(t0)和t1时刻的角速度增量获得t1时刻的四元素其中， Φ＝|Φ|；
(3.3)、根据t1时刻的四元素Q(t1)和t2时刻的角速度增量获得t2时刻的四元素Q(t2)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q0 q1 q2 q3]，根据Q(k)计算如下：
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(4)中差分GPS数据包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(VE，VN，VU)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(5)中根据捷联惯性系统误差方程，选取的状态向量XINS为13阶，具体表示如下：
其中：δL为纬度误差；
δλ为经度误差；
δve、δvn分别为东速误差和北速误差；
φe、φn和φu分别为三个姿态误差角；
εx、εy和εz为激光陀螺的零位；
和为加速度计零位。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(6)中修正后的东北天向比力值fn'表示如下：
其中：φ×为反对称阵，
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
fb(ti)为当前时刻的速度增量；
ΔT为系统采样间隔时间。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，对所述步骤(7)中计算的重力异常粗值δg采用数字滤波器进行滤波处理，以提高数据精度。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，数字滤波采用FIR和IIR低通滤波器，截止频率小于0.01Hz；或者采用正反卡尔曼滤波器。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，若船舶在行驶过程中存在形成网格的交叉点，则进行测线网平差处理方法，消除仪器的系统误差和半系统误差，提高测量精度。
在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(4)中采用PPP技术将GPS测量获得的船舶运动信息进行差分处理。
本发明与现有技术相比的优点在于：
(1)、本发明采用kalman滤波器，估计惯性导航解算中的姿态误差，对姿态矩阵和垂向比力分量进行修正，得到更为精确的比力信息，再由PPP技术得提供的位置、速度和高度信息计算重力各改正项，最后经低通滤波得到沿航线的重力异常信息，本发明可用于海洋重力仪数据的处理，抗干扰能力强，具有较高的数据处理精度，可以作为高精度海洋重力仪重力异常提取的数据处理方法，特别适用于高精度捷联式动基座海洋重力仪的重力异常提取。
(2)、本发明在初始对准阶段，采用了惯性凝固假设，进行动基座初始对准，可以有效减小由于船体晃动带来的姿态误差，提高初始对准精度。
(3)、本发明采用PPP技术处理航次中的GPS数据，克服了针对于差分GPS技术受限于基站与移动站的距离的问题，从而获得较高精度的差分数据。
(4)、本发明采用组合导航技术获取天向比力信息，利用PPP技术处理整个航次的GPS数据，经过改正计算后，通过数字低通方式去除高频噪声得到高精度的重力异常信号。
(5)、本发明可以采用数字低通滤波器对计算得到的重力异常粗值去噪处理，截止频率小于0.01Hz，进一步提高了重力异常信号的精度；此外如果船舶在行驶过程中存在形成网格的交叉点，可以采用测线网平差方法处理航次中的系统误差和半系统误差，以提高重力异常的测量精度。
附图说明
图1为本发明适用于捷联式海洋重力仪的重力异常数据处理方法的工作流程图；
图2为本发明载体坐标系、导航坐标系和实际数学平台坐标系之间关系示意图；
图3为本发明所针对的捷联式海洋重力仪的硬件安装示意图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步详细的描述：
如图1所示为本发明适用于捷联式海洋重力仪的重力异常数据处理方法的工作流程图，本发明重力异常数据处理方法具体包括如下步骤：
(1)定义载体坐标系(b系)为右前上坐标系，记为oxbybzb，原点位于重力仪质心，xb、yb和zb分别指向重力仪右、前、上。将捷联式重力仪安装于测量船舶，重力仪Y轴指向船头即前向，此时，重力仪XYZ轴分别指向船舶的右前上方向，如图3所示为本发明所针对的捷联式重力仪的硬件安装示意图。
将捷联式海洋重力仪安装于船舱内，系统的Y轴指向船首，捷联式海洋重力仪记录船舶的速度增量和运动角速度增量，并投影至载体坐标系下得到速度增量fb和角速度增量
(2)采用了惯性凝固假设，进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵所述姿态转移矩阵为对准时间段内最后时刻k时刻的姿态转移矩阵。
对载体坐标系(b系)到实际数学平台坐标系(n系)的姿态转移矩阵进行初始化设定；海洋重力数据处理中，导航坐标系通常选为地理坐标系。如图2所示为本发明载体坐标系、导航坐标系和实际数学平台坐标系之间关系示意图；
该方法用到的坐标系定义如下：
a)经线地球坐标系e：原点位于地心，oze轴沿地球自转轴方向，oxe轴位于赤道平面内，从地心指向重力仪所在点经线，oye轴在赤道平面内，oxe、oye、oze轴构成右手坐标系。
b)经线地心惯性坐标系i：在对准起始时刻t0时刻将经线地球坐标系oxeyeze惯性凝固后形成的坐标系。
c)导航坐标系n’：原点位于捷联式重力仪中心，ox轴指向东(E)，oy轴指向北(N),oz轴指向天(U)。
d)实际数学平台坐标系n：坐标系Ox1y1z1，近似指向东北天，与理想导航坐标n’系之间存在失准角，例如，水平失准角(φe、φn)为0.005°，方位失准角φu为0.08°。
e)载体系b系：坐标系Oxbybzb，原点位于重力仪质心，xb、yb和zb指向重力仪右前上。
f)载体惯性坐标系ib0：在对准起始时刻t0时刻将载体坐标系oxbybzb经惯性凝固后的坐标系。
动基座初始对准算法中，姿态阵分散成4个矩阵求取。设对准点的纬度为L，则姿态转移矩阵通过如下方法获得：
其中：
其中：为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ωie为地球自转角速度，t0为对准时间段内的初始对准时刻，tk为对准时间段内的任意时刻；
为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：
式中：q0 q1 q2 q3为对准数据段最后时刻k的四元素。
载体坐标系到载体惯性坐标系的姿态转移矩通过如下方法获得：
(a)、初始对准时刻t0，载体坐标系到载体惯性坐标系的姿态转移矩阵表示如下：
其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t0)＝[1000]；
(b)、根据t0时刻的四元素Q(t0)和t1时刻的角速度增量获得t1时刻的四元素其中， Φ＝|Φ|；
(c)、根据t1时刻的四元素Q(t1)和t2时刻的角速度增量获得t2时刻的四元素Q(t2)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q0 q1 q2 q3]，根据Q(k)计算如下：
为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：
其中：g为地球重力值，分别计算Vi(tk1)和Vi(tk2)的值，tk1和tk2分别是对准时段内的两个时刻；
Δtk＝tk-t0，
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
fb(ti)为当前时刻的速度增量。
(3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵获得船舶导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，其中ti为当前时刻；具体方法如下：
根据导航过程中起始时刻，即第k+1时刻的速度增量fb(tk+1)、角速度增量和所述姿态转移矩阵获得导航过程中船舶在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)；
根据船舶在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，以及第k+2时刻的速度增量fb(tk+2)、角速度增量获得的船舶在第k+2时刻的姿态转移矩阵东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)；
根据船舶在第k+2时刻的姿态转移矩阵东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)，以及第k+3时刻的加速度增量fb(tk+3)、角速度增量获得的船舶在第k+3时刻的姿态转移矩阵东速ve(tk+3)、北速vn(tk+3)、纬度lat(tk+3)和经度lon(tk+3)。
依次类推，获得船舶在导航过程中当前时刻的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)。按照上述方法可以获得船舶在导航过程中各时刻的姿态转移矩阵东速ve(t)、北速vn(t)、纬度lat(t)和经度lon(t)。
(4)、根据船舶测量期间记录的GPS移动站数据和精密星历，采用waypoint软件中的PPP技术，获得差分GPS信息，包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(VE，VN，VU)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。获得差分GPS数据后，进行步骤(5)。
(5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，以及步骤(4)中差分GPS数据中的当前时刻的东速VE，北速VN、经度λ和纬度L，得到相应东速的差值ve(ti)-VE、北速的差值vn(ti)-VN、经度的差值lon(ti)-λ和纬度的差值lat(ti)-L，作为卡尔曼滤波器的观测量进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计。
根据kalman滤波器估计修正各项参数误差，特别是位置、速度和姿态误差，修正分为开环和闭环修正两种方式。完成后，进入步骤(6)。
选取的状态向量XINS为13阶，具体表示如下：
根据测量环境，忽略部分非主要误差参数，采用了如下的捷联惯性系统的误差方程
式中，δL为纬度误差；
δλ为经度误差；
δve、δvn分别为东、北速误差；
φe、φn和φu分别为三个姿态误差角，通常情况下，φ为小量；
εx、εy和εz为激光陀螺的零位；
和为加速度计零位；
Tij(i＝1,2,3；j＝1,2,3)为姿态阵的元素。
(6)、根据步骤(5)中得到的当前时刻姿态误差，修正东北天向比力值，得到修正后的东北天向比力值fn'，完成修正后进入步骤(7)。
其中：φ×为反对称阵，
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
fb(ti)为当前时刻的速度增量。
(7)、计算重力异常粗值δg，公式如下：
式中：gb为码头处的重力基准值，为设定值；
fu为fn'中的天向比力值；
为码头处的天向比力初值；
au为天向运动加速度，根据GPS提供的高度二次差分得到；
δaE为厄特弗斯改正；
δaF为自由空间改正；
γ0为正常重力改正；
δgdrift为零点漂移改正，零点漂移改正为每个航次重力测量时不同时间在同一点的观测值变化的改正。
其中计算重力异常各改正项，包括厄特弗斯改正、天向运动加速度、正常重力改正、自有空间改正、零点漂移改正，各改正项计算公式如下。厄特弗斯：
天向运动加速度：
可以根据GPS提供的伪距、载波相位和多普勒频移观测值，以及其单差、双差组建观测方程，利用最小二乘法，求出载体位置、速度和加速度。然后利用差分GPS数据，采用位置差分、速度差分或者载波相位差分等方法计算天向运动加速度。
正常重力：
γ0＝9.780327(1+0.0053024sin2L-0.0000058sin22L)
自由空间：
零点漂移改正：
重力测量的零点漂移率可按线性化近似计算，有
其中，C为此次测量的零点漂移变化率，和分别为前后校基准点处的天向比力值，和分别为前后校基准点处的重力场，t1和t0分别对应f1 u和的观测时间。
则零点漂移改正值为
δgdrift＝C·(t-t0)
式中，C为此次测量的零点漂移变化率，Δti为第i个测点的测量时刻与在基准点前校时刻的时间差。
完成计算后进入步骤(8)。
(8)采用数字滤波器对δg进行滤波，得到高精度重力异常信号，采用滤波器可采用FIR和IIR滤波器，截止频率小于0.01Hz；也可采用正反kalman滤波器。
(9)若船舶在行驶过程中存在形成网格的交叉点，则进行测线网平差处理方法，消除仪器的系统误差和半系统误差，提高测量精度。
以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。
本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (10)

1.一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：包括如下步骤：
(1)、将捷联式海洋重力仪安装于船舱内，捷联式海洋重力仪实时记录船舶的速度增量和角速度增量，并投影至载体坐标系下得到速度增量fb和角速度增量
(2)、进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵所述姿态转移矩阵为对准时间段内最后时刻k时刻的姿态转移矩阵；
(3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵获得船舶导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，其中ti为当前时刻；
(4)、将GPS测量获得的船舶运动信息进行差分处理，获得差分GPS数据；
(5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，以及步骤(4)中差分GPS数据中的当前时刻的东速VE，北速Vn、经度λ和纬度L，得到相应东速的差值ve(ti)-VE、北速的差值vn(ti)-VN、经度的差值lon(ti)-λ和纬度的差值lat(ti)-L，作为卡尔曼滤波器的观测量，进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计；
(6)、根据步骤(5)中得到的当前时刻的姿态误差修正东北天向比力值，得到修正后的东北天向比力值fn'；
(7)、计算重力异常粗值δg，公式如下：
其中：gb为码头处的重力基准值；
fu为fn'中的天向比力值；
为码头处的天向比力初值；
au为天向运动加速度；
δaE为厄特弗斯改正；
δaF为自由空间改正；
γ0为正常重力改正；
δgdrift为零点漂移改正；
所述步骤(6)中修正后的东北天向比力值fn'表示如下：
其中：φ×为反对称阵，
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
fb(ti)为当前时刻的速度增量；
ΔT为系统采样间隔时间。
2.根据权利要求1所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：获得船舶导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)的具体方法如下：
根据导航过程中起始时刻，即第k+1时刻的速度增量fb(tk+1)、角速度增量和所述姿态转移矩阵获得导航过程中船舶在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，根据船舶在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，以及第k+2时刻的速度增量fb(tk+2)、角速度增量获得的船舶在第k+2时刻的姿态转移矩阵东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)，依次类推，获得船舶在导航过程中当前时刻的姿态转移矩阵东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)。
3.根据权利要求1所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：所述姿态转移矩阵通过如下方法获得：
其中：
其中：为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ωie为地球自转角速度，t0为对准时间段内的初始对准时刻，tk为对准时间段内的任意时刻；
为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：
式中：q0、q1、q2、q3为对准数据段最后时刻k的四元素；
为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：
其中：g为地球重力值，分别计算Vi(tk1)和Vi(tk2)的值，tk1和tk2分别是对准时段内的两个时刻；
Δtk＝tk-t0，
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
fb(ti)为当前时刻的速度增量。
4.根据权利要求3所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：所述载体坐标系到载体惯性坐标系的姿态转移矩通过如下方法获得：
(3.1)、初始对准时刻t0，载体坐标系到载体惯性坐标系的姿态转移矩阵表示如下：
其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t0)＝[1 0 0 0]；
(3.2)、根据t0时刻的四元素Q(t0)和t1时刻的角速度增量获得t1时刻的四元素其中，Φ＝|Φ|；
(3.3)、根据t1时刻的四元素Q(t1)和t2时刻的角速度增量获得t2时刻的四元素Q(t2)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q0 q1 q2 q3]，根据Q(k)计算如下：
5.根据权利要求1～4之一所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：所述步骤(4)中差分GPS数据包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(VE，VN，VU)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。
6.根据权利要求1～4之一所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：所述步骤(5)中根据捷联惯性系统误差方程，选取的状态向量XINS为13阶，具体表示如下：
其中：δL为纬度误差；
δλ为经度误差；
δve、δvn分别为东速误差和北速误差；
φe、φn和φu分别为三个姿态误差角；
εx、εy和εz为激光陀螺的零位；
和为加速度计零位。
7.根据权利要求1～4之一所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：对所述步骤(7)中计算的重力异常粗值δg采用数字滤波器进行滤波处理，以提高数据精度。
8.根据权利要求7所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：所述数字滤波采用FIR和IIR低通滤波器，截止频率小于0.01Hz；或者采用正反卡尔曼滤波器。
9.根据权利要求1～4之一所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：若船舶在行驶过程中存在形成网格的交叉点，则进行测线网平差处理方法，消除仪器的系统误差和半系统误差，提高测量精度。
10.根据权利要求1～4之一所述的一种适用于捷联式海洋重力仪的重力异常数据处理方法，其特征在于：所述步骤(4)中采用PPP技术将GPS测量获得的船舶运动信息进行差分处理。
CN201610885518.7A 2016-10-10 2016-10-10 一种适用于捷联式海洋重力仪的重力异常数据处理方法 Active CN106405670B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610885518.7A CN106405670B (zh) 2016-10-10 2016-10-10 一种适用于捷联式海洋重力仪的重力异常数据处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610885518.7A CN106405670B (zh) 2016-10-10 2016-10-10 一种适用于捷联式海洋重力仪的重力异常数据处理方法

Publications (2)

Publication Number Publication Date
CN106405670A CN106405670A (zh) 2017-02-15
CN106405670B true CN106405670B (zh) 2018-06-19

Family

ID=59228939

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610885518.7A Active CN106405670B (zh) 2016-10-10 2016-10-10 一种适用于捷联式海洋重力仪的重力异常数据处理方法

Country Status (1)

Country Link
CN (1) CN106405670B (zh)

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107589464B (zh) * 2017-09-07 2018-05-15 中国石油大学(华东) 一种卫星测高重力数据与船测重力数据融合方法
CN108279440B (zh) * 2018-02-01 2019-08-23 中国自然资源航空物探遥感中心 一种航空重力测网交叉点非遍历搜索方法
CN109001829B (zh) * 2018-07-12 2019-11-05 中国人民解放军国防科技大学 一种捷联式水下动态重力测量仪
CN111123381B (zh) * 2018-11-01 2022-04-12 北京自动化控制设备研究所 一种用于平台式重力仪减小水平加速度影响的方法
CN109682397B (zh) * 2018-12-18 2021-01-29 上海航天控制技术研究所 一种不受历史数据影响快速收敛的地面静态对准方法
CN111812737B (zh) * 2020-06-17 2021-05-11 东南大学 水下导航与重力测量一体化系统
CN111735442A (zh) * 2020-06-17 2020-10-02 东南大学 一种水下重力无源导航系统
CN111722302A (zh) * 2020-06-29 2020-09-29 宁夏大学 用于auv载重力仪的垂直加速度改正方法
CN111722295B (zh) * 2020-07-04 2021-04-23 东南大学 一种水下捷联式重力测量数据处理方法
CN112346140B (zh) * 2020-10-15 2022-10-14 北京航天控制仪器研究所 一种捷联式海洋重力仪主机减摇装置
CN112487604B (zh) * 2020-10-27 2022-08-16 青岛海洋地质研究所 海洋重力仪输出数据长时间非线性漂移补偿方法
CN112415634B (zh) * 2020-10-27 2021-12-07 青岛海洋地质研究所 基于卫星重力异常信息的动态重力仪零位漂移补偿方法
CN112762927B (zh) * 2020-12-18 2021-09-10 中国人民解放军战略支援部队信息工程大学 水下动态重力数据采集半实物仿真方法及系统
CN115166856B (zh) * 2022-07-12 2024-05-28 中国自然资源航空物探遥感中心 一种无人船重磁测量方法、系统、设备及计算机可读存储介质

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7739045B2 (en) * 2006-05-31 2010-06-15 Honeywell International Inc. Rapid self-alignment of a strapdown inertial system through real-time reprocessing
US8024119B2 (en) * 2007-08-14 2011-09-20 Honeywell International Inc. Systems and methods for gyrocompass alignment using dynamically calibrated sensor data and an iterated extended kalman filter within a navigation system
CN103604442A (zh) * 2013-11-14 2014-02-26 哈尔滨工程大学 应用于捷联惯导系统在线标定的可观测性分析方法
CN103557876B (zh) * 2013-11-15 2016-01-20 山东理工大学 一种用于天线跟踪稳定平台的捷联惯导初始对准方法
CN103852799A (zh) * 2014-02-25 2014-06-11 中国人民解放军92859部队 一种基于ppp技术的船载重力测量方法

Also Published As

Publication number Publication date
CN106405670A (zh) 2017-02-15

Similar Documents

Publication Publication Date Title
CN106405670B (zh) 一种适用于捷联式海洋重力仪的重力异常数据处理方法
CN107655476B (zh) 基于多信息融合补偿的行人高精度足部导航方法
CN104181574B (zh) 一种捷联惯导系统/全球导航卫星系统组合导航滤波系统及方法
CN109813311B (zh) 一种无人机编队协同导航方法
CN107314718B (zh) 基于磁测滚转角速率信息的高速旋转弹姿态估计方法
CN105737828B (zh) 一种基于强跟踪的相关熵扩展卡尔曼滤波的组合导航方法
CN108051866B (zh) 基于捷联惯性/gps组合辅助水平角运动隔离的重力测量方法
CN102829785B (zh) 基于序列图像和基准图匹配的飞行器全参数导航方法
JP5068531B2 (ja) 測定及び記憶された重力傾度を用いて慣性航法測定値の精度を改善する方法及びシステム
CN108594283B (zh) Gnss/mems惯性组合导航系统的自由安装方法
CN106443827B (zh) 一种用于动基座重力仪的动态精度评估方法
CN103743395B (zh) 一种惯性重力匹配组合导航系统中时间延迟的补偿方法
CN106500693B (zh) 一种基于自适应扩展卡尔曼滤波的ahrs算法
EP3460399B1 (en) Methods, apparatuses, and computer programs for estimating the heading of an axis of a rigid body
CN104698485B (zh) 基于bd、gps及mems的组合导航系统及导航方法
EP1582840A1 (en) Inertial navigation system error correction
CN108225308A (zh) 一种基于四元数的扩展卡尔曼滤波算法的姿态解算方法
Rios et al. Fusion filter algorithm enhancements for a MEMS GPS/IMU
CN106932804A (zh) 天文辅助的惯性/北斗紧组合导航系统及其导航方法
CN104697520B (zh) 一体化无陀螺捷联惯导系统与gps系统组合导航方法
CN110849360B (zh) 面向多机协同编队飞行的分布式相对导航方法
CN110440830A (zh) 动基座下车载捷联惯导系统自对准方法
CN103017787A (zh) 适用于摇摆晃动基座的初始对准方法
CN103278165A (zh) 基于剩磁标定的磁测及星光备份的自主导航方法
CN105988129A (zh) 一种基于标量估计算法的ins/gnss组合导航方法

Legal Events

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