基于本征模函數(shù)能量的外測數(shù)據(jù)誤差分析方法

劉一涵， 黃楷宸， 柏業(yè)超*

(1.南京大學(xué) 電子科學(xué)與工程學(xué)院，江蘇 南京 210023； 2.南京大學(xué) 人工智能學(xué)院，江蘇 南京 210023)

在武器裝備試驗獲得大量的試驗測量數(shù)據(jù)后，如何利用這些數(shù)據(jù)建立數(shù)學(xué)模型，分析測量數(shù)據(jù)誤差并給出可靠的數(shù)據(jù)處理結(jié)果是武器裝備系統(tǒng)和測控系統(tǒng)進(jìn)行精度分析與鑒定的重要環(huán)節(jié)[1]。受外測設(shè)備測量環(huán)境、測量手段、測量元素等因素影響，且導(dǎo)彈飛行瞬態(tài)(如助推器脫落、級間點火和級間分離等)彈道特征點附近噪聲相對較大、大量連續(xù)系統(tǒng)誤差難以識別，使得外彈道測量數(shù)據(jù)誤差復(fù)雜[2]，隨機(jī)誤差大于測量的量化誤差且存在非平穩(wěn)時變特性，可模型化的誤差帶來了具有復(fù)雜特性和不確定性的模型誤差，影響了數(shù)據(jù)處理精度。張峰等[3]采用樣條函數(shù)建立了主動段運動模型，但是對于波動較大的彈道數(shù)據(jù)采用多項式或者樣條函數(shù)描述彈道的方法不能保證軌道參數(shù)的逼近精度，會帶來模型誤差。李蟬等[4]提出了采用自適應(yīng)抗差濾波方法實時計算再入彈道數(shù)據(jù)，該方法適用于非平穩(wěn)測量數(shù)據(jù)處理，但是需要得到精確的測量數(shù)據(jù)和隨機(jī)誤差的統(tǒng)計特性，否則可能導(dǎo)致濾波發(fā)散，得到與實際情況相背離的結(jié)果。姚尚等[5]利用小波變換方法在時域和頻域同時分析數(shù)據(jù)，根據(jù)分解的小波系數(shù)表現(xiàn)出的不同性質(zhì)識別信號和噪聲。李書興[6]研究了經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)方法在外彈道測量數(shù)據(jù)方面的應(yīng)用，但是外彈道測量數(shù)據(jù)中不可避免地會帶有各種噪聲，特別是無線電測量設(shè)備受到環(huán)境電磁干擾時，采用EMD分析測量數(shù)據(jù)，經(jīng)常出現(xiàn)模態(tài)混疊現(xiàn)象，導(dǎo)致分解后的本征模分量時-頻分布混亂，失去實際物理意義。

本文提出了基于本征模態(tài)函數(shù)(Intrinsic Mode Function，IMF)能量拐點的外測數(shù)據(jù)誤差分析方法。首先，采用EMD方法將外測數(shù)據(jù)分解成不同頻率成分的本征模；然后，根據(jù)外測數(shù)據(jù)隨機(jī)誤差、系統(tǒng)誤差和真實數(shù)據(jù)的頻率特征，采用IMF能量拐點方法將分解得到的IMF分成高頻隨機(jī)誤差、混合信息和有效信息共3個集合；接著將高頻隨機(jī)誤差集合的本征模直接濾掉，對混合信息集合的IMF采用本文改進(jìn)的反雙曲正弦函數(shù)閾值函數(shù)作為調(diào)節(jié)因子進(jìn)行小波濾噪，保留有效信息集合；最后，分別采用小波重構(gòu)和EMD重構(gòu)得到濾掉隨機(jī)誤差的外測有效數(shù)據(jù)。本文方法綜合利用了EMD和小波變換方法的優(yōu)點，避免了單一使用一種方法的局限，有效地分離了低頻和次低頻誤差，避免了非平穩(wěn)測量數(shù)據(jù)近似成平穩(wěn)測量數(shù)據(jù)處理過程中引入的模型誤差，適合非線性、多頻段測量數(shù)據(jù)誤差分析。

1 外測數(shù)據(jù)特點

