基于非機動單站的目標運動狀態無偏估計

林建軍，王驤予涵，班曉軍*，黃顯林

（1.哈爾濱工業大學 航天學院，黑龍江 哈爾濱 150000；2.哈爾濱工程大學 航天與建筑工程學院，黑龍江 哈爾濱 150001）

1 引言

與有源定位相比，無源定位不需要向對方發射探測信號，具有隱蔽性良好、生存能力強的優點［1］。近年來，機載單站快速定位技術是發展最快的無源偵察定位技術之一，因其作用距離遠、機動性強的特點，在軍事上存在著巨大的應用價值。根據定位原理的不同，該技術可以分為到達角（Angle of Arrival，AOA），DOA，TDOA 等。隨著陣列天線及其信號處理技術的發展，AOA測量的硬件平臺成本降低，這使得基于AOA 的目標跟蹤技術具有更廣的適用范圍［2］。

三維到達角（Three-Dimension Angle of Arrival，3D-AOA）目標跟蹤是通過測量目標與運動站之間的方位角、俯仰角，進而對目標進行運動狀態估計的技術，其測量方程具有非線性。常見的基于角度信息的目標運動分析方法可以分為偽線性濾波法與非線性濾波法［3］，其中，偽線性濾波法由于其計算量小、不易發散的優點而倍受關注。文獻［4］在二維條件下對方位角進行偽線性處理，提出了適用于勻速運動目標的偽線性濾波器（Pseudo-Linear Estimator，PLE），相比于擴展卡爾曼濾波器（EKF）具有更高的穩定性，但該方法是一種有偏估計。為了解決PLE 的有偏性問題，文獻［5］提出了一種輔助變量法（IV），但沒有閉式解且對初值敏感；文獻［6］考慮了觀測方程噪聲特性提出了WIV，獲得了接近最大似然估計（MLE）的估計性能。通過添加約束條件求解偽線性方程的極小值，文獻［7］給出了漸近無偏的估計器，文獻［8］將該方法應用到二維勻速目標跟蹤上，文獻［9-10］則進一步考慮觀測平臺的位置誤差，對該方法進行了優化。在單站條件下，目標運動的可觀測信息量少，易出現目標運動信息不完全可觀的問題，若不能有效機動無法完成對目標的有效跟蹤［11］，可以通過增加傳感器數量［12］、引入速率先驗條件［13］、固定角度探測范圍［14］等方式，使傳感器在靜止條件下也能對目標進行有效跟蹤。然而，基于偽線性模型的目標運動狀態估計研究主要集中在小距離以及二維條件下，對遠距離三維場景下的目標定位跟蹤問題研究較少，且要求觀測平臺具有較為復雜的機動特性。在實際環境中，偵察機若處于安全狀態，通常會以巡航狀態工作。而借助于目標識別技術，能夠對目標型號進行有效識別，并給出相應的先驗速率參數。

本文在眾多偽線性測量研究的基礎上，針對巡航狀態下的目標，提出了一種觀測平臺非機動條件下基于三維到達角的目標運動狀態無偏估計法。該方法借助目標速率先驗的條件，以相對運動模型對系統進行建模，將偽線性最小二乘方法應用到該模型中；針對偽線性最小二乘存在有偏性的問題，通過添加約束條件對偽線性最小二乘法進行優化求解，達到消除模型有偏性的作用，從而得到一種具有漸進無偏特性的目標運動狀態估計方法。為了保證算法的運行效率，算法的解為閉式解，而根據測量噪聲特性對算法提出約束條件，使得問題轉換為帶約束的最優化問題，能夠在犧牲偽線性最小二乘方法一定的運行效率后大大提高算法的目標跟蹤精度。

2 基于非機動單站的3D-AOA建模

2.1 運動模型

