炮彈起始擾動姿軌耦合建模與微多普勒特征分析

宮志華 李開明 段鵬偉 陳春江

①(中國人民解放軍63861部隊 白城 137001)

②(空軍工程大學信息與導航學院 西安 710077)

1 引言

依據彈道學理論，炮彈目標運動包括質心運動和繞心轉動，據此構建的炮彈運動模型就是炮彈質心運動和繞心轉動的動力學方程[1-4]。炮彈繞心運動，亦稱為角運動，包括自旋、進動、章動、翻滾和擺動等，在武器試驗鑒定領域，這些角運動表征參數的獲取對武器系統的戰術性能研判、故障判定、射表編擬和目標辨識等具有重要意義。

2000年，美國海軍研究院的Chen[5]將目標質心平動以外的旋轉、振動等微小運動稱為微動(micromotion)，并將微動對雷達回波產生的附加頻率調制稱為微多普勒效應(micro-Doppler effect)。在高頻區，目標回波可被等效為一系列散射中心電磁散射的合成，通過各散射中心在空間的變化可近似反映目標的運動情況[6]。炮彈角運動使彈體上某些強散射中心發生空間變化，對雷達回波信號形成的頻率調制，也屬于典型的微多普勒效應，其產生的微多普勒特征是炮彈目標狀態判定的重要依據。因此，炮彈角運動的動力學模型是炮彈目標角運動參數與測量雷達回波信號映射模型的理論基礎，是描述炮彈角運動特征量與雷達信號之間關系的依據，可以反映炮彈角運動規律與目標回波之間的內在聯系，是炮彈目標特征參數辨識的重要參考。

近年來，雷達目標微多普勒效應研究是國內外雷達領域學術界和工程界的研究熱點之一[7-12]。文獻[10]提出一種基于物理光學法和等效電磁流法的擴展性微動目標回波模擬方法，并分析了圓錐與帶翼彈頭的進動特性；文獻[11]針對相位污染的缺損回波提取了空間碎片目標的微動特征；文獻[12]基于正交頻分復用信號外輻射源雷達，對帶旋翼低空目標進行回波建模與微動參數提取；文獻[13]給出了飛行動目標及其微動電磁散射特征分析的系統方法，借助軟件建模平臺建立了目標多自由度運動模型，對目標的雷達回波信號進行了仿真計算分析。實際應用場景中，對微多普勒效應進行建模分析時結合目標環境與運動特征能夠獲得更好的解釋性，從而對目標特性分析和識別更有價值[14]。同時，現有的對彈道目標的微動回波建模大多針對外太空飛行的彈道導彈目標展開[15]，彈道導彈目標的彈道中段通常在平流層，受大氣影響可以忽略，而對于炮彈、火箭彈等戰術類導彈、彈丸目標，其飛行基本在大氣層內，受空氣擾動影響明顯。因此，其動力學模型與彈道導彈存在明顯差異，常規的彈道導彈微動回波建模不再適用[1-4,6]。特別是，炮彈目標在彈道初始直線段受到起始擾動的影響[16]，其特殊的角運動對目標回波的調制作用與彈道導彈目標微動回波存在明顯差異。

基于以上分析，本文針對炮彈類目標，重點分析炮彈目標角運動機理，在此基礎上結合目標運動學方程建立目標角運動回波信號模型，從炮彈目標運動的物理原理出發，揭示目標微動與回波信號之間的內在關系，為有效分析炮彈目標的雷達測量信息、準確辨識炮彈目標運動參數提供參考依據。

2 炮彈起始擾動分析

在炮彈飛行過程中，根據彈丸在空中飛行的實際彈道軌跡，將彈道劃分為直線彈道段和曲線彈道段。相較于飛行全彈道，剛出炮口的一小段初始彈道，如圖1中紅色粗線段所示，彈丸的初速度較大，彈體結構偏心、質量偏心以及風等強迫擾動項對彈丸的作用還不顯著，即彈道傾角速度≈0、彈道傾角加速度≈0，一般而言彈道傾角下降量小于 30′，此時的彈道可近似為直線[1-4]，通常把這段離炮口不遠的一段彈道稱為“直線段”[2]。

