基于近似观测和最小熵约束的SAR稀疏自聚焦方法

熊旭颖, 李 根, 马彦恒

(陆军工程大学石家庄校区无人机工程系, 河北 石家庄 050003)

0 引 言

压缩感知和稀疏重构理论表明,当成像场景稀疏或可以被稀疏表示时,可以利用少量的雷达回波数据获取高分辨和高质量的合成孔径雷达(synthetic aperture radar,SAR)图像,稀疏SAR成像技术能够有效降低雷达接收机的数据采集与存储的压力,相关的成像方法受到了广泛关注[1-4]。

传统的稀疏SAR成像方法需要根据平台的运动轨迹构建测量矩阵。对于机载SAR成像,由于惯导设备精度有限,平台的运动轨迹通常难以精确测量,运动误差的存在将导致成像理论模型和测量矩阵失配,这种失配将产生严重的散焦现象。因此,基于回波数据的自聚焦是稀疏SAR成像过程中必不可少的环节。Onhon等人在2012年提出了一种稀疏驱动的自聚焦方法[5],该方法将误差相位引入测量矩阵，构建了一种稀疏自聚焦模型,并将模型的迭代求解分为SAR场景重建和误差相位估计两步,能够校正一维和二维的相位误差,该方法也被称为基于两步优化的稀疏自聚焦方法。基于两步优化的自聚焦思想,Yang等人提出了能够估计观测位置误差的压缩感知SAR成像方法[6]。Chen等人提出了一种基于参数化稀疏表示的SAR运动补偿方法[7],进一步提高了成像质量。以上这些方法的自聚焦模型均是在时域构建测量矩阵,这种测量矩阵又被称为精确观测算子。当成像场景较大时,精确观测算子将产生巨大的运算量和内存占用量[8-9],限制了稀疏SAR成像方法的应用。近年来,中科院电子所吴一戎院士团队研究了基于近似观测算子的稀疏微波成像方法[7，10],成功地将传统的匹配滤波类SAR成像算法和压缩感知理论相结合,极大地减小了测量矩阵的尺寸,使稀疏SAR成像方法能够用于大场景成像,相关的成像算法受到了广泛的关注[11-13]。对于运动误差条件下稀疏SAR成像,Li等人将近似观测算子引入稀疏自聚焦模型中,用于替代精确观测算子,极大地降低了稀疏自聚焦方法的运算量[14]。Pu等人基于近似观测算子构建了一种基于幅相误差联合校正的稀疏自聚焦模型,用于校正空变的相位误差,在全采样数据下有效估计了三维的运动误差偏离[15]。

以上这些基于两步优化的稀疏自聚焦方法均采用最大似然估计器(maximum likelihood estimator,MLE)估计相位误差,在该模型中,误差相位在二维时域(回波数据域)引入,难以获取误差相位的初始解。在基于两步优化的稀疏自聚焦方法中,合适的误差相位初始解对算法性能至关重要,当回波数据采样率较低且重建散射点数量较多时,MLE收敛缓慢并且容易陷入误差较大的局部最优解,导致自聚焦失败。本文基于近似观测算子构建了一种新的稀疏自聚焦模型,在该模型中,误差相位在聚焦图像的傅里叶变换域引入,因此相位梯度自聚焦(phase gradient autofocus,PGA)[16]等经典匹配滤波类自聚焦方法可用于提供误差相位的有效初始解。所提模型基于两步迭代的思想采用坐标下降法进行求解,为加快算法收敛速度,快速迭代阈值收缩(fast iterative shrinkage thresholding,FIST)[17]被用于重建场景散射系数,同时引入了场景散射系数的最小熵约束,这一附加的先验信息在加快算法迭代收敛速度的同时,使迭代结果更接近全局最优解。

1 传统的稀疏自聚焦模型

假设成像区域中共含有PQ个散射点,在一个合成孔径时间内,SAR接收机采集到的离散二维回波数据可表示为

exp[-j4πR0(k,ηm)/λ]，

m=1,2，…，M;n=1,2，…，N

(1)

式中：ηm和τn分别表示离散的方位和距离采样点；M和N分别表示距离和方位向上的采样点数；g(k)表示第k个散射点的复散射系数；R0(k,ηm)表示成像区域中的第k个散射点在方位采样时刻ηm的瞬时斜距；c和λ分别表示光速和发射信号波长；p(τ)表示发射的基带线性调频信号,具体为

p(τ)=rect(τ/Tr)exp(jπKrτ2)

