侯博文，王炯琦,2，周萱影，李 冬，何章鳴,2

(1.國防科技大學(xué)理學(xué)院，長沙 410072； 2.北京控制工程研究所，北京 100086；3.中國人民解放軍91550部隊94分隊，大連 116023)

0 引言

高精度彈道跟蹤數(shù)據(jù)是分析和鑒定航天飛行器的制導(dǎo)系統(tǒng)、再入系統(tǒng)、測控系統(tǒng)及其它分系統(tǒng)的重要基礎(chǔ)。然而，無論采用何種彈道跟蹤體制對飛行器進(jìn)行跟蹤測量，跟蹤或測量設(shè)備的沖擊、振動、系統(tǒng)故障、環(huán)境干擾，或操作人員的誤判等都會造成過失誤差，引起異常值。測量數(shù)據(jù)集合中偏離正常數(shù)據(jù)所呈現(xiàn)趨勢的數(shù)據(jù)點(diǎn)稱為野值[1]。野值的存在會嚴(yán)重影響彈道跟蹤數(shù)據(jù)的精度，甚至歪曲測量體系真相，在需要高精度彈道跟蹤數(shù)據(jù)的場合是不允許的。因而，必須在彈道跟蹤數(shù)據(jù)預(yù)處理中，檢測、識別野值，并剔除或作必要的修正。

根據(jù)野值是否連續(xù)，可分為孤立型和斑點(diǎn)型野值[2]。根據(jù)彈道跟蹤數(shù)據(jù)處理模式又可分為事后處理和實時處理。不同情況下跟蹤測量數(shù)據(jù)的野值特征、表現(xiàn)形式、變化規(guī)律各不相同。目前已經(jīng)有很多學(xué)者提出了不同的野值剔除方法[3-10]。傳統(tǒng)的野值剔除方法[3](包括萊特準(zhǔn)則、羅曼諾夫斯基準(zhǔn)則、格拉布斯準(zhǔn)則及肖維勒準(zhǔn)則等)對于彈道數(shù)據(jù)的統(tǒng)計特性具有約束性，導(dǎo)致這些方法在使用時受限；自適應(yīng)最小二乘法[4]對于孤立型野值具有較好的剔除效果，但沒有討論關(guān)于斑點(diǎn)型野值處理效果；中值濾波差分法[7]可有效剔除斑點(diǎn)型野值，但剔除效果不穩(wěn)定，會出現(xiàn)數(shù)據(jù)失真或缺失的情況；外推擬合法[8]可有效剔除孤立型和斑點(diǎn)型野值，但斑點(diǎn)型野值剔除后會出現(xiàn)數(shù)據(jù)缺失；插值剔除法[9]適用于光學(xué)測量設(shè)備，在其他情況下效果并不顯著。

本文對上述文獻(xiàn)中提到的各種野值剔除方法進(jìn)行了綜述，從理論上分析了各種方法的優(yōu)缺點(diǎn)，基于各野值剔除方法建立了相應(yīng)的識別準(zhǔn)則和剔除步驟，并通過仿真分析了其針對孤立型野值和斑點(diǎn)型野值的剔除能力，對比了各類野值剔除方法效果，為實際工程應(yīng)用提供了理論與技術(shù)支撐。

1 事后數(shù)據(jù)野值剔除方法

彈道跟蹤數(shù)據(jù)的事后野值剔除是指在導(dǎo)彈或火箭飛行試驗后對彈道測量數(shù)據(jù)的野值剔除。

假設(shè)測量數(shù)據(jù)序列為y=[y1y2…yn]Τ，其可以表示為

(1)

由于飛行器的飛行數(shù)據(jù)一般是具有趨勢的數(shù)據(jù)，因此事后野值剔除主要分為兩步，即趨勢剔除(包括真實信號特征提取、系統(tǒng)誤差補(bǔ)償?shù)?和野值剔除。這里主要介紹野值剔除的方法。假設(shè)趨勢剔除后的數(shù)據(jù)為x=[x1x2…xn]Τ，趨勢剔除后的測量數(shù)據(jù)殘差序列為v=[v1v2…vn]Τ。工程中常用的有以下幾種事后野值處理準(zhǔn)則。

