《目標(biāo)跟蹤系統(tǒng)中的濾波方法》課件第8章

第8章RTS平滑及分段融合方法

8.1引言8.2RTS平滑算法8.3基于分段RTS平滑的凸組合航跡融合算法8.4小結(jié)

8.1引言

多傳感器信息融合系統(tǒng)的結(jié)構(gòu)可以分為集中式和分布式。集中式融合需要將所有局部傳感器量測數(shù)據(jù)都傳送到中心節(jié)點，網(wǎng)絡(luò)傳輸負(fù)載大，對中心節(jié)點的處理能力要求高。在分布式融合系統(tǒng)中，每個傳感器均有自己的處理器，在各自節(jié)點進行預(yù)處理，然后將結(jié)果送到中心節(jié)點，進行融合處理。由于融合中心的主要任務(wù)是對各局部航跡進行融合，所以這種融合結(jié)構(gòu)也稱為航跡融合[1]。航跡融合網(wǎng)絡(luò)負(fù)載小，對中心節(jié)點的處理能力要求不高，所以航跡融合算法歷來是人們研究的熱點。最易于實現(xiàn)的航跡融合算法是凸組合融合算法[2]，該算法不考慮局部航跡之間的相關(guān)性，直接進行線性加權(quán)平均，涉及到的參數(shù)少，算法簡單，運算量小。其它的航跡融合算法還有Bar-Shalom-Campo融合算法、最優(yōu)分布式估計融合、協(xié)方差交叉法以及聯(lián)邦濾波器融合算法等[1-3]。也有學(xué)者采用模糊邏輯的方法進行航跡融合[4-5]，并且取得了較好的整體融合性能。但是，這些算法涉及到的參數(shù)多，運算量大。文獻[6]在不考慮融合實時性的前提下，應(yīng)用RTS(Rauch-Tung-Striebel)平滑算法，取得了較好的濾波效果。Willsky等人提出了多尺度自回歸估計融合算法[7]，該算法綜合了卡爾曼濾波和RTS平滑算法，提高了多尺度估計融合性能。本章考慮到RTS平滑算法[8-10]對噪聲的抑制效果要比單純的卡爾曼濾波結(jié)果好，所以將RTS平滑算法應(yīng)用到融合過程中，以提高航跡融合性能。

RTS平滑算法的逆向平滑過程從最后一個濾波結(jié)果反向遞推，逐步向前平滑，這就意味著獲得了整個目標(biāo)航跡跟蹤過程之后，才使用RTS平滑算法對系統(tǒng)航跡進行平滑，平滑步驟嚴(yán)重滯后，融合結(jié)果的實時性難以得到保證。為了避免輸出結(jié)果大幅度延遲，本章提出了分段RTS平滑算法。該算法思想是對需要平滑的時間段內(nèi)的濾波結(jié)果分段，然后逐段逆向平滑，由于分段長度可大可小，從而在不同程度上改善融合過程中應(yīng)用RTS平滑造成的實時性欠佳問題。仿真結(jié)果表明，分段平滑算法的分段長度對整個航跡融合結(jié)果的誤差影響不大，并且比傳統(tǒng)航跡融合算法效果好。將分段RTS平滑過程應(yīng)用到凸組合航跡融合過程中，能夠較好地克服由于過程噪聲使得凸組合航跡融合算法性能降低的問題。在融合過程中，針對局部傳感器有無額外計算能力，結(jié)合實施分段平滑過程的時機，本章提出了基于分段RTS平滑的先平滑再融合和先融合再平滑兩種改進的凸組合航跡融合算法。仿真結(jié)果表明，這兩種算法的航跡融合性能均高于凸組合航跡融合算法和最優(yōu)航跡融合算法。

