CN106168682B - 一种基于旋转重力场的运动目标体监测方法 - Google Patents

Abstract

本发明提供了一种基于旋转重力场的运动目标体监测方法，包括以下步骤：(1)在一个固定或者移动的重力基站上，按一定时间间隔，观测由质量为M的单个运动物体引起的剩余重力矢量及其梯度张量，得到每个时刻的剩余重力矢量g及其梯度张量T；(2)将所观测到的剩余重力矢量g及其梯度T按照zyz型欧拉角旋转，得到旋转之后的数据集g_z′(t)和T′(t)；(3)根据旋转之后的数据，利用解析公式，精确计算单个运动目标体的坐标位置、质量和速度矢量。本发明可以实时地监测运动目标体的轨迹和质量，有效提高了定位精度，并且不受基站位置的影响，不仅具有较高的民用价值，也具有较大的国防应用价值。

Description

一种基于旋转重力场的运动目标体监测方法

技术领域

本发明涉及一种勘探地球物理领域的目标体定位方法，特别涉及一种基于旋转重力数据的运动目标体的监测方法。

背景技术

在地球重力学中，重力梯度张量除用于勘探地球地质结构构造，寻找油气、矿产资源外，还有其他多种应用。例如，利用卫星重力梯度，高精度恢复中短波的地球重力场；利用高精度重力梯度张量，为舰艇或巡航导弹进行导航；2010年，Majid Beiki and LaustB.Pedersen(Majid Beiki and Laust B.Pedersen，Eigenvector analysis of gravitygradient tensor to locate geologic bodies，Geophysics,2010)应用重力梯度张量特征值与特征向量、垂直方向的重力场计算了地下单一异常体的位置与质量，并被用于重力梯度张量剖面数据的定性解释。2014年，Lockerbie(N A Lockerbie，The location ofsubterranean voids using tensor gravity gradiometry，Class.Quantum Grav.2014)基于2个不同测点处的重力梯度张量特征向量唯一确定了单一目标体的位置，并应用于地下掩体范围的大致圈定。

现有技术中，必须基于2个不同测点处的重力梯度张量确定单一目标体位置测量，且只能应用于地下掩体范围的大致圈定。因此，有必要设计一种基于一个观测基站的旋转重力矢量和重力梯度张量进行单个运动物体定位和跟踪的解析方法。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供了一种基于旋转重力场的运动目标体监测方法，本发明抛开了传统面积性数据采集方式，给出了利用一个观测站所测得的重力矢量和重力梯度张量进行旋转监测单个运动目标体的坐标位置和质量的核心算法，可以实时地监测运动目标体的轨迹和质量。

本发明的技术方案为：

一种基于旋转重力场的运动目标体监测方法，包括以下步骤：

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

根据观测目标确定重力观测基站的位置，并在重力观测基站安装重力仪和重力梯度仪；

步骤2、计算空间直角坐标系(x,y,z)【笛卡尔直角坐标系】下单个运动物体引起的剩余重力矢量g(t)和剩余重力梯度张量T(t)：

定义空间直角坐标系(x,y,z)；利用重力仪和重力梯度仪按一定时间间隔记录空间直角坐标系(x,y,z)下重力观测基站的重力矢量g_total(t)和重力梯度张量T_total(t)，根据公式(1)计算质量为M的单个运动物体引起的剩余重力矢量g(t)和剩余重力梯度张量T(t)：

上式中，T_total(t₀)和g_total(t₀)分别为运动目标体对重力基站不产生影响时观测到的重力矢量及重力梯度张量，t₀称为时间基点；

步骤3、旋转重力场，计算局部旋转坐标系(x′,y′,z′)下的剩余重力矢量和剩余重力梯度张量：

将计算得到的剩余重力矢量g(t)和剩余重力梯度张量T(t)按照zyz型欧拉角旋转，令剩余重力矢量g(t)的方向为旋转后的局部旋转坐标系为(x′,y′,z′)，其中z′的方向为旋转空间直角坐标系使得T_x′y′(t)的绝对值最小，得到旋转后的局部旋转坐标系中的g_z′(t)和T′(t)：

