一種基於對數似然估計的無源定位跟蹤新方法與流程
2023-05-24 04:51:11 1
本發明涉及自動控制和信息技術領域，可用於在軍事和民用領域中對高機動運動目標的實時定位跟蹤，具體涉及一種基於對數似然估計的無源定位跟蹤新方法。
背景技術：
無源定位技術是一種定位設備本身不發射信號，僅僅是依靠被動地接收輻射源的信息來實現定位的技術。它可以利用未知位置的輻射源的輻射信息，確定出該輻射源的類型、空間和地理位置；或者利用已知地理位置的輻射源來確定航行中物體的空間和地理位置，這也是導航和制導定位中的一項重要技術手段。與有源定位技術相比，無源定位技術具有作用距離遠、隱蔽接收、不易被對方發覺等優點，是現代一體化防空系統、機載對地對海攻擊以及對付隱身目標的遠程預警系統的重要組成部分，對於提高系統在電子戰環境下的生存能力和作戰能力具有重要作用，因此無源定位跟蹤技術一直是研究的熱點和難點。
無源定位是通過觀測站接收來自目標的無線電信號，並從信號中挖掘出用於定位的觀測量。一般的觀測量包括到達時間(Time of Arrival，TOA)，接受信號強度(Received Signal Strength，RSS)，到達時間差(Time Difference of Arrival，TDOA)，到達方位角和俯仰角(Angle of Arrival，AOA)[8,9],到達頻率差(FrequencyDifference of Arrival，FDOA)等。根據上述觀測信息均能夠建立關於目標位置或速度與觀測站位置之間的(非線性)方程，再通過優化求解該方程即可獲得關於目標位置或速度的參數信息。近些年來，基於上述觀測量的目標定位算法已相繼提出，其中包括Taylor級數迭代算法，總體最小二乘(Total Least Squares，TLS)算法,約束加權最小二乘(Constrained Weighted Least Squares，CWLS)算法，約束總體最小二乘(Constrained Total Least Squares，CTLS)算法，結構總體最小二乘(Structured Total Least Squares，STLS)算法。然而，上述算法大都需要迭代運算，這除了帶來較複雜的運算量外，還會出現迭代發散和局部收斂等問題。
技術實現要素：
為了克服上述現有技術的不足，本發明的目的是提供一種基於對數似然估計的無源定位跟蹤新方法，解決了傳統無源定位跟蹤方法中迭代運算出現迭代發散、收斂於局部解以及運算量大的問題。本發明針對多基站測角無源定位跟蹤技術，利用對數似然估計將多站測角無源定位非線性觀測方程進行偽線性化處理，用最大似然估計求出目標的位置並給出了目標定位算法的閉式解，接著利用「當前」統計機動模型和卡爾曼濾波，在線實時計算目標位置、速度和加速度，實現對目標的實時精確定位跟蹤。通過理論證明和仿真實驗證明本發明提出閉式解的定位性能均能夠達到克拉美羅下限(Cramér-Rao Bound，CRB)，從而驗證文中理論分析的有效性。
一種基於對數似然估計的無源定位跟蹤新方法，其特徵在於，包括如下步驟：
1)根據基站測角與目標位置之間數學關係，建立目標定位跟蹤測量方程，測量噪聲服從高斯分布，同時確定目標狀態變量；
在直角坐標系下，假設觀測站i的位置分別為(xi,yi,zi)，i＝1,2,…,N，飛行目標的空間位置在k時刻為(xk,yk,zk)，觀測站i測得飛行目標在k時刻的方位角為Ai,k和俯仰角為Ei,k，在存在測量噪聲情況下，Ai,k和Ei,k可以用下式(1)、式(2)表示：
這裡定義和以式(3)、式(4)表示：
和表示觀測站i測得方位角和俯仰角的真值，ni,k和ei,k是互不相關的測量噪聲序列，服從均值為0的高斯分布，其方差分別由下式給出：
重新定義測量序列，用向量ak表示為式(5)，
ak＝[A1,k,A2,k,…,AN,k,E1,k,E2,k,…,EN,k]T (5)
同時測量真值序列用向量表示為式(6)，
根據式(2)、(5)、(6)，可以得出下式(7),
這裡
bk＝[n1,k,n2,k,…,nN,k,e1,k,e2,k,…,eN,k]T
由於噪聲序列ni,k和ei,k是互不相關的零均值高斯隨機變量，方差分別為和噪聲序列bk是2N維向量，均值為0，方差矩陣由式(8)給出：
這裡
為了利用含噪聲項向量ak估計飛行目標在k時刻的位置(xk,yk,zk)，定義如下向量sk，向量sk由式(9)表示：
向量元素(xk,yk,zk)給出了飛行器在k時刻的位置，給出了k時刻飛行器在(x,y,z)方向上的飛行速度，給出了k時刻飛行器在(x,y,z)方向上的飛行加速度，
假定sk的最大似然估計和它的誤差協方差矩陣pk/k-1先於ak被計算出來，和pk/k-1定義如下形式向量，由式(10)、式(11)表示，
2)利用新的量測序列ak更新sk的最大似然估計和誤差協方差矩陣pk/k-1。sk的更新估計將用表示，更新誤差協方差矩陣用pk/k表示，根據最大似然估計原理，建立求解目標狀態變量的對數似然函數，通過求導可以計算出目標位置的全局最優解；
因為最大似然估計是漸近的高斯函數分布，因此我們可以假定估計值的概率密度函數是均值為sk，方差為pk/k-1的高斯函數，的概率密度函數定義由式(12)表示：
這裡ρ1為：
因為在式(7)中，bk是均值為0，方差為Rk的高斯高斯隨機向量，所以向量ak是均值為方差為Rk的高斯隨機向量，因此ak的概率密度函數g(ak)由式(13)給出，
這裡ρ2為：
在給定和ak的條件下，sk的似然函數l(sk)表示如式(14)：
由此，負對數似然函數L(sk)經計算等於式(15)：
這裡常數項被省略，在給定和ak的條件下，為了得到sk的最大似然估計，必須最小化-L(sk)，由式(16)表示：
3)為了便於快速實時濾波，利用對數似然估計將多站測角無源定位非線性觀測方程進行偽線性化處理，給出了目標位置計算的閉式解；
根據(13)式和(14)式，飛行目標位置坐標狀態變量xk，yk，zk與測量值Ai,k和Ei,k存在非線性形式，為了得到狀態變量與測量之間的線性關係，我們將在線性化這兩個等式，
對(13)式在進行泰勒級數展開，忽略掉高次項，得式(17)：
這裡
符號和表示估計的觀測站與目標的方位角和水平距離，觀測站位於(xi,yi,zi),為從(17)式估計目標在k時刻的位置，將(17)式進行數學變換得到下式：
這裡
對(14)式在進行泰勒級數展開，忽略掉高次項，得到下式：
其中
和同上，表示估計的觀測站與目標的距離和方位角。為從(19)式估計目標在k時刻的位置，將(19)式進行數學變換，得到下面的關係式(20)：
其中
通過對量測Ai,k和Ei,k重複上面提到的線性化過程，可以得到下面偽測量向量，由式(21)表示：
向量也可以表示成矩陣形式，如式(22)，
這裡
將上式計算結果代入(15)式，可以得到-L(sk)另外一個表達式(23)，
因此，為了得到狀態向量sk在k時刻的最大似然估計，必須最小化-L(sk)；
4)通過偽線性化處理，給出了目標位置狀態變量真實值與預測值之間的關係式，在此基礎上建立目標定位跟蹤狀態方程；
為了求解sk的估計值求-L(sk)對sk的梯度▽(-L(sk))，令▽(-L(sk))＝0，得到式(24)，
因為
所以(23)式求出的是全局最小值，在(23)式中Gk由式(25)表示，
矩陣Gk叫增益矩陣，經變換後得到下面的等式(26)，
的誤差協方差矩陣pk/k由式(27)給出
由上述推導結果可以看出，在給定pk/k-1和ak時，和pk/k計算過程如下：
第一步：計算
第二步：計算
第三步：計算pk/k＝(I-GkHk)pk/k-1；
5)運用「當前」統計模型和卡爾曼濾波，結合測量方程和狀態方程，建立機動目標定位跟蹤的遞推關係式，實現基於測角的機動目標無源定位跟蹤；
在上述中，當計算和pk/k，假定和它的誤差協方差矩陣pk/k-1已經給定，因此在k+1時刻計算和pk+1/k+1時，必須事先給定和pk+1/k的最大似然估計，在這部分將給出如何利用和pk/k估計和pk+1/k值，運動目標從k時刻到k+1時刻機動模型採用「當前」統計模型，即式(28)所表示，
這裡狀態轉移矩陣F由式(29)表示，
其中
U由
U＝[u1 u2 u3 u4 u5 u6 u7 u8 u9]T (30)
其中
這裡，F表示狀態轉移矩陣，T表示採樣間隔，(ax,ay,az)表示(x,y,z)方向機動，表示隨機加速度的均值(加速據)時間常數的倒數，W表示狀態噪聲，服從零均值高斯分布，方差為Q，U表示加速度係數矩陣。
是加速度的方差，由估計理論，可以求出和pk+1/k，由式(32)、式(33)表示，
本發明的有益效果是：
1)運用最大似然估計及對數似然函數偽線性化處理，使得目標位置計算具有閉式解，而且計算結果具有全局最優解，避免了傳統無源定位方法中迭代運算出現迭代發散、收斂於局部解的問題；
2)在利用對數似然估計將多站測角無源定位非線性觀測方程進行偽線性化處理的基礎上，運用「當前」統計模型和卡爾曼濾波算法，實現了在線實時的目標定位跟蹤，解決了傳統迭代計算複雜度高的問題；
3)通過理論證明和仿真實驗驗證，本發明提出閉式解的定位性能均能夠達到克拉美羅下限(Cramér-Rao Bound，CRB)，從而驗證本方法的有效性。
附圖說明
圖1為本發明的第一種仿真位置估計曲線；
圖2為本發明的第一種仿真速度估計曲線；
圖3為本發明的第一種仿真加速的估計曲線；
圖4為本發明的第二種仿真位置估計曲線；
圖5為本發明的第二種仿真速度估計曲線；
圖6為本發明的第二種仿真加速度估計曲線；
圖7為本發明的位置估計RMS與Cramér Rao Bound下限。
具體實施方式
以下結合實施例對本發明進一步敘述。
一種基於對數似然估計的無源定位跟蹤新方法，其特徵在於，包括如下步驟：
1)根據基站測角與目標位置之間數學關係，建立目標定位跟蹤測量方程，測量噪聲服從高斯分布，同時確定目標狀態變量；
在直角坐標系下，假設觀測站i的位置分別為(xi,yi,zi)，i＝1,2,…,N，飛行目標的空間位置在k時刻為(xk,yk,zk)，觀測站i測得飛行目標在k時刻的方位角為Ai,k和俯仰角為Ei,k，在存在測量噪聲情況下，Ai,k和Ei,k可以用下式(1)、式(2)表示：
這裡定義和以式(3)、式(4)表示：
和表示觀測站i測得方位角和俯仰角的真值，ni,k和ei,k是互不相關的測量噪聲序列，服從均值為0的高斯分布，其方差分別由下式給出：
重新定義測量序列，用向量ak表示為式(5)，
同時測量真值序列用向量表示為式(6)，
根據式(2)、(5)、(6)，可以得出下式(7),
這裡
bk＝[n1,k,n2,k,…,nN,k,e1,k,e2,k,…,eN,k]T
由於噪聲序列ni,k和ei,k是互不相關的零均值高斯隨機變量，方差分別為和噪聲序列bk是2N維向量，均值為0，方差矩陣由式(8)給出：
E[bkbl]＝δklRk (8)
這裡
為了利用含噪聲項向量ak估計飛行目標在k時刻的位置(xk,yk,zk)，定義如下向量sk，向量sk由式(9)表示：
向量元素(xk,yk,zk)給出了飛行器在k時刻的位置，給出了k時刻飛行器在(x,y,z)方向上的飛行速度，給出了k時刻飛行器在(x,y,z)方向上的飛行加速度，
假定sk的最大似然估計和它的誤差協方差矩陣pk/k-1先於ak被計算出來，和pk/k-1定義如下形式向量，由式(10)、式(11)表示，
2)利用新的量測序列ak更新sk的最大似然估計和誤差協方差矩陣pk/k-1。sk的更新估計將用表示，更新誤差協方差矩陣用pk/k表示，根據最大似然估計原理，建立求解目標狀態變量的對數似然函數，通過求導可以計算出目標位置的全局最優解；
因為最大似然估計是漸近的高斯函數分布，因此我們可以假定估計值的概率密度函數是均值為sk，方差為pk/k-1的高斯函數，的概率密度函數定義由式(12)表示：
這裡ρ1為：
因為在式(7)中，bk是均值為0，方差為Rk的高斯高斯隨機向量，所以向量ak是均值為方差為Rk的高斯隨機向量，因此ak的概率密度函數g(ak)由式(13)給出，
這裡ρ2為：
在給定和ak的條件下，sk的似然函數l(sk)表示如式(14)：
由此，負對數似然函數L(sk)經計算等於式(15)：
這裡常數項被省略，在給定和ak的條件下，為了得到sk的最大似然估計，必須最小化-L(sk)，由式(16)表示：
3)為了便於快速實時濾波，利用對數似然估計將多站測角無源定位非線性觀測方程進行偽線性化處理，給出了目標位置計算的閉式解；
根據(13)式和(14)式，飛行目標位置坐標狀態變量xk，yk，zk與測量值Ai,k和Ei,k存在非線性形式，為了得到狀態變量與測量之間的線性關係，我們將在線性化這兩個等式，
對(13)式在進行泰勒級數展開，忽略掉高次項，得式(17)：
這裡
符號和表示估計的觀測站與目標的方位角和水平距離，觀測站位於(xi,yi,zi),為從(17)式估計目標在k時刻的位置，將(17)式進行數學變換得到下式：
這裡
對(14)式在進行泰勒級數展開，忽略掉高次項，得到下式：
其中
和同上，表示估計的觀測站與目標的距離和方位角。為從(19)式估計目標在k時刻的位置，將(19)式進行數學變換，得到下面的關係式(20)：
其中
通過對量測Ai,k和Ei,k重複上面提到的線性化過程，可以得到下面偽測量向量，由式(21)表示：
向量也可以表示成矩陣形式，如式(22)，
這裡
將上式計算結果代入(15)式，可以得到-L(sk)另外一個表達式(23)，
因此，為了得到狀態向量sk在k時刻的最大似然估計，必須最小化-L(sk)；
4)通過偽線性化處理，給出了目標位置狀態變量真實值與預測值之間的關係式，在此基礎上建立目標定位跟蹤狀態方程；
為了求解sk的估計值求-L(sk)對sk的梯度▽(-L(sk))，令▽(-L(sk))＝0，得到式(24)，
因為
所以(23)式求出的是全局最小值，在(23)式中Gk由式(25)表示，
矩陣Gk叫增益矩陣，經變換後得到下面的等式(26)，
的誤差協方差矩陣pk/k由式(27)給出
由上述推導結果可以看出，在給定pk/k-1和ak時，和pk/k計算過程如下：
第一步：計算
第二步：計算
第三步：計算pk/k＝(I-GkHk)pk/k-1；
5)運用「當前」統計模型和卡爾曼濾波，結合測量方程和狀態方程，建立機動目標定位跟蹤的遞推關係式，實現基於測角的機動目標無源定位跟蹤；
在上述中，當計算和pk/k，假定和它的誤差協方差矩陣pk/k-1已經給定，因此在k+1時刻計算和pk+1/k+1時，必須事先給定和pk+1/k的最大似然估計，在這部分將給出如何利用和pk/k估計和pk+1/k值，運動目標從k時刻到k+1時刻機動模型採用「當前」統計模型，即式(28)所表示，
這裡狀態轉移矩陣F由式(29)表示，
其中
U由
U＝[u1 u2 u3 u4 u5 u6 u7 u8 u9]T (30)
其中
這裡，F表示狀態轉移矩陣，T表示採樣間隔，(ax,ay,az)表示(x,y,z)方向機動，a表示隨機加速度的均值(加速據)時間常數的倒數，W表示狀態噪聲，服從零均值高斯分布，方差為Q，U表示加速度係數矩陣。
是加速度的方差，由估計理論，可以求出和pk+1/k，由式(32)、式(33)表示，
實施例
為了驗證算法的可行性，我們作以下三種仿真分析。
第一種：設目標在平面上作勻加速運動，起始坐標在原點，起始速度v0＝2m/s，加速度a0＝0.2m/s2。角度測量誤差服從均值為0，方差為0.001o的高斯分布。仿真中取T＝1s，機動常數為10，兩個觀測站的位置o1(0,50)、o2(0,100)。圖(1)給出了跟蹤曲線結果，圖(2)給出了速度估計結果，圖(3)給出了加速度估計結果。由圖可以看出位置、速度，加速度估計值和理論值很接近。位置的均方根誤差0.032，速度的均方根誤差為0.0897，加速度均方根誤差0.176。
第二種：設目標在平面上作s型運動，起始坐標在原點。其運動曲線由下面方程決定：
x＝10sin(5tpi/180)
初始條件同上。圖(4)給出了跟蹤曲線結果，圖(5)給出了速度估計結果，圖(6)給出了加速度估計結果。由圖可以看出位置、速度，估計值和理論值很接近，加速度估計值與理論值相差25％，在一定誤差範圍內，這種估計結果可以接受。經過計算位置的均方根誤差0.16，速度的均方根誤差為0.267，加速度均方根誤差0.47。
在經過100次蒙特卡羅仿真後，雖然每次估計的曲線不同，但統計結果基本相同，尤其是位置估計相對誤差小於4％。這說明採用上述不僅可以獲得高精度目標定位精度，而且提供了精度較高的速度和加速度估計。但是，當目標運動非常複雜時，加速度的估計精度可能下降。
第三種，仿真條件兩觀測站的位置坐標分別為：(-40km，0，0)，(-40km，0，0)，採樣時間間隔為1s，方位角和俯仰角的觀測誤差為5mrad，目標的初始位置是(-10km,50km,10km)。運動軌跡分三個階段：
運動階段1：勻速直線運動,各方向分速度為：vx＝0.2m/s，vy＝0.1m/s，vz＝0m/s，運行時間為1500s；
運動階段2：勻速圓周運動,角速度為0.157rard/s，向心加速度為74m/s2，線速度為471m/s，圓半徑為3km，運行時間為500s；
運動階段3：勻速直線運動,各方向分速度為:vx＝0.2m/s，vy＝0.2m/s，vz＝0.2m/s，運行時間為1000s。
本文經過100次蒙特卡羅仿真實驗，實驗結果如圖7所示，由實驗結果可以看出，利用本文方法對目標位置進行計算，其位置估計均方根誤差(RMS)漸近達到了克拉美羅下限(Cramér-Rao Bound，CRB)。