基于ADMM 算法的機載太赫茲SAR成像高頻振動補償

郝昭昕 張啟明 孫進平 王彥平

（1.北京航空航天大學電子信息工程學院，北京 100191；2.北方工業(yè)大學信息學院，北京 100144）

1 引言

太赫茲波的頻率為0.1～100 THz，具有頻率高、波長短的特點，因此與傳統(tǒng)的微波合成孔徑雷達（Synthetic Aperture Radar，SAR）相比，太赫茲SAR 具有高成像分辨率、良好的穿透性和高幀頻成像能力等優(yōu)異性能，擅長對低速和微動目標進行檢測和識別［1-3］。因此隨著太赫茲雷達系統(tǒng)的逐步成熟，太赫茲SAR 在無損檢測、視頻SAR、安防安檢等領(lǐng)域均得到了廣泛的應(yīng)用［4-8］。但與此同時，太赫茲SAR 也存在對平臺的高頻振動十分敏感的問題。在實際的機載太赫茲SAR 成像過程中，由于氣流以及載機自身飛行特性的影響，載機的飛行軌跡并不穩(wěn)定，自身會存在高頻振動，使成像結(jié)果出現(xiàn)散焦現(xiàn)象和鬼影目標。因此太赫茲SAR 成像對運動補償算法提出了更高的要求。載機的非理想運動可以分為高頻振動和低頻運動兩部分。其中低頻運動分量造成的相位誤差可以通過傳統(tǒng)的運動補償方法進行補償，如傳感器測量數(shù)據(jù)補償以及自聚焦［9-10］。由于現(xiàn)有的傳感器測量精度不足，無法準確測量高頻振動的數(shù)據(jù)，因此傳統(tǒng)的運動補償算法難以補償高頻振動誤差，必須采用先進的信號處理方法估計和補償太赫茲SAR 成像中的高頻誤差。

目前對于太赫茲SAR 高頻振動相位誤差的估計與補償已取得了一定的研究成果。平臺的高頻振動可以建模為多分量簡諧運動之和［11］，一些研究人員因此提出了基于單通道回波數(shù)據(jù)的振動相位誤差補償算法。其中文獻［12-14］提出使用離散正弦調(diào)頻變換并結(jié)合智能優(yōu)化算法進行振動參數(shù)估計，文獻［15］使用正弦調(diào)頻傅里葉變換估計振動參數(shù)，文獻［16］提出使用分數(shù)階傅里葉變換并結(jié)合準最大似然法和隨機采樣一致性原理估計信號的瞬時調(diào)頻率，然后進一步估計高頻振動參數(shù)。這些方法一般假定高頻振動包含2 個簡諧運動分量，當振動的簡諧運動分量數(shù)增加時，其參數(shù)估計準確性會變差，嚴重時甚至導(dǎo)致方法失效。此外也有一些無須對振動進行建模的振動補償方法，但也具有一定局限性，如文獻［17-19］提出了基于雙通道回波相位干涉的振動估計補償方法，然而在高頻振動幅度較大、頻率較高時，該方法的估計精度下降。

交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）是一種針對可分解凸優(yōu)化問題的常用算法，通常用于解決存在兩個優(yōu)化變量的等式約束優(yōu)化類問題，具有處理速度快、收斂性能好的優(yōu)點。文獻［20］將ADMM 算法用于稀疏孔徑ISAR 成像中，利用ADMM 解決稀疏約束欠定問題，通過稀疏孔徑數(shù)據(jù)重構(gòu)ISAR 圖像，并在重構(gòu)過程中利用最小熵準則對各脈沖的相位誤差進行估計和補償。該算法具有聚焦性能良好，運算效率高的優(yōu)點。

本文將ADMM 應(yīng)用到太赫茲SAR 成像中，提出了一種新的高頻振動誤差補償算法。首先通過回波信號的方位向快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）得到粗聚焦的圖像，并將振動補償問題轉(zhuǎn)化為使圖像的l1范數(shù)最小化的優(yōu)化問題，然后引入輔助變量將目標函數(shù)構(gòu)造成ADMM 算法所要求的形式，并利用圖像熵最小化原則在迭代過程中引入了相位誤差更新步驟，按順序?qū)D像、輔助變量、拉格朗日乘子和相位誤差四個變量分別求解并循環(huán)迭代。最后利用估計得到的相位誤差進行補償和成像，以獲得聚焦良好的圖像。