1.1 萊特準(zhǔn)則

萊特準(zhǔn)則是最常用、最簡單的判別準(zhǔn)則，該方法以測量次數(shù)充分大為前提[11]。

原理：對于某一個測量序列，假設(shè)各測量值只含有隨機(jī)誤差，則根據(jù)隨機(jī)誤差的正態(tài)分布規(guī)律，其殘差落在3σ以外的概率不超過0.3%。

準(zhǔn)則：以3σ為野值判斷閾值，其中σ為測量殘差v的標(biāo)準(zhǔn)差。

步驟如下：

1) 計算剔除趨勢后數(shù)據(jù)均值及標(biāo)準(zhǔn)差為

(2)

2) 計算各個時刻數(shù)據(jù)的殘差，有

(3)

3) 將殘差逐一與3σ進(jìn)行比較，若

vi≥3σ

(4)

1.2 羅曼諾夫斯基準(zhǔn)則

羅曼諾夫斯基準(zhǔn)則又稱為t檢驗準(zhǔn)則，對于一批獨(dú)立等精度測量數(shù)據(jù)中存在單個異常值的情況，是一種有效的識別方法[12]。

原理：根據(jù)抽樣分布定理，按照t分布的實際誤差分布范圍來檢測野值。

準(zhǔn)則：根據(jù)樣本確定的檢測統(tǒng)計量與給定顯著水平下確定的t檢驗系數(shù)進(jìn)行比較，確定野值。

步驟如下：

1) 根據(jù)抽樣分布定理，確定服從t分布的隨機(jī)變量，有

(5)

1.3 格拉布斯準(zhǔn)則

格拉布斯準(zhǔn)則又稱極大殘差檢驗[14]，是檢測服從正態(tài)分布的單變量數(shù)據(jù)集中野值的方法。

原理：將測量值由大到小排序，依次進(jìn)行野值剔除。

準(zhǔn)則：確定構(gòu)造統(tǒng)計量的分布，根據(jù)預(yù)設(shè)的顯著性水平剔除野值。

步驟如下：

1) 將xi按大小順序排列成順序統(tǒng)計量x(i)，滿足

x(1)≤x(2)≤…≤x(n)

(6)

2) 構(gòu)造統(tǒng)計量

(7)

3) 取定顯著性水平α，有

P(gn≥G(n,α))=α

(8)

其中格拉布斯確定了g分布的臨界值[11]，有

(9)

1.4 肖維勒準(zhǔn)則

肖維勒準(zhǔn)則以正態(tài)分布為前提，是一種剔除實驗數(shù)據(jù)中野值的有效方法[15]。

原理：設(shè)在一個n次的等精度測量中，不出現(xiàn)vi>a的誤差，那么概率P{vi>a}接近于0。當(dāng)n足夠大時，由大數(shù)定律，概率P≈m/n，其中m是vi>a出現(xiàn)的次數(shù)。因此可認(rèn)為

m/n?P{vi>a}→0

(10)

或

m?nP{vi>a}→0

(11)

式中：m為整數(shù)。因此有

nP{vi>a}≤1/2

(12)

即最低限度應(yīng)有

P{vi>a}=1-erf(k)=1/(2n)

(13)

式中：k=a/σ，k值可查表[16]獲得。

準(zhǔn)則：確定一個概率范圍，其中心為一個標(biāo)準(zhǔn)正態(tài)分布的均值，在該范圍之外的點(diǎn)即可判斷為野值。

步驟如下：

1)計算數(shù)據(jù)殘差vi。

2)剔除野值，如滿足

vi>kσ

(14)

1.5 奇偶提取法

奇偶提取法，即將數(shù)據(jù)按照奇偶位置分為若干組(集)進(jìn)行野值剔除[17]。

原理：原數(shù)據(jù)樣本標(biāo)準(zhǔn)差誤差較大，將原數(shù)據(jù)拆分成小樣本數(shù)據(jù)集后，每個小樣本數(shù)據(jù)集中野值個數(shù)減少，標(biāo)準(zhǔn)差更接近于真值。