圖1 彈道軌跡與直線彈道段示意圖Fig.1 Schematic diagram of ballistic and the straight line segment

在直線彈道段，炮彈發射后效期結束瞬間，僅受起始擾動的影響，炮彈的動力學方程可用復攻角運動的2階變系數齊次微分方程表示[17]：

起始擾動[4]，即當初始s=0時，由發射瞬間隨機產生的復攻角?0和復攻角速度引起的彈道攻角變化。采用“系數凍結法”求解式(1)，則復攻角?的解可寫為

由上述分析可知，由起始擾動引起的炮彈攻角變化是由兩個角頻率分別為ω1和ω2的快、慢圓運動復合成的，快圓運動是陀螺力矩作用的結果，慢圓運動是外力矩作用的結果。對于動態穩定良好的炮彈目標來說，必須滿足λ1<0,λ2<0，則攻角的兩圓運動的半徑不斷縮小，即模態振幅K1,K2的值不斷變小。此時，每個圓運動都成為半徑遞減的螺線，隨著攻角的不斷縮小，兩個圓運動逐漸合并，趨于單圓運動，再逐漸衰減而消失，復攻角的矢端將在復平面上畫出不斷縮小的外擺線或橢圓曲線。

3 雷達信號建模

3.1 坐標系解析

通過上述對炮彈在直線段僅受初始擾動影響的繞心運動特性分析，如圖2所示，可將炮彈角運動分解為3個運動的疊加：即繞彈軸矢量OA的自旋、彈軸繞陀螺動量矩矢量G(即章動軸)旋轉的章動(快圓運動)和陀螺動量矩矢量G繞速度矢量V1(即進動軸)旋轉的進動(慢圓運動)。

圖2 炮彈直線段角運動空間幾何關系Fig.2 Spatial geometric relation of projectile angular motion in the straight line segment

進一步，結合彈道學理論[2]和測量雷達相關坐標系，建立炮彈目標直線段角運動空間幾何關系，以雷達對炮彈目標的尾追式跟蹤為例，設炮彈目標自旋角速度為，章動角為?，章動角速度為ω1，進動角為ψ，進動角速度為ω2。

圖2中引入了以下4種坐標系，包括全局坐標系(即雷達坐標系)Q-UV W、炮口坐標系O1-xyz、參考坐標系O-XY Z和彈軸坐標系O-ξηζ。炮口坐標系以炮口中心為原點O1，水平軸O1x為射擊面與炮口水平面的交線，順時針射向為正，鉛直軸O1y在射擊面內并與水平軸O1z相垂直；雷達坐標系以站址中心為原點Q，水平軸QU在水平面內指向正北方向，QV為鉛垂線方向并與水平軸QW相垂直；參考坐標系以炮彈質心為原點O，與雷達坐標系始終平行；彈軸坐標系原點在炮彈質心上，Oξ軸沿彈軸向前為正，Oη軸垂直于彈軸向上為正，Oζ軸由右手法則確定。

設發射初始時刻，自旋軸OA、章動軸G和進動軸V1重合，參考坐標系原點與炮口坐標系原點重合，雷達坐標系Q-UV W、炮口坐標系O-xyz、參考坐標系O-XY Z和彈軸坐標系O-ξηζ之間的幾何關系如圖3所示。

圖3 初始時刻4種坐標系空間幾何關系Fig.3 Spatial geometric relation of the four coordinate systems as the projectile is stationary

此時，炮口坐標系Oz軸和彈軸坐標系Oζ軸重合，Oξη和Oxy共面，并繞Oζ軸逆時針旋轉θ0角，θ0稱為射角，即彈丸飛離炮口瞬間的彈軸矢量線與炮口水平面間的夾角(位于垂直面內)；炮口坐標系Oy軸和參考坐標系OY軸重合，Oxz和OXZ共面，并繞OY軸逆時針旋轉α0角，α0稱為射向角，即發射前炮口坐標系和參考坐標系水平面重合時，起始方向存在的夾角(位于水平面內)。則彈軸坐標系到參考坐標系的坐標變換矩陣可表示為