本文的結(jié)構(gòu)安排如下：第2 節(jié)建立了具有高頻振動誤差的機載太赫茲SAR 信號模型，并分析了高頻振動誤差對回波相位的影響；第3 節(jié)介紹了基于ADMM 算法的高頻振動補償方法的原理，并給出了方法的實現(xiàn)步驟和流程圖；第4 節(jié)進行了點目標仿真實驗和分布式場景仿真實驗，并將實驗結(jié)果與已有的太赫茲SAR 振動補償方法進行對比，證明了本文方法的有效性；第5節(jié)為結(jié)論。

2 具有高頻振動誤差的太赫茲SAR 信號模型

對于機載SAR，由于載機自身飛行特性以及氣流的影響，平臺在飛行時會產(chǎn)生周期性的振動。平臺的振動接近于多個不同頻率的簡諧運動分量的疊加，因此可以將該周期性振動建模為

其中t為慢時間，L為簡諧運動分量數(shù)，Al、fl和φl分別為第l個簡諧運動分量的振幅、頻率和初相。高頻振動須滿足的條件為

其中fv為平臺振動的最大頻率分量，Ts為合成孔徑時間。以直升機平臺為例，其主諧振分量的頻率通常在10 Hz 到30 Hz 范圍內(nèi)［11］。太赫茲SAR 合成孔徑時間約為0.2 s，因此直升機的振動可以認為是高頻振動。

對于工作在正側(cè)視模式的機載太赫茲SAR，平臺的高頻振動誤差主要分為沿航向誤差和沿視線向誤差兩部分，其中對回波相位造成影響的主要是沿視線向的部分，而沿航向部分的誤差對回波相位的影響通常可以忽略［21］。因此，本文中討論的高頻振動誤差為沿視線向的誤差。

工作在正側(cè)視模式下的機載太赫茲SAR 成像幾何如圖1 所示。飛機在h的高度上以速度V沿y軸方向飛行。成像區(qū)域內(nèi)的一個散射點P與雷達之間的斜距為rc(τ，t)，其中τ為快時間。假設(shè)載機沿視線向的高頻振動為rv(t)，無高頻振動時點P與雷達天線相位中心的最短斜距為r0，零多普勒時間為tp，則點P與雷達天線相位中心之間的距離為

圖1 太赫茲SAR成像幾何Fig.1 Imaging geometry of THz-SAR

其中

用一階泰勒展開對式（4）近似，可以得到

其中rm(t)為與快時間無關(guān)的瞬時斜距，km(t)τ為快時間引入的距離徙動誤差。

對于脈沖體制雷達，脈沖寬度一般較小，可以采用“停-走-停”近似處理，忽略快時間對斜距的影響。但若雷達發(fā)射的信號為線性調(diào)頻連續(xù)波（Linear Frequency Modulated Continuous Wave，LFMCW），須根據(jù)成像質(zhì)量的要求，對比快時間造成的距離偏移誤差與成像的距離向分辨率，以判斷是否可以忽略快時間對斜距的影響。點P對應(yīng)的距離偏移誤差ΔR為［22］

其中Tp為脈沖寬度，Br為信號帶寬，f0為中心頻率。根據(jù)本文中仿真實驗的參數(shù)設(shè)置（見表1），可以得到ΔR的最大值約為0.0027 m，遠小于距離分辨率。因此可以忽略快時間對斜距歷程的影響，點P與雷達天線相位中心之間的距離可以近似為

表1 仿真實驗中的系統(tǒng)參數(shù)設(shè)置Tab.1 System parameters for Simulation

假設(shè)雷達發(fā)射的LFMCW信號的調(diào)頻率為Kr，wr和wa分別為距離向和方位向的窗函數(shù)，則發(fā)射信號為

雷達接收到的點目標P的回波信號為

其中λ為載波波長。可以得到，由于高頻振動的存在，回波信號的相位受到了如下指數(shù)項的調(diào)制，需要對其進行相位誤差補償

將式（12）中的回波信號表示為矩陣形式，記為S∈CN×M，其中N和M分別是方位向和距離向的采樣點數(shù)。通過對信號S進行方位向FFT，可以得到粗聚焦的SAR圖像X∈CN×M，即

其中F1∈CN×N為傅里葉矩陣，N∈CN×M為噪聲矩陣，E1=diag[exp(jφ)]為相位誤差矩陣，diag(·)表示用向量構(gòu)成對角矩陣，φ∈RN代表相位誤差。將式（14）中的矩陣的所有列按順序排成一列，構(gòu)成相應(yīng)的向量x、s、n，可以得到

其中F=IM?F1，E=IM?E1，IM∈CM×M為單位矩陣，“?”表示求克羅內(nèi)克積。此時式（15）可以改寫為

