基于模糊控制交互式多模型粒子濾波的靜電機動目標跟蹤

付巍，鄭賓

(中北大學 計算機與控制工程學院，山西 太原030051)

0 引言

靜電探測技術是一種通過測量目標自身的靜電場獲取目標運動狀態信息的探測技術［1］。基于靜電探測技術的探測器具有被動式、反低空、反隱身的優勢，可被應用于探測低空飛行的目標。目前已經設計出能夠應用于探測空中目標的靜電探測器主要有短路軸向式探測器［2］、電極掃描式探測器［3］、球形探測器［4］和旋轉式三維靜電探測器［5］等。這些探測器所得到的探測方程包含探測器與目標之間的相對距離和速度等真實信息，并且這些探測方程都是非線性的，探測數據中包含著非高斯分布的噪聲，因此，需要利用目標跟蹤算法根據靜電探測器當前采樣時刻的測量信息和前一個采樣時刻的目標運動狀態的預測值，估算出目標當前采樣時刻的運動狀態，以實現對機動目標的有效跟蹤。

在機動目標跟蹤中，目標經常有多種運動狀態，單一模型的目標跟蹤算法不能滿足對機動目標跟蹤性能的要求。交互式多模型算法假定目標的運動狀態可用有限個運動模型描述，并為每個運動模型分配一個濾波算法和概率，根據當前采樣時刻的測量信息在不同的運動模型之間進行交互和更新，實現對機動目標的跟蹤［6］。交互式多模型算法對機動目標跟蹤的效果主要依靠各運動模型概率和運動模型間的轉換概率矩陣［7］。如果概率較大的運動模型與當前采樣時刻目標的實際運動狀態相一致，跟蹤的誤差就小，反之就大。交互式多模型算法的模型轉換概率矩陣是先驗矩陣，即在算法開始前就預先指定的，如果能夠根據目標的實際運動狀態情況實時調整模型轉換概率矩陣每一列的值，使與目標當前運動狀態最接近的運動模型在混合產生這一采樣時刻的初始狀態向量里占有更大的比重，交互式多模型算法將會更有效地實現對機動目標的跟蹤。為此，提出一種基于模糊控制的交互式多模型算法，即利用模糊控制方法來實現對模型轉換概率矩陣的實時調整，使與目標當前運動狀態最接近的運動模型在混合產生這一采樣時刻的初始狀態向量里占有更大的比重。

為了在交互式多模型算法中實現對目標運動狀態的更新，并解決探測方程非線性和探測數據含有非高斯噪聲的問題，常采用粒子濾波算法。粒子濾波的效果在很大程度上取決于重要性密度函數的選取，重要性密度函數選擇不好，常會導致濾波發散或較差的濾波效果［8］。為了提高粒子濾波算法的精度，減小算法更新時間，提出利用中心差分擴展卡爾曼濾波算法得出基本粒子濾波的建議分布函數，通過這個建議分布函數來代替基本粒子濾波算法的重要性密度函數，實現對目標運動狀態的更新。

1 靜電探測系統數學模型

靜電探測器采用旋轉式三維靜電探測器，該探測器有x、y、z 3 個方向的感應電極。假設目標靜電場的三維電場強度分量Ex、Ey、Ez，在感應電極上產生感應電量Qx、Qy、Qz，Qx、Qy、Qz隨時間的變化率就是電流，因此，測量3 個方向上的感應電極的電流變化，就能夠確定目標靜電場的三維電場強度分量的變化，通過目標靜電場的三維電場強度分量可以計算得到目標與靜電探測器之間的方位角和俯仰角。

利用旋轉式三維靜電探測器探測飛行目標的方式如圖1所示。當目標與探測器的距離較遠，即距離大于目標尺寸的5 ～6 倍時，可以認為目標為點目標［9］。圖1中探測器的架設高度為h，被探測目標的初始坐標為(x0，y0，z0)，并以沿x 軸、y 軸和z 軸的速度vx、vy、vz從P 點向探測器飛來。

圖1 旋轉式三維靜電探測器探測低空目標示意圖Fig.1 Schematic diagram of the rotating three-dimensional electrostatic detector detecting a low-altitude target

在旋轉式三維靜電探測器3 個方向的感應電極上產生的電場強度的三維分量分別為

式中:ax，ay和az分別為沿x，y 和z 軸的單位方向向量;R1為目標距探測器電極的距離;K 為電場強度系數。

則旋轉式三維靜電探測器與目標之間的方位角