在發射瞬間，由于受到初始擾動的影響，彈軸出現章動和進動，章動角為?，進動角為ψ，如圖2所示。不考慮章動角和進動角初始相位，則彈軸坐標系向參考坐標系轉換的坐標變換矩陣可表示為

可見，∠AOV1值 即是攻角值δ。式(5)用于下文驗證回波建模的正確性。

3.2 回波信號建模

其中，根據Euler-Rodrigues繞矢量軸旋轉公式，各旋轉矩陣可展開為[5,18]

設雷達發射單頻連續波信號，其表達式為

其中，f0為載頻。則目標回波可表示為

將發射信號作為參考信號，與目標回波做相參處理，得到回波信號為

3.3 微多普勒分析

對式(10)中的相位項Φm(t)關于時間t求導，可得到回波的多普勒頻率為

由式(12)可以看出，由彈軸在空間的角運動引起的微多普勒表達式較為復雜[19]。進一步分析可得如下結論[20]：即雷達微多普勒頻率隨時間的變化具有周期性，彈頂散射點運動周期Ttop為章動周期Tu和進動周期Tc的最小公倍數；尾翼散射點運動周期Tbom為自旋周期Ts、進動周期Tc和章動周期Tu的最小公倍數，即有如下關系式成立：

其中，k1,k2,k3∈N，N為自然數集合。

4 仿真分析

為驗證第3節炮彈角運動對測量雷達回波信號的映射關系模型構建是否準確，分別選取典型高速自旋炮彈和低速自旋尾翼炮彈進行雷達微動信號建模和時頻分析。采用單頻連續波雷達測量炮彈目標，設載頻f0=10 GHz。在炮口坐標系下，射向角α0=5°，射角θ=53°。考慮到實際中炮彈的回波將以頂點和底面圓環邊緣點的回波為主，尾翼炮彈的回波將以頂點和尾翼散射點的回波為主，仿真中分別設置3個強散射點，分別為彈頂散射點1個、彈底邊緣散射點或尾翼散射點2個。

4.1 高速旋轉炮彈分析

如圖2所示，在彈軸坐標系O-ξηζ下，設彈體有3個強散射點：彈頂散射點A和彈底邊緣散射點B1,B2，仿真參數設置如表1所示。

表1 仿真參數設置(高速旋轉炮彈)Tab.1 Simulation parameters setting (high-speed rotating projectile)

根據式(12)可推算得，彈頂散射點運動周期理論值Ttop為0.50 s，彈底邊緣散射點運動周期理論值Tbom為0.50 s。考慮章動角和進動角衰減情況，仿真結果分別如圖4-圖6所示。

圖4 快圓和慢圓不衰減情況(高速旋轉炮彈)Fig.4 In the case of neither the fast nor the slow circular motion decays (high-speed rotating projectile)

圖4(a)和圖4(b)為運動補償前的目標微多普勒理論值和Gabor變換時頻圖像，圖4(d)和圖4(e)為進行運動補償后，得到的目標微多普勒理論值和Gabor變換時頻圖像。由圖4(d)和圖4(e)可以明顯看出，運動補償后彈底邊緣散射點由彈體自旋、章動和進動的共同作用產生微多普勒頻率調制，且包絡幅度主要受章動和進動的影響，彈頂散射點的微多普勒僅由章動和進動產生。由于炮彈自旋頻率高，時頻圖像出現自旋邊帶，且由于章動的調制，各散射點均出現鋸齒波動。由圖4(c)可見，當炮彈角運動的快圓和慢圓角度幅值均不衰減時，由彈頂散射點體現的攻角變化曲線呈現出明顯周期性，且幅度較為平穩，體現了快圓和慢圓運動在空間上的耦合，符合炮彈運動攻角的變化規律。對圖4(e)的時頻圖做自相關處理，其結果如圖4(f)所示，可見，自相關處理后得到的時頻圖像中的微多普勒信號周期為0.50 s，與理論值一致；進一步應用骨架提取和Hough變換，也可以從時頻圖像中提取微多普勒信號的周期，經過5,10,20次蒙特卡羅試驗獲得的彈頂和彈底邊緣微動周期均值如表2所示。

