曲志昱 王超然 孫 萌

(哈爾濱工程大學 哈爾濱 150001)

(先進船舶通信與信息技術工業(yè)和信息化部重點實驗室 哈爾濱 150001)

1 引言

星載無源定位系統(tǒng)因其具有探測距離遠、監(jiān)控范圍廣、不受地域國界影響、隱蔽性能好等優(yōu)勢，在電子對抗、偵察監(jiān)視等領域[1–3]發(fā)揮著越來越重要的作用。近年來，隨著星載電子偵察技術的發(fā)展，其定位體制日趨多樣化[4]，定位目標也從最初的地面固定目標延伸至空域運動目標[5]。傳統(tǒng)星載時頻差無源定位算法應用于高速運動目標時，常忽略目標的運動速度，這會導致頻差定位結(jié)果出現(xiàn)較大偏差[6]，另外，高速運動的目標一般都位于高空而非地表，但傳統(tǒng)星載定位普遍利用假設高程的球面模型或橢球面模型，即通過假設目標高程為0的方式將地球表面視作定位面進而實現(xiàn)目標定位[7]。因此，對于搭載于高速運動平臺(如戰(zhàn)機、導彈等)上的輻射源目標而言，這些算法將不再適用。

文獻[8]提出一種基于高程值估計的3星時頻差運動輻射源定位方法，使用解析法實現(xiàn)了動目標的測速和定位，但精度較低，且隨著輻射源高程的增加，定位、測速的誤差會越來越大，文獻[9]研究一種基于擴展卡爾曼濾波(Extended Kalman Filter,EKF)的時頻差聯(lián)合定位跟蹤方法，改善了高速目標的速度帶來的測量誤差大的問題，但仍將目標高程假設為已知量參與跟蹤計算；文獻[10]提出一種基于時差頻差的多站多次定位跟蹤濾波方法，但該算法應用于星載3維定位時至少需要4個觀測站，這將大大提升探測成本和衛(wèi)星組網(wǎng)難度；文獻[11]給出一種高精度測向輔助的3星時差定位算法，該算法通過引入方位角信息，使定位目標不再局限于地面輻射源，但由于缺乏對目標速度的直接觀測，因此算法對動目標定位跟蹤的效果較差。針對上述問題，本文在3星時頻差定位系統(tǒng)基礎上，提出了一種利用主星1維測向信息的改進迭代擴展卡爾曼濾波(Iterative Extended Kalman Filter,IEKF)融合跟蹤算法，算法在EKF方法中加入迭代，后采用LM(Levenberg-Marquardt)方法[12]對迭代過程的狀態(tài)更新進行優(yōu)化，并給出了迭代終止條件。仿真結(jié)果表明，算法可以實現(xiàn)對未知高程運動目標的跟蹤、定位和測速，且性能優(yōu)于傳統(tǒng)EKF和IEKF算法。

2 時差頻差測向融合跟蹤模型

3星融合體制定位系統(tǒng)由1顆主星和2顆輔星構成，具體形式如圖1所示，其中，α為地球自轉(zhuǎn)時間角，f為真近心角(真近點角)，Ω為升交點赤經(jīng)，w為近地點幅角，i為軌道傾角。系統(tǒng)的主星可通過干涉儀測得目標輻射源到達主星的方位角信息(Direction Of Arrival,DOA)，并可以測得輻射源信號到達該衛(wèi)星與兩個輔星的到達時間差(Time Difference Of Arrival,TDOA)與到達頻率差(Frequency Difference Of Arrival,FDOA)信息。

圖1 TDOA-FDOA-DOA跟蹤定位系統(tǒng)示意圖

設置衛(wèi)星軌道為圓軌道，定義星體坐標系原點為衛(wèi)星質(zhì)心，X軸指向衛(wèi)星前進方向，Z軸指向地球質(zhì)心，Y軸與另外兩軸成右手關系。為簡化討論，假設衛(wèi)星姿態(tài)角為0，天線陣面安裝角也為0，即主星的天線陣面坐標系與定義的主星星體坐標系重合。

