李美蘭, 吳寧寧, 王雄志, 黃嘉智,, 韓永金

(1.上海無(wú)線電設(shè)備研究所,上海 201109;2.上海目標(biāo)識(shí)別與環(huán)境感知工程技術(shù)研究中心,上海 201109)

0 引言

隨著高能裝藥、高密度裝填技術(shù)在現(xiàn)代高性能火炮中的應(yīng)用,彈丸在膛內(nèi)的發(fā)射環(huán)境越來(lái)越復(fù)雜。火炮連續(xù)發(fā)射過(guò)程中,彈丸在膛內(nèi)做高速運(yùn)動(dòng),產(chǎn)生的摩擦和高溫會(huì)造成火炮身管內(nèi)部的燒蝕,導(dǎo)致彈丸受力不平衡。彈丸在發(fā)射過(guò)程中所承受的橫向過(guò)載有時(shí)能達(dá)到縱向過(guò)載的2倍以上,嚴(yán)重時(shí)會(huì)導(dǎo)致彈丸卡膛、引信瞎火或早炸等故障。因此彈丸膛內(nèi)運(yùn)動(dòng)姿態(tài)研究對(duì)于引信設(shè)計(jì)、彈載機(jī)構(gòu)設(shè)計(jì)、強(qiáng)度校核、卡膛故障分析、身管壽命評(píng)價(jià)等具有重要意義。

目前獲取彈丸膛內(nèi)高速運(yùn)動(dòng)參數(shù)的測(cè)試方法主要有高速攝影測(cè)量法[1-2]和激光杠桿測(cè)量法[3-4]等。高速攝影測(cè)量法由于空間分辨率有限,在測(cè)量精度上難以滿(mǎn)足要求。激光杠桿測(cè)量法對(duì)反射鏡面的大小有要求,對(duì)于尖頭彈丸,其頂部反射鏡的鏡面太小,可能無(wú)法接收和反射激光信號(hào)。增大激光光束直徑可以保證光斑面積足夠大,但又會(huì)降低光束質(zhì)量。同時(shí)彈丸在膛內(nèi)做高速運(yùn)動(dòng)使空氣壓縮,局部空氣的折射率發(fā)生改變,也會(huì)導(dǎo)致測(cè)量精度偏低。

文獻(xiàn)[5]針對(duì)旋轉(zhuǎn)機(jī)械狀態(tài)監(jiān)測(cè)過(guò)程中,因無(wú)線通信設(shè)備間相對(duì)位置發(fā)生變化而產(chǎn)生多普勒效應(yīng)的問(wèn)題進(jìn)行了研究,對(duì)轉(zhuǎn)速與多普勒頻移之間的關(guān)系進(jìn)行了仿真,但是沒(méi)有考慮旋轉(zhuǎn)線速度與無(wú)線通信設(shè)備發(fā)射信號(hào)垂直時(shí)的情況。文獻(xiàn)[6]提出了基于光學(xué)杠桿測(cè)量系統(tǒng)和毫米波雷達(dá)的彈丸膛內(nèi)姿態(tài)與縱向運(yùn)動(dòng)聯(lián)合測(cè)試方法,分析了彈丸膛內(nèi)擺動(dòng)和炮口振動(dòng)的頻譜分布和時(shí)頻特性。但是該文獻(xiàn)未指明炮彈的類(lèi)型和徑向速度范圍。文獻(xiàn)[7]將剛體運(yùn)動(dòng)學(xué)知識(shí)與多普勒效應(yīng)相結(jié)合,建立了微多普勒雷達(dá)回波信號(hào)模型,通過(guò)仿真驗(yàn)證了其準(zhǔn)確性。但是該文獻(xiàn)只考慮了一個(gè)假想的簡(jiǎn)單仿真環(huán)境,并沒(méi)有針對(duì)實(shí)際信號(hào)進(jìn)行分析,也沒(méi)有從信號(hào)里提取出真實(shí)的擺動(dòng)參數(shù)。

