基于格拉斯曼流形的狀態跟蹤方法研究*

趙 成，蘇嶺東，馬祥林

（1.杭州邦友安派智能科技有限公司，浙江 杭州 310012；2.國網徐州供電公司，江蘇 徐州 221000；3.常州致科自動化科技有限公司，江蘇 常州 213001）

0 引言

近年來，格拉斯曼流形在醫療科學、天文學和機器人控制等領域得到了廣泛運用，特別是在雷達信號處理過程中，m個基點觀察k個目標，格拉斯曼流形Gk，m-k可以作為一個新的目標群采樣空間[1]。而對于觀測過程中難以避免的噪聲干擾問題，國外學者Y.Chikuse 完成了很多基礎性工作并進行了有益探索。文獻[2]中提出一種將狀態空間模型置于格拉斯曼流行上，在后驗跟蹤模型的基礎上遞推估計的方法。但是，該方法將觀測模型和噪聲都限制在格拉斯曼流形上，使得估計值與實際值存在較大的偏差。本文中嘗試利用粒子濾波良好的濾波性能，提出了一種基于格拉斯曼流形的粒子濾波（Grass-Mann Manifolds-Paticle Filter，GM-PF）算法來估計流形上的隱馬爾科夫過程。

安排如下：首先，簡要介紹必要的背景知識，包括流形上的粒子濾波的概念；其次，提出詳細的動態模型和估計方法；再次，進行仿真并分析和討論；最后，簡要總結和討論今后的研究方向。

1 背景知識

1.1 格拉斯曼流形及其上的概率分布

格拉斯曼流形Gk，m-k可以看作是所有點為k維子空間v的空間集合，即在Rm上包含原點的k維超平面。每個在Gk，m-k上的k維子空間v都對應一個唯一的m×m正交投影矩陣P，冪為k。假如存在 在Vk，m上k列 的m×k矩 陣Y，則YY′=P；Pk，m-k表示全部冪為k的m×m正交映射矩陣的集合。因此，在流形Pk，m-k上進行的統計分析等于在格拉斯曼流形Gk，m-k上進行分析。Pk，m-k的隨機矩陣P分布的密度方程為：

pFq為超幾何函數矩陣。在文獻[3]中的分布上略加修改即可得到Pk，m-k的矩陣Langevin 分布，記為L(p)(m，k;B)。為了簡化分析，這里假設B的秩為k，則B的譜分解為：

其中：

其中在B上限制如下條件：trB=b被固定或者rankB=b＜m。

對于分布模式ΓΓ′，當B為半正定時，λi可能影響分布的形式。Chikuse and Watson 在文獻[4]中定義并討論了式（1）的分布在矩陣參數為多項式或超幾何函數的情況下的一般性。

為了更好地理解格拉斯曼流形上的概率分布，這里引入不變測度的概念。定義矩陣X∈Rmm，其點集X∈Rmm，關于矩陣變量的實值函數f(x)可視為矩陣元上面的多重積分。M為一黎曼流形，而Ξ是作用在M上的可測子集。那么，如果對于每一個可測量S?M，有μ(S)=μ(TS)=μ(ST)，T∈Ξ，μ是關于可測子集的不變測度。若令μ(M)=1，則μ是單位不變測度。

與通常定義在Rn上的Lebesgue 測度dX不同，Gk，m-k上的概率分布被定義表示為[dX]。格拉斯曼流形上的概率分布的數學表達形式為，其中為關于不變測度微分形式[dX]的矩陣變量的密度函數[5]。

1.2 基于格拉斯曼流形的粒子濾波

在經典的歐式空間信號處理過程中，一個典型的濾波框架由離散時間的隱馬爾科夫過程xt和帶噪聲觀測的yt=g(xt)+nt構成，其中g為映射，nt通常設為高斯白噪聲。概括來說，就是通過在t時刻關于未知狀態xt的一些條件函數和觀測變量y1，…，yt來估計E(f(xt)|y1，…，yt)。

