基于虛擬慢時間的雙基地ISAR成像算法

2019-05-25 02:12:04史林韓寧宋祥君王立兵崔東輝

航空學報 2019年5期

史林，韓寧，宋祥君，王立兵，崔東輝

1. 陸軍工程大學石家莊校區電子與光學工程系，石家莊 050003 2. 中國人民解放軍32181部隊，石家莊 050003 3. 中國人民解放軍63961部隊，北京 100010 4. 中國人民解放軍78616部隊，成都 610214

雙基地逆合成孔徑雷達(Inverse Synthetic Aperture Radar, ISAR)是空間目標監視的一種重要手段，可實現目標的多姿態觀測，當發射站靠前布置時還可增大成像系統作用距離。雙基地ISAR具有良好的“四抗性能”[1]，且成像不受目標運動方向限制，已逐漸成為國內外的研究熱點[2-8]。因雙基地配置固有的“三大同步”問題，目前針對雙基地ISAR的研究還處于較為開放的研究階段。國內外對雙基地ISAR的研究也主要集中在運動補償[9-10]、成像平面確定[11-13]、三維干涉成像[14-16]以及方位定標[17-18]等方面。

因雙基地雷達系統固有的“三大同步”問題，反向投影(Back Projection, BP)、極坐標格式算法(Polar Format Algorithm, PFA)等對雙基地雷達的系統配置精度、目標坐標位置和同步時鐘誤差敏感的算法難以在雙基地ISAR成像中應用。距離-多普勒(Range Doppler, RD)算法物理意義明確、對收發同步精度的要求較低，被廣泛應用于雙基地ISAR成像仿真及實測數據處理。目前，大部分針對雙基地ISAR的研究都假定成像期間雙基地角是恒定不變的，但在實際成像過程中，雙基地角是個時變量，雙基地角時變會造成二維ISAR成像在方位向的散焦，影響圖像的聚焦度。針對此問題，文獻[19]利用粒子群優化算法估計雙基地角時變帶來的高次項系數，然后構造補償相位項，實現自聚焦，但該算法中使用的粒子群優化算法是一種智能優化算法，只能以一定的概率估計出最優補償系數，算法無法保證每次都收斂到最優解。文獻[20]將稀疏分解算法引入雙基地ISAR成像的自聚焦過程中，基于冗余基的高分辨特性，估計出高次項的系數而后完成補償。為了提高補償系數的精度，需擴大冗余基的個數，使得該算法的運算量較大，且該算法的最終成像效果受算法正則參數的影響較大，正則參數選取不當可能會造成無法對目標進行成像。文獻[21]研究了雙基地ISAR成像系統中越分辨單元徙動帶來的影響及相應的處理算法。該算法基于廣義Keystone變換，消除越距離單元走動；基于最大圖像對比度準則，估計等效旋轉中心，進行越多普勒單元徙動校正；依據圖像畸變角度，通過對距離單元內像素進行移位操作實現圖像的畸變校正。文獻[22]基于圖像旋轉相關度最大準則，提出了一種更為精確的雙基地ISAR等效旋轉中心估計算法。其采用與文獻[21]相同的越多普勒單元徙動校正和圖像畸變校正算法。文獻[21-22]中均假定完全精確已知目標空間位置和對應的雙基地角。實際系統中，此假定并不成立。利用先驗信息得到的衛星位置信息，是有誤差的。此誤差會導致每個周期雙基地角和累積轉角均產生誤差。若對應的誤差值較大，會影響后續越多普勒單元徙動校正和圖像畸變校正。需進一步考慮雙基地角、累計轉角存在誤差時相應的處理算法。

針對以上問題，本文聚焦于雙基地角時變下，實際成像系統中多普勒向散焦及圖像畸變問題，提出一種基于虛擬慢時間采樣消除高次項影響，而后通過非均勻傅里葉變換完成方位壓縮得到目標二維ISAR像的新型成像算法。

1 雙基地ISAR成像模型

圖1 雙基地ISAR成像原理模型Fig.1 Imaging principle model for bisatic ISAR

