一種估計艦空導彈脫靶量的方法

裴 喆，史 震

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

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

1 脫靶量數學模型

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

圖1 脫靶量數學模型

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

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

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

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

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

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

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

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

圖2 無野值多普勒頻率a1系數的擬和結果

為驗證方法的可行性和有效性，現利用蒙特卡洛方法進行仿真。給定 Vr、tρ、ρ、ti的一組設定值，取 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)序列的小波系數a1進行最小二乘擬合，得到ρ和 tρ的估計值。由于Fd(i)都是ρ和 tρ的非線性函數，須采用非線性最小二乘擬合[5]。仿真1 000次的統計結果如表1所示，可以看出：

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

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

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

表1 a1系數擬合統計結果

3 多普勒頻率數據有野值時的脫靶量估計

3.1 野值識別

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

圖3 有野值多普勒頻率a1系數的擬和結果

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

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

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

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

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

±3σ。[7]

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

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

3.2 野值去除

3.2.1 野值去除順序

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

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

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

3.2.2 野值替換方法

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

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

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

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

3.3 算法仿真

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

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

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

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

第4步：按3.2節中的方法去除部分野值，對野值替換后生成的新fd(i)數據循環執行第2步。

通過設置不同野值進行仿真后結果表明，k的取值范圍主要與野值分布有關。后半區野值較少時，全區替換時k 一般在0＜k≤1范圍內取值；后半區野值較多時，后半區k的取值范圍能達到 0＜k＜2，將后半區大部分野值替換后，再在0＜k≤1范圍內取值對全區替換即可。

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

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

4 結束語

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

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

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

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

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

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

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

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