1 引 言

现有海洋重力异常场分辨率已经达到2′×2′,进一步插值精细化处理可达1′×1′,海洋重力仪的测量精度达到1 mGal(1 Gal=10^-2 m/s²),这给高精度的重力辅助导航提供了可能性^{[1, 2, 3, 4]}。自主水下航行器重力辅助惯性导航作为一种利用地球重力场特征信息结合匹配算法对惯性导航系统位置信息进行初次估计，而后采用导航数学模型方程结合一定的滤波算法对惯性导航系统位置信息进行进一步精化修正的新技术，近年来无论是理论还是实践应用都得到了快速的发展。目前关于重力匹配算法的文献较多，研究也比较透彻，按照匹配算法设计原理可分为序列相关匹配方法和递推滤波方法两种,序列相关匹配方法主要分为最近等值线迭代算法(ICCP)和相关极值分析算法两大类^{[5, 6, 7, 8, 9, 10, 11]}。文献^[7]给出了一种基于平均平方差最小准则构造差分降权相关目标函数模型的相关极值分析算法，有效修正了重力辅助导航的误差。文献^[9]在0.2′×0.2′重力辅助异常数据库的基础上，提出将匹配位置误差作为观测量，用卡尔曼滤波对惯性导航系统进行最优估计的最近等值线迭代算法(ICCP)匹配算法。重力异常图一般是以规则格网形式预先存储在匹配计算机中，如何获取高精度、高分辨率的重力异常图是进行准确匹配的关键所在。高斯样条函数具有局部支撑的良好特性，用其作为样条基函数对计算区域重力异常进行二维整体逼近，能在保证重力异常局部特性不失真的前提下获取计算区域重力异常的统一解析式等优点。文献^[12]给出了一种基于快速傅里叶技术建立连续局部重力场的重力辅助导航算法，定位误差小于350 m。文献^[13]提出了一种基于斐波那契数列利用二维高斯样条函数逼近局部重力异常场的方法，得到高精度逼近局部离散格网数据的局部重力异常场连续解析模型，连续解析模型逼近绝对误差均值达到0.000 69 mGal。文献^[15]给出了一种基于余弦变换的Parker正演重力梯度方法。

自主水下航行器惯性导航系统中加速度计无法区分运动加速度和重力加速度，需要通过其他手段将比力中的重力加速度扣除，称为重力补偿。如文献^{[5, 8, 9, 10, 12]}所述，目前传统自主水下航行器重力辅助导航系统模型的构建过程中，认为比力中的重力加速度为标准重力模型，未考虑重力异常扰动矢量、标准重力模型误差(标准重力模型误差由测量位置维度误差解算引起)和厄特弗斯修正计算值3类误差影响(本文所述的传统模型就是指不考虑3类误差影响，认为比力中的重力加速度为标准重力模型)。这些误差都将引起导航误差，在惯性元件精度不高的系统中，惯性元件带来的系统误差较大，相比之下，正常重力模型足够精确，这些系统误差可以忽略。但在精度要求较高的场合(重力辅助导航中陀螺仪零偏小于0.001 deg/h，加速度计零偏小于0.01 mg的情况下)，尤其是随着惯性元件精度的不断提高，其带来的系统误差与上述3类误差带来的系统误差处在同一个量级，甚至前者可能低于后者，系统误差不能被忽略，有必要采用精度更高的系统状态方程模型对正常重力矢量进行补偿，应在初次重力匹配算法计算经纬度结果的基础上，对标准重力模型误差、重力扰动矢量和厄特弗斯修正计算值进行补偿，而后再进行整个系统的导航滤波计算，目前并未有这方面的研究^{[15, 16, 17, 18, 19, 20]}。因此，本文基于二维高斯样条局部重力异常场提出了一种考虑重力扰动矢量、标准重力模型误差、厄特弗斯修正计算值的精化自主水下航行器重力辅助导航模型构建算法。

2 高斯基函数插值

高斯样条基函数G(x)=e^(－x2/a2)，a为局部支撑参数，设局部重力异常图等分格网点为{(x_i,y_j,z_i,j)|i=1,2,…,m;j=1,2,…,n}，Δx、Δy为格网点分辨率。这样符合上述x、y平面的二维高斯基函数可写成矩阵形式

