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

本发明涉及一种适用于捷联式海洋重力仪的重力异常数据处理方法，属于航空航天高精度惯性元件的测试技术领域。

背景技术：

目前，重力测量在资源勘探、固体潮监测等民用和科学研究方面应用广泛。动基座重力测量系统的出现使得高效、大范围的重力测量成为可能。对比于平台式重力仪，捷联式重力仪具有体积小、功耗低、结构简单、可靠性高等优点。捷联式动基座海洋重力仪重力异常提取的数据处理方法是捷联式重力仪应用于海洋重力测量的关键技术。

要求重力测量系统的精度优于1mGal(1mGal＝10^-5m/s²＝1μg)。这对于重力异常提取的数据处理方法而言，提出了很高的要求。由于重力数据的敏感性，国外成熟产品的软件已经封存为程序包，购买国外的仪器并不能知晓其详细的处理过程和方法。目前，国内对复杂海况下的捷联式重力仪重力异常的提取方法并未形成好的解决方案。

技术实现要素：

本发明的目的在于克服现有技术的上述不足，提供一种适用于捷联式海洋重力仪的重力异常数据处理方法，该方法可用于海洋重力仪数据的处理，抗干扰能力强，具有较高的数据处理精度，可以作为高精度海洋重力仪重力异常提取的数据处理方法，特别适用于高精度捷联式动基座海洋重力仪的重力异常提取。

本发明的上述目的主要是通过如下技术方案予以实现的：

一种适用于捷联式海洋重力仪的重力异常数据处理方法，包括如下步骤：

(1)、将捷联式海洋重力仪安装于船舱内，捷联式海洋重力仪实时记录船舶的速度增量和角速度增量，并投影至载体坐标系下得到速度增量f^b和角速度增量

(2)、进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵所述姿态转移矩阵为对准时间段内最后时刻k时刻的姿态转移矩阵；

(3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵获得船舶导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)，其中t_i为当前时刻；

(4)、将GPS测量获得的船舶运动信息进行差分处理，获得差分GPS数据；

(5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)，以及步骤(4)中差分GPS数据中的当前时刻的东速V_E，北速V_N、经度λ和纬度L，得到相应东速的差值v_e(t_i)-V_E、北速的差值v_n(t_i)-V_N、经度的差值lon(t_i)-λ和纬度的差值lat(t_i)-L，作为卡尔曼滤波器的观测量进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计；

(6)、根据步骤(5)中得到的当前时刻的姿态误差修正东北天向比力值，得到修正后的东北天向比力值f^n'；

(7)、计算重力异常粗值δg，公式如下：

其中：g_b为码头处的重力基准值；

f^u为f^n'中的天向比力值；

为码头处的天向比力初值；

a^u为天向运动加速度；

δa_E为厄特弗斯改正；

δa_F为自由空间改正；

γ₀为正常重力改正；

δg_drift为零点漂移改正。

在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，获得船舶导航过程中当前时刻的姿态转移矩阵东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)的具体方法如下：

根据导航过程中起始时刻，即第k+1时刻的速度增量f^b(t_k+1)、角速度增量和所述姿态转移矩阵获得导航过程中船舶在第k+1时刻的姿态转移矩阵东速v_e(t_k+1)、北速v_n(t_k+1)、纬度lat(t_k+1)和经度lon(t_k+1)，根据船舶在第k+1时刻的姿态转移矩阵东速v_e(t_k+1)、北速v_n(t_k+1)、纬度lat(t_k+1)和经度lon(t_k+1)，以及第k+2时刻的速度增量f^b(t_k+2)、角速度增量获得的船舶在第k+2时刻的姿态转移矩阵东速v_e(t_k+2)、北速v_n(t_k+2)、纬度lat(t_k+2)和经度lon(t_k+2)，依次类推，获得船舶在导航过程中当前时刻的姿态转移矩阵东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)。

在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，姿态转移矩阵通过如下方法获得：

其中：

其中：为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ω_ie为地球自转角速度，t₀为对准时间段内的初始对准时刻，t_k为对准时间段内的任意时刻；

为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：

式中：q₀ q₁ q₂ q₃为对准数据段最后时刻k的四元素；

为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：

其中：g为地球重力值，分别计算Vⁱ(t_k1)和Vⁱ(t_k2)的值，t_k1和t_k2分别是对准时段内的两个时刻；

Δt_k＝t_k-t₀，

为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；

f^b(t_i)为当前时刻的速度增量。

在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，载体坐标系到载体惯性坐标系的姿态转移矩通过如下方法获得：

(3.1)、初始对准时刻t₀，载体坐标系到载体惯性坐标系的姿态转移矩阵表示如下：

其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t₀)＝[1000]；

(3.2)、根据t₀时刻的四元素Q(t₀)和t₁时刻的角速度增量获得t₁时刻的四元素其中，Φ＝|Φ|；

