基于稀疏重构的全极化SAR联合多维重建

孙 豆 路东伟 邢世其*杨 潇 李永祯 王雪松

①(国防科技大学电子信息系统复杂电磁环境效应国家重点实验室 长沙 410073)

②(国防科技大学电子科学系 长沙 410073)

1 引言

合成孔径雷达(Synthetic Aperture Radar,SAR)具有全天时全天候对目标场景进行观测的优势，日益成为目标识别的关键传感器[1,2]。极化是电磁波的一个重要本质属性，深刻反映目标形状、结构和姿态等物理特征，对于提升目标识别能力具有重要作用[3,4]。近年来，SAR越来越朝着全极化和三维的方向发展[5–8]。三维信息使得雷达对目标的描述越来越精细，极化使得雷达对目标散射特性的描述更加完整，将全极化信息和三维散射信息相结合可以得到人造目标更为全面的重建结果，有利于提高目标分类识别精度。

梳理现有的文献资料，将极化三维重建相关的研究工作进行分类。按照成像方法的不同，相关研究可以分为两类。一类是基于傅里叶变换的成像方法[9,10]，这类方法是传统二维SAR成像方法的三维扩展，具有快捷高效的优势，但是这类方法只能单通道独立成像，无法进行多通道联合同时成像。另一类是基于稀疏重构的成像方法[11–13]，由于稀疏成像本质是一种基于模型匹配的方法，因此这类方法的成像结果分辨率高，且对旁瓣和噪声的抑制效果好，此外，稀疏重构类的方法可以进行多通道联合处理。按照成像方式的不同，相关研究可以分为两类。一类是先进行二维成像再获取高度维信息[9,11,14]，这类方法可以使用现有的二维成像方法，只需进一步对高度信息进行估计，其优势在于计算复杂度低，效率高，但是成像结果受二维图像配准质量的影响。另一类是直接三维成像[10,13,15]，即同时得到目标的三维信息，这类方法的优势在于可以保证数据之间的关联性，不需要进行额外的图像配准，但是直接三维成像会面临巨大的计算压力。按照极化信息的提取方式的不同，相关研究也可以分为两类。一类是各极化通道独立处理[16,17]，即先每个通道分别成像，再提取极化散射矩阵。这种方式的优势在于快捷，适用性广，但是由于忽略了各个极化通道数据之间的关联性，不能保证提取的极化散射矩阵准确。另一类是所有极化通道联合处理[11,18]，即所有通道联合成像，同时得到极化散射矩阵，这种联合的方式保证了极化信息提取的准确度，可以有效克服独立处理中散射中心失配的问题，但是联合处理的计算复杂度较大。

现有文献所采用的方法都是对上述3种分类的不同组合，但是还没有方法将稀疏成像、直接三维成像和极化联合处理这三者相结合。实际上，相比于基于傅里叶变换的成像，稀疏成像在分辨率、抗噪性等方面有着更多的优势[19,20]。相比于直接三维成像，先二维后一维的成像方法很可能破坏散射的一致性，特别是每个二维成像结果中散射点的数目和位置估计可能不同，这种不一致性会使得高度信息的提取不准确。此外，由于每个极化通道观测相同的区域，各通道数据高度相关，且之间存在着丰富的互补信息，三维信息和极化信息同时联合获取不破坏数据之间的关联性，才能充分利用这些互补信息来提高估计精度和分辨率。因此，考虑到稀疏成像、直接三维成像和极化联合处理这三种方式各自的优势，将它们结合进行极化多维重建可以得到质量、精度更高的重建结果。

