郝昭昕 張啟明 孫進(jìn)平 王彥平

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

1 引言

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

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

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

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

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

2 具有高頻振動(dòng)誤差的太赫茲SAR 信號(hào)模型

對(duì)于機(jī)載SAR，由于載機(jī)自身飛行特性以及氣流的影響，平臺(tái)在飛行時(shí)會(huì)產(chǎn)生周期性的振動(dòng)。平臺(tái)的振動(dòng)接近于多個(gè)不同頻率的簡(jiǎn)諧運(yùn)動(dòng)分量的疊加，因此可以將該周期性振動(dòng)建模為

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

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

對(duì)于工作在正側(cè)視模式的機(jī)載太赫茲SAR，平臺(tái)的高頻振動(dòng)誤差主要分為沿航向誤差和沿視線向誤差兩部分，其中對(duì)回波相位造成影響的主要是沿視線向的部分，而沿航向部分的誤差對(duì)回波相位的影響通?？梢院雎裕?1］。因此，本文中討論的高頻振動(dòng)誤差為沿視線向的誤差。

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

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

其中

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

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

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

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

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

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

雷達(dá)接收到的點(diǎn)目標(biāo)P的回波信號(hào)為

其中λ為載波波長(zhǎng)?？梢缘玫?，由于高頻振動(dòng)的存在，回波信號(hào)的相位受到了如下指數(shù)項(xiàng)的調(diào)制，需要對(duì)其進(jìn)行相位誤差補(bǔ)償

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

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

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

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

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

3 基于ADMM 的相位誤差補(bǔ)償方法

3.1 ADMM 算法簡(jiǎn)介

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

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

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

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

3.2 基于ADMM 的相位誤差估計(jì)與補(bǔ)償

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

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

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

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

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

使用最小熵準(zhǔn)則估計(jì)相位誤差，即

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

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

其中

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

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

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

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

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

圖2 基于ADMM的振動(dòng)誤差補(bǔ)償流程圖Fig.2 Flowchart of the proposed vibration error compensation method based on ADMM

4 仿真實(shí)驗(yàn)

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

4.1 點(diǎn)目標(biāo)仿真實(shí)驗(yàn)

點(diǎn)目標(biāo)設(shè)置為3×3的點(diǎn)陣，相鄰點(diǎn)的方位、距離向間距均為5 m。載機(jī)的高頻振動(dòng)設(shè)置為4 個(gè)不同頻率的簡(jiǎn)諧運(yùn)動(dòng)分量之和

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

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

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

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

圖3 點(diǎn)目標(biāo)成像結(jié)果Fig.3 Imaging results of point targets

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

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

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

4.2 分布式場(chǎng)景仿真實(shí)驗(yàn)

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

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

5 結(jié)論

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