本發(fā)明涉及一種勘探地球物理領(lǐng)域的多目標體定位方法，特別涉及一種基于重力矢量及其梯度張量的多運動目標體的監(jiān)測方法。

背景技術(shù)：

航空器的飛行軌跡監(jiān)測、水下潛艇的質(zhì)量估計及其潛航軌跡探測、彈道導(dǎo)彈末端高精度跟蹤是目前民用、軍事領(lǐng)域的研究難點。目前，雖然存在一些方法技術(shù)來定位運動目標體的軌跡，如聲學(xué)雷達和聲吶、電磁激光波法和磁測技術(shù)，但是這些技術(shù)易受人文干擾影響，具有不穩(wěn)定性。而質(zhì)量或者密度是運動目標體的固有屬性，不隨電磁和聲波環(huán)境而改變。為此，人們越來越重視基于重力場或重力梯度數(shù)據(jù)探測的方法來探測和估計水下或者空中運動目標體的位置和剩余質(zhì)量。其中，2010年和2013年，武凜等人利用面積性的重力梯度數(shù)據(jù)提出了一組解析公式來定位單個水下運動目標體；2013年，華中科技大學(xué)閆足等人利用重力梯度的差分比值定位單個水下物體，并于2015年提出了一種基于重力矢量和垂直重力梯度的單目標定位解析公式。最近，湯井田等人于2017年提出了一種基于重力矢量和重力梯度張量的單運動目標體定位方法和一種基于旋轉(zhuǎn)重力場的單運動目標體監(jiān)測方法。但是，現(xiàn)有的技術(shù)中，基于重力矢量及其梯度張量只能確定單一目標體位置和質(zhì)量，無法實現(xiàn)多運動目標體的定位。因此，有必要設(shè)計一種進行多個運動目標體定位和跟蹤的方法。

技術(shù)實現(xiàn)要素：

本發(fā)明所解決的技術(shù)問題是，針對現(xiàn)有技術(shù)的不足，提供一種基于重力矢量及其梯度張量的多運動目標體監(jiān)測方法，利用重力觀測基站觀測剩余重力矢量和剩余重力梯度張量，以確定多個運動目標體的位置和質(zhì)量，并跟蹤其運動軌跡。

本發(fā)明的技術(shù)方案為：

一種基于重力矢量及其梯度張量的多運動目標體監(jiān)測方法，包括以下步驟：

步驟1、布置重力觀測基站：

設(shè)空間中存在n個運動目標體，根據(jù)運動目標體的方位確定m個重力觀測基站的位置，并在m個重力觀測基站上安裝重力儀和重力梯度儀；其中m≥n；

步驟2、計算重力觀測基站上多個運動目標體引起的剩余重力矢量g(t)和剩余重力梯度張量t(t)：

定義空間直角坐標系(x,y,z)；利用重力儀和重力梯度儀按一定時間間隔記錄空間直角坐標系(x,y,z)下第i個重力觀測基站的重力矢量和重力梯度張量(即運動目標體對重力觀測基站不產(chǎn)生影響時觀測到的重力矢量及重力梯度張量)，根據(jù)公式(1)計算第i個重力觀測基站上運動目標體引起的剩余重力矢量gⁱ(t)和剩余重力梯度張量tⁱ(t)：

上式中，t0為時間基點；

步驟3、由gⁱ(t)的垂直分量和消除質(zhì)量的影響，獲得只包含位置信息的第i個重力觀測基站對應(yīng)的非線性方程組，其中i＝1,2,…,m：

其中，

式中，為第j個運動目標體與第i個重力觀測基站之間的距離，δx^ij(t)＝x^j(t)-aⁱ,δy^ij(t)＝y(tǒng)^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、根據(jù)步驟3得到m個非線性性方程組，將它們組合成一個非線性性方程組，求解組合后的非線性性方程組，得到多個運動目標體在t時刻的質(zhì)心坐標位置(x(t),y(t),z(t))n×3。

進一步地，將所述步驟4得到的每個運動目標體在t時刻的質(zhì)心坐標位置帶入公式(4)，通過線性方程組求解得到第j個運動目標體的剩余質(zhì)量mj(剩余質(zhì)量一旦獲得便不改變)，j＝1,2,…,n；

進一步地，所述步驟4中，將由單站單目標定位算法所計算的位置作為初始值，利用標準的levenberg-marquardt(中文為列文伯格-馬夸爾特法，簡稱lm算法)算法求解組合后的非線性性方程組，得到多個運動目標體在t時刻的質(zhì)心坐標位置(x(t),y(t),z(t))n×3。

進一步地，所述步驟4具體步驟為：

將多個運動目標體看成一個整體，通過單站單目標定位算法計算每個重力觀測基站上估計的該整體的位置，形成位置信息矩陣從中分別為n個運動目標體自動挑選一組坐標作為其位置初始值，形成初始位置矩陣x0(t)＝(x0(t),y0(t),z0(t))n×3；