此时局部旋转坐标系中剩余重力矢量的垂直分量g_z′(t)与剩余重力矢量g(t)的模g(t)相等，即g_z′(t)＝g(t)，且g_x′(t)＝g_y′(t)＝0，其中g_x(t)、g_y(t)、g_z(t)为剩余重力矢量g(t)的三个分量；

步骤4、由局部旋转坐标系中剩余重力矢量的垂直分量g_z′(t)及剩余重力梯度张量T′(t)各分量的值T_x′x′(t)，T_y′y′(t)，T_z′z′(t)，根据公式(3)、(4)和(5)计算t时刻在局部旋转坐标系下重力观测基站与运动目标体质心的距离，即矢径d_z′(t)：

当T_x′x′(t)＝T_y′y′(t)时，

当T_x′x′(t)＞T_y′y′(t)时，

当T_x′x′(t)<T_y′y′(t)时，

步骤5、由t时刻在局部旋转坐标系下重力观测基站与运动目标体质心的距离d_z′(t)，经过坐标转换，得到t时刻运动目标体在空间直角坐标系(x,y,z)【笛卡尔直角坐标系】下任意时刻的质心坐标位置：

上式中，x(t),y(t),z(t)为t时刻重力观测基站的坐标位置，G为引力常数，x_c(t),y_c(t),z_c(t)为t时刻运动目标体的质心坐标位置。

运动目标体t时刻的质量M(t)为：

所述步骤3中，重力场旋转的准则是基于如公式(8)和(9)所示的zyz型欧拉角的旋转方式：

(x′,y′,z′)^T＝R(x,y,z)^T (8)

旋转后的局部旋转坐标系中的T′(t)计算公式为：

T′(t)＝RT(t)R^T (10)

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

由运动目标体的质心坐标位置(运动轨迹)，根据公式(11)计算运动目标体相对于重力观测基站的运动速度V＝(V_x,V_y,V_z)：

其中，t_i+1为第i+1个时刻，t_i为第i个时刻；

假设重力基站的移动速度为V_s，则运动目标体的绝对速度为：

V_绝对速度＝V_s+V (12)。

本发明原理为：

假设目标体的体积为V，剩余密度为ρ，笛卡尔坐标系(垂直向下为z轴正方向)，观察点(x,y,z)的重力位满足的积分方程为：

剩余重力矢量g满足下列积分方程：

剩余重力梯度张量T满足下列积分方程：

其中T_xx+T_yy+T_zz＝0,(x_c,y_c,z_c)∈V,r＝|(x,y,z)-(x_c,y_c,z_c)|，V为目标体的体积。

当观测点与运动目标体距离较远时，目标体便可以被看出为一个质点，则公式(14)和(15)可以写成如下形式：

其中，x,y,z为观测点位置，G为引力常数，x_c,y_c,z_c为目标体重心位置，M为目标体质量。

根据公式(16)和(17)之间的相互关系，可以得到以下三种定位方法，为了便于说明将此类定位方法称之为普通重力场定位方法：

当g_z≠0，可以采用公式(18)所示的普通重力场定位方法1：

当g_x≠0，可以采用公式(19)所示的普通重力场定位方法2：

当g_y≠0，可以采用公式(20)所示的普通重力场定位方法3：

上述定位方法是基于运动球体或质点这一前提条件提出来的，在运动物体相对测点不能看成为质点的情况下，会导致识别出错误的结果；另外，公式(18)～(20)中存在多个分母项，当这些分母项符号发生变化时便会增加定位的计算误差。为此，在上述定位方法的基础上提出了局部旋转重力场定位方法。下面将详细讨论局部旋转重力场定位算法的原理和应用效果。

