裴 喆，史 震

（1.哈爾濱工程大學 自動化學院，哈爾濱 150001；2.92941 部隊，遼寧 葫蘆島 125001）

在靶場進行艦空導彈攔截靶標（靶機或靶彈）飛行試驗時，脫靶量是必須測量的參數(shù)，因為它是衡量導彈在彈道末端制導與控制性能的關(guān)鍵指標。脫靶量測量一般采用兩種方法：一種是設(shè)法測量導彈和靶標的外彈道數(shù)據(jù)，利用遭遇段的空間位置參數(shù)重建彈靶運動軌跡，外推計算脫靶量[1]。目前主要采用地面雷達或光電經(jīng)緯儀進行外彈道測量，雷達測量精度較低，光電經(jīng)緯儀測量精度高但易受氣象條件影響。另一種則是設(shè)法測量遭遇段彈靶之間的多普勒頻率，根據(jù)脫靶量與多普勒頻率之間的數(shù)學關(guān)系估計脫靶量[2-3]。目前主要采用遙測或靶機專用脫靶量測量設(shè)備進行多普勒頻率測量，靶機專用脫靶量測量設(shè)備精度較高，但重復測量成本高，且無法裝在靶彈上。遙測是靶場主要的內(nèi)彈道測量方法，技術(shù)成熟、工作可靠，遙測數(shù)據(jù)中包含引信探測的多普勒頻率。實際彈靶遭遇時，由于靶標體目標效應(yīng)的影響，多普勒頻率經(jīng)常包含許多孤立野值和成片野值，而如何去除成片野值在公開發(fā)表的文獻中還沒有見到通用有效的方法[4]。本文首先討論利用無野值的多普勒頻率估計脫靶量的方法，然后研究利用小波分解識別并去除野值的方法。

1 脫靶量數(shù)學模型

導彈和靶標遭遇過程中，由于相對速度大，遭遇時間短，可以認為相對速度保持不變。脫靶量數(shù)學模型如圖1所示，圖中rV為導彈相對靶標的運動速度，ρ為脫靶量（ρ垂直于 Vr），ti為多普勒頻率采樣時間，tρ為脫靶時間。

圖1 脫靶量數(shù)學模型

對于不同的ti(i=1,2,3,…,n)，多普勒頻率Ed(i)與脫靶量ρ的數(shù)學關(guān)系可由式(1)表示[2]：

式中：λ為引信天線波長；Vr可由引信開機前的導引頭遙測數(shù)據(jù)得出。由于引戰(zhàn)配合需要，引信在 tρ到來前引爆戰(zhàn)斗部，tρ時沒有Ed(i)數(shù)據(jù)，tρ和ρ只能由已測的Ed(i)數(shù)據(jù)根據(jù)式(1)計算得到。實際遙測測量的多普勒頻率fd(i)可由式(2)表示：

式中：ε (i)表示隨機誤差；εy(i)表示野值。

利用fd(i)估計ρ 需要盡量消除 ε (i)和 εy(i)，從而得到與Ed(i)非常接近的曲線。

2 多普勒頻率數(shù)據(jù)僅含隨機誤差時的脫靶量估計

先討論fd(i)僅含隨機誤差ε (i)的情況，即εy(i)=0。多數(shù)情況下 ε (i)近似于零均值高斯白噪聲，一般采用擬合或濾波的方法去除 ε (i)。小波變換可同時進行時域和頻域分析，具有可變的時頻分辨率，通過小波多尺度分解可將信號分解為反映整體趨勢的低頻分量和反映噪聲細節(jié)的高頻分量。

對fd(i)利用db4 小波進行4層分解，得到低頻系數(shù)a1、a2、a3、a4和高頻系數(shù)d1、d2、d3、d4，再對各層系數(shù)進行重構(gòu)。分析各層系數(shù)發(fā)現(xiàn)：從第1層到第4層，高頻系數(shù)將噪聲按頻率降低的順序分別提取出來；低頻系數(shù)所含噪聲越來越少，多普勒頻率趨勢越來越明顯。a1與fd(i)曲線如圖2所示，由于d1提取了白噪聲最尖銳的部分，a1不但沒有丟失Ed(i)的有用信息，而且比fd(i)光滑許多。對a1進行最小二乘擬合，使其與擬合曲線fd' (i)殘差序列的平方和 S (i)最小，表達式如式(3)所示：

