1 引 言

惯性导航系统(inertial navigation system) 的定位误差随时间积累发散，无法长时间保持高精度，为了满足潜艇长期隐蔽潜航的要求，有必要借助另外一种无源的辅助导航手段对INS位置进行校正。重力匹配辅助惯性导航是一种利用地球重力场特征信息对惯性导航系统位置漂移误差进行修正的导航新技术^[1,2]。重力匹配辅助导航系统由惯性导航系统、重力图、重力传感器、匹配计算机等组成。载体在运动过程中，重力传感器可准实时测得重力数据，同时根据INS的位置信息，可以从重力图中读取重力数据。将这2种数据输入匹配解算计算机，利用匹配定位软件进行解算，即可求得最佳匹配位置^{[3, 4, 5, 6]}。利用匹配位置信息对INS系统进行校正，即可起到抑制INS定位误差，延长INS重调周期的作用。重力匹配辅助惯性导航系统中代表性产品包括美国Lockheed Martin公司研制的由通用重力模块(UGM)所组成的重力导航系统(NGS)和贝尔实验室研制的重力辅助导航系统(GAINS)。

重力异常图一般是以规则格网形式预先存储在匹配计算机中，如何获取高精度、高分辨率重力异常图是进行准确匹配的关键所在。目前国内外主要采用逆Stokes公式、逆Vening-Meinesz公式对Geosat、Seasat、ERS-1和T/P等卫星测高数据进行反演获取海洋格网重力异常模型^{[7, 8, 9, 10]}，分辨率达到了2′×2′。为了进一步提高重力异常图分辨率，许多学者采用二维插值的方法对原重力异常图进行精化^{[11, 12]}，二维插值常用方法有最近邻法、最小曲率法、距离倒数加权法、多项式回归法、Kriging法、改进的Shepard方法以及径向基函数法等^{[13, 14]}。采用函数逼近对重力异常图进行解析重构使重力异常图的应用不受分辨率限制，在重力异常函数逼近方面，相关学者也做出了大量有价值的研究^{[15, 16, 17, 18]}。文献[15]对各种函数模型逼近及统计模型逼近的特性进行了分析，提出由函数模型逼近与统计逼近组合起来的综合逼近法，克服了2种逼近方法的不足，达到了较好的效果。文献[16]系统地建立了地球重力场逼近中的B样条表示理论和算法。文献[17]则采用快速傅里叶变换方法对局部地球重力场模型进行逼近。文献[18]提出了1种利用Fourier级数展开表示的局部地球重力场模型，该方法获得的最终解析表达式是连续的。综合来看，B样条函数逼近及傅里叶级数逼近两种方法都具有很高的精度，但由于B样条函数逼近最终解析表达式不统一，因而不适用于连续场匹配导航，而傅里叶级数逼近涉及大量参数计算，限制了连续场匹配算法实时性。高斯样条函数在医学影像、数字地形构建、计算机图形学及地球物理等领域有较为广泛的应用，文献[19]讨论了一维高斯样条在协方差函数代数确定中的应用，并对平稳随机过程逼近中一维高斯样条最优参数确定进行了研究。本文考虑到高斯函数具有局部支撑的良好特性(其函数特性与B样条函数非常接近)，采用高斯函数作为样条基函数对计算区域重力异常进行二维整体逼近，因此能在保证重力异常局部特性不失真的前提下获取计算区域重力异常的统一解析式，同时考虑到高斯样条函数逼近局部重力异常的相关影响因素，采用斐波那契寻优方法获取样条函数最优参数^[20]。采用高斯样条函数逼近局部重力异常运算简单且最终解析式是统一的，因此可以较好的满足了连续场匹配导航的要求，具有一定的应用价值。

设高斯样条函数表达式为G(x)=exp(－x²/a²)。高斯函数衰减很快，它可以看作局部支撑的基函数，但它的局部支撑区间受到了参数a的控制，局部支撑区间随a增大而增大。在一维函数y=f(x)的X方向上以Δx等距采样获得离散点集合{(x_i,y_i)|i=1,…,n}，则在X方向上的高斯样条函数表达式为，在已知离散点处满足插值条件：{L(x_i)=f(x_i)|i=1,…,n}，则根据插值条件有如下线性方程：

点击浏览公式