将笛卡尔直角坐标系进行旋转。(x,y,z)为使用者定义的笛卡尔坐标系，(x′,y′,z′)为旋转后的坐标系。假设两个坐标系有相同的坐标原点，可以互相旋转，坐标系的轴或者向量的分量可以通过一个3D的旋转矩阵联系起来：

(x′,y′,z′)^T＝R(x,y,z)^T (21)

这里R为旋转矩阵，通常用欧拉角(φ,θ,ψ)来表示。采用地球物理中常用的zyz型欧拉角，首先关于z轴旋转φ，然后关于y′旋转角度θ，最后关于z′旋转ψ。这个旋转矩阵可以用下式表示：

R是一个正交矩阵，它的逆矩阵就是它的转置矩阵，(x,y,z)^T＝R^T(x′,y′,z′)^T

将旋转矩阵应用到任意向量或者向量操作，则有：

将计算得到的剩余重力矢量g和剩余重力梯度张量T按照zyz型欧拉角旋转，令剩余重力矢量g的方向为旋转后的局部旋转坐标系为(x′,y′,z′)，其中z′的方向为旋转空间直角坐标系使得T_x′y′的绝对值最小，得到旋转后的局部旋转坐标系中的g_z′和T′：

此时局部旋转坐标系中剩余重力矢量的垂直分量g_z′与剩余重力矢量g的模g相等，即g_z′＝g，且g_x′＝g_y′＝0，其中g_x、g_y、g_z为剩余重力矢量g的三个分量。

在旋转之后的局部参考系中，结合重力场一阶导数和重力梯度之间的相互关系，利用g_x′,g_y′,g_z′,T_z′z′和式(18)可以得到第一组目标体坐标位置与质量的计算公式：

同理，利用g_x′,g_y′,g_z′,T_x′x′和式(19)可以得到第二组目标体坐标位置与质量的计算公式：

同理，利用g_x′,g_y′,g_z′,T_y′y′和式(20)可以得到第三组目标体坐标位置与质量的计算公式：

经过比较，得到一套比较理想的方案，即：

当T_x′x′＝T_y′y′时，

当T_x′x′>T_y′y′时，

当T_x′x′<T_y′y′时，

经过坐标转换便可以得到目标体在笛卡尔直角坐标系下的位置：

由物体质心轨迹，根据公式(32)计算物体相对于重力基站的运动速度V＝(V_x,V_y,V_z)：

其中，t_i+1为第i+1个时刻，t_i为第i个时刻。

假设重力基站的移动速度为V_s，运动物体的绝对速度为：

V_绝对速度＝V_s+V (33)

有益效果：

本发明抛开了传统面积性数据采集方式，给出了利用一个观测站所测得的重力矢量和重力梯度张量进行旋转监测单个运动目标体的坐标位置和质量的核心算法，可以实时地监测运动目标体的位置和质量。通过观测重力梯度张量和重力矢量随时间的变化，即可求出物体的运动轨迹。

本发明解决了普通定位方法受质点或球体假设的限制，以及重力场或梯度在符号变化的邻域内数值离散的问题，有效提高了定位精度；并且运动目标体可以位于地表面及以下，如地下洞穴及其中运动的汽车、人员，常规潜艇及核潜艇，航空母舰及各种水面舰艇等；也可以位于空气中，如飞机、导弹，特别是隐身飞机。由于质量是物体的固有属性，地球重力场及其梯度不受电磁波、声波、热等因素的影响，因此，目前的隐身方法相对于重力梯度是通明的。依据本发明方法，在理论上，只要布置合适的观测台网(固定台网+移动台站)，在适当的先验信息下，就可以监测特定区域一些特定目标的运动状态，并根据其质量，评估其载弹量及其型号。配合其他监测手段，可以获得更详细的目标信息。因此，本发明方法不仅具有较高的民用价值，也具有较大的国防应用价值。

附图说明：

图1为本发明利用一个重力观测基站进行单体运动轨迹识别的示意图。