衛(wèi)星定位中，描述輻射源坐標時常采用WGS-84系(ECEF坐標系)，但參數(shù)的測量卻是在衛(wèi)星的天線陣面上進行的，因此為方便分析，應通過坐標系轉(zhuǎn)換將輻射源和兩顆輔星的位置轉(zhuǎn)換至主星的天線陣面坐標系下(本文中即主星星體坐標系)。ECEF系至主星的星體坐標系間轉(zhuǎn)換方式為

其中，(xe,ye,ze)為輻射源或衛(wèi)星在ECEF坐標系下的坐標，(xt,yt,zt)為其在主星的星體坐標系下的坐標，(xe0,ye0,ze0)為主星在ECEF坐標系下的坐標，而(xt0,yt0,zt0)為主星在其本體坐標系下的坐標(0,0,0)，因此可在計算中省略，R1為ECEF坐標系到ECI坐標系的轉(zhuǎn)換矩陣，忽略對其影響較小的章動和歲差，可由仿真時刻儒勒日計算出地球自轉(zhuǎn)時間角α得 到；R2為ECI坐標系到主星的星體坐標系的轉(zhuǎn)換矩陣，由軌道根數(shù)f,Ω,w,i得到。由式(1)可知，星體坐標系至ECEF坐標系的轉(zhuǎn)換關系為

設置輻射源為一高程未知的等高程勻速運動目標，在k時刻，其狀態(tài)方程為

其待定的位置和速度分別為(xk,yk,zk)和(vxk,vyk,vzk)，則該時刻輻射源運動狀態(tài)矢量Xk=[xk,yk,zk,vxk,vyk,vzk]T,A為狀態(tài)轉(zhuǎn)移矩陣，定義為

其中，Δt為采樣周期，wk-1為相應的過程噪聲，其對應協(xié)方差為Qk-1。

方程組式(7)前3式為測量方程，其中，c為光速，tn,k和fn,k分別表示輔星n在k時刻與主星的時差(TDOA)、頻差(FDOA)觀測量；θk為主星的1維干涉儀測得的k時刻到達角觀測量，為簡化討論，這里假設該角為來波方向與ECEF系的X軸方向的夾角，eX與eY為主星的星體系坐標軸的單位矢量，eX=[1,0,0]T,eY=[0,1,0]T。在跟蹤過程中，輻射源始終做等高程勻速運動，因此其速度和位置矢量應滿足正交性，由此得出方程組第4式。dtk,dfk,dθk,dOk為相互獨立的高斯白噪聲，其中，dtk,dfk,dθk分別代表時差、頻差和角度的觀測誤差，dOk為目標位置與速度的正交約束誤差。

3 改進IEKF跟蹤算法

3.1 濾波初值計算方法

工程上，擴展卡爾曼濾波器(EKF)是一種應用廣泛的非線性濾波器，利用EKF對輻射源目標的位置和速度進行濾波跟蹤是解決此類問題的常用方式，對于EKF而言，一個理想的濾波初值可以帶來更小的截斷誤差和更快的收斂速度，初始狀態(tài)選擇不當會導致濾波器發(fā)散，因此濾波初值選取的優(yōu)劣對濾波結(jié)果有著直接影響。由定位方程組式(7)可知，時差和角度觀測方程與目標速度無關，不會由于運動目標速度引入誤差。因此，在研究運動目標時，首先融合TDOA和DOA兩種定位信息來得到濾波初值。

濾波初值的計算在主星的星體坐標系下進行，首先設主星星體坐標系下輻射源初始位置為ub,t=[xb,t,yb,t,zb,t]，3 顆衛(wèi)星初始位置為sb,n=[xb,n,yb,n,zb,n]，其中n=0時表示主星初始位置，其在星體坐標系中的坐標為(0,0,0),n=1,2時分別表示兩顆輔星初始位置。設主星與目標輻射源的距離為r0，輔星與目標輻射源的距離為rn,n=1,2，則有

在主星的星體系下，有如式(9)測量方程