表2 Hough變換提取的時頻圖像周期(高速旋轉炮彈)Tab.2 Extracted period of time-frequency image with Hough Transform (high-speed rotating projectile)

從表2的統計結果可以看出，經過20次蒙特卡羅試驗得到時頻圖像中兩個微動信號的周期都約為0.49 s，同時可以發現，隨著試驗次數的增加，微動參數提取結果的均值逐漸逼近真實值。此時，相對應的彈頂散射點運動周期估計值約為0.49 s，彈底散射點運動周期估計值約為0.49 s，與理論值及自相關法提取的微動周期十分接近，且對比可以看出自相關法提取的微動周期更加精確，后續將主要采用自相關法進行高速自旋炮彈的微動周期提取。

進一步，考慮快圓角度幅值快速衰減的情況，設衰減規律為?′=?·exp(-3.5t)，由圖5(a)和圖5(b)可見，運動補償后，當章動幅值衰減時，彈底邊緣散射點自旋微多普勒僅被進動所調制，此時彈頂散射點的微多普勒僅由進動產生，各散射點的鋸齒波動均消失。同樣，對圖5(b)的時頻圖采用自相關法可得：彈頂散射點運動周期Ttop為0.50 s，彈底散射點運動周期Tbom為0.50 s，與理論值一致；由圖5(c)可見，由彈頂散射點體現的攻角幅值曲線經章動角幅值衰減后，攻角幅值整體呈現衰減趨勢，但均值整體接近進動角，符合炮彈運動攻角變化規律。

圖5 快圓衰減而慢圓不衰減情況(高速旋轉炮彈)Fig.5 In the case of the fast circular motion decays and the slow circular motion remains (high-speed rotating projectile)

當同時考慮快圓和慢圓角度幅值衰減的情況下，快圓衰減規律同上，設慢圓衰減規律為ψ′=ψ·exp(-0.05t)，由圖6(a)和圖6(b)可見，當進動角和章動角同時衰減時，由于進動角幅值相對較大，彈底邊緣散射點和彈頂散射點的微多普勒主要由自旋和進動產生，各散射點鋸齒波動均消失，進動衰減使彈頂和彈底散射點的微多普勒包絡幅值在不斷減小。對圖6(b)的時頻圖采用自相關法可得，彈頂散射點運動周期Ttop為0.50 s，彈底散射點運動周期Tbom為0.50 s，與理論值一致；由圖6(c)可見，由彈頂散射點體現的攻角幅值曲線經章動角幅值衰減后，同理攻角幅值整體呈現衰減趨勢，但上下包絡的中心表現為以進動角為初值的指數衰減形式，符合炮彈運動攻角變化規律。

圖6 快圓和慢圓均在衰減情況(高速旋轉炮彈)Fig.6 In the case of both the fast and the slow circular motion decays (high-speed rotating projectile)

4.2 旋轉尾翼炮彈分析

如圖2所示，在彈軸坐標系O-ξηζ下，設彈體有3個強散射點：彈頂散射點A和尾翼散射點B1,B2，仿真參數設置如表3所示。

表3 仿真參數設置(旋轉尾翼炮彈)Tab.3 Simulation parameters setting (rotating tail projectile)

根據式(12)可計算得到彈頂散射點運動周期的理論值Ttop為2.00 s，尾翼散射點運動周期的理論值Tbom為2.00 s。考慮章動角和進動角衰減情況，仿真結果分別如圖7-圖9所示。

圖7 快圓和慢圓不衰減情況(旋轉尾翼炮彈)Fig.7 In the case of neither the fast nor the slow circular motion decays (rotating tail projectile)