最后由fd' (i)根據(jù)式(1)計算出ρ。

圖2 無野值多普勒頻率a1系數(shù)的擬和結(jié)果

為驗證方法的可行性和有效性，現(xiàn)利用蒙特卡洛方法進行仿真。給定 Vr、tρ、ρ、ti的一組設(shè)定值，取 Vr=1 200 m/s、tρ=30.000 s、ρ=10 m、ti=29.992 s時引信起爆。由式(1)得到Fd(i)，再模擬生成 ε (i)～N(0，1 0002)，按式(2)合成fd(i)序列。然后對fd(i)序列的小波系數(shù)a1進行最小二乘擬合，得到ρ和 tρ的估計值。由于Fd(i)都是ρ和 tρ的非線性函數(shù)，須采用非線性最小二乘擬合[5]。仿真1 000次的統(tǒng)計結(jié)果如表1所示，可以看出：

a）利用a1進行擬合精度很高，其中一次的擬合曲線fd' (i)如圖2所示，fd' (i)與Fd(i)幾乎重合。ρ的估計均值與設(shè)定值僅差約0.01 m，估計標準差僅約為0.24 m。統(tǒng)計還表明ρ的估計值近似服從正態(tài)分布，因此ρ的估計精度能夠控制在（?0.75 m～+0.75 m）范圍內(nèi)，完全滿足測量要求。

b）由于對實際數(shù)據(jù)進行擬合時事先并不知道真值，仿真中設(shè)定了3組不同的初值，結(jié)果表明擬合算法對不同初值有很好的收斂性。

因此，當多普勒頻率數(shù)據(jù)僅含隨機誤差時，利用該擬合算法估計脫靶量是可行而且有效的。

表1 a1系數(shù)擬合統(tǒng)計結(jié)果

3 多普勒頻率數(shù)據(jù)有野值時的脫靶量估計

3.1 野值識別

式(2)中的εy(i)表示野值，一批數(shù)據(jù)中有部分數(shù)據(jù)與其余數(shù)據(jù)相比明顯不一致的稱為野值，或稱離群值、異常值[6]，一般分為孤立型和成片型。由于彈靶遭遇時靶標的體目標效應(yīng)影響，實際的引信多普勒頻率數(shù)據(jù)經(jīng)常含有野值。直接對有野值fd(i)的a1系數(shù)進行擬合，結(jié)果的局部圖如圖3所示。擬合曲線fd' (i)與Fd(i)偏離較大，ρ估計值也與設(shè)定值有較大偏差，因此必須對野值進行識別、去除。

圖3 有野值多普勒頻率a1系數(shù)的擬和結(jié)果

通過隨機加入不同類型、不同數(shù)量、不同位置、不同幅值的野值 εy(i)對fd(i)進行模擬，利用db4小波進行4層分解，將重構(gòu)后的各層系數(shù)與無野值多普勒頻率小波分解后的各層系數(shù)對比分析發(fā)現(xiàn)：

a）孤立野值和白噪聲主要影響d1的幅值，因此a1能反映成片野值的輪廓，定位比較準確，a1局部圖如圖3所示。而且a1與擬合曲線fd' (i)的殘差序列能夠反映成片野值的大小。

b）成片野值主要影響d2、d3，它們局部幅值的增大說明成片野值的存在。如果設(shè)置一定門限可用于成片野值的識別，且這種識別方法與擬合曲線fd' (i)無關(guān)。

由此可總結(jié)出fd(i)有無野值的識別準則：

a）孤立野值和兩點成片野值可直接利用fd(i)與擬合曲線fd' (i)的殘差序列識別，門限值一般取

±3σ。[7]