本文以平穩空間目標為研究對象，其成像原理模型如圖1所示。圖中：T為發射站雷達；R為接收站雷達；TR為基線；L為長度為；Rt0、Rr0分別為觀測起始時刻，目標相位中心距發射站和接收站的距離；Rt1、Rr1分別為目標運動到下一時刻，目標相位中心距發射站和接收站的距離；V為目標運動方向矢量；C為目標上的任一散射點；在觀測起始時刻，散射點C在收發雙站雷達和目標質心O確定的平面內的投影，記為E；RtC0、RrC0分別為散射點C在觀測起始時刻，距發射站、接收站雷達的距離。

在觀測起始時刻，以目標質心O為原點，建立慣性坐標系：在觀測起始時刻，y軸正方向定義為目標雙基地角平分線延長線方向；在目標雙基地角平分線與目標軌道曲線構成的平面內，將y軸的法線定義為x軸，x軸正方向定義為目標運動方向。該坐標軸指向，不隨目標的運動變化。此坐標系下，目標散射點距離的變化可分解為平動和相對轉動兩部分。O′由目標質心O移動得到。為便于分析目標轉動情況，以O′為原點，建立以雙基地角平分線延長線方向為y′軸的坐標系x′O′y′，該坐標系的y′軸隨著雙基地角平分線指向變化而改變；x′軸則是此刻目標雙基地角平分線與目標軌道曲線構成的平面內y′軸的法向。因此，在觀測時間內x′O′y′坐標系與xOy坐標系間的相對轉動角反映了目標的相對轉動情況。假定雷達發射的線性調頻信號為

(1)

假雙基地雷達理想同步，且成像期間雙基地角恒定不變，在中頻采樣后通過數字下變頻得到基頻信號為

(2)

式中：σC為散射點C的散射系數；c為真空中的光速；R(tm)可表示為

RC(tm)=Rref(tm)+Rrot(tm)≈

(3)

式中：Rref(tm)為目標散射中心的平動分量；Rrot(tm)為散射點C的轉動分量；θC為散射點C的方位矢量與xOy坐標系中x軸正向的夾角；ψ(tm)為成像期間雙基地角平分線的轉動角度；β為成像期間的雙基地角。

結合式(2)和式(3)可以看出，雙基地ISAR回波中散射點到收發雙站的距離依然可以分解為平動項與轉動項兩部分。對式(2)進行脈沖壓縮利用數字脈壓完成距離維成像，而后進行運動補償，在包絡對齊之后，將平動和轉動導致的相位項統一建模，進行相位補償實現自聚焦，通過方位維壓縮得到二維ISAR像[4,18]。

2 雙基地角時變對ISAR成像的影響機理

在前文推導式(3)的過程中，假定了成像期間雙基地角恒為β，當雙基地角隨慢時間變化時，式(3)可改寫為

(4)

式中：β(tm)表示成像期間雙基地角隨慢時間tm變化。在較短的相干處理時間(Coherent Processing Interval，CPI)內，目標處于遠場位置的條件下，雙基地角隨慢時間近似成線性變化，可以用一階泰勒展開進行近似[4]

β(tm)=β0+Δβtm

(5)

式中：β0為零時刻的雙基地角；Δβ為雙基地角在零時刻的一階導數。

將式(5)代入cos(β(tm)/2)進行展開，并忽略二次以上的高次項得到

K0+K1tm

(6)

式中：

(7)

(8)

在下文中，將K0、K1分別稱為虛擬慢時間系數K0和虛擬慢時間系數K1。實際的雙基地雷達系統中，收發雙站的位置是固定的，可以根據收發雙站的位置信息和目標的軌道信息，依據幾何位置關系，獲得雙基地角信息，進而獲得虛擬慢時間系數的值。

在較短的CPI時間內，平穩目標旋轉一個微小的角度ψ(tm)=ωtm，ω為雙基地角的平均轉動角速度,此時sinψ(tm)和cosψ(tm)可近似為sinψ(tm)≈ωtm，cosψ(tm)≈1。對目標回波作平動補償后，目標平動項被去除，目標依然可等效為轉臺目標。將式(6)代入式(4)中的目標轉動項進行化簡可得

