CN108415096B - 基于牛顿迭代法的水下重力梯度目标探测方法 - Google Patents

Abstract

本发明公开了一种基于牛顿迭代法的水下重力梯度目标探测方法，利用水下潜器的重力梯度仪对航行的区域进行实时测量，获取实测的重力梯度值；利用实测的重力梯度值去减去正常的重力梯度值，得到重力梯度异常值；将重力梯度异常值作为重力梯度模型非线性方程组的输入，求解非线性方程组获目标体的质心相对于探测点的相对位置；由目标体的质心相对于探测点的相对位置，确定目标体的质心与探测点之间的距离；最后将上述距离和万有引力常量代入重力梯度垂直分量公式，获取目标体的质量；还能根据需要推测出密度异常体的大致物质种类。具有防伪装，抗干扰特性，探测方法简单通俗、使用方便快速、大大减少了人工的干预，实现了自动探测。

Description

基于牛顿迭代法的水下重力梯度目标探测方法

技术领域

本发明属于水下目标探测技术领域，具体涉及一种利用重力场相关辅助的无源导航目标探测方法。

背景技术

现今对于水下导航的方式有很多种，例如有GPS卫星导航、天文导航、无线电导航等，基于重力梯度导航(此处指竖直方向重力梯度)不需要接收外部的信息，同时也不向外部辐射信息进而难以进行干扰等诸多优点，是严格意义上的无源导航，具有十分重要和可行的现实意义，利用重力场相关辅助导航避障已经成为了现代导航避障技术研究的热点。

利用重力梯度仪进行辅助导航开始于20世纪60年代，从最初的研究集中在实时检测重力异常以提高惯性导航系统的精度，发展为基于重力梯度仪的高精度高分辨率的重力图研究，利用高精度的重力图和重力测量仪器便可以对惯性导航系统进行辅助修正。同时，重力梯度作为重力位的二阶导数，经过前人多番验证，重力梯度具有比重力本身更高的分辨率，且同一位置有多个分量。很多情况下，由于重力位与重力梯度都受测量点到质心距离和密度的影响，其中影响最大的一般都是垂直分量。所以我们可以选择单一分量(垂直分量)来进行水下导航避障。

目前，重力梯度张量测量在国外已经由理论研究迈入到十分成熟并进入商业应用阶段，而在我国还处于起步阶段，但重力梯度张量数据的处理和解释已经得到了开展。除了用重力梯度仪直接观测重力梯度张量，还可由重力异常数据经过一定的算法求得重力梯度张量，更直接的办法是直接用模型正演计算简单模型研究复杂地质体。国内、外学者已经研究和提出了众多重力梯度的求取算法，例如样条函数插值法，差商方法、频率域中的傅里叶变换法和余弦变换法以及有限单元方法等。

目前有AGTI算法对目标进行探测，探测质量难以确定，并且AGTI算法的质量误差会带入到距离探测中，同时，AGTI算法需要重力梯度全张量或5个张量进行计算，计算量大效率不高，同时它的最大缺陷是质量的探测误差很大，它的质量探测与积分范围有关，积分范围越大，质量的误差会越小，而积分越大，计算时间越长。这样计算时间越长，要等好长时间才能探测出目标体信息，它的实时性就会很差，其次他的距离R的探测依赖于质量，导致质量的误差会引入到距离中，那么它的位置探测也会有很大误差。

发明内容

本发明要解决的技术问题是：针对目前水下潜器对目标探测时精度不高、方法复杂的问题，提供一种基于牛顿迭代法的水下重力梯度目标探测方法，可以准确高效获取目标质量和位置，并大致判断物质种类而不依赖于形状，具有防伪装，抗干扰特性，探测方法简单通俗、使用方便快速、大大减少了人工的干预，实现了自动探测。不仅仅适用于战时水下水下潜器(潜艇)的探测，而且还可适用于水下机器人的目标探测。

为解决上述技术问题，本发明采用如下技术方案：

一种基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于包括如下主要步骤：

S1：利用水下潜器的重力梯度仪对航行的区域进行实时测量，获取实测的力梯度值；

S2：利用实测的重力梯度值去减去正常的重力梯度值，得到重力梯度异常，即为目标体或密度异常体引起的重力梯度值；

S3：将步骤S2得到的目标体或密度异常体引起的重力梯度值作为重力梯度模型非线性方程组的输入，求解非线性方程组获目标体的质心相对于探测点的相对位置；由目标体的质心相对于探测点的相对位置，确定目标体的质心与探测点之间的距离；

