一種基于微波光子超高分辨雷達機翼振動參數估計方法

符吉祥 邢孟道* 徐 丹 王安樂(西安電子科技大學雷達信號處理國家重點實驗室 西安 710071)

②(西安電子科技大學信息感知技術協同創新中心 西安 710071)

③(空軍預警學院 武漢 430019)

1 引言

得益于全天時全天候和超遠距離探測的優勢，雷達被廣泛應用到軍事和民事的應用中，如戰場監視和遙感測繪等。逆合成孔徑雷達，因為其能對非合作目標進行探測和成像，獲得目標的形狀、尺寸、姿態和運動信息，進而可以對非合作目標進行識別，在國防安全領域有著重要應用，受到各個國家的高度重視和廣泛研究。隨著雷達技術和需求的不斷提高，雷達向高分辨、多功能和多任務方向發展。然而，傳統微波雷達受限于微波器件的“電子瓶頸”[1]，無法產生大帶寬信號，使得雷達分辨率受限。近年來，微波光子技術的提出打破了這一瓶頸。微波光子技術將光子技術應用到雷達信號的產生和調制階段，充分利用了光子的大帶寬低損耗和抗電磁干擾的優勢，使得雷達能夠發射超大帶寬的高質量信號。美國自然雜志評價其為照亮雷達未來的技術[2,3]。目前，已經研制成功的微波光子雷達所能產生的最大信號帶寬為12 GHz[4]，這大大超過傳統的微波雷達。超高的分辨率，使得雷達能夠對目標進行更加精細的探測和感知，如目標微振動的觀測和提取，多頻段同時多任務觀測等。然而，分辨率的大幅度提升也給雷達成像和參數估計帶來了新的挑戰。一是高分辨率會帶來嚴重的越距離單元徙動[5]；二是運動誤差補償的精度提高，為此需要研究新的雷達成像和參數估計方法。

圖1 微波光子實測飛機目標1維距離像Fig.1 One dimensional range profile of airplane measured by microwave photonic radar

在對10 G帶寬的微波光子雷達實測民航飛機的數據處理過程中，我們發現，飛機目標在轉向過程中，因為目標狀態的改變，機翼與飛機的主體存在運動不一致的現象，如圖1所示的飛機目標的高分辨1維距離像，圖1(a)是飛機整體包絡對齊后的1維距離像包絡，機身的包絡如圖1(b)所示，可以發現，機身散射點包絡變化較為平穩，部分機翼的包絡如圖1(c)所示，可以發現機翼回波的包絡線存在明顯的抖動現象。機翼的這種非平穩抖動使得傳統的成像方法無法對其進行聚焦成像。機翼的振動相對于平動以及機身相對于雷達的轉動來說，其具有幅度小頻率高的特點，所以可以將其視為飛機的微動特征。利用微動特征可以評估飛機的運動狀態以及飛機的安全狀態，進一步為飛機目標的識別和監測作為依據。目前微動參數提取的相關方法多是應用在彈道類目標上，因為彈頭散射結構較為簡單，散射點較少，所以其可以利用基于時頻分析的方法進行微動參數估計[6,7]。具體做法是，首先利用包絡對齊平動補償[8,9]提取單個散射點回波，再通過對其進行時頻分析獲得其時間頻率分布曲線，最后通過對時頻曲線進行相關處理即可得到彈頭的微動參數。但是對于飛機目標來說，一方面，機翼結構較彈頭要復雜得多，純凈的單個散射點回波較難獲取，而其他散射點的回波會對參數估計造成嚴重干擾；另一方面，因為機翼回波的信噪比較低，直接采用時頻分析很難從其時頻分布圖中提取時頻曲線。基于此，本文提出一種基于子孔徑序列圖的振動參數估計方法。首先將機翼的振動建模為單擺運動，然后通過分析信號波數譜表達式，得到短時間機翼散射點聚焦位置和單擺運動參數之間的關系，進一步通過機翼子孔徑成像結果序列圖中散射點位置的變化來反演單擺運動的參數。最后提出一種修正的極坐標格式算法(Modified Polar Format Algorithm , MPFA)，通過非均勻傅里葉變換和估計的振動參數對振動的機翼進行高效的聚焦成像，并通過構造成像結果熵值優化函數對振動參數進行精估計。因為是利用圖像序列圖來進行擺動參數的反演，這充分地利用了2維相干積累增益，使得在低信噪比下的參數估計成為可能。另外相比于傳統的時頻分析方法，一方面，序列圖能夠提供距離和多普勒兩維的信息；另一方面，序列圖中各個散射點之間是解耦的，沒有交叉項干擾的影響。

2 機翼振動建模分析和振動參數估計

2.1 運動模型分析

根據力學分析可知，機翼的振動與單擺運動存在一致性，所以將其建模為單擺運動。如圖2所示，以飛機整體所在平面作為XOY平面，X軸指向為在參考時刻沿著機身方向指向機頭的方向，Z軸為飛機平面的法線方向。雷達視線與Z軸的夾角為φ0，其在XOY平面上的投影向量與Y軸之間的夾角為β0。機翼沿著Z軸振動，振動過程中，機翼與XOY平面之間的夾角為?，最大夾角為?0，根據單擺轉角公式可知，夾角隨慢時間正弦變化，即