我們可以將振動補償問題轉(zhuǎn)化為如下尋優(yōu)問題

其中，‖·‖1、‖·‖2分別指l1范數(shù)和l2范數(shù)，α是正則化系數(shù)。

3 基于ADMM 的相位誤差補償方法

3.1 ADMM 算法簡介

ADMM 是一種求解可分離凸優(yōu)化問題的重要算法，具有處理速度快、收斂性能好的優(yōu)勢，在機器學習等領(lǐng)域具有廣泛應(yīng)用。ADMM 的目標函數(shù)應(yīng)當為包含兩組可分離自變量，且兩組自變量之間存在線性等式約束，其形式為

其中，ρ為懲罰因子。此時可以得到目標函數(shù)的擴展拉格朗日表達式為

其中u∈CN是拉格朗日乘子。ADMM 算法將求解分為三步，每一步都將其中兩個變量固定，并更新第三個變量的值

其中k表示迭代次數(shù)。ADMM算法要求不斷重復(fù)這三步，直至收斂。

3.2 基于ADMM 的相位誤差估計與補償

式（17）所示優(yōu)化問題中僅包含一組自變量x，此時引入輔助變量z∈CN，將式（17）變?yōu)?/p>

式（22）符合式（18）的形式，可以使用ADMM 算法求解?？梢缘玫綌U展拉格朗日表達式為

分別求擴展拉格朗日表達式關(guān)于x、z的一階導(dǎo)數(shù)，令導(dǎo)數(shù)為0，可以求得令目標函數(shù)值最小的x、z的值，則式（21）被代替為

其中sign(·)表示取符號函數(shù)。事實上，式（24）中相位誤差矩陣E也是未知的，需要進行求解。E可以通過使圖像熵最小來求解。圖像熵是一種圖像聚焦質(zhì)量的評價標準，熵指的是體系的混亂程度，圖像熵越小，代表圖像越清晰，圖像的聚焦程度越好［23-24］。對第k次迭代得到的x(k)進行重構(gòu)，生成的圖像X(k)∈CN×M對應(yīng)的圖像熵為

式（25）和式（26）為二維圖像X(k)中所有像素求和的形式，可以轉(zhuǎn)換為向量x(k)的所有像素求和的形式，以與ADMM算法的形式相對應(yīng)

使用最小熵準則估計相位誤差，即

計算圖像熵關(guān)于第n個脈沖對應(yīng)相位誤差φn的導(dǎo)數(shù)

其中Re(·)表示取實部，*表示取共軛，

其中

其中vec(·)表示將矩陣的所有列按順序排成一列。因此可以得到

令導(dǎo)數(shù)為0，可以得到

其中angle(·)表示取相位。

按照順序執(zhí)行式（24）和式（34），即完成了一次對圖像x、輔助變量z、拉格朗日乘子u和相位誤差φ四個變量的求解，這就完成了本文提出的基于ADMM 算法的相位誤差估計方法的一次迭代。當?shù)_到預(yù)設(shè)次數(shù)或收斂時，最終得到的φ即為式（13）中相位的估計。利用估計得到的φ構(gòu)造相位補償信號并與式（12）的回波信號共軛相乘，即可補償高頻振動產(chǎn)生的相位誤差。

在相位誤差補償后即可使用常規(guī)的成像算法（如ωk 算法）進行成像。本文方法的流程如圖2 所示。在參數(shù)初始化中，z、φ、u中的所有元素均可初始化為0，ρ可設(shè)置為1，α可設(shè)置為μvar(s)，其中μ取0.01，var(·)為求方差算子。

圖2 基于ADMM的振動誤差補償流程圖Fig.2 Flowchart of the proposed vibration error compensation method based on ADMM

4 仿真實驗

為了驗證本文提出的相位誤差補償方法的有效性，本節(jié)對具有高頻振動誤差的太赫茲SAR 進行相位誤差補償仿真實驗，分別模擬了點目標和分布式場景的回波信號，然后進行相位誤差補償和成像。仿真實驗中的系統(tǒng)參數(shù)設(shè)置如表1所示。

4.1 點目標仿真實驗

點目標設(shè)置為3×3的點陣，相鄰點的方位、距離向間距均為5 m。載機的高頻振動設(shè)置為4 個不同頻率的簡諧運動分量之和

其中Al，fl和φl分別為第l個簡諧運動分量的振幅、頻率和初相，數(shù)值如表2所示。

表2 高頻振動參數(shù)設(shè)置Tab.2 Parameters of the high-frequency vibration

