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

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

0 引言

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

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

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

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

1 彈丸膛內(nèi)微運動的回波信號建模

1.1 理論模型

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

圖1 膛內(nèi)彈丸運動示意圖

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

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

式中:?是彈丸做簡諧運動時的位移;B是振動的最大幅度;ωb是振動的頻率;φ是振動運動的初始相位。

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

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

擺動變換矩陣可以表示為

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

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

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

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

其中

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

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

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

1.2 模型仿真分析

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

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

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

圖3 膛內(nèi)彈丸運動規(guī)律曲線

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

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

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

2 時頻分析方法

2.1 短時傅里葉變換

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

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

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

2.2 Wigner-Ville分布

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

式中:*表示共軛運算。

2.3 Hilbert-Huang變換

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

a)IMF的極值點個數(shù)和過零點個數(shù)相等或者最多相差一個;

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

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

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

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

其中

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

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

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

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

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

3 多普勒信號的時頻分析

3.1 時頻聚集性分析

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

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

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

表1 不同C值信號的時頻方法M值比較

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

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

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

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

圖6 C=5時的回波信號時頻分布圖

圖7 C=30時的回波信號時頻分布圖

圖8 C=50時的回波信號時頻分布圖

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

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

3.2 時頻方法的誤差分析

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

圖10 三種時頻分析方法的頻率誤差對比

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

表2 三種時頻分析方法估計的不同C值信號瞬時頻率誤差對比

3.3 實際膛內(nèi)彈丸回波信號分析

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

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

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

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

4 結(jié)論

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

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

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

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