因為目標飛行高度z0比探測器架設高度h 大很多，所以R1≈R. 旋轉式三維靜電探測器與目標之間的俯仰角

旋轉式三維靜電探測器在每一個采樣時刻可以同時得到當前時刻目標與靜電探測器之間的方位角和俯仰角等角度信息。利用探測器測得的角度信息，結合目標跟蹤算法就能夠估算出目標當前的運動參數，如位置、速度等信息，實現對機動目標的跟蹤。

2 基于模糊控制的交互式多模型粒子濾波(FIMMPF)算法

2.1 交互式多模型粒子濾波(IMMPF)算法及改進

假設動態系統狀態空間模型為

交互式多模型算法的核心思想是假定目標的運動狀態可用有限個運動模型描述，各運動模型之間的轉換概率由如下一階馬爾可夫鏈約束:

式中:Ψij表示運動模型轉換概率矩陣;rk表示機動運動目標在第k 個采樣時刻所屬的運動模型，其值域為{1，…，M}，M 表示系統包含運動模型的數目。

在第k 個采樣時刻系統處于第rk=i 個運動模型的概率記為

利用中心差分擴展卡爾曼濾波算法對基本粒子濾波的建議分布進行構造，形成改進的IMMPF 算法，其運算流程如下:

1)計算運動模型混合概率。計算運動模型混合概率的公式為

2)計算基于運動模型的混合狀態初值和混合初始狀態估計誤差協方差矩陣。設粒子點的個數為N，而運動模型個數為M，這些粒子點在第k 個采樣時刻第i 個運動模型產生的狀態估計向量和狀態估計誤差協方差矩陣可以表示為Xni(k)，Pni(k)，其中:i=1，…，M;n=1，…，N.

混合所有運動模型的相應的粒子點，則第j 個運動模型的混合狀態初值及初始狀態估計協方差矩陣為

3)用中心差分擴展卡爾曼濾波算法進行濾波。利用第k+1 采樣時刻的測量數據Z(k+1)，對每個運動模型中的N 個粒子點，通過中心差分擴展卡爾曼濾波算法進行濾波，得到第k+1 個采樣時刻的每一個粒子點的運動狀態向量的更新值，包括更新的狀態向量(k+1)，更新的狀態向量估計的誤差協方差矩陣(k+1)，似然函數Lnj(k +1)和每個粒子點的權重(k+1).

采用中心差分擴展卡爾曼濾波算法進行建議分布的構造，其步驟如下:

預測狀態向量

預測協方差矩陣

式中:Q(k)為狀態方程噪聲的協方差矩陣。

殘差協方差矩陣:

令

則

式中:H(k)為觀測方程狀態轉移矩陣;R(k)為觀測方程噪聲的協方差矩陣。

令

應用中心差分理論

式中u 表示中心差分的步長，則

濾波增益矩陣

濾波輸出

誤差協方差矩陣

計算每個粒子點在第k +1 個采樣時刻第j 個運動模型的的似然函數:

通過公式wnj(k +1)=wnj(k)Lnj(k +1)計算權重，規則化權重

4)更新每個運動模型的概率。

5)加權求和。

在IMMPF 的第3)步計算過程中，利用中心差分擴展卡爾曼濾波算法對基本粒子濾波進行建議分布構造，實現對目標運動狀態的更新，可大幅減小狀態更新時間，且跟蹤精度較高。

2.2 模糊推理機

從(24)式可以看出，各個運動模型概率的更新是用運動模型的似然函數計算得到的，而似然函數表示的是在當前采樣時刻運動模型與目標實際運動的一致性程度。因此，運動模型概率也同樣反映了在當前采樣時刻運動模型與目標實際運動的一致性程度，模型概率越大，表示這個運動模型與目標實際運動越接近。用新更新得到的運動模型概率μi(k+1)作為模糊推理機的輸入，根據模糊推理機的結果計算第k+1 個采樣時刻運動模型間的轉換概率矩陣Ψij(k +1)，能夠反映當前采樣時刻運動模型與目標實際運動的一致性程度。

為了避免出現隸屬函數值過小的情況，同時考慮到計算各個運動模型的似然函數采用的是高斯概率密度函數的形式，因此，采用高斯形隸屬函數計算每個模型概率的隸屬度。運動模型概率的隸屬函數如圖2所示。圖2中的NB(負大)、PB(正大)、ZO(0)、PS(正小)、PM(正中)表示模糊子集，利用運動模型概率所屬的模糊子集可以確定模型轉換概率矩陣Ψij.