其中，tn,0分別表示輔星n在初始時刻與主星的時差(TDOA)觀測量；θ0為主星的1維干涉儀測得的初始時刻到達角觀測量(來波方向與主星星體坐標系X軸夾角)，dt0,dθ0分別代表初始時刻時差和角度的觀測誤差。

由方程組式(9)中的第1式可得

聯(lián)立定位方程組式(9)第2式與式(10)可得

其中，n=1或2，進而有

將式(12)和方程組式(9)第2式代入式(8)，可得到關于r0的一元二次方程

由式(13)解得r0的兩個根，去除鏡像模糊值后，將得到的解代入式(12)和定位方程組式(9)第2式，解出目標在主星星體系下的初始定位坐標，將該位置坐標由式(2)坐標轉(zhuǎn)換至ECEF坐標系下，即得到濾波初值u0。對于協(xié)方差矩陣，分別對式(2)及定位方程組式(9)的兩式求解全微分，可得到初始位置估計的協(xié)方差矩陣

其中，d iag為以相應元素為對角線元素的對角陣，σt表示時差的測量精度，σθ表示方位角的測量精度，矩陣M為

目標運動狀態(tài)此時未知，故將初始速度置為0，可得濾波初始值為(,0,0,0)，初始協(xié)方差矩陣為

其中，vmax為目標速度可能達到的最大值。

3.2 基于融合定位體制的EKF算法

卡爾曼濾波是一種可以對線性系統(tǒng)狀態(tài)進行最優(yōu)估計的算法，但融合體制3星定位系統(tǒng)是一個非線性系統(tǒng)[13]，為解決這一問題，本文首先給出了融合定位系統(tǒng)的擴展卡爾曼濾波(Extended Kalman Filter,EKF)方法，其基本思想是將非線性系統(tǒng)線性化，對非線性函數(shù)的泰勒展開式進行1階線性化截斷，忽略其高階項，從而將非線性問題轉(zhuǎn)化為線性，然后進行卡爾曼濾波。

首先設系統(tǒng)狀態(tài)方程和觀測方程分別為

其中，A是狀態(tài)轉(zhuǎn)移矩陣，狀態(tài)轉(zhuǎn)換函數(shù)h指代k個時間采樣點下狀態(tài)量Xk與觀測量Zk之間的函數(shù)關系，wk-1為過程噪聲，Rk為測量噪聲。

EKF過程如下：

目標狀態(tài)預測

預測協(xié)方差矩陣

卡爾曼增益矩陣計算

更新目標狀態(tài)

更新協(xié)方差矩陣

3.3 改進的迭代卡爾曼濾波(IEKF)算法

EKF算法應用在非線性系統(tǒng)時，會產(chǎn)生2次項以上的截斷誤差，尤其在線性模型誤差較大時，可能導致濾波結(jié)果誤差較大。IEKF算法在EKF的基礎上，在測量更新階段加入多次迭代，重復對觀測方程進行線性化，最終得到更為準確的濾波結(jié)果，由文獻[14]可定義其代價函數(shù)

IEKF的迭代過程可視為利用高斯-牛頓方法迭代求解式(24)的極小值點，但由于線性化帶來的截斷誤差等因素，往往只能求出其非零極值點，即等價于將高斯-牛頓方法應用于小殘差問題的迭代優(yōu)化。對于此類問題，IEKF算法存在著不穩(wěn)定性，主要表現(xiàn)為對狀態(tài)的觀測更新不能保障誤差的一致減小，對協(xié)方差陣估計值比真實值偏低，進而影響對觀測信息的有效利用。為解決這一問題，本文算法利用L-M方法調(diào)整預測協(xié)方差陣，以保證算法具有全局收斂性，算法核心是在每次迭代過程中使用阻尼因子μi對協(xié)方差預測矩陣進行修正，以修正后的協(xié)方差矩陣進行觀測的迭代更新[15–17]。具體過程如下：

目標狀態(tài)預測

預測協(xié)方差矩陣

迭代開始

L-M方法修正協(xié)方差矩陣

卡爾曼增益矩陣計算

更新目標狀態(tài)

更新協(xié)方差矩陣

