基于MCIS技術(shù)的相位差變化率單站無源定位

邢懷璽,張宇暉,陳游,*,周一鵬,何文波

1. 空軍工程大學(xué) 航空工程學(xué)院,西安 710038 2. 空軍裝備部 駐成都第三軍事代表室,成都 610000

單站無源定位技術(shù)由于避免了多站定位的時間同步和數(shù)據(jù)融合等問題,具有良好的機(jī)動性和靈活性,對于提高飛機(jī)的突防能力、生存能力和隱蔽攻擊能力都有重要的現(xiàn)實意義,一直是無源定位領(lǐng)域的研究熱點[1-5]。最早的單站無源定位利用輻射源來波方向(DOA)實現(xiàn)交叉定位[6-8],但這一方法受限于來波方向的測量精度,利用長基線干涉儀(LBI)可以精確獲得非合作輻射源信號到達(dá)陣元間的相位差及相位差變化率。許耀偉和孫仲康[9]提出了利用干涉儀測量相位差變化率實現(xiàn)對目標(biāo)的測距,研究表明,采用相位差變化率定位較只測向定位大幅提高了定位精度和速度。相位差變化率是相位差對時間的導(dǎo)數(shù),因此不需要對相位差的解模糊運算,同時可以消除天線間通道幅度相位不一致帶來的相位差系統(tǒng)偏差。王強等[10]提出了僅利用相位差變化率單站無源定位方法,只需要兩陣元構(gòu)成的長基線干涉儀測量來波的相位差變化率,便可實現(xiàn)對三維空間中固定輻射源的無源定位,定位系統(tǒng)設(shè)備相對簡單且成本低。

基于相位差變化率的單站無源定位實質(zhì)上是一個非線性最優(yōu)估計問題,一般很難直接獲得目標(biāo)位置的解析解。目前定位算法主要分為迭代法和搜索法兩類,其中迭代法應(yīng)用最多的是非線性最小二乘法和擴(kuò)展卡爾曼濾波(EKF)以及性能更好的容積卡爾曼濾波(CKF)、不敏卡爾曼濾波(UKF)算法等。文獻(xiàn)[11]采用非線性最小二乘估計的方法,獲取多個時刻的相位差變化率迭代遞推目標(biāo)的位置；單月暉等[12]將EKF和修正協(xié)方差擴(kuò)展卡爾曼濾波(MVEKF)方法用于測相位差變化率的無源定位中,驗證了MVEKF方法有良好的收斂性能；郭福成等[13]首次將MVEKF應(yīng)用于無源定位跟蹤中,無需計算測量方程的修正函數(shù),適用范圍更廣；文獻(xiàn)[14]利用一種基于改進(jìn)的變尺度最小斜度不敏卡爾曼濾波(SMSS-UKF)估計目標(biāo)位置,在低信噪比條件下具有較好的適應(yīng)性,但這一類方法非常依賴于初始化的結(jié)果,在狀態(tài)估計誤差較大的情況下,尤其是濾波初期,受測量噪聲影響較大,估計過程中協(xié)方差易出現(xiàn)病態(tài), 致使濾波結(jié)果不穩(wěn)定。此外,最大似然估計(ML)[15]是一種漸進(jìn)無偏的估計算法,被廣泛用于定位問題中,可以達(dá)到克拉美羅下界(CRLB)。迭代方法也可以用來解決ML問題,常用的方法有牛頓法、共軛梯度法等。這些方法對初始估計同樣很敏感,在高測量噪聲下可能會導(dǎo)致收斂至局部最優(yōu)或發(fā)散。另一個解決方案是將ML問題松弛為凸優(yōu)化問題,如半定規(guī)化松弛(SDP)[16]和二階錐規(guī)劃松弛(SOCP)[17]。搜索類算法優(yōu)勢在于當(dāng)缺少目標(biāo)位置先驗信息也可以通過構(gòu)建代價函數(shù),搜索代價函數(shù)極值點獲得目標(biāo)位置。文獻(xiàn)[10]在沒有目標(biāo)位置先驗信息的前提下,用ML代替最大后驗估計構(gòu)造似然函數(shù),將目標(biāo)位置估計轉(zhuǎn)化為一個非線性非凸優(yōu)化問題,然后利用網(wǎng)格搜索法得到全局的最優(yōu)解,算法穩(wěn)定性較高。然而,要想實現(xiàn)較高的定位精度,需要足夠大的網(wǎng)格密度,無法避免算法的運算量隨著網(wǎng)格分辨率增加呈指數(shù)增長這一固有矛盾,同時難以實現(xiàn)工程化的實時應(yīng)用。