图2为单运动目标体单站定位结果图，其中实线为运动目标体实际航行轨迹，空心圆为常规普通重力场定位算法所识别的结果，三角形为本发明局部旋转重力场定位算法所识别的结果；图2(a)-(c)分别为运动目标体x、y、z轴定位结果，图2(d)为运动目标体质量识别结果。

图3为加入5％高斯噪声运动目标体单站定位结果图，其中图3(a)-(c)中实线为运动目标体实际航行轨迹，三角形为本发明局部旋转重力场定位算法所识别的结果；图3(d)为相对误差对比图，其相对误差计算公式为：相对误差＝|(估计值-真实值)/真实值|×100。

具体实施方式：

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

本发明涉及的定位方法包括以下步骤：

步骤1、基站布置设计：确定观测目标及仪器灵敏度范围，合理设计观测台站的位置；在固定的台站或者移动的基站上安装高精度的重力仪和重力梯度仪。

步骤2、异常重力场和重力梯度张量的计算：利用重力仪和重力梯度仪实时记录基站的重力场与重力梯度张量，以便实时计算出重力场与重力梯度张量的异常变化，即剩余重力矢量值g和重力梯度值T；将所计算的剩余重力矢量值g和重力梯度值T按照zyz型欧拉角旋转，得到旋转之后的数据集g_z′(t)和T′(t)。

步骤3、实时坐标与剩余重量计算与监测：整理基站采集数据，利用公式(1)～(12)所示的计算方法，实时确定单一移动目标体的坐标、质量和运动速度等。

以下为本发明水下移动单体轨迹定位的实例。

设单一棱柱体的质量为1.6×10⁷kg，长170m，宽12.8m，高11m。

运动轨迹为的抛物型方程，y＝0.1875x²+100，z＝200，Δt＝1,2,...,24。运动单体起始点x坐标为-1000米，重力观察基站位于(0m，-100m，0m)处，共记录24个时刻的重力梯度张量和垂直重力场。从图2可以看出，在距离测点较远的区域，普通重力场定位方法可以准确地识别棱柱体的轨迹和质量，但是当物体运动到(0m，100m，200m)附近区域时，追踪识别的结果明显偏离了物体的轨迹，这主要有以下两个方面的原因：(1)在观测点附近目标体不能看成为质点；(2)公式中的分母项在该区域符号发生了变化，由于数值离散不恰当增加计算误差。而本发明局部旋转重力场定位方法解决了上述两个问题，识别棱柱体的轨迹和质量吻合了真实运动轨迹和质量。

为了检验本发明的稳定性和可靠性，有必要考虑噪声影响。棱柱体运动轨迹为如图3实线，目标体起始点x坐标为-1000m，观测间隔为1/8分钟，共记录12个时刻的重力矢量及其梯度。观测点坐标为(0m，-100m，0m)，在所测量的重力矢量及其梯度张量中分别加入5％的高斯随机噪声。其计算结果见图3和表1。如图3和表1所示，加入5％的噪声，本发明算法同样可以定位运动目标体的轨迹和质量，所识别的x、y、z和质量M的平均相对误差分别为：3.32％、2.58％、1.85％和3.48％，最大相对误差分别为：6.26％、6.37％、5.97％和7.51％，其中绝对误差D＝估计结果-真实值，相对误差R＝|(估计值-真实值)/真实值|×100。其中x方向的最大绝对偏差为43.55m，小于85m，没有超出物体的边界；y和z方向的最大绝对偏差分别为10.37m和8.87m，个别点偏离了物体的边界，但偏差较小，均在允许的误差范围以内。当运动目标体距离测量点较远时，偏差有所增大，这是因为此时重力梯度张量各分量的值较小，易受噪声的影响。上述分析说明该方法稳定性强，具有一定的抗噪能力。

表1加入5％高斯噪声运动目标体单站定位结果误差表

分析表明，采用本发明所提供的基于局部旋转重力场的运动目标体监测方法，可以实时地监测具有复杂运动轨迹的目标体的位置和质量，比在直角坐标系下目标识别算法更具有适用性，不仅具有较高的民用价值，也具有较大的国防应用价值。