为了充分利用数据之间的关联性带来的额外信息，实现极化散射矩阵和目标三维信息的同时获取，本文提出基于稀疏重构的全极化联合多维重建方法。首先建立各极化通道的直接三维稀疏重建模型，然后对所有通道设置联合稀疏约束，将全极化三维重建建模为多通道联合稀疏重构问题。考虑到上述稀疏重建问题规模大、计算压力大，本文通过数据插值对模型进行简化，并结合三维快速傅里叶变换、共轭梯度法和牛顿迭代法给出了一种高效的模型求解方法。最后，对三维极化散射矩阵进行极化分解，得到包含目标三维位置信息以及散射类型的多维重建结果。本文提出的全极化联合多维重建方法保证了所有极化通道中散射中心的数目和三维位置一致，有效地克服了散射中心失配的问题，确保了极化特征提取的准确性。此外，该方法大大降低了多维联合处理面临的计算压力，能够在不改变结果精度的情况下有效地获取目标多维重建结果。基于仿真数据和电磁计算数据的实验验证了本文方法的有效性。

2 单通道信号模型

假设雷达距离场景足够远，位于相对于场景中心方位角ϕ和俯仰角θ的位置，并发射中心频率为fc，带宽为BW的观测信号，使用平面波模型，接收信号可以表示为

其中，G(k,0,0)为反射率函数的三维傅里叶变换在(k,0,0)处的值。

以角度(ϕ,θ) 对和(k,0,0)旋转，得到旋转后的(x,y,z)和 (kx,ky,kz)。根据傅里叶变换的旋转特性，式(3)变为

由式(4)可知，反射率函数是波数域采样G(kx,ky,kz)的三维傅里叶逆变换。

3 联合稀疏模型建立

对于每个极化通道，需要被重建的成像场景可以通过式(6)获得

其中，gl(x,y,z;ϕ,θ)为第l(l1,2,···,L)个极化通道的场景反射率函数，Gl(kx,ky,kz)是第l个极化通道的波数域采样。

在三维重建空间中，定义N个位置为候选散射中心

通常，这些位置是从均匀的矩形网格中选择的。基于这些位置，定义M ×N的字典矩阵

其中，m表示波数域采样的索引，n表示C中N个位置的索引。

根据字典矩阵A，将式(6)写成矩阵形式，可得

其中，N ×1的向量βl是要重建的第l个极化通道的三维场景，M ×1的向量bl表示第l个极化通道的波数域采样。

每个极化通道单独三维重建就是要解决下面的直接三维稀疏成像模型

对式(10)中的直接三维稀疏成像模型进行求解则可以同时得到目标的三维信息。相比于先二维稀疏成像再第三维稀疏重建的方法，直接三维稀疏成像避免了二维图像稀疏支撑集的不一致。

对于每个极化通道，由于其观测场景是相同的，且采样方式不变，因此，各极化通道共享相同的字典矩阵A。基于此，建立联合成像模型

其中β[β1···βl···βL]是联合稀疏解矩阵，‖β‖0‖|β1|+···+|βl|+···+|βL|‖0表示稀疏解矩阵β中不为零的行数。

在联合成像模型中，每个极化通道中βl的支撑集是相同的，只是不同极化通道的βl的幅度不同。联合稀疏约束项‖β‖0‖|β1|+···+|βl|+···+|βL|‖0保证了不同极化通道的重建结果共享相同数量和位置的散射中心。联合稀疏为散射系数的估计提供了额外的约束条件，使得在重建过程中可以保留更多有价值的信息。

4 全极化联合多维重建方法

式(11)中的数学模型不是一个凸优化问题，而是一个NP难问题。为了用数值方法求解该模型，首先对目标函数进行放松，定义混合规范

其中βl(i)表示βl中的第i个元素。混合范数ℓ2,p沿不同极化通道的方向计算ℓ2范数，即计算β的张量，并沿散射位置的方向计算ℓp范数。

结合式(12)中的混合范数，式(11)中的模型被松弛为式(13)的优化问题

其中µ是控制稀疏性的正则化参数。联合多维重建就是对式(13)中的模型进行求解。

为了求解式(13)，本文在Cetin等人[21]处理思想的基础上，提出了一种快速迭代方法。定义代价函数

当代价函数J取最小值时，可以得到稀疏重建的结果。计算J对β的偏导数，得到

其中D(β)是N ×N的对角矩阵，b[b1···bl···bL]。