當目標以巡航狀態工作時，運動狀態近似呈勻速直線運動。即使目標出現小機動其狀態方程成為時變系統，根據分段線性定常系統（Piece-Wise Const Systerm，PWCS）理論，仍可以認為目標在短時間內速度不變［15］。當觀測平臺以勻速運動形式巡航時，二者之間的相對運動仍保持勻速。在該場景下，定義觀測平臺的位置為po=[pox，poy，poz]T，速度為vo=[vox，voy，voz]T，目標的位置為pt=[ptx，pty，ptz]T，速度為vt=[vtx，vty，vtz]T。以目標與觀測平臺之間的相對運動作為研究對象，假設觀測平臺固定于點O，目標做勻速直線運動，其先驗速率為vc。如圖1 所示，以正東方向為X軸，正北方向為Y軸，天向為Z軸，建立右手系直角坐標系。目標與觀測站連線與Y軸夾角為方位角θ，與該連線的X-Y平面投影夾角為俯仰角φ，目標的位置與速度分別用p=pt-po=[px，py，pz]T，v=vt-vo=[vx，vy，vz]T表示。

圖1 東北天坐標系Fig.1 North-East-Up（ENU）coordinate system

tk表示k時刻（k=0，1，…），目標與觀測站的距離為rk，那么初始時刻距離為r0。忽略過程噪聲的影響，任意時刻tk的目標位置可表示為：

2.2 三維AOA 量測模型

三維空間中，k時刻方位角、俯仰角的測量方程為：

式中：假設nk，mk為服從高斯分布的測量噪聲，均值為0，方差分別為，測量方位角與俯仰角之間相互獨立，由此可得測量方程協方差矩陣

不失一般性，假設采樣時間間隔為T，將式（1）代入式（3）中進行代數變換：

然而，實際應用過程中無法獲得真值，因此用測量值代替，而會產生新的誤差向量。為了形式簡潔，假設初始時刻的角測量為真值，用·表示含有噪聲影響的變量。

經過推導可以得到噪聲向量εk為：

2.3 可觀性分析

不考慮測量噪聲，對方位角的測量方程進行偽線性處理。當k=1，2，…，k，k≥2 時，有：

定義偽觀測矩陣：

式中N=[I2×2，02×1]T。

令M=HTθ Hθ，則有：

經推導得到M的行列式值為：

3 三維定位算法

3.1 偽線性最小二乘法

利用觀測信息重新構建新觀測量的方法稱為偽線性化，而基于偽線性化觀測方程得到最小二乘解的方法稱為偽線性最小二乘法（Pseudolinear Least Square，PLE）。在上一節中，忽略了方位角與俯仰角的測量噪聲，證明了系統的可觀測性。然而，偽線性化處理使得量測方程中的噪聲特性發生變化，這會影響最小二乘法的估計精度。

將1 到k時刻的偽線性觀測方程寫成矩陣形式：

根據先驗速率vc，進一步可以得到目標的初始距離、位置矢量、目標與觀測平臺距離：

式中e0=[cosφ0sinθ0，cosφ0cosθ0，sinφ0]T。

與ε均受角測量噪聲的影響，從而具有相關性，這使得該估計方法并不是無偏估計，估計期望偏差為E{PLE-b}=-E{()-1}≠0。

3.2 約束總體最小二乘法

3.2.1 算法原理

基于偽線性觀測方程的PLE 方法是一種有偏估計。當角測量噪聲增大或者距離增大時算法性能會惡化，因此，本文提出了一種具有無偏性的約束總體最小二乘法（Constrain Total Least Square，CTLS），使得算法能夠適用于遠距離場景。PLE 本質上是一種以觀測方程均方誤差最小的方法，其損失函數表示為：

為了解決其有偏性問題，定義增廣矩陣，

那么，可將損失函數J重新表示為：

由于觀測噪聲的存在，測量矩陣與觀測量存在誤差，可將其分解為真值與誤差部分。于是，增廣矩陣可以表示為：

更進一步，將損失函數J展開，并對等式兩邊取期望，得到：