Claims (5)

1.一种基于旋转重力场的运动目标体监测方法，其特征在于，包括以下步骤：

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

根据运动目标体确定重力观测基站的位置，并在重力观测基站安装重力仪和重力梯度仪；

步骤2、计算空间直角坐标系(x,y,z)下运动目标体引起的剩余重力矢量g(t)和剩余重力梯度张量T(t)：

定义空间直角坐标系(x,y,z)；利用重力仪和重力梯度仪按一定时间间隔记录空间直角坐标系(x,y,z)下重力观测基站的重力矢量g_total(t)和重力梯度张量T_total(t)，根据公式(1)计算质量为M的运动目标体引起的剩余重力矢量g(t)和剩余重力梯度张量T(t)：

$\begin{array}{c}T\left(t\right)=T_{\mathrm{total}}\left(t\right)-T_{\mathrm{total}}\left(t_0\right)\\ g\left(t\right)=g_{\mathrm{total}}\left(t\right)-g_{\mathrm{total}}\left(t_0\right)\end{array}---\left(1\right)$

上式中，T_total(t₀)和g_total(t₀)分别为运动目标体对重力观测基站不产生影响时观测到的重力矢量及重力梯度张量，t₀称为时间基点；

步骤3、旋转重力场，计算局部旋转坐标系(x′,y′,z′)下的剩余重力矢量和剩余重力梯度张量：

将计算得到的剩余重力矢量g(t)和剩余重力梯度张量T(t)按照zyz型欧拉角旋转，令剩余重力矢量g(t)的方向为旋转后的局部旋转坐标系为(x′,y′,z′)，其中z′的方向为旋转空间直角坐标系使得T_x′y′(t)的绝对值最小，得到旋转后的局部旋转坐标系中的g_z′(t)和T′(t)：

$T^{\&prime;}\left(t\right)=\left[\begin{array}{ccc}{T}_{x^{\&prime;}{x}^{\&prime;}}\left(t\right)& 0& T_{x^{\&prime;}{z}^{\&prime;}}\left(t\right)\\ 0& T_{y^{\&prime;}{y}^{\&prime;}}\left(t\right)& T_{y^{\&prime;}{z}^{\&prime;}}\left(t\right)\\ T_{z^{\&prime;}{x}^{\&prime;}}\left(t\right)& T_{z^{\&prime;}{y}^{\&prime;}}\left(t\right)& T_{z^{\&prime;}{z}^{\&prime;}}\left(t\right)\end{array}\right]---\left(2\right)$

此时局部旋转坐标系中剩余重力矢量的垂直分量g_z′(t)与剩余重力矢量g(t)的模g(t)相等，即g_z′(t)＝g(t)，且g_x′(t)＝g_y′(t)＝0，其中g_x(t)，g_y(t)，g_z(t)为剩余重力矢量g(t)的三个分量；

步骤4、由局部旋转坐标系中剩余重力矢量的垂直分量g_z′(t)及剩余重力梯度张量T′(t)各分量的值T_x′x′(t)，T_y′y′(t)，T_z′z′(t)，根据公式(3)、(4)和(5)计算t时刻在局部旋转坐标系下重力观测基站与运动目标体质心的距离，即矢径d_z′(t)：

当T_x′x′(t)＝T_y′y′(t)时，

$d_{z^{\&prime;}}\left(t\right)=-\frac{2g_{z^{\&prime;}}\left(t\right)}{T_{z^{\&prime;}{z}^{\&prime;}}\left(t\right)}---\left(3\right)$

当T_x′x′(t)＞T_y′y′(t)时，

$d_{z^{\&prime;}}\left(t\right)=\frac{g_{z^{\&prime;}}\left(t\right)}{T_{x^{\&prime;}{x}^{\&prime;}}\left(t\right)}---\left(4\right)$

当T_x′x′(t)<T_y′y′(t)时，