外測測量設(shè)備包括各式雷達(dá)和光電經(jīng)緯儀等。這些設(shè)備的測量元件可分為單脈沖雷達(dá)、相控陣?yán)走_(dá)、紅外傳感器和可見光傳感器等，測量元素主要有距離、角度、速度，測量平臺有固定站、車載平臺、艦載平臺等，因而設(shè)備誤差也是復(fù)雜多樣的。其復(fù)雜性不僅表現(xiàn)在不同的測量體制間，也表現(xiàn)在同一測量體制、不同測量手段的測量數(shù)據(jù)，例如多傳感器光電經(jīng)緯儀，雖然中波紅外傳感器和可見光傳感器都是光學(xué)測量，但是由于安裝的位置和成像的原理等不同，測量數(shù)據(jù)誤差特性也不相同。但是即使采用同一測量手段，不同時間進(jìn)場的測量設(shè)備其數(shù)據(jù)誤差表現(xiàn)也不盡相同。外測數(shù)據(jù)隨機(jī)誤差具有不可觀測性、強(qiáng)相關(guān)性和非平穩(wěn)時變特性。從工程技術(shù)來看，隨機(jī)誤差有時大于測量的量化誤差。從測量時間段看，各個時間區(qū)域方差大小是有變化的。系統(tǒng)誤差只有部分可以通過標(biāo)校測定和修正，這樣的系統(tǒng)誤差數(shù)值較大，即使實際上可以建模的系統(tǒng)誤差采用最可靠的誤差修正模型，在彈道不同測量弧段也存在著不同的誤差修正殘余量，還有些系統(tǒng)誤差不能標(biāo)?；蛘卟荒芙⑿拚Ｐ?，這些系統(tǒng)誤差往往具有潛伏性，不易識別，很難與隨機(jī)誤差分離。

文獻(xiàn)[1]定義了外測測量數(shù)據(jù)的數(shù)學(xué)模型：

y(t)=f(t)+s(t)+e(t)

(1)

式中：y(t)為測量數(shù)據(jù)；f(t)為真實數(shù)據(jù)；s(t)為系統(tǒng)誤差；e(t)為隨機(jī)誤差。系統(tǒng)誤差反映了測量數(shù)據(jù)的準(zhǔn)確度,隨機(jī)誤差反映了測量數(shù)據(jù)的精密度。要獲得高精度的數(shù)據(jù),必須準(zhǔn)確地分離系統(tǒng)誤差和隨機(jī)誤差。

文獻(xiàn)[2]根據(jù)外測數(shù)據(jù)在不同的頻帶中都具有能譜分布，測量數(shù)據(jù)系統(tǒng)誤差一般可分為常值誤差、線性漂移誤差、周期性誤差和復(fù)雜規(guī)律變化誤差等，因此定義測量數(shù)據(jù)模型為

y(t)=P(t)+B(t)+ε(t)

(2)

式中：y(t)為測量數(shù)據(jù)；P(t)為真實數(shù)據(jù)；B(t)為漂移誤差；ε(t)為隨機(jī)誤差。B(t)中所含的頻率成分比真實數(shù)據(jù)P(t)快,但與隨機(jī)誤差ε(t)相比要慢，因此又可稱為次低頻誤差。

2 本文方法

2.1 外測數(shù)據(jù)3個集合{fIMF(i)}與{E(i)}的能量關(guān)系

文獻(xiàn)[2]根據(jù)大量實例認(rèn)為式(2)中P(t)和B(t)之間、B(t)和ε(t)之間通常都存在一個頻帶間隔，提出三頻帶信號分離的樣條函數(shù)方法[2]。該方法分析非平穩(wěn)數(shù)據(jù)會因所選的節(jié)點數(shù)目、位置和待估參數(shù)的增多而導(dǎo)致信號估計的精度不高。

EMD的基本思想為將一個頻率不規(guī)則的波轉(zhuǎn)換為多個單一頻率的波+殘波的形式[7-9]。

EMD算法采用“篩分”的方法將IMF層層篩選出來,將外彈道測量數(shù)據(jù)分解成不同時間尺度的IMF和一個殘差信號之和,詳細(xì)步驟如文獻(xiàn)[6]所示。

外彈道測量數(shù)據(jù)經(jīng)EMD后得到n個IMF，各個IMF分量記作{fIMF(i)}，{fIMF(i)}的能量{E(i)}的計算公式為

(3)