其中，?0表示振動角，fv表示振動頻率，φ0表示振動初相，表示方位慢時間。在飛機沿著圖2中曲線軌跡運動的時候，飛機整體還有相對于雷達的平動和轉動。對于機身上的散射點p1，因為其不存在振動，所以其相對于雷達的運動可以表示為

其中，R0表示飛機旋轉中心O到雷達的初始距離，vr和ar分別表示飛機相對于雷達運動的徑向速度和徑向加速度。轉動分量是隨散射點位置發生空變的，其形式為(平面波近似)

其中，(xp1，yp1)表示散射點p1參考時刻的坐標位置，θ表示飛機相對于雷達轉動的角度，其是慢時間的函數，即θ(tm)=ωtm,ω表示轉速。對于機翼上的散射點，其不僅包含以上的運動，還包含機翼振動對其的運動調制。通過運動分解可以得到機翼上散射點p2各個時刻的坐標位置如式(5)

圖2 飛機運動幾何模型Fig.2 Geometry model of airplane

其中，(xp2，yp2)表示散射點p2參考時刻的坐標位置。所以散射點p2到雷達的斜距變化為

假設雷達發射線性調頻信號，使用去斜方式接收回波，經過剩余視頻項補償之后，機翼散射點雷達回波信號的波數譜表達式為[10]

2.2 振動參數估計方法

首先通過粗成像對機翼機身回波進行分離。因為機翼存在周期性振動而機身只存在平穩的運動，所以在粗成像結果中，機翼和機身上的散射點散焦情況存在明顯的差異，即機翼散射點散焦成線性分布而機身散射點散焦呈現塊狀分布。另外，從粗成像的輪廓中也可以輕易地將機翼和機身進行區分。如圖4所示，圖4(a)是飛機整體粗成像結果，可以發現，盡管飛機散焦嚴重，但是仍能夠從輪廓中輕易區分機身和機翼，圖4(b)是機身散射點塊狀散焦結果，圖4(c)是機翼散射點線性散焦結果。最后通過加窗濾波將成像結果中的機身回波濾出，從而實現機身機翼回波分離。

對機身回波使用keystone算法[13]校正其線性距離走動，得到

其中，sinc(x)=sin(x)/(x),表示徑向波數中心。觀察式(9)可以發現，keystone變換之后，包絡已經對齊，然而散射點回波存在和其距離坐標yp線性相關的方位向二次相位。通過對不同距離單元信號進行調頻率估計并對調頻率進行線性擬合即可得到轉速ω，然后再構造調頻率補償函數

圖3 轉臺模型Fig.3 Turntable model

其中，r表示式(9)中各個距離單元的位置。對式(9)補償之后直接方位向FFT即可得到聚焦圖像

從式(11)中可以發現，成像結果中散射點的聚焦位置相對于其在XOY中的坐標發生了旋轉和縮放，旋轉角度為方位角β0，縮放因子為sinφ0。所以可以從定標圖像中，通過計算機身相對水平方向的夾角來估計β0，通過測量定標圖像中機身長度與寬度與飛機真實的尺寸的比值來估計φ0。

對分離后的機翼部分回波進行平動補償，將其旋轉中心和振動中心補償為機翼中心點O′，如圖3(b)所示，再對其進行子孔徑成像，子圖像的距離和多普勒如式(12)

圖4 飛機粗成像結果Fig.4 Coarse imaging results of airplane

其中，Δyq的值可由以下方程組解得

觀察機翼散射點回波，發現其距離和方位是解耦的，利用極坐標格式算法的思想，本文提出一種修正的極坐標格式算法，即

將式(17)代入到式(8)中的機翼回波后得到

對式(18)進行2維傅里葉變換即可得到距離的機翼圖像。為了提高算法的效率，使用快速非均勻傅里葉變換(Non?Uniform Fast Fourier Transform,NUFFT)[14,15]高效地實現插值成像的操作。

其中，E[·]表示計算熵值操作，表示圖像矩陣，矩陣大小為M×N,si為 其第i個 元素。表示圖像的總能量。所以構造的最小熵優化函數如式(20)

其中，wing(kr，tm)表示機翼波數域回波，表示利用參數對進行非均勻快速傅里葉變換。使用粗估計的振動參數作為初值，利用式(20)進行參數搜索，得到振動參數的精確估計。這里采用的搜索策略是利用坐標下降法的思想，每次只對1個參數進行搜索，其他參數固定不變，循環迭代，直至相鄰兩次迭代的圖像熵值變化小于門限停止迭代。算法的流程圖如圖5所示。

3 仿真實驗及實測數據驗證

本文使用Ka波段雷達參數對飛機散射點模型進行參數估計，為了簡化仿真，將飛機兩側機翼振動參數設置為相同的，仿真參數如表1所示。飛機模型長35 m，寬30 m，如圖6所示。雷達視線方位角為80°，俯仰角為85°。