圖7(a)和圖7(b)為未進行運動補償的自旋尾翼炮彈微多普勒理論值和Gabor變換時頻圖像，圖7(d)和圖7(e)為進行運動補償后得到的微多普勒理論值和Gabor變換時頻圖像。由圖7(d)和圖7(e)可見，尾翼散射點的微多普勒同時被自旋、章動和進動調制，曲線包絡清晰可見，由于章動幅度偏小，其微多普勒受自旋和進動的影響較大；彈頂散射點的微多普勒主要由章動和進動產生。在快圓和慢圓角度幅值均不衰減的情況下，由圖7(c)可見，由彈頂散射點體現的攻角變化曲線體現了快圓和慢圓角運動在空間上的合成，符合炮彈運動攻角周期性的變化規律。同樣，對圖7(e)的時頻圖做自相關處理，可以得到時頻圖像中的微多普勒信號周期為2.00 s，如圖7(f)所示；進一步，應用骨架提取和Hough變換從時頻圖像中提取時頻圖像周期，同樣經過多次蒙特卡羅試驗進行微動周期的提取，表4分別給出5,10,20次蒙特卡羅試驗獲得的彈頂和尾翼的微動周期估計的均值。

表4 Hough變換提取的時頻圖像周期(旋轉尾翼炮彈)Tab.4 Extracted period of time-frequency image with Hough Transform (rotating tail projectile)

從表4的統計結果可以看出，經過20次蒙特卡羅試驗得到時頻圖像中兩個微動信號的周期估計值分別約為2.00 s和1.99 s，相對應的彈頂散射點運動周期Ttop約 為2.00 s，彈底散射點運動周期Tbom約為1.99 s，與理論值及自相關法提取的微動周期十分接近，且可以看出自相關法提取的微動周期更加精確，后續將主要采用自相關法進行自旋尾翼炮彈的微動周期提取。

在考慮快圓角度幅值快速衰減的情況下，衰減規律為?′=?·exp(-1.5t)，由圖8(a)和圖8(b)可見，當章動角幅值衰減后，尾翼散射點的微多普勒主要受自旋和進動的影響，彈頂散射點的微多普勒主要受進動的影響。對圖8(b)的時頻圖同樣采用自相關法可得，彈頂散射點運動周期Ttop為2.00 s，尾翼散射點運動周期Tbom為2.00 s，與理論值一致；由圖8(c)可見，由彈頂散射點體現的攻角幅值曲線經章動角幅值衰減過后，攻角幅值整體呈現衰減趨勢，但均值整體接近進動角，符合炮彈運動攻角變化規律。

圖8 快圓衰減而慢圓不衰減情況(旋轉尾翼炮彈)Fig.8 In the case of the fast circular motion decays and the slow circular motion remains (rotating tail projectile)

在同時考慮快圓和慢圓角度幅值均衰減的情況下，快圓衰減規律同上，慢圓衰減規律為ψ′=ψ·exp(-0.05t)，由圖9(a)和圖9(b)可見，當進動角和章動角幅值同時衰減后，由于進動角幅值相對較大，尾翼散射點的微多普勒包絡主要被自旋和進動調制，彈頂散射點微多普勒主要由進動產生，進動幅值的衰減使得彈頂和尾翼散射點的微多普勒包絡幅值在不斷減小。對圖9(b)的時頻圖采用自相關法可得，彈頂散射點運動周期為2.00 s，尾翼散射點運動周期為2.00 s，與理論值一致；由圖9(c)可見，由彈頂散射點體現的攻角幅值經章動角幅值衰減過后，攻角幅值整體呈現衰減趨勢，但上下包絡的中心表現為以進動角為初值的指數衰減形式，符合炮彈運動攻角變化規律。

圖9 快圓和慢圓均在衰減情況(旋轉尾翼炮彈)Fig.9 In the case of both the fast and the slow circular motion decays (rotating tail projectile)

5 實例驗證

首先給出某型炮彈和某型末制導炮彈(帶4個尾翼)尾追式跟蹤得到的直線段回波時頻分析仿真結果，考慮到實際炮彈目標的瞬時速度變化，圖10和圖11分別給出瞬時徑向速度隨時間的變化，仿真采用連續波雷達，載頻與實際裝備完全相同。

