一種改進(jìn)粒子濾波算法實(shí)現(xiàn)的多徑參數(shù)估計(jì)

國 強(qiáng)，劉雪萌，周 凱

(哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江 哈爾濱 150001)

隨著全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System，GNSS)的應(yīng)用領(lǐng)域日益廣泛，人們對其定位精度的要求越來越高。但由于衛(wèi)星信號在傳播過程中會經(jīng)障礙物產(chǎn)生反射[1]，致使接收到的信號往往由直達(dá)信號與多徑信號兩部分組成，進(jìn)而引起多徑誤差。多徑誤差具有位置上不相關(guān)及時(shí)間上不確定特性，使其很難通過差分技術(shù)消除，因此它是影響全球?qū)Ш叫l(wèi)星系統(tǒng)定位精度的主要誤差源之一[2]。

目前，對于全球?qū)Ш叫l(wèi)星系統(tǒng)中抗多徑干擾研究的方法大致可分為3類：一是基于接收機(jī)天線的設(shè)計(jì)[3-5]；二是基于接收機(jī)跟蹤環(huán)路的改進(jìn)[6-7]；三是基于數(shù)據(jù)處理的多徑抑制[8-9]。前兩者在實(shí)現(xiàn)過程中均面臨一定困難：基于天線端的設(shè)計(jì)，會增加天線成本，并不是一個經(jīng)濟(jì)的選擇；而基于接收機(jī)跟蹤環(huán)路的改進(jìn)則需要消耗較多的硬件資源。因此，基于數(shù)據(jù)處理的多徑抑制技術(shù)得到了廣泛應(yīng)用。數(shù)據(jù)處理技術(shù)的核心是信號參數(shù)估計(jì)，也正是筆者所研究的內(nèi)容。