近年來,蒙特卡洛(MC)技術(shù)[18-19]廣泛應(yīng)用于信號處理領(lǐng)域,其中重要性抽樣是計算多維積分強有力的工具。文獻(xiàn)[20]提出一種基于重要性抽樣的最大似然估計方法,利用重要性抽樣(IS)技術(shù)實現(xiàn)了DOA的ML估計,大大降低了ML估計的計算量；文獻(xiàn)[21]采用馬爾科夫鏈蒙特卡洛(MCMC)方法求解似然函數(shù)的全局極大值,得到時差-頻差聯(lián)合估計,估計精度優(yōu)于傳統(tǒng)的互模糊函數(shù)(CAF)算法,而且克服了迭代算法的初值依賴和收斂問題,同時又避免了網(wǎng)格搜索大量的計算過程；文獻(xiàn)[22]利用一種采用蒙特卡洛重要性抽樣(MCIS)方法近似多維積分的方法解決到達(dá)時間差(TDOA)的源定位問題,通過概率形式選取合適的樣本,以樣本均值代替目標(biāo)位置估計,取得理想的效果。

基于此,本文提出一種基于相位差變化率的高精度,低復(fù)雜度的單站無源定位方法。首先給出高斯測量噪聲下的目標(biāo)輻射源位置的最大似然估計,然后引入Pincus定理導(dǎo)出多維積分下的全局最優(yōu)解,最后通過MCIS技術(shù)可以得到一個近似的全局解。這一過程需要目標(biāo)輻射源位置粗略的初始估計值。一般的迭代方法本身也需要初值估計,而最小二乘算法能夠根據(jù)定位方程直接解算目標(biāo)位置。所以本文利用加權(quán)最小二乘(WLS)算法[23-24]獲得目標(biāo)位置封閉形式解作為本文定位方法的初值估計。本文將相位差變化率測量方程在目標(biāo)位置估計初值處一階泰勒級數(shù)展開得到線性化方程,構(gòu)造一個高斯分布的PDF作為重要性函數(shù)。以重要性函數(shù)的概率分布抽取足夠多的樣本,根據(jù)大數(shù)定律,以樣本估計的期望值作為全局最優(yōu)解,從而得到輻射源位置。結(jié)果表明,MCIS方法能夠提供最優(yōu)的性能,不需要大量計算,在一定的噪聲水平下定位誤差逼近CRLB,并且對初始估計誤差的敏感性較低。

1 基于相位差變化率的定位原理

1.1 定位模型

假設(shè)單個觀測平臺上垂直機(jī)身安裝一組二陣元的一維單基線干涉儀,在一段觀測時間內(nèi),干涉儀多次測量來波方向的相位差變化率,定位場景如圖1所示。

空中觀測平臺的位置速度和姿態(tài)角信息可以從慣導(dǎo)系統(tǒng)或GPS中獲取。假設(shè)固定目標(biāo)輻射源的位置狀態(tài)為xT=[xT,yT,zT]T,在t1～tn

圖1 單基線測相位差變化率的定位模型示意圖Fig.1 Schematic diagram of positioning model for measuring change rate of phase difference with single baseline

(1)

式中:κ=2πdf/c,f為輻射源信號頻率,d為基線長度,c為光速;Lk=rk/rk。相應(yīng)的tk時刻的相位差變化率為

(2)

(3)

(4)

將K次測量的相位差變化率用向量形式表示

(5)

1.2 僅測相位差變化率的三維可觀測性分析

選機(jī)腹下機(jī)身軸和機(jī)翼軸的交點O作為參考原點,沿著機(jī)身軸方向并指向機(jī)頭的方向為參考X軸,OY軸與OX軸在同一水平面上且指向左側(cè)機(jī)翼方向,按照右手關(guān)系作OZ軸垂直于OXY平面且指向上方,建立三維直角坐標(biāo)系O-XYZ,如圖2所示。