調(diào)用函數(shù)[x,fval,flag]＝fsolve(@fun,x0(t),option)，其中fsolve是matlab軟件的一條函數(shù)，用于求解非線性方程組，fun為要求解的非線性方程組，即步驟4中組合后的非線性性方程組，x0(t)為求解過程的初始值，即初始位置矩陣，option為優(yōu)化算法的選項設(shè)定，其中優(yōu)化算法為levenberg-marquardt算法；x＝(x(t),y(t),z(t))n×3為返回的解，fval為將x代入非線性方程組的驗證結(jié)果，即將x代入非線性方程組左邊的數(shù)學(xué)式，計算其結(jié)果，再計算該結(jié)果與非線性方程組右邊的1的差值，將其絕對值的最大值記為fval；flag為終止迭代的標志，若fval小于給定的閾值，則flag＝1，將此時的x＝(x(t),y(t),z(t))n×3作為求解得到的運動目標體在t時刻的質(zhì)心坐標位置，否則flag＝0，從中重新挑選坐標作為位置初始值，形成初始位置矩陣x0(t)，重復(fù)上述步驟。

依次求解各個運動目標體在每個時刻的質(zhì)心坐標位置，得到其運動軌跡。

進一步地，所述步驟4中，單站單目標定位算法為：

由剩余重力矢量gⁱ(t)的水平分量和垂直分量以及剩余重力梯度張量tⁱ(t)的水平分量和將多個運動目標體看成一個整體，根據(jù)公式(5)計算第i個重力觀測基站上估計的該整體的位置：

進一步地，所述運動目標體位于太空、地表或水面以下、地表或水面以上，或位于地表面或水面。

進一步地，所述重力觀測基站固定在地表面或地表以上某個位置，或位于船只、飛機或衛(wèi)星等移動平臺上。

本發(fā)明原理為：

(1)單運動目標體定位算法

假設(shè)運動目標體的體積為λ，剩余密度為ρ，第i個重力觀測基站p(aⁱ,bⁱ,cⁱ)上的剩余重力位函數(shù)φⁱ滿足下列積分方程：

式中為運動目標體與重力觀測基站p(aⁱ,bⁱ,cⁱ)之間的距離，δxⁱ＝x-aⁱ，δyⁱ＝y(tǒng)-bⁱ，δzⁱ＝z-cⁱ，x,y,z為運動目標體的位置，g＝6.67×10^-11kg^-1m³s^-2為引力常數(shù)。

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

如果重力觀測基站與運動目標體之間的距離足夠大，運動目標體可以被看成一個質(zhì)點，此時，公式(7)和公式(8)可以寫成如下形式：

使用公式(9)的重力矢量和公式(10)中的水平重力梯度張量，則在第i個重力觀測基站上所計算的運動目標體的質(zhì)心坐標和剩余質(zhì)量可從公式(11)獲得。

(2)多運動目標體定位算法

假設(shè)空間中存在n個運動目標體，其剩余密度為ρj(j＝1,2,…,n)。根據(jù)公式(9)和(10)，可以得到第i個重力觀測基站上的剩余重力矢量gⁱ和重力梯度張量tⁱ：

式中，為運動目標體o(x^j,y^j,z^j)與重力觀測基站p(aⁱ,bⁱ,cⁱ)之間的距離，δx^ij＝x^j-aⁱ，δy^ij＝y(tǒng)^j-bⁱ，δz^ij＝z^j-cⁱ；x^j,y^j,z^j和mj為運動目標體o(x^j,y^j,z^j)的位置和質(zhì)量。

此時，可以被寫成如下形式：

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

結(jié)合公式(13)和(14)，用來表示和

設(shè)公式(15)和(16)可以改寫為：

公式(17)依然受到運動目標體質(zhì)量mj的影響，為了消除質(zhì)量的影響，考慮公式(14)，將公式(17)除以便可得到歸一化非線性方程：

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

同理，將公式(17)分別除以和可得：

其中，

本發(fā)明利用標準的levenberg-maquardt算法求解公式(18)和公式(19)所組成的非線性方程組。眾所周知，對于這種迭代的方法，初始值的選擇直接影響了計算結(jié)果的可靠性。本發(fā)明給出了一種數(shù)值穩(wěn)定且高精度的單站單目標定位算法，并在單站單目標定位算法的基礎(chǔ)上給出了一個挑選位置初始值的可行的方案：將多個運動目標體看成一個整體，通過單站單目標定位算法計算每個重力觀測基站上估計的該整體的位置，形成位置信息矩陣從中分別為n個運動目標體自動挑選一組坐標作為其位置初始值，形成初始位置矩陣x0(t)＝(x0(t),y0(t),z0(t))n×3；進而基于levenberg-marquardt(中文為列文伯格-馬夸爾特法，簡稱lm算法)算法求解非線性性方程組，得到的運動目標體在t時刻的質(zhì)心坐標位置：令fun為公式(18)和公式(19)所組成的非線性方程組，調(diào)用函數(shù)[x,fval,flag]＝fsolve(@fun,x0(t),option)，當fval小于給定的閾值時，將此時的x＝(x(t),y(t),z(t))n×3作為求解得到的運動目標體在t時刻的質(zhì)心坐標位置。