本文針對(duì)高速攝影法和光學(xué)杠桿法在彈丸膛內(nèi)運(yùn)動(dòng)姿態(tài)測(cè)試中存在的問(wèn)題,利用雷達(dá)微多普勒測(cè)試技術(shù),對(duì)靶場(chǎng)試驗(yàn)獲取的回波數(shù)據(jù)進(jìn)行分析,建立彈丸膛內(nèi)運(yùn)動(dòng)的微多普勒模型,采用多種時(shí)頻分析方法分析彈丸回波信號(hào)的頻譜分布和時(shí)頻特性,為膛內(nèi)彈丸的運(yùn)動(dòng)姿態(tài)研究提供參考。

1 彈丸膛內(nèi)微運(yùn)動(dòng)的回波信號(hào)建模

1.1 理論模型

線膛炮身管內(nèi)壁上刻有陽(yáng)線和陰線,彈丸在線膛炮身管內(nèi)做旋轉(zhuǎn)前進(jìn)運(yùn)動(dòng)。為了研究彈丸在膛內(nèi)的姿態(tài),需要簡(jiǎn)化模型,在這里只考慮彈丸的自旋和振動(dòng)兩種微運(yùn)動(dòng)。圖1為彈丸運(yùn)動(dòng)示意圖。

圖1 膛內(nèi)彈丸運(yùn)動(dòng)示意圖

oxyz坐標(biāo)系是以雷達(dá)質(zhì)心所在位置為坐標(biāo)原點(diǎn)建立的坐標(biāo)系。OXYZ坐標(biāo)系是以彈丸的質(zhì)心為原點(diǎn)建立的坐標(biāo)系,并隨著彈丸的運(yùn)動(dòng)平移,平移后的坐標(biāo)系記為OiXiYiZi。彈丸在雷達(dá)坐標(biāo)系中的方位角與俯仰角分別為α和β。因?yàn)閺椡柙谔艃?nèi)是緊貼著身管壁的,可以認(rèn)為彈丸在繞著其自身的中心軸旋轉(zhuǎn),即做自旋運(yùn)動(dòng)。彈丸旋轉(zhuǎn)角速度ωr與彈丸的平動(dòng)速度v之間的關(guān)系可以表示為

式中:θ是線膛炮的纏角;d是線膛炮口徑。將彈丸在膛內(nèi)的振動(dòng)簡(jiǎn)化成彈丸繞中心軸的簡(jiǎn)諧運(yùn)動(dòng),運(yùn)動(dòng)方程為

式中:?是彈丸做簡(jiǎn)諧運(yùn)動(dòng)時(shí)的位移;B是振動(dòng)的最大幅度;ωb是振動(dòng)的頻率;φ是振動(dòng)運(yùn)動(dòng)的初始相位。

彈丸質(zhì)心初始位置O0在雷達(dá)坐標(biāo)系oxyz中的坐標(biāo)向量R0可表示為

式中:R0為彈丸質(zhì)心初始位置O0到雷達(dá)坐標(biāo)系坐標(biāo)原點(diǎn)o的距離。設(shè)彈丸上微運(yùn)動(dòng)的參考點(diǎn)M在彈丸初始OXYZ坐標(biāo)系中的坐標(biāo)為(X0,Y0,Z0),t時(shí)刻彈丸自旋運(yùn)動(dòng)的旋轉(zhuǎn)矩陣可以表示為

擺動(dòng)變換矩陣可以表示為

t時(shí)刻雷達(dá)到彈丸上的點(diǎn)Mi的距離可表示為

式中:γ為雷達(dá)視線與彈丸縱向速度方向的夾角。由于R0?vt且R0?OiMi,對(duì)應(yīng)的α和β很小,即余弦值接近于1。設(shè)雷達(dá)發(fā)射的電磁波頻率為fc,則彈丸上散射點(diǎn)Mi產(chǎn)生的回波信號(hào)可以表示為

式中:Ai為點(diǎn)Mi產(chǎn)生的回波信號(hào)的幅度;c是電磁波傳播的速度。

由式(7)可知,RMi(t)決定著彈丸回波多普勒的特征。式(6)可簡(jiǎn)化為

其中

正弦函數(shù)可以利用泰勒級(jí)數(shù)展開(kāi),表達(dá)式為