(2)

式中：Tr表示发射信号的脉冲持续时间；Kr表示距离调频率。

在考虑噪声的情况下,式(1)可矩阵化表示为

s=Ξg+n0

(3)

式中：s∈CMN×1表示由2维回波组成的1维回波向量；Ξ=HA∈CMN×PQ表示精确观测矩阵,其中H∈CMN×MN为降采样矩阵,是对角元素仅含有0或1的对角矩阵,A∈CMN×PQ为测量矩阵,其具体形式见文献[1]；g∈CPQ×1为散射点组成的1维向量；n0为1维的加性高斯白噪声向量。

式(3)是基于精确观测的压缩感知模型,当成像场景稀疏时,散射系数向量g可以通过求解如下所示的L1范数正则化问题获得：

(4)

式中：γ表示正则化参数,由成像场景的稀疏度决定。

当SAR回波数据存在非空变的相位误差时,常规的稀疏自聚焦模型在二维时域引入误差相位[5],表示为

(5)

其中，Ψ为相位误差矩阵,具体表示为

(6)

式中:φ(m)表示第m个方位采样点的误差相位。

文献[5]将式(5)的求解转化为两步优化问题,并采用坐标下降法进行求解,具体的求解过程如算法1所示。

算法 1 基于精确观测的稀疏自聚焦算法初始化: i=0, g0=(Ψ0Ξ)Hs且Ψ0=Ψ0|ϕ(m)=0步骤 1 gi+1=argming12s-ΨiΞg22+γg1{}步骤 2 Ψi+1=argminΨ12s-ΨΞgi+122+γgi+11{}步骤 3 i=i+1,然后回到步骤1。当gi+1-gi22gi22小于预设的门限值时,算法停止。

算法1中的步骤1在已知相位误差的基础上,对场景散射系数进行重建,是一个L1范数正则化问题,可采用迭代软阈值(iterative soft thresholding,IST)[18]等方法进行求解;步骤2在已知场景散射系数的基础上估计误差相位,可以通过最小化均方误差,利用MLE进行求解,则第m个方位采样点的误差相位为

(7)

其中，[Ξ]m和[s]m分别表示第m个方位采样点的观测子矩阵和回波数据。

算法1中的稀疏自聚焦模型需要将2维的散射场景和回波数据拉成1维向量,精确观测矩阵Ξ的维度为MN×PQ,这将产生巨大的极大内存占用和运算量,限制了其在大场景成像中的应用。近似观测算子的提出极大地降低了观测矩阵的维度,提高了稀疏重构效率[11]。接下来,以线性调频变标算法为例,简要说明近似观测算子的构造过程。

图1给出了线性调频变标算法的成像流程图,其中，FFT表示快速傅里叶变换，IFFT表示逆快速傅里叶变换。用Γ(·)算子表示该成像过程,则Γ(·)可以矩阵化表示为

图1 线性调频变标算法成像流程图Fig.1 The flow chart of chirp scaling algorithm for the linear frequency modulation

(8)

从式(8)中可以看出Γ(·)是仅含矩阵乘法和矩阵点乘运算的线性可逆算子,原始回波S可通过其逆运算获得,即

(9)

其中，*表示矩阵的共轭。

在式(9)中，回波数据通过匹配滤波成像的逆过程获得,Γ-1(·)被称为近似观测算子。计算Γ(·)和Γ-1(·)所需的维度与原始回波数据一致,极大地降低了观测算子的内存占用和运算量。

为降低稀疏自聚焦模型中观测矩阵Ξ的维度,可通过引入近似观测算子构造如下式所示的稀疏自聚焦模型：

(10)

式中：L是一个仅包含0和1的二进制随机矩阵,表示回波数据的欠采样矩阵。设ξa表示方位向上的回波数据抽取率,ξr表示距离向上的回波数据抽取率,则L共有Nξa个非零列,每一列有Mξr个1。Φ表示相位误差矩阵,具体形式为

Φ=diag[ejφ(1),ejφ(2),…,ejφ(M)]

(11)

式(10)中基于近似观测的稀疏自聚焦模型同样可以采用算法1进行求解,尽管该模型利用近似观测算子降低了内存占用量,但误差相位依旧在回波数据域引入,难以采用常规的自聚焦方法获取有效的误差相位初始解。误差相位需要从0开始采用MLE进行迭代更新,当采样率较低且重建的散射点数量较多时,该模型的迭代过程容易陷入误差较大的局部最优解。