如果選擇噪聲方差為高斯密度函數，就是人們熟知的卡爾曼濾波器（Kalman Filter），映射g為線性映射，p(xt|xt-1)為馬爾科夫過程[6]。在這種情況下，大量統計得到的后驗分布可以接近真實分布。但是，出現非線性非高斯噪聲時，常規的遞歸濾波方法難以處理，使得需尋求一種近似計算的方法。一種常用的近似方法將局部線性化，即擴展卡爾曼濾波（Extended Kalman Filter，EKF）。另一種基于序貫重要性采樣的近似方法，即粒子濾波（Particle Filter）。在重要性分布上抽取N個樣本（即“粒子”），根據密度p(xt|xt-1)和p(yt|yt)遞歸更新歸一化重要性權值，則后驗期望E(f(xt)|y1，…，yt)的估計可由得到。顯然，重要性函數的選擇直接決定變量的估計和濾波器的性能[7-9]。

然而，在格拉斯曼流形框架下粒子濾波方法同樣通過計算加權均值樣本來接近期望的函數E(f)。設點位于格拉斯曼流形Gk，m-k上，目前基本上有兩種方法計算他們的均值。外蘊均值法將流形中嵌入一個歐式空間，在歐式空間上計算算術平均值，然后將結果點映射回流形上。但是，假如最接近的點不唯一，外蘊均值法會產生多個均值或者均值為空。內蘊均值法充分利用流形的幾何結構[10]，可以有效解決由于指數映射導致的低效問題，但同樣存在缺陷，需要解決相同尺度下最近點的最小化問題（如Karcher 均值法和質心法[11]）。本文采用外蘊均值法，內蘊均值計算法是今后工作的重點。

2 模型選擇

2.1 矩陣的概率分布

一個m×k隨機矩陣X的正態分布如下：

M為均值矩陣，Ω為協方差矩陣，上的概率密度。dX為Lebesgue 測度，tr(·)表示跡。

在格拉斯曼流形上一個m×k隨機矩陣X和參數矩陣B的von Mises-Fisher 分布（也被稱為Langebin 分布）為[12]：

其中C為標準化常數，[dX]為在格拉斯曼流形上的單位不變測度，p(X)[dX]=1[dX]是相應的均勻分布。

矩陣的von Mises-Fisher 分布與歐幾里得空間的正態分布密切相關。假如一個始于矩陣正態密度的隨機矩陣X，參數，限制條件X′X=I，易得到由此產生的密度Gmpd(X;Ω-1G)。所以，將選擇不同的噪聲能量對比GM-PF 算法和GM-PMTA 算法比較優劣。

2.2 GM-PF 動態模型

為了解決在統計信號處理過程中的矩陣變量的估計和跟蹤問題，這里介紹一種基于格拉斯曼流形的動態空間模型。假定離散時間的隱馬爾科夫過程Xt位于格拉斯曼流形Gk，m-k上，則矩陣的von Mises-Fisher 分布為

設觀測Yt由Yt=Xt+Nt得到，Nt為均值為零的高斯矩陣，參數，所以觀測模型為：

2.3 后驗跟蹤算法的狀態空間模型

與上文不同，文獻[2]中的狀態空間模型的觀測被限制于格拉斯曼流形上，包含如下兩個von Mises-Fisher 分布。Xt為在格拉斯曼流形上的隱馬爾科夫過程：

觀測模型為：

注意到觀測誤差很有可能是量測噪聲干擾的結果，所以這里優先考慮2.2 節提出的動態模型。但是，文獻[2]中提出的后驗模型跟蹤算法可以作為比較的基礎，給定觀測噪聲來仿真兩種模型，設β=Ω-1。

3 仿真實驗

3.1 算 法

結合2.2 節和2.3 節提出的不同狀態空間模型得到詳細的濾波算法。這種算法需要在mk? 空間上的點Y映射到Gk，m-k上，根據文獻[2]中的定義，通過對Y進行極分解得到X，即。

根據2.2 節的模型，假設有N個粒子，初始化權值，i-1，…，N。當t≥1 時粒子濾波過程如下所示[5，13-14]。

算法1：格拉斯曼流形上的粒子濾波

1.得到觀測Yt～N(Yt;Xt，Ω，∑)

2.fori=1，…，N

4.是否需要重采樣