假設(shè)觀測平臺水平直線勻速運動,對于靜止目標(biāo),記k時刻,以O(shè)Y方向為基準(zhǔn),機(jī)載觀測平臺和輻射源之間的方位角為βk,以平面O-XY為基準(zhǔn),機(jī)載觀測平臺和輻射源之間的俯仰角為εk;觀測平臺到輻射源徑向距離在水平面的投影距離為rpk。

所以k時刻相位差用方位角和俯仰角表示為

(6)

相位差變化率進(jìn)一步表示為

(7)

(8)

(9)

圖2 觀測平臺與輻射源相對位置示意圖Fig.2 Schematic diagram of relative position between observing station and radiation source

(10)

(11)

其中：rpk=(xT-xok)2+(yT-yok)2。

將式(8)～式(11)代入式(7)可得

f(xk,yk,zk)

(12)

2 重要性抽樣定位方法

無源定位實質(zhì)上是一個非線性最優(yōu)估計問題,在基于相位差變化率的單站無源定位中,相位差變化率與目標(biāo)位置存在非線性的約束關(guān)系,通過多個時刻觀測值獲取目標(biāo)位置狀態(tài)信息。目前的研究中以非線性最小二乘(NLS)、EKF等其他改進(jìn)的非線性跟蹤濾波算法為主體,但這類算法的性能嚴(yán)重依賴初始狀態(tài)估計,初始階段定位誤差較大會出現(xiàn)發(fā)散的現(xiàn)象,同時受測量噪聲影響大,導(dǎo)致濾波結(jié)果不穩(wěn)定。除此之外,也可以通過構(gòu)建代價函數(shù),通過網(wǎng)格搜索代價函數(shù)最小值獲取目標(biāo)位置信息,ML就是其中之一。

目標(biāo)輻射源位置ML可表示為矢量形式

(13)

式中：Q=E(eeT)為測量誤差的協(xié)方差。

2.1 ML估計的全局解

ML問題是一種非線性、多維的極值優(yōu)化問題,需要全局極值的多維搜索,因此其計算量巨大。因此在實際的操作中往往將其降維處理,這里將ML問題轉(zhuǎn)化為一維優(yōu)化問題,這里介紹Pincus定理。

(14)

(15)

在實際的測量中,需確定有足夠大的范圍保證x在范圍邊界以內(nèi)。重要性抽樣的關(guān)鍵在于重要性函數(shù)的選取,合適的重要性函數(shù)可以減少計算量,大大提高算法的計算速度,較差的重要性函數(shù)則會影響仿真時間,同樣影響優(yōu)化結(jié)果的質(zhì)量。這里,定義-λf(x)的歸一化函數(shù)為

(16)

(17)

利用式(17)計算x的最佳估計值,實際上是求解一個多維多重積分問題,直接計算非常復(fù)雜且困難。由大數(shù)定律可知,當(dāng)樣本數(shù)量足夠多時,樣本均值的方差越小,樣本均值逼近真實積分值。因此針對上述問題可以通過p(x)生成多個x的樣本,用樣本均值代替式(17)的復(fù)雜積分。

2.2 蒙特卡洛重要性抽樣技術(shù)

蒙特卡洛重要性抽樣技術(shù)是解決多維積分的強有力工具。蒙特卡洛方法建立簡單的概率統(tǒng)計模型可以求解某個隨機(jī)變量的數(shù)字特征,通過合適的抽樣方法抽取樣本,求取樣本的統(tǒng)計量,從而得到參數(shù)估計。顯而易見,蒙特卡洛方法的關(guān)鍵在于抽樣環(huán)節(jié),抽樣過程中并不能對服從任意概率分布的隨機(jī)變量直接進(jìn)行抽樣,然而重要性抽樣可以把難以直接抽樣的概率分布,轉(zhuǎn)換為從易于抽取樣本的重要函數(shù)進(jìn)行抽樣,簡單方便,且可以降低估計方差。其基本思想是通過一個相對簡單的分布函數(shù)生成一定數(shù)量的隨機(jī)數(shù),并將其加權(quán)平均來近似計算目標(biāo)分布函數(shù)的數(shù)學(xué)期望。下面以蒙特卡洛積分為例介紹重要性抽樣過程：