ΔRC(tm)≈ 2(xiωtm+yi)(K0+K1tm)=

(9)

式中：(xi,yi)為散射點C在xOy坐標系中的坐標。

經過理想的運動補償(包絡對齊和初相校正)后，脈沖壓縮后散射點的回波式(2)可表示為

(10)

將式(9)代入PC，可得

(11)

忽略常數項2yiK0，則第n個距離單元的對應的相位多項式信號為

(12)

式中：Ln為第n個距離單元內散射點的個數；i為第n個距離單元內散射點對應的下標。位于第n個距離單元內散射點的Ln個散射點的縱向距離(y1=y2=…=yLn)相等，因此，式(12)可以重寫為

(13)

從式(13)可以看出，因成像期間雙基地角隨慢時間的變化，式中第1項將會引起圖像畸變，且畸變量與距離坐標成正比，式中第2項導致成像所需的轉動相位項出現了二次高階項，若不對其進行補償則會造成散射點方位向的散焦，進而會影響圖像的聚焦度。

3 虛擬慢時間成像算法

3.1 初次相位補償

為了消除圖像畸變，同時為消除式(13)中高次項的影響，首先構建初次補償相位項:

(14)

yn=(n-nc)Δy

(15)

式中：nc為等效旋轉中心對應的縱向距離下標。因此，為構造初次相位補償項φ1，需要確定散射點相對目標等效旋轉中心的縱向定標距離yn、虛擬慢時間系數K0和虛擬慢時間系數K1。其中，虛擬慢時間系數可以通過雙基地角得到。散射點縱向距離的定標量需要確定圖像的等效旋轉中心的縱向位置。等效旋轉中心的位置的估計在3.4節進行分析。

已獲得初次補償相位φ1的基礎上，將式(14)與式(13)相乘，完成初次相位補償后可得

(16)

3.2 非均勻傅里葉變換實現方位壓縮

完成初次相位補償后，由于雙基地角時變的影響，回波數據中包含與散射點方位向坐標有關的高次項。此時基于傅里葉變換進行方位壓縮，將引起方位向散焦。為消除高次項的影響，可采用距離瞬時多普勒(Range Instantaneous Doppler, RID)成像算法，來提高成像質量。RID成像算法可分為兩大類。第一類是基于時頻分析的成像算法。基于Wigner-Ville分布(Wigner-Ville Distribution, WVD)類算法屬于雙線性變換，需要在時頻聚集性和交叉項之間進行平衡。基于參數化的時頻分布算法，如基于自適應線性調頻分解類的算法，無交叉項影響，但運算量大，且對噪聲和初值選擇比較敏感，應用受限。另外一類是基于參數估計的高次相位補償算法。該類算法通過分數階傅里葉變換(Fractional Fourier Transform, FrFT)等參數估計算法，估計信號的高階項系數，補償高次相位的影響。但此類算法需要估計每個距離單元相位多項式信號的參數，運算量大，且成像效果依賴與參數估計的精度。空間目標脈沖壓縮后回波信噪比低，參數估計精度不高，此類算法應用受限。

雙基地ISAR成像期間，基于空間目標軌道的先驗信息和成像幾何關系，可獲得相對精確的雙基地角信息，進一步估計相應的時變信息。定義虛擬慢時間τm滿足如下關系式:

(17)

定義式(17)為虛擬慢時間。將式(7)和式(8)代入式(17)可得每個慢時間tm對應的虛擬慢時間τm。將式(17)代入式(16)可得

(18)

圖2 虛擬慢時間采樣前后的數據平面圖Fig.2 Data plane before and after virtual slow time sampling

式(18)中相位項僅包含虛擬慢時間τm的一次項。利用式(17)對回波數據完成虛擬慢時間映射前后的數據平面如圖2所示。圖2(a)和圖2(b)的橫軸分別為慢時間和虛擬慢時間，縱軸均為快時間。對比圖2(a)和圖2(b)，回波數據完成虛擬慢時間映射后，相鄰虛擬慢時間之間不再是固定不變的時間間隔。經過虛擬慢時間映射，去掉了慢時間二次項影響，等效對方位向回波實現了虛擬非均勻采樣，可利用非均勻傅里葉變換完成方位壓縮得到目標的二維ISAR像。

