專(zhuān)利名稱：卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法
技術(shù)領(lǐng)域：
本發(fā)明涉及機(jī)動(dòng)目標(biāo)濾波跟蹤問(wèn)題，具體涉及一種卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有 機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法。
背景技術(shù)：
目標(biāo)跟蹤濾波方法在探測(cè)、識(shí)別與跟蹤等多個(gè)領(lǐng)域有著重要現(xiàn)實(shí)意義和廣闊的應(yīng) 用背景，其中機(jī)動(dòng)目標(biāo)跟蹤濾波問(wèn)題更是研究熱點(diǎn)之一。經(jīng)過(guò)眾多學(xué)者的不懈努力提出或 發(fā)展了許多相關(guān)的理論，如卡爾曼理論及其擴(kuò)展方法、小波分析、粒子濾波、多模型概率估 計(jì)等等，其中以卡爾曼理論及其擴(kuò)展的濾波方法最為人們接收，應(yīng)用也是最為廣泛，但它們 和許多其他方法一樣，仍存在著對(duì)參數(shù)依賴性強(qiáng)、易發(fā)散等問(wèn)題，從一定程度上限制了這些 方法的應(yīng)用。經(jīng)驗(yàn)?zāi)B(tài)分解(EmpiricalMode Decomposition-EMD)方法是于 1998 年由美 國(guó)國(guó)家工程學(xué)院院士黃鍔提出，該方法是將原信號(hào)分解為從高頻到低頻的固有模態(tài)函數(shù) (Intrinsic Mode Function-IMF)，且分解出的各個(gè)IMF的頻率高低是由信號(hào)本身的頻率特 性決定的，具有自適應(yīng)的特性，因此適用于噪聲形式未知的非平穩(wěn)信號(hào)濾波。在應(yīng)用EMD方 法時(shí)，有一個(gè)無(wú)法回避的問(wèn)題，即在用極點(diǎn)擬合上下包絡(luò)的時(shí)候，在兩個(gè)端點(diǎn)附近就會(huì)出現(xiàn) 包絡(luò)失真問(wèn)題，造成計(jì)算上的不準(zhǔn)確，這種現(xiàn)象被稱為“邊界效應(yīng)”。而且隨著求取IMF階次 的升高，迭代的次數(shù)不斷增加，“邊界效應(yīng)”還會(huì)傳播到信號(hào)內(nèi)部，進(jìn)一步影響EMD分解的質(zhì) 量，使分解出的IMF失去了實(shí)際的物理意義。而在對(duì)目標(biāo)進(jìn)行跟蹤濾波時(shí)所關(guān)心的是目標(biāo) 最新的狀態(tài)，恰恰是信號(hào)序列的邊界值，正是這個(gè)原因限制了 EMD在機(jī)動(dòng)目標(biāo)跟蹤濾波中 的應(yīng)用。
發(fā)明內(nèi)容
本發(fā)明的目的是克服EMD方法中的“邊界效應(yīng)”以及強(qiáng)機(jī)動(dòng)目標(biāo)卡爾曼濾波的發(fā) 散問(wèn)題，提供了一種卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法。卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法，它的具體過(guò)程如下步驟一、獲得當(dāng)前采樣時(shí)刻的量測(cè)數(shù)據(jù)序列，并利用卡爾曼預(yù)測(cè)方程，獲得下一采 樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)，并將所述下一采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)和當(dāng)前量測(cè)數(shù)據(jù) 序列組合，獲得組合后的數(shù)據(jù)序列；步驟二、在步驟一獲得的組合后的數(shù)據(jù)序列中的每相鄰兩個(gè)數(shù)據(jù)之間均勻分布多 個(gè)插值點(diǎn)，生成濾波序列；步驟三、利用EMD方法對(duì)步驟二生成的濾波序列進(jìn)行分解，獲得包含噪聲的IMF， 在濾波序列中剔除包含噪聲的IMF，然后獲得當(dāng)前時(shí)刻的濾波值；步驟四、將當(dāng)前時(shí)刻的濾波值作為當(dāng)前濾波結(jié)果顯示；步驟五、將當(dāng)前時(shí)刻的濾波值作為系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)，在下一采樣時(shí)刻來(lái) 臨時(shí)，結(jié)合下一采樣時(shí)刻的量測(cè)數(shù)據(jù)，并利用卡爾曼方程組得到用于下一采樣時(shí)刻的計(jì)算
5參數(shù)，然后，返回執(zhí)行步驟一，將下一采樣時(shí)刻作為當(dāng)前采樣時(shí)刻，實(shí)現(xiàn)機(jī)動(dòng)目標(biāo)的跟蹤。本發(fā)明的積極效果本發(fā)明有效地克服了 EMD方法中“邊界效應(yīng)”以及針對(duì)強(qiáng)機(jī)動(dòng)目標(biāo)卡爾曼濾波的 發(fā)散問(wèn)題，具有較高的跟蹤精度。
圖1為本發(fā)明的流程圖；圖2為步驟三的流程圖；圖3為步驟三一的流程圖；圖4 為應(yīng)用EKF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系X軸向的目標(biāo)速度信息的濾波效果圖； 圖5為應(yīng)用KF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系X軸向的目標(biāo)速度信息的濾波效果 圖；圖6為應(yīng)用EMD-KF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系X軸向的目標(biāo)速度信息的 濾波效果圖；圖7為應(yīng)用EKF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系Y軸向的目標(biāo)速度信 息的濾波效果圖；圖8為應(yīng)用KF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系Y軸向的目標(biāo)速 度信息的濾波效果圖；圖9為應(yīng)用EMD-KF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系Y軸向 的目標(biāo)速度信息的濾波效果圖；圖10為應(yīng)用EKF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐標(biāo)系 Z軸向的目標(biāo)速度信息的濾波效果圖；圖11為應(yīng)用KF方法在30秒時(shí)間段內(nèi)測(cè)得的地理坐 標(biāo)系Z軸向的目標(biāo)速度信息的濾波效果圖；圖12為應(yīng)用EMD-KF方法在30秒時(shí)間段內(nèi)測(cè)得 的地理坐標(biāo)系Z軸向的目標(biāo)速度信息的濾波效果圖。
具體實(shí)施例方式具體實(shí)施方式
一本實(shí)施方式的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo) 跟蹤方法，它的具體過(guò)程如下步驟一、獲得當(dāng)前采樣時(shí)刻的量測(cè)數(shù)據(jù)序列，并利用卡爾曼預(yù)測(cè)方程，獲得下一采 樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)，并將所述下一采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)和當(dāng)前量測(cè)數(shù)據(jù) 序列組合，獲得組合后的數(shù)據(jù)序列；步驟二、在步驟一獲得的組合后的數(shù)據(jù)序列中的每相鄰兩個(gè)數(shù)據(jù)之間均勻分布多 個(gè)插值點(diǎn)，生成濾波序列；步驟三、利用EMD方法對(duì)步驟二生成的濾波序列進(jìn)行分解，獲得包含噪聲的IMF， 在濾波序列中剔除包含噪聲的IMF，然后獲得當(dāng)前時(shí)刻的濾波值；步驟四、將當(dāng)前時(shí)刻的濾波值作為當(dāng)前濾波結(jié)果顯示；步驟五、將當(dāng)前時(shí)刻的濾波值作為系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)，在下一采樣時(shí)刻來(lái) 臨時(shí)，結(jié)合下一采樣時(shí)刻的量測(cè)數(shù)據(jù)，并利用卡爾曼方程組得到用于下一采樣時(shí)刻的計(jì)算 參數(shù)，然后，返回執(zhí)行步驟一，將下一采樣時(shí)刻作為當(dāng)前采樣時(shí)刻，實(shí)現(xiàn)機(jī)動(dòng)目標(biāo)的跟蹤。在系統(tǒng)運(yùn)行過(guò)程中，若接收到來(lái)自外界的結(jié)束控制命令，則結(jié)束當(dāng)前計(jì)算過(guò)程。本發(fā)明綜合運(yùn)用了卡爾曼濾波(Kalman Filter-KF)理論和經(jīng)驗(yàn)?zāi)B(tài)分解的方法， 利用卡爾曼濾波方法對(duì)EMD濾波序列進(jìn)行端點(diǎn)沿拓，并將EMD的濾波值代入卡爾曼遞推方 程中，有效地克服了 EMD方法中“邊界效應(yīng)”以及針對(duì)強(qiáng)機(jī)動(dòng)目標(biāo)卡爾曼濾波的發(fā)散問(wèn)題， 具有較高的跟蹤精度。本發(fā)明繼承了 EMD方法的自適應(yīng)特性，因此對(duì)目標(biāo)的機(jī)動(dòng)方式、噪聲 的分布形式等條件并不敏感，具有較強(qiáng)的魯棒性。
具體實(shí)施方式
二本實(shí)施方式是對(duì)實(shí)施方式一的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟一所述內(nèi)容的過(guò)程如下在采樣時(shí)刻t(n)，n = 1,2,…，濾波系統(tǒng)獲取當(dāng)前量測(cè)數(shù)據(jù)序列，且在各個(gè)采樣 時(shí)刻所獲取的量測(cè)數(shù)據(jù)序列的長(zhǎng)度均為N ；用{y(i)，i = 1，2，…，N}表示t (n)時(shí)刻的量測(cè)數(shù)據(jù)序列，其中y (N)為t (n)時(shí)刻 采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，y(N-l)為t(n-l)時(shí)刻采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，…， y(l)為t(n-N+l)時(shí)刻采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)；所述卡爾曼預(yù)測(cè)方程為 其中，n為當(dāng)前采樣時(shí)刻的序數(shù)，即當(dāng)前采樣時(shí)刻為t(n)；F(n+1, n)為狀態(tài)轉(zhuǎn)移矩陣，它為已知量；| 為用于預(yù)測(cè)的輸入量， iO + lIn)為預(yù)測(cè)結(jié)果輸出量；在首次執(zhí)行步驟一時(shí)，令| 77)等于系統(tǒng)當(dāng)前時(shí)刻的系統(tǒng)狀態(tài)測(cè)量數(shù)據(jù)y(N)； 在非首次執(zhí)行步驟一時(shí)，令iO I ？2)等于上次執(zhí)行步驟五時(shí)獲得的所述用于下一采樣時(shí)刻 的計(jì)算參數(shù)；通過(guò)卡爾曼預(yù)測(cè)方程，并根據(jù)I n)，獲得下一采樣時(shí)刻t (n+l)的系統(tǒng)狀態(tài)預(yù) 測(cè)數(shù)據(jù)+ ；將所述下一時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)義(W + 1 |作)和原始量測(cè)數(shù)據(jù)序列組合，令 {y⑴，y⑵，...，y(N)}表示原始量測(cè)數(shù)據(jù)序列，則組合后的數(shù)據(jù)序列為具體實(shí)施方式
三本實(shí)施方式是對(duì)實(shí)施方式一或二的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解 有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟二所述內(nèi)容的具體過(guò)程為令
表示步驟一獲得的組合后的數(shù)據(jù)序列， 其中y(l)、y(2)、……y(N)為步驟一獲得的當(dāng)前量測(cè)數(shù)據(jù)序列{y (i)，i = 1，2，…，N}中 的N個(gè)系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)+ 11…為步驟一獲得的下一采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)；在
中的每相鄰兩個(gè)數(shù)據(jù)之間分別均勻分布j 個(gè)插值點(diǎn)，j為正整數(shù)，生成濾波序列{X⑴，i = 1,2,…，nX(j+l)+l}。
具體實(shí)施方式
四本實(shí)施方式是對(duì)實(shí)施方式一至三中的任意一種卡爾曼濾波與經(jīng) 驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟三所述內(nèi)容的具體過(guò)程為步驟三一、令{x(i)，i = 1，2，…，nX (j+l)+l}表示步驟二獲得的濾波序列，利用 EMD方法對(duì)該濾波序列進(jìn)行分解，得到多個(gè)IMF，然后在所述多個(gè)IMF中獲得含有噪聲信號(hào) 的 IMF ；步驟三二、在濾波序列{x⑴，i = 1,2,…，nX (j+l)+l}中減去所述含有噪聲信 號(hào)的IMF，進(jìn)而獲得當(dāng)前時(shí)刻的濾波值。IMF為固有模態(tài)函數(shù)，滿足以下兩個(gè)條件的函數(shù)稱為固有模態(tài)函數(shù)條件一過(guò)零點(diǎn)和極值點(diǎn)個(gè)數(shù)相差不超過(guò)一；條件二 上下包絡(luò)線關(guān)于x軸對(duì)稱。
其中，條件二不可能完全滿足，只能近似滿足。
具體實(shí)施方式
五本實(shí)施方式是對(duì)實(shí)施方式四的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī) 結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟三一所述內(nèi)容的具體過(guò)程為步驟三一一、令{x(i)，i = 1，2，…，nX (j+l)+l}表示步驟二獲得的濾波序列，并 將該濾波序列作為當(dāng)前的待變換信號(hào)；初始化1 = 0;步驟三一二、令1 = 1+1 ；獲取當(dāng)前的待變換信號(hào)的極大值點(diǎn)序列和極小值點(diǎn)序 列，通過(guò)擬合獲得所述極大值點(diǎn)序列的包絡(luò)曲線以及所述極小值點(diǎn)序列的包絡(luò)曲線；步驟三一三、獲得極大值點(diǎn)序列的包絡(luò)曲線與極小值點(diǎn)序列的包絡(luò)曲線的均值曲 線，然后開(kāi)始進(jìn)行迭代計(jì)算，并初始化&為待變換信號(hào)，所述迭代計(jì)算的公式為hk = hk_rmk, k = 1,2,...其中，hk為第k次迭代計(jì)算的結(jié)果，hn為第k_l次迭代計(jì)算的結(jié)果，k為迭代次 數(shù)，mk為、的極大值點(diǎn)序列的包絡(luò)曲線與、的極小值點(diǎn)序列的包絡(luò)曲線的均值曲線；在所述迭代計(jì)算的過(guò)程中，當(dāng)hk滿足IMF的定義條件時(shí)即停止迭代，并令Cl等于 當(dāng)前的迭代計(jì)算的結(jié)果，其中Cl為原始濾波序列的第1個(gè)IMF分量；步驟三一四、計(jì)算r1: 步驟三一五、判斷ri的極值點(diǎn)個(gè)數(shù)是否小于或等于2 若是，則執(zhí)行步驟三一六； 否則，返回執(zhí)行步驟三一二；步驟三一六、當(dāng)前巧即為信號(hào)殘余；在1個(gè)IMF分量Cl、c2. . . . Cl中，獲得含有噪 聲信號(hào)的IMF分量。
具體實(shí)施方式
六本實(shí)施方式是對(duì)實(shí)施方式五的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī) 結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟三一六中所述的在1個(gè)IMF分量Cl、c2. . . . Cl 中，獲得含有噪聲信號(hào)的IMF分量的具體過(guò)程為步驟A1、 步驟A2、計(jì)算當(dāng)前Dif(p)，并判斷當(dāng)前Dif(p)所包含的極點(diǎn)個(gè)數(shù)是否小于預(yù)設(shè)的 閾值若是，則執(zhí)行步驟A4 ；否則，執(zhí)行步驟A3 ；步驟A3、令p = p+1，然后返回步驟A2 ；步驟A4、在1個(gè)IMF分量ci、c2. . . . Cl中，含有噪聲信號(hào)的IMF分量為Cl、c2. . . . cp。需要指出的是，由于噪聲分布和目標(biāo)機(jī)動(dòng)形式的不同，包含有噪聲信號(hào)的IMF的 階次也是不同的。甚至對(duì)于噪聲分布和機(jī)動(dòng)形式相同的目標(biāo)，由于采樣數(shù)據(jù)長(zhǎng)度的限制，包 含有噪聲信號(hào)的IMF的階次也是在不斷變化的，因此在濾波過(guò)程中需要引入決策機(jī)制加以 判斷。
具體實(shí)施方式
七本實(shí)施方式是對(duì)實(shí)施方式五的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī) 結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟三一六中所述的在1個(gè)IMF分量Cl、c2. . . . Cl 中，獲得含有噪聲信號(hào)的IMF分量的具體過(guò)程為
分別依次判斷1個(gè)IMF分量C1、C2. ...Cl的信號(hào)能量是否低于預(yù)設(shè)閾值，并將信號(hào) 能量低于預(yù)設(shè)閾值的IMF分量作為含有噪聲信號(hào)的IMF。
具體實(shí)施方式
八本實(shí)施方式是對(duì)實(shí)施方式四的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī) 結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟三二所述內(nèi)容的具體過(guò)程為在濾波序列{x⑴，i = 1,2,…，nX (j+l)+l}中減去所述含有噪聲信號(hào)的IMF, 獲得一個(gè)新序列，則該新序列的第j+2個(gè)點(diǎn)的值即為當(dāng)前時(shí)刻的濾波值。
具體實(shí)施方式
九本實(shí)施方式是對(duì)實(shí)施方式一至八中的任意一種卡爾曼濾波與經(jīng) 驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法的進(jìn)一步說(shuō)明，步驟五所述內(nèi)容的具體過(guò)程為令t(n)表示當(dāng)前時(shí)刻，則將t(n)時(shí)刻的濾波值作為系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)代 入如下的卡爾曼方程組中，即將步驟三中所得到的濾波值替換一下方程組中的1 n)，并 結(jié)合t(n+l)時(shí)刻的測(cè)量值z(mì)(n+l)進(jìn)行卡爾曼濾波估計(jì) x{n + \\n + \) = x{n + \\ n) + K{n + l)[z{n +1) - H(n + \)x{n +11 n)]K (n+l) = P (n+11 n) HT (n+1) [H (n+1) P (n+11 n) HT (n+1) +R (n+1) ]P (n+11 n) = F (n+1，n) P (n | n) FT (n+1，n) +G (n+1，n) Q (n) GT (n+1，n)P (n+11 n+1) = P (n+11 n) _K (n+1) H (n+1) P (n+11 n)= [I-K (n+l)H(n+l)]P(n+l |n)即可獲得用于下一采樣時(shí)刻t (n+1)的計(jì)算參數(shù):+11 +1)，該參數(shù)用于求取 t(n+l)時(shí)刻的預(yù)測(cè)值i(> + 2| + l)。以下結(jié)合數(shù)字仿真算例闡述本發(fā)明的
具體實(shí)施例方式數(shù)字仿真是評(píng)判濾波跟蹤算法性能的重要測(cè)試手段。本仿真環(huán)境及相關(guān)參數(shù)來(lái)自 某雷達(dá)對(duì)運(yùn)動(dòng)目標(biāo)的試驗(yàn)?zāi)M數(shù)據(jù)。運(yùn)動(dòng)目標(biāo)作錐形螺旋機(jī)動(dòng)。圓柱螺旋線的數(shù)學(xué)方程為 其中，xt、yt和zt分別為運(yùn)動(dòng)目標(biāo)的空間三維坐標(biāo)，t為時(shí)間變量，I為目標(biāo)的螺 旋機(jī)動(dòng)半徑，《為其螺旋機(jī)動(dòng)的角速度，\為螺旋線軸線方向的速度。由于目標(biāo)的機(jī)動(dòng)能力和大氣密度、自身速度等因素有關(guān)，所以在不同高度的最大 螺旋加速度a_是不同的，數(shù)值由下式給出 其中，Ht為目標(biāo)所處高度，g為重力加速度。
本實(shí)施例中目標(biāo)螺旋機(jī)動(dòng)角速度ω = 3rad/S，初始螺旋線軸線方向的速度λ ^ = 500m/s, ξ由aMX= ω2ξ實(shí)時(shí)求出，彈道傾角在-50° -80°變化，雷達(dá)可以獲取地理坐標(biāo) 系下目標(biāo)的三維速度信息，且都含有均值為0，方差為10的隨機(jī)噪聲，采樣周期0. 1秒。仿 真中將EMD-KF方法與當(dāng)前被人們廣為接受的卡爾曼方法及擴(kuò)展卡爾曼方法的濾波效果做 了對(duì)比，采用標(biāo)準(zhǔn)差作為評(píng)價(jià)指標(biāo)，使得結(jié)果更加具有說(shuō)服力。在在采樣時(shí)刻t(n)開(kāi)始跟蹤濾波，其中η = 1,2,…；將t (η)時(shí)刻的 測(cè)量值和之前η-1個(gè)采樣時(shí)刻的測(cè)量值按時(shí)間順序排成數(shù)據(jù)序列{y(i)，i = 1， 2，…，N}，其中y(N)為t(n)時(shí)刻采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，y (N-I)為t(n_l) 時(shí)刻采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，…，y(l)為t(n-N+l)時(shí)刻采樣獲得的系統(tǒng) 狀態(tài)量測(cè)數(shù)據(jù)。根據(jù)測(cè)量值y(N)以及設(shè)定的狀態(tài)轉(zhuǎn)移矩陣F(n+1，η),由卡爾 曼理論的預(yù)測(cè)方程xO + 11 w) = F(n + l,n)x(n 可以推斷出下一時(shí) 刻系統(tǒng)的狀態(tài)+ …，^> + 11打)與之前的量測(cè)數(shù)據(jù)構(gòu)成新的數(shù)據(jù)序列 {少(1)，少(2)，·..,ΛΝ),+11 )}，此時(shí)數(shù)據(jù)序列長(zhǎng)度為n+l，本實(shí)施例中，η = 5。在新數(shù)據(jù)序列的每相鄰兩個(gè)數(shù)據(jù)之間均勻分布j個(gè)插值點(diǎn)，仿真中采用線性插 值，j = 10，構(gòu)成了可用于濾波的序列Ix⑴，i = 1,2,…，nX (j+l)+l}。應(yīng)用EMD方法對(duì)序列|x(i)，i = 1，2，…，nX (j+l)+l}進(jìn)行分解，其分解過(guò)程如 下a)、通過(guò)下式迭代得到第一個(gè)IMF，公式中的x(t)是需要被變換的原信號(hào)，此處 x(t)等于{x(i)，i = 1,2,…，nX(j+l)+l}，獲取該信號(hào)的極大值點(diǎn)和極小值點(diǎn)的序列，用 擬合的方法構(gòu)造極值序列的包絡(luò)曲線，Hl1是它的上下包絡(luò)線的均值，其余Hli是?。。、迳舷掳?絡(luò)線的均值。 hk滿足IMF的定義即可停止本次迭代，得到第一個(gè)IMF分量Cl。b)、令Γι = x(t)_Cl，對(duì)Γι重復(fù)上述迭代過(guò)程依次得到C2, C3- cn,即 X(I)-Cl^rv c)、當(dāng)&的極值點(diǎn)小于或等于2時(shí)，則停止對(duì)χ (t)的分解。采用頻率判斷的方法，用原信號(hào)減去含有噪聲的IMF后，得到最終的濾波信號(hào)，該信號(hào)序列中第j+2個(gè)點(diǎn)的濾波值即為當(dāng)前時(shí)刻的濾波值。將步驟三中所得到的濾波值作為 系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)趙《| )，并結(jié)合t(n+l)時(shí)刻的測(cè)量數(shù)據(jù)z(n+l)進(jìn)行卡爾曼濾波 估計(jì)χ(η +11= F(η +1, )χ{η | ) P (η+11 η+1) = P (η+11 η) -K (η+1) H (η+1) P (η+11 η)= [I-K(η+1)H(η+1)]P(η+1 |η)經(jīng)過(guò)迭代運(yùn)算，可以得到參數(shù)+11 W +1)， 在t(n+l)采樣時(shí)刻，用η+1代替η 代入+ l丨“)=+| “)中，即通過(guò) x(n + 2\n + l) = F(n + 2,n + l)x(n +1丨“+1)獲得t (n+2)采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)
數(shù)據(jù)+ 2\n + \)，之后反復(fù)進(jìn)行上述步驟，直到跟蹤濾波過(guò)程的結(jié)束。為了體現(xiàn)本發(fā)明所提設(shè)計(jì)方法的有效性和優(yōu)越性，經(jīng)過(guò)Monte Carlo隨機(jī)仿真，采 用標(biāo)準(zhǔn)差作為評(píng)價(jià)指標(biāo)，將EMD-KF方法的濾波效果與當(dāng)前被人們廣為接受的卡爾曼方法 (KF)及擴(kuò)展卡爾曼方法(EKF)做了對(duì)比。表1為三種方法的濾波效果比較。表 1 應(yīng)用EKF方法，在30秒時(shí)間段內(nèi)，在地理坐標(biāo)系X軸向、Y軸向和Z軸向獲得的目 標(biāo)速度信息的濾波效果圖分別為圖4、圖7和圖10 ；應(yīng)用KF方法，在30秒時(shí)間段內(nèi)，在地理坐標(biāo)系X軸向、Y軸向和Z軸向獲得的目 標(biāo)速度信息的濾波效果圖分別為圖5、圖8和圖11 ；應(yīng)用EMD-KF方法，在30秒時(shí)間段內(nèi)，在地理坐標(biāo)系X軸向、Y軸向和Z軸向獲得 的目標(biāo)速度信息的濾波效果圖分別為圖6、圖9和圖12 ；從表1及圖4至圖12中可以看出，卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)合的機(jī)動(dòng)目標(biāo)跟 蹤方法(EMD-KF)的濾波效果要優(yōu)于KF方法，與EKF相當(dāng)。而且由于EMD濾波的自適應(yīng)特 性，卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)合的濾波方法可以不需要事先獲知噪聲的統(tǒng)計(jì)特性，對(duì) 噪聲的分布、目標(biāo)的機(jī)動(dòng)形式也并不敏感，因此比EKF具有更好的適用性。另外，隨著模型 誤差或目標(biāo)機(jī)動(dòng)能力的進(jìn)一步加大，KF和EKF的濾波效果都出現(xiàn)了明顯的降低甚至發(fā)散， 而卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)合濾波方法的效果變化不大，體現(xiàn)了較強(qiáng)的魯棒性。
權(quán)利要求
卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法，其特征在于它的具體過(guò)程如下步驟一、獲得當(dāng)前采樣時(shí)刻的量測(cè)數(shù)據(jù)序列，并利用卡爾曼預(yù)測(cè)方程，獲得下一采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)，并將所述下一采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)和當(dāng)前量測(cè)數(shù)據(jù)序列組合，獲得組合后的數(shù)據(jù)序列；步驟二、在步驟一獲得的組合后的數(shù)據(jù)序列中的每相鄰兩個(gè)數(shù)據(jù)之間均勻分布多個(gè)插值點(diǎn)，生成濾波序列；步驟三、利用EMD方法對(duì)步驟二生成的濾波序列進(jìn)行分解，獲得包含噪聲的IMF，在濾波序列中剔除包含噪聲的IMF，然后獲得當(dāng)前時(shí)刻的濾波值；步驟四、將當(dāng)前時(shí)刻的濾波值作為當(dāng)前濾波結(jié)果顯示；步驟五、將當(dāng)前時(shí)刻的濾波值作為系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)，在下一采樣時(shí)刻來(lái)臨時(shí)，結(jié)合下一采樣時(shí)刻的量測(cè)數(shù)據(jù)，并利用卡爾曼方程組得到用于下一采樣時(shí)刻的計(jì)算參數(shù)，然后，返回執(zhí)行步驟一，將下一采樣時(shí)刻作為當(dāng)前采樣時(shí)刻，實(shí)現(xiàn)機(jī)動(dòng)目標(biāo)的跟蹤。
2.根據(jù)權(quán)利要求1所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟一所述內(nèi)容的過(guò)程如下在采樣時(shí)刻t(n)，n = 1,2,…，濾波系統(tǒng)獲取當(dāng)前量測(cè)數(shù)據(jù)序列，且在各個(gè)采樣時(shí)刻 所獲取的量測(cè)數(shù)據(jù)序列的長(zhǎng)度均為N ；用{y⑴，i = 1,2,…，N}表示t(n)時(shí)刻的量測(cè)數(shù)據(jù)序列，其中y (N)為t (n)時(shí)刻采 樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，y(N-l)為t(n-l)時(shí)刻采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，…， y(l)為t(n-N+l)時(shí)刻采樣獲得的系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)； 所述卡爾曼預(yù)測(cè)方程為 關(guān)(n +11 n) = +1, n)x(n I 打)， 其中，n為當(dāng)前采樣時(shí)刻的序數(shù)，即當(dāng)前采樣時(shí)刻為t (n)； F(n+1, n)為狀態(tài)轉(zhuǎn)移矩陣，它為已知量為用于預(yù)測(cè)的輸入量，+ 為預(yù)測(cè)結(jié)果輸出量；在首次執(zhí)行步驟一時(shí)，令| 等于系統(tǒng)當(dāng)前時(shí)刻的系統(tǒng)狀態(tài)測(cè)量數(shù)據(jù)y(N)；在非 首次執(zhí)行步驟一時(shí)，令^> | 等于上次執(zhí)行步驟五時(shí)獲得的所述用于下一采樣時(shí)刻的計(jì) 算參數(shù)；通過(guò)卡爾曼預(yù)測(cè)方程，并根據(jù)| n)，獲得下一采樣時(shí)刻t (n+1)的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù) mx{n + \\n)；將所述下一時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)+ 勸和原始量測(cè)數(shù)據(jù)序列組合，令{y(l)， y(2)，…，y(N)}表示原始量測(cè)數(shù)據(jù)序列，則組合后的數(shù)據(jù)序列為
3.根據(jù)權(quán)利要求1所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟二所述內(nèi)容的具體過(guò)程為令{><1),少(2),...,少(#),^> + 1| )}表示步驟一獲得的組合后的數(shù)據(jù)序列，其中 y(l)、y(2)、……y(N)為步驟一獲得的當(dāng)前量測(cè)數(shù)據(jù)序列{y (i)，i = 1，2，…，N}中的N個(gè)系統(tǒng)狀態(tài)量測(cè)數(shù)據(jù)，^> + n)為步驟一獲得的下一采樣時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)； 在{少(1),少(2),...,少0)，^> + 11")}中的每相鄰兩個(gè)數(shù)據(jù)之間分別均勻分布j個(gè)插值點(diǎn)，j為正整數(shù)，生成濾波序列{X⑴，i = 1,2,…，nX(j+l)+l}。
4.根據(jù)權(quán)利要求1所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟三所述內(nèi)容的具體過(guò)程為步驟三一、令{x(i)，i = 1,2,…，nX (j+l)+l}表示步驟二獲得的濾波序列，利用EMD 方法對(duì)該濾波序列進(jìn)行分解，得到多個(gè)IMF，然后在所述多個(gè)IMF中獲得含有噪聲信號(hào)的 IMF ；步驟三二、在濾波序列Ix⑴，i = 1,2,…，nX (j+l)+l}中減去所述含有噪聲信號(hào)的 IMF，進(jìn)而獲得當(dāng)前時(shí)刻的濾波值。
5.根據(jù)權(quán)利要求4所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟三一所述內(nèi)容的具體過(guò)程為步驟三一一、令{x⑴，i = 1,2,…，nX(j+l)+l}表示步驟二獲得的濾波序列，并將該 濾波序列作為當(dāng)前的待變換信號(hào)；初始化1=0;步驟三一二、令1 = 1+1 ；獲取當(dāng)前的待變換信號(hào)的極大值點(diǎn)序列和極小值點(diǎn)序列，通 過(guò)擬合獲得所述極大值點(diǎn)序列的包絡(luò)曲線以及所述極小值點(diǎn)序列的包絡(luò)曲線；步驟三一三、獲得極大值點(diǎn)序列的包絡(luò)曲線與極小值點(diǎn)序列的包絡(luò)曲線的均值曲線 然后開(kāi)始進(jìn)行迭代計(jì)算，并初始化&為待變換信號(hào)，所述迭代計(jì)算的公式為hk = hH-nik, k = 1,2,…其中，hk為第k次迭代計(jì)算的結(jié)果，hn為第k-1次迭代計(jì)算的結(jié)果，k為迭代次數(shù)，mk 為、的極大值點(diǎn)序列的包絡(luò)曲線與、的極小值點(diǎn)序列的包絡(luò)曲線的均值曲線；在所述迭代計(jì)算的過(guò)程中，當(dāng)hk滿足IMF的定義條件時(shí)即停止迭代，并令Cl等于當(dāng)前 的迭代計(jì)算的結(jié)果，其中Cl為原始濾波序列的第1個(gè)IMF分量； 步驟三一四、計(jì)算r1: 步驟三一五、判斷巧的極值點(diǎn)個(gè)數(shù)是否小于或等于2 若是，則執(zhí)行步驟三一六；否則， 返回執(zhí)行步驟三一二；步驟三一六、當(dāng)前巧即為信號(hào)殘余；在1個(gè)IMF分量Cl、c2. . . . Cl中，獲得含有噪聲信 號(hào)的IMF分量。
6.根據(jù)權(quán)利要求5所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟三一六中所述的在1個(gè)IMF分量Cl、(v"Cl中，獲得含有噪聲信號(hào)的IMF分 量的具體過(guò)程為 步驟A2、計(jì)算當(dāng)前Dif(p)，并判斷當(dāng)前Dif(p)所包含的極點(diǎn)個(gè)數(shù)是否小于預(yù)設(shè)的閾 值若是，則執(zhí)行步驟A4 ；否則，執(zhí)行步驟A3 ；步驟A3、令p = p+1，然后返回步驟A2 ；步驟A4、在1個(gè)IMF分量Cl、c2. . . . Cl中，含有噪聲信號(hào)的IMF分量為Cl、c2. . . . cp。
7.根據(jù)權(quán)利要求5所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟三一六中所述的在1個(gè)IMF分量Cl、c2. . . Cl中，獲得含有噪聲信號(hào)的IMF 分量的具體過(guò)程為分別依次判斷1個(gè)IMF分量ci、c2. ...Cl的信號(hào)能量是否低于預(yù)設(shè)閾值，并將信號(hào)能量 低于預(yù)設(shè)閾值的IMF分量作為含有噪聲信號(hào)的IMF。
8.根據(jù)權(quán)利要求4所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟三二所述內(nèi)容的具體過(guò)程為在濾波序列{x⑴，i = 1,2,…，nX (j+l)+l}中減去所述含有噪聲信號(hào)的IMF，獲得 一個(gè)新序列，則該新序列的第j+2個(gè)點(diǎn)的值即為當(dāng)前時(shí)刻的濾波值。
9.根據(jù)權(quán)利要求1所述的卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法， 其特征在于步驟五所述內(nèi)容的具體過(guò)程為令t(n)表示當(dāng)前時(shí)刻，則將t(n)時(shí)刻的濾波值作為系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)代入如 下的卡爾曼方程組中，即將步驟三中所得到的濾波值替換一下方程組中的| n)，并結(jié)合 t(n+l)時(shí)刻的測(cè)量值z(mì)(n+l)進(jìn)行卡爾曼濾波估計(jì) =[I-K(n+l)H(n+l)]P(n+l|n)即獲得用于下一采樣時(shí)刻t (n+1)的計(jì)算參數(shù) x(n +11 n +1)。
全文摘要
卡爾曼濾波與經(jīng)驗(yàn)?zāi)B(tài)分解有機(jī)結(jié)合的機(jī)動(dòng)目標(biāo)跟蹤方法，涉及機(jī)動(dòng)目標(biāo)濾波跟蹤領(lǐng)域，解決了無(wú)法克服EMD方法中的“邊界效應(yīng)”及強(qiáng)機(jī)動(dòng)目標(biāo)卡爾曼濾波的發(fā)散問(wèn)題。該方法為一、獲得當(dāng)前采樣時(shí)刻的量測(cè)數(shù)據(jù)序列，利用卡爾曼預(yù)測(cè)方程，獲得下一時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)數(shù)據(jù)，將其與當(dāng)前序列組合；二、在組合序列中插值生成濾波序列；三、利用EMD方法對(duì)濾波序列進(jìn)行分解，剔除包含噪聲的IMF，獲得當(dāng)前濾波值；四、將該濾波值作為當(dāng)前濾波結(jié)果顯示；五、將該濾波值作為系統(tǒng)當(dāng)前狀態(tài)的后驗(yàn)估計(jì)，在下一時(shí)刻結(jié)合量測(cè)數(shù)據(jù)及卡爾曼方程組，獲得計(jì)算參數(shù)，然后返回一，將下一時(shí)刻作為當(dāng)前時(shí)刻，實(shí)現(xiàn)機(jī)動(dòng)目標(biāo)的跟蹤。本發(fā)明可用于機(jī)動(dòng)目標(biāo)濾波跟蹤領(lǐng)域。
文檔編號(hào)G06F17/14GK101894097SQ20101023506
公開(kāi)日2010年11月24日 申請(qǐng)日期2010年7月23日 優(yōu)先權(quán)日2010年7月23日
發(fā)明者張淼, 沈毅, 趙振昊 申請(qǐng)人:哈爾濱工業(yè)大學(xué)