基于小波包算法的压缩传感SAR成像方法

时 燕 陈迪荣

(北京航空航天大学数学与系统科学学院 北京 100191)

1 引言

合成孔径雷达(Synthetic Aperture Radar,SAR)具有全天候、全天时、高分辨率等特点，目前已经得到了广泛的应用。合成孔径雷达的高分辨率必然带来高数据率和高数据量，这对数据传送带宽和数据存储带来较大压力。因此，减少雷达数据的数据量是很有必要的。

近年来，Donoho,Emmanuel Candès等人[1－5]提出了压缩传感理论，可以利用低于Nyquist采样率的数据精确重构稀疏信号。近年来已引起图像处理、雷达应用等方面的关注及应用。Baraniuk等人[6]首次提出将压缩传感应用于雷达成像中，分析了其应用的可行性。将压缩传感应用于雷达成像中能够有效地减少数据量，但压缩传感应用的前提是稀疏信号(或有稀疏表示的信号)。文献[7]提出基于小波稀疏表示的压缩感知 SAR成像算法(Compressive Sensing Imaging Approach Based on Wavelet Sparse Representation,WCS)，该算法利用小波的稀疏表示功能，利用压缩传感模型进行SAR成像。该算法在严重降采样情况下对具有较多细节的 SAR场景成像效果不是很好。这很可能源于小波分解只对低频部分分解，而不对高频部分分解。而小波包分解是对低频和高频部分都进行分解。所以小波包分解可能可以更好地表达较复杂的SAR场景。

本文提出基于小波包算法的压缩传感 SAR成像方法，通过对 SAR图像进行分类，对多幅同类SAR图像进行小波包分解，获得一个综合小波包树，对综合小波包树应用小波包算法，得到稀疏表示该类SAR场景的小波包最优基，并通过求解l1范数最小化问题，得到SAR图像。仿真实验结果表明，在严重降采样下，本文方法仍能够得到效果较好的SAR图像。本文的算法虽然比WCS算法的计算复杂度大一些，但在降采样情况下成像效果好于WCS算法。

2 传统的SAR成像原理

SAR成像就是从回波信号中提取地表后向散射系数的2维分布。SAR回波信号可以看作是发射脉冲与地表的后向散射系数卷积的结果，成像处理可看成一个解卷积的过程。

回波信号的模型为

其中σk为散射点k的后向散射系数，N表示散射点的个数。f0为线性调频脉冲信号的载频，rk(η)为散射点k在脉冲发射时刻η与雷达平台间的距离，c为光速，Tp为脉冲重复间隔，Kr为线性调频脉冲信号的时间宽度，nx(τ)为复基带回波信号中由弱散射中心回波信号合成的等效加性噪声。

传统的 SAR成像算法主要是应用匹配滤波相位相乘等来解卷积，比如距离多普勒(Range Doppler,RD)算法[8]。

3 基于小波包最优基稀疏表示的压缩传感SAR成像

3.1 压缩传感SAR成像原理

将回波信号的卷积形式写成如下矩阵形式：

式(4)中Φ是NrNa×N的chirp矩阵(Nr表示距离向采样数，Na表示方位向采样数)，其显示表示如下：

其中

SAR成像的过程可以看作一个线性系统，通过回波信号sR和 chirp矩阵Φ来恢复场景反射系数σ。

稀疏度的定义：若向量s中只有K个元素非零，其它元素均为零，则称向量s的稀疏度为K。

稀疏性的定义：信号中大部系数能量较小，只有少量系数的能量较大。稀疏度越小，稀疏性越大。

若σ的稀疏度为K，当K足够小时，且l0范数最小化问题式(5)具有唯一最优解，则恢复σ就是求解l0范数最小化问题。表示零范数，表示向量σ中非零元素的个数。

l0范数最小化问题需要穷举σ中所有非零值的所有种排列可能，是一个NP难问题。

Candès等人[2]指出：如果要精确重构稀疏度为K的信号σ，测量矩阵Φ必须满足2K阶限制等距(Restricted Isometry Property,RIP)条件[4,5]。当测量矩阵满足RIP性质时，l0范数最小化问题式(5)可以与l1范数最小化问题式(6)等价。l1范数最小化问题是一个凸优化问题。

