符吉祥 邢孟道* 徐 丹 王安樂(西安電子科技大學(xué)雷達(dá)信號(hào)處理國家重點(diǎn)實(shí)驗(yàn)室 西安 710071)

②(西安電子科技大學(xué)信息感知技術(shù)協(xié)同創(chuàng)新中心 西安 710071)

③(空軍預(yù)警學(xué)院 武漢 430019)

1 引言

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

圖1 微波光子實(shí)測飛機(jī)目標(biāo)1維距離像Fig.1 One dimensional range profile of airplane measured by microwave photonic radar

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

2 機(jī)翼振動(dòng)建模分析和振動(dòng)參數(shù)估計(jì)

2.1 運(yùn)動(dòng)模型分析

根據(jù)力學(xué)分析可知，機(jī)翼的振動(dòng)與單擺運(yùn)動(dòng)存在一致性，所以將其建模為單擺運(yùn)動(dòng)。如圖2所示，以飛機(jī)整體所在平面作為XOY平面，X軸指向?yàn)樵趨⒖紩r(shí)刻沿著機(jī)身方向指向機(jī)頭的方向，Z軸為飛機(jī)平面的法線方向。雷達(dá)視線與Z軸的夾角為φ0，其在XOY平面上的投影向量與Y軸之間的夾角為β0。機(jī)翼沿著Z軸振動(dòng)，振動(dòng)過程中，機(jī)翼與XOY平面之間的夾角為?，最大夾角為?0，根據(jù)單擺轉(zhuǎn)角公式可知，夾角隨慢時(shí)間正弦變化，即

其中，?0表示振動(dòng)角，fv表示振動(dòng)頻率，φ0表示振動(dòng)初相，表示方位慢時(shí)間。在飛機(jī)沿著圖2中曲線軌跡運(yùn)動(dòng)的時(shí)候，飛機(jī)整體還有相對(duì)于雷達(dá)的平動(dòng)和轉(zhuǎn)動(dòng)。對(duì)于機(jī)身上的散射點(diǎn)p1，因?yàn)槠洳淮嬖谡駝?dòng)，所以其相對(duì)于雷達(dá)的運(yùn)動(dòng)可以表示為

其中，R0表示飛機(jī)旋轉(zhuǎn)中心O到雷達(dá)的初始距離，vr和ar分別表示飛機(jī)相對(duì)于雷達(dá)運(yùn)動(dòng)的徑向速度和徑向加速度。轉(zhuǎn)動(dòng)分量是隨散射點(diǎn)位置發(fā)生空變的，其形式為(平面波近似)

其中，(xp1，yp1)表示散射點(diǎn)p1參考時(shí)刻的坐標(biāo)位置，θ表示飛機(jī)相對(duì)于雷達(dá)轉(zhuǎn)動(dòng)的角度，其是慢時(shí)間的函數(shù)，即θ(tm)=ωtm,ω表示轉(zhuǎn)速。對(duì)于機(jī)翼上的散射點(diǎn)，其不僅包含以上的運(yùn)動(dòng)，還包含機(jī)翼振動(dòng)對(duì)其的運(yùn)動(dòng)調(diào)制。通過運(yùn)動(dòng)分解可以得到機(jī)翼上散射點(diǎn)p2各個(gè)時(shí)刻的坐標(biāo)位置如式(5)

圖2 飛機(jī)運(yùn)動(dòng)幾何模型Fig.2 Geometry model of airplane

其中，(xp2，yp2)表示散射點(diǎn)p2參考時(shí)刻的坐標(biāo)位置。所以散射點(diǎn)p2到雷達(dá)的斜距變化為

假設(shè)雷達(dá)發(fā)射線性調(diào)頻信號(hào)，使用去斜方式接收回波，經(jīng)過剩余視頻項(xiàng)補(bǔ)償之后，機(jī)翼散射點(diǎn)雷達(dá)回波信號(hào)的波數(shù)譜表達(dá)式為[10]

2.2 振動(dòng)參數(shù)估計(jì)方法

首先通過粗成像對(duì)機(jī)翼機(jī)身回波進(jìn)行分離。因?yàn)闄C(jī)翼存在周期性振動(dòng)而機(jī)身只存在平穩(wěn)的運(yùn)動(dòng)，所以在粗成像結(jié)果中，機(jī)翼和機(jī)身上的散射點(diǎn)散焦情況存在明顯的差異，即機(jī)翼散射點(diǎn)散焦成線性分布而機(jī)身散射點(diǎn)散焦呈現(xiàn)塊狀分布。另外，從粗成像的輪廓中也可以輕易地將機(jī)翼和機(jī)身進(jìn)行區(qū)分。如圖4所示，圖4(a)是飛機(jī)整體粗成像結(jié)果，可以發(fā)現(xiàn)，盡管飛機(jī)散焦嚴(yán)重，但是仍能夠從輪廓中輕易區(qū)分機(jī)身和機(jī)翼，圖4(b)是機(jī)身散射點(diǎn)塊狀散焦結(jié)果，圖4(c)是機(jī)翼散射點(diǎn)線性散焦結(jié)果。最后通過加窗濾波將成像結(jié)果中的機(jī)身回波濾出，從而實(shí)現(xiàn)機(jī)身機(jī)翼回波分離。