假設(shè)一個復(fù)雜積分

(18)

式中：p(θ)為θ的概率密度函數(shù)。如果q(θ)為積分域上的概率密度函數(shù),那么

(19)

(20)

根據(jù)大數(shù)定律,式(20)的積分被近似為

(21)

2.3 測相位差變化率定位的重要性函數(shù)選擇

對于重要性函數(shù)的選擇,這里希望能夠找到一種近似于式(16)的概率密度函數(shù)來直接計算積分?？紤]到相位差變化率的測量誤差服從高斯分布,同時包含目標(biāo)的位置信息,除此之外,高斯分布因為只需要均值和協(xié)方差,容易構(gòu)造,樣本的抽取比較簡單,所以本文選擇高斯分布的概率密度函數(shù)作為重要性函數(shù)。

由式(4)可得

(22)

(23)

(24)

為了簡化公式,令

所以式(23)可以改寫成

(25)

因為是關(guān)于xT的線性函數(shù),可以將K次觀測的線性方程用矩陣形式表示

e=HxT-F

(26)

式中：

根據(jù)式(13)得到目標(biāo)位置的最大似然估計近似為

(27)

因此,Pincus定理構(gòu)造重要性函數(shù)為

q(xT)=

(28)

λ的大小關(guān)系著抽樣區(qū)域的大小,為了構(gòu)造高斯分布的概率密度函數(shù),這里變換似然函數(shù)：

(HxT-F)TQ-1(HxT-F)=[xT-

(HTQ-1H)-1HTQ-1F]T·(HTQ-1H)·[xT-(HTQ-1H)-1HTQ-1F]+FTQ-1F-FTQ-1H(HTQ-1H)-1HTQ-1F=

(29)

其中：

(30)

R=(HTQ-1H)-1

(31)

ζ=FTQ-1F-FTQ-1H(HTQ-1H)-1HTQ-1F

(32)

這里因為ζ項與xT無關(guān),所以q(xT)分子中常數(shù)項exp(-λ1ζ)與分母對應(yīng)項相互抵消。重要性函數(shù)q(xT)可以改寫為

q(xT)=

(33)

所以式(33)進(jìn)一步簡化為

(34)

從式(34)中可以看出高斯分布的協(xié)方差與密切相關(guān),λ1的大小決定高斯分布重要性函數(shù)的幅度和寬度,當(dāng)λ1很大時,重要性函數(shù)很窄,容易生成大量樣本接近高斯分布均值,當(dāng)均值明顯偏離輻射源真實位置時,會使最終的定位結(jié)果產(chǎn)生偏差。參考文獻(xiàn)[19],這里取λ1=0.5。第4節(jié)中討論了參數(shù)λ1對MCIS算法定位精度的影響。

3 重要性抽樣定位

3.1 輻射源位置估計

在確定重要性函數(shù)之后,可以通過逆變換采樣的方法選擇樣本集,這里不作贅述。繼而很容易獲得目標(biāo)輻射源位置信息的最大似然估計,可以通過足夠多的樣本均值來近似計算,即

(35)

λf(xTm)]

(36)

顯然似然函數(shù)是一個多峰函數(shù),存在局部極值。為了明確區(qū)分全局最優(yōu)和局部極值,使樣本均值更加趨近于真實值,在這里對其指數(shù)歸一化處理得到標(biāo)準(zhǔn)化的重要性系數(shù)。

(37)

從式(30)中可以看出,參數(shù)λ的選擇對于標(biāo)準(zhǔn)化重要性系數(shù)ω′(xTm)估計的準(zhǔn)確度有著重要影響。根據(jù)Pincus定理,當(dāng)λ→∞時,式(16)的似然函數(shù)呈現(xiàn)n維的Dirac-delta函數(shù)形狀,使全局最大值和其他極值差異明顯。實際上,對于本文的定位模型,輻射源真實位置附近區(qū)域只存在單一的極值點,所以λ的取值不必太大,只需要保證滿足相應(yīng)的計算要求。經(jīng)過實驗分析,當(dāng)λ=0.2時已經(jīng)能夠滿足計算要求。第4節(jié)也給出了參數(shù)λ對MCIS算法估計性能的影響。