只將最高的3次冪保留下來(lái),簡(jiǎn)化得到RMi(t)的表達(dá)式為

將式(2)代入式(16),可得對(duì)應(yīng)的點(diǎn)Mi產(chǎn)生 的多普勒頻率fD為

1.2 模型仿真分析

假設(shè)線膛炮身管的口徑是130 mm,彈丸在膛內(nèi)運(yùn)動(dòng)的時(shí)間為5 ms,最大速度約為1 000 m/s,纏度為1/20。根據(jù)內(nèi)彈道學(xué)方程組,可得到符合實(shí)際膛內(nèi)彈丸縱向運(yùn)動(dòng)規(guī)律的速度關(guān)系曲線,如圖2所示。

圖2 膛內(nèi)彈丸縱向運(yùn)動(dòng)規(guī)律曲線

使用C=k B來(lái)定量表示擺動(dòng)幅度的大小,稱(chēng)為C值。其中:k為理論仿真與現(xiàn)實(shí)信號(hào)能量的調(diào)節(jié)系數(shù),由觀察點(diǎn)的位置及雷達(dá)信號(hào)發(fā)射方向決定;B是簡(jiǎn)諧運(yùn)動(dòng)方程中的最大幅度,由觀察點(diǎn)的位置及雷達(dá)信號(hào)發(fā)射方向決定。同時(shí)C值也定量地表征了微多普勒能量占比的大小,C越大則擺動(dòng)幅度就越大,微多普勒能量占比就越大。采用Matlab進(jìn)行仿真,只考慮一個(gè)觀察點(diǎn),假設(shè)k=1/(50 000π),根據(jù)式(17)得到不同C值條件下彈丸在膛內(nèi)的運(yùn)動(dòng)規(guī)律仿真曲線,如圖3所示。在縱向運(yùn)動(dòng)規(guī)律曲線上加載的小波動(dòng)正是微多普勒的特征,隨著C值的增大,微多普勒特征更加明顯。

圖3 膛內(nèi)彈丸運(yùn)動(dòng)規(guī)律曲線

根據(jù)雷達(dá)工作原理,多普勒頻率fD與雷達(dá)目標(biāo)速度vc的關(guān)系可表示為

根據(jù)彈丸的運(yùn)動(dòng)規(guī)律,仿真得到雷達(dá)的回波信號(hào)如圖4所示。隨著時(shí)間的增大,波形越來(lái)越密,即頻率越來(lái)越高。

圖4 雷達(dá)回波信號(hào)仿真結(jié)果

2 時(shí)頻分析方法

2.1 短時(shí)傅里葉變換

短時(shí)傅里葉變換(STFT)是一種常用的線性時(shí)頻分析方法,采用分段的思想在時(shí)域上將完整的信號(hào)分成多段,每一段信號(hào)表示某一時(shí)刻的信號(hào)特征。信號(hào)x(t)的STFT表達(dá)式為

式中:w(t)是窗函數(shù)。為了兼顧時(shí)間分辨率和頻率分辨率,窗函數(shù)需要選擇一個(gè)合適的長(zhǎng)度。根據(jù)前人研究的經(jīng)驗(yàn)及實(shí)際處理情況,得到雷達(dá)測(cè)膛內(nèi)彈丸運(yùn)動(dòng)速度的最佳窗長(zhǎng)WD的計(jì)算公式為

式中:MC是模式修正參數(shù),與發(fā)射的微波性質(zhì)有關(guān);ve是彈丸運(yùn)動(dòng)最大速度的估計(jì)值。

2.2 Wigner-Ville分布

Wigner-Ville分布(WVD)是一種非線性的時(shí)頻分析方法,基本思路是對(duì)信號(hào)的自相關(guān)函數(shù)做傅里葉變換。WVD在一定程度上克服了STFT的海森堡測(cè)不準(zhǔn)問(wèn)題,但是它對(duì)包含兩個(gè)及以上頻率的信號(hào)進(jìn)行分析時(shí)會(huì)存在交叉項(xiàng),干擾正確的頻率信息。信號(hào)x(t)的WVD變換公式為

式中:*表示共軛運(yùn)算。

2.3 Hilbert-Huang變換