圖5 算法流程圖Fig.5 Flow chart of the proposed algorithm

表1 仿真雷達參數Tab.1 Simulation radar parameters

圖6 飛機散射點模型Fig.6 Scatter model of airplane

圖7 機身聚焦成像定標結果Fig.7 Focused and scaled result of airplane body

圖8 不同時刻機翼子孔徑成像結果Fig.8 Sub?aperture imaging results of airplane wing at different time

圖7是機身聚焦成像定標結果，從圖7中可以看出機身主軸和Y軸存在夾角，通過之前的分析可知，其夾角大小和雷達視線的方位角是相等的，通過機頭標記點A的坐標位置，可得雷達視線的方位，通過測量飛機寬度BC與實際中機身寬度之間的比例，可以得到對俯仰角的估計，即，其中W0為飛機的真實寬度。

圖8是機翼不同子孔徑成像結果，可以發現其方位多普勒存在周期性的變化。以圖8(a)中用圓圈標記的散射點1,2,3為例，分別提取其不同子孔徑方位多普勒和距離的變化，提取的多普勒變化曲線如圖9(a)所示，可以發現曲線存在明顯的階梯跳躍，這是因為子孔徑圖像多普勒分辨率有限，不能反映多普勒的連續變化，為了降低跳變對參數估計的影響，對提取的多普勒曲線進行平滑濾波處理，得到的曲線如圖9(b)所示。為了進行振動角的估計，還需要對曲線進行去均值處理，處理結果如圖9(c)所示。通過測量多普勒曲線波峰之間的時間差即可估計出振動頻率，即。提取的散射點1, 2, 3的距離變化曲線如。通過式(16)可以得到初相的估計值。最后對振動角和初相參數進行精搜索，機翼回波熵值優化函數的代價曲面如圖11所示，可以發現，熵值是平滑變化的，且是凸的，所以使用本文所提搜索策略能夠快速得到參數的精確估計結果。最終振動參數的估計結果如表2所示。機翼區域精細化處理前后的對比圖像如圖12所示。

圖9 提取的不同散射點的多普勒變化曲線Fig.9 Extracted Doppler curves of different scatterers

圖10 提取的不同散射點的距離變化曲線Fig.10 Extracted range curves of different scatterers

圖11 圖像熵值優化函數代價曲面Fig.11 Cost surface of entropy optimization function

表2 仿真數據機翼振動參數估計結果Tab.2 Vibration parameters estimation result of simulation data

為了進一步驗證本文方法的有效性，下面對實測數據進行處理，實測機身成像定標結果如圖13所示，通過測量機身軸線與Y軸之間的夾角，可以得到方位角的估計，即。觀測飛機的機型為A320，通過公開資料可知其機身寬度為3.95 m，通過測量定標圖像中機身寬度，可以得到俯仰角的估計。圖14是振動較為明顯的機翼不同子孔徑的距離多普勒成像結果，通過提取圖中用圓圈標記的散射點圖10所示，通過在距離變化曲線上選取兩個時間點，聯立式(14)和式(15)求解得到振動角的估計值1的距離和多普勒變化，如圖15所示。最小熵優化函數的代價平面如圖16所示，振動參數估計結果如表3所示。機翼區域精細化處理前后的對比圖像如圖17所示，從圖中用圓圈標記的散射點的聚焦情況可以看出，使用本文方法對機翼區域精細化處理之后，散射點的聚焦質量有所改善。

圖12 機翼區域精細化處理前后結果對比Fig.12 Before and after accurate processing of airplane wings

圖13 實測數據機身聚焦成像定標結果Fig.13 Focused and scaled airplane body result of the measured data

圖15 提取的散射點的距離和多普勒變化曲線Fig.15 Range and Doppler curves of extracted scatter

圖16 圖像熵值優化函數代價曲面Fig.16 Cost surface of image entropy optimization function

表3 實測數據機翼振動參數估計結果Tab.3 Vibration parameters estimation result of measured data

4 結束語

本文針對機翼在飛機非平穩運動的狀態下振動的現象，提出一種利用微波光子超高分辨雷達估計機翼振動參數的方法。首先通過對平穩運動的機身進行聚焦成像和定標處理，從中估計雷達視線角。然后通過對機翼振動進行建模，推導了機翼回波的解析式，并分析了振動對成像的影響。針對振動導致的機翼成像結果散焦的問題，本文首先通過子孔徑序列成像提取散射點距離和多普勒曲線，再聯合雷達視線角參數進行振動參數的初步估計，然后提出了修正的極坐標格式算法，構造了最小熵優化函數對振動參數進行精確估計。利用估計的振動參數，修正的極坐標格式算法能夠對機翼進行聚焦成像。本文對仿真數據和飛機實測數據進行了試驗，試驗結果驗證了本文方法的有效性。

圖17 機翼區域精細化處理前后結果對比Fig.17 Before and after accurate processing of airplane wings