基于交替方向乘子法的Capon层析SAR成像方法

刘 慧，郭馨宇，郭子夜，程碧辉

（北京建筑大学电气与信息工程学院，北京 102616）

0 引 言

合成孔径雷达层析成像技术（TomoSAR）是一种多基线干涉测量技术，通过多次航过目标，在不同高度形成基线，构成第二合成孔径，反演垂直于视线方向即斜距垂向的目标后向散射系数［1-2］。通过谱估计的方法可以得到目标散射体斜距垂向的位置，同时还可以得到后向散射体的散射强度信息，实现斜距垂向的合成孔径成像［3］。

傅里叶变换法是最早被应用于层析SAR 三维成像的方法［4］，通过对同一距离-方位单元的数据序列进行傅里叶逆变换来获取高度向上的后向散射系数，实现三维成像，但此算法受限于奈奎斯特采样定理对均匀采样的要求。根据瑞利准则［5］，斜距垂向上的分辨率与合成孔径大小成反比，这就需要保证较多且较长基线，即足够多的航过次数。但是多次航过成本高且难以实现，实际中航过间隔甚至超过了奈奎斯特采样定理的限制，且非均匀基线会导致斜距垂向的低分辨率、高旁瓣问题。不同于冰原、森林等自然地貌，城市场景的三维成像需要高分辨率二维影像的支撑［6］。如何从多次航过且不满足奈奎斯特采样定理的高分辨率SAR影像中反演出斜距垂向的后向散射系数，实现城市场景的三维成像是层析SAR 成像算法研究的关键。

目前，有两种层析SAR 成像的主流算法。一种是压缩感知［7-11］（Compressive Sensing,CS）算法，另一种是谱估计算法［12-13］。压缩感知能够突破Nyquist 采样定理的限制，使用稀疏的测量信号获取信号离散采样值，实现稀疏信号的重构［14］。在层析SAR 成像问题中，待恢复的后向散射系数在斜距垂向上呈稀疏分布，因此可采用压缩感知算法实现层析SAR 的三维成像。层析SAR 三维成像的压缩感知算法主要有凸优化算法［8］以及迭代贪婪算法［11］等。CS算法根据可解析的散射体数量进行三维重构，但是受散射体数量的限制［7］。

而谱估计算法不考虑散射体的个数，直接根据散射体的功率谱强弱判别为目标或噪声。其中基于波束形成理论的一类算法，从噪声抑制的角度出发，通过加权处理，使信号输出噪声方差最小［15］的约束，计算获得最佳加权，进而实现对斜距垂向信号的估计。文献［16-17］采用波束形成的思想，在二维影像较少的情况下，通过谱估计，恢复斜距垂向上的后向散射系数，从而实现TomoSAR三维成像。此类基于谱估计的非参数化波束形成方法也可以看作是一种光谱分析信号估计器，非常适合处理分布式散射体的成像问题，它的分辨率很大程度上取决于真实孔径，即基线的总跨度，能够改善傅里叶变换成像的低分辨、高旁瓣、成像效果受系统误差影响大等缺点，具有超分辨能力［18］。

基于波束形成的谱估计算法主要有标准Capon 波束形成算法（Standard Capon Beamformer,SCB）［18］、双约束鲁棒Capon 波束形成算法（Doubly Constrained Robust Capon Beamforming,DCRCB）［19］。波束形成的谱估计算法性能优劣主要可以从旁瓣高度、主瓣宽度及稳健性三个方面来评价:低旁瓣可以有效抑制干扰，较窄的主瓣宽度可以提高目标的分辨能力，高稳健性可以减小各种失配对信号反演结果造成的影响。文献［12］采用SCB 算法解决了层析SAR 成像的问题，能够得到高分辨率的同时还拥有更强的抗干扰能力，但是稳健性较差；但在实际场景中，受导向向量和噪声协方差矩阵误差的影响，该方法鲁棒性变差，信号处理性能严重下降。文献［19］提出DCRCB 算法，在SCB 算法的基础上，提高了鲁棒性。