若σ在某个正交基Ψ下的展开系数为x,x具有很好的稀疏性(即稀疏度小，只有少量元素为零)。则上述问题可转化为优化问题式(7)。

根据压缩传感理论，大多数随机矩阵以高概率满足RIP条件，可以精确重构原稀疏信号。RIP性质并不容易验证，通常用不相干性来间接验证。

文献[6]中提到雷达发射的 chirp波形经过时间平移和频率调制而形成的字典与场景反射系数的稀疏表示基(时间、频率、时频基)是不相干的。文献[9]通过分析 chirp矩阵和随机高斯矩阵的性质，表明chirp矩阵与随机高斯矩阵的M-RIP性质类似。文献[10]也通过大量试验证明，chirp矩阵可以以高概率重构K稀疏信号。所以chirp矩阵Φ可以以高概率恢复满足稀疏条件的场景反射系数σ。

当场景反射系数σ具有较好的稀疏性时，可直接用 chirp矩阵Φ来恢复场景。当场景较为复杂，细节较多时，场景反射系数不具有满足重构条件的稀疏性时，需要一个能够稀疏表示该复杂场景的正交基Ψ使之满足重构条件，应用ΦΨ来恢复场景。

3.2 小波包算法选择稀疏表示基

正交小波变换只对信号的低频部分做分解，而对高频部分，即信号的细节部分不再继续分解，所以小波变换能够表征一大类以低频信息为主要成分的信号，但它不能很好地分解和表示包含大量细节信息的信号，如遥感图像。与小波变换不同的是，小波包变换可以对高频部分提供更精细的分解，对包含大量中、高频信息的信号能够更好的时频局部化分析[11]。SAR场景主要表现地物地貌，包含大量细节信息。所以小波包能够很好地表示包含大量细节信息的SAR复杂场景。

小波包分解树是一个二叉树(如图1)，其中小波包树中的子二叉树是一组小波包基。一个小波包树相当于一个小波包库。压缩传感需要的是稀疏的信号，需要从中选择能够稀疏表示的小波包基。

图1 小波包树(3层小波包分解)Fig.1 Wavelet package decomposition tree (3-layer decomposition)

本文对SAR场景图像进行分类，将具有类似特征的场景图像归类(根据图像内容分类，例如陆地类SAR图像、海岸类SAR图像等)。用小波包算法找出稀疏表示基。主要步骤如下：

第1步：选择L张同类场景图像作为样本，图像展开为1维向量，长度为N。对每个样本进行J层小波包分解，得到L个小波包分解树；

第2步：计算每个样本的小波包树的第k层，第m个节点 β=(m,k),0≤m＜2k,0≤k≤J的系数的稀疏度；

第 4步：应用小波包算法[11]对综合小波包树以平均稀疏度sβ作为信息代价函数值，求出最优树。具体算法为：从综合小波包树的最底层开始，标记最底层的节点，将两个孩子节点的稀疏度之和与它们的父节点的稀疏度进行比较。如果比父节点的稀疏度大，则标记父节点，去掉孩子节点的标记。反之，用两个孩子节点的稀疏度之和代替父节点的稀疏度，依次自底向顶。最终做标记的节点即为最优树，最优树对应的即是小波包最优基。

本文选择了多张128×128陆地类SAR图像作为样本，由于db3小波的压缩性能较好，所以本文选择db3小波3层小波包分解，根据上述算法，得到小波包最优树，如图2所示。该树对应的该类的小波包最优基。仿真用SAR场景图像如图3所示。

图2 陆地类小波包最优树Fig.2 Best wavelet package tree for land scene

图3 SAR场景图像(tier)Fig.3 SAR scene (tier)

图4是SAR场景图像(图3)的灰度值的直方图，该场景的图像的灰度值相对较为平均，稀疏度大，稀疏性小。该场景在小波基下(如图5)和小波包最优基下(如图6)的系数都主要集中在0附近，稀疏度小，稀疏性较强。从图中可以看出，小波包最优基下的系数分布更“陡”一些，系数更集中，稀疏性更好。

从表1中的数据也可以看出，小波系数和小波包稀疏的分布都集中在(0,0.1]范围内，在(0,0.1]范围内，小波包最优基下的系数分布比小波系数多，而在(0.1,1.0]范围内，小波包最优基下的系数分布比小波系数少。同样体现了小波包最优基下的系数比小波基下的系数更集中，稀疏度小，稀疏性更好。

