膛內旋轉彈丸微運動建模研究

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

(1.上海無線電設備研究所,上海 201109;2.上海目標識別與環境感知工程技術研究中心,上海 201109)

0 引言

隨著高能裝藥、高密度裝填技術在現代高性能火炮中的應用,彈丸在膛內的發射環境越來越復雜。火炮連續發射過程中,彈丸在膛內做高速運動,產生的摩擦和高溫會造成火炮身管內部的燒蝕,導致彈丸受力不平衡。彈丸在發射過程中所承受的橫向過載有時能達到縱向過載的2倍以上,嚴重時會導致彈丸卡膛、引信瞎火或早炸等故障。因此彈丸膛內運動姿態研究對于引信設計、彈載機構設計、強度校核、卡膛故障分析、身管壽命評價等具有重要意義。

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

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

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

1 彈丸膛內微運動的回波信號建模

1.1 理論模型

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

圖1 膛內彈丸運動示意圖

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

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

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

彈丸質心初始位置O0在雷達坐標系oxyz中的坐標向量R0可表示為

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

擺動變換矩陣可以表示為

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

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

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

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

其中

正弦函數可以利用泰勒級數展開,表達式為

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

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

1.2 模型仿真分析

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

圖2 膛內彈丸縱向運動規律曲線

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

圖3 膛內彈丸運動規律曲線

根據雷達工作原理,多普勒頻率fD與雷達目標速度vc的關系可表示為

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

圖4 雷達回波信號仿真結果

2 時頻分析方法

2.1 短時傅里葉變換

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

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

式中:MC是模式修正參數,與發射的微波性質有關;ve是彈丸運動最大速度的估計值。

2.2 Wigner-Ville分布

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

式中:*表示共軛運算。

2.3 Hilbert-Huang變換

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

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

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

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

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

可得到第i個IMF分量的解析函數

其中

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

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

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

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

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

3 多普勒信號的時頻分析

3.1 時頻聚集性分析

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

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

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

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

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

圖5 回波信號EMD分解結果

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

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

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

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

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

圖9是采用三種時頻分析方法提取的信號時頻關系曲線。可以看出,WVD和H HT回波信號微多普勒特征的提取能力明顯優于STFT。

圖9 三種方法得到的時頻關系曲線

3.2 時頻方法的誤差分析

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

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

表2是不同C值的信號在三種時頻分析方法處理下,估計的信號頻率平均誤差、方差和最大誤差。可以看出在估計精度方面,HHT明顯優于WVD和STFT,WVD又明顯優于STFT,其中使用STFT得到的結果已經將微多普勒特征丟失了。

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

3.3 實際膛內彈丸回波信號分析

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

圖11 完整的雷達回波信號

圖12 截取的待處理雷達回波信號

圖13 三種方法處理的雷達回波信號時頻關系圖

4 結論

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

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

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

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