付巍，鄭賓

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

0 引言

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

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

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

1 靜電探測系統(tǒng)數(shù)學模型

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

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

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

在旋轉(zhuǎn)式三維靜電探測器3 個方向的感應(yīng)電極上產(chǎn)生的電場強度的三維分量分別為

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

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

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

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

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

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

假設(shè)動態(tài)系統(tǒng)狀態(tài)空間模型為

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

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

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

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

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

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

混合所有運動模型的相應(yīng)的粒子點，則第j 個運動模型的混合狀態(tài)初值及初始狀態(tài)估計協(xié)方差矩陣為

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

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

預(yù)測狀態(tài)向量

預(yù)測協(xié)方差矩陣

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

殘差協(xié)方差矩陣:

令

則

式中:H(k)為觀測方程狀態(tài)轉(zhuǎn)移矩陣;R(k)為觀測方程噪聲的協(xié)方差矩陣。

令

應(yīng)用中心差分理論

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

濾波增益矩陣

濾波輸出

誤差協(xié)方差矩陣

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

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

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

5)加權(quán)求和。

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

2.2 模糊推理機

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

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

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

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

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

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

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

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

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

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

式中:

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

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

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

3 仿真實例

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

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

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

算法參數(shù)設(shè)置如下:對于交互式多模型算法，主要是運動模型的選擇和設(shè)置，在本次仿真實驗中，選擇3 種運動模型，3 種運動模型的角速度分別設(shè)置為ω1=0 rad/s，ω2=0.115 rad/s，ω3= -0.046 rad/s.3 種運動模型的初始模型概率為μ1=μ2=μ3=1/3;模型間的初始轉(zhuǎn)換概率矩陣都是Ψij=［0.33，0.33，0.33;0.33，0.33，0.33;0.33，0.33，0.33］. 設(shè)置4 個模糊子集，這4 個模糊子集分別為NB(負大)、ZO(0)、PS(正小)、PB(正大)，它們的中心分別為0，1/3，2/3，1;模糊函數(shù)選擇高斯形隸屬函數(shù)。在改進的粒子濾波算法中，應(yīng)用了中心差分擴展卡爾曼濾波算法，取中心差分的步長無跡粒子濾波算法的粒子數(shù)取為N =100. 分別用3 種目標跟蹤算法對其進行跟蹤，蒙特卡羅仿真次數(shù)為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 種跟蹤算法經(jīng)過100 次蒙特卡羅仿真后得到的相對距離，相對速度的均方根誤差RMSE 和運算時間。其中表示目標與探測器之間距離的均方根誤差分別為目標與探測器之間在水平和垂直距離上的均方根誤差，t 表示運行100 次蒙特卡羅仿真的時間，t'表示在一次完整的仿真中，每完成一次狀態(tài)向量更新所使用的時間。

圖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中的數(shù)據(jù)可以看出，由于模糊控制方法在每一次狀態(tài)更新時，實時調(diào)整模型轉(zhuǎn)換概率矩陣，使與目標當前運動狀態(tài)最接近的運動模型在混合產(chǎn)生這一采樣時刻的初始狀態(tài)向量里占有更大的比重，因此，F(xiàn)IMMPF 的目標跟蹤精度最高，其次是IMMUPF，而IMMPF 的精度最低。從計算時間上，IMMUPF 在完成每一次目標狀態(tài)向量更新的時間是最長的，其次是IMMPF，而FIMMPF 所需的時間最短。這是因為中心差分擴展卡爾曼濾波算法只在計算殘差協(xié)方差矩陣時用中心差分公式完成非線性觀測方程傳遞函數(shù)偏導(dǎo)數(shù)的計算，而其他部分與擴展卡爾曼濾波算法相同，因此，中心差分擴展卡爾曼濾波算法的運算時間很短。

4 結(jié)論

研究了基于旋轉(zhuǎn)式三維靜電探測器的機動目標跟蹤問題，建立了探測低空飛行目標的數(shù)學模型。針對測量方程為非線性的特點和機動目標跟蹤問題，為了提高濾波算法的精度，減少計算時間，提出了一種新的適用于靜電機動目標跟蹤的FIMMPF.理論分析與仿真結(jié)果表明，F(xiàn)IMMPF 比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］畢軍建. 對空彈藥引信用靜電矢量探測技術(shù)研究［D］. 北京:北京理工大學，2003.BI Jun-jian. Electrostatic vector probing technology research of air ammunition fuze［D］. Beijing:Beijing Institute of Technology，2003. (in Chinese)

［3］代方震. 基于電極掃描原理的被動式靜電探測技術(shù)研究［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］陳曦.被動式地面靜電探測技術(shù)研究［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.