對于處理非線性系統(tǒng)的擴展卡爾曼濾波、不敏卡爾曼濾波[11-12]、高斯厄米特卡爾曼濾波[13-14]、容積卡爾曼濾波[15]和粒子濾波[16-17]，對應(yīng)的RTS平滑算法有擴展卡爾曼RTS平滑算法、不敏卡爾曼RTS平滑算法[18]、高斯RTS平滑的通用形式[19-20](包括高斯厄米特RTS平滑算法、容積卡爾曼RTS平滑算法)和粒子RTS平滑算法[21-22]等。本章對這些處理非線性系統(tǒng)的平滑算法也做了論述。這些平滑算法都可以結(jié)合航跡融合算法進行分段融合，并能夠獲得較好的融合效果。

8.2RTS平滑算法

8.2.1卡爾曼濾波RTS平滑

最優(yōu)平滑的目的是在接收到T時刻的量測值后，計算k時刻狀態(tài)向量xk的后驗邊緣分布，即p(xk|y1:T)，此處T>k。濾波和平滑之間的差別在于，最優(yōu)濾波器僅僅使用直到時刻k之前的量測值進行估計，而最優(yōu)平滑器不僅使用k時刻之前的量測，而且使用k時刻之后的量測值。在獲得濾波后驗狀態(tài)分布之后，用貝葉斯最優(yōu)平滑方程計算從k時刻到T時刻中各時刻的邊緣后驗分布。根據(jù)馬爾可夫特性，在狀態(tài)xk+1給定的條件下，狀態(tài)向量xk與yk+1:T相互獨立，由此可得p(xk|xk+1，y1:T)=p(xk|xk+1，y1:k)。使用貝葉斯公式，在xk+1和y1:T給定的條件下，狀態(tài)向量xk的概率分布可以表示為(8-1)在y1:T給定的條件下，xk和xk+1的聯(lián)合分布可以計算為(8-2)其中，在k+1時刻平滑后，狀態(tài)向量的概率為p(xk+1|y1:T)。通過對xk+1求積分，獲得在y1:T條件下xk的邊緣分布，即(8-3)因此，計算平滑分布p(xk|y1:T)的后向遞推方程通過下列貝葉斯平滑方程給出(8-4)(8-5)其中，p(xk|y1:k)表示在時刻k濾波后的概率分布，p(xk+1|y1:k)是k+1時刻的預(yù)測概率分布，若將狀態(tài)向量看做為離散分布，則在上式的積分運算可以被替換為累加和。

RTS平滑器可以用來求解平滑過程，該方法獲得近似解為(8-6)為了獲得RTS平滑器的方程，首先給出多元高斯分布的定義及其中的兩個引理。

若隨機變量x∈Rn是具有均值為m∈Rn、協(xié)方差為P∈Rn×n的高斯分布，則其概率密度函數(shù)為(8-7)其中|P|為矩陣P的行列式。

引理8.1

多元高斯變量的聯(lián)合密度。如果隨機變量x∈Rn和y∈Rm具有如下概率密度(8-8)則x，y的聯(lián)合概率密度和y的邊緣分布概率密度為(8-9)(8-10)

引理8.2

多元高斯變量的條件密度。如果隨機變量x,y具有如下聯(lián)合高斯概率密度(8-11)則x，y的邊緣和條件概率密度如下(8-12)(8-13)(8-14)定理8.1RTS平滑器中，離散固定區(qū)間RTS平滑器的后向遞推方程如下證明：根據(jù)多元高斯分布中的引理8.1，給定y1:k時，xk和xk+1的聯(lián)合高斯分布為(8-20)其中(8-21)根據(jù)狀態(tài)向量的馬爾可夫特性，可得(8-22)于是借助于引理8.2，可得條件概率分布(8-23)其中給定所有量測數(shù)據(jù)條件下xk和xk+1的聯(lián)合概率分布為

(8-27)

其中(8-28)由此通過引理8.2，xk的邊緣分布為(8-29)(8-30)(8-31)其中證畢。8.2.2高斯RTS平滑算法的通用形式