外測數(shù)據(jù)序列的總能量記為E0，其計算公式為

(4)

對各層{fIMF(i)}能量{E(i)}歸一化處理，其計算公式為

(5)

文獻(xiàn)[2]從理論上證明了式(2)中P(t)與B(t)有分頻界點、P(t)+B(t)與ε(t)有分頻界點并且采用分頻方法對運載火箭的外彈道跟蹤測量信號進(jìn)行分析處理,成功地分離了測量信號中的軌道、漂移誤差和隨機(jī)誤差[2]。

根據(jù)外測數(shù)據(jù)特點和測量數(shù)據(jù)的數(shù)學(xué)模型以及文獻(xiàn)[2]，根據(jù)IMF能量集合{E(i)}，將外測數(shù)據(jù)EMD后得到的{fIMF(i)}劃分為高頻隨機(jī)誤差集合A、混合信息集合B和有效信息集合(含有低頻的系統(tǒng)誤差和真實數(shù)據(jù))C，即

A={fIMF(i),i=1,…,m-1}
B={fIMF(i),i=m,…,p-1}
C={fIMF(i),i=p,…,n}

集合A含有高頻的隨機(jī)誤差，其{fIMF(i)}隨著分解層數(shù)的增加，高頻隨機(jī)誤差成分減少，因此fIMF(i)的能量逐漸降低。

集合B含有較高頻的真實數(shù)據(jù)，部分次高頻和較高頻的系統(tǒng)誤差、次高頻和較高頻隨機(jī)誤差，其{fIMF(i)}隨著分解層數(shù)的增加，頻率成分增加。

劃分能量拐點的準(zhǔn)則為

if((EIMF(i+1)-EIMF(i))>0&(EIMF(i+2)-EIMF(i+1))<0),
k=i+1

(6)

2.2 基于IMF能量拐點的測量數(shù)據(jù)成分集合判斷算法

基于IMF能量拐點的測量數(shù)據(jù)成分集合判斷的具體算法如下。

① 輸入外彈道測量序列x(t)，i=1。

② 對x(t)進(jìn)行 EMD分解，得到不同的{fIMF(i)}。

③ 根據(jù)式(3)計算各分量的能量Ei。

④ 根據(jù)式(5)計算各分量歸一化能量εi。

⑤ 計算相鄰分量的能量差，即δ=εi+εi-1。

⑥ 按式(6)尋找分量能量拐點。

2.3 混合信息集合的處理

基于改進(jìn)的EMD和小波變換濾噪算法思想是對測量數(shù)據(jù){x(t)}(t=1,2,…,N)進(jìn)行EMD，采用本文的能量拐點計算方法判定3個集合的臨界點，將高頻噪聲集合直接濾除，有效信息集合保留，混合信息采用db4小波分解，根據(jù)各層小波系數(shù)的標(biāo)準(zhǔn)方差確定濾噪閾值。

采用db4小波對混合信息集合和有限信息集合的各{fIMF(i)}進(jìn)行分解，一般分解層數(shù)為3或4層即可，具體分解過程見文獻(xiàn)[3]。設(shè){fIMF(i)}k時刻在第j層的低頻系數(shù)為ai,k、高頻系數(shù)為di,k[10]。

混合信息中噪聲和有效數(shù)據(jù)沒有明確的界限，采用硬閾值函數(shù)和軟閾值函數(shù)濾噪時，若將絕對值比設(shè)定閾值大的小波系數(shù)直接置為0[11]，會刪除掉一些有用的信息，從而降低信號重構(gòu)后的精度[12]。本文采用反雙曲正弦函數(shù)作為閾值函數(shù)的調(diào)節(jié)因子。

(7)

因此本文的閾值函數(shù)為

(8)

(9)

式中:σj為第j層小波系數(shù)的標(biāo)準(zhǔn)方差。

sinh是奇函數(shù)且關(guān)于原點對稱,值域區(qū)間為(-∞，∞)，是在定義域區(qū)間內(nèi)非線性單調(diào)上升的連續(xù)函數(shù)。采用反雙曲正弦函數(shù)作為閾值函數(shù)調(diào)節(jié)因子，則當(dāng)|dj|=Tj時，Dj=0，Dj連續(xù)；當(dāng)|dj|→∞時，Dj=|dj|。其克服了硬閾值函數(shù)Gibbs振蕩現(xiàn)象，也解決了軟閾值函數(shù)當(dāng)|dj|→∞時Dj存在的與Tj的固定偏差[13-14]，其估計的小波系數(shù)小于真實的小波系數(shù)絕對值，導(dǎo)致重構(gòu)后的信號出現(xiàn)細(xì)節(jié)缺陷問題。

