引信多普勒頻率的線性穩(wěn)健擬合

裴 喆

（92941部隊，遼寧葫蘆島125000）

防空導(dǎo)彈前向探測無線電引信在彈目接近時實時測量多普勒頻率。多普勒頻率不僅是導(dǎo)彈飛行試驗分析引信啟動和引戰(zhàn)配合結(jié)果的重要參數(shù)，還可作為光測輔助手段估算導(dǎo)彈脫靶量[1-4]。由于體目標(biāo)效應(yīng)影響多普勒頻率經(jīng)常包含誤差（尤其野值）[2-7]，給結(jié)果分析帶來困難。文獻[1]研究了基于非線性最小二乘的多普勒頻率擬合以及利用擬合殘差去除野值的方法；文獻[2]研究了非線性最小二乘擬合和小波分解相結(jié)合的方法，可去除野值比例達1/3，但野值識別和替換步驟有些復(fù)雜。為此，本文提出一種多普勒頻率的線性穩(wěn)健擬合方法，該方法的抗野值能力更強，步驟卻更簡單。

1 引信多普勒頻率線性轉(zhuǎn)換

1.1 引信多普勒頻率數(shù)學(xué)模型

防空導(dǎo)彈與目標(biāo)遭遇時間短，可認(rèn)為彈目相對速度Vr保持不變，彈目遭遇過程如圖1所示。

圖1 彈目遭遇過程Fig.1 Missile-target encounter process

圖1 中：ρ 為脫靶量（ρ 垂直于Vr）；tρ為脫靶時間。則多普勒頻率理論值為[1-3]：

式（1）中：λ 為引信工作波長；Vr通常可由導(dǎo)引頭測量得到。

由于引信近程探測的體目標(biāo)效應(yīng)影響，引信多普勒頻率實測值fd(t)的中后段會出現(xiàn)起伏，這種起伏可用加法性誤差ε(t)表示，即：

取Vr=1 000 m/s 、tρ=20.000 0 s 、ρ=10 m ，在Fd(t)中后段加入白噪聲誤差ε(t)～N(0,σ20)合成模擬的實測值fd(t)。σ0=800 Hz，F(xiàn)d(t)和fd(t)如圖2所示。

圖2 多普勒頻率理論值Fd(t)和實測值fd(t)曲線Fig.2 Theoretical Doppler Fd(t) and actual Doppler fd(t)

多普勒頻率擬合的目的就是濾除fd(t)中的ε(t)，使擬合曲線盡量逼近Fd(t)。因?qū)椧龖?zhàn)配合延時需要，引信須在tρ到來前引爆戰(zhàn)斗部。所以，ρ 和tρ為未知量，但tρ和ρ 可由fd(t)擬合估計得到。一次彈目遭遇時ρ 和tρ是唯一的，因而ρ 和tρ的估計精度可用來衡量fd(t)的擬合精度。

1.2 引信多普勒頻率線性轉(zhuǎn)換

把式（3）等號左側(cè)函數(shù)記為Fx(t)，即：

再令K1=-1/ρ，K2=tρ/ρ，則得Fx(t)關(guān)于時間t 的一次線性函數(shù)為：

式（5）中：K1和K2為多項式系數(shù)，因而式（4）即為多普勒頻率線性轉(zhuǎn)換表達式。

同理，將fd(t)線性轉(zhuǎn)換為fx(t)，則fx(t)只與已知量fd(t)、λ 和Vr有關(guān)，而系數(shù)K1、K2只與未知量ρ、tρ有關(guān)。將圖2 中的Fd(t) 和fd(t) 分別轉(zhuǎn)換為Fx(t) 和fx(t)，結(jié)果如圖3所示。

圖3 多普勒頻率線性轉(zhuǎn)換的Fx(t)和fx(t)曲線Fig.3 Fx(t) and fx(t) of Doppler linear transform

對fx(t)線性穩(wěn)健擬合估計得到K1和K2，計算出ρ 和tρ，再利用式（1）即可計算fd(t)的擬合結(jié)果。

2 線性穩(wěn)健擬合方法

2.1 線性穩(wěn)健擬合及其權(quán)重函數(shù)