令τ=[τ0,τ1,…,τM-1]，則補償系數矩陣的基可表示為

(19)

基于式(19)，構建用于方位壓縮的補償系數矩陣Ψ:

Ψ=[γ0,γ1,…,γM-1]

(20)

利用式(20)對初次相位補償后的數據進行方位壓縮，即可得到目標的二維ISAR像。由于非均勻傅里葉變換采用的不是完全正交的基函數，無法采用Cooley和Tukey算法[23]。若直接進行矩陣運算，計算虛擬非均勻采樣信號sn(τm)與每一個頻率樣本點的相關值，需要4N次的實數加法和(4N-2)次的矩陣乘法。與N個頻率點進行相關運算，總的運算量為O(N2)數量級。為了減少運算量，采用文獻[24]中的Goertzel算法。這種算法總的計算量是(2N+4)次實數相乘和(4N-2)次實數相加，將運算復雜度降低到O(N)的數量級。若對W個距離單元進行方位壓縮，計算量為O(WN)。

3.3 等效旋轉中心估計

(21)

式中：C_imag為圖像對比度；I(x,y)為復圖像的幅度；A(I(x,y))為圖像在整個成像平面上的幅度平均。首先假定等效旋轉中心位置，并以此中心構造補償相位項φ1，進行初次相位補償后，進行虛擬慢時間采樣、方位向非均勻方位壓縮，得到ISAR二維像，計算圖像對比度。然后，更換假定的等效旋轉中心位置，重復以上步驟。當等效旋轉中心位置是實際旋轉中心時，圖像對比度最大，圖像聚焦效果最好。為了減小搜索范圍，考慮到等效旋轉中心位置一般與強散射點有關，因此，選擇峰值較大的距離單元兩側進行搜索，既可以快速找到等效旋轉中心，又減少了運算量。

3.4 虛擬慢時間系數參數估計

前文的分析中，假定根據目標的軌道信息和收發雙站的位置的先驗信息，精確獲得雙基地角的信息。實際的系統中，雙基地角的估計值和真實值之間存在著一定的誤差。因此，本節分析在雙基地角存在估計誤差的情況下的虛擬慢時間系數K0和K1參數估計算法和魯棒性。

假定系統測量的雙基地角為

(22)

(23)

(24)

3.5 算法流程

綜合以上分析，本文提出的虛擬慢時間ISAR成像算法流程如圖3所示。

圖3 所提算法的成像流程Fig.3 Imaging process of proposed algorithm

在圖3中用虛線框內內容為本文所提算法的關鍵步驟。具體步驟為

步驟1對雙基地ISAR回波進行脈沖壓縮、包絡對齊、相位校正，得到平動校正后的一維距離像序列。

步驟3假定等效旋轉中心位置，按照式(14) 和式(15)構造初次相位補償項，補償一維距離像。

步驟4按照式(17)進行虛擬慢時間采樣，并完成非均勻傅里葉變換，得到ISAR圖像，計算圖像對比度。

步驟5假定新的等效旋轉中心位置，重復步驟3和步驟4并將此次圖像的對比度與上次的對比，若對比度變大，存儲二維ISAR圖像，如此循環，直到遍歷完畢可能的等效旋轉中心位置。

步驟6遍歷結束，存儲的即為圖像對比度最大的二維ISAR圖像，輸出圖像矩陣。

4 仿真實驗

4.1 典型散射點目標成像結果分析

仿真雙基地雷達參數及目標的初始軌道根數如表1和表2所示。積累脈沖數為512。回波基于文獻[25]的算法模擬生成。通過收發雙站及目標的位置信息和成像幾何關系，獲得雙基地角和等效累積轉角。

表1 仿真雙基地雷達參數Table 1 Simulation parameters of bistatic radar

表2 初始軌道根數Table 2 Initial two line elements

