基于本征模函數能量的外測數據誤差分析方法

劉一涵， 黃楷宸， 柏業超*

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

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

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

1 外測數據特點

外測測量設備包括各式雷達和光電經緯儀等。這些設備的測量元件可分為單脈沖雷達、相控陣雷達、紅外傳感器和可見光傳感器等，測量元素主要有距離、角度、速度，測量平臺有固定站、車載平臺、艦載平臺等，因而設備誤差也是復雜多樣的。其復雜性不僅表現在不同的測量體制間，也表現在同一測量體制、不同測量手段的測量數據，例如多傳感器光電經緯儀，雖然中波紅外傳感器和可見光傳感器都是光學測量，但是由于安裝的位置和成像的原理等不同，測量數據誤差特性也不相同。但是即使采用同一測量手段，不同時間進場的測量設備其數據誤差表現也不盡相同。外測數據隨機誤差具有不可觀測性、強相關性和非平穩時變特性。從工程技術來看，隨機誤差有時大于測量的量化誤差。從測量時間段看，各個時間區域方差大小是有變化的。系統誤差只有部分可以通過標校測定和修正，這樣的系統誤差數值較大，即使實際上可以建模的系統誤差采用最可靠的誤差修正模型，在彈道不同測量弧段也存在著不同的誤差修正殘余量，還有些系統誤差不能標校或者不能建立修正模型，這些系統誤差往往具有潛伏性，不易識別，很難與隨機誤差分離。

文獻[1]定義了外測測量數據的數學模型：

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

(1)

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

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

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

(2)

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

2 本文方法

2.1 外測數據3個集合{fIMF(i)}與{E(i)}的能量關系

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

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

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

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

(3)

外測數據序列的總能量記為E0，其計算公式為

(4)

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

(5)

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

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

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

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

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

劃分能量拐點的準則為

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

(6)

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

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

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

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

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

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

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

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

2.3 混合信息集合的處理

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

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

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

(7)

因此本文的閾值函數為

(8)

(9)

式中:σj為第j層小波系數的標準方差。

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

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

① 對x(t)進行EMD。

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

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

⑥ 有效信息與x1(t)重構后得到去噪后的有效外測數據。

3 實驗分析

3.1 分析案例

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

將北斗測量數據經坐標轉換得到待分析艦載雷達測量數據真值[15-16]，將真值和實際測量數據作差，得到殘差數據，此時殘差數據包括測量數據的隨機誤差和系統誤差。然后采用表1中4種方法分離殘差數據中系統誤差和隨機誤差。

表1 實驗方法

按照表1的4種方法，分離了隨機誤差和系統誤差，用濾掉隨機誤差的測量數據計算目標彈道軌跡，并用與北斗定位結果的彈道差異來分析方法的有效性。

采用本文方法結合方位角殘差數據說明具體分析過程。

圖1為方位角殘差數據圖。圖1中數據起伏較大，說明數據隨機性強，顯然是非平穩數據，若視為平穩數據濾波，必然會影響數據精度。

圖1 方位角殘差數據圖

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

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

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

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

圖2 IMF1～IMF9及其頻譜圖

圖3 IMF1～IMF9時域振幅

圖4 IMF1～IMF9頻域幅度圖

表2 IMF能量和能量變化

圖5 IMF能量和EMD分解層數圖

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

有效信息集合是低頻和次低頻的系統誤差。

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

3.2 4種方法得到的結果分析

4種方法分離得到的系統誤差和隨機誤差如圖6所示。

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

采用4種方法分離系統誤差和隨機誤差后，用測量數據補償系統誤差，計算目標彈道軌跡，得到與北斗定位結果的彈道差別均方根，如表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倍。顯然將本文方法分離得到的系統誤差補償到測量數據，計算得到的彈道軌跡與真值差別最小、精度最高，特別是23.7 s數據，方法3和方法4分離的結果是估高了隨機誤差影響，降低了系統誤差估值，系統誤差補償不夠，導致計算的坐標與真值坐標差值很大，方法2分離的結果是估高了系統誤差數值，方法1分離的結果是系統誤差和隨機誤差共同作用引起殘差數據跳變，對測量數據補償的系統誤差值最合理，因此與真值彈道差別最小。

4 結束語

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

圖6 4種方法分離得到的系統誤差和隨機誤差