準(zhǔn)則：分組后的數(shù)據(jù)按照萊特準(zhǔn)則進(jìn)行野值剔除。

步驟如下：

1) 奇偶序列提取。對于數(shù)據(jù)x1,x2,…,xn按照奇偶位置分成兩組，分別記為x1-O,x1-E；再進(jìn)行第二次分組，得到x2-OO,x2-OE,x2-EO,x2-EE；以此類推，將數(shù)據(jù)x1,x2,…,xn分成2m組。

(15)

(16)

(17)

(18)

以此類推，根據(jù)需要選定分組次數(shù)m(根據(jù)精度設(shè)定)，分成2m個小樣本數(shù)據(jù)集，對每組都進(jìn)行上述求均值和方差的過程。

3) 利用萊特判別準(zhǔn)則進(jìn)行野值檢驗和剔除。

4) 數(shù)組整合。

經(jīng)過步驟1)～3)后，將剔除野值后的分組數(shù)據(jù)進(jìn)行整合，得到完整的數(shù)據(jù)。

需要注意的是，在每組小樣本數(shù)據(jù)集里可能存在第1個數(shù)據(jù)就是野值的情況，此時用第2個數(shù)據(jù)替代。

1.6 差分輔助中值濾波

在中值濾波的基礎(chǔ)上借差分思想對數(shù)據(jù)進(jìn)行野值剔除。

原理：對測量數(shù)據(jù)進(jìn)行一次差分，可得

(ωi+1-ωi)+(δi+1-δi)

(19)

考慮野值差分結(jié)果為(δi+1-δi)，對于斑點(diǎn)型野值，若野值點(diǎn)幅度相近，在差分?jǐn)?shù)據(jù)中，野值點(diǎn)數(shù)據(jù)將只保留在當(dāng)前窗口內(nèi)的第1個野值數(shù)據(jù)，其余野值將會被差分消除。因此，斑點(diǎn)型野值會變?yōu)楣铝⑿鸵爸?。孤立型野值可直接利用中值濾波方法進(jìn)行剔除。

準(zhǔn)則：對一次中值濾波后的觀測數(shù)據(jù)進(jìn)行一次差分，對差分后的數(shù)據(jù)再進(jìn)行中值濾波，積分還原后進(jìn)行野值判斷。

步驟如下：

1) 對原始測量數(shù)據(jù)x進(jìn)行中值濾波得到xmed，剔除外彈道測量數(shù)據(jù)中的孤立型野值。

2) 對xmed進(jìn)行一次差分后，再進(jìn)行一次中值濾波剔除野值。

3) 對差分后的數(shù)據(jù)進(jìn)行積分，得到原測量數(shù)據(jù)。

4) 檢測剔除后的數(shù)據(jù)是否還存在野值，如有，重復(fù)步驟1)～3)。

2 實時數(shù)據(jù)野值剔除方法

彈道跟蹤數(shù)據(jù)的實時野值剔除是指在導(dǎo)彈或火箭飛行試驗過程中，對彈道測量數(shù)據(jù)的野值進(jìn)行實時檢測與剔除[4]。

由于實時數(shù)據(jù)的獲取會有延遲，所以實時數(shù)據(jù)野值剔除也主要分為兩步，即數(shù)據(jù)外推和野值剔除，此外，還要求起始段數(shù)據(jù)無野值情況。

2.1 五點(diǎn)線性預(yù)報野值剔除法

五點(diǎn)線性預(yù)報法分為差分檢驗、線性預(yù)報。首先用求一階差分和四階差分的方法對數(shù)據(jù)進(jìn)行合理性檢驗，發(fā)現(xiàn)并剔除明顯的野值點(diǎn)，并按五點(diǎn)線性預(yù)報公式補(bǔ)點(diǎn)。

原理：利用最小二乘法對數(shù)據(jù)進(jìn)行預(yù)報，與實測數(shù)據(jù)做對比，從而實現(xiàn)野值剔除。

準(zhǔn)則：將測量數(shù)據(jù)殘差與預(yù)先確定的野值檢測門限做對比并進(jìn)行檢測。