2 基于近似观测和最小熵约束的稀疏自聚焦

2.1 新的稀疏自聚焦模型

(12)

其中

(13)

在式(12)的稀疏自聚焦模型中,相位误差Φ在聚焦图像的傅里叶变换域引入,常规的匹配滤波类自聚焦方法可以用于提供误差相位的初始解。

2.2 基于两步优化的模型求解方法

基于算法1中两步优化的思想,式(12)中所提的稀疏自聚焦模型可以采用坐标下降法进行求解,求解过程的伪代码如算法2所示,主要内容包括误差相位初始解Φ0的获取,以及G和Φ的重建方法,具体的推导过程如下。

算法 2 基于近似观测和最小熵约束的稀疏自聚焦算法初始化: i=0, Imax,V0=0, t0=1,Φ0=PGA(Γ(L☉S))步骤 1 基于FIST的场景散射系数重建Gi+1=argminGJ(G, Φi)步骤 1.1 Vi+1←Ψλ,δ(Vi-δFη(Φ∗iΩ(L☉[Ω-1(Φi·F-1(Gi))-S])))步骤 1.2 ti+1=1+1+4t2i2步骤 1.3 Gi+1=Vi+1+ti-1ti+1(Vi+1-Vi)步骤 2 相位误差估计Φi+1=argminΦJ(Gi+1, Φ)步骤 2.1 ϕi(m)=angle([Ω(S)]m·[F-1η(Gi+1☉(1+log|Gi+1|2))]Hm)步骤 2.2 Φi=diag[ejϕi(1),ejϕi(2),…,ejϕi(M)]步骤 2.3 i=i+1,然后重复步骤1～步骤5,当i=Imax时,算法停止。

2.2.1 误差相位初始解的获取

采用坐标下降法交替求解G和Φ时,仅能保证迭代收敛到局部最优解。因此,合适的Φ0能够减少算法的迭代次数并使收敛结果更接近全局最优解。在所构建的稀疏自聚集模型中,PGA和最小熵等方法均可以用于提供Φ0。但在稀疏采样的情况下,点目标聚焦后的点散布函数会出现较高的副瓣电平,大量高副瓣电平的存在会使聚焦后的SAR图像存在大量的欠采样噪声,导致图像熵值增加,因此最小熵方法不适用于稀疏采样条件下Φ0的估计。稀疏采样对强散射点的主瓣影响较小,在PGA方法中,随着迭代次数的增加,窗长度逐渐减小,欠采样产生的高副瓣电平可以被滤除,PGA方法更适用于稀疏采样条件下Φ0的估计。

因此,在算法2中,Φ0由PGA方法对稀疏采样的匹配滤波成像结果Γ(L⊙S)进行自聚焦得到。与全采样条件下的自聚焦过程相同,应用PGA方法时,仅需选取约10%的平均能量较大的距离6单元进行自聚焦,最大迭代次数可固定为10次,具体的算法实现过程可参考文献[17]。

2.2.2G和Φ的交替重建

在算法2的步骤1中,利用给定的相位误差Φ来重建图像,则式(12)的优化问题可以表示为

(14)

式(14)是一个L1范数正则化问题,相较于IST方法,FIST方法[16]具有更快的收敛速度,因此算法1结合FIST方法进行求解。在步骤1中,Ψγ,δ(·)表示软阈值算子,δ为搜索步长,通常取δ=1来保证算法收敛,参数γ的选择取决于成像场景的稀疏度K0。本文取γ=|Vi|(K0),其中,|Vi|(K0)是矩阵Vi所有元素幅度排序(由大到小)第K0的幅度值。

在步骤2中,利用给定的场景散射系数Gi来估计相位误差Φi,优化问题可以表示为

(15)

在该优化问题中,Gi是常数,则有

(16)

基于最小方差和MLE求得第m个观测位置的误差相位为

(17)

式中：SΩ=Ω(S)表示相位校正的回波数据;[X]m表示矩阵X的第m行。

为进一步加快求解过程中的收敛速度,本文利用成像场景稀疏的先验信息,在重建相位误差Φ时采用了见成像场景的最小熵约束,如步骤2.1所示,具体推导过程见第2.3节。

2.3 成像场景的最小熵约束

在稀疏SAR成像中,成像场景中散射强度较弱的背景信息可以看作是噪声,因此SAR成像场景可以看作是稀疏的。在相同稀疏度的条件下,运动误差将导致图像散焦,使稀疏SAR图像的熵增加。本节将通过对重建的稀疏SAR场景进行最小熵约束来进一步提高误差相位的重建精度,从而减少迭代次数,提高算法稳定性。