$d_{z^{\&prime;}}\left(t\right)=\frac{g_{z^{\&prime;}}\left(t\right)}{T_{y^{\&prime;}{y}^{\&prime;}}\left(t\right)}---\left(5\right)$

步骤5、由t时刻在局部旋转坐标系下重力观测基站与运动目标体质心的距离d_z′(t)，经过坐标转换，得到t时刻运动目标体在空间直角坐标系(x,y,z)下任意时刻的质心坐标位置：

${\{x_c\left(t\right),y_c\left(t\right),z_c\left(t\right)\}}^T={\{x\left(t\right),y\left(t\right),z\left(t\right)\}}^T-{\hat{z}}^{\&prime;}\&CenterDot;d_{z^{\&prime;}}\left(t\right)---\left(6\right)$

上式中，{x(t),y(t),z(t)}为t时刻重力观测基站的坐标位置，{x_c(t),y_c(t),z_c(t)}为t时刻运动目标体的质心坐标位置。

2.根据权利要求1所述的基于旋转重力场的运动目标体监测方法，其特征在于：运动目标体t时刻的质量M(t)为：

$M\left(t\right)=|d_{z^{\&prime;}}^2\left(t\right)\&CenterDot;g_{z^{\&prime;}}\left(t\right)\&CenterDot;\left(\frac{1}{-G}\right)|---\left(7\right)$

上式中，G为引力常数。

3.根据权利要求1或2所述的基于旋转重力场的运动目标体监测方法，其特征在于：所述步骤3中，重力场旋转的准则是基于如公式(8)和(9)所示的zyz型欧拉角的旋转方式：

(x′,y′,z′)^T＝R(x,y,z)^T (8)

$R=\left(\begin{array}{ccc}\cos \mathrm{\&psi;}\cos \mathrm{\&theta;}\cos \mathrm{\&phi;}& \cos \mathrm{\&psi;}\cos \mathrm{\&theta;}\sin \mathrm{\&phi;}+\sin \mathrm{\&psi;}\cos \mathrm{\&phi;}& \cos \mathrm{\&psi;}\sin \mathrm{\&theta;}\\ -\sin \mathrm{\&psi;}\cos \mathrm{\&theta;}\cos \mathrm{\&phi;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}& -\sin \mathrm{\&psi;}\cos \mathrm{\&theta;}\sin \mathrm{\&phi;}+\cos \mathrm{\&psi;}\cos \mathrm{\&phi;}& -\sin \mathrm{\&psi;}\sin \mathrm{\&phi;}\\ -\sin \mathrm{\&theta;}\cos \mathrm{\&phi;}& -\sin \mathrm{\&theta;}\sin \mathrm{\&phi;}& \cos \mathrm{\&theta;}\end{array}\right)---\left(9\right)$

旋转后的局部旋转坐标系中的T′(t)计算公式为：

T′(t)＝RT(t)R^T (10)。

4.根据权利要求1所述的基于旋转重力场的运动目标体监测方法，其特征在于：重力观测基站固定在地表面或地表以上某个位置，或位于移动平台上。

5.根据权利要求1所述的基于旋转重力场的运动目标体监测方法，其特征在于：由运动目标体的质心坐标位置，根据公式(11)计算运动目标体相对于重力观测基站的运动速度V＝(V_x,V_y,V_z)：

$V_x=\frac{x_c\left(t_{i+1}\right)-x_c\left(t_i\right)}{t_{i+1}-t_i}$

$V_y=\frac{y_c\left(t_{i+1}\right)-y_c\left(t_i\right)}{t_{i+1}-t_i}---\left(11\right)$

$V_z=\frac{z_c\left(t_{i+1}\right)-z_c\left(t_i\right)}{t_{i+1}-t_i}$

其中，t_i+1为第i+1个时刻，t_i为第i个时刻；

假设重力基站的移动速度为V_s，则运动目标体的绝对速度为：

V_绝对速度＝V_s+V (12)。