在點目標仿真實驗中，將文獻［16］中基于分數(shù)階傅里葉變換（Fractional Fourier Transform，F(xiàn)rFT）并結(jié)合準最大似然法和隨機采樣一致性原理的參數(shù)估計方法作為對比方法，以下簡稱FrFT方法。FrFT方法首先使用FrFT獲得信號的瞬時調(diào)頻率，然后使用離散傅里葉變換得到振動頻率的粗估計，再根據(jù)最小二乘回歸得到振動幅度和初相的粗估計，最后使用準最大似然法搜索得到更精確的估計結(jié)果。根據(jù)估計得到的參數(shù)構(gòu)造相位誤差補償信號，可以補償高頻振動誤差。當高頻振動的簡諧運動分量數(shù)不多于2時，F(xiàn)rFT 方法可以準確地估計振動參數(shù)；但當高頻振動包含更多簡諧運動分量時，由于待估計參數(shù)變多，F(xiàn)rFT方法的估計準確性將會下降，甚至失效。

根據(jù)表1 與表2 中的參數(shù)模擬了點目標的回波信號并添加了信噪比（Signal-to-Noise Ratio，SNR）為5 dB 的高斯白噪聲。分別使用FrFT 方法以及本文方法進行振動補償，補償后的成像結(jié)果分別如圖3（c）和圖3（d）所示，圖3（a）與圖3（b）分別為無高頻振動與高頻振動未補償時的參考圖像。可以看到，F(xiàn)rFT 方法補償后的成像結(jié)果存在明顯的散焦現(xiàn)象，而本文方法有效抑制了鬼影目標的產(chǎn)生，得到了聚焦良好的成像結(jié)果。

圖3 點目標成像結(jié)果Fig.3 Imaging results of point targets

為了對成像質(zhì)量進行更準確的評價，以中間點為例，對點目標的峰值旁瓣比（Peak Sidelobe Ratio，PSLR）和積分旁瓣比（Integrated Sidelobe Ratio，ISLR）進行比較。兩種方法的補償結(jié)果和無振動條件下的成像結(jié)果的方位向響應(yīng)如圖4 所示，方位向響應(yīng)的PSLR 和ISLR 列于表3 中?？梢钥闯觯現(xiàn)rFT方法的成像結(jié)果散焦現(xiàn)象嚴重，無法提取PSLR 與ISLR，而本文方法的方位向響應(yīng)接近理想無振動誤差情況下的成像結(jié)果，證明了本文方法的有效性。

表3 中間點方位向響應(yīng)的PSLR與ISLRTab.3 PSLR and ISLR of azimuth responses of center point

圖4 中間點目標的方位向響應(yīng)Fig.4 Azimuth responses of the center point target

4.2 分布式場景仿真實驗

使用聚焦良好的高分辨率SAR 圖像產(chǎn)生分布式場景的仿真回波，產(chǎn)生回波數(shù)據(jù)的方法與文獻［21，25-26］相同。仿真實驗中使用的參數(shù)參照表1與表2，噪聲的SNR 為0 dB。無高頻振動時的成像結(jié)果如圖5（a）所示，存在高頻振動但未進行振動補償?shù)某上窠Y(jié)果如圖5（b）所示。分別使用FrFT 方法與本文方法進行相位誤差補償，成像結(jié)果分別如圖5（c）和圖5（d）所示。可以看出，F(xiàn)rFT方法補償后的成像結(jié)果具有明顯的散焦現(xiàn)象，而本文方法可以有效抑制鬼影目標的產(chǎn)生，使補償后的成像結(jié)果聚焦。

圖5 分布式場景成像結(jié)果Fig.5 Imaging results of the distributed imaging scene

5 結(jié)論

載機的高頻振動會使機載太赫茲SAR 的成像結(jié)果出現(xiàn)鬼影目標和散焦現(xiàn)象?，F(xiàn)有的高頻振動相位誤差估計與補償方法大多需要對振動誤差進行建模，如建模為2個簡諧運動之和的形式，當實際的振動誤差簡諧運動分量數(shù)增加時，其估計準確度將會受到嚴重影響。同時，傳統(tǒng)的自聚焦方法難以補償高頻振動誤差。因此，需要一種無須對振動進行建模的高頻振動補償方法。本文提出了一種基于ADMM 的太赫茲SAR 高頻振動補償方法，可以在無須建模的情況下補償高頻振動的相位誤差，且在低SNR 下的補償結(jié)果也滿足聚焦成像的要求。對點目標與分布式場景的仿真實驗結(jié)果證明了本文方法的有效性。