S4：最后将上述距离和万有引力常量代入重力梯度垂直分量公式，获取目标体的质量,完成水下目标探测。

进一步的，在步骤S4之后，还可以进行步骤S5：根据目标体的质量和前述测得的目标体或密度异常体引起的重力梯度值，计算出目标体的平均密度，根据平均密度大致判断目标体物质种类。

进一步的，所述目标体为海底山峰、敌方潜艇、礁石、沉船之类的至少一种物质，当它们连成一体或重叠时将他们看成一个整体处理。

进一步的，步骤S3中，将重力梯度模型非线性方程组化简，选取三个张量求解，得到目标体距离探测点在三维坐标系中各个方向的相对距离。

进一步的，步骤S3中，选取xx,yy,zz三个张量，分别是重力位对x方向的二阶导数，y方向的二阶导数，zz方向的二阶导数。

进一步的，步骤S3中，非线性方程组化简如下：

其中x,y,z是目标体质心与测量点在x,y,z三个方向的相对距离，l是相邻测量点的步长， Γ_xx ^(t)是t时刻xx方向的重力梯度值，Γ_yy ^(t)是t时刻yy方向的重力梯度值，Γ_zz ^(t)是t时刻zz方向的重力梯度值；通过求解上述方程组，便可以得到x,y,z的值，获取目标体与探测点在三个方向的相对距离。

进一步的，步骤S3中，根据以下公式获取目标体质心与测量点的距离：

R代表目标体质心与测量点的距离，x,y,z是目标体质心与测量点在x,y,z三个方向的相对距离。

进一步的，步骤S3中，采用经典的牛顿迭代法求解非线性方程组。

进一步的，步骤S4中，目标体的质量m通过如下重力梯度垂直分量公式计算得到：

其中G为万有引力常量，Γ_zz是由目标体或密度异常体引起的重力梯度值，Δz为目标体z方向到探测点的相对距离。

本发明通过分析重力梯度的基本原理，设计一种重力梯度模型，将水下目标探测问题转化为数学问题，对目标在三维坐标系上相对位置的求取转化为重力梯度模型非线性方程组的求解；非线性方程组的求解采用经典的牛顿迭代法，再将求的值带入重力梯度垂直分量来反求质量，完成水下目标的探测。

基于牛顿迭代法的水下重力梯度目标探测方法，是首先通过重力梯度仪测得重力梯度值，此处说的重力梯度值均为密度异常体造成的重力梯度值。再带入化简得到一个非线性方程组，那么目标的质心相对于探测点的相对位置变成了求解非线性方程组的问题。对于非线性方程组的问题，采用经典的牛顿迭代法，可以很快的借助计算机求解得到，三维坐标系的相对位置确定后，距离也就可以确定，万有引力常量已知，最后由重力梯度垂直分量公式，将已知量和已经求取的量带入，就可以获取目标体的质量。进而完成了目标探测的任务。

相对于现有技术，本发明有益效果在于：

(1)本发明与其他水下目标探测方法相比，属于无源目标探测，具有隐蔽安全，精确有效，稳定性好，方法简单等优点。

(2)本发明先探测相对位置，再探测质量，避免了其他方法因质量探测的误差引入位置探测从而导致位置探测精度不足的问题。

(3)本发明专注于目标的总质量和相对位置的探测，这两个参数对水下潜器非常重要，还可以由重力梯度值来大致判断何种物质，关注的是总质量、平均密度，质心与测量点的距离，与他的形状无关，具有防伪装，抗干扰，还可以应用于避障的和水下矿物质勘探。

(4)本发明使用方法简便，可以快速有效的进行水下目标探测，与其他的重力梯度目标探测相比，具有简单通俗，使用方便，快速的目的，全张量中只用三个张量就可以进行水下目标探测，计算量少，且大大减少了人工的干预，对相对位置的探测转化为数学上求解非线性方程组的问题，减少了人工的工作量，实现了自动探测的目的。

附图说明

图1为探测目标的三维坐标系原理图。

图2为本发明方法目标探测的流程图。

图3为本发明的实测示例图。

图4-6依次为本发明Γ_xx,Γ_yy,Γ_zz与距离的变化关系图。

具体实施方式

探测目标的三维坐标系原理图如图1所示，根据万有引力定律，在如图1空间中两个物体，质量分别为m,M的两物体，坐标分别为(x₁,y₁,z₁),(x₂,y₂,z₂)，之间的距离为R。

两物体之间的万有引力：

