一种基于重力矢量及其梯度张量的多运动目标体监测方法与流程

本发明涉及一种勘探地球物理领域的多目标体定位方法，特别涉及一种基于重力矢量及其梯度张量的多运动目标体的监测方法。

背景技术：

航空器的飞行轨迹监测、水下潜艇的质量估计及其潜航轨迹探测、弹道导弹末端高精度跟踪是目前民用、军事领域的研究难点。目前，虽然存在一些方法技术来定位运动目标体的轨迹，如声学雷达和声呐、电磁激光波法和磁测技术，但是这些技术易受人文干扰影响，具有不稳定性。而质量或者密度是运动目标体的固有属性，不随电磁和声波环境而改变。为此，人们越来越重视基于重力场或重力梯度数据探测的方法来探测和估计水下或者空中运动目标体的位置和剩余质量。其中，2010年和2013年，武凛等人利用面积性的重力梯度数据提出了一组解析公式来定位单个水下运动目标体；2013年，华中科技大学闫足等人利用重力梯度的差分比值定位单个水下物体，并于2015年提出了一种基于重力矢量和垂直重力梯度的单目标定位解析公式。最近，汤井田等人于2017年提出了一种基于重力矢量和重力梯度张量的单运动目标体定位方法和一种基于旋转重力场的单运动目标体监测方法。但是，现有的技术中，基于重力矢量及其梯度张量只能确定单一目标体位置和质量，无法实现多运动目标体的定位。因此，有必要设计一种进行多个运动目标体定位和跟踪的方法。

技术实现要素：

本发明所解决的技术问题是，针对现有技术的不足，提供一种基于重力矢量及其梯度张量的多运动目标体监测方法，利用重力观测基站观测剩余重力矢量和剩余重力梯度张量，以确定多个运动目标体的位置和质量，并跟踪其运动轨迹。

本发明的技术方案为：

一种基于重力矢量及其梯度张量的多运动目标体监测方法，包括以下步骤：

步骤1、布置重力观测基站：

设空间中存在n个运动目标体，根据运动目标体的方位确定m个重力观测基站的位置，并在m个重力观测基站上安装重力仪和重力梯度仪；其中m≥n；

步骤2、计算重力观测基站上多个运动目标体引起的剩余重力矢量g(t)和剩余重力梯度张量t(t)：

定义空间直角坐标系(x,y,z)；利用重力仪和重力梯度仪按一定时间间隔记录空间直角坐标系(x,y,z)下第i个重力观测基站的重力矢量和重力梯度张量(即运动目标体对重力观测基站不产生影响时观测到的重力矢量及重力梯度张量)，根据公式(1)计算第i个重力观测基站上运动目标体引起的剩余重力矢量gⁱ(t)和剩余重力梯度张量tⁱ(t)：

上式中，t0为时间基点；

步骤3、由gⁱ(t)的垂直分量和消除质量的影响，获得只包含位置信息的第i个重力观测基站对应的非线性方程组，其中i＝1,2,…,m：

其中，

式中，为第j个运动目标体与第i个重力观测基站之间的距离，δx^ij(t)＝x^j(t)-aⁱ,δy^ij(t)＝y^j(t)-bⁱ,δz^ij(t)＝z^j(t)-cⁱ；x^j(t),y^j(t)和z^j(t)为t时刻第j个运动目标体在x、y和z轴上的坐标位置；aⁱ,bⁱ和cⁱ第i个重力观测基站在x、y和z轴上的坐标位置；并且，当n＝2时，gl(t)＝g(t)^-1；当n>2时，gl(t)＝g(t)^t[g(t)g(t)^t]^-1；

步骤4、根据步骤3得到m个非线性性方程组，将它们组合成一个非线性性方程组，求解组合后的非线性性方程组，得到多个运动目标体在t时刻的质心坐标位置(x(t),y(t),z(t))n×3。

进一步地，将所述步骤4得到的每个运动目标体在t时刻的质心坐标位置带入公式(4)，通过线性方程组求解得到第j个运动目标体的剩余质量mj(剩余质量一旦获得便不改变)，j＝1,2,…,n；

进一步地，所述步骤4中，将由单站单目标定位算法所计算的位置作为初始值，利用标准的levenberg-marquardt(中文为列文伯格-马夸尔特法，简称lm算法)算法求解组合后的非线性性方程组，得到多个运动目标体在t时刻的质心坐标位置(x(t),y(t),z(t))n×3。

进一步地，所述步骤4具体步骤为：

将多个运动目标体看成一个整体，通过单站单目标定位算法计算每个重力观测基站上估计的该整体的位置，形成位置信息矩阵从中分别为n个运动目标体自动挑选一组坐标作为其位置初始值，形成初始位置矩阵x0(t)＝(x0(t),y0(t),z0(t))n×3；