首先假設(shè)濾波結(jié)果的高斯近似分布為(8-32)此外，假定在時刻k+1進行平滑后的近似高斯分布為(8-33)根據(jù)狀態(tài)向量的馬爾可夫特性，在給定xk+1的條件下，xk和yk+1|T相互獨立，因此有(8-34)借助于貝葉斯公式，在xk+1和y1:T條件下，xk的條件概率密度函數(shù)可計算為(8-35)也就是說，p(xk|xk+1，y1:T)的高斯近似形式由以下兩項構(gòu)成：

(1)構(gòu)造聯(lián)合分布p(xk，xk+1|y1:k)的高斯近似形式；

(2)以xk+1為條件，計算上述聯(lián)合分布，可得p(xk|xk+1,y1:T)的高斯近似形式。

這里需要采用假設(shè)密度濾波的方法，通過近似計算分布p(xk，xk+1|y1:k)的前二階矩。注意到因為該分布獨立于未來量測值yk+1:T，因此近似也僅僅依賴于濾波結(jié)果和系統(tǒng)動態(tài)模型。

如果我們假定濾波結(jié)果是高斯分布p(xk|y1:k)=N(xk|mk|k，Pk|k)，則近似推導(dǎo)過程如下。首先計算高斯積分為(8-36)(8-37)(8-38)下列聯(lián)合分布的高斯近似可以構(gòu)造為(8-39)式(8-35)的除法結(jié)果仍然是高斯分布，因此根據(jù)高斯分布的計算規(guī)則可得(8-40)(8-41)(8-42)其中如果定義平滑增益為(8-43)則上述方程可寫為(8-44)(8-45)給定所有觀測數(shù)據(jù)條件下，xk和xk+1的聯(lián)合分布現(xiàn)在變?yōu)?8-46)其中(8-47)(8-48)xk的邊緣分布為(8-49)其中(8-50)(8-51)因此，該通用的固定區(qū)間平滑器以濾波結(jié)果mT|T，PT|T為起點，按照k=T－1，T－2，…，0的次序做如下步驟的計算：

(1)預(yù)測。根據(jù)式(8-36)~式(8-38)計算預(yù)測均值mk+1|k，預(yù)測協(xié)方差Pk+1|k和互協(xié)方差Ck，k+1。

(2)平滑。計算平滑增益Gk、平滑均值mk|T和協(xié)方差Pk|T如下(8-52)(8-53)(8-54)上述方法可以更一般化到非加性噪聲的情形，此時，狀態(tài)轉(zhuǎn)移方程和量測方程如下(8-55)(8-56)其中，qk－1和rk是服從高斯分布的高斯噪聲。用更為一般的形式替換式(8-36)~(8-38)中的積分(8-57)(8-58)(8-59)需要指出的是，在不敏變換情況下該方法等價于在UKF中使用的擴維方法。實際上，隨機變量qk也可以被替換為非高斯分布，只需要在式(8-57)~式(8-59)中將高斯分布N(qk|0，Qk)用非高斯噪聲分布函數(shù)進行替換即可。8.2.3不敏卡爾曼濾波RTS平滑

結(jié)合不敏變換和高斯RTS平滑的通用形式，具有加性噪聲形式的不敏RTS平滑算法步驟如下。

Step1：構(gòu)造西格瑪點集:(8-60)(8-61)(8-62)其中，λ=α2(n+κ)－n，n為狀態(tài)向量維數(shù)，運算符[]i表示矩陣的第i列。

Step2：借助狀態(tài)轉(zhuǎn)移方程計算傳播后的西格瑪點集:(8-63)

Step3：計算預(yù)測均值mk+1|k、預(yù)測協(xié)方差矩陣Pk+1|k和互協(xié)方差矩陣Dk+1:(8-64)(8-65)(8-66)其中權(quán)值W(m)i表示第i個西格瑪點對應(yīng)的權(quán)值，W(c)i表示第i個西格瑪點對應(yīng)的協(xié)方差矩陣的權(quán)值，為方便起見,常令W(m)i=W(c)i，其定義見UKF濾波一節(jié)。