圖4給出了仿真用的目標三維和二維散射點模型。雙基地角是雙基地ISAR成像的重要參數，為此，給出了某一可見觀測時間段內的雙基地角變化曲線，如圖5所示，并由此確定所選擇的仿真成像弧段。從圖中可以看出，在觀測時間內，只有在曲線頂點附近的小部分時間內雙基地角是非線性變化的(在基線的中心位置附近)，大多數時間內雙基地角是線性變化的。選取特定的512個周期回波數據作為成像段數據，如圖5中加粗線段所示，對應的積累時間為5.12 s。該觀測時間內雙基地角變化約為6.36°，放大后可看出，此期間雙基地角與慢時間近似成線性關系。等效單基地雷達對應的成像累積轉角隨積累脈沖數的變化情況如圖6所示。期間目標等效累積轉角約為4.26°，并且累積轉角與慢時間成線性關系。

圖4 三維和二維散射點模型Fig.4 Three-dimensinal and two-dimensinal scattering point model

圖5 雙基地角隨觀測時間變化曲線Fig.5 Variation curve of bistatic angles with observation time

圖6 目標等效累積轉角隨脈沖個數的變化曲線Fig.6 Variation curve of equivalent cumulative angle of target with pulse number

圖7(a)和圖7(b)為RD成像算法得到的二維ISAR像(基于最大互相關法完成包絡對齊；利用相位梯度自聚焦(Phase Gradient Autofocusing，PGA)算法完成初相校正)。圖7(a)為定標前的ISAR像，圖7(b)為定標后的ISAR像。由于雙基地時變角的影響，圖7(a)和圖7(b)中的二維ISAR像是“歪斜”的，圖像產生了畸變，并在方位向上存在散焦現象。

依然基于最大互相關法完成包絡對齊；利用PGA算法完成初相校正，完成平動運動補償。采用基于圖像對比度最大的搜索方法搜索等效旋轉中心。等效旋轉中心估計曲線如圖8所示，圖像對比度最大時對應第500個距離采樣單元。以等效旋轉中心位置為中心點，進行初次相位補償，并進行虛擬慢時間采樣，構造補償系數矩陣，完成非均勻傅里葉變換，得到ISAR二維圖像如圖7(c)和圖7(d)所示。圖7(c)為定標前的ISAR像，圖7(d) 為定標后的ISAR像。從圖7(c)和圖7(d) 可以看出，利用所提算法可以正確生成目標的二維ISAR像，驗證了所提算法的有效性。相比于圖7(a)和圖7(b)，圖7(c)和圖7(d)利用所提成像算法生成的ISAR像其聚焦度優于RD成像算法生成的ISAR像。從圖7(c)和圖7(d)可以看出，所提算法通過初次相位補償，亦可以有效校正“歪斜”項，消除圖像畸變。圖7(d)中定標后的圖像形狀與散射點模型一致，可等效為原圖形圍繞等效選擇中心旋轉了一定角度得到瞬時像，有利于后期目標的正確識別。

為了定量分析圖像聚集度的變化程度，分別計算圖7(b)和圖7(d) 定標后圖像的對比度和方位向3 dB寬度均值。圖7(b)和圖7(d)定標后圖像對比度分別為16.28和24.05，方位向3 dB寬度均值分別為0.331 m和0.229 9 m。據此可以看出，基于所提成像算法獲得的ISAR成像，圖像的對比度有明顯提升。

圖9為采用與圖7(c)和圖7(d)采用相同的等效旋轉中心估計和初次相位補償操作后，然后基于偽WVD(Pseudo-Wigner-Ville Distribution, PWVD)時頻分析算法得到t=2.5 s時ISAR圖像。圖9(a)是未定標的ISAR圖像，圖9(b)表示定標后的ISAR圖像。雖然基于PWVD得到圖像分辨率更高，但對于存在兩個散射點的距離單元，得到的二維像中出現了交叉項，會影響目標識別。若距離單元中存在更多散射點，則會出現更為嚴重的交叉項干擾，影響目標識別。在此方面，相比基于PWVD的算法，所提算法具有優勢。

圖7 采用RD算法和所提成像算法的ISAR成像結果Fig.7 ISAR imaging results using RD and proposed algorithms