(3.3)、根据t₁时刻的四元素Q(t₁)和t₂时刻的角速度增量获得t₂时刻的四元素Q(t₂)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q₀ q₁ q₂ q₃]，根据Q(k)计算如下：

在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(4)中差分GPS数据包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(V_E，V_N，V_U)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。

在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(5)中根据捷联惯性系统误差方程，选取的状态向量X_INS为13阶，具体表示如下：

其中：δL为纬度误差；

δλ为经度误差；

δv_e、δv_n分别为东速误差和北速误差；

φ_e、φ_n和φ_u分别为三个姿态误差角；

ε_x、ε_y和ε_z为激光陀螺的零位；

和为加速度计零位。

在上述适用于捷联式海洋重力仪的重力异常数据处理方法中，步骤(6)中修正后的东北天向比力值f^n'表示如下：

其中：φ×为反对称阵，

为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；

f^b(t_i)为当前时刻的速度增量；

Δ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系)为右前上坐标系，记为ox_by_bz_b，原点位于重力仪质心，x_b、y_b和z_b分别指向重力仪右、前、上。将捷联式重力仪安装于测量船舶，重力仪Y轴指向船头即前向，此时，重力仪XYZ轴分别指向船舶的右前上方向，如图3所示为本发明所针对的捷联式重力仪的硬件安装示意图。

将捷联式海洋重力仪安装于船舱内，系统的Y轴指向船首，捷联式海洋重力仪记录船舶的速度增量和运动角速度增量，并投影至载体坐标系下得到速度增量f^b和角速度增量

(2)采用了惯性凝固假设，进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵所述姿态转移矩阵为对准时间段内最后时刻k时刻的姿态转移矩阵。

对载体坐标系(b系)到实际数学平台坐标系(n系)的姿态转移矩阵进行初始化设定；海洋重力数据处理中，导航坐标系通常选为地理坐标系。如图2所示为本发明载体坐标系、导航坐标系和实际数学平台坐标系之间关系示意图；

该方法用到的坐标系定义如下：

a)经线地球坐标系e：原点位于地心，oz_e轴沿地球自转轴方向，ox_e轴位于赤道平面内，从地心指向重力仪所在点经线，oy_e轴在赤道平面内，ox_e、oy_e、oz_e轴构成右手坐标系。

b)经线地心惯性坐标系i：在对准起始时刻t₀时刻将经线地球坐标系ox_ey_ez_e惯性凝固后形成的坐标系。

c)导航坐标系n’：原点位于捷联式重力仪中心，ox轴指向东(E)，oy轴指向北(N),oz轴指向天(U)。

d)实际数学平台坐标系n：坐标系Ox₁y₁z₁，近似指向东北天，与理想导航坐标n’系之间存在失准角，例如，水平失准角(φ_e、φ_n)为0.005°，方位失准角φ_u为0.08°。

e)载体系b系：坐标系Ox_by_bz_b，原点位于重力仪质心，x_b、y_b和z_b指向重力仪右前上。

f)载体惯性坐标系i_b0：在对准起始时刻t₀时刻将载体坐标系ox_by_bz_b经惯性凝固后的坐标系。

动基座初始对准算法中，姿态阵分散成4个矩阵求取。设对准点的纬度为L，则姿态转移矩阵通过如下方法获得：

其中：

其中：为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ω_ie为地球自转角速度，t₀为对准时间段内的初始对准时刻，t_k为对准时间段内的任意时刻；

为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：

式中：q₀ q₁ q₂ q₃为对准数据段最后时刻k的四元素。

载体坐标系到载体惯性坐标系的姿态转移矩通过如下方法获得：

(a)、初始对准时刻t₀，载体坐标系到载体惯性坐标系的姿态转移矩阵表示如下：

其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t₀)＝[1000]；

(b)、根据t₀时刻的四元素Q(t₀)和t₁时刻的角速度增量获得t₁时刻的四元素其中，Φ＝|Φ|；

(c)、根据t₁时刻的四元素Q(t₁)和t₂时刻的角速度增量获得t₂时刻的四元素Q(t₂)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q₀ q₁ q₂ q₃]，根据Q(k)计算如下：

为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：

其中：g为地球重力值，分别计算Vⁱ(t_k1)和Vⁱ(t_k2)的值，t_k1和t_k2分别是对准时段内的两个时刻；

Δt_k＝t_k-t₀，

为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；

f^b(t_i)为当前时刻的速度增量。

(3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵获得船舶导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)，其中t_i为当前时刻；具体方法如下：

根据导航过程中起始时刻，即第k+1时刻的速度增量f^b(t_k+1)、角速度增量和所述姿态转移矩阵获得导航过程中船舶在第k+1时刻的姿态转移矩阵东速v_e(t_k+1)、北速v_n(t_k+1)、纬度lat(t_k+1)和经度lon(t_k+1)；