從理論上看，若E［J］能夠取最小值，即殘差均方差E[εTε]取最小時，v為最優解。但是，由于測量誤差的存在，v出現偏差。定義：

將噪聲項約束為常數，即設vTΩv為任意常數c，此時將不會影響E［J］取到極值時的v，選擇c=1。

于是，問題可重新被描述為：

拉格朗日極值法可處理帶約束條件的極值問題，構造函數：

對方程求導并令其導數為0，得：

于是，v*，λ*即為矩陣（P，Ω）的最小廣義特征向量與最小特征值，改進的漸進無偏估計b*為：

為了提高算法的運行效率，P，Ω可以由遞推形式獲得。

k時刻記矩陣P，Ω為P（k），Ω(k)，那么k+1時刻有：

同理可得：

算法步驟如表1 所示。

表1 約束總體最小二乘法算法步驟Tab.1 Steps of constrained total least squares algorithm

3.2.2 算法無偏性證明

將約束條件代入式（18）可得：

從理論上看，f(v，λ)的最小值也就是E[J]的最小值為1。對式（23）兩邊同時左乘vT，此時λ=1。將λ=1 代入式（23），并將Pv展開：

根據隨機矩陣理論，隨機矩陣的極限也是它的漸進期望，于是當測量數據量足夠大時有：

此時向量v∈null（A），等價于：

定義b的真值為b0，顯然其滿足Hb0=-y。根據可觀性理論可知，H為列滿秩矩陣，因此b只有唯一解，即：b=b0。那么，當測量數據量足夠大時，b依概率收斂于真值。

4 仿真分析

4.1 性能指標

克拉美羅下界（Cramer-Rao Lower Bound，CRLB）能夠給出狀態估計的理論最小方差。基于非機動條件下的觀測平臺對非機動目標進行運動狀態估計時，角度測量誤差與觀測平臺的位置誤差是目標估計精度的重要影響因素，而本文不考慮觀測平臺位置誤差的影響。定義狀態觀測量κ=[θ，φ]T，觀測站與目標之間距離為r，待估計向量p=[px，py，pz]T。

克拉美羅下界等于費舍爾矩陣I的逆，即：

于是，狀態估計的位置誤差與距離誤差的理論下界相同，可由式（35）得到：

其中：M表示蒙特卡洛仿真次數，，xk分別表示k時刻狀態后驗估計與真值。

時均誤差可用于衡量算法在一段時間的估計性能，其定義如下：

式中：U=N-L+1，N為跟蹤過程的總時刻數，L為開始記錄誤差的時刻數。

η用于衡量算法的提升精度，其方程如下：

式中：Index表示待評價指標參數，下標new，old分別表示初始方法與改進方法。

4.2 3D-AOA 的非機動單站運動仿真

4.2.1 不同角測量噪聲下的仿真對比

仿真場景設計如下：采樣時間為100 s，采樣周期為0.2 s，觀測站起始點坐標原點O，速度為[0，250，0]T，目標初始位置為[100，100，10]T，速度為[-340，80，-1.21]T，相對先驗速率vc為380 m/s。蒙特卡洛仿真次數為200，忽略過程噪聲。在光電跟蹤系統的角測量誤差水平能達到6 mrad［16］。因此，將觀測角噪聲分別假設為方差為0.1°～0.3°的高斯白噪聲，且獨立不相關。以相對距離誤差（Relative Distance Error，RDE）與絕對位置誤差（Absolue Position Error，APE）作為評價算法的定位精度的評價指標，分別為：

仿真場景下，觀測站與目標運動的位置態勢如圖2 所示。仿真對比了PLE，IV，WIV，EKF 以及約束總體最小二乘法。其中，EKF 設置初值b0|-1=03×1，初始狀態誤差的協方差P0|-1=10I3×3。

圖2 觀測站-目標運動態勢Fig.2 Observation station -target motion position