3.2 算法框架

步驟1通過加權(quán)最小二乘(WLS)算法獲取初始目標(biāo)位置估計,計算系數(shù)矩陣H和F。

步驟3設(shè)置樣本數(shù)量M,利用逆變換采樣定理(Probability Integral Transformation Theorem)生成樣本,生成M個獨立且服從[0,1]均勻分布的隨機(jī)向量um,對于每一個um,在q(xT)找到最接近的值,對應(yīng)的xTm為生成的樣本。

步驟4計算標(biāo)準(zhǔn)化重要性系數(shù)ω′(xTm)并根據(jù)式(38)估計目標(biāo)位置

(38)

3.3 定位誤差的CRLB

因為CRLB清楚地刻畫了參數(shù)測量誤差下的定位誤差及其分布,所以推導(dǎo)相位差變化率定位誤差的CRLB作為定位方法性能的評價指標(biāo)。根據(jù)CRLB的定義：

(39)

(40)

因此,可得

(41)

4 仿真分析

4.1 實驗1：MCIS算法估計輻射源位置

4.2 實驗2：MCIS參數(shù)對定位性能的影響

后面的討論將建立在實驗1模擬場景下進(jìn)行。同時,為了更好地衡量算法的定位精度,將目標(biāo)位置估計的均方根誤差(RMSE)作為指標(biāo)來衡量定位的精度。RMSE定義為

(42)

4.2.1 參數(shù)λ和λ1對定位精度的影響

改變參數(shù)λ和λ1,進(jìn)行200次蒙特卡洛實驗,討論MCIS的定位精度。圖4和圖5為不同參數(shù)λ和λ1設(shè)置下的定位誤差的RMSE。從圖4可以看出,當(dāng)λ≥0.2時定位結(jié)果趨于穩(wěn)定,說明無需達(dá)到理論上的無窮大,很容易滿足本文定位模型的要求。圖5中λ1在[0,5]區(qū)間內(nèi)改變,發(fā)現(xiàn)λ1的大小對定位精度幾乎沒有影響。也進(jìn)一步驗證了實驗1中λ和λ1的選取是合理有效的。

圖4 λ對MCIS定位精度的影響Fig.4 Influence of λ on MCIS positioning accuracy

圖5 λ1對MCIS定位精度的影響Fig.5 Influence of λ1 on MCIS positioning accuracy

4.2.2 抽樣數(shù)量對定位精度的影響

圖6為200次蒙特卡洛試驗下,不同樣本數(shù)量定位誤差的RMSE,顯然抽樣數(shù)量在10～100時始終能夠達(dá)到3 km以內(nèi)的定位誤差,說明樣本數(shù)量不是影響定位精度的主要因素。這是因為構(gòu)造的重要性采樣函數(shù)和似然函數(shù)相似度很高,生成樣本均服從某一個確定的高斯分布。即使生成不同數(shù)量樣本,抽樣值在目標(biāo)真實位置周圍的分布規(guī)律是相同的,所以生成樣本數(shù)量對定位結(jié)果影響不大。因此,可以選擇較少的樣本達(dá)到一個理想的定位結(jié)果,進(jìn)一步減小計算的復(fù)雜度,提高定位的收斂速度。

圖6 樣本數(shù)量對MCIS定位精度的影響Fig.6 Influence of sample number on MCIS positioning accuracy

4.3 實驗3：不同觀測條件下的定位性能

觀測平臺飛行速度保持300 m/s不變,圖7為不同的相位差變化率測量間隔下的定位相對誤差R隨觀測時長的變化曲線,圖中縱坐標(biāo)表示相對于觀測平臺和目標(biāo)歐氏距離的相對誤差?？梢?當(dāng)相位差變化率測量間隔一定,增加觀測時長可以提高定位精度,當(dāng)觀測時間T在30 s左右時能夠達(dá)到R=5%的定位誤差；此外,在相同的觀測時長下,測量間隔越小,也就是數(shù)據(jù)率越高時,收斂速度越快且定位精度越高。這是因為觀測時間的延長和數(shù)據(jù)率的提高使得系數(shù)矩陣維度增加,會在一定程度上修正IS估計帶來的誤差,但會導(dǎo)致算法的計算量增加的弊端。