当∇βJ0时，则得到了散射系数的估计。考虑到式(15)中的2AHA+µpD(β)可近似看作Hessian矩阵，因此可以采用近似高斯迭代法求解∇βJ0，得到迭代式：

其中∆n+1是在迭代中的步长，为了使算法快速收敛，在迭代过程中采用变步长∆n+1(∆n)0.9。当小于预设阈值时，算法退出，此时就得到了所有极化通道的三维重建结果。为了进一步提取目标极化特性，还可对式(16)得到的极化散射矩阵进行Cameron分解[22]，最终就得到了包含极化信息和目标三维信息的多维重建结果。

实际上，式(16)的迭代计算还存在一些问题。由于本文方法进行三维直接重建，同时获取目标的全极化三维信息，不可避免地面临计算量巨大的问题，即式(16)中的字典矩阵维度很大，无法直接存储计算。具体地，字典矩阵A是M ×N的，N与场景的大小和分辨率有关,M与波数域的采样点数有关。假设要重建的场景大小为5 m×5 m×5 m，每个维度的分辨率是0.05 m，那么三维重建结果的维度是100×100×100，则N1×106。假设频率的采样点数为100，方位向的采样点数为100，俯仰向的采样点数为100，则M1×106。在这种情况下，由于M ×N1×1012，字典矩阵A需要的存储空间大小约为3GB。除了要存储巨大的字典矩阵之外，每次迭代都不可避免地要进行字典矩阵的乘法运算，如计算Aβl，需要M ×N1×1012次乘法和加法计算，因此直接对式(16)进行迭代的计算复杂度很高。为了解决上述计算规模大的问题，本文借鉴文献[12]中的思路，结合三维非均匀快速傅里叶变换(3-D Non-Uniform Fast Fourier Transform,3-D NUFFT)和共轭梯度法求解式(16)中的部分。

由式(5)可知，kx,ky,kz与f,θ,ϕ之间的关系不是线性的，因此kx,ky,kz处于非均匀网格上。对于字典矩阵虽然x,y,z处于均匀的网格上，但由于kx,ky,kz是非均匀的，因此字典矩阵不可以被看作三维傅里叶变换矩阵。为了可以使用三维快速傅里叶变换(3-D Fast Fourier Transform,3-D FFT)代替字典矩阵，就需要在稀疏重构前对每个极化通道的波数域采样进行预处理，使得波数域采样处于均匀的网格上。

这里使用3-D NUFFT对非均匀网格kx,ky,kz上的每个极化通道的波数域采样Gl(kx,ky,kz)进行插值，插值后的波数域采样处于均匀网格上。对应于插值后的波数域采样，字典矩阵则更新为式(16)更新为

表1总结了本文所提出的全极化联合多维重建方法的具体步骤。

该联合多维重建方法实现了极化散射矩阵和目标三维信息的同时获取，充分利用了数据之间的关联性带来的额外信息。实际上，本文方法的关键在于反映了联合稀疏性的矩阵包含了所有极化通道和全部三维的散射系数。联合稀疏性为散射系数的估计提供了额外的约束条件，保证了不同极化通道、不同维度的稀疏支撑集是一致的。不管散射中心是不是在所有极化通道都有散射信息，联合稀疏对4个极化通道进行联合，同时估计散射中心在各个极化通道的散射信息，有效地克服了独立处理和三维分步成像中散射中心失配的问题，进而保证了获取的极化散射矩阵的准确性。此外，同时获取极化和三维信息不可避免地要面对计算规模大的问题。本文方法结合3-D FFT和共轭梯度法对模型进行简化，避免了巨大字典矩阵的存储，大大降低了计算压力，使联合多维重建可以在不改变解的精度的前提下高效地进行。

表1 全极化联合多维重建方法的步骤Tab.1 Steps of full polarization joint multi-dimensional reconstruction method

5 实验结果与分析

本节基于仿真数据和电磁计算数据开展实验，验证本文提出方法的有效性。实验在一台配备Intel Core I5-6500 CPU和12 GB RAM的计算机上的MATLAB R2016b中进行。