圖3 比較了各算法在100 s 時，不同角度測量標準差條件下的RDE 與APE。表2 與表3 比較了不同算法相比PLE 的精度。從圖3（a）可以看出，當測量噪聲為0.1°時，各算法均能有效收斂到10%以內的相對距離誤差。其中，約束總體最小二乘法與IV，WIV 方法誤差接近。隨著測量噪聲逐漸增大到0.2°，各算法性能受到不同程度的影響。PLE 與CTLS 的估計誤差均保持平穩上升，但PLE 的收斂精度明顯低于CTLS；EKF，IV，WIV 均出現較大的波動，性能不夠穩定，EKF 的估計精度反而有所提升，這可能與初值P的選擇是否恰當有關。當測量噪聲增大至0.3°時，EKF 與PLE 方法在100 s 時的相對距離誤差達到30% 以上，此時無法有效工作，而IV 與WIV 則出現發散傾向，RDE 與APE 均高于其他算法。其中，WIV 的絕對位置誤差超過100 km，但是CTLS 在100 s 時的RDE 仍能達到10%以內。根據表2 與表3，在100 s 時CTLS 的提升效果最為穩定與明顯，在角測量標準差為0.1°，0.15°，0.2°，0.25°，0.3°時，相對距離的估計精度相比于PLE 分別提高了70%，80.8%，82.9%，85.5%，85.9%；絕對位置的估計精度相比于PLE 分別提高了 65.4%，75.5%，78.3%，81.3%，81.8%，精度的提升效果隨著角測量標準差的增大逐步提升，算法的優勢逐漸明顯。

表2 不同角測量標準差下的RDE 及性能提升效果Tab.2 RDE and performance improvement effects under different angle measurement standard deviations（%）

表3 不同角測量標準差下的APE 及性能提升效果Tab.3 APE and performance improvement effects under different angle measurement standard deviations

圖3 不同角度測量標準差下的算法性能比較Fig.3 Comparison of algorithm performance at different angles of measurement standard deviation

為了進一步精確評估算法的精度與收斂性能，計算50～100 s 的平均時間偏差，如表4 所示。可以看出，在0.1°，0.2°，0.3°時，CTLS 的相對距離誤差與絕對距離誤差均為最低，這說明CTLS 相比于其他算法具有較高的收斂速度與定位精度。

表4 不同角測量誤差RDE 的平均時間偏差Tab.4 Time-ave deviation of RDE under different angle measurement standard deviations

4.2.2 不同初始距離條件下的仿真對比

在速度一定的條件下，目標起始點會影響角度觀測角的變化，因此進一步對初始距離的影響進行探究。取單位起始點位置矢量為[100，100，5]T（km），設置不同起始點，分別為起始點位置矢量的0.5～2 倍，間隔為0.25 倍，其他條件相同，角測量標準差均設置為0.1°。圖4 為算法在不同初始距離條件下的RDE 與APE。由圖4 可知，當取0.5 倍單位初始距離時，除了EKF外，其他算法的相對距離誤差與絕對距離誤差均十分接近；當初始距離逐漸提高至一倍初始距離時，各算法性能逐漸下降并有所區別，100 s 時按RDE 排序為PLE＞EKF＞IV＞WIV＞CTLS，APE 為EKF＞PLE＞IV＞WIV＞CTLS，各算法100 s 時均能達到10%以內的相對距離誤差；當初始距離增加到2 倍初始距離時，各種算法誤差大大提高，其中IV，WIV 出現較大波動，當初始距離為1.75～2 倍單位初始距離之間時，二者的APE 均大于PLE 與EKF，已超過100 km，對目標的運動信息估計作用已失效。而PLE，EKF，CTLS 的誤差均以較小幅度增加，CTLS 的誤差曲線明顯低于PLE，EKF。表5 與表6 比較了不同初始距離條件下不同算法相比PLE 的精度。CTLS 在初始距離分別為0.5 倍、1 倍、1.5 倍、2倍初始距離時，相對距離的估計精度相比PLE 分別提高33.3%，71.2%，80.9%，80.8%，絕對位置估計精度相比PLE 分別提高20%，65.4%，75.6%，75.8%，精度提升效果隨著初始距離的增大逐步提升，且提升效果始終最佳，性能穩定。