Hilbert-Huang變換(HHT)突破了傅里葉變換的理論條件限制。它的實(shí)現(xiàn)需要兩個(gè)步驟:首先對(duì)信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)得到固有模態(tài)函數(shù)(IMF),然后對(duì)IMF分量進(jìn)行希爾伯特變換。固有模態(tài)分量必須滿(mǎn)足兩個(gè)條件:

a)IMF的極值點(diǎn)個(gè)數(shù)和過(guò)零點(diǎn)個(gè)數(shù)相等或者最多相差一個(gè);

b)在任意時(shí)刻,由極大值構(gòu)成的上包絡(luò)線和極小值構(gòu)成的下包絡(luò)線關(guān)于橫軸對(duì)稱(chēng),即均值為0。

經(jīng)過(guò)模態(tài)分解后可以得到各個(gè)IMF分量的時(shí)頻關(guān)系,從中找到包含有效信息最多的分量進(jìn)行希爾伯特變換。對(duì)第i個(gè)IMF分量ci(t)做希爾伯特變換,表達(dá)式為

式中:h(t)為希爾伯特變換響應(yīng)函數(shù);?表示卷積運(yùn)算。

可得到第i個(gè)IMF分量的解析函數(shù)

其中

式中:ai(t)為幅值函數(shù);φi(t)為相位函數(shù)。

由相位函數(shù)φi(t)可以得到瞬時(shí)頻率

則ci(t)的希爾伯特譜可表示為

式中:n為IMF分量的總數(shù);Re·（）用于求解復(fù)數(shù)函數(shù)的實(shí)數(shù)部分。

相比前兩種方法,HHT方法獲得的瞬時(shí)頻率信息更精確,但是計(jì)算量偏大。

3 多普勒信號(hào)的時(shí)頻分析

3.1 時(shí)頻聚集性分析

使用時(shí)頻聚集性衡量指標(biāo)M值來(lái)客觀評(píng)價(jià)各種時(shí)頻分析方法的性能。M值的表達(dá)式為

式中:fTED(t,ω)為對(duì)應(yīng)時(shí)頻分析方法的分布函數(shù),其中ω為信號(hào)角頻率。M值越大,表示時(shí)頻聚集性越好,反之則時(shí)頻聚集性越差。為了真實(shí)反映時(shí)頻分析方法的性能,對(duì)不同條件下的信號(hào)進(jìn)行處理。表1是C值不同的情況下,各種時(shí)頻方法的M值比較結(jié)果。

由表1可知:STFT方法的時(shí)頻聚集性隨著C值的增大而變差;相比較而言,WVD方法和HHT方法的時(shí)頻聚集性雖然在C值小的情況下較差,但是在C值大的情況下優(yōu)于STFT方法。

表1 不同C值信號(hào)的時(shí)頻方法M值比較

對(duì)于H HT方法,需要先進(jìn)行EMD分解,再選擇合適的IMF分量進(jìn)行分析。圖5是信號(hào)經(jīng)過(guò)EMD分解得到的c1(t),c2(t)兩個(gè)分量的時(shí)間與幅值和時(shí)間與頻率的關(guān)系圖?？芍?第一個(gè)分量c1(t)所占的能量大,并且時(shí)頻關(guān)系符合整個(gè)信號(hào)時(shí)頻關(guān)系的走向趨勢(shì),所以主要對(duì)c1(t)分量進(jìn)行H HT變換分析。

圖5 回波信號(hào)EMD分解結(jié)果

圖5 回波信號(hào)EMD分解結(jié)果(續(xù))

圖6～圖8是C值分別為5,30,50情況下,STFT、WVD和H HT三種方法的時(shí)頻分布三維圖?？梢悦黠@看出,隨著C值的增大,HHT和WVD的時(shí)頻分布中代表微多普勒特征的微小波動(dòng)越來(lái)越明顯,而STFT只能獲取多普勒的普通輪廓,不能從中分析得到微多普勒的細(xì)節(jié)信息。因此,STFT不適用于膛內(nèi)彈丸的運(yùn)動(dòng)姿態(tài)研究。

圖6 C=5時(shí)的回波信號(hào)時(shí)頻分布圖