调用函数[x,fval,flag]＝fsolve(@fun,x0(t),option)，其中fsolve是matlab软件的一条函数，用于求解非线性方程组，fun为要求解的非线性方程组，即步骤4中组合后的非线性性方程组，x0(t)为求解过程的初始值，即初始位置矩阵，option为优化算法的选项设定，其中优化算法为levenberg-marquardt算法；x＝(x(t),y(t),z(t))n×3为返回的解，fval为将x代入非线性方程组的验证结果，即将x代入非线性方程组左边的数学式，计算其结果，再计算该结果与非线性方程组右边的1的差值，将其绝对值的最大值记为fval；flag为终止迭代的标志，若fval小于给定的阈值，则flag＝1，将此时的x＝(x(t),y(t),z(t))n×3作为求解得到的运动目标体在t时刻的质心坐标位置，否则flag＝0，从中重新挑选坐标作为位置初始值，形成初始位置矩阵x0(t)，重复上述步骤。

依次求解各个运动目标体在每个时刻的质心坐标位置，得到其运动轨迹。

进一步地，所述步骤4中，单站单目标定位算法为：

由剩余重力矢量gⁱ(t)的水平分量和垂直分量以及剩余重力梯度张量tⁱ(t)的水平分量和将多个运动目标体看成一个整体，根据公式(5)计算第i个重力观测基站上估计的该整体的位置：

进一步地，所述运动目标体位于太空、地表或水面以下、地表或水面以上，或位于地表面或水面。

进一步地，所述重力观测基站固定在地表面或地表以上某个位置，或位于船只、飞机或卫星等移动平台上。

本发明原理为：

(1)单运动目标体定位算法

假设运动目标体的体积为λ，剩余密度为ρ，第i个重力观测基站p(aⁱ,bⁱ,cⁱ)上的剩余重力位函数φⁱ满足下列积分方程：

式中为运动目标体与重力观测基站p(aⁱ,bⁱ,cⁱ)之间的距离，δxⁱ＝x-aⁱ，δyⁱ＝y-bⁱ，δzⁱ＝z-cⁱ，x,y,z为运动目标体的位置，g＝6.67×10^-11kg^-1m³s^-2为引力常数。

重力矢量gⁱ及其梯度张量tⁱ可以写成如下形式：

如果重力观测基站与运动目标体之间的距离足够大，运动目标体可以被看成一个质点，此时，公式(7)和公式(8)可以写成如下形式：

使用公式(9)的重力矢量和公式(10)中的水平重力梯度张量，则在第i个重力观测基站上所计算的运动目标体的质心坐标和剩余质量可从公式(11)获得。

(2)多运动目标体定位算法

假设空间中存在n个运动目标体，其剩余密度为ρj(j＝1,2,…,n)。根据公式(9)和(10)，可以得到第i个重力观测基站上的剩余重力矢量gⁱ和重力梯度张量tⁱ：

式中，为运动目标体o(x^j,y^j,z^j)与重力观测基站p(aⁱ,bⁱ,cⁱ)之间的距离，δx^ij＝x^j-aⁱ，δy^ij＝y^j-bⁱ，δz^ij＝z^j-cⁱ；x^j,y^j,z^j和mj为运动目标体o(x^j,y^j,z^j)的位置和质量。

此时，可以被写成如下形式：

式中aⁱ＝[1…1](1×n)，

结合公式(13)和(14)，用来表示和

设公式(15)和(16)可以改写为：

公式(17)依然受到运动目标体质量mj的影响，为了消除质量的影响，考虑公式(14)，将公式(17)除以便可得到归一化非线性方程：

当n＝2时，gl＝g^-1；当n>2时，gl＝g^t(gg^t)^-1。

同理，将公式(17)分别除以和可得：

其中，

本发明利用标准的levenberg-maquardt算法求解公式(18)和公式(19)所组成的非线性方程组。众所周知，对于这种迭代的方法，初始值的选择直接影响了计算结果的可靠性。本发明给出了一种数值稳定且高精度的单站单目标定位算法，并在单站单目标定位算法的基础上给出了一个挑选位置初始值的可行的方案：将多个运动目标体看成一个整体，通过单站单目标定位算法计算每个重力观测基站上估计的该整体的位置，形成位置信息矩阵从中分别为n个运动目标体自动挑选一组坐标作为其位置初始值，形成初始位置矩阵x0(t)＝(x0(t),y0(t),z0(t))n×3；进而基于levenberg-marquardt(中文为列文伯格-马夸尔特法，简称lm算法)算法求解非线性性方程组，得到的运动目标体在t时刻的质心坐标位置：令fun为公式(18)和公式(19)所组成的非线性方程组，调用函数[x,fval,flag]＝fsolve(@fun,x0(t),option)，当fval小于给定的阈值时，将此时的x＝(x(t),y(t),z(t))n×3作为求解得到的运动目标体在t时刻的质心坐标位置。