步驟如下：

1) 對跟蹤數(shù)據(jù)進(jìn)行一階差分，其表達(dá)式為

Δ1yi=yi+1-yi

(20)

用四階差分檢驗法進(jìn)行初始檢驗，找出一組合理點(diǎn)。

2) 數(shù)據(jù)四階差分值

Δ4yj=yj-4-4yj-3+6yj-2-4yj-1+yj

(21)

式中：j≥5。取門限值M1(經(jīng)驗值為17σ，σ為測量數(shù)據(jù)的精度)，判斷Δ4yj≤M1是否成立，若是，則為一組合理點(diǎn)；否則，令j=j+1，繼續(xù)進(jìn)行四階差分檢驗。

(22)

2.2 α-β-γ濾波器野值剔除法

原理：利用α-β-γ濾波器對數(shù)據(jù)進(jìn)行預(yù)報，與實測數(shù)據(jù)做對比，從而實現(xiàn)野值剔除。

準(zhǔn)則：將測量數(shù)據(jù)殘差與預(yù)先確定的野值檢測門限做對比并進(jìn)行檢測。

步驟如下：

(23)

式中：u，s為濾波器的中間變量。

2) 從第4點(diǎn)開始按照式(24)、(25)遞推計算，有

(24)

(25)

3) 對測量數(shù)據(jù)序列進(jìn)行一步預(yù)測，若實測值與預(yù)測值之差的絕對值小于某門限值，此數(shù)據(jù)為合理值，反之為野值。

(26)

式中：M為門限值，一般為正常訓(xùn)練數(shù)據(jù)預(yù)測殘差精度的3倍。

2.3 基于穩(wěn)健估計的自適應(yīng)門限法

由于數(shù)據(jù)中存在的野值會影響樣本標(biāo)準(zhǔn)差的確定，結(jié)合穩(wěn)健估計原理，實時確定野值檢測門限，從而實現(xiàn)跟蹤數(shù)據(jù)的野值實時剔除[18]。

(27)

式中：ψ(·)為影響函數(shù)；β為待定參數(shù)。適當(dāng)選擇ψ(y)和β可對數(shù)據(jù)中異常值加以抑制。

2)β公式的推導(dǎo)。選擇Huber的ψH(y0)函數(shù)作為影響函數(shù)，有

(28)

(29)

因為異常值不多且可用CΗσ代替，故有

(30)

(31)

因此取不同的CH有不同的β，在實時數(shù)據(jù)處理過程中一般取CH=1.5，此時β=0.778 5。又有

(32)

于是有

(33)

(34)

式中：∑1為滿足Δyi≤CHσ的觀測數(shù)據(jù)求和；NH為觀測數(shù)據(jù)中Δyi>CHσ的數(shù)據(jù)個數(shù)；n為窗口大小，通常取16～30。由于實測數(shù)據(jù)是動態(tài)的，故采用滑動數(shù)據(jù)窗口。

準(zhǔn)則：將測量數(shù)據(jù)殘差與實時更新的野值檢測門限做對比，實現(xiàn)野值檢測和剔除。

步驟：

2) 判斷是否滿足Δyi≤CHσ，若滿足，則yj+1為合理值；否則，為野值，用預(yù)測值代替。

2.4 基于匹配測元的方法

彈道跟蹤數(shù)據(jù)中存在匹配關(guān)系的測元，如某一個測元是另一個測元的微分(或積分)，根據(jù)某一測元的正常數(shù)據(jù)來剔除另一測元測量數(shù)據(jù)的野值。

(35)

式(35)為匹配測元之間的關(guān)系，可改寫為

(36)

在實際測量數(shù)據(jù)中，由于系統(tǒng)誤差和隨機(jī)誤差，有

(37)

(38)

準(zhǔn)則：測量和匹配數(shù)據(jù)的殘差與設(shè)定的野值檢測門限做對比，實現(xiàn)野值檢測和剔除。

步驟如下：

1) 選取連續(xù)n個點(diǎn)的匹配測元數(shù)據(jù)，計算

(39)

并進(jìn)行統(tǒng)計，有

(40)