表5 不同初始距離下的RDE 及性能提升效果Tab.5 RDE and performance improvement effects under different initial distances（%）

表6 不同初始距離下的APE 及性能提升效果Tab.6 APE and performance improvement effects under different initial distances

圖4 不同初始距離的算法性能比較Fig.4 Comparison of algorithm performance for different initial distances

在50～100 s 內，不同初始距離條件下的時均RDE 與APE 如表7 所示。通過比較可以發現，當初始距離較近，取到0.5 倍單位初距時，各算法的性能接近，IV，WIV 與CTLS 能夠收斂到時均RDE 為1%，APE 為1 km 以內。而初始距離增加到單位初始距離的1 倍與2 倍時，其他算法均在CTLS 性能參數的1.5 倍以上。其中，在初始距離為單位初始距離的2 倍時，其他算法的50～100 s 時均RDE 均在60%以上。

表7 不同初距RDE 平均時間偏差Tab.7 Time-ave deviation of RDE under different initial distances

最后，以EKF 的一次蒙特卡洛仿真時間作為單位時間，各算法的相對運行時間如表8 所示。可以看出，CTLS 的算法運行速度與IV 相當，而WIV 的運行效率遠低于其他算法，且實驗過程中其定位并不明顯優于IV，與IV 一樣受測量噪聲與初始距離誤差的影響較大，具有性能不穩定的缺點。

表8 算法相對運行時間對比Tab.8 Comparison of algorithm's relative running time（s）

通過對比算法在不同測量誤差與不同初始距離下的仿真性能，可知本文提出的CTLS 相比于PLE，EKF 具有更高的定位精度，估計效果均較為穩定，不易出現較大波動；而相比于同樣具有無偏特性的IV，WIV，當測量噪聲較小或初始距離較近時，CTLS 均能達到十分接近克拉美羅下界的估計效果，但是當測量噪聲較高或初始距離較遠時，則具有更強的適應性。

5 結論

本文研究了百公里級別的單站無源定位問題，在運動平臺非機動條件下，利用先驗速率的已知條件，提出了一種漸進無偏的基于3D-AOA的目標跟蹤算法。以相對運動模型為研究對象，對系統的量測方程進行了偽線性處理，并對該系統的可觀性進行了分析。然后，介紹了基于偽線性測量方程的PLE 方法，針對PLE 的有偏問題，提出了一種具有漸進無偏性能的CTLS，對該算法原理進行了推導，并證明了算法的無偏性。最后，通過仿真分析在0.1°，0.2°，0.3°測量標準差時發現，CTLS 在50～100 s 內的平均RDE 分別為6%，12%，21%，平均APE 為9，19，35 km。同時，還比較了在不同初始距離條件下的算法性能，在50 km 左右的仿真場景中各算法性能接近，時均RDE 均能達到5%以內；而當距離逐漸提高到200 km 級別時，其他算法逐漸失效，而CTLS 在100 s 內已能明顯看出收斂趨勢，時均APE 為30%，明顯優于其他算法，100 s 的APE已能達到10%以內；比較各算法相對PLE 方法的精度提升效果發現，CTLS 的精度提升效果最好、最穩定。仿真實驗結果表明，CTLS 與其他具有無偏性的算法相比，對角測量誤差與初始距離的抗干擾能力最強，且定位精度高，性能穩定，不易發生較大波動，同時兼具收斂速度快的優點。而且算法的運行速度與EKF 以及PLE 同一個量級，運行效率較高，因此特別適合用于遠距離場景、角測量誤差較大時的單站定位跟蹤問題。