3個集合的處理算法如下。

① 對x(t)進(jìn)行EMD。

② 判定高頻噪聲集合、混合信息集合和有效信息集合。

③ 對混合信息集合的m～p-1的fIMF(i)進(jìn)行db4小波變換。

⑥ 有效信息與x1(t)重構(gòu)后得到去噪后的有效外測數(shù)據(jù)。

3 實驗分析

3.1 分析案例

以某艦載雷達(dá)的一次測量數(shù)據(jù)為例進(jìn)行分析。雷達(dá)固連在船上，測量平臺受海況和船的航行狀況影響，測量平臺有顛簸、搖擺、振動和航行等復(fù)雜運動，測量數(shù)據(jù)是誤差成分復(fù)雜的非平穩(wěn)測量序列。

將北斗測量數(shù)據(jù)經(jīng)坐標(biāo)轉(zhuǎn)換得到待分析艦載雷達(dá)測量數(shù)據(jù)真值[15-16]，將真值和實際測量數(shù)據(jù)作差，得到殘差數(shù)據(jù)，此時殘差數(shù)據(jù)包括測量數(shù)據(jù)的隨機(jī)誤差和系統(tǒng)誤差。然后采用表1中4種方法分離殘差數(shù)據(jù)中系統(tǒng)誤差和隨機(jī)誤差。

表1 實驗方法

按照表1的4種方法，分離了隨機(jī)誤差和系統(tǒng)誤差，用濾掉隨機(jī)誤差的測量數(shù)據(jù)計算目標(biāo)彈道軌跡，并用與北斗定位結(jié)果的彈道差異來分析方法的有效性。

采用本文方法結(jié)合方位角殘差數(shù)據(jù)說明具體分析過程。

圖1為方位角殘差數(shù)據(jù)圖。圖1中數(shù)據(jù)起伏較大，說明數(shù)據(jù)隨機(jī)性強(qiáng)，顯然是非平穩(wěn)數(shù)據(jù)，若視為平穩(wěn)數(shù)據(jù)濾波，必然會影響數(shù)據(jù)精度。

圖1 方位角殘差數(shù)據(jù)圖

圖2為IMF1～I(xiàn)MF9及其頻譜圖。從圖2可以看出，由殘差數(shù)據(jù)經(jīng)EMD得到的IMF在時域上特點不同，看不到與時間有關(guān)的周期性成分存在。IMF1 存在很多高頻信號，且幅值較大。IMF2信號主要集中在頻段0.7～1.8 MHz。IMF1～I(xiàn)MF4低頻分量較多，表現(xiàn)為低頻信號。IMF5～I(xiàn)MF9頻譜中低頻分量的幅值遠(yuǎn)高于高頻分量的幅值。IMF8和IMF9可以視作沒有高頻分量，相較于IMF7低頻幅值更大。

圖3為IMF1～I(xiàn)MF9時域振幅圖。顯然IMF8和IMF9振幅較大，余下的IMF振幅在時域上很難分辨。圖4為IMF1～I(xiàn)MF9頻域幅度圖。IMF信號從高頻到低頻排列，IMF4、IMF7和IMF8呈現(xiàn)明顯的周期性。

將分解得到的IMF系數(shù)代入式(3)和式(5)，計算各層IMF的能量和能量變化，如表2所示，歸一化的各層IMF的能量和EMD分解層數(shù)，如圖5所示。

如果殘差數(shù)據(jù)為白噪聲數(shù)據(jù)則經(jīng)EMD會從第2個IMF開始，IMF能量呈線性關(guān)系，從圖5可以看出，殘差數(shù)據(jù)不是白噪聲數(shù)據(jù)。分解層數(shù)為4和7時是能量曲線的拐點，這與表2中δ4<0,δ7<0，一致，顯然IMF4和IMF7為能量拐點。

圖2 IMF1～I(xiàn)MF9及其頻譜圖

圖3 IMF1～I(xiàn)MF9時域振幅