本文提出一种改进的波束形成优化算法，在DCRCB 的基础上，结合L1 范数的约束函数，构建交替方向乘子法［20］（Alternating Direction Method of Multipliers,ADMM）的代价函数，将DCRCB 恢复的后向散射系数进行进一步稀疏优化，实现层析SAR 的三维成像。ADMM 算法以增广拉格朗日算法为基础，将较为复杂的全局求解问题转换为两个或多个更易求解的简单局部子问题。ADMM 算法在迭代中，各子问题可分别完成稀疏重构和降噪运算，被分离的局部子问题代数式都较为简单，均能较容易地求出确定的解，且不必对其进行收敛运算与约束操作。因此，ADMM 算法具有重建精度高的优势。

本文在第1 节中构建了层析SAR 模型及反演层析SAR成像的波束形成理论；第2节介绍基于压缩感知和谱估计方法的层析SAR 成像方法，并提出改进的层析SAR成像算法；第3节利用仿真参数和山西运城地区的8 通道机载阵列干涉SAR 真实数据进行实验验证，通过实验对比，验证了本文提出算法的有效性和优势；第4节对本文进行总结。

1 TomoSAR成像原理

1.1 TomoSAR成像模型

层析SAR 成像几何模型如图1 所示。图1 中，x为SAR 平台的航过方向，即方位向；r为距离向；s为斜距垂向；y为地距向；z垂直于x-y平面。θ为主影像的观测视角。SAR 平台经过L次航过，得到L幅单视复数（Single Look Complex,SLC）图像，经过复图像配准、去斜［21］以及相位误差补偿等处理后，影像的同一距离-方位单元构成一个长度为L的序列g={}g0,g2,…,gL-1，其中g0表示参考影像，gl(l=1,2,…,L-1)表示辅影像。根据文献［22］，沿斜距垂向上的散射点和影像gl(l=0,1,…,L-1)之间的关系表示为

图1 层析SAR成像几何模型

式中，γ(s)为观测目标沿斜距垂向的连续域后向散射，Δs为斜距垂向的后向散射信号的分布范围，为斜距垂向的空间观测频率，b⊥l为SAR 平台第l次航过时与主影像构成的垂直基线，λ为入射波长，r为斜距。层析SAR 成像就是已知观测序列g={}g0,g2,…,gL-1，反演后向散射系数γ(s)，从而将斜距垂向叠掩的各散射体分离并获知其散射强度以及位置的技术，实现SAR三维成像。

如果沿斜距垂向将式（1）进行离散化，则式（1）可以近似写为如式（2）所示矩阵形式:

式中，向量gL×1=[g0,g2,…,gL-1]T为对目标的L次观测数据，M是沿斜距垂向的离散化个数，A是一个L×M阶的测量矩阵，n=[n1,n2,…,nL]T为噪声向量，γM×1=[γ1,γ2,…,γM]T为待估计的后向散射系数向量。

1.2 基于波束形成估计理论的TomoSAR信号模型

多基线层析SAR 系统的各条航迹沿斜距垂向排列，各个航迹的天线中心可以看作是沿着斜距垂向构成的一个相控阵天线阵列，因此上述成像问题可以采用波达方向（Direction-of-Arrival,DOA）估计理论［23］来解决。

由式（2）所表示的线性方程组构成的SAR 层析成像信号模型，矢量γ由M个样本组成，这些样本是沿斜距垂向位置分别为的散射体的未知连续复反射率。矢量n为加性噪声。测量矩阵A由导向矢量组成，每个矢量维度为L，其中包含由位置给定的干涉相位信息，对于每一已知的sm，有

式中，kl=2πξl(l=1,2,…,L-1)，kl(l=1,2,…,L-1)还可以写为

式中，kl(l=1,2,…,L-1)代表参考卫星和第l个辅星之间的双向垂直波数。其中，dl(l=1,2,…,L-1)为z方向上参考卫星和第l个辅星之间的基线，r是目标到参考卫星的斜距，θ代表入射角度，λ是载波波长。

由于维度L通常远低于样本数M，式（2）中的线性方程有无穷多个解，这意味着未知数个数M比式（2）中的方程个数L多。此外，加性噪声和乘性噪声的存在增加了求解的不确定性。因此，为了确保对未知矢量γ求解唯一性，必须对式（2）中的测量矩阵A施加一些约束。