圖10 某型炮彈直線段回波時頻分析仿真結果Fig.10 Simulated time-frequency image of projectile in the straight line segment

圖11 某型末制導炮彈時頻分析仿真結果Fig.11 Time-frequency image of terminal guided projectile in the straight line segment

在上述理論和仿真分析的基礎上，進一步結合實測數據進行驗證。在試驗鑒定中，采用某型連續波雷達尾追式跟蹤測量，分別得到炮彈和末制導炮彈在初始直線彈道段的兩幅時頻圖像，分別如圖12和圖13所示。

圖13 某型末制導炮彈直線段回波真實時頻圖像Fig.13 Time-frequency image of terminal guided projectile in the straight line segment

從圖12可以看出，炮彈目標真實回波的時頻分析結果整體呈現衰減形式，這是由于目標遠離雷達的速度引起的，且對比圖10可以看出，實際目標回波受到遮擋效應的影響，時頻圖像主要表現為彈底散射點的多普勒，與圖10中間亮線的變化趨勢一致。彈底邊緣散射點的微多普勒由自旋、章動和進動共同作用產生，其微多普勒包絡主要受章動和進動的影響。由于炮彈的結構尺寸小、自旋頻率高，進動的調制形成比較明顯的包絡，而章動的調制對自旋微多普勒邊帶形成更為微小的鋸齒波動。隨著攻角逐漸變小，進動調制包絡也逐漸衰減。但整體而言，由于炮彈底面半徑較小且彈底散射點微動幅度小，其微多普勒受章動和進動的調制效果不明顯。

圖13給出末制導炮彈直線段真實回波的時頻圖像，同樣與圖11的仿真結果變化趨勢相同。但實際應用場景中，由于末制導炮彈自旋頻率低，尾翼散射點微多普勒受自旋、章動和進動的復合調制，同時進動的調制形成明顯的包絡，而章動幅值相對較小，章動調制效果不明顯。隨著攻角變小，調制包絡也逐漸衰減。由于帶尾翼的炮彈彈底散射點微動幅度較大，受進動調制效果明顯可見。

通過仿真分析與實際回波時頻分析結果的進一步對比可以發現，兩者變化的整體趨勢是一致的，證明了本文的建模方法的正確性和可行性。

6 結論

本文通過分析炮彈目標在直線段運動過程中彈軸和速度矢量之間角運動的變化機理，構建了炮彈角運動對雷達測量信息的投影變換模型，通過選取典型炮彈目標進行仿真分析和實測數據驗證得出如下結論：(1)炮彈目標受空氣阻力的影響，在隨彈道前進過程中，微動受到自旋、進動和章動的復合調制，且進動角和章動角的初值存在差異，隨大氣干擾而逐步衰減，與彈道導彈目標的微動建模存在明顯不同，兩者在彈道術語(如進動、章動等含義)方面的界定存在明顯差異；(2)由于炮彈自旋頻率較高，章動衰減快、幅度小，進動產生的微多普勒效果比較明顯，體現在雷達時頻圖像中進動特征比較明顯；末制導炮彈自旋頻率較低，章動衰減快、幅度小，由自旋和進動產生的微多普勒效果比較明顯，體現在雷達時頻圖像中自旋和進動特征明顯。實測數據驗證了所建立的基于炮彈目標角運動回波模型是正確有效的，符合炮彈運動規律。

以上研究結論為炮彈類目標的微動回波建模探索了新的途徑，為尋求彈道目標與炮彈目標微動特征差異進行了嶄新的嘗試，為后續利用實際雷達測量信息辨識炮彈目標運動特性、氣動特性、結構特性和參數匹配辨識提供了參考依據。

在實際工程應用中，由于受到雷達信號形式、雷達照射角度、目標散射特性、結構特性、信噪比和自然環境、氣象條件等多種因素的共同影響，基于各型炮彈目標的角運動雷達回波信號建模更為復雜，下一步需要持續開展相關研究，為炮彈這一類特殊目標的微動特征提取與識別提供參考。