圖7 C=30時(shí)的回波信號(hào)時(shí)頻分布圖

圖8 C=50時(shí)的回波信號(hào)時(shí)頻分布圖

圖9是采用三種時(shí)頻分析方法提取的信號(hào)時(shí)頻關(guān)系曲線?？梢钥闯?WVD和H HT回波信號(hào)微多普勒特征的提取能力明顯優(yōu)于STFT。

圖9 三種方法得到的時(shí)頻關(guān)系曲線

3.2 時(shí)頻方法的誤差分析

分析理論仿真得到的速度信號(hào)與通過(guò)三種時(shí)頻分析方法得到的信號(hào)的時(shí)頻關(guān)系,可以得到頻率的誤差曲線,如圖10所示。

圖10 三種時(shí)頻分析方法的頻率誤差對(duì)比

表2是不同C值的信號(hào)在三種時(shí)頻分析方法處理下,估計(jì)的信號(hào)頻率平均誤差、方差和最大誤差?？梢钥闯鲈诠烙?jì)精度方面,HHT明顯優(yōu)于WVD和STFT,WVD又明顯優(yōu)于STFT,其中使用STFT得到的結(jié)果已經(jīng)將微多普勒特征丟失了。

表2 三種時(shí)頻分析方法估計(jì)的不同C值信號(hào)瞬時(shí)頻率誤差對(duì)比

3.3 實(shí)際膛內(nèi)彈丸回波信號(hào)分析

圖11是130 mm口徑的榴彈炮炮射實(shí)驗(yàn)獲得的雷達(dá)回波實(shí)測(cè)信號(hào)曲線。圖12是截取的待處理雷達(dá)回波信號(hào)曲線。圖13是采用STFT、WVD和HHT方法處理后得到的雷達(dá)回波信號(hào)時(shí)頻分布局部放大圖。從圖13中可以看出:STFT處理得到的信號(hào)輪廓時(shí)頻關(guān)系效果很好,說(shuō)明STFT適用于求解膛內(nèi)彈丸的縱向運(yùn)動(dòng)速度;采用WVD得到的信號(hào)時(shí)頻關(guān)系輪廓最模糊,可見(jiàn)其既不適合求解縱向運(yùn)動(dòng)規(guī)律,也很難提取到微多普勒特征;采用HHT處理的結(jié)果存在明顯的微多普勒特征,但是由于實(shí)際信號(hào)存在噪聲,導(dǎo)致波動(dòng)的規(guī)律性不強(qiáng)。綜合考慮回波信號(hào)分析結(jié)果,采用HHT進(jìn)行信號(hào)微多普勒特征提取的效果較好,WVD效果最差,STFT會(huì)將信號(hào)微多普勒特征丟失。

圖11 完整的雷達(dá)回波信號(hào)

圖12 截取的待處理雷達(dá)回波信號(hào)

圖13 三種方法處理的雷達(dá)回波信號(hào)時(shí)頻關(guān)系圖

4 結(jié)論

針對(duì)膛內(nèi)彈丸的微運(yùn)動(dòng)姿態(tài)測(cè)試問(wèn)題,本文建立了雷達(dá)監(jiān)測(cè)彈丸在膛內(nèi)運(yùn)動(dòng)的模型,分別用STFT、WVD和H HT三種時(shí)頻分析方法處理理論仿真信號(hào)和實(shí)測(cè)信號(hào),對(duì)處理后的結(jié)果進(jìn)行了對(duì)比與誤差分析。得出三點(diǎn)結(jié)論:

a)STFT處理縱向速度的結(jié)果明顯優(yōu)于WVD和HHT,但是隨著微多普勒能量的增大,其時(shí)頻聚集性變差;

b)相比STFT和WVD,HHT處理得到的瞬時(shí)頻率信息更多,能檢測(cè)到微多普勒信息,隨著微多普勒能量的增大,它的優(yōu)勢(shì)更加明顯;

c)對(duì)于微多普勒能量占比大的信號(hào),HHT處理得到的信號(hào)精確度更高,但對(duì)于信噪比低的信號(hào),HHT精確度不高,需要進(jìn)一步優(yōu)化。