通常认为，后向散射矢量γ、噪声n和观测向量g为均值为零的复高斯随机过程，其自相关矩阵［17］可表示为

式中（·）+代表埃尔米特共轭，{E}是期望算子，因子N0是白噪声功率的功率谱密度。

在层析SAR 中，通常后向散射矢量γ是不相关的，因此相关矩阵Rγ为对角矩阵，则其主对角线上的元素所构成的矢量b为后向散射矢量γ的二阶统计量，即功率谱，该功率谱可用来表征层析SAR 三维成像斜距垂向上的后向散射。为求解出后向散射功率谱矢量b，可基于DOA估计理论的波束形成算法，利用导向矩阵A以及信号和噪声统计的先验信息［24］解决该问题。

2 TomoSAR三维成像方法

2.1 压缩感知层析SAR成像方法

基于压缩感知的层析SAR 成像主要解决的数学问题如式（8）所示:

式中AL×M为测量矩阵，由基线、斜距、波长等已知参数确定，观测所得的SAR 影像为gL×1，γM×1为待重建的稀疏信号，即斜距垂向的后向散射。

本文采用压缩感知中的贪婪算法开展层析SAR 三维成像的研究。正交匹配追踪算法［11］（Orthogonal Matching Pursuit,OMP）是贪婪算法中的经典算法。该算法采用迭代的方法得到待恢复的稀疏信号。算法初始化γ(0)=0，在第n次迭代时，用测量矩阵AL×M与残差r(n)=g-A·γ(n)进行施密特正交，然后找出前K（K表示后向散射的稀疏度）个最大元素对应的索引值添加到索引集Λ(n)。设第n次迭代满足结束条件，则恢复稀疏向量如式（9）所示:

OMP算法的伪代码迭代流程如表1所示。

表1 正交匹配追踪算法（OMP）伪代码

2.2 基于Capon波束形成谱估计层析SAR成像方法

谱估计算法中的波束形成技术获取层析SAR成像所需的后向散射系数，从噪声抑制的角度出发，通过加权处理，使信号输出噪声方差最小E[σ2]=h+mRghm的约束，计算获得最佳加权，进而实现对斜距垂向信号的估计。其设计原理是让滤波器对感兴趣方位的信号无失真地输出。本文基于式（7）研究了Capon 波束形成谱估计的算法，其最佳加权的约束方程为

该约束方程可由拉格朗日因子法求解，有

式（11）中的自相关矩阵Rg可由样本协方差矩阵代替，J代表观测数目。滤波器的输出功率构成的对角线矩阵为后向散射γ的自相关矩阵。由此可得出式（5）的Capon波束形成谱估计结果，具体的可写为

由于Capon滤波器要求数据协方差矩阵可逆，对角加载法是针对式（7）所示数据协方差矩阵Rg的秩不足时，提高Capon鲁棒性而普遍采用的一种方法［13］。一般来说，理想模型化的数据自相关矩阵Rg是不等于样本数据协方差矩阵G的。由于理想信号的波达方向和真实信号的波达方向之间存在的差异，当预期响应和实际响应之间存在差异时，Rγ和对角加载矩阵可以解释为转向矢量误差造成的差异［24］。

式中U和Γ通过特征分解G=UΓU+得到，U的列包含G的特征向量，对角矩阵Γ的主对角线聚集了各向量的特征值，因子对应于从下式获得的拉格朗日常数:

2.3 基于DCRCB谱估计的交替方向乘子法

对于如下的分布式凸优化问题［20］:

式中，f1(x)为损失函数，f2(y)为正则项，C，D，E为常系数矢量。对于对偶变量(x,y)，由于其由各自独立的函数约束，因此可以采用分解协调的形式，将小局部子问题的解通过协调以找到大局或全局问题的解。交替乘子法正是基于这样的思想而求解上述凸优化问题的。交替乘子法以增广的拉格朗日函数作为基础，融合双重分解和增广拉格朗日方法的优点，用于约束优化。具体的对于上述凸优化问题，交替乘子法构建如下的代价函数:

