本發(fā)明涉及一種勘探地球物理領域的目標體定位方法，特別涉及一種基于旋轉重力數據的運動目標體的監(jiān)測方法。
背景技術：
：在地球重力學中，重力梯度張量除用于勘探地球地質結構構造，尋找油氣、礦產資源外，還有其他多種應用。例如，利用衛(wèi)星重力梯度，高精度恢復中短波的地球重力場；利用高精度重力梯度張量，為艦艇或巡航導彈進行導航；2010年，MajidBeikiandLaustB.Pedersen(MajidBeikiandLaustB.Pedersen，Eigenvectoranalysisofgravitygradienttensortolocategeologicbodies，Geophysics,2010)應用重力梯度張量特征值與特征向量、垂直方向的重力場計算了地下單一異常體的位置與質量，并被用于重力梯度張量剖面數據的定性解釋。2014年，Lockerbie(NALockerbie，Thelocationofsubterraneanvoidsusingtensorgravitygradiometry，Class.QuantumGrav.2014)基于2個不同測點處的重力梯度張量特征向量唯一確定了單一目標體的位置，并應用于地下掩體范圍的大致圈定。現有技術中，必須基于2個不同測點處的重力梯度張量確定單一目標體位置測量，且只能應用于地下掩體范圍的大致圈定。因此，有必要設計一種基于一個觀測基站的旋轉重力矢量和重力梯度張量進行單個運動物體定位和跟蹤的解析方法。技術實現要素：本發(fā)明所解決的技術問題是，針對現有技術的不足，提供了一種基于旋轉重力場的運動目標體監(jiān)測方法，本發(fā)明拋開了傳統(tǒng)面積性數據采集方式，給出了利用一個觀測站所測得的重力矢量和重力梯度張量進行旋轉監(jiān)測單個運動目標體的坐標位置和質量的核心算法，可以實時地監(jiān)測運動目標體的軌跡和質量。本發(fā)明的技術方案為：一種基于旋轉重力場的運動目標體監(jiān)測方法，包括以下步驟：步驟1、布置重力觀測基站：根據觀測目標確定重力觀測基站的位置，并在重力觀測基站安裝重力儀和重力梯度儀；步驟2、計算空間直角坐標系(x,y,z)【笛卡爾直角坐標系】下單個運動物體引起的剩余重力矢量g(t)和剩余重力梯度張量T(t)：定義空間直角坐標系(x,y,z)；利用重力儀和重力梯度儀按一定時間間隔記錄空間直角坐標系(x,y,z)下重力觀測基站的重力矢量gtotal(t)和重力梯度張量Ttotal(t)，根據公式(1)計算質量為M的單個運動物體引起的剩余重力矢量g(t)和剩余重力梯度張量T(t)：T(t)=Ttotal(t)-Ttotal(t0)g(t)=gtotal(t)-gtotal(t0)---(1)]]>上式中，Ttotal(t0)和gtotal(t0)分別為運動目標體對重力基站不產生影響時觀測到的重力矢量及重力梯度張量，t0稱為時間基點；步驟3、旋轉重力場，計算局部旋轉坐標系(x′,y′,z′)下的剩余重力矢量和剩余重力梯度張量：將計算得到的剩余重力矢量g(t)和剩余重力梯度張量T(t)按照zyz型歐拉角旋轉，令剩余重力矢量g(t)的方向為旋轉后的局部旋轉坐標系為(x′,y′,z′)，其中z′的方向為旋轉空間直角坐標系使得Tx′y′(t)的絕對值最小，得到旋轉后的局部旋轉坐標系中的gz′(t)和T′(t)：T′(t)=Tx′x′(t)0Tx′z′(t)0Ty′y′(t)Ty′z′(t)Tz′x′(t)Tz′y′(t)Tz′z′(t)---(2)]]>此時局部旋轉坐標系中剩余重力矢量的垂直分量gz′(t)與剩余重力矢量g(t)的模g(t)相等，即gz′(t)＝g(t)，且gx′(t)＝gy′(t)＝0，其中gx(t)、gy(t)、gz(t)為剩余重力矢量g(t)的三個分量；步驟4、由局部旋轉坐標系中剩余重力矢量的垂直分量gz′(t)及剩余重力梯度張量T′(t)各分量的值Tx′x′(t)，Ty′y′(t)，Tz′z′(t)，根據公式(3)、(4)和(5)計算t時刻在局部旋轉坐標系下重力觀測基站與運動目標體質心的距離，即矢徑dz′(t)：當Tx′x′(t)＝Ty′y′(t)時，dz′(t)=-2gz′(t)Tz′z′(t)---(3)]]>當Tx′x′(t)＞Ty′y′(t)時，dz′(t)=gz′(t)Tx′x′(t)---(4)]]>當Tx′x′(t)<Ty′y′(t)時，dz′(t)=gz′(t)Ty′y′(t)---(5)]]>步驟5、由t時刻在局部旋轉坐標系下重力觀測基站與運動目標體質心的距離dz′(t)，經過坐標轉換，得到t時刻運動目標體在空間直角坐標系(x,y,z)【笛卡爾直角坐標系】下任意時刻的質心坐標位置：{xc(t),yc(t),zc(t)}T={x(t),y(t),z(t)}T-z^′·dz′(t)---(6)]]>上式中，x(t),y(t),z(t)為t時刻重力觀測基站的坐標位置，G為引力常數，xc(t),yc(t),zc(t)為t時刻運動目標體的質心坐標位置。運動目標體t時刻的質量M(t)為：M(t)=|dz′2(t)·gz′(t)·(1-G)|---(7).]]>所述步驟3中，重力場旋轉的準則是基于如公式(8)和(9)所示的zyz型歐拉角的旋轉方式：(x′,y′,z′)T＝R(x,y,z)T(8)R=cosψcosθcosφcosψcosθsinφ+sinψcosφcosψsinθ-sinψcosθcosφ-cosψsinφ-sinψcosθsinφ+cosψcosφ-sinψsinφ-sinθcosφ-sinθsinφcosθ---(9)]]>旋轉后的局部旋轉坐標系中的T′(t)計算公式為：T′(t)＝RT(t)RT(10)重力基站固定在地表面或地表以上某個位置，或位于船只、飛機或衛(wèi)星等移動平臺上。由運動目標體的質心坐標位置(運動軌跡)，根據公式(11)計算運動目標體相對于重力觀測基站的運動速度V＝(Vx,Vy,Vz)：Vx=x′(ti+1)-x′(ti)ti+1-tiVy=y′(ti+1)-y′(ti)ti+1-tiVz=z′(ti+1)-z′(ti)ti+1-ti---(11)]]>其中，ti+1為第i+1個時刻，ti為第i個時刻；假設重力基站的移動速度為Vs，則運動目標體的絕對速度為：V絕對速度＝Vs+V(12)。本發(fā)明原理為：假設目標體的體積為V，剩余密度為ρ，笛卡爾坐標系(垂直向下為z軸正方向)，觀察點(x,y,z)的重力位滿足的積分方程為：φ=-G∫VρdVr---(13)]]>剩余重力矢量g滿足下列積分方程：gz=-G∫∫∫vρz-zcr3dv,gx=-G∫∫∫vρx-xcr3dv,gy=-G∫∫∫vρy-ycr3dv---(14)]]>剩余重力梯度張量T滿足下列積分方程：Txx=-G∫∫∫vρr2-3(x-xc)2r5dv,Txy=-G∫∫∫vρ-3(x-xc)(y-yc)r5dv,Txz=-G∫∫∫vρ-3(x-xc)(z-zc)r5dv,Tyy=-G∫∫∫vρr2-3(y-yc)2r5dv,Tyz=-G∫∫∫vρ-3(y-yc)(z-zc)r5dv,Tzz=-G∫∫∫vρr2-3(z-zc)2r5dv---(15)]]>其中Txx+Tyy+Tzz＝0,(xc,yc,zc)∈V,r＝|(x,y,z)-(xc,yc,zc)|，V為目標體的體積。當觀測點與運動目標體距離較遠時，目標體便可以被看出為一個質點，則公式(14)和(15)可以寫成如下形式：gz=-GMz-zcr3,gx=-GMx-xcr3,gy=-GMy-ycr3---(16)]]>Txx=-GMr2-3(x-xc)2r5,Txy=-GM-3(x-xc)(y-yc)r5,Txz=-GM-3(x-xc)(z-zc)r5,Tyy=-GMr2-3(y-yc)2r5,Tyz=-GM-3(y-yc)(z-zc)r5,Tzz=-GMr2-3(z-zc)2r5---(17)]]>其中，x,y,z為觀測點位置，G為引力常數，xc,yc,zc為目標體重心位置，M為目標體質量。根據公式(16)和(17)之間的相互關系，可以得到以下三種定位方法，為了便于說明將此類定位方法稱之為普通重力場定位方法：當gz≠0，可以采用公式(18)所示的普通重力場定位方法1：Dxz=gxgz,Dyz=gygz,A=[Dxz2+Dyz2+1]12,dz=A2-3A2gzTzz,dx=Dxzdz,dy=Dyzdz,M=|dz2·gz·A3(1-G)|---(18)]]>當gx≠0，可以采用公式(19)所示的普通重力場定位方法2：Dzx=gzgx,Dyx=gygx,A=[Dzx2+Dyx2+1]12,dx=A2-3A2gxTxx,dy=Dyxdx,dz=Dzxdx,M=|dx2·gz·A3(1-GDzx)|.---(19)]]>當gy≠0，可以采用公式(20)所示的普通重力場定位方法3：Dzy=gzgy,Dxy=gxgy,A=[Dzy2+Dxy2+1]12,dy=A2-3A2gyTyy,dx=Dxydy,dz=Dzydy,M=|dy2·gz·A3(1-GDzy)|.---(20)]]>上述定位方法是基于運動球體或質點這一前提條件提出來的，在運動物體相對測點不能看成為質點的情況下，會導致識別出錯誤的結果；另外，公式(18)～(20)中存在多個分母項，當這些分母項符號發(fā)生變化時便會增加定位的計算誤差。為此，在上述定位方法的基礎上提出了局部旋轉重力場定位方法。下面將詳細討論局部旋轉重力場定位算法的原理和應用效果。將笛卡爾直角坐標系進行旋轉。(x,y,z)為使用者定義的笛卡爾坐標系，(x′,y′,z′)為旋轉后的坐標系。假設兩個坐標系有相同的坐標原點，可以互相旋轉，坐標系的軸或者向量的分量可以通過一個3D的旋轉矩陣聯系起來：(x′,y′,z′)T＝R(x,y,z)T(21)這里R為旋轉矩陣，通常用歐拉角(φ,θ,ψ)來表示。采用地球物理中常用的zyz型歐拉角，首先關于z軸旋轉φ，然后關于y′旋轉角度θ，最后關于z′旋轉ψ。這個旋轉矩陣可以用下式表示：R=cosψcosθcosφcosψcosθsinφ+sinψcosφcosψsinθ-sinψcosθcosφ-cosψsinφ-sinψcosθsinφ+cosψcosφ-sinψsinφ-sinθcosφ-sinθsinφcosθ---(22)]]>R是一個正交矩陣，它的逆矩陣就是它的轉置矩陣，(x,y,z)T＝RT(x′,y′,z′)T將旋轉矩陣應用到任意向量或者向量操作，則有：▿′=R▿---(23)]]>將計算得到的剩余重力矢量g和剩余重力梯度張量T按照zyz型歐拉角旋轉，令剩余重力矢量g的方向為旋轉后的局部旋轉坐標系為(x′,y′,z′)，其中z′的方向為旋轉空間直角坐標系使得Tx′y′的絕對值最小，得到旋轉后的局部旋轉坐標系中的gz′和T′：T′=Tx′x′0Tx′z′0Ty′y′Ty′z′Tz′x′Tz′y′Tz′z′---(24)]]>此時局部旋轉坐標系中剩余重力矢量的垂直分量gz′與剩余重力矢量g的模g相等，即gz′＝g，且gx′＝gy′＝0，其中gx、gy、gz為剩余重力矢量g的三個分量。在旋轉之后的局部參考系中，結合重力場一階導數和重力梯度之間的相互關系，利用gx′,gy′,gz′,Tz′z′和式(18)可以得到第一組目標體坐標位置與質量的計算公式：A=[(gx′gz′)2+(gy′gz′)2+1]12=1,dz′=A2-3A2gz′Tz′z′=-2gz′Tz′z′,dx′=gx′gz′·dz′=0,dy′=gy′gz′·dz′=0,M=|dz′2·gz′·(1-G)|.---(25)]]>同理，利用gx′,gy′,gz′,Tx′x′和式(19)可以得到第二組目標體坐標位置與質量的計算公式：A=[(gz′gx′)2+(gy′gx′)2+1]12=∞,dx′=A2-3A2gx′Tx′x′=gx′Tx′x′=0,dy′=gy′gx′·gx′Tx′x′=0,dz′=gz′gx′·gx′Tx′x′=gz′Tx′x′,M=|dz′2·gz′·(1-G)|.---(26)]]>同理，利用gx′,gy′,gz′,Ty′y′和式(20)可以得到第三組目標體坐標位置與質量的計算公式：A=[(gz′gy′)2+(gx′gy′)2+1]12=∞,dy′=A2-3A2gy′Ty′y′=gy′Ty′y′=0,dx′=gx′gy′·gy′Ty′y′=0,dz′=gz′gy′·gy′Ty′y′=gz′Ty′y′,M=|dz′2·gz′·(1-G)|.---(27)]]>經過比較，得到一套比較理想的方案，即：當Tx′x′＝Ty′y′時，dz′=-2gz′Tz′z′M=|dz′2·gz′·(1-G)|---(28)]]>當Tx′x′>Ty′y′時，dz′=gz′Tx′x′M=|dz′2·gz′·(1-G)|---(29)]]>當Tx′x′<Ty′y′時，dz′=gz′Ty′y′M=|dz′2·gz′·(1-G)|---(30)]]>經過坐標轉換便可以得到目標體在笛卡爾直角坐標系下的位置：(xc,yc,zc)T=(x,y,z)T-z^′·dz′---(31)]]>由物體質心軌跡，根據公式(32)計算物體相對于重力基站的運動速度V＝(Vx,Vy,Vz)：Vx=xc(ti+1)-xc(ti)ti+1-tiVy=yc(ti+1)-yc(ti)ti+1-tiVz=zc(ti+1)-zc(ti)ti+1-ti---(32)]]>其中，ti+1為第i+1個時刻，ti為第i個時刻。假設重力基站的移動速度為Vs，運動物體的絕對速度為：V絕對速度＝Vs+V(33)有益效果：本發(fā)明拋開了傳統(tǒng)面積性數據采集方式，給出了利用一個觀測站所測得的重力矢量和重力梯度張量進行旋轉監(jiān)測單個運動目標體的坐標位置和質量的核心算法，可以實時地監(jiān)測運動目標體的位置和質量。通過觀測重力梯度張量和重力矢量隨時間的變化，即可求出物體的運動軌跡。本發(fā)明解決了普通定位方法受質點或球體假設的限制，以及重力場或梯度在符號變化的鄰域內數值離散的問題，有效提高了定位精度；并且運動目標體可以位于地表面及以下，如地下洞穴及其中運動的汽車、人員，常規(guī)潛艇及核潛艇，航空母艦及各種水面艦艇等；也可以位于空氣中，如飛機、導彈，特別是隱身飛機。由于質量是物體的固有屬性，地球重力場及其梯度不受電磁波、聲波、熱等因素的影響，因此，目前的隱身方法相對于重力梯度是通明的。依據本發(fā)明方法，在理論上，只要布置合適的觀測臺網(固定臺網+移動臺站)，在適當的先驗信息下，就可以監(jiān)測特定區(qū)域一些特定目標的運動狀態(tài)，并根據其質量，評估其載彈量及其型號。配合其他監(jiān)測手段，可以獲得更詳細的目標信息。因此，本發(fā)明方法不僅具有較高的民用價值，也具有較大的國防應用價值。附圖說明：圖1為本發(fā)明利用一個重力觀測基站進行單體運動軌跡識別的示意圖。圖2為單運動目標體單站定位結果圖，其中實線為運動目標體實際航行軌跡，空心圓為常規(guī)普通重力場定位算法所識別的結果，三角形為本發(fā)明局部旋轉重力場定位算法所識別的結果；圖2(a)-(c)分別為運動目標體x、y、z軸定位結果，圖2(d)為運動目標體質量識別結果。圖3為加入5％高斯噪聲運動目標體單站定位結果圖，其中圖3(a)-(c)中實線為運動目標體實際航行軌跡，三角形為本發(fā)明局部旋轉重力場定位算法所識別的結果；圖3(d)為相對誤差對比圖，其相對誤差計算公式為：相對誤差＝|(估計值-真實值)/真實值|×100。具體實施方式：以下結合附圖和具體實施方式對本發(fā)明作進一步的說明。本發(fā)明涉及的定位方法包括以下步驟：步驟1、基站布置設計：確定觀測目標及儀器靈敏度范圍，合理設計觀測臺站的位置；在固定的臺站或者移動的基站上安裝高精度的重力儀和重力梯度儀。步驟2、異常重力場和重力梯度張量的計算：利用重力儀和重力梯度儀實時記錄基站的重力場與重力梯度張量，以便實時計算出重力場與重力梯度張量的異常變化，即剩余重力矢量值g和重力梯度值T；將所計算的剩余重力矢量值g和重力梯度值T按照zyz型歐拉角旋轉，得到旋轉之后的數據集gz′(t)和T′(t)。步驟3、實時坐標與剩余重量計算與監(jiān)測：整理基站采集數據，利用公式(1)～(12)所示的計算方法，實時確定單一移動目標體的坐標、質量和運動速度等。以下為本發(fā)明水下移動單體軌跡定位的實例。設單一棱柱體的質量為1.6×107kg，長170m，寬12.8m，高11m。運動軌跡為的拋物型方程，y＝0.1875x2+100，z＝200，Δt＝1,2,...,24。運動單體起始點x坐標為-1000米，重力觀察基站位于(0m，-100m，0m)處，共記錄24個時刻的重力梯度張量和垂直重力場。從圖2可以看出，在距離測點較遠的區(qū)域，普通重力場定位方法可以準確地識別棱柱體的軌跡和質量，但是當物體運動到(0m，100m，200m)附近區(qū)域時，追蹤識別的結果明顯偏離了物體的軌跡，這主要有以下兩個方面的原因：(1)在觀測點附近目標體不能看成為質點；(2)公式中的分母項在該區(qū)域符號發(fā)生了變化，由于數值離散不恰當增加計算誤差。而本發(fā)明局部旋轉重力場定位方法解決了上述兩個問題，識別棱柱體的軌跡和質量吻合了真實運動軌跡和質量。為了檢驗本發(fā)明的穩(wěn)定性和可靠性，有必要考慮噪聲影響。棱柱體運動軌跡為如圖3實線，目標體起始點x坐標為-1000m，觀測間隔為1/8分鐘，共記錄12個時刻的重力矢量及其梯度。觀測點坐標為(0m，-100m，0m)，在所測量的重力矢量及其梯度張量中分別加入5％的高斯隨機噪聲。其計算結果見圖3和表1。如圖3和表1所示，加入5％的噪聲，本發(fā)明算法同樣可以定位運動目標體的軌跡和質量，所識別的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，個別點偏離了物體的邊界，但偏差較小，均在允許的誤差范圍以內。當運動目標體距離測量點較遠時，偏差有所增大，這是因為此時重力梯度張量各分量的值較小，易受噪聲的影響。上述分析說明該方法穩(wěn)定性強，具有一定的抗噪能力。表1加入5％高斯噪聲運動目標體單站定位結果誤差表分析表明，采用本發(fā)明所提供的基于局部旋轉重力場的運動目標體監(jiān)測方法，可以實時地監(jiān)測具有復雜運動軌跡的目標體的位置和質量，比在直角坐標系下目標識別算法更具有適用性，不僅具有較高的民用價值，也具有較大的國防應用價值。當前第1頁1 2 3