裴 喆

（92941部隊(duì)，遼寧葫蘆島125000）

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

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

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

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

圖1 彈目遭遇過(guò)程Fig.1 Missile-target encounter process

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

式（1）中：λ 為引信工作波長(zhǎng)；Vr通?？捎蓪?dǎo)引頭測(cè)量得到。

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

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

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

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

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

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

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

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

同理，將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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 擬合方法仿真驗(yàn)證

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

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

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

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

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

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

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

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

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

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

圖6 fd(t)含野值時(shí)的第一次穩(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所示?？梢?jiàn)，fxn2(t)非常逼近Fx(t)。 ρ 和tρ的最終估計(jì)值分別為9.945 m和19.999 5 s，比第一次穩(wěn)健擬合的估計(jì)精度進(jìn)一步提高；同時(shí)，fd(t)最終擬合值fdn2(t)也非常逼近Fd(t)，滿足引信啟動(dòng)和引戰(zhàn)配合分析需要。

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

4 結(jié)論

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