文獻(xiàn)[10]采用一種基于高斯和近似的擴(kuò)展切片高斯混合濾波器實(shí)現(xiàn)直接信號時(shí)延和多徑相關(guān)輸出的參數(shù)估計(jì)，但其算法復(fù)雜度較高。隨著壓縮感知技術(shù)的興起，一種基于壓縮感知技術(shù)的GPS多徑估計(jì)方法被提出[11]，而最新研究顯示傳統(tǒng)的壓縮感知以稀疏結(jié)構(gòu)為先驗(yàn)信息來進(jìn)行信號恢復(fù)，使得壓縮感知的性能并不高。文獻(xiàn)[12]針對雷達(dá)系統(tǒng)參考信道中存在多徑干擾的問題，提出一種基于改進(jìn)恒模算法的參考通道多徑干擾抑制算法，該算法與傳統(tǒng)的盲均衡方法相比，具有較好的收斂性，多徑抑制性能有所提高。在信號特征提取領(lǐng)域，針對多徑信號產(chǎn)生的影響，一種基于時(shí)變AR模型的多徑信號抑制算法被提出[13]，該算法最終可精確地估計(jì)出多徑時(shí)延及幅度因子，進(jìn)而實(shí)現(xiàn)多徑信號的逐級抑制。由于大多數(shù)多徑抑制算法僅限于高斯噪聲，導(dǎo)致在非高斯環(huán)境下這些算法性能降低，為了解決這一問題有學(xué)者采用基于相關(guān)熵準(zhǔn)則的自適應(yīng)濾波多徑抑制算法[14]，提高多徑抑制性能。近年來，基于粒子濾波(Particle Filter，PF)算法的參數(shù)估計(jì)越來越普遍，王憲等[15]提出一種基于粒子濾波的GPS多徑估計(jì)，雖然解決了算法運(yùn)算量過大的問題，但是隨著時(shí)間的增加粒子退化、粒子多樣性降低現(xiàn)象仍然存在。因此基于PF的改進(jìn)算法大量涌現(xiàn)：在狀態(tài)和參數(shù)聯(lián)合估計(jì)過程中采用雙重粒子濾波算法[16]可避免粒子退化現(xiàn)象，但其計(jì)算量相對較大；文獻(xiàn)[17]中提出基于遺傳粒子濾波(Genetic Algorithm-Particle Filter，GA-PF)算法實(shí)現(xiàn)的多徑跟蹤，可以緩解粒子濾波算法中粒子多樣性降低的現(xiàn)象，但由于遺傳算法自身存在著過早收斂的問題，導(dǎo)致多徑跟蹤的效果并不理想；利用粒子群優(yōu)化算法(Particle Swarm Optimization，PSO)使得采樣粒子集向似然函數(shù)較大的區(qū)域移動[18]，該算法雖然減小了參數(shù)估計(jì)所需的粒子數(shù)目但卻沒有解決粒子退化問題；基于擴(kuò)展卡爾曼濾波(Extended Kalman Filter，EKF)改進(jìn)的粒子濾波[19]，利用EKF產(chǎn)生粒子濾波的重要性密度函數(shù)，然而，由于實(shí)際上很多系統(tǒng)的狀態(tài)方程是非線性的，使得EKF在線性化程中忽略了泰勒展開式的高階項(xiàng)[20]，產(chǎn)生截?cái)嗾`差；基于差分進(jìn)化(Differential Evolution，DE)改進(jìn)粒子濾波的多徑估計(jì)算法[21]，在差分進(jìn)化算法中縮放因子F與交叉因子CR的選取嚴(yán)重影響算法性能，使得差分進(jìn)化算法的尋優(yōu)能力和收斂速度相互矛盾，參數(shù)估計(jì)存在一定的誤差。

為克服上述問題，筆者提出一種無跡差分粒子濾波算法(Unscented Improved Differential Evolution-Particle Filter，UIDE-PF)。該算法首先將無跡卡爾曼波(Unscented Kalman Filter，UKF)算法“嵌入”到粒子濾波算法重要性采樣階段，得到合理的采樣分布，避免了粒子退化問題。然后，采用改進(jìn)差分進(jìn)化(Improved Differential Evolution，IDE)算法對重要性采樣后的粒子進(jìn)行自適應(yīng)變異、自適應(yīng)交叉以及選擇操作更新粒子，從而提高粒子多樣性，進(jìn)而減小參數(shù)估計(jì)結(jié)果的波動幅度及均方根誤差。

1 系統(tǒng)模型建立

全球?qū)Ш叫l(wèi)星系統(tǒng)在進(jìn)行定位時(shí)，其接收機(jī)接收到的信號通?？杀硎緸?/p>

r0(t)=A0c(t-τ0)cos(ω0t+φ) ，

(1)

其中，A0表示直達(dá)信號幅度；φ表示直達(dá)信號相位延遲；c(t)表示擴(kuò)頻碼序列；τ0表示直達(dá)信號時(shí)延；ω0表示載波相位；式中忽略導(dǎo)航電文。然而，實(shí)際接收信號還包含多徑及噪聲干擾，此時(shí)信號模型表示為

r1(kTs)=A0c(kTs-k0Ts)cos(ω0kTs+φ0)+

(2)

式(2)表示數(shù)字化后中頻信號，M表示多徑信號的路數(shù)；Ai、kiTs、φi分別表示第i路多徑信號的幅值、相對時(shí)延、相對相位延遲，其中，k表示第k個采樣點(diǎn)，Ts表示采樣周期；k0=τ0/Ts，表示τ0對應(yīng)采樣點(diǎn)數(shù)；ki=τi/Ts，表示τi對應(yīng)采樣點(diǎn)數(shù)；n(kTs)表示高斯白噪聲。對于靜態(tài)環(huán)境而言，一般是指接收機(jī)周圍環(huán)境變化相對較小，如平原、村落等地。此時(shí)接收機(jī)接收信號以雙徑模型為例[18]，具體形式表示如下：

r1(kTs)=A0c(kTs-k0Ts)cos(ω0kTs+φ0)+

A1c(kTs-k0Ts-k1Ts)cos(ω0kTs+φ0+φ1)+n(kTs) 。

(3)

(4)

Y=[R1(A0，Ai，η，ki，d1)，R2(A0，Ai，η，ki，d2)，…，RL(A0，Ai，η，ki，dL)] ，

(5)

其中，L表示相關(guān)器個數(shù)。

為了將無跡差分粒子濾波算法用于多徑參數(shù)估計(jì)，需要利用p(xk|xk-1)描述參數(shù)的演化特性[15]：

p(xi|xi-1，xi-2，…，x1)=p(xi|xi-1) 。

(6)

根據(jù)式(6)得到系統(tǒng)具體模型，其狀態(tài)方程與觀測方程如下：

(7)

其中，H表示狀態(tài)轉(zhuǎn)移矩陣；ω表示過程噪聲；R表示觀測矩陣；v表示觀測噪聲；Xn表示待估計(jì)參數(shù)。筆者的主要目的是用觀測向量Y估計(jì)出狀態(tài)向量Xn。

2 基于改進(jìn)算法的多徑參數(shù)估計(jì)

2.1 改進(jìn)差分進(jìn)化算法

差分進(jìn)化(Differential Evolution，DE)算法是一種智能優(yōu)化算法，標(biāo)準(zhǔn)差分進(jìn)化算法主要包括以下步驟[22]：

(1) 初始化種群。其表達(dá)式為

(8)

(2) 變異過程。在差分進(jìn)化算法中變異策略表達(dá)式為

(9)

(3) 交叉過程。該過程的主要目的是產(chǎn)生一個試驗(yàn)向量，其交叉過程表達(dá)式為

(10)

其中，CR一般在0～1之間。CR選取決定著目標(biāo)向量xj，i和變異向量vj，i在試驗(yàn)向量中所占比例的大小。

(11)

其中，f(·)表示適應(yīng)度函數(shù)。由于筆者的主要目的在于實(shí)現(xiàn)信號參數(shù)估計(jì)，因此可將參數(shù)估計(jì)過程中直達(dá)信號幅度、多徑信號幅度、本地估計(jì)偏差、多徑相對時(shí)延4種參數(shù)估計(jì)結(jié)果的絕對誤差之和作為適應(yīng)度函數(shù)，并保留適應(yīng)度函數(shù)較小的個體。

(5) 邊界條件處理。對于超出范圍的個體按如下表達(dá)式處理：

(12)

由于差分進(jìn)化算法中的F和CR往往是固定值，不能滿足算法各個階段的需求。因此筆者提出一種自適應(yīng)策略，使F和CR動態(tài)變化，改進(jìn)后的F和CR表達(dá)式如下：

(13)

其中，G表示最大進(jìn)化次數(shù)，g表示當(dāng)前進(jìn)化次數(shù)。由式(13)可以看出，在進(jìn)化前期g取值較小，此時(shí)獲得較大的F和較小的CR，有利于提高算法的全局搜索能力；在進(jìn)化后期g取值較大，此時(shí)獲得較小的F和較大的CR，有利于局部搜索并提高收斂速度。

2.2 多徑參數(shù)估計(jì)的實(shí)現(xiàn)過程

盡管粒子濾波算法可以應(yīng)用在任何非高斯非線性系統(tǒng)中，但其存在的粒子退化、多樣性降低問題仍不可避免。因此，針對該問題采用無跡差分粒子濾波算法實(shí)現(xiàn)信號參數(shù)估計(jì)，基本原理如下：

(1) 初始化。由于多徑信號的幅度小于直達(dá)信號的幅度，且歸一化幅度均在0～1之間；多徑信號延時(shí)一般在0～1.5碼片之間，即各參數(shù)的取值范圍如下：

(14)

其中，A0、A1分別表示直達(dá)信號、多徑信號幅度；τ表示多徑相對時(shí)延。按式(14)產(chǎn)生NP個粒子。

(2) 重要性采樣。利用UKF算法為每個粒子計(jì)算其均值及協(xié)方差，并不斷更新粒子。

① 利用無跡變換(Unscented Transformation，UT)變化計(jì)算2n+1個Sigma點(diǎn)集：

(15)

其中，P表示初始狀態(tài)協(xié)方差矩陣；i表示對應(yīng)的采樣點(diǎn)；n表示狀態(tài)參數(shù)的維數(shù)，文中n=4。

② 計(jì)算采樣點(diǎn)的相應(yīng)權(quán)值：

(16)

其中，m表示均值；c表示協(xié)方差；上標(biāo)表示對應(yīng)的采樣點(diǎn)；λ表示縮放比例；α和β表示常數(shù)。

③ 對2n+1個Sigma點(diǎn)集進(jìn)行下一步預(yù)測：

X(i)(k+1|k)=f[k，X(i)(k|k)]，i=1，…，2n+1 ，

(17)

其中，f(·)表示狀態(tài)函數(shù)。

④ 對步驟③中的Sigma點(diǎn)集加權(quán)求均值進(jìn)而得到狀態(tài)向量的一步預(yù)測結(jié)果，并計(jì)算協(xié)方差矩陣：

(18)

⑤ 根據(jù)式(18)再次使用無跡變換，得到新的Sigma點(diǎn)集如下：

(19)

⑥ 將式(19)代入式(7)中的觀測方程，得到一步預(yù)測觀測值：

Z(i)(k+1|k)=h[X(i)(k+1|k)] ，

(20)

其中，h表示觀測函數(shù)。

⑦ 將步驟⑥得到的觀測值通過加權(quán)平均，求得系統(tǒng)預(yù)測均值及協(xié)方差：

(21)

⑧ 計(jì)算卡爾曼增益：

(22)

⑨ 更新狀態(tài)向量及協(xié)方差：

(23)

(3)重采樣。采用IDE算法代替?zhèn)鹘y(tǒng)重采樣過程，得到最優(yōu)粒子并計(jì)算其歸一化權(quán)重：

(24)

其中，i表示樣本中第i個粒子；k表示迭代次數(shù)。

(4)根據(jù)步驟(1)～(3)可以得到采用無跡差分粒子濾波算法實(shí)現(xiàn)多徑參數(shù)估計(jì)的結(jié)果，其表達(dá)式如下：

(25)

根據(jù)以上對無跡差分粒子濾波算法原理的分析，給出其流程框圖，如圖1所示。

圖1 無跡差分粒子濾波算法流程框圖

圖2描述了在第500次迭代時(shí)采用無跡差分粒子濾波算法(UIDE-PF)與標(biāo)準(zhǔn)粒子濾波算法(PF)在參數(shù)估計(jì)過程中粒子權(quán)重的大小。從圖2(a)可以看出，采用標(biāo)準(zhǔn)粒子濾波算法在參數(shù)估計(jì)過程中，大權(quán)重的粒子不斷被復(fù)制，小權(quán)重的粒子逐漸被淘汰，導(dǎo)致粒子多樣性逐漸降低；而從圖2(b)可以看出，采用無跡差分粒子濾波算法在該過程中幾乎不存在粒子權(quán)重相等的現(xiàn)象。由此可見，無跡差分粒子濾波算法相比于標(biāo)準(zhǔn)粒子濾波算法可有效避免粒子多樣性降低問題。

(a) 采用PF算法得到的粒子權(quán)重

3 仿真結(jié)果及分析

為了驗(yàn)證無跡差分粒子濾波算法在參數(shù)估計(jì)過程中的有效性，將采用該算法得到的參數(shù)估計(jì)結(jié)果與無跡卡爾曼粒子濾波(Unscented Kalman Particle Filter，UPF)算法、PF算法及文獻(xiàn)[23]中的算法進(jìn)行對比分析。為了實(shí)驗(yàn)的嚴(yán)謹(jǐn)性，4種算法參數(shù)保持一致且與文獻(xiàn)[23]中算法參數(shù)設(shè)置相同：粒子數(shù)M=40；過程噪聲協(xié)方差Q=10-4eye(4)；各支路與即時(shí)碼之間的相關(guān)間距為±0.5 Tc，±0.3 Tc，±0.1 Tc；狀態(tài)轉(zhuǎn)移矩陣F=eye(4)；各參數(shù)初始狀態(tài)X0=[0.5 0.7 0.2 0.4]T。其中，UKF算法中狀態(tài)維數(shù)L=4；縮放比例參數(shù)λ=α2L-L，采樣點(diǎn)分布狀態(tài)由α來確定，經(jīng)過多次實(shí)驗(yàn)得出當(dāng)α=0.001時(shí)算法性能較好；而權(quán)系數(shù)β通常是非負(fù)數(shù)，其取值通常為β=2。IDE算法中Fmax=0.6，F(xiàn)min=0.2，CRmax=0.8，CRmin=0.2。仿真過程進(jìn)行500次迭代，仿真結(jié)果如圖3所示。

圖3描述了各參數(shù)的估計(jì)結(jié)果及其對應(yīng)的均方根誤差，其中擴(kuò)展卡爾曼粒子濾波(Extended Kalman Particle Filter，EPF)算法所代表的為文獻(xiàn)[23]中所提及的算法。從圖中可以看出文獻(xiàn)[23]中所采用的算法雖然解決了EKF算法對初值敏感的問題，但其線性化過程引起的誤差，導(dǎo)致估計(jì)結(jié)果波動幅度及均方根誤差相對較大；而UIDE-PF算法則有效減小了參數(shù)估計(jì)結(jié)果的波動幅度及均方根誤差。

(a) 直達(dá)信號幅度估計(jì)

圖4描述了采用上述估計(jì)結(jié)果實(shí)現(xiàn)多徑信號重構(gòu)，并去除多徑干擾后的相關(guān)函數(shù)值。從圖中可以看出采用UIDE-PF算法進(jìn)行多徑抑制后得到的相關(guān)函數(shù)曲線更平滑，且最接近真實(shí)值。同時(shí)，從圖的形狀上來看，采用UIDE-PF算法得到的相關(guān)函數(shù)曲線圖在碼延時(shí)為0時(shí)，對應(yīng)的幅值恰好為曲線的峰值點(diǎn)且附近沒有發(fā)生畸變。

(a) 采用PF算法實(shí)現(xiàn)多徑抑制

表1給出了4種不同算法每次迭代所消耗的時(shí)間。由于UPF算法與UIDE-PF算法在參數(shù)估計(jì)過程中需要對每個粒子進(jìn)行均值及協(xié)方差的計(jì)算，導(dǎo)致兩種算法時(shí)間消耗較長。但是對于以時(shí)速為30 m/s行駛的汽車而言，要使接收機(jī)定位精度達(dá)到2 m左右，即算法完成一次運(yùn)行時(shí)間不應(yīng)超過66.67 ms。而筆者所提出的算法運(yùn)行一次大約需15 ms，滿足接收機(jī)在高速環(huán)境下的定位需求。

表1 4種不同算法每次迭代所消耗的時(shí)間

4 結(jié)束語

筆者提出了一種無跡差分粒子濾波算法。該算法可有效地解決粒子濾波算法在參數(shù)估計(jì)過程中存在的粒子退化及多樣性降低的問題。相對于文中提到的粒子濾波算法、無跡粒子濾波算法及對比文獻(xiàn)中的算法，筆者所提算法在滿足了接收機(jī)對定位實(shí)時(shí)性要求的情況下，提高了參數(shù)估計(jì)精度，減小了估計(jì)結(jié)果均方根誤差。接下來，將針對如何減小算法的運(yùn)算量，進(jìn)一步提高參數(shù)估計(jì)精度展開研究。