在三维坐标系中位于(x₁,y₁,z₁)处质量为m的质点对位于(x₂,y₂,z₂)测量点的产生的重力位为：

此处G为万有引力常量。对上式求一阶导可得：

对式(1)求二阶导数得：

由m＝δv，其中m表示质量，δ表示密度，v表示体积。对于体积，可以得到如果密度分布不均匀，则有其中δ_i是第i个立方体的密度。

故(6)～(14)可以变成以下形式：

为与重力参数区分，

令Γ_xx＝g_xx,Γ_xy＝g_xy,Γ_xz＝g_xz,Γ_yy＝g_yy,Γ_yx＝g_yx,Γ_yz＝g_yz，

Γ_zz＝g_zz,Γ_zx＝g_zx,Γ_zy＝g_zy。同时由公式(1～15)和公式(1～23)，可以很明显的看出有以下关系：

Γ_xy＝Γ_yx，Γ_yz＝Γ_zy，Γ_xz＝Γ_zx (24)

此外，在无源空间内，引力的散度与旋度都为零，故有：

▽×g＝0 (25)

▽代表拉普拉斯算子，g代表重力梯度，上式代表重力梯度的散度即：

Γ_xx+Γ_yy+Γ_zz＝0 (26)

对于水下目标探测，其基本理论是属于反演，知道g,g_x,g_y,g_z,Γ_xx,Γ_xy,Γ_xz,Γ_yy,Γ_yz,Γ_zz这些量，来推导出探测目标体的质量m，密度ρ，位置等信息。因为我们求得质量是目标的总质量，密度是平均密度，位置是目标的质心，则有以下关系：

Δx＝x₁-x₂；Δy＝y₁-y₂；Δz＝z₁-z₂ (27)

上式中(x₁,y₁,z₁)代表测量点的坐标，(x₂,y₂,z₂)为目标体质心的坐标，Δx代表测量点与目标体质心的x方向的距离，Δy代表测量点与目标体质心的y方向的距离，Δz代表测量点与目标体质心的y方向的距离。

代入公式(3)～(14)可以推导出以下关系：

目标体探测重力梯度模型方程组的建立：

由于不是每一个区域都有重力梯度基准图进行对比：有些区域无法通过测量，有些是因为基准图年代久远未进行校准。所以此处假设一种没有基准图数据进行目标探测。

目前重力梯度或者重力基准图主要通过机载或航天设备载重力梯度仪或者重力仪获取，它们距离地面都有一定的距离，其中因为各种干扰的存在，精度相对来说较低，国外的测量仪器较比国内更为先进，得到的精度相较国内高，但基本对我国保密，目前大部分公开的精度超过1km。精度较高的数据由于在军事方面应用，故目前还处于保密阶段。对水下潜器而言，在水下航行的范围较广，在执行任务时，很多情况下，航行的区域并非我国控制的区域，要提前获得该区域的高精度数据非常困难。其次，基于基准图的重力梯度目标探测还要求水下潜器静止，这对水下潜器的运动状态要求较高。

综上考虑，本发明选择一种新的目标探测方法，这种方法既可以利用重力梯度保证水下潜器的隐蔽性，又能去除高精度重力及重力梯度基准图的要求，由此引入基于重力梯度均值比的三次测量法。

在介绍本发明的目标探测算法之前，首先，要指出这个目标探测算法中的目标的定义。这里所谓的目标可能是大型礁石或者一首沉船，只关心它的质心位置、总质量与平均密度，故只是将其看成一块整体。此处目标是判断水下潜器的航迹中障碍物的位置，以起到辅助导航避障的目的。

本发明采用重力梯度均值比的测量方法流程如图2所示。

在采用该方法进行目标探测时，水下潜器一定是保持在直线上，因为如果不在一条直线上，采样间距会不一致，这对该方法进行目标探测不适用。实际航行中，水下潜器也是基本做匀速直线运动的，对本文的方法是非常适用的。

在上图中，潜器的航是匀速直线运动，等时间间隔采样三个重力梯度值，分别记为Γ^(t),Γ^(t+1),Γ^(t+2)。只要水下潜器的速度是匀速，且保持直线运动，那么在采样时间一定时，潜器的距离也是相同的，那么这样得到的三个重力梯度值也是等间隔的。

其次两次采样求均值可以将干扰降低噪声的影响，因为每次采样都会有干扰存在，求均值会将两次的干扰进行综合。其次两次的均值一定是两个采样点之间的均值一定是某一点所对应的重力梯度值，由此来进行目标探测可以降低扰动对探测结果的影响。