當代價函數(shù)隨迭代次數(shù)增加而小于一定限度時，即認為收斂已至極值點附近，迭代終止條件可寫為

其中，ε是預先設置的門限值。

整個過程的具體流程圖如圖2所示。

圖2 改進IEKF流程示意圖

4 仿真結(jié)果與分析

利用融合定位體制對輻射源進行定位時，衛(wèi)星的編隊方式將對定位結(jié)果產(chǎn)生影響[18]，設3顆衛(wèi)星軌道高度均為800 km，星間距100 km；目標初始位置為(28.73°E,8.61°N)，高程為10 km，目標運動狀態(tài)方程如式(5)所示，做速度為400 m/s勻速直線運動。觀測時長為100 s，觀測時間間隔為1 s，時差測量誤差為20 ns，多普勒頻率測量誤差為3 Hz，測向誤差為0.25°。圖3、圖4分別給出了3星同軌和2星同軌1星異軌時的融合定位體制的GDOP圖，圖中的誤差單位均為km。由圖中可見在融合定位體制中，同軌3星的定位誤差情況優(yōu)于2星同軌1星異軌情況，這是由于本算法的主要誤差來自測角誤差，而3星同軌的編隊方式在主星星體系下聯(lián)立時差方程求解輻射源位置時所受測角誤差影響相對三角構型要小得多，因此雖然對于時差測量來說等腰三角構型的誤差情況要優(yōu)于同軌3星構型[19]，但綜合考慮總體誤差，3星同軌編隊方式依然更優(yōu)。

圖3 3星同軌GDOP分布

圖4 2星同軌1星異軌GDOP分布

綜合考慮衛(wèi)星傳感器的測量參數(shù)誤差和覆蓋范圍情況，本文的算法更適合3星同軌編隊的應用場景。圖5—圖8仿真分析了3星TDOA/FDOA系統(tǒng)的EKF跟蹤方法、3星TDOA/DOA系統(tǒng)的EKF跟蹤方法以及本文的融合跟蹤系統(tǒng)EKF方法對于未知高程的運動目標位置、速度的跟蹤性能，并將融合跟蹤系統(tǒng)的EKF方法與改進IEKF方法進行了比較。

圖5 目標定位跟蹤結(jié)果

圖6 目標高程跟蹤結(jié)果

圖7 目標位置估計誤差

圖8 目標速度估計誤差

從圖5、圖6可以看出，在未知高程的情況下，本文所提融合定位算法可以實現(xiàn)運動輻射源目標的跟蹤定位。從圖7可以看出，本文融合定位系統(tǒng)的定位精度優(yōu)于3星TDOA/DOA定位系統(tǒng)，且改進IEKF的融合跟蹤方法定位精度略高于EKF算法。而由于對目標高程估計誤差大，3星TDOA/FDOA定位性能明顯不如其他算法，且其誤差會隨著目標高程增加而變得更大，4種算法收斂速度都較為緩慢，這是時差或頻差的測量方程相對于目標位置的強非線性所導致的，在跟蹤時加入迭代能一定程度上緩解這一問題。由圖8可以看出，由于缺乏對速度參數(shù)的直接觀測，3星TDOA/DOA定位系統(tǒng)對于速度的跟蹤誤差較大，而加入FDOA觀測量后，對速度的跟蹤性能明顯提升。而由于定位初值偏差較大，3星TDOA/FDOA定位系統(tǒng)收斂速度明顯不如本文提出的融合定位系統(tǒng)。

5 結(jié)束語

本文針對未知高程的運動目標跟蹤問題提出了一種基于3星TDOA/FDOA/DOA的融合定位方法，并針對傳統(tǒng)EKF方法定位跟蹤精度低、收斂慢的問題，提出了改進的IEKF融合跟蹤方法。首先利用3星時差信息和主星1維測向信息得到濾波初值及其協(xié)方差陣，再利用改進的IEKF結(jié)合3星頻差信息進行濾波跟蹤。仿真結(jié)果表明，在未知高程信息的情況下，算法可以實現(xiàn)運動目標的高精度位置跟蹤，并能獲得目標的高精度速度估計。