易知A是非奇异的，即上式有唯一解，解方程即得到高斯样条函数一维插值解析式。参数a对插值影响很大，当a过小时，由于局部支撑区间很小，使得插值曲线刚度小，容易出现锯齿状；当a过大时，由于局部支撑区间很大，插值曲线刚度大，插值曲线过于平滑甚至使插值解析式在已知点的插值误差很大。由上述论述可以知道，a太大或太小均失去插值意义，因此有必要对参数a进行寻优。由上式可以看出插值由A与b唯一确定，而A由矩阵的阶数和a决定，而在求解方程的过程中存在数值计算误差，如A阶数过高或a过大，对A的求逆存在较大误差。下面计算对不同阶数和a值时A矩阵求逆数值计算误差，设I为单位矩阵，B=inv(A)A，误差计算公式如式(2)，误差分析如图 1所示。

点击浏览公式

图 1 A矩阵求逆误差分析图 Fig. 1 The matrix inversion error analysis of A

图选项

由A矩阵求逆误差分析图可以看出，当A矩阵阶数小于5或a值小于5时，A矩阵求逆运算误差较小。考虑到实际应用中对局部重力场逼近时用到的格网区域远大于5×5，这里仅对a∈[1,5]范围内进行一维变量寻优。

斐波那契数列方法是一种经典的一维寻优方法，它能以较少的迭代次数获取精度较高的最优解，且只需计算和比较函数值，算法简单易于实现，因此采用斐波那契数列方法对参数a进行寻优。为了保证插值函数与已知数据点保形性良好，采用式(3)准则函数进行寻优。

点击浏览公式

设原函数曲线为f(x)=10sin(x)+10,x∈[0,6]，以Δx=0.2在该曲线上进行等距采样获得已知离散点集，由已知数据点通过准则函数斐波那契数列寻优获得a最优值为2.8909。为说明对参数a寻优的有效性。分别令a=1、2、2.890 9、4、4.5对Δx=0.1等距离散点处进行插值并与该点原函数值相比，插值误差如图 2所示。由仿真误差分析图可以看出，当a值较小时插值曲线刚度小，在区域边缘龙格现象明显；当a较大时，曲线刚度变大，虽龙格现象得到有效抑制但整体误差变大。综合来看，取a为最优值时的插值误差及龙格现象均得到有效抑制。一维高斯样条函数插值误差数值分析结果如表 1所示，从误差分析结果可以看出a取最优值时其插值绝对误差均值及绝对误差均方差均为5组试验中的最小，从而说明对参数a进行寻优是有效的。最后采用B样条函数对已知数据进行仿真试验，从图 2及表 1可以表明采用B样条函数进行插值边缘误差较大，总体精度较好，但其精度远低于最优高斯样条函数插值。

图 2 一维高斯样条函数插值仿真误差分析图 Fig. 2 The analysis of 1-D Gauss spline interpolation simulation errors

图选项

下面在一维高斯样条函数插值基础上讨论二维高斯样条插值。已知XY平面上等分离散格网点集合{(x_i,y_j)|i=1,…,m;j=1,…,n}及其对应函数值z_i,j(i=1,…,m;j=1,…,n)，则其X方向一维高斯样条函数解析式为L_x=Span{G_x((x－x_i)/Δx)}，Y方向一维高斯样条函数解析式为L_y=Span{G_y((y－y_j)/Δy)}，XY平面的二维高斯样条函数可以写成L_x与L_y的张量积形式：

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

易知X、Y均为非奇异矩阵，即上式有唯一解，解方程即得到高斯样条函数二维插值解析式。由上一节论述可知参数a对插值影响很大，因此有必要在二维高斯样条函数插值前对X,Y方向的a_x、a_y进行寻优。同上考虑到求逆计算误差，仅需对a_x∈[1,5],a_y∈[1,5]区间内进行寻优，理论上必须对a_x、a_y进行二维联合寻优。考虑到二维准则函数构建复杂且需进行二重积分影响算法速度，将二维寻优模型解耦为X方向的a_x寻优模型及Y方向的a_y寻优模型。解耦方法为：将格网点划分为n行，对每行分别进行一维高斯样条函数插值，将这n次一维插值准则函数之和作为二维高斯样条函数插值X方向的a_x寻优准则函数；将格网点划分为m列，同样对每列分别进行一维高斯样条函数插值，将这m次一维插值准则函数之和作为二维高斯样条函数插值Y方向的a_y寻优准则函数；二维高斯样条函数插值准则函数解耦后表达式如式(9)、式(10)所示。根据准则函数解耦模型我们即可采用斐波那契数列方法对a_x、a_y分别进行寻优。

点击浏览公式

点击浏览公式