普通最小二乘擬合給予每個測量數(shù)據(jù)的權(quán)重均為1，而擬合目標(biāo)是使殘差平方和最小，因而擬合結(jié)果受野值（也稱異常值、粗差）的影響較大。穩(wěn)健擬合就是為抵抗野值發(fā)展起來的。目前，線性穩(wěn)健擬合理論較為成熟，主要有包括L估計、R估計和M估計等多種方法[8-9]。本文采用M估計法，它將經(jīng)典極大似然估計的思想推廣用于穩(wěn)健擬合[8，10]，其基本思想是采用迭代加權(quán)最小二乘法估計多項式系數(shù)，根據(jù)權(quán)重函數(shù)確定殘差的權(quán)重大小，殘差越大，權(quán)重越小，從而達到抵抗野值的目的[11-12]。

選擇合適的權(quán)重函數(shù)是M估計法的重要環(huán)節(jié)，比較典型的函數(shù)有Huber、Bisquare、Andrew 等。根據(jù)fx(t)的誤差特性，本文采用Bisquare權(quán)重函數(shù)（也稱雙平方權(quán)重函數(shù)），即[8，13]：

式（6）中：y 為歸一化殘差；cB為調(diào)節(jié)常數(shù)（通常取4.685）[8]，對絕對值大于cB的殘差，其權(quán)重均為0。

求解基于M 估計的線性穩(wěn)健擬合需要使用迭代加權(quán)最小二乘法[14]，比較簡便的方法是直接采用Matlab、R 等常用數(shù)學(xué)軟件的內(nèi)置函數(shù)，本文采用Matlab軟件的robustfit函數(shù)[15]。

2.2 多普勒頻率線性穩(wěn)健擬合方法及步驟

由式（4）將多普勒頻率實測值fd(t)線性變換成fx(t)，對fx(t)進行線性穩(wěn)健擬合的方法及步驟如下：

1）對fx(t)進行第1 次穩(wěn)健擬合，得到系數(shù)K1、K2的第1 次估計值以及擬合直線fxn1(t)；由K1、K2計算ρ、tρ的第1 次估計值，再利用式（1）計算fd(t)擬合值fdn1(t)。

2）計算擬合殘差εxn(t)=fx(t)-fdn1(t)，再計算εxn(t)的中位數(shù)絕對離差MAD，公式為[16-17]：

式（7）中，MED 為中位數(shù)。

MED 是穩(wěn)健擬合中重要的位置測度，MAD 是重要的尺度測度，他們的抗野值能力都很強。為此，基于MAD 計算穩(wěn)健擬合標(biāo)準(zhǔn)差σxn=MAD/0.674 5[16]；最后，將擬合殘差絕對值＞3σxn的數(shù)據(jù)fx(tk)判別為野值，并將其修正為對應(yīng)的擬合值fxn1(tk)，得到新的fx(t)。

3）對新的fx(t)進行第2 次穩(wěn)健擬合，得到系數(shù)K1、K2的第2 次估計值以及擬合直線fxn2(t)，再計算ρ、tρ的最終估計值以及fd(t)的最終擬合值fdn2(t)。

3 擬合方法仿真驗證

3.1 僅含隨機誤差時的穩(wěn)健擬合結(jié)果

圖2 中，fd(t)數(shù)據(jù)僅含隨機誤差，曲線起伏較小。圖3 中，fx(t)為其線性轉(zhuǎn)換數(shù)據(jù)。僅執(zhí)行步驟1），對fx(t)進行穩(wěn)健擬合得到擬合直線fxn1(t)如圖4 所示?？梢?，fxn1(t)幾乎與理論值Fx(t)重合。同時,得到ρ 和tρ的估計值分別為10.025 m 和20.000 2 s，生成的擬合值fdn1(t)也非常逼近理論值Fd(t)（圖略）。

圖4 fd(t)僅含隨機誤差時的穩(wěn)健擬合結(jié)果Fig.4 Robust fitting result when fd(t) contains random errors