采用匹配滤波原理成像时,相位补偿后的聚焦SAR图像可以表示为

G=Fη(ΦSΩ)

(18)

根据式(18)可得,二维散射矩阵G的第m行第n列元素G(m,n)可以表示为

(19)

用于估计误差相位的图像熵表示为

(20)

(21)

EG对相位误差φ(l)的偏导数表示为

(22)

由于图像熵EG是相位误差的函数,熵值最小时,熵对每一点相位误差的偏导数为0,即

(23)

由于|G(m,n)|2=G(m,n)*G(m,n),所以

(24)

根据式(19)可知,G(m,n)相对于φ(l)的偏导数为

(25)

将式(25)代入式(24),再将式(24)代入式(22)可得

(26)

令式(26)等于0,可得误差相位φ(l)的最大似然估计为

(27)

式(27)是基于图像最小熵的误差相位最大似然估计方法,对比式(27)和式(17)可以发现,在稀疏自聚焦框架中,求解误差相位时,对重建的场景散射系数Gi+1乘以(1+ln|Gi+1|2)的权系数即可实现最小熵约束,即步骤2.1。

对于任意的SAR成像模式,在获取稀疏采样的原始数据L⊙S和频域成像算子Γ(·)后,采用算法2即可实现稀疏自聚焦成像。

3 实测数据验证

本节通过机载SAR的实测数据来验证所提算法的有效性,选用文献[14]中基于近似观测的稀疏自聚焦方法作为参考算法,采用图像熵来评估SAR图像的聚焦质量,在相同的稀疏度下,熵越小图像质量越高。除此之外,取全采样条件下的PGA自聚焦结果作为全局最优解,用于衡量稀疏采样条件下算法的自聚焦效果。与参考算法相比,在每一步的迭代中,所提算法2需要额外进行一次用于最小熵约束的矩阵点乘运算和步骤2.1中用于加快收敛速度的矩阵加法运算,相对于迭代过程中近似观测算子的运算,这些附加的运算量可以忽略。因此,所提方法和参考算法单次迭代的运算量是近似一致的,算法效率可通过收敛所需的迭代次数来体现。

为便于直观比较算法的收敛速度和精度,本文算法和对比算法中的迭代次数均固定为30,为便于比较图像的聚焦效果,所有图片的最大亮度一致。实测数据对应的机载SAR参数如表1所示,采用线性调频变标算子作为近似观测算子。本文的算法模型能够用于距离和方位二维稀疏采样回波数据的自聚焦成像,为便于分析回波采样率对算法收敛性的影响,实验过程中仅在方位向上对实测数据进行抽取。方位向上的回波数据抽取率为ξa,ξa=1表示全采样回波数据。在表1对应的实测数据中,分辨率约为0.35 m,全采样回波数据的方位向过采样率为1.2(即PRF与奈奎斯特采样率的比值,其中奈奎斯特采样率为回波信号的方位向带宽)。因此,在稀疏采样的条件下,回波信号相对奈奎斯特采样的欠采样率为1.2ξa,即当ξa=0.5时,欠采样率为0.6。成像场景的相对稀疏度定义为α=K0/PQ,表示重建的散射点占场景所有散射点的比例。

表1 机载SAR参数Table 1 Parameters of airborne SAR

图2给出了全采样和方位稀疏采样(ξa=0.5)的实测数据成像结果,其中图2(a)为成像场景的谷歌地图。对比图2(b)和图2(c)可以看出,在全采样数据下,PGA可以实现理想的运动补偿。而从图2(d)中可以看出,回波数据的方位向稀疏采样使PGA聚焦后的图像存在大量欠采样噪声,严重降低图像质量。

图2 实测SAR数据成像结果Fig.2 Imging results of real SAR data