根據2.3 節中的模型，用χt來表示Xt的估計。文獻[2]將χ0=X0，但在這里將χ0直接隨機均勻抽取。第一個觀測量Y1～Gmpd(Y1;X1β)，則。當t≥2 時，算法過程如下。

算法2：后驗模型跟蹤算法

3.fori=1，…，收斂

3.2 仿真結果分析

將算法1 用于2.2 節中的模型，其中算法仿真時間步為125，過程噪聲方差0.01，觀測噪聲方差0.25，每個參數都固定實驗30 次取平均值。圖1 為在不同粒子數的情況下算法的最小均方誤差。如圖1 所示，隨著粒子數的加大，基于格拉斯曼流形的粒子濾波方法的誤差逐漸減小，但當粒子數超過100 時，GM-PF 方法付出了巨大的計算代價，而性能提高并不明顯。因此，為了兼顧濾波精度和效率，下文仿真中都將設置N=100。

圖1 在G1，3上，粒子數N=25 到N=1 000 時，GM-PF 算法的MMSE

在格拉斯曼流形上的馬爾科夫過程過程Xt，觀測Yt通過式（7）產生。最小均方差（Minimum Mean Square Error，MMSE）狀態估計根據算法1 得到。圖2 和圖3 顯示了一個在格拉斯曼流形G1，3上典型的HMM 過程的跟蹤效果，，兩種算法的跟蹤性能比較。圖2、圖3 比較可以發現，兩種算法都能成功跟蹤隱馬爾科夫過程Xt的真值，但是后驗模型跟蹤算法的狀態估計更為平滑，粒子濾波算法的MMSE 估計更接近真實Xt曲線。此外發現，算法1 計算所需時長是算法2 的10 倍，但是算法1可得到更低的均方誤差。

為了將算法2 應用于最大后驗（MAP）估計，將觀測Yt映射ε(Yt)到格拉斯曼流形上產生“偽觀測”。與映射前對比發現，這種簡單的映射χt=ε(Yt)也能夠產生一定程度的“降噪”效果。

圖2 在G1，3上相同隱馬爾科夫過程和噪聲條件下，GM-PF 濾波算法實驗結果

圖3 在G1，3上相同隱馬爾科夫過程和噪聲條件下，GM-PMTA 方法實驗結果

圖4 根據2.2 節的模型在G1，3上對比幾種不同方法的實驗結果。首先固定過程參數α，對比各種方法在不同噪聲情況下的估計誤差，其中Yt和Xt之間的距離用Frobenius 范數處理，過程參數α=10I。仿真時間設為100 時間步，粒子數N=100，ε(Yt)在格拉斯曼流形上為GM-PMTA 算法產生合適的觀測量。圖4 中的實驗結果都是同一參數30 次實驗取均值的結果，Ω=σ2I，β=σ-2I。其中，GM-PF為格拉斯曼流形粒子濾波后得到的觀測誤差，GMPGMA 就是后驗模型跟蹤算法（PMTA）濾波后得到的觀測誤差，Projected Observed 為映射后觀測到的誤差，GM Observed 為在流形G1，3上未進行濾波觀測到的誤差。

在圖4 中可以清楚看出，GM-PF 方法勝過GM-PMTA 方法，特別是在噪聲增大的時候（相同的情況同樣出現在2.3 節中的模型）。當噪聲增大到一定程度后，GM-PMTA 方法的性能迅速下降，甚至會比僅僅將狀態估計映射ε(Yt)所產生的降噪效果還要差。

圖4 不同噪聲能量σ2 情況下，在G1，3上不同方法的MMSE變化曲線

4 結語

本文介紹了一種矩陣變量的動態模型和與之匹配的粒子濾波算法。在低維空間上的仿真表明，相比較確定性模型跟蹤算法，該方法顯著提高了濾波性能，但付出了較高的計算代價。為了兼顧粒子濾波方法的精度和效率，未來研究的重點放在利用流形內部的幾何結構進行內蘊均值計算，或者通過傳遞矩陣并行傳遞來替代指數映射，避免反復指數映射帶來的大量計算，在保證較高精度的情況下提高效率。另外，著重研究一種在矩陣von Mises-Fisher 分布上高效的采樣方法，使得這種基于格拉斯曼流形的粒子濾波方法能夠解決高維空間的狀態估計問題。