图4 SAR场景图像的灰度值直方图Fig.4 Gray histogram for SAR scene

图5 小波系数分布图Fig.5 Distribution of the wavelet coefficients for the SAR scene

图6 小波包最优基下的系数分布图Fig.6 Distribution of the wavelet package best basis coefficients for the SAR scene

表1 tier场景灰度值、小波系数、小波包最优基系数的稀疏度对比(%)Tab.1 Sparsity Compared from scene gray,wavelet coefficients and wavelet package best basis coefficients for tier scene (%)

3.3 小波包最优基压缩传感SAR成像

小波包最优基压缩传感(Compressive Sensing Based on Wavelet Package Best Basis,WPCS)SAR成像的过程为：将同类型SAR图像根据3.2节中的算法训练得到的小波包最优基作为该类 SAR图像的稀疏表示基Ψ，以ΦΨ为观测矩阵，回波信号sR为观测值，通过压缩传感恢复算法求得场景σ。

WPCS是以小波包最优基作为稀疏表示基Ψ,WCS是以小波基作为稀疏表示基Ψ。虽然小波包最优基的计算量要比直接小波基的计算量大，但是，小波包最优基是对一类图像通用的，只要对每类图像训练一次将该类的小波包最优基存储于计算机内，在成像时，直接使用存储的合适的最优基即可，并不比使用小波基复杂。

4 仿真结果

本文通过对SAR图像用chirp信号进行该场景的回波仿真，用仿真回波数据进行成像处理，比较RD算法、WCS算法[7]、WPCS算法的成像效果。表2是仿真雷达参数。本文中的降采样是对方位向进行降采样，距离向采样不变。方位向采样率决定方位向分辨率。

图7(a)、图7(b)、图7(c)分别为距离多普勒算法全采样、50%采样、30%采样处理回波数据的结果。图7(d)、图7(e)、图7(f)为WCS全采样、50%采样、30%采样处理成像结果。图 7(g)、图 7(h)、图 7(i)是WPCS全采样、50%采样、30%采样处理成像的结果。从图 7中可以看出，RD算法的成像结果图像质量依次下降，30%采样图像发生严重混叠，已模糊不清。可以很明显地看出WCS和WPCS在3种采样情况下均比RD算法成像的效果好。

表2 仿真雷达参数Tab.2 The parameters of radar simulation

在50%采样情况下，RD算法成像有轻微混叠，RD算法图像质量下降较为明显，WCS算法和WPCS算法的图像质量基本没有下降，几乎不影响视觉效果。30%采样的RD算法成像由于混叠严重，已经很模糊了，而WCS和WPCS成像，仍能保持图像的主要特征，无模糊成像。从图像中可以看出30%采样的成像结果中，比较图7(f)和图7(i)，两幅图中从红色方框内的细节中可以看出，图 7(f)中红色方框内的细节已经模糊，看不清物体的轮廓了，图 7(i)中红色方框内仍能看出物体的大致轮廓。图7(i)比图7(f)细节保留得更完整。

表3是图7中RD成像、WCS成像，WPCS成像 3种成像结果的相关度、信噪比(SNR)、峰值信噪比(PSNR)。

图7 SAR图像(tier)成像处理(陆地类)Fig.7 SAR scene (tier) imaging processing (land scene)

表3 图7中各种成像算法的相关度、SNR和PSNRTab.3 The correlation,SNR and PSNR of various imaging algorithms in Fig.7

最小化平均平方误差(Mean Squared Error,MSE)的定义如下：

式中：xi为原始图像的第 i个像素；为处理后的图像的第i个像素。

信噪比(Signal to Noise Ratio,SNR)的定义如下：

信噪比的单位是dB。

图像的峰值信噪比(PSNR)的定义：该指标一般应用在灰度在 0～255的单色图像的降噪评估中。