b）成片野值可利用d2、d3的局部幅值識別，門限值可通過與無野值多普勒頻率的d2、d3幅值統(tǒng)計比較得到。

c）成片野值也可利用a1與擬合曲線fd' (i)的殘差序列識別，多次仿真發(fā)現(xiàn)無野值多普勒頻率的a1與擬和曲線fd' (i)的殘差序列中幾乎沒有連續(xù)三點數(shù)據(jù)幅值都超過1.5σ的。因此，識別兩點以上成片野值的門限值可取±1.5σ。

3.2 野值去除

3.2.1 野值去除順序

確定野值去除順序應(yīng)首先分析野值的特殊性，通過對不同類型、數(shù)量、位置、幅值的野值進行模擬發(fā)現(xiàn)：若幅值相同、分布位置相同，則成片野值相比孤立野值對擬合影響大；若幅值相同、類型相同，則位于后半?yún)^(qū)的野值比前半?yún)^(qū)的對擬合影響大，這是因為多普勒頻率是非線性的，離靶標較遠時變化較緩，較近時急劇下降。另外，由于靶標體目標效應(yīng)逐漸增強，通常后半?yún)^(qū)野值數(shù)量更多，幅值更大，且起伏劇烈。由于這些特性，擬合曲線fd' (i)對Fd(i)有時有較大偏離（尤其后半?yún)^(qū)），如圖3所示。

野值有正有負、有大有小，導致fd' (i)對Fd(i)有不同的偏離方向和不同的偏離程度，因此野值去除順序必須保證fd' (i)向Fd(i)逐步收斂。偏離的方向有全區(qū)上升、全區(qū)下降、前降后升、前升后降4種情況。全區(qū)上升時必須先去除正野值，后去除負野值；全區(qū)下降時與之相反。由野值的特殊性決定，前降后升時必須先去除后半?yún)^(qū)正野值和前半?yún)^(qū)負野值，最后去除其他野值；前升后降時與之相反。為表征fd' (i)偏離Fd(i)的方向和程度，將a1與fd' (i)的殘差序列定義為偏離因子，因為該序列能夠較準確地反映成片野值的位置和幅值。偏離因子初始值為0，它的正負對應(yīng)野值的正負，它的大小對應(yīng)偏離的程度。將偏離因子大小根據(jù)成片野值不同的長度、幅值、位置進行分檔，長度越長、幅值越大，則偏離因子越大，且后半?yún)^(qū)的偏離因子比前半?yún)^(qū)的大。再將殘差序列按時間等分為幾個（一般6個即可）分區(qū)，計算各分區(qū)偏離因子以及全區(qū)偏離因子總和，這樣就可以識別偏離方向及程度了。全區(qū)上升或下降表現(xiàn)為全區(qū)偏離因子總和較大，且各分區(qū)偏離因子大部分為正或大部分為負。前降后升或前升后降表現(xiàn)為全區(qū)偏離因子總和較小，且前半?yún)^(qū)的偏離因子與后半?yún)^(qū)的偏離因子正負相反。

由此可總結(jié)出野值去除的順序：先確定擬合曲線的偏離方向和程度；然后對比前后半?yún)^(qū)，若后半?yún)^(qū)野值多則先去除，少則全區(qū)同時去除；無論哪個區(qū)域都須先去除較長成片野值，再去除較短成片野值，最后去除孤立野值。

3.2.2 野值替換方法

野值替換是指在野值點fd(j)處的擬合曲線fd' (j)數(shù)據(jù)增加適當?shù)恼{(diào)整值（一般為k?σ），生成新的fd(j)數(shù)據(jù)：

孤立野值替換直接使用式(4)。成片野值替換需調(diào)整為式(5)：

因為如果成片野值的替換數(shù)據(jù)相同會影響再次擬和的結(jié)果，所以應(yīng)在式(4)基礎(chǔ)上增加與成片野值相同長度的零均值高斯白噪聲ε (j,j+m)～N(0,σ2)，m+1為成片野值長度。