Step4：計算平滑增益Gk、平滑均值msk|k和協(xié)方差Psk|k:(8-67)(8-68)(8-69)上述步驟的計算是從最后一個時刻的濾波結(jié)果msT|T=mT|T，PsT|T=PT|T開始的，遞推時刻序列為k=T－1，…，0。

假定濾波結(jié)果的近似均值和協(xié)方差為

p(xk|y1:k)≈N(xk|mk|k，Pk|k)

(8-70)

時刻k+1的平滑分布已知且近似高斯分布:

p(xk+1|y1:T)≈N(xk+1|msk+1|k+1，Psk+1|k+1)

(8-71)

對最優(yōu)平滑估計的基于不敏變換的近似方法推導(dǎo)如下。

Step1：給出xk和xk+1的聯(lián)合分布的高斯近似如下(8-72)上式中變量可以通過式(8-64)~式(8-66)求解，也就是利用不敏變換求解具有加性噪聲的狀態(tài)估計過程。

Step2：由于式(8-72)是高斯分布的，借助于高斯分布的計算規(guī)則，xk的條件分布給定如下其中

Step3：給定所有觀測數(shù)據(jù)，xk，xk+1的聯(lián)合分布為(8-77)其中(8-78)(8-79)于是xk的邊緣分布為(8-80)其中(8-81)(8-82)采用類似的方法，結(jié)合擴維不敏卡爾曼濾波的思路，易于獲得擴維的不敏RTS平滑算法。8.2.4高斯厄米特RTS平滑算法

在加性高斯RTS平滑算法的通用形式中使用高斯厄米特近似方法，可以獲得高斯厄米特RTS平滑算法，算法過程如下：

Step1：用一維單位西格瑪點的笛卡爾乘積構(gòu)造多維單位西格瑪點，即(8-83)Step2：構(gòu)造西格瑪點如下(8-84)

Step3：將西格瑪點帶入狀態(tài)轉(zhuǎn)移方程，計算西格瑪點的狀態(tài)轉(zhuǎn)移值(8-85)

Step4：用一維權(quán)值的乘積構(gòu)造多維權(quán)值(8-86)其中每一個ik取值為1，…，p。

Step5：計算預(yù)測均值mk+1|k、預(yù)測協(xié)方差Pk+1|k以及互協(xié)方差Dk+1，即(8-87)(8-88)(8-89)

Step6：增益Gk，狀態(tài)向量估計msk|k和協(xié)方差矩陣Psk|k計算如下(8-90)(8-91)(8-92)8.2.5容積卡爾曼濾波RTS平滑

在加性高斯RTS平滑算法的通用形式中使用3階球形容積近似方法，可以獲得加性噪聲條件下的容積卡爾曼濾波RTS平滑器，該算法過程如下:

Step1：構(gòu)造西格瑪點(8-93)其中單位西格瑪點定義為(8-94)Step2：將西格瑪點代入狀態(tài)轉(zhuǎn)移方程，可得(8-95)

Step3：計算預(yù)測均值mk+1|k、預(yù)測協(xié)方差矩陣Pk+1|k和互協(xié)方差矩陣Dk+1：(8-96)(8-97)(8-98)

Step4：計算增益Gk、均值msk|k和協(xié)方差矩陣Psk|k：(8-99)(8-100)(8-101)在非加性形式的高斯RTS平滑算法的通用形式中使用3階球形容積近似方法，可以獲得非加性噪聲條件下的容積卡爾曼濾波RTS平滑器，該算法過程如下：

Step1：由n′=n+nq維擴維隨機變量(xTk，qTk)T構(gòu)造西格瑪點集(8-102)其中，。Step2：通過狀態(tài)轉(zhuǎn)移方程傳播西格瑪點集(8-103)