圖8驗證了不同運動速率下定位的RMSE,觀測平臺作直線飛行。當(dāng)平臺速度為100 m/s時收斂速度較慢,觀測至60 s左右時才能達(dá)到R=10%的定位精度,觀測平臺速度越大,定位效果越好,尤其是當(dāng)速度高于400 m/s左右時,定位精度最高,而且觀測至10 s左右時就達(dá)到R=10%定位誤差,收斂速度較快。

圖7 不同測量間隔下定位性能Fig.7 Positioning performance at different measurement intervals

圖8 不同平臺速度下的定位性能Fig.8 Positioning performance at different platform speeds

對比分析延長觀測時間和增大觀測平臺的速度對算法的定位性能的影響,兩種情況都會延長觀測平臺的位移。所以,在本文的定位模型中,增大觀測平臺位移能夠有效提高算法的定位性能。這是由于單次測量相位差變化率對應(yīng)一個定位曲面,多次測量后只有經(jīng)過目標(biāo)輻射源的定位曲面才能交于一點,前后觀測平臺位移越大,使定位曲面的差異性越明顯,利于算法在迭代過程中快速排除虛假點,有效提高收斂速度和估計精度,定位結(jié)果的可信度更高。但是延長觀測時間與定位收斂速度之間相互沖突,因此可以考慮通過增大飛行速度提升定位性能。

從圖10中可以看出觀測平臺作機(jī)動時,MCIS算法的定位性能得到改善,定位精度更高且收斂速度更快。產(chǎn)生這一情況的原因是觀測平臺姿態(tài)發(fā)生變化的過程中,干涉儀基線矢量方向相對于目標(biāo)輻射源產(chǎn)生旋轉(zhuǎn),這種旋轉(zhuǎn)產(chǎn)生的“大”相位差變化率明顯改善了定位精度,這足以說明觀測平臺作機(jī)動轉(zhuǎn)彎時可觀測性更強。

圖9 直線運動和機(jī)動運動軌跡圖Fig.9 Trajectories of straight and maneuvering motion

圖10 直線運動和機(jī)動運動下的定位性能Fig.10 Positioning performance in straight and maneuvering motion

4.4 實驗4：測量噪聲對定位精度的影響

圖11 不同測量噪聲下幾種定位算法的RMSE對比Fig.11 Comparison of RMSE of several positioning algorithms with different measurement noise

4.5 實驗5：算法對初始估計的敏感性分析

圖12 初始估計誤差對算法定位性能的影響Fig.12 Influence of initial estimation error on algorithm positioning performance

4.6 實驗6：算法復(fù)雜度分析

本文除了將定位精度和收斂速度作為衡量本文算法的性能標(biāo)準(zhǔn)之外,也對算法的復(fù)雜度進(jìn)行了對比分析。表1對比了本文算法與EKF算法、NLS算法以及網(wǎng)格搜索(GS)的計算復(fù)雜度。K表示觀測次數(shù);D表示目標(biāo)狀態(tài)維度;M表示樣本數(shù)量;Ngrid表示網(wǎng)格密度。

分析表1中的各算法的計算復(fù)雜度,EKF算法和NLS算法復(fù)雜度與迭代次數(shù)成正比,采用網(wǎng)格搜索(GS)求解最大似然估計的計算復(fù)雜度與網(wǎng)格密度密切相關(guān),而且計算量隨估計量維度D的增加呈指數(shù)型爆炸增長。這里選擇樣本數(shù)量M=20,觀測次數(shù)K=100,網(wǎng)格密度Ngrid=25,計算復(fù)雜度比,本文提出的MCIS算法計算復(fù)雜度上明顯優(yōu)于其他幾種算法,GS算法計算量最大,耗時最長,可以通過減小網(wǎng)格密度或觀測次數(shù)降低其他幾種算法的復(fù)雜度,但必然會影響定位精度。

表1 算法復(fù)雜度對比Table 1 Complexity assessment of existing algorithms

5 結(jié) 論