圖8 等效旋轉中心估計曲線Fig.8 Curve of equivalent rotation center estimation

圖9 RID算法成像結果(t=2.5 s)Fig.9 Imaging results using RID algorithm (t=2.5 s)

4.2 魯棒性分析

假定目標距離信息含有[-5 m, 5 m]均勻分布的隨機誤差，與文獻[21]進行了仿真對比驗證。文獻[21]中的角度誤差累積分布如圖10所示。雙基地角誤差近似在[-0.018°, 0.018°]內均勻分布，累積轉角誤差近似在[-0.003 2°, 0.003 2°]內均勻分布，雙基地角誤差近在[-20°, 110°]內不規則分布。

圖10 各類角的誤差累積分布Fig.10 Cumulative distribution of errors at different angles

在此誤差條件下，所提算法能有效成像，獲得目標的正確形狀，消除方位向散焦(見圖11(c))。采樣文獻[21]中的算法，部分距離單元方位向無法完成有效壓縮(見圖11(a))，且進行歪斜校正后無法獲得目標真實形狀(見圖11(b))。這是由于在實際系統中，完全精確已知目標空間位置的假定是不成立的。包含誤差的先驗信息，會導致每個周期雙基地角和累積轉角均產生誤差(見圖10(a)和10(b))，且較小的測距誤差會導致很大圖像畸變角誤差，(見圖10(c))。這些誤差影響了后續越分辨單元徙動校正和圖像畸變校正。本文通過LSE算法，估計虛擬慢時間系數。在每一個雙基地角均存在隨機誤差的情況下(見圖10(a)和10(b))，基于均方誤差最小的約束，亦可以精確的估計出相應系數，可以有效成像(見圖11(c))，所提算法更為魯棒。

圖11 誤差條件下的成像結果Fig.11 Imaging results under error conditions

為了分析雙基地角誤差對成像算法性能的影響，將式(22)重寫為

(25)

在4.1節選用的成像段中，K0和K1的值分別為0.765 6和-0.007 5，其絕對值相對較小。但對于ISAR成像，需要考察目標空間位置的變化與波長的比擬程度，誤差相位的數量級要求一般與λ/8相比擬。可以通過分析其相對變化值，來定量分析參數的誤差及對成像的影響。因此，進一步定義參數的相對誤差：

(26)

假定測量誤差εβ(tm)服從[-Δε/2,Δε/2]的均勻分布，Δε/2代表εβ(tm)誤差的最大值。為了驗證Δε對參數估計的影響，將Δε從0°～4°，依照步長0.2°步進，進行500次Monte Carlo仿真。圖12中，Δε的值較大時，εK0和εK1的均值和方差依然很小。在Δε=4°時，εK0的均值和方差分別為0.052 9%和0.039 1%，εK1的均值和方差分別為0.304 8%和0.226 3%。在存在一定量的雙基地角測量誤差時，通過LSE估計，能獲得高精度的K0和K1的估計值。

在Δε=4°時，計算圖像對比度，仿真500次，表3給出用RD算法和所提算法得到的圖像對比度的平均值。表3還給出了散射點的距離向3 dB主瓣寬度和方位向3 dB主瓣寬度平均值。

觀察表3可以看出，所提算法得到的圖像對比度高于RD成像算法生成的圖像對比度，距離向3 dB主瓣寬度在兩種成像算法下并無明顯變化，方位向3 dB主瓣寬度的聚焦改善明顯。這與本文所提算法可提高方位向聚焦度的理論分析一

圖12 Δε對K0和K1估計精度的影響Fig.12 Effects of Δε on estimation accuracy of K0 and K1

致，證明了所提算法的有效性和魯棒性。需要說明的是，距離壓縮和方位壓縮時使用了Hamming窗。由于存在誤差，表3中距離向和方位向3 dB寬度較理論寬度值有所增大(加Hamming窗后，3 dB距離向和方位向寬度的理論值為相應分辨率乘以系數1.3，距離向和方位向分辨率理論值分別為0.200 1×1.3=0.261 m，0.176 8×1.3=0.229 8 m)[26]。