根据船舶在第k+1时刻的姿态转移矩阵东速v_e(t_k+1)、北速v_n(t_k+1)、纬度lat(t_k+1)和经度lon(t_k+1)，以及第k+2时刻的速度增量f^b(t_k+2)、角速度增量获得的船舶在第k+2时刻的姿态转移矩阵东速v_e(t_k+2)、北速v_n(t_k+2)、纬度lat(t_k+2)和经度lon(t_k+2)；

根据船舶在第k+2时刻的姿态转移矩阵东速v_e(t_k+2)、北速v_n(t_k+2)、纬度lat(t_k+2)和经度lon(t_k+2)，以及第k+3时刻的加速度增量f^b(t_k+3)、角速度增量获得的船舶在第k+3时刻的姿态转移矩阵东速v_e(t_k+3)、北速v_n(t_k+3)、纬度lat(t_k+3)和经度lon(t_k+3)。

依次类推，获得船舶在导航过程中当前时刻的姿态转移矩阵东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)。按照上述方法可以获得船舶在导航过程中各时刻的姿态转移矩阵东速v_e(t)、北速v_n(t)、纬度lat(t)和经度lon(t)。

(4)、根据船舶测量期间记录的GPS移动站数据和精密星历，采用waypoint软件中的PPP技术，获得差分GPS信息，包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(V_E，V_N，V_U)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。获得差分GPS数据后，进行步骤(5)。

(5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速v_e(t_i)、北速v_n(t_i)、纬度lat(t_i)和经度lon(t_i)，以及步骤(4)中差分GPS数据中的当前时刻的东速V_E，北速V_N、经度λ和纬度L，得到相应东速的差值v_e(t_i)-V_E、北速的差值v_n(t_i)-V_N、经度的差值lon(t_i)-λ和纬度的差值lat(t_i)-L，作为卡尔曼滤波器的观测量进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计。

根据kalman滤波器估计修正各项参数误差，特别是位置、速度和姿态误差，修正分为开环和闭环修正两种方式。完成后，进入步骤(6)。

选取的状态向量X_INS为13阶，具体表示如下：

根据测量环境，忽略部分非主要误差参数，采用了如下的捷联惯性系统的误差方程

式中，δL为纬度误差；

δλ为经度误差；

δv_e、δv_n分别为东、北速误差；

φ_e、φ_n和φ_u分别为三个姿态误差角，通常情况下，φ为小量；

ε_x、ε_y和ε_z为激光陀螺的零位；

和为加速度计零位；

T_ij(i＝1,2,3；j＝1,2,3)为姿态阵的元素。

(6)、根据步骤(5)中得到的当前时刻姿态误差，修正东北天向比力值，得到修正后的东北天向比力值f^n'，完成修正后进入步骤(7)。

其中：φ×为反对称阵，

为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；

f^b(t_i)为当前时刻的速度增量。

(7)、计算重力异常粗值δg，公式如下：

式中：g_b为码头处的重力基准值，为设定值；

f^u为f^n'中的天向比力值；

为码头处的天向比力初值；

a^u为天向运动加速度，根据GPS提供的高度二次差分得到；

δa_E为厄特弗斯改正；

δa_F为自由空间改正；

γ₀为正常重力改正；

δg_drift为零点漂移改正，零点漂移改正为每个航次重力测量时不同时间在同一点的观测值变化的改正。

其中计算重力异常各改正项，包括厄特弗斯改正、天向运动加速度、正常重力改正、自有空间改正、零点漂移改正，各改正项计算公式如下。厄特弗斯：

天向运动加速度：

可以根据GPS提供的伪距、载波相位和多普勒频移观测值，以及其单差、双差组建观测方程，利用最小二乘法，求出载体位置、速度和加速度。然后利用差分GPS数据，采用位置差分、速度差分或者载波相位差分等方法计算天向运动加速度。

正常重力：

γ₀＝9.780327(1+0.0053024sin²L-0.0000058sin²2L)

自由空间：

零点漂移改正：

重力测量的零点漂移率可按线性化近似计算，有

其中，C为此次测量的零点漂移变化率，和分别为前后校基准点处的天向比力值，和分别为前后校基准点处的重力场，t₁和t₀分别对应f₁^u和的观测时间。

则零点漂移改正值为

δg_drift＝C·(t-t₀)

式中，C为此次测量的零点漂移变化率，Δt_i为第i个测点的测量时刻与在基准点前校时刻的时间差。

完成计算后进入步骤(8)。

(8)采用数字滤波器对δg进行滤波，得到高精度重力异常信号，采用滤波器可采用FIR和IIR滤波器，截止频率小于0.01Hz；也可采用正反kalman滤波器。

(9)若船舶在行驶过程中存在形成网格的交叉点，则进行测线网平差处理方法，消除仪器的系统误差和半系统误差，提高测量精度。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。