若ΔL(t1)>3σ，則R(t1)為野值點(diǎn)，其中σ為數(shù)據(jù)精度。記野值點(diǎn)個數(shù)為l，則有

(41)

2) 以n點(diǎn)為窗口滑動，按照步驟1)進(jìn)行數(shù)據(jù)檢測。

3 仿真驗證

針對實際彈道跟蹤測量存在的孤立型野值和斑點(diǎn)型野值，分別就事后處理和實時處理2種模式，進(jìn)行野值剔除的仿真計算。為了更好地描述每種方法的性能，定義野值剔除率為

(42)

3.1 事后數(shù)據(jù)野值剔除仿真

用計算機(jī)仿真出200個服從均值為0，方差為2的剔除趨勢后的數(shù)據(jù)點(diǎn)，在第96～105個點(diǎn)間加入常值為10的野值(斑點(diǎn)型野值)，在第50、70、130、150個點(diǎn)分別加入常值為-15、-10、15、12的野值(孤立型野值)。為更好地體現(xiàn)各方法的性能，在第151～160個點(diǎn)間加入了10sin(i)的斑點(diǎn)型時變野值，i為時間。分別用第2節(jié)中的6種方法進(jìn)行野值剔除，顯著性水平設(shè)置為0.01，進(jìn)行了10次仿真，取某一次仿真結(jié)果如圖1所示。圖1(a)～(d)中紅線為野值檢測閾值，藍(lán)線為數(shù)據(jù)殘差；圖1(e)～(f)中紅線為剔除野值后數(shù)據(jù)，藍(lán)線為剔除野值前數(shù)據(jù)。

10次仿真的數(shù)據(jù)野值剔除率和誤檢點(diǎn)數(shù)統(tǒng)計結(jié)果見表1。

表 1 事后數(shù)據(jù)野值檢測結(jié)果

差分輔助中值濾波法是直接將野值剔除，故無法判斷野值檢測結(jié)果。從圖1(f)可看出，該方法在一定程度上可實現(xiàn)野值剔除，但該方法受數(shù)據(jù)本身波動影響較大，當(dāng)原始數(shù)據(jù)出現(xiàn)較大波動時，野值剔除后的數(shù)據(jù)也會出現(xiàn)波動，效果不穩(wěn)定。

綜合表1中準(zhǔn)則，結(jié)合仿真結(jié)果，得到各種方法對比結(jié)果見表2。

3.2 實時數(shù)據(jù)野值剔除仿真

設(shè)定二次多項式仿真彈道數(shù)據(jù)

y(t)=0.4×t2+5×t+220

(43)

圖1 事后數(shù)據(jù)野值剔除結(jié)果Fig.1 Results of off-line data outlier elimination

方法統(tǒng)計量門限性能適用范圍萊特準(zhǔn)則xi-x3σ 簡單易行,大數(shù)據(jù)量情況下對孤立型野值剔除效果好 數(shù)據(jù)量較大,無時變型野值,需要快速剔除野值羅曼諾夫斯基準(zhǔn)則xi-xK(n,α)σi 野值剔除效果較好,但需要逐點(diǎn)計算均值、方差,計算復(fù)雜度高 數(shù)據(jù)量在4～20的較少數(shù)據(jù),無時變型野值格拉布斯準(zhǔn)則x(n)-xσG(n,α)σ 野值剔除效果最好 數(shù)據(jù)量在20～100的較少數(shù)據(jù),存在時變型野值肖維勒準(zhǔn)則xi-xkσ 野值剔除效果較好,計算復(fù)雜度相對較低 數(shù)據(jù)量大,無時變型野值,斑點(diǎn)型野值較少,無時變型野值奇偶提取法xi-x3σn-O/n-E 能有效地剔除斑點(diǎn)型野值,存在虛警情況 數(shù)據(jù)量大,無時變型野值,孤立型野值較少,存在時變型野值差分輔助中值濾波—— 具有一定剔除野值效果,且可處理斑點(diǎn)型野值,但處理效果不穩(wěn)定 數(shù)據(jù)波動小,存在時變型野值