式中，μ为拉格朗日乘子向量，向量维度取决于线性方程组Cx+Dy=E的维度，c为交替乘子法的正则化因子。

论文引入交替方向乘子法解决Capon 谱估计的层析SAR 成像最优化问题。层析SAR 成像的分布式凸优化问题构建如下:

式中，ρ是正则化因子，通常取值为0.5～1.0，引入Capon 谱估计的层析SAR 成像结果作为中间变量，=(m=1,2,…,M)，由式（18）可以构建对于层析SAR 成像问题的交替方向乘子法（ADMM）优化函数:

通过迭代的方法即可估计出最优的成像结果。

2.4 交替方向乘子法的处理流程

交替方向乘子法通过交替迭代运算实现对式（18）所表示的优化问题的求解。具体如下:

式中，γk+1表示第k+1 次迭代时对γ估计结果，表示第k+1 次迭代时对重新估计的结果，重新估计可以有效抑制噪声，μk+1表示第k+1 次迭代时拉格朗日乘子向量更新，c为常数，通常可在0.03～0.05 范围内取值。γk+1与可通过对式（18）求偏导并将偏导置零来求解得到上述运算的解析表达式如下:

式中，I表示单位矩阵，S(·)(·)表示软门限函数，对于任意向量e与实标量门限ζ，有

式中，ui=Sζ(ei)，其中ei和ui分别表示向量e和u的第i个元素。

交替方向乘子算法通过式（20）进行迭代更新，当满足迭代结束条件

则输出=γk+1，其中ε为极小数为稀疏优化后的后向散射系数。具体的交替方向乘子算法流程如图2 所示。图2 中，首先进行初始化，然后通过式（20）进行迭代更新，最后输出最优化后向散射系数。

图2 交替方向乘子算法流程图

2.5 算法复杂度

压缩感知OMP 算法包括选择和更新两个环节。首先是索引集Jn的选择，即表1迭代过程的第一步，涉及排序算法，时间复杂度为O(KLM)；K表示稀疏度，L、M分别表示测量矩阵的维度。压缩感知OMP 算法第n次更新迭代的环节涉及到最小二乘法求解问题，时间复杂度为O(nN)。OMP 最多在2K次迭代终止，因此，OMP 总运行复杂度为O(KLM)。

本文提出的基于DCRCB 谱估计的交替方向乘子算法，在迭代过程中的运算量主要集中在更新估计向量γ时的矩阵求逆过程。从式（20）可以看出，当A和c固定时，待求逆的矩阵为一常数矩阵，这使得求逆运算可从循环中分离出来且仅为1次，其运算量为O(LM2)次浮点操作。可见两种算法的运算复杂度都与测量矩阵维度有关，在成像过程中，受基线和瑞利分辨率的限制，M的取值过大时，对成像结果的影响不大，因此两种算法的计算复杂度可视为同一量级。表2 列出涉及两种算法的时间复杂度。

表2 两种算法的时间复杂度

3 实验结果

为验证提出算法的成像结果，下面分别使用仿真数据和真实数据进行实验，其中真实数据采用中国科学院空天信息创新研究院在《雷达学报》官方网站公开的机载阵列干涉SAR 系统获取的山西运城地区的8通道数据［25］，进行大范围城区建筑物层析SAR 三维重建。实验环境的参数如表3所示。

表3 机载阵列干涉SAR系统参数

3.1 层析SAR仿真实验

利用表2所示的机载阵列干涉SAR系统参数，我们仿真了如图3 所示的非均匀基线分布时4 种算法的成像，结果如图4 所示。从图4（a）、（b）和（c）中我们不难看出，压缩感知算法和以Capon、DCRCB 为代表的波束形成算法在非均匀基线成像时对建筑物的顶、底以及墙体各面的还原效果都不佳，点云间断且有丢失现象，两种波束形成算法还有大量的噪点存在；而图4（d）所示的本文算法对建筑物形状具有较强的还原能力，点云形状完整且噪声抑制能力更优，同时成像结果更加清晰，这说明本文算法的层析SAR 成像方法在真实数据的应用场景中更加具有优势。

图3 非均匀基线分布

图4 仿真实验成像结果