从表3中可以看出，全采样情况下，WPCS算法的相关度(0.8686),SNR(5.7422),PSNR(22.0874)和 WCS算法的相关度(0.8676),SNR(5.7289),PSNR(22.0741)均比 RD 算法的相关度(0.8153),SNR(4.4745),PSNR(20.8197)高；WPCS算法的相关度比RD算法的高6.54%，比WCS算法高0.12%；WPCS算法的SNR,PSNR比RD算法的高28.33%,6.09%；比WCS算法高0.23%,0.06%。50%采样情况下，WPCS算法的相关度比RD算法的高6.35%,比 WCS算法的高 4.12%；WPCS算法的 SNR,PSNR比RD算法的高22.31%,4.76%；比WCS算法高18.83%,4.11%。30%采样情况下，WPCS算法的相关度比RD算法的高8.64%，比WCS算法高5.89%，WPCS算法的SNR,PSNR比RD算法的高36.24%,7.07%，比WCS算法的高25.23%,5.26%。

图8是海岸类SAR图像成像处理的仿真结果。对于海岸类图像经过 3.2节中的算法得到的小波包最优基与陆地类小波包最优基一样，是如图2所示的小波包最优基树。图8(a)-图8(c)成像质量依次下降，50%采样相干斑较多，海上目标已经被混淆了；30%采样时混叠严重，轮廓模糊。图8(d)-图8(f)和图8(g)-图8(i)图像的大致轮廓都能够很好地保持，但图 8(f)中存在一块较大斑点(红色方框中)，而图8(i)图像的质量效果要好于图8(f)。

在降采样情况下 WPCS算法的成像效果好于WCS算法。

5 结论

本文通过对不同场景类型的 SAR图像在小波基下和小波包最优基下的稀疏度进行分析比较，小波包最优基下的系数比小波基下的系数更具有稀疏性。

本文通过对SAR图像仿真成像，分别用传统的成像算法(距离多普勒算法，RD算法)、WCS算法、WPCS算法进行仿真成像。成像结果表明WCS算法和WPCS算法比RD算法成出的图像质量高。而在降采样情况下，WPCS算法比WCS算法成出的图像效果更好一点。对于较复杂的场景，WPCS算法的优势更明显一些。

[1]Donoho D L.Compressed sensing[J].IEEE Transactions on Information Theory,2006,52(4): 1289-1306.

[2]Candès E,Romberg J,and Tao T.Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[J].IEEE Transactions on Information Theory,2006,52(2): 489-509.

[3]Candès E.Compressive sampling[C].In: Proceedings of International Congress of Mathematicians,Zürich,Switzerland,European Mathematical Society Publishing House,2006:1433-1452.

[4]Candès E and Tao T.Near optimal signal recovery from random projections: universal encoding strategies?[J].IEEE Transactions on Information Theory,2006,52(12):5406-5425.

[5]Candès E and Tao T.Decoding by linear programming[J].IEEE Transactions on Information Theory,2005,51(12):4203-4215.

[6]Baraniukc Richard C and Steeghs Philippe.Compressive radar imaging[C].IEEE Radar Conference,Waltham,MA,April 17-20,2007: 128-133.

[7]王伟伟,廖桂生,吴孙勇,等.基于小波稀疏表示的压缩感知SAR成像算法研究[J].电子与信息学报,2011,33(6):1440-1446.Wang Wei-wei,Liao Gui-sheng,Wu Sun-yong,et al..A compressive sensing imaging approach based on wavelet sparse representation[J].Journal of Electronics & Information Technology,2011,33(6): 1440-1446.

[8]Ian G Cumming,Frank H Wong,洪文,et al..合成孔径雷达成像—算法与实现[M].北京: 电子工业出版社,2007:154-169.Ian G Cumming,Frank H Wong,Hong W,et al..Digital Processing of Synthetic Aperture Radar Data: Algorithms and Implementation[M].Beijing: Electronic Industry Press,2007: 154-169.

[9]Lorne Applebaum,Stephen D Howard,Stephen Searle,et al..Chirp sensing codes: deterministic compressed sensing measurements for fast recovery[J].Applied and Computational Harmonic Analysis,2009,26(2): 283-290.

[10]Tropp J,Wakin M,Duarte M,et al..Random filters for compressive sampling and reconstruction[C].Proc.IEEE ICASSP,2005.

[11]孙延奎.小波分析及其应用[M].北京: 机械工业出版社,2005:245-259.Sun Yan-kui.Wavelet Analysis and Applications[M].Beijing:Machinery Industry Press,2005: 245-259.