式中，L_x=[G_x((x－x₁)/Δx)G_x((x－x₂)/Δx…G_x((x－x_m)/Δx)]；L_y=[G_y((y－y₁)/Δy)G_y((y－y₂)/Δy…G_x((y－y_n)/Δy)]；简记为L_x=Inte{G_x((x－x_i)/Δx)}，L_y=Inte{G_y((y－y_i)/Δy)}

根据插值条件

{L(x_i,y_i)=f(x_i,y_i) i=1,2,…，m;j=1,2,…，n}，

将矩阵组合，则有如下线性方程成立

式中，Z=(Z_ij)_m×n；Z_ij=L(x_i,y_i)；

可见，高斯基插值函数存在两个重要参数需要确定：①插值系数矩阵λ；②局部支撑参数a。无论是一维高斯基函数还是二维高斯基函数，如果X、Y均为非奇异矩阵，则式(2)有唯一解，对X、Y矩阵求逆，解方程即得到插值系数矩阵，将插值系数矩阵λ代入式(1)便得到该局部重力异常基准图二维高斯基函数逼近解析式。研究证明，局部支撑参数a的取值对高斯基函数逼近效果影响很大，a过大或过小都将导致很大的逼近误差。当a过小时，由于局部支撑区间很小，使得插值曲线刚度小，容易出现锯齿状；当a过大时，由于局部支撑区间很大，插值曲线刚度大，插值曲线过于平滑甚至使插值解析式在已知点的插值误差很大，所以a太大或太小均会失去插值的意义。因此有必要对参数a进行寻优，结合插值模型精度评估方法，解算参数a的最优值。文献^[13]中的二维高斯函数插值局部寻优方法对a进行寻优，建立局部连续重力场。

3 基于高斯插值的精化算法 3.1 速度误差方程

取地理坐标系n为导航坐标系；b表示载体坐标系；C_bⁿ表示由载体坐标系向导航坐标系转换。假设姿态误差角φ是小角度，在考虑因重力扰动矢量、标准重力模型误差和厄特弗斯修正计算值的情况下，重力辅助导航速度误差方程的具体推导叙述如下：

在不考虑任何误差的情况下，认为根据导航仪表所示载体位置为，经过重力匹配算法初次计算得到的实际位置为；根据比力方程，速度理想方程由下述方程(3)给出

考虑各种误差的情况下，实际速度计算值应由下述方程确定

式中，v^c=vⁿ+δvⁿ；ω^c_ie=ωⁿ_ie+δωⁿ_ie；ω^c_en=ωⁿ_en+δωⁿ_en；(I+[δA])+^b；

[δK_A]=diag[δK_Ax δK_Ay δK_Az]；φ_E、φ_N、φ_U为姿态误差角；δK_Ai和δA_i(i=x,y,z)分别为加速度计的刻度系数误差和安装角误差。

式(4)与式(3)相减，略去二阶小量可得线性近似的捷联惯导误差方程为

式中，δg=[_m(x,y)－E_C(x,y)+；δfⁿ_sf=Cⁿ_b{(δA+[δK_A])f^b_sf+^b}；_m(x,y)为前述高斯插值计算得到的连续重力异常场表达式；g_s(x)、E_C(x,y)分别为后文叙述的标准重力模型表达式和厄特弗斯修正计算值。

3.2 状态方程

除速度误差方程外，本文所提导航系统的其他误差方程与其他文献中的误差方程相同^{[20, 21, 22, 23]}，定义SINS的状态变量如下，共15维状态矢量X_SINS=[φ^T δv^nT δp^T ε^{bT bT}]^T，其中，φ^T=[φ_E φ_N φ_U]；δv^nT =[δvⁿ_Eδvⁿ_Nδvⁿ_U]；δp^T=[L λ h]分别为惯导系统姿态角、速度、位置误差。将加速度计误差模型近似为随机常值零偏加白噪声，陀螺误差模型近似为随机常值漂移加白噪声，这样陀螺仪漂移误差ε^bT、加速度偏移误差Δ^bT方程见式(6)

于是可得到基于高斯插值的导航状态方程

式中，系统噪声w_SINS=[w_gx w_gyw _gz w_ax w_ay w_az]^T为6维系统干扰白噪声;系统噪声分配矩阵为G_SINS是15×6矩阵；状态矩阵F_SINS是15×15矩阵。

3.3 量测方程

本文系统量测值取重力仪量测重力与惯导指示位置对应连续场重力图重力之差。设定此时系统的量测方程Z_g可由式(8)计算得出

式中，g₀=[0 0 978 049]^T，单位为mm/s²。

设导航仪表所示载体位置为(L+ΔL,λ+Δλ)；对应的实际位置为(L,λ)；_m(L+ΔL,λ+Δλ)为在重力图上根据载体测量位置(L+ΔL,λ+Δλ)读出的重力异常值；_s(L+ΔL)为根据载体纬度测量位置信息L+ΔL利用标准重力公式(9)计算得到的重力值；(L,λ)为重力传感器输出的实际重力值；E_C(L+ΔL,λ+Δλ)为利用式(10)计算得到的厄特弗斯修正计算值。

厄特弗斯修正计算值的计算公式为^[1]

式中，vⁿ_E和vⁿ_N为载体的东向和北向速度；L为载体当地的大地纬度；Ω=7.292 115 67×10^－5rad/s表示地球自转角速率。

由于绘制成的重力图值等于实际重力测量值经过厄特弗斯修正而后减去标准重力值，所以对于实际位置点(L,λ)，有^[3]

则(L,λ)=_m(L,λ)－E_C(L,λ)+g_s(L)，将其代入方程式(8)，经线性化处理后可得系统量测方程

式中，V_gravity为重力传感器噪声、重力异常图制作噪声和线性化误差总和，在实际计算中，重力异常梯度可以采用9点拟合、全平面拟合等方法计算。

测量矩阵为

式中，Z_gline为式(8)计算得到的重力测量误差值；H_line为1×15维矩阵，在矩阵H_line中，H_line(1,4)=；其余元素均为0。V_gravityline为重力传感器的噪声、重力异常图制作噪声和线形化误差总和。

4 仿真试验与分析

本文采用第二炮兵工程大学惯性导航实验室基于Windows 2000操作系统、Matlab 7.0软件自行开发研制的惯性导航轨迹生成软件，使用代码为real-trace和inertial navigation-trace的两个函数调用相应的子库函数，仿真生成需要的真实航迹和惯导指示航迹^[23]。而重力异常数据使用渤海地区经度范围118°E-119°E，纬度范围37°N-38°N的真实数据，该重力异常场的三维图如图 1所示。重力图分辨率为1′×1′，重力仪的实时量测数据是利用真实航迹在重力数据库中的采样值叠加0.4 mGal的随机量测噪声获取，陀螺常值漂移取0.01°/h，随机噪声漂移为0.01°/h，加速度计随机常值零偏1×10^－4 g，随机噪声漂移为1×10^－5 g。初始状态：速度10 m/s、经度118.02°E、纬度37.24°N、惯导测量经纬度误差0.02°、对准水平姿态误差角2′、对准航向误差角10′。运动载体保持俯仰角和横滚角10°、航向角30°的姿态，0.005 m/s²加速度每隔10 min加速一次，3 min采样匹配滤波一次，行进了70 min，得到24个测量点^[23]。在捷联惯导的情况下进行仿真，对比分析了ICCP算法初次估计解算，局部连续场情况下(文献^[13]进行高斯插值，最优局部支撑参数a_best=1.535 7)传统重力辅助导航模型(认为比力中的重力加速度为标准重力模型值，未考虑重力异常扰动矢量、标准重力模型值误差和厄特弗斯修正计算值3类误差影响建立重力辅助导航系统模型)和本文的基于高斯插值的精化重力辅助导航模型构建算法的效果(分别称其为tradition-ICCP和Gauss-compensation-ICCP算法，简记为TRAI和GACOI算法)。

图 2(a)-(c)分别为TRAI和GACOI算法导航经、纬度位置误差；航向、俯仰、横滚姿态误差角；东、北、天向速度误差分布图。如图 2所示；TRAI算法导航位置误差在125～330 m、平均为210 m，水平姿态误差在0.2′～0.75′、平均为0.50′，航向姿态误差在4′～9′、平均误差为7.5′，速度误差在0.027～0.085 m/s、平均误差为0.045 m/s。GACOI算法导航能够保持在100～200 m、平均110 m的位置误差，0.01′～0.25′、平均0.18′的水平姿态误差，2′～6′、平均2.5′的航向姿态误差，0.002～0.05 m/s、平均0.016 m/s的速度误差。可见GACOI算法将传统重力辅助惯性导航算法的位置导航精度提高了1倍左右，姿态角导航精度、速度导航精度提高了1～2倍。新算法由于引入了先进的局部支撑模型参数确定方法，同时考虑重力扰动矢量、标准重力模型误差、厄特弗斯修正计算值，提高了导航精度，达到了近乎高精度导航的水平。

5 结 论

本文首先对高斯样条优化插值以及所涉及的局部支撑参数求解问题进行研究探讨，给出了局部连续重力异常场的构建方法，而后基于此给出了考虑重力扰动矢量误差、标准重力值误差、厄特弗斯修正计算误差的高斯插值精化导航算法，最终提高了重力辅助导航系统的导航精度，实现了精化导航计算的目的。本文算法可以为重力辅助惯性导航算法的实践化提供一定参考，然而在重力辅助导航系统观测量为非线性的情况下，如何基于局部连续重力异常场进行相应的重力辅助导航非线性滤波算法设计还有待进一步研究。