圖2 運動模型概率的隸屬函數Fig.2 Membership functions of motion model probability

模型轉換概率矩陣Ψij是M×M 維矩陣，它能夠混合產生所有運動模型的模型混合概率。為了實現根據目標實際運動情況實時調整模型轉換概率矩陣，設計模糊推理機的運行規則，設Ψmij表示第m 條規則下模型轉換概率矩陣。假定運動模型個數等于4，則模糊推理機的運行規則如下所示。

規則0:如果μ1是ZO 并且μ2是ZO 并且μ3是ZO 并且μ4是ZO，那么Ψij=Ψ0ij.

規則1:如果μ1是PS 并且μ2是NB 并且μ3是NB 并且μ4是NB，那么Ψij=Ψ1ij.

規則2:如果μ2是PS 并且μ1是NB 并且μ3是NB 并且μ4是NB，那么Ψij=Ψ2ij.

規則3:如果μ3是PS 并且μ1是NB 并且μ2是NB 并且μ4是NB，那么Ψij=Ψ3ij.

規則4:如果μ4是PS 并且μ1是NB 并且μ2是NB 并且μ3是NB，那么Ψij=

當模糊推理機判斷模型概率屬于以上5 條規則之一時，模型轉換概率矩陣Ψij就等于相應的，對的定義如下:

式中:

當模糊推理機判斷模型概率不屬于以上5 條規則之一時，模型轉換概率矩陣Ψij由加權平均法決定:

式中:Am表示第m 條模糊規則在總輸出中所占的權重，按取小法得出Am的值。

式中:∧表示取小的符號;Fmi(μi)表示輸入到模糊推理機中的模型概率的隸屬函數值。

3 仿真實例

選擇單個機動運動的跟蹤問題進行仿真實驗，將FIMMPF 與文獻［10］IMMPF、文獻［11］IMMUPF進行仿真比較，觀察這3 種目標跟蹤算法在精度和實時性方面的性能。

圖3 機動目標跟蹤場景示意圖Fig.3 Schematic diagram of maneuvering target tracking

圖4 測量得到的方位角變化曲線Fig.4 Changing curve of measured azimuth

算法參數設置如下:對于交互式多模型算法，主要是運動模型的選擇和設置，在本次仿真實驗中，選擇3 種運動模型，3 種運動模型的角速度分別設置為ω1=0 rad/s，ω2=0.115 rad/s，ω3= -0.046 rad/s.3 種運動模型的初始模型概率為μ1=μ2=μ3=1/3;模型間的初始轉換概率矩陣都是Ψij=［0.33，0.33，0.33;0.33，0.33，0.33;0.33，0.33，0.33］. 設置4 個模糊子集，這4 個模糊子集分別為NB(負大)、ZO(0)、PS(正小)、PB(正大)，它們的中心分別為0，1/3，2/3，1;模糊函數選擇高斯形隸屬函數。在改進的粒子濾波算法中，應用了中心差分擴展卡爾曼濾波算法，取中心差分的步長無跡粒子濾波算法的粒子數取為N =100. 分別用3 種目標跟蹤算法對其進行跟蹤，蒙特卡羅仿真次數為100，將它們的相對距離誤差和相對速度誤差曲線進行比較，如圖5～圖8所示。

圖9～圖11 分別是用3 種算法跟蹤時模型概率的變化曲線圖。

圖5 目標與探測器水平相對距離誤差Fig.5 Horizontal relative distance estimation error between target and detector

圖6 目標與探測器垂直相對距離誤差Fig.6 Veritcal relative distance estimation error between target and detector

圖7 目標與探測器水平相對速度誤差Fig.7 Horizontal relative velocity estimation error between target and detector

表1給出了3 種跟蹤算法經過100 次蒙特卡羅仿真后得到的相對距離，相對速度的均方根誤差RMSE 和運算時間。其中表示目標與探測器之間距離的均方根誤差分別為目標與探測器之間在水平和垂直距離上的均方根誤差，t 表示運行100 次蒙特卡羅仿真的時間，t'表示在一次完整的仿真中，每完成一次狀態向量更新所使用的時間。

圖8 目標與探測器垂直相對速度誤差Fig.8 Veritcal relative velocity estimation error between target and detector

圖9 IMMPF 中3 種模型概率變化曲線Fig.9 The probability curves of three models of IMMPF

圖10 IMMUPF 中3 種模型概率變化曲線Fig.10 The probability curves of three models of IMMUPF

圖11 FIMMPF 中3 種模型概率變化曲線Fig.11 The probability curves of three models of FIMMPF