k的正負與野值的正負相反，野值為正時，調(diào)整值為負；野值為負時，調(diào)整值為正。k的大小取決于擬和曲線fd' (i)偏離Fd(i)的程度，可用仿真統(tǒng)計的方法確定k的取值范圍。

3.3 算法仿真

為驗證野值識別和去除方法的有效性，并確定k的取值范圍，編制了仿真程序，程序流程分為：

第1步：模擬生成如式(2)所示的fd(i)數(shù)據(jù)。其中野值 εy(i)由人工隨機設(shè)置，成片野值長度最長為7點，幅值最大5倍σ，數(shù)量不大于全部數(shù)據(jù)的二分之一。

第2步：利用db4 小波對fd(i)進行4層分解，并重構(gòu)各層系數(shù)。對a1進行非線性最小二乘擬合，得到擬合曲線fd' (i)、擬合殘差序列及ρ的估計值。

第3步：按3.1節(jié)中的準則對fd(i)進行野值識別，若有野值，執(zhí)行第4步；若無野值，則程序結(jié)束，ρ即為所求脫靶量。

第4步：按3.2節(jié)中的方法去除部分野值，對野值替換后生成的新fd(i)數(shù)據(jù)循環(huán)執(zhí)行第2步。

通過設(shè)置不同野值進行仿真后結(jié)果表明，k的取值范圍主要與野值分布有關(guān)。后半?yún)^(qū)野值較少時，全區(qū)替換時k 一般在0＜k≤1范圍內(nèi)取值；后半?yún)^(qū)野值較多時，后半?yún)^(qū)k的取值范圍能達到 0＜k＜2，將后半?yún)^(qū)大部分野值替換后，再在0＜k≤1范圍內(nèi)取值對全區(qū)替換即可。

為提高替換精度，可將 k?σ的大小按0.25·σ分檔，即0.25·σ、0.5·σ、0.75·σ、σ，具體選檔需根據(jù)擬合曲線偏離程度確定。

仿真還表明，如果野值數(shù)量小于全部數(shù)據(jù)的三分之一，后半?yún)^(qū)野值較少時，野值去除后的脫靶量估計精度能控制在（?1 m～+1 m）范圍內(nèi)；后半?yún)^(qū)野值較多時，也能控制在（?1.5 m～+1.5 m）范圍內(nèi)，仍能滿足測量要求，表明算法是合理有效的。但如果野值數(shù)量大于全部數(shù)據(jù)的三分之一，擬合曲線有時會發(fā)散。

4 結(jié)束語

本文研究了一種利用小波分解和非線性最小二乘擬和相結(jié)合估計艦空導彈脫靶量的方法。通過仿真驗證可知，對于同時含有成片野值和孤立野值的引信多普勒頻率數(shù)據(jù)，如果野值數(shù)量小于全部數(shù)據(jù)的三分之一，則利用該方法估計的脫靶量誤差在（?1.5 m～+1.5 m）范圍內(nèi)，滿足測量要求，但野值較多時的去除方法仍需加以研究。

[1]袁俊.由雷達測量數(shù)據(jù)求解脫靶量算法研究[D].南京:南京理工大學,2003.

[2]林松.用導彈“視頻遙測”信號推算導彈脫靶量[J].上海航天,2002(3):35-38.

[3]魏國華,吳嗣亮,王菊,等.脫靶量測量技術(shù)綜述[J].系統(tǒng)工程與電子技術(shù),2004,26(6):768-772.

[4]胡紹林,孫國基.靶場外測數(shù)據(jù)野值點的統(tǒng)計診斷技術(shù)[J].宇航學報,1999,20(2):68-74.

[5]張光澄.非線性最優(yōu)化計算方法[M].北京:高等教育出版社,2005:159-168.

[6]張德然.統(tǒng)計數(shù)據(jù)中異常值的檢驗方法[J].統(tǒng)計研究,2003(5):53-55.

[7]費業(yè)泰.誤差理論與數(shù)據(jù)處理[M].5 版.北京:機械工業(yè)出版社,2008:43-49.