3.2 层析SAR三维成像方法验证

运城地区的实验场景的Google Earth光学影像和SAR影像分别如图5（a）和（b）所示。

图5 实验场景的SAR图像及其对应的光学影像

为了验证算法，论文针对图5中红线所示的方位向，即影像的第2 000 行进行了层析SAR 三维成像反演。图6 给出了分别采用CS 算法、Capon 算法、DCRCB 算法和本文提出算法的成像结果，并对实验结果对比分析。

图6 运城地区真实数据SAR层析成像结果

从图6 可以看出，CS 算法将目标看作一组稀疏向量，需要提前设定稀疏度进行成像，但在实际场景中，强噪声可能也会替代目标，导致目标散射体丢失；也可能由于噪声的影像，会导致散射体的错位或偏移。Capon 算法、DCRCB 算法和本文算法是波束形成算法，从信号处理的角度出发，只单纯对所有散射点的回波信号进行滤波，本质上是从成像角度解决层析SAR 第三维的成像反演问题，并不会和CS 算法一样依赖于稀疏度。Capon算法由于其测量矩阵A所具有的正交性，功率谱也呈现较尖锐的形状，能够将旁瓣抑制得很低，但鲁棒性较差，在实际场景中极易收到干扰导致误差和噪声严重。DCRCB 算法在Capon 算法的基础上增强了鲁棒性，但同时就会带来分辨率的损失。本文提出的算法则进一步提高分辨率的同时，并降低旁瓣，可提升层析SAR 成像结果的分辨率及抗干扰能力。

图7 给出了第2 000 行第460 列的分辨单元的成像结果。可以看出，CS 算法和Capon 算法功率谱呈现较尖锐的形状，鲁棒性较差，但分辨率较高。而DCRCB 算法和本文算法增强了鲁棒性，本文算法相对DCRCB 算法降低旁瓣的同时更好地提高了分辨率。

图7 某一距离-方位单元成像结果（第2 000行第460列）功率强度

3.3 山西运城层析SAR三维成像结果

图8 展示了基于上述4 种算法对山西运城某地区的层析SAR 三维重建结果，其中，图8（a）和（b）所示的CS 算法和Capon 算法噪点过多，建筑立面分辨情况不佳，图8（c）所示的DCRCB 算法分辨率有所提升，但噪点过多现象依然存在，最后本文提出的基于交替方向乘子法的Capon 谱估计算法如图8（d）所示，在减少噪点的同时最大程度重建楼体轮廓，在不采用任何降噪手段处理时也能得到清晰成像结果。

图8 运城地区实验场景层析SAR三维重建结果

为了进一步解决噪声问题，本文采用DBSCAN方法提取有效点云［26］，利用点云的密度信息，提取其中密度较大的点云聚类。图9 给出了本文所提及的4 种算法经过DBSCAN 提取点云之后对整个建筑群的三维重建结果。从图9的处理结果来看，本文提出的基于交替方向乘子法的Capon 估计算法对后续的去噪方法没有强依赖关系，而其他3种成像方法必须经过去噪处理才能获得建筑的清晰成像结果。

图9 经DBSCAN降噪后的运城地区实验场景层析SAR三维重建结果

4 结束语

本文提出了一种基于交替方向乘子法的Capon 谱估计算法进行层析SAR 三维成像，并大规模重建了山西运城某城区的建筑物群，同时在噪点情况、分辨率、主瓣以及旁瓣等方面，与现有算法包括CS 算法、Capon 谱估计算法以及DCRCB 谱估计算法进行了SAR 层析成像结果的对比分析。通过实验，表明了基于传统CS 算法的层析SAR 成像在散射点选择问题上的一些弊端，而本文中提及的3 种基于波束形成的层析SAR 成像算法则从信号处理的角度，对所有散射点进行全面的反演工作，但这类算法需要目标信号的导向向量及噪声干扰的协方差矩阵，这一限制使其在真实场景中的性能大幅度降低，本文在对导向向量加以约束来提高其鲁棒性的基础上，通过L1范数，进一步提升此类算法的分辨率，减少噪声，获得更优成像结果，目的是突出信号处理算法的强还原能力。