從3 種目標跟蹤算法得到的相對距離誤差曲線圖、相對速度誤差曲線圖、3 種模型概率變化曲線圖和表1中的數據可以看出，由于模糊控制方法在每一次狀態更新時，實時調整模型轉換概率矩陣，使與目標當前運動狀態最接近的運動模型在混合產生這一采樣時刻的初始狀態向量里占有更大的比重，因此，FIMMPF 的目標跟蹤精度最高，其次是IMMUPF，而IMMPF 的精度最低。從計算時間上，IMMUPF 在完成每一次目標狀態向量更新的時間是最長的，其次是IMMPF，而FIMMPF 所需的時間最短。這是因為中心差分擴展卡爾曼濾波算法只在計算殘差協方差矩陣時用中心差分公式完成非線性觀測方程傳遞函數偏導數的計算，而其他部分與擴展卡爾曼濾波算法相同，因此，中心差分擴展卡爾曼濾波算法的運算時間很短。

4 結論

研究了基于旋轉式三維靜電探測器的機動目標跟蹤問題，建立了探測低空飛行目標的數學模型。針對測量方程為非線性的特點和機動目標跟蹤問題，為了提高濾波算法的精度，減少計算時間，提出了一種新的適用于靜電機動目標跟蹤的FIMMPF.理論分析與仿真結果表明，FIMMPF 比IMMUPF 和IMMPF 更適用于機動運動的靜電目標跟蹤。

References)

［1］郝曉輝，虞健飛，崔占忠. 直升機靜電場研究［J］. 兵工學報，2012，33(5):583 -587 HAO Xiao-hui，YU Jian-fei，CUI Zhan-zhong. Research on electrostatic field of helicopter［J］. Acta Armamentarii，2012，33(5):583 -587. (in Chinese)

［2］畢軍建. 對空彈藥引信用靜電矢量探測技術研究［D］. 北京:北京理工大學，2003.BI Jun-jian. Electrostatic vector probing technology research of air ammunition fuze［D］. Beijing:Beijing Institute of Technology，2003. (in Chinese)

［3］代方震. 基于電極掃描原理的被動式靜電探測技術研究［D］.北京:北京理工大學，2004.DAI Fang-zhen. Based on the principle of passive electrode scanning electrostatic detection technology research［D］. Beijing:Beijing Institute of Technology，2004. (in Chinese)

［4］陳曦.被動式地面靜電探測技術研究［D］. 北京:北京理工大學，2005.CHEN Xi. Research on the ground passive electrostatic detection technology［D］. Beijing:Beijing Institute of Technology，2005.(in Chinese)

［5］付巍.低空目標靜電探測與跟蹤算法［D］. 北京:北京理工大學，2009.FU Wei. The low altitude target detection and tracking based on electrostatic detectors［D］. Beijing:Beijing Institute of Technology，2009. (in Chinese)

［6］Boers Y，Driessen J N. Interacting multiple model particle filter［J］. IEEE Proceedings of Radar Sonar and Navigation，2003，150(5):334 -349.

［7］Sanjeev M，Simon M. A tutorial on particle filter for online nonlinear/non-Gaussian Bayesian tracking［J］. IEEE Trans on Signal Processing，2002，50(2):174 -188.

［8］Hu Z T，Pan Q. An improved Particle filtering algorithm based on observation inversion optimal sampling［J］. Journal of Central South University of Technology，2009，16(5):815 -820.

［9］劉尚合，武占成. 靜電放電及危害防護［M］. 北京:北京郵電大學出版社，2004:35 -36.LIU Shang-he，WU Zhan-cheng. Electrostatic and discharge and damage protection［M］. Beijing:Beijing University of Posts and Telecommunications Press，2004:35 -36. (in Chinese)

［10］曹陽，趙明富，羅彬彬，等. 機載空間光通信平臺的交互多模型粒子濾波跟蹤算法［J］. 紅外與激光工程，2012，41(11):3065 -3068.CAO Yang，ZHAO Ming-fu，LUO Bin-bin，et al. Airbone platform＇s tracking algorithm for free space optical communication based on IMMPF methods［J］. Infrared and Laser Engineering，2012，41(11):3065 -3068. (in Chinese)

［11］Zahir M，Abdelaziz O. Comparision of interactive multiple model partide filter and interactive multiple model unscented particle filter for tracking multiple manoeuviring targets in sensors array［J］. IEEE Proceedings of Cybernetic Intelligent Systems，2010，41(1):1 -6.