将得到的每个运动目标体在t时刻的质心坐标位置带入公式(12)和公式(13)，求解线性方程组，便可得到运动目标体的剩余质量mj。

有益效果：

本方法给出了一种利用重力观测基站的重力矢量及其梯度张量监测多个运动目标体的坐标位置(运动轨迹)和质量的算法，可以实时地监测运动目标体的位置和质量，精度高，为下一步军用或者民用成果转化奠定了基础。通过观测剩余重力矢量和剩余重力梯度张量随时间的变化，即可求出多个运动目标体的运动轨迹。本发明中运动目标体可以位于外太空、地表面及以下，如地下洞穴及其中运动的汽车、人员，常规潜艇及核潜艇，航空母舰及各种水面舰艇等；也可以位于空气中，如飞机、导弹等，适用范围广。

附图说明

图1为本发明利用多个重力观测基站进行多目标体运动轨迹识别的示意图，其中oi表示第i个观测点，sj表示第j个目标体，rij表示其间距离；

图2为无噪时双运动目标体双站定位结果图，其中实线为运动目标体a和b的实际航行轨迹，空心三角形为运动目标体a的定位结果，空心圆为运动目标体b的定位结果，空心五角星为重力观测基站。图2(a)分别为运动目标体x-y轴定位结果，图2(b)为运动目标体x-z轴定位结果。

图3为加入方差为0，标准差为0.03e的高斯噪声双运动目标双站定位结果图，其中实线为运动目标体a和b的实际航行轨迹，空心三角形为运动目标体a的定位结果，空心圆为运动目标体b的定位结果，空心五角星为重力观测基站。图3(a)分别为运动目标体x-y轴定位结果，图3(b)为运动目标体x-z轴定位结果。

具体实施方式

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明公开了一种基于重力矢量及其梯度张量的多运动目标体监测方法，包括以下步骤：

(1)重力观测基站布置设计：假设空间中存在n个运动目标体，根据观测目标确定观测基站的位置，并在m个重力观测基站上安装重力仪和重力梯度仪；

(2)计算重力观测基站上水下或空中多个运动目标体引起的剩余重力矢量g(t)和剩余重力梯度张量t(t)；

(3)由每一个重力观测基站上测量的剩余重力矢量g(t)的水平分量gx(t)、gy(t)和垂直分量gz(t)，以及剩余重力梯度张量t(t)的水平分量txx(t)、tyy(t)和txy(t)，将空中或者水下多个目标体看成一个整体，根据公式(5)初步计算每一个重力观测基站上估计的多个运动目标体整体的位置；

(4)实时坐标与质量的计算与监测：整理基站采集的数据，以步骤(3)估计的运动目标体的整体的位置为初始值，利用公式(6)～(20)所示的计算方法，实时确定多运动目标体的坐标和质量等。

实施例1：

本实施例为进行水下多运动目标体轨迹定位的实例。假设重力梯度仪的观测精度为0.01e，运动目标体a和b的质量分别为4.32×10⁶kg和4.61×10⁶kg，运动速度为-2.315千米每小时，目标体a的运动轨迹为方程：x＝115.75-(i-1)×δt×2315，y＝9.0625×10^-1×x²+30，z＝10；目标体b的运动轨迹为方程：x＝115.75-(i-1)×δt×2315，y＝0.75×x，z＝50，21个时间间隔，时间间隔δt为7.5秒，i＝1,2,k,21。图2展示了运动目标体的运动轨迹与重力观测基站的位置。运动目标体a和b起始点分别为(115.75m,42.56m,10m)和(115.75m,86.81m,50m)，重力观测基站位置分别为(50m,90m,0m)和(50m,10m,0m)处(用符号“☆”标示)。从图2可以看出，本发明算法估计的运动目标体的轨迹与实际位置完全吻合，所识别a和b的x、y、z的平均相对误差分别为：0.1％、0.4％、0.4％和1.4％、1.2％、0.6％。其中相对误差r＝|(估计值-真实值)/真实值|×100。说明我们的算法能够精准的实时监测双目标体的运动轨迹，适用性强。

为了检验本发明的稳定性和可靠性，有必要考虑噪声影响。在上述模型中加入方差为0，标准差为0.03e(3倍仪器观测精度)的高斯噪声，其计算结果见图3。如图3所示，加入方差为0，标准差为0.03e的高斯噪声，本发明算法同样可以定位运动目标体的轨迹，所识别a和b的x、y、z的平均相对误差分别为：0.7％、1.5％、4.7％和2.3％、2.5％、3.2％。上述分析说明该方法稳定性强，具有一定的抗噪能力。

分析表明，采用本发明所提供的基于重力矢量及其梯度张量的多运动目标体监测方法，可以实时地监测具有复杂运动轨迹的多个运动目标体的位置，不仅具有较高的民用价值，也具有较大的国防应用价值。