由于重力梯度有6个梯度分量，也就是可以通过三次采样建立6个方程，建立的坐标系是立体的三维坐标系，而定位则是求出三维坐标系的三个坐标值。这样也就是6个方程组求三个未知数的问题。而对于三个未知数，通过数学解方程的思想，只需要3个方程就可以解出三个未知数。而对于这种非线性方程组的求解问题，目前有很多种方法，本文采用较为经典的牛顿迭代法。

步骤1：首先我们通过重力梯度仪获取一系列的重力梯度数据值，假设潜器的水下速度为10m/s，每隔50米进行一次采样。现采用如下方法进行目标探测：

假设连续三个采样的重力梯度值(也即重力异常梯度值)分别为Γ^(t)，Γ^(t+1)，Γ^(t+2)。

步骤2、构建重力梯度模型：

对相邻时刻的重力梯度值求均值，得到如下等式：

再将两者相除

之所以相除就是为了消去质量m，由于之前关于目标探测的最大的问题，就是密度分布不均匀，形状不一定规则，而上面的三重积分则是求质量。通过上面方法建立方程，方程中只含有探测的质心坐标。因重力梯度张量有6个方向，这里面是立体坐标系有三个未知数，必须至少三个方程才方便求解，由于xx,yy,zz方向特征明显，便于求解，故本文选择这三个方向的张量建立非线性方程组。

通过上述分析，我们建立如下方程：

其中，l是相邻测量点的步长。

将上面的等式进行化简处理得到本发明的重力梯度模型：

在本发明的重力梯度模型中，重力梯度的值都可以通过重力梯度仪测得，进而可以得到GGVRxx,GGVRyy,GGVRzz。从而建立非线性方程组。

步骤3、非线性方程组的求解：

设定初始值，利用牛顿迭代法进行求解，具体过程参照图2：

非线性方程组F(x,y,z)＝0，三个方程分别为f₁(x,y,z)＝0，f₂(x,y,z)＝0，

f₃(x,y,z)＝0，x＝(x,y,z)^T∈R³，求解方案如下：

步骤3.1：设定初始值x⁽⁰⁾，以及误差限E；

步骤3.2：由公式x^(k+1)＝x^(k)-F'(x^(k))^-1F(x^(k))，将x⁽⁰⁾代入，可以求得x⁽¹⁾，x⁽²⁾...x⁽ⁿ⁾；

其中F'(x)^-1是F(x)的雅克比矩阵的的逆矩阵，F(x)的雅克比矩阵如下：

步骤3.3：判断，由Δx^(k)＝x^(k+1)-x^(k)与误差限E进行对比，如果|Δx^(k)|＜E，则停止迭代，否则重复步骤3.2；

步骤3.4：如果一直为大于误差限E，则重新设定初始值x⁽⁰⁾，以及误差限E，进行上述的步骤，迭代计算。

步骤4、目标体的质量探测：

目标体的质量m可以通过如下公式计算得到：

上式中G是万有引力常量，R由三个方向的距离平方和开方获取，z方向的相对距离Δz由非线性方程组求得，将各个量带入，即可求得质量。

下面选取一个具体测量实施例来对本发明的方法进行说明：

目标体选择为长74.9米，宽8.4米，高5.3米，密度2.8g/cm³，总质量9336.7334吨的长方体。图3均以目标体质心正上方建立如图1所示的三维坐标系，其中z轴垂直向下，航线从左至右方向进行实验。目标正上方300m为坐标原点，z轴垂直向下，航线方向为x轴正方向，建立三维坐标系。

(1)探测重力梯度：在如图3所示的目标体的活动水域内，潜器在目标体质心上方100米的水平面内作实时测量，首先我们通过重力梯度仪获取一系列的重力梯度数据值，假设潜器的水下速度为10m/s，每隔50米进行一次采样。现采用如下方法进行目标探测：由此得出实测重力梯度信号，并记录数据。图4-6依次是xx,yy,zz方向的张量随x正方向的变化曲线图，图4-6中，横坐标的起点x＝-1000，终点是x＝0，单位m(米)；纵坐标是重力梯度值，单位为E，1E＝1×10^-9/s²＝1mGal/(10km)。

(2)构建重力梯度模型：

通过获得的重力梯度数据，可将每一个测量点(即水下潜器所处位置)测得的由目标引起的重力梯度异常带入下面的公式，通过下面的方程组求解得出目标质心方位的相对位置坐标。

其中x,y,z是三个方向的相对距离，l是相邻测量点的步长，