為進行比較，用蒙特卡洛法對fd(t)進行非線性最小二乘擬合[1-2，18]、fx(t)進行普通線性最小二乘擬合、以及fx(t)進行穩(wěn)健擬合（2 種權(quán)重函數(shù)）分別仿真1 000次。每次加入不同的白噪聲誤差ε(t)～N(0,)，σ0仍取800 Hz。 ρ 和tρ估計值的均值分別為10.000 m 和20.000 0 s，ρ 和tρ估計值的標(biāo)準(zhǔn)差（分別用σρ和σtρ表示）統(tǒng)計結(jié)果如表1所示。

表1 不同擬合方法ρ 和tρ 估計值的標(biāo)準(zhǔn)差統(tǒng)計結(jié)果Tab.1 Statistics results of standard errors of ρ and tρ in different fitting methods

分析表1 可知：僅含隨機誤差時，4 種擬合精度均很高，σρ均小于0.15 m 0，σtρ均小于0.3 ms；但穩(wěn)健擬合精度相對更高，尤其采用Bisquare 權(quán)重函數(shù)的本文方法。這是因為由式（6）可知，Bisquare函數(shù)對絕對值不大于cB的中間段殘差也進行權(quán)重分配，越靠近分布中心的殘差其權(quán)重越接近1。

3.2 包含多個野值時的穩(wěn)健擬合結(jié)果

主要討論成片型（或斑點型）野值[2，10]，且野值主要分布于曲線一側(cè)的情形，因為這種情形對普通最小二乘擬合最為不利。在圖2 所示fd(t)數(shù)據(jù)加入11 段成片野值（每段4～7個），野值數(shù)量比例為40%，加野值的fd(t)及其非線性擬合結(jié)果fdn(t)如圖5所示，可見擬合值fdn(t)與理論值Fd(t)偏離較大，此時ρ 和tρ估計值分別為11.650 m 和20.001 4 s。如果利用擬合殘差εdn(t)=fd(t)-fdn(t)識別野值，則會導(dǎo)致部分正常數(shù)據(jù)被誤判為野值。可見，普通最小二乘擬合受野值影響較大，或者說擬合穩(wěn)健性不強。

圖5 fd(t)含野值時的非線性擬合結(jié)果Fig.5 Nonlinear fitting result when fd(t) contains outliers

將加野值的fd(t)線性變換成fx(t)，對fx(t)先執(zhí)行步驟1）進行第一次穩(wěn)健擬合得到擬合直線fxn1(t)，如圖6 所示，fxn1(t)沒有很逼近Fx(t)，但偏離并不大，說明第一次擬合的穩(wěn)健性就很好。此時，ρ 和tρ的第一次估計值分別為9.790 m 和19.999 0 s，相比fd(t)非線性擬合的估計精度明顯提高。

圖6 fd(t)含野值時的第一次穩(wěn)健擬合結(jié)果Fig.6 First robust fitting result when fd(t) contains outliers

繼續(xù)執(zhí)行步驟2）和3），野值修正后的新fx(t)及其第二次穩(wěn)健擬合直線fxn2(t)，如圖7所示?？梢?，fxn2(t)非常逼近Fx(t)。 ρ 和tρ的最終估計值分別為9.945 m和19.999 5 s，比第一次穩(wěn)健擬合的估計精度進一步提高；同時，fd(t)最終擬合值fdn2(t)也非常逼近Fd(t)，滿足引信啟動和引戰(zhàn)配合分析需要。

圖7 野值修正后的穩(wěn)健擬合結(jié)果Fig.7 Robust fitting result after correcting outliers

4 結(jié)論

本文研究了引信多普勒頻率的線性穩(wěn)健擬合方法，將常用的多普勒頻率非線性擬合轉(zhuǎn)換為線性穩(wěn)健擬合。研究表明：該方法充分利用了線性穩(wěn)健擬合抗野值能力強的特點，多普勒頻率野值比例為40%時的擬合精度仍很高；同時，該方法僅需2次線性穩(wěn)健擬合和1次野值修正，步驟簡單、編程實現(xiàn)容易。方法可用于飛行試驗引信啟動和引戰(zhàn)配合結(jié)果分析，也可用于導(dǎo)彈脫靶量估算。