式中：y(t)為位置數(shù)據(jù)；t為時間。共設(shè)200個采樣點(diǎn)，采樣時間間隔為0.05 s，野值與事后處理仿真設(shè)定相同，共進(jìn)行10次仿真，取某一次仿真結(jié)果如圖2、3所示。圖2中紅線為檢測門限，藍(lán)線為預(yù)測殘差；圖3中紅線為剔除野值后數(shù)據(jù)，藍(lán)線為剔除野值前數(shù)據(jù)。

圖2 實時數(shù)據(jù)野值檢測結(jié)果Fig.2 Results of on-line data outlier detection

圖3 實時數(shù)據(jù)野值剔除結(jié)果Fig.3 Results of on-line data outlier elimination

由于實時野值剔除方法對計算效率有較高的要求，因此幾種方法的時間復(fù)雜度和仿真運(yùn)算見表3。

綜合以上幾種實時數(shù)據(jù)處理方法，并結(jié)合圖2、3及表3可得各種方法對比結(jié)果，見表4。

表3 實時數(shù)據(jù)野值檢測結(jié)果

表4 實時數(shù)據(jù)野值剔除方法效果對比

4 結(jié)束語

野值剔除結(jié)果會直接影響彈道跟蹤數(shù)據(jù)處理的精度。本文結(jié)合彈道數(shù)據(jù)事后處理和實時處理2種模式，對常見的，用于彈道跟蹤數(shù)據(jù)的野值剔除方法進(jìn)行了綜述。通過理論分析和仿真計算表明：不同的野值剔除方法對野值的剔除效果不同。在事后數(shù)據(jù)野值剔除中，羅曼諾夫斯基準(zhǔn)則對斑點(diǎn)型野值和孤立型野值的剔除效果較好，但需要以犧牲計算效率為代價；肖維勒準(zhǔn)則對于斑點(diǎn)型野值剔除效果較好，且計算效率高；萊特準(zhǔn)則在數(shù)據(jù)量較大情況下對斑點(diǎn)型野值剔除效果較好；格拉布斯準(zhǔn)則對于孤立型野值和斑點(diǎn)型野值剔除效果最好；奇偶提取法對于斑點(diǎn)型野值和時變型野值剔除效果較好；差分輔助中值濾波法在一定程度上具有剔除野值的作用，但效果不穩(wěn)定。在實時數(shù)據(jù)野值剔除中，五點(diǎn)線性預(yù)報法比α-β-γ濾波法好，自適應(yīng)門限方法比固定門限方法好，主要體現(xiàn)在虛警點(diǎn)少，剔除效率高。匹配測元方法的野值剔除效果僅次于基于自適應(yīng)門限的五點(diǎn)線性預(yù)報法，且虛警點(diǎn)最少，但時間復(fù)雜度較高，實時性相對較差。當(dāng)孤立型、斑點(diǎn)型、時變型野值都存在時，基于自適應(yīng)門限的五點(diǎn)線性預(yù)報法能較好地實現(xiàn)野值剔除功能。

綜合比較可知：在事后處理中，對孤立型野值，羅曼諾夫斯基準(zhǔn)則和格拉布斯準(zhǔn)則可有效剔除；對斑點(diǎn)型野值，格拉布斯準(zhǔn)則的剔除效果最好；對時變型野值，格拉布斯準(zhǔn)則和奇偶提取法都能實現(xiàn)一定程度的剔除?？紤]到各種方法的適用性，當(dāng)彈道數(shù)據(jù)量較大時，對彈道數(shù)據(jù)的孤立型野值剔除也可選擇萊特準(zhǔn)則；對斑點(diǎn)型野值和時變型野值可選擇奇偶提取法；在實時處理中，各種方法都能有效剔除孤立型野值，出于實時性要求，選擇五點(diǎn)線性預(yù)報法最佳；基于自適應(yīng)門限的五點(diǎn)線性預(yù)報法對斑點(diǎn)型野值的剔除效果最好，效率最高；α-β-γ濾波法對時變型野值的剔除效果最好；如果3種野值同時存在時，采用基于自適應(yīng)門限的五點(diǎn)線性預(yù)報法最合適。