將得到的每個運動目標體在t時刻的質(zhì)心坐標位置帶入公式(12)和公式(13)，求解線性方程組，便可得到運動目標體的剩余質(zhì)量mj。

有益效果：

本方法給出了一種利用重力觀測基站的重力矢量及其梯度張量監(jiān)測多個運動目標體的坐標位置(運動軌跡)和質(zhì)量的算法，可以實時地監(jiān)測運動目標體的位置和質(zhì)量，精度高，為下一步軍用或者民用成果轉(zhuǎn)化奠定了基礎(chǔ)。通過觀測剩余重力矢量和剩余重力梯度張量隨時間的變化，即可求出多個運動目標體的運動軌跡。本發(fā)明中運動目標體可以位于外太空、地表面及以下，如地下洞穴及其中運動的汽車、人員，常規(guī)潛艇及核潛艇，航空母艦及各種水面艦艇等；也可以位于空氣中，如飛機、導(dǎo)彈等，適用范圍廣。

附圖說明

圖1為本發(fā)明利用多個重力觀測基站進行多目標體運動軌跡識別的示意圖，其中oi表示第i個觀測點，sj表示第j個目標體，rij表示其間距離；

圖2為無噪時雙運動目標體雙站定位結(jié)果圖，其中實線為運動目標體a和b的實際航行軌跡，空心三角形為運動目標體a的定位結(jié)果，空心圓為運動目標體b的定位結(jié)果，空心五角星為重力觀測基站。圖2(a)分別為運動目標體x-y軸定位結(jié)果，圖2(b)為運動目標體x-z軸定位結(jié)果。

圖3為加入方差為0，標準差為0.03e的高斯噪聲雙運動目標雙站定位結(jié)果圖，其中實線為運動目標體a和b的實際航行軌跡，空心三角形為運動目標體a的定位結(jié)果，空心圓為運動目標體b的定位結(jié)果，空心五角星為重力觀測基站。圖3(a)分別為運動目標體x-y軸定位結(jié)果，圖3(b)為運動目標體x-z軸定位結(jié)果。

具體實施方式

以下結(jié)合附圖和具體實施方式對本發(fā)明作進一步的說明。

本發(fā)明公開了一種基于重力矢量及其梯度張量的多運動目標體監(jiān)測方法，包括以下步驟：

(1)重力觀測基站布置設(shè)計：假設(shè)空間中存在n個運動目標體，根據(jù)觀測目標確定觀測基站的位置，并在m個重力觀測基站上安裝重力儀和重力梯度儀；

(2)計算重力觀測基站上水下或空中多個運動目標體引起的剩余重力矢量g(t)和剩余重力梯度張量t(t)；

(3)由每一個重力觀測基站上測量的剩余重力矢量g(t)的水平分量gx(t)、gy(t)和垂直分量gz(t)，以及剩余重力梯度張量t(t)的水平分量txx(t)、tyy(t)和txy(t)，將空中或者水下多個目標體看成一個整體，根據(jù)公式(5)初步計算每一個重力觀測基站上估計的多個運動目標體整體的位置；

(4)實時坐標與質(zhì)量的計算與監(jiān)測：整理基站采集的數(shù)據(jù)，以步驟(3)估計的運動目標體的整體的位置為初始值，利用公式(6)～(20)所示的計算方法，實時確定多運動目標體的坐標和質(zhì)量等。

實施例1：

本實施例為進行水下多運動目標體軌跡定位的實例。假設(shè)重力梯度儀的觀測精度為0.01e，運動目標體a和b的質(zhì)量分別為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可以看出，本發(fā)明算法估計的運動目標體的軌跡與實際位置完全吻合，所識別a和b的x、y、z的平均相對誤差分別為：0.1％、0.4％、0.4％和1.4％、1.2％、0.6％。其中相對誤差r＝|(估計值-真實值)/真實值|×100。說明我們的算法能夠精準的實時監(jiān)測雙目標體的運動軌跡，適用性強。

為了檢驗本發(fā)明的穩(wěn)定性和可靠性，有必要考慮噪聲影響。在上述模型中加入方差為0，標準差為0.03e(3倍儀器觀測精度)的高斯噪聲，其計算結(jié)果見圖3。如圖3所示，加入方差為0，標準差為0.03e的高斯噪聲，本發(fā)明算法同樣可以定位運動目標體的軌跡，所識別a和b的x、y、z的平均相對誤差分別為：0.7％、1.5％、4.7％和2.3％、2.5％、3.2％。上述分析說明該方法穩(wěn)定性強，具有一定的抗噪能力。

分析表明，采用本發(fā)明所提供的基于重力矢量及其梯度張量的多運動目標體監(jiān)測方法，可以實時地監(jiān)測具有復(fù)雜運動軌跡的多個運動目標體的位置，不僅具有較高的民用價值，也具有較大的國防應(yīng)用價值。