圖4 IMF1～I(xiàn)MF9頻域幅度圖

表2 IMF能量和能量變化

圖5 IMF能量和EMD分解層數(shù)圖

高頻隨機(jī)誤差集合={IMF1，IMF2，IMF3，IMF4}，混合集合={IMF5，IMF6，IMF7}，有效集合={IMF8，IMF9}。

有效信息集合是低頻和次低頻的系統(tǒng)誤差。

對混合信息集合的混合信息集合元素采用本文新的閾值函數(shù)由式(7)～式(9)計算，分別小波濾噪，經(jīng)重構(gòu)后得到濾噪后的混合信息和混合信息中的隨機(jī)誤差，然后將混合信息的隨機(jī)誤差與高頻隨機(jī)誤差集合元素重構(gòu)得到了該測量數(shù)據(jù)的隨機(jī)誤差序列?；旌闲畔V噪后進(jìn)行重構(gòu)得到了混合信息中的系統(tǒng)誤差，并與其有效信息集合的元素重構(gòu)數(shù)據(jù)相加得到測量數(shù)據(jù)的系統(tǒng)誤差序列。

3.2 4種方法得到的結(jié)果分析

4種方法分離得到的系統(tǒng)誤差和隨機(jī)誤差如圖6所示。

從圖6(a)、圖6(b)和圖6(d)中的系統(tǒng)誤差和隨機(jī)誤差界限分明，系統(tǒng)誤差數(shù)值明顯小于隨機(jī)誤差，圖6(c)中系統(tǒng)誤差和隨機(jī)誤差混疊，系統(tǒng)誤差數(shù)值小于隨機(jī)誤差，圖6(d)中系統(tǒng)誤差曲線光順，系統(tǒng)誤差和隨機(jī)誤差沒有混疊。

采用4種方法分離系統(tǒng)誤差和隨機(jī)誤差后，用測量數(shù)據(jù)補(bǔ)償系統(tǒng)誤差，計算目標(biāo)彈道軌跡，得到與北斗定位結(jié)果的彈道差別均方根，如表3所示。

表3 彈道差別表 單位：m

從表3可以看出，本文方法與3個方法相比，在X、Y、Z軸3個方向精度都有所提高。X方向方法4、方法3、方法2的誤差分別是本文方法的1.59倍、1.28倍、1.09倍；Y方向方法4、方法3、方法2的誤差分別是本文方法的1.52倍、1.24、1.09倍；Z方向方法4、方法3、方法2的誤差分別是本文方法的1.52倍、1.30倍、1.06倍。顯然將本文方法分離得到的系統(tǒng)誤差補(bǔ)償?shù)綔y量數(shù)據(jù)，計算得到的彈道軌跡與真值差別最小、精度最高，特別是23.7 s數(shù)據(jù)，方法3和方法4分離的結(jié)果是估高了隨機(jī)誤差影響，降低了系統(tǒng)誤差估值，系統(tǒng)誤差補(bǔ)償不夠，導(dǎo)致計算的坐標(biāo)與真值坐標(biāo)差值很大，方法2分離的結(jié)果是估高了系統(tǒng)誤差數(shù)值，方法1分離的結(jié)果是系統(tǒng)誤差和隨機(jī)誤差共同作用引起殘差數(shù)據(jù)跳變，對測量數(shù)據(jù)補(bǔ)償?shù)南到y(tǒng)誤差值最合理，因此與真值彈道差別最小。

4 結(jié)束語

外測數(shù)據(jù)誤差分析是外彈道參數(shù)處理的關(guān)鍵環(huán)節(jié)，利用基于IMF能量拐點的測量數(shù)據(jù)成分集合判斷方法和基于反雙曲正弦函數(shù)調(diào)節(jié)的小波去噪方法，有效地分離了測量數(shù)據(jù)的系統(tǒng)誤差和隨機(jī)誤差。經(jīng)數(shù)據(jù)驗證，相比較本文新定義的閾值函數(shù)小波濾噪、EMD濾噪和多項式濾噪方法，本文方法與真值彈道差別最小，提高了定位精度。但是混合信息的小波濾噪如何設(shè)置最優(yōu)閾值需要進(jìn)一步解決。

圖6 4種方法分離得到的系統(tǒng)誤差和隨機(jī)誤差