Step3：計算預(yù)測均值mk+1|k、預(yù)測協(xié)方差矩陣Pk+1|k和互協(xié)方差矩陣Dk+1(8-104)(8-105)(8-106)

Step4：計算增益Gk、均值msk|k和協(xié)方差矩陣Psk|k，(8-107)(8-108)(8-109)8.3基于分段RTS平滑的凸組合航跡融合算法

8.3.1分段RTS平滑算法

假設(shè)某一動態(tài)系統(tǒng)的總時長為k=1，2，…，N。分段RTS平滑算法就是對N個觀測值分段進行前向濾波和逆向平滑，假設(shè)分段長度為L，且1≤L≤N，則分段RTS平滑算法的步驟如下：

Step1：輸出變量清空，即xs=[]，Ps=[]。

Step2：局部濾波結(jié)果變量清空，即M=[]，P=[]，分段平滑計數(shù)器清零，即count=0。

Step3：如果對目標(biāo)的跟蹤沒有結(jié)束，則繼續(xù)向下處理；否則輸出結(jié)果xs，Ps，整個分段RTS平滑算法結(jié)束；

Step4：接收當(dāng)前濾波結(jié)果并保存在變量M，P中，并采用計數(shù)器count計數(shù)。若count=L，則針對局部濾波結(jié)果變量M，P采用常規(guī)RTS平滑算法進行逆平滑過程，平滑結(jié)果保存在輸出變量xs，Ps中；轉(zhuǎn)到Step2。

算法中，x，P既可以是局部節(jié)點的濾波狀態(tài)及方差陣，也可以是中央節(jié)點對各個局部節(jié)點濾波結(jié)果融合之后的狀態(tài)值及方差陣；xs和Ps表示分段RTS平滑之后的狀態(tài)值及方差陣。

如果定義分段RTS平滑算法的響應(yīng)時間為從數(shù)據(jù)采樣開始到獲得第一個平滑結(jié)果為止的時間，則有

tr=(tsampling+tfiltering+tsmoothing)L

(8-110)^其中，tsampling表示平均采樣時間，tfiltering表示平均濾波時間，tsmoothing表示平均逆向平滑時間，L表示RTS平滑算法分段的長度。當(dāng)L≥N時，分段RTS平滑算法退化為常規(guī)RTS平滑算法；當(dāng)L<N時，分段RTS平滑算法的響應(yīng)時間比不分段的RTS平滑算法的響應(yīng)時間短。

融合中心接收到L個濾波數(shù)據(jù)后，立刻進行逆向平滑過程，不必等待后續(xù)數(shù)據(jù)的到來，這樣做可以降低逆向平滑的滯后時間。分段長度L取值越小，融合算法的實時性就越好。8.3.2算法描述及分析

在凸組合航跡融合過程中應(yīng)用分段RTS平滑算法，針對局部節(jié)點有無額外計算能力，結(jié)合實施逆向平滑過程的時機，將其分為兩種：

(1)先平滑再融合(SmoothingFirstandFusingNext，SFFN)，即先對局部航跡進行分段RTS逆向平滑，然后再利用凸組合算法進行融合。

(2)先融合再平滑(FusingFirstandSmoothingNext，F(xiàn)FSN)，即先采用凸組合算法融合局部航跡，再對融合后的航跡進行分段RTS逆向平滑。

這兩種方法的融合過程分別如圖8.1(a)和圖8.1(b)所示。圖8.1基于RTS平滑的航跡融合過程在工程實踐過程中，如果局部節(jié)點有額外的計算能力，就將逆向平滑過程放在局部節(jié)點處理，即采用SFFN；如果局部節(jié)點沒有額外的計算能力，就將逆向平滑過程放在中心節(jié)點處理，即采用FFSN。由于在中心節(jié)點實施逆向平滑過程使用的數(shù)據(jù)是融合后的，經(jīng)過的處理比在局部節(jié)點多，因此從理論上講，采用SFFN算法比采用FFSN算法效果好。仿真結(jié)果與該結(jié)論一致。