對(duì)機(jī)身回波使用keystone算法[13]校正其線性距離走動(dòng)，得到

其中，sinc(x)=sin(x)/(x),表示徑向波數(shù)中心。觀察式(9)可以發(fā)現(xiàn)，keystone變換之后，包絡(luò)已經(jīng)對(duì)齊，然而散射點(diǎn)回波存在和其距離坐標(biāo)yp線性相關(guān)的方位向二次相位。通過對(duì)不同距離單元信號(hào)進(jìn)行調(diào)頻率估計(jì)并對(duì)調(diào)頻率進(jìn)行線性擬合即可得到轉(zhuǎn)速ω，然后再構(gòu)造調(diào)頻率補(bǔ)償函數(shù)

圖3 轉(zhuǎn)臺(tái)模型Fig.3 Turntable model

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

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

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

圖4 飛機(jī)粗成像結(jié)果Fig.4 Coarse imaging results of airplane

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

觀察機(jī)翼散射點(diǎn)回波，發(fā)現(xiàn)其距離和方位是解耦的，利用極坐標(biāo)格式算法的思想，本文提出一種修正的極坐標(biāo)格式算法，即

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

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

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

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

3 仿真實(shí)驗(yàn)及實(shí)測數(shù)據(jù)驗(yàn)證

本文使用Ka波段雷達(dá)參數(shù)對(duì)飛機(jī)散射點(diǎn)模型進(jìn)行參數(shù)估計(jì)，為了簡化仿真，將飛機(jī)兩側(cè)機(jī)翼振動(dòng)參數(shù)設(shè)置為相同的，仿真參數(shù)如表1所示。飛機(jī)模型長35 m，寬30 m，如圖6所示。雷達(dá)視線方位角為80°，俯仰角為85°。

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

表1 仿真雷達(dá)參數(shù)Tab.1 Simulation radar parameters

圖6 飛機(jī)散射點(diǎn)模型Fig.6 Scatter model of airplane

圖7 機(jī)身聚焦成像定標(biāo)結(jié)果Fig.7 Focused and scaled result of airplane body

圖8 不同時(shí)刻機(jī)翼子孔徑成像結(jié)果Fig.8 Sub?aperture imaging results of airplane wing at different time

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

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

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

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

圖11 圖像熵值優(yōu)化函數(shù)代價(jià)曲面Fig.11 Cost surface of entropy optimization function

表2 仿真數(shù)據(jù)機(jī)翼振動(dòng)參數(shù)估計(jì)結(jié)果Tab.2 Vibration parameters estimation result of simulation data

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

圖12 機(jī)翼區(qū)域精細(xì)化處理前后結(jié)果對(duì)比Fig.12 Before and after accurate processing of airplane wings

圖13 實(shí)測數(shù)據(jù)機(jī)身聚焦成像定標(biāo)結(jié)果Fig.13 Focused and scaled airplane body result of the measured data

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

圖16 圖像熵值優(yōu)化函數(shù)代價(jià)曲面Fig.16 Cost surface of image entropy optimization function

表3 實(shí)測數(shù)據(jù)機(jī)翼振動(dòng)參數(shù)估計(jì)結(jié)果Tab.3 Vibration parameters estimation result of measured data

4 結(jié)束語

本文針對(duì)機(jī)翼在飛機(jī)非平穩(wěn)運(yùn)動(dòng)的狀態(tài)下振動(dòng)的現(xiàn)象，提出一種利用微波光子超高分辨雷達(dá)估計(jì)機(jī)翼振動(dòng)參數(shù)的方法。首先通過對(duì)平穩(wěn)運(yùn)動(dòng)的機(jī)身進(jìn)行聚焦成像和定標(biāo)處理，從中估計(jì)雷達(dá)視線角。然后通過對(duì)機(jī)翼振動(dòng)進(jìn)行建模，推導(dǎo)了機(jī)翼回波的解析式，并分析了振動(dòng)對(duì)成像的影響。針對(duì)振動(dòng)導(dǎo)致的機(jī)翼成像結(jié)果散焦的問題，本文首先通過子孔徑序列成像提取散射點(diǎn)距離和多普勒曲線，再聯(lián)合雷達(dá)視線角參數(shù)進(jìn)行振動(dòng)參數(shù)的初步估計(jì)，然后提出了修正的極坐標(biāo)格式算法，構(gòu)造了最小熵優(yōu)化函數(shù)對(duì)振動(dòng)參數(shù)進(jìn)行精確估計(jì)。利用估計(jì)的振動(dòng)參數(shù)，修正的極坐標(biāo)格式算法能夠?qū)C(jī)翼進(jìn)行聚焦成像。本文對(duì)仿真數(shù)據(jù)和飛機(jī)實(shí)測數(shù)據(jù)進(jìn)行了試驗(yàn)，試驗(yàn)結(jié)果驗(yàn)證了本文方法的有效性。

圖17 機(jī)翼區(qū)域精細(xì)化處理前后結(jié)果對(duì)比Fig.17 Before and after accurate processing of airplane wings