5.1 重建性能分析

选取三面角，偶极子，30°二面角和45°二面角作为仿真目标，进行全极化联合三维重建，并从目标类型依赖性、极化散射矩阵估计精度、噪声敏感度这3个方面分析本文方法的性能。表2给出了4个仿真目标的信息，其中三面角，偶极子和45°二面角的稀疏支撑集不一致，30°二面角的稀疏支撑集一致。回波数据由MATLAB生成，使用理想点散射模型，表3记录了仿真实验时的各项参数。

表2 仿真目标信息Tab.2 Information of simulated targets

表3 目标的仿真参数Tab.3 Simulation parameters of simulated targets

图1给出了仿真目标的4个极化通道的成像结果，其中左列是HH,HV,VH和VV极化通道下的二维投影结果，右列是HH,HV,VH和VV极化通道下的三维结果。以HH极化通道成像结果的峰值为0 dB，图1中所有子图显示成像结果的幅度阈值为–20 dB。观察HH通道的成像结果有3个点目标，对应三面角、偶极子、30°二面角；HV通道和VH通道的成像结果有两个点目标，对应30°二面角和45°二面角；VV通道的成像结果有两个点目标，对应三面角和30°二面角。各个极化通道的成像结果和目标的极化散射矩阵可以对应上，目标极化散射矩阵分量是0的通道的成像结果没有该目标。

图1 各个极化通道的仿真目标成像结果Fig.1 Each polarization channel’s imaging results of simulated targets

图2给出了仿真目标的多维重建结果，其中图2(a)为全极化三维成像结果，图2(b)为Cameron分解结果。对比图2(a)的全极化三维成像结果和图1中的各个极化通道的三维成像结果，可以看出图2(a)融合了所有极化通道的结果，包含了全部的四个目标，对应三面角、偶极子、30°二面角和45°二面角。图2(b)中不同的图标表示由Cameron分解得到的不同的散射类别。对比图2(b)中的结果和表2中的目标信息可知，各散射类型判别结果与实际情况相符合，验证了本文方法对极化散射矩阵提取的有效性。

为了进一步分析本文提出的联合多维重建方法对极化散射矩阵的估计精度，表4给出了仿真目标的极化散射矩阵估计结果，并同时列出了各个极化通道分别独立多维重建的估计结果进行对比。变型前的结果是多维重建得到的极化散射矩阵原始值。为了和目标真实的极化散射矩阵进行对比分析，我们对估计得到的每个散射矩阵进行归一化处理，将幅度和相位提到矩阵外，得到变形后的极化散射矩阵估计值。对于变形后的极化散射矩阵估计值，蓝色部分表示各个目标的幅度估计结果，红色部分表示各个目标的极化散射矩阵估计结果，指数部分表示与当前场景有关的相对相位，不影响结果的分析。从表4中联合多维重建方法得到的变型后的结果可以看出，各个目标的散射幅度相近，这与表2中目标真实的幅度情况一致。由于NUFFT插值和3-D FFT会带来增益的变化，因此幅度的估计结果不是1。与表2中目标真实的极化散射矩阵对比，可以看出对于三面角、偶极子和45°二面角，联合多维重建方法得到的极化散射矩阵的估计值和真实值完全一致，对于30°二面角，联合多维重建方法得到的极化散射矩阵估计结果基本上和真实值一致，交叉极化通道的估计结果略微大于真实值。观察表4中独立多维重建法得到的变型后的结果，可以看出极化散射矩阵的估计结果并不准确。具体地，30°二面角的散射幅度和其他目标的散射幅度差异很大，这与表2中目标真实的幅度情况不一致。此外，对于30°二面角，独立多维重建方法得到的极化散射矩阵估计结果与真实值差异很大，交叉极化通道的估计结果远大于真实值。因此，各个极化通道分别独立多维重建不能保证极化散射矩阵提取的准确性，而联合多维重建方法则可以得到较准确的极化散射矩阵的估计结果。