为同时补偿运动误差和抑制欠采样噪声,分别采用本文和文献[14]中的稀疏自聚焦方法对该回波数据进行处理,为比较算法对场景散射点的重构能力,固定ξa为0.5,相对稀疏度α分别设置为0.1和0.5,自聚焦结果分别如图3和图4所示。从图3(a)和图3(b)中可以看出,当重建的散射点数量较少时(α=0.1),两种算法都可以较好地估计和补偿运动误差并抑制欠采样噪声,图3(c)表明,两种算法在30次迭代内都收敛到了全局最优解附近,而所提方法具有更快的收敛速度,仅需约7～8次迭代即可收敛。对于多数SAR图像,仅重建少数强散射点将导致图像信息损失。在欠采样的条件下,先验信息的利用程度将极大的影响散射点数量的重建能力,从图4(a)和图4(b)中可以看出,参考方法仅利用了成像场景稀疏的先验信息,当重建的散射点数量大量增加时,场景的稀疏性将变弱,在欠采样的条件下重建的成像场景将有较大误差,进而导致参考方法不能准确地估计和补偿运动误差,致使迭代结果依旧存在明显的方位散焦。由于利用了场景最小熵的附加先验信息,且提供了更精确的初始误差相位,所提方法在欠采样条件下能够有效重建更多的散射点。观察图4(c)可以发现,随着迭代次数的增加,参考方法逐渐收敛到误差较大的局部最优解,而所提方法则可以快速收敛到全局最优解附近,具有良好的自聚焦效果,极大地减少了图像信息损失。

图3 相对稀疏度α=0.1时的稀疏自聚焦结果Fig.3 Sparse autofocuing results with relative sparsity α=0.1

图4 相对稀疏度α=0.5时的稀疏自聚焦成像结果Fig.4 Sparse autofocuing results with relative sparsity α=0.5

在本文方法中,构建了一种新的稀疏自聚焦模型,使得PGA方法能够被用于提供相位误差的初始解,最小熵(minimum entropy, ME)约束能被用于为相位误差的最大似然估计提供更多的先验信息,同时所提方法采用FIST方法重构场景散射信息,进一步加快了收敛速度。为验证PGA提供初始解和ME约束在自聚焦算法中的重要性,在本文的稀疏自聚焦模型的基础上分别构建了不含PGA和ME约束的FIST+MLE算法,含有PGA提供初始解的PGA+FIST+MLE算法,以及含有PGA和最小熵约束的PGA+FIST+ME+MLE算法。其中FIST+MLE方法的成像效果可近似看作为文献[14]中的算法,但事实上,受益于FIST方法的高效率,FIST+MLE的收敛速度会好于文献[14]中的IST+MLE方法。图5和图6分别给出了方位采样率和相对稀疏度对FIST+MLE,PGA+FIST+MLE和PGA+FIST+ME+MLE 3种算法收敛速度影响的实测数据分析结果。可以看出,本文方法在不同的采样率和相对稀疏度下都有很好的稳定性。从图5(b)和图5(c)中可以看出,含有PGA提供初始解的两种算法的初始图像熵较小,这表明在相对较高的采样率下,PGA方法可以提供有效的初始解,能够减少算法所需的迭代次数。同时可以看出,在PGA提供初始解的基础上,随着方位采样率的升高,自聚焦所需的迭代次数迅速减小。从图5(a)中可以看出,当采样率较低时,欠采样噪声将降低PGA方法对初始解的估计精度,此时自聚焦算法所需的迭代次数较多,ME约束能够明显提高收敛速度,并避免算法陷入误差较大的局部最优解。从图6中可以看出,PGA提供初始解以及ME约束对提高算法收敛速度避免陷入误差较大的局部最优解的作用非常显著,随着场景重建散射点数量的增多,PGA+FIST+MLE方法将逐渐收敛到误差较大的局部最优解,而PGA+FIST+ME+MLE方法的收敛结果更接近全局最优解。

图5 α=0.4时,采样率对收敛速度的影响Fig.5 Influence of the sampling rate on convergence rate with α=0.4

图6 ξa=0.5时,相对稀疏度对收敛速度的影响Fig.6 Influence of the relative sparsity on convergence rate with ξa=0.4

4 结束语

基于近似观测的稀疏SAR成像方法将匹配滤波类成像算法和稀疏成像方法结合在一起,极大提高了稀疏成像方法的效率,在此基础上本文提出了一种基于近似观测和最小熵约束的稀疏自聚焦方法,该方法将匹配滤波类成像算法中常用的PGA和最小熵等自聚焦方法结合到SAR的稀疏自聚焦成像中。实测数据成像结果表明,与常规的稀疏自聚焦方法相比,所提方法进一步利用了成像场景的最小熵约束,且为误差相位提供了更精确的初始解,能够在稀疏采样的条件下,有效增大散射点的重建数量,减小迭代次数,使迭代结果更接近全局最优解。本文所构建的稀疏自聚焦模型仅能校正非空变的相位误差,下一步将进一步研究空变相位误差下的稀疏自聚焦方法。