基于非機動單站的目標(biāo)運動狀態(tài)無偏估計

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

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

1 引言

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

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

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

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

2.1 運動模型

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

圖1 東北天坐標(biāo)系Fig.1 North-East-Up（ENU）coordinate system

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

2.2 三維AOA 量測模型

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

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

不失一般性，假設(shè)采樣時間間隔為T，將式（1）代入式（3）中進(jìn)行代數(shù)變換：

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

經(jīng)過推導(dǎo)可以得到噪聲向量εk為：

2.3 可觀性分析

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

定義偽觀測矩陣：

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

令M=HTθ Hθ，則有：

經(jīng)推導(dǎo)得到M的行列式值為：

3 三維定位算法

3.1 偽線性最小二乘法

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

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

根據(jù)先驗速率vc，進(jìn)一步可以得到目標(biāo)的初始距離、位置矢量、目標(biāo)與觀測平臺距離：

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

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

3.2 約束總體最小二乘法

3.2.1 算法原理

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

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

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

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

更進(jìn)一步，將損失函數(shù)J展開，并對等式兩邊取期望，得到：

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

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

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

拉格朗日極值法可處理帶約束條件的極值問題，構(gòu)造函數(shù)：

對方程求導(dǎo)并令其導(dǎo)數(shù)為0，得：

于是，v*，λ*即為矩陣（P，Ω）的最小廣義特征向量與最小特征值，改進(jìn)的漸進(jìn)無偏估計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展開：

根據(jù)隨機矩陣?yán)碚摚S機矩陣的極限也是它的漸進(jìn)期望，于是當(dāng)測量數(shù)據(jù)量足夠大時有：

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

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

4 仿真分析

4.1 性能指標(biāo)

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

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

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

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

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

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

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

式中：Index表示待評價指標(biāo)參數(shù)，下標(biāo)new，old分別表示初始方法與改進(jìn)方法。

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

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

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

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

圖2 觀測站-目標(biāo)運動態(tài)勢Fig.2 Observation station -target motion position

圖3 比較了各算法在100 s 時，不同角度測量標(biāo)準(zhǔn)差條件下的RDE 與APE。表2 與表3 比較了不同算法相比PLE 的精度。從圖3（a）可以看出，當(dāng)測量噪聲為0.1°時，各算法均能有效收斂到10%以內(nèi)的相對距離誤差。其中，約束總體最小二乘法與IV，WIV 方法誤差接近。隨著測量噪聲逐漸增大到0.2°，各算法性能受到不同程度的影響。PLE 與CTLS 的估計誤差均保持平穩(wěn)上升，但PLE 的收斂精度明顯低于CTLS；EKF，IV，WIV 均出現(xiàn)較大的波動，性能不夠穩(wěn)定，EKF 的估計精度反而有所提升，這可能與初值P的選擇是否恰當(dāng)有關(guān)。當(dāng)測量噪聲增大至0.3°時，EKF 與PLE 方法在100 s 時的相對距離誤差達(dá)到30% 以上，此時無法有效工作，而IV 與WIV 則出現(xiàn)發(fā)散傾向，RDE 與APE 均高于其他算法。其中，WIV 的絕對位置誤差超過100 km，但是CTLS 在100 s 時的RDE 仍能達(dá)到10%以內(nèi)。根據(jù)表2 與表3，在100 s 時CTLS 的提升效果最為穩(wěn)定與明顯，在角測量標(biāo)準(zhǔn)差為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%，精度的提升效果隨著角測量標(biāo)準(zhǔn)差的增大逐步提升，算法的優(yōu)勢逐漸明顯。

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

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

圖3 不同角度測量標(biāo)準(zhǔn)差下的算法性能比較Fig.3 Comparison of algorithm performance at different angles of measurement standard deviation

為了進(jìn)一步精確評估算法的精度與收斂性能，計算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 不同初始距離條件下的仿真對比

在速度一定的條件下，目標(biāo)起始點會影響角度觀測角的變化，因此進(jìn)一步對初始距離的影響進(jìn)行探究。取單位起始點位置矢量為[100，100，5]T（km），設(shè)置不同起始點，分別為起始點位置矢量的0.5～2 倍，間隔為0.25 倍，其他條件相同，角測量標(biāo)準(zhǔn)差均設(shè)置為0.1°。圖4 為算法在不同初始距離條件下的RDE 與APE。由圖4 可知，當(dāng)取0.5 倍單位初始距離時，除了EKF外，其他算法的相對距離誤差與絕對距離誤差均十分接近；當(dāng)初始距離逐漸提高至一倍初始距離時，各算法性能逐漸下降并有所區(qū)別，100 s 時按RDE 排序為PLE＞EKF＞IV＞W(wǎng)IV＞CTLS，APE 為EKF＞PLE＞IV＞W(wǎng)IV＞CTLS，各算法100 s 時均能達(dá)到10%以內(nèi)的相對距離誤差；當(dāng)初始距離增加到2 倍初始距離時，各種算法誤差大大提高，其中IV，WIV 出現(xiàn)較大波動，當(dāng)初始距離為1.75～2 倍單位初始距離之間時，二者的APE 均大于PLE 與EKF，已超過100 km，對目標(biāo)的運動信息估計作用已失效。而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%，精度提升效果隨著初始距離的增大逐步提升，且提升效果始終最佳，性能穩(wěn)定。

表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 內(nèi)，不同初始距離條件下的時均RDE 與APE 如表7 所示。通過比較可以發(fā)現(xiàn)，當(dāng)初始距離較近，取到0.5 倍單位初距時，各算法的性能接近，IV，WIV 與CTLS 能夠收斂到時均RDE 為1%，APE 為1 km 以內(nèi)。而初始距離增加到單位初始距離的1 倍與2 倍時，其他算法均在CTLS 性能參數(shù)的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 相當(dāng)，而WIV 的運行效率遠(yuǎn)低于其他算法，且實驗過程中其定位并不明顯優(yōu)于IV，與IV 一樣受測量噪聲與初始距離誤差的影響較大，具有性能不穩(wěn)定的缺點。

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

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

5 結(jié)論

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