为分析本文方法对噪声的敏感度，仿真中给波数域采样G(kx,ky,kz)增加噪声，在图像域信噪比(Signal to Noise Ratio,SNR)为13 dB,18 dB,23 dB时分别进行三维联合重建，并记录极化散射矩阵估计结果如表5所示。从表5中变型后的结果可以看出，当SNR为18 dB和23 dB时，蓝色部分表示的各个目标的散射幅度相近，这与表2中目标真实的幅度情况一致。此外，红色部分表示的各个目标的极化散射矩阵估计结果也与表2中目标真实的极化散射矩阵基本一致。当SNR=13 dB时，可以看出各个目标的散射幅度估计结果差异略大，特别是30°二面角和45°二面角的散射幅度估计结果相差0.09。此外，与表2中目标真实的极化散射矩阵对比，可以看出对于三面角、偶极子和45°二面角，其极化散射矩阵的估计值和真实值基本一致，但对于30°二面角，受噪声影响，其极化散射矩阵的估计值和真实值略有差异。

图2 仿真目标全极化联合多维重建结果Fig.2 Full polarization joint multi-dimensional reconstruction results of simulated targets

表4 仿真目标的极化散射矩阵估计结果Tab.4 Polarization scattering matrix estimation results of simulated targets

表5 不同SNR下仿真目标的极化散射矩阵估计结果Tab.5 Polarization scattering matrix estimation results of simulated targets under different SNR

根据上述仿真分析可得，本文方法的性能不受目标类型影响，不论目标的稀疏支撑集是否一致，本文方法都可以较准确地估计出极化散射矩阵。相比较于各个极化通道分别独立多维重建不能保证极化散射矩阵提取的准确性，本文的联合多维重建方法对极化散射矩阵的提取准确且估计精度高。此外，本文方法具有一定的抗噪性，在SNR=13 dB时，极化散射矩阵的估计结果开始变差。

5.2 典型目标重建结果

使用Slicy模型的电磁计算数据进行仿真实验。该模型由三面角、二面角、圆柱等组成，图3给出了其三维CAD模型及尺寸信息。回波数据由CST微波工作室生成，采用弹跳射线法。表6记录了Slicy仿真实验中的各项参数。

图3 Slicy的CAD模型Fig.3 CAD model of Slicy

表6 Slicy的仿真参数Tab.6 Simulation parameters of Slicy

图4给出了Slicy模型的4个极化通道的成像结果，其中左列是HH,HV,VH和VV极化通道下的二维投影结果，右列是HH,HV,VH和VV极化通道下的三维结果。以HH极化通道成像结果的峰值为0 dB，图4中所有子图显示成像结果的幅度阈值为–30 dB，可以看出成像结果稀疏，旁瓣基本都被去除掉，验证了本文基于稀疏重构的重建方法在旁瓣抑制上的有效性。观察两个主极化通道，其成像结果相似，都可清晰地看到Slicy模型的外形轮廓，Slicy模型中的两个帽型结构，一个三面角和两个二面角都被重建出来，成像结果的位置和尺寸与图3中的模型吻合。观察两个交叉极化通道，目标的散射相对较弱，轮廓不太清晰，有些部分没有被重建出来，且两个交叉极化通道的结果略有不同，Slicy模型的一个帽型结构的强度有些差异。

图4 各个极化通道的Slicy成像结果Fig.4 Each polarization channel’s imaging results of Slicy

图5给出了Slicy模型的多维重建结果，其中图5(a)为全极化三维成像结果，图5(b)图为Cameron分解结果。对比图5(a)的全极化三维成像结果和图4中的各个极化通道的三维成像结果，可以看出图5(a)和图4中的主极化通道的结果基本一致，不存在散射中心失配的问题，这也就验证了本文的联合重建方法的有效性，即保证了不同极化通道、不同维度的稀疏支撑集一致。图5(b)中不同的图标表示由Cameron分解得到的不同的散射类别。对比图5(b)中的结果和图3中的模型可知，各散射类型判别结果基本与实际情况相符合，验证了本文方法对极化散射矩阵提取的有效性。但由于Cameron分解自身的局限性，散射类型判别存在一定的不确定性，Slicy模型中的二面角和帽型结构没法分辨。