Γ_xx ^(t)是t时刻xx方向的重力梯度值，Γ_yy ^(t)是t时刻yy方向的重力梯度值，Γ_zz ^(t)是t时刻zz方向的重力梯度值。

(3)方程组求解,计算目标相对于水下潜器的方位，得出目标质心方位的相对位置坐标。

通过求解上述方程，便可以得到x，y，z的值。目标体距离质心的距离就可以由上面求解出来的值然如其中，坐标系定义如下：把海平面当成水平面，置x轴和y轴于此平面内，使z轴垂直向下以目标体正上方300米处为原点。

(4)计算目标体的质量

目标的质量m，在多次测量的情况下，将说测得的重力梯度垂直分量可以通过如下公式计算得到：

上式中G是万有引力常量，R由三个方向的距离平方和开方获取，z方向的相对距离Δz由非线性方程组求得，将各个量带入，即可求得质量。对于目标质量的测量，随着离目标体越近，探测范围的增大，测量点越多，最后得到的目标质量误差会逐渐减少，与真实值越来越近，探测的误差会越小。多次测量求均值也是降低误差的一种方式。

(5)目标体的物质种类大致确定

在获取目标体质量后，还可以获取物质种类的大致信息。对于目标体距离测量点的相对位置相等，仅仅密度不同时，重力梯度值和密度成正比，海水的密度已知，海水造成的重力梯度值也是可以测量的，密度异常体的密度也可以测量，那么由如下关系就可以反推密度异常体的平均密度，进而由密度值来推测出大概为何种物质。

综上所述，对于载有高精度重力梯度仪的水下潜器，通过实施精确测量航迹中的重力梯度值以及运用本发明提出的牛顿迭代法的水中重力梯度目标探测，可以对目标距离测量点的相对位置以及目标体质量进行准确估计。

Claims (7)

1.一种基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于包括如下主要步骤：

S1：利用水下潜器的重力梯度仪对航行的区域进行实时测量，获取实测的重力梯度值；

S2：利用实测的重力梯度值去减去正常的重力梯度值，得到重力梯度异常值，即为目标体或密度异常体引起的重力梯度值；

S3：将步骤S2得到的目标体或密度异常体引起的重力梯度值作为重力梯度模型非线性方程组的输入，求解重力梯度模型非线性方程组获目标体的质心相对于探测点的相对位置；由目标体的质心相对于探测点的相对位置，确定目标体的质心与探测点之间的距离；

步骤S3中，将重力梯度模型非线性方程组化简，选取三个张量求解，得到目标体距离探测点在三维坐标系中各个方向的相对距离；重力梯度模型非线性方程组化简如下：

其中x,y,z是目标体质心与测量点在x,y,z三个方向的相对距离，l是相邻测量点的步长，Γ_xx ^(t)是t时刻xx方向的重力梯度值，Γ_yy ^(t)是t时刻yy方向的重力梯度值，Γ_zz ^(t)是t时刻zz方向的重力梯度值；通过求解上述方程组，便可以得到x,y,z的值，获取目标体与探测点在三个方向的相对距离；

S4：最后将上述距离和万有引力常量代入重力梯度垂直分量公式，获取目标体的质量,完成水下目标探测。

2.根据权利要求1所述的基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于：

在步骤S4之后，还可以进行步骤S5：根据目标体的质量和前述测得的目标体或密度异常体引起的重力梯度值，计算出目标体的平均密度，根据平均密度大致判断目标体物质种类。

3.根据权利要求1所述的基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于：所述目标体为海底山峰、敌方潜艇、礁石、沉船中的至少一种物质，当它们连成一体或重叠时将它们看成一个整体处理。

4.根据权利要求1所述的基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于：步骤S3中，选取xx,yy,zz三个张量，分别是重力位对x方向的二阶导数，y方向的二阶导数，zz方向的二阶导数。

5.根据权利要求1所述的基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于：步骤S3中，根据以下公式获取目标体质心与测量点的距离：

R代表目标体质心与测量点的距离，x,y,z是目标体质心与测量点在x,y,z三个方向的相对距离。

6.根据权利要求5所述的基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于：步骤S3中，采用经典的牛顿迭代法求解非线性方程组。

7.根据权利要求5所述的基于牛顿迭代法的水下重力梯度目标探测方法，其特征在于：步骤S4中，目标体的质量m通过如下重力梯度垂直分量公式计算得到：

其中G为万有引力常量，Γ_zz是由目标体或密度异常体引起的重力梯度值，Δz为目标体z方向到探测点的相对距离。