为方便仿真验证，直接采用MATLAB里的高斯合成曲面函数peaks函数组合成局部连续重力异常场作为仿真用真实背景场，该重力异常场解析表达式如式(11)，其3D示意图如图 3所示。

点击浏览公式

图 3 局部连续重力异常场3D示意图 Fig. 3 The 3D map of local continuous gravity anomaly field

图选项

以分辨率0.2×0.2在该曲面上进行等距采样获取31×31标准离散格网点集，该区域重力异常变化范围为-51.185~86.1819mGal，在已知离散格网点集的基础上采用斐波那契数列方法以10^－2精度对二维高斯样条函数插值解耦准则函数分别进行寻优得到a_x的最优值为3.5391及a_y的最优值为2.5272。按表 2给定的5组a_x、a_y值对仿真区域0.1×0.1等距离散格网点进行高斯样条函数二维插值，插值结果与对应点的函数值之差作为该点的插值误差，5组试验插值误差图如图 4所示。

图 4 二维样条函数逼近仿真误差分析图 Fig. 4 The analysis of 2D spline approximation simulation errors

图选项

从前4组仿真试验插值误差区域分布来看，当a_x、a_y值较小时，插值曲面刚度小，较大误差都分布在边缘很小的一块区域内(约1～2原始格网点宽)，在内部区域的插值精度都非常高，且a_x、a_y值越小边缘效应明显，因此在构建整个海区的高精度重力异常基准图时，有必要剔除边缘区域的异常误差数据。由后2组仿真试验插值误差图可以看出，当a_x、a_y值较大时，插值曲面刚度大，虽边缘效应得到有效抑制，但在内部区域误差过大甚至导致插值算法发散。由第3组仿真试验插值误差图可以看出，当a_x、a_y值为斐波那契数列寻优最优值时，插值误差及边缘效应均得到有效抑制。综合来看，由于已知数据沿X方向与Y方向的数据特征不同，a_x、a_y一般是不等的；a_x、a_y取过大或过小均可能导致较大误差，因此在插值前对a_x、a_y进行寻优是非常有必要的。5组试验的误差分析结果如表 3所示，从误差分析结果也可以看出a_x、a_y取最优值时其插值绝对误差均值为5组试验中的最小，从而说明在二维高斯样条函数插值时对参数a进行解耦寻优是有效的。最后设计试验6对B样条函数二维逼近进行仿真，仿真误差分析见图 4(f)及表 3。从试验结果可以看出，采用B样条对格网数据进行函数逼近时精度较好，但与最优二维高斯样条逼近有相当大的差距，从而进一步说明了最优二维高斯样条逼近的有效性。

为了验证该方法在实际重力异常数据处理中的有效性，取范围为λ:122°E~126°E,φ:20°N~24°N，分辨率为2′×2′的卫星测高反演重力异常数据作为仿真用真实重力异常基准图进行二维高斯样条逼近计算。由于在除格网点外其它位置真实重力异常并不知道，为了方便进行逼近误差分析，采用该方域4′×4′格网处重力异常数据作为已知数据点进行二维高斯样条函数逼近，获得该范围局部重力异常基准图解析表达式，在此基础上计算2′×2′格网点处重力异常逼近值，然后与已知的2′×2′格网处重力异常值进行比对分析。如图 5(a)及图 5(b)所示，二维高斯样条函数逼近的局部重力异常基准图解析式能较好的描述这一范围的重力异常基准图。从图 5(c)、图 5(d)及对应的误差分析可以得出，试算区域4′×4′格网点处二维高斯逼近计算值绝对误差均值为0.229 6

图 5 局部重力异常基准图二维高斯逼近仿真试验 Fig. 5 The simulation of local gravity anomaly map’s approximation with 2D Gauss spline

图选项

图 6 局部重力异常基准图B样条逼近仿真试验误差分析< Fig. 6 Error analysis for the simulation of local gravity anomaly map’s approximation with B-spline

图选项

对局部离散格网重力异常场进行解析重构是研究连续场重力匹配导航的前提条件。为了适应重力匹配导航的实际应用，本文提出了利用高斯样条函数采用张量积算子对离散重力异常场进行逼近。分析表明该方法能获得局部范围内统一解析式且计算简单。需要强调的是，重力场解析重构的最终精度与原始数据分辨率及数据质量是密切相关的，而本文提出的重构方法无法消除由原始数据分辨率及其精度带来的影响。本文最后对仿真数据及实际数据进行了试算，结果表明高斯样条逼近法所获得的精度能够满足匹配导航的要求，如何在高斯样条逼近模型的基础上展开重力匹配算法解析建模有待进一步研究。