图5 Slicy全极化联合多维重建结果Fig.5 Full polarization joint multi-dimensional reconstruction results of Slicy

使用卫星模型的电磁计算数据进行仿真实验，图6给出了其三维CAD模型，卫星中心位于坐标原点，太阳能两帆板外侧的距离为8.21 m，星体最长尺寸3.55 m。回波数据由电磁计算软件[23]生成，采用积分方程快速计算方法。表7记录了卫星仿真实验中的各项参数。

图7给出了卫星的4个极化通道的成像结果，其中左列是HH,HV,VH和VV极化通道下的二维投影结果，右列是HH,HV,VH和VV极化通道下的三维结果。以HH极化通道成像结果的峰值为0 dB，图7中所有子图显示成像结果的幅度阈值为–30 dB，可以看出成像结果稀疏，卫星的轮廓很清晰，旁瓣基本都被去除掉了，进一步验证了本文基于稀疏重构的重建方法在旁瓣抑制上的有效性。观察两个主极化通道，都可清晰地看到卫星模型的外形轮廓，卫星主体的柱状结构和卫星两侧的太阳能帆板都被重建出来，与图6中的模型吻合。不过两个主极化通道的成像结果略有不同，HH极化通道中卫星两个帆板的散射强度明显大于VV极化通道。观察两个交叉极化通道，其成像结果相似，目标的散射相对较弱，轮廓不太清晰，卫星两侧的太阳能帆板没有被重建出来，只有卫星主体柱状结构被重建出来。

图6 卫星的CAD模型Fig.6 CAD model of satellite

图8给出了卫星的多维重建结果，其中图8(a)为全极化三维成像结果，图8(b)为Cameron分解结果。对比图8(a)的全极化三维成像结果和图7中的各个极化通道的三维成像结果，可以看出图8(a)和图7中的HH极化通道的结果基本一致，因此不存在散射中心失配的问题，进一步验证了本文联合重建方法的有效性。图8(b)中不同的图标表示由Cameron分解得到的不同的散射类别，可以看出卫星的两个帆板上的散射类型为圆柱体，两个帆板下沿以及帆板与卫星主体连接处的散射类型为平板/三面角，卫星主体上部的散射类型为圆柱体，卫星主体中部比较复杂，散射类型较多，卫星主体下部的散射类型为平板/三面角和帽型/二面角的组合。由于Cameron分解自身的局限性，二面角和帽型结构，三面角和平板没法分辨。对比图8(b)中的结果和图6中的模型，各散射类型判别结果基本与实际情况相符合，进一步验证了本文方法对极化散射矩阵提取的有效性。

表7 卫星的仿真参数Tab.7 Simulation parameters of satellite

图7 各个极化通道的卫星成像结果Fig.7 Each polarization channel’s imaging results of satellite

图8 卫星全极化联合多维重建结果Fig.8 Full polarization joint multi-dimensional reconstruction results of satellite

6 结论

本文提出一种基于稀疏重构的全极化联合多维重建方法。该方法对所有极化通道和所有维度的散射系数进行了联合稀疏约束，不仅保证了不同极化通道、不同维度的稀疏支撑集一致，并且充分利用了数据之间的关联性带来的额外信息。此外，该方法解决了同时获取极化和三维信息面临的计算规模大的问题，使联合多维重建可以在不改变解的精度的前提下高效地进行。实验结果表明本文方法的性能不受目标类型影响，具有一定的抗噪性，可以实现极化散射矩阵和目标三维信息同时获取，有效地克服了散射中心失配的问题，得到了准确的目标多维重建结果。此外，该方法保持了稀疏重构的优势，有效地抑制了旁瓣，得到的重建结果分辨率高。