假定各個局部節(jié)點的數(shù)據(jù)采樣是同步的，并且忽略數(shù)據(jù)在分布式傳感器網(wǎng)絡(luò)中的傳輸時間，那么航跡融合系統(tǒng)的響應(yīng)時間可以采用下式估算(8-111)其中，tr1、tr2、tr3分別表示SFFN、FFSN和CCF(ConvexCombinationtrack-to-trachFusion)三種算法的響應(yīng)時間，

tsampling表示平均采樣時間，tfiltering表示平均濾波時間，tfusion表示平均融合時間，tsmoothing表示平均逆向平滑時間，L表示RTS平滑算法分段的長度。當(dāng)L≥N時，分段RTS平滑算法退化為常規(guī)RTS平滑算法，故此時前兩種算法的響應(yīng)時間并未減少；當(dāng)1≤L<N，尤其是當(dāng)L<<N時，分段RTS算法的響應(yīng)時間能夠大大縮短。由以上分析可知，當(dāng)分段長度選取得比較小時，分段平滑算法可以有效地縮短融合結(jié)果響應(yīng)時間，從而確保航跡融合的實時性。在使用SFFN和FFSN算法過程中，如果分段RTS平滑算法分段長度過長，則融合結(jié)果滯后，無法滿足融合實時性要求；如果分段長度過短，則融合性能有可能達不到最優(yōu)。這就需要在實時性和精確性之間做折中。由于最優(yōu)分段長度難于提前估算，因此使用過程中只能憑借經(jīng)驗選取合適的分段長度。本章實驗一和實驗三選取的分段長度值為1，試驗結(jié)果表明，雖然所選取的分段長度并不能保證估計結(jié)果的精度達到理想值，但是同樣能夠大幅度提高航跡融合性能。8.3.3仿真實驗及結(jié)果分析

假設(shè)某系統(tǒng)狀態(tài)變量為x=[ξ

ξ]T，其離散狀態(tài)方程和量測方程為·(8-112)其中，i=1，2表示傳感器1和傳感器2。過程噪聲Q和量測噪聲R

i是相互獨立的零均值高斯白噪聲，過程噪聲、量測噪聲和狀態(tài)初始值互不相關(guān)。實驗中取x0=[0

10]T，P0=[101;110]，R1=5，R2=15，T=1，仿真步數(shù)N=100步。實驗一算法性能分析

采用分段長度L=1和過程噪聲Q=10的基于分段RTS平滑的凸組合融合算法進行1000次MonteCarlo仿真，對比SFFN、FFSM、CCF和OF(Optimaltrack-to-trachFusion)四種融合算法的均方根誤差。

實驗結(jié)果如圖8.2(a)和圖8.2(b)所示。結(jié)果表明，在融合過程中應(yīng)用分段長度為1的RTS平滑過程后，融合效果明顯得到改善，狀態(tài)向量兩個分量的誤差大幅度降低。從實驗結(jié)果來看，所提出的兩種方法比CCF和OF算法的誤差性能好。圖8.2融合算法誤差性能對比實驗二分段長度對算法性能影響分析

令過程噪聲Q=10，依次改變RTS平滑算法的分段長度L=1，2，3，…，N，分別采用不同的分段長度進行MonteCarlo仿真100次，計算不同分段長度仿真結(jié)果的均方根誤差。當(dāng)使用平滑尺寸參數(shù)k進行仿真時，其均方根誤差的計算式為(8-113)實驗結(jié)果如圖8.3(a)和圖8.3(b)所示，橫坐標(biāo)表示分段的長度。實驗結(jié)果表明，分段RTS平滑算法的分段長度對融合結(jié)果影響不大，隨著分段長度的增加，融合誤差趨向常數(shù)。而且，不論對分段RTS平滑算法的分段長度如何選擇，基于分段RTS航跡融合方法的效果總是優(yōu)于C