1.   引言

层析SAR成像(Tomographic SAR, TomoSAR)是一种先进的SAR干涉测量技术，该技术能够监测并获取目标区域的高分辨率雷达影像，其中城市和森林区域是层析SAR成像技术的重要研究对象^[1,2]。不同于传统的合成孔径雷达(Synthetic Aperture Radar, SAR)只能进行二维成像，层析SAR成像将合成孔径雷达的合成孔径原理扩展用于三维成像的高程方向，能够重建散射体的三维信息并反演高程向剖面，其可以有效解决在二维SAR成像中处于同一散射单元内的目标散射点与雷达间斜距相等时存在的叠掩效应^[3-5]。相比于干涉合成孔径雷达(Interferometric SAR, InSAR)，层析SAR成像技术不仅可以获得目标散射体的高程信息，同时还可以获得散射体在高程向上的分布，能够完全恢复真实三维场景^[6]。基于以上的优点，TomoSAR成为现在最受关注且最具潜力的三维成像方法之一。

目前国际上主流的层析SAR成像方法主要包括多航过成像与阵列SAR两种，其体制不同但原理和信号处理方法是类似的。但多航过成像随着航过数的增加，成像的时间会大大增加，时效性较差；而阵列SAR系统的成本、重量与复杂度也会随着通道数的增加而急剧增加。因此，目前国内外的层析SAR研发方向都是着力于减少通道/航过数。德国宇航局(Deutsches zentrum für Luf- tund Raumfahrt, DLR)已经做到了3～5轨的层析SAR成像^[6]，还在继续向下减少通道/航过数。我国提出的“微波视觉”^[4]等新理论也在努力减少通道数。

层析SAR成像的信号处理方法研究始于20世纪80年代，在数十年的研究和发展中，其成像效果也在不断改善。谱估计方法是层析SAR成像的经典方法，包括CAPON算法^[7,8]、MUSIC算法^[9]等，但谱估计方法CAPON和MUSIC需要估计空间的协方差矩阵，而这个过程会降低方位-距离单元高程向的分辨率^[10,11]。稀疏信号处理方法，尤其是压缩感知(Compressed Sensing, CS)^[12-14]作为近年来的研究热点，也被应用到层析SAR成像中，该方法在相同基线数量下可以有效提高高程向的分辨率。文献[15]针对目标区域的小波域中的稀疏表征，同时结合相应多基线测量的二阶统计量，实现目标三维成像。基于压缩感知的层析SAR成像方法要求场景是点目标并且具有一定的稀疏性，压缩感知的重构算法主要包括贪婪算法和凸优化算法两类，贪婪算法一般计算速度较快但是在低信噪比和低采样数条件下重建精度较低，而凸优化算法的重建精度较高却计算效率较低。

城市和森林区域是层析SAR的重要应用对象。森林区域具有复杂的电磁散射特性，森林在高程向不具有典型的稀疏特性。同时城市区域中的建筑分布在高程向上屋顶部分通常具有稀疏特性，但地面上由于存在大量人工设施，例如花坛等，导致城市区域的地面不具有典型的稀疏特性。这为将稀疏信号处理应用于城市和森林区域提出了一定的挑战。

KL变换(Karhunen Loeve Transform, KLT)是一种数据衍生正交变换，它充分利用数据的相关性得到目标的自适应表征形式，在MRI等场景中已经得到了充分应用^[16]。应用KL变换处理一个信号，其离散展开式为

其中，$X(t)$为输入信号，${Z}_{n}={\displaystyle \int X(t){e}_{n}(t){\rm{d}}t}$为展开式的系数，${e_n}(t)$为信号的特征函数，特征函数需要根据输入信号来确定函数的类型。特征函数和信号自相关函数$E\left\{ {{X_k}{X_l}} \right\}$有以下关系：

此外在KL变换中，其展开式的系数${Z_n}$的范数为特征值${\lambda _n}$，表示相对应的特征函数${e_n}(t)$的能量值，且秩为非零特征值${\lambda _n}$的个数。同时，KL变换在均方误差准则下与其他正交变换相比效果更佳，是失真最小的一种变换，因此本文引入KL变换对目标的高程向进行处理，以获得更佳的成像结果。基于KL变换的重构问题可以转换为谱正则化的矩阵恢复问题，通过低秩矩阵的模型进行求解。由于目标区域在相邻方位-距离单元的高程向分布具有较强相关性，本文针对低航过数和低通道数的情况提出了一种基于稀疏和低秩结构的层析SAR成像方法，利用目标区域的低秩结构特性，引入KL变换对目标的方位-距离单元的高程向进行处理，表征其低秩结构特性，以此获得低航过数和低通道数情况下高分辨层析SAR成像结果。本文提出的基于KL变换的层析SAR成像方法不仅可以处理稀疏点目标场景，也可以处理分布式目标场景，其优化模型比传统压缩感知方法更为复杂，传统的压缩感知重构算法不再适用，本文采用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)优化算法，将复杂的原优化问题分解为若干相对简单的子问题，通过优化变量交替投影的方式进行算法求解，其计算复杂度与传统压缩感知在同一量级。

本文结构如下：第2节给出层析SAR成像的信号模型；第3节介绍基于稀疏和低秩结构的层析SAR成像方法；第4节利用仿真和实测数据的实验结果说明本文算法能够更好地重建三维目标场景；第5节对全文进行总结。

2.   层析SAR成像信号模型

层析SAR成像沿每一个轨道飞行获得一副聚焦的二维SAR复数图像，设${g_n}$是第n次获取的二维SAR图像，它可以理解为高程向散射单元的傅里叶变换的非均匀采样结果^[17]：

其中，$x$代表方位-距离二维平面位置，$s$代表高程向高度，$\gamma (x,s)$代表高程向后向散射系数，$\lambda$代表载波波长，$l$代表观测点到场景中心的斜距，${b}_{n}$代表第n条基线的高程向高度，$\varepsilon$代表噪声。通过对式(3)进行离散化处理，可以将式(3)记为${\boldsymbol{g}} =$${\boldsymbol{R}}(\gamma ) + \varepsilon$，其中${\boldsymbol{R}}$表示层析SAR成像的观测算子。

通过将场景的后向散射系数$\gamma (x,s)$进行重新排列，可以分析目标在相邻方位-距离单元的高程向分布的相关性，构造矩阵

其中，矩阵的列表示给定高度的方位-距离平面成像结果，也就是方位-距离平面二维成像结果按照方位向展开为一维向量，而矩阵的行表示给定方位-距离单元的高程向成像结果，矩阵的空间意义为三维网格每个格点的复散射信息。

以城市中的建筑为例，一个建筑物的外墙高度一致，所以一个高度的建筑物可以看作一个秩为1的成像结果矩阵，目标区域内建筑物的高度差异越复杂，则成像结果矩阵的秩越大。而对于森林区域，一种类型的树木其高程向分布基本一致，不同类型的树木其高程向分布则有一定差异，由此树木种类越丰富，则成像结果中的秩越大。在常见的城市、森林场景中，典型目标如建筑物、树木等，其种类较为一致，所以成像结果矩阵可预期具有先验低秩性。因为矩阵的行线性相关，所以矩阵${\boldsymbol{\varGamma}}$的秩$r \le \min (N,L)$。由矩阵的相关性，对矩阵${\boldsymbol{\varGamma}}$进行KL分解，可表示为以下3个矩阵的乘积

由此$\gamma (x,s)$可表示为r个高程基函数的加权线性组合

其中，${\rho }_{i}(x)$是矩阵${\boldsymbol{U\varSigma}}$的行，高程基函数${\nu_i}(s)$是矩阵${\boldsymbol{V}}$的列。

3.   基于ADMM的层析SAR成像方法

本文在稀疏信号处理的基础上利用相邻方位-距离单元的高程向的相关性，引入KL变换来表征其低秩结构特性，提高目标区域层析SAR成像能力。根据低秩矩阵的重构理论，矩阵${\boldsymbol{\varGamma}} \in {\mathbb{R}^{N \times L}}$的秩小于等于$\min (N,L)$，可通过其量测值${\boldsymbol{g}} = {\boldsymbol{R}}({\boldsymbol{\varGamma}} )$利用正则化方法进行高精度重建。由此形成如下稀疏与低秩结构相结合的层析SAR成像模型^[18,19]：

其中，${\boldsymbol{R}}$代表层析SAR成像的观测算子，$\psi ({\boldsymbol{\varGamma}} ) =$${\left\| {{\boldsymbol{\varGamma}} {\boldsymbol{\varPhi}} } \right\|_{{l_1}}}$, ${\boldsymbol{\varPhi}}$为小波变换矩阵^[15]，$\psi ({\boldsymbol{\varGamma}} )$表示矩阵中行向量小波系数的${l_1}$范数，用以表征矩阵${\boldsymbol{\varGamma}}$的稀疏特性，$\varphi ({\boldsymbol{\varGamma}} ) = {({\left\| {\boldsymbol{\varGamma}} \right\|_p})^p} = \displaystyle\sum\nolimits_{i = 1}^{\min \left\{ {N,L} \right\}} {\sigma _i^p}$表示矩阵${\boldsymbol{\varGamma}}$的Schatten-P范数，${\boldsymbol{\varGamma}} = {\boldsymbol{U\varSigma}} {{\boldsymbol{V}}^*}$表示矩阵${\boldsymbol{\varGamma}}$的KL分解， ${\boldsymbol{\varSigma}} = {\rm{diag}}\left( {\left[ {{\sigma _0},{\sigma _1},\cdots,{\sigma _{r - 1}}} \right]} \right)$, $\varphi ({\boldsymbol{\varGamma}} )$用以表征模矩阵低秩结构特性，${\lambda _1}$和${\lambda _2}$表示模型的正则化参数。本文方法和CS方法都需要利用目标场景的先验信息，对于CS而言，对于城市的建筑一般假设在高程向具有2～3个散射点，而森林树冠、城市地面等由于不具有典型稀疏性，其散射点个数往往较多，需要设置合适的矩阵秩。在理论上，根据稀疏信号处理理论与阵列信号处理理论，三维重建所需通道数/轨数最少也应大于目标的自由度，自由度在本文中体现为成像结果矩阵的秩^[6]。

由于式(7)同时使用到了稀疏和低秩结构，本文采用ADMM方法对其进行求解，式(7)可以转化为

其中，${\boldsymbol{K}},{\boldsymbol{S}}$都是辅助变量，${\beta }_{1},{\beta }_{2}$为正则化参数。当${\beta }_{1},{\beta }_{2}\to \infty$时，式(8)等价于式(7)的解。对于式(8)，本文使用如下的3步交替最小化方法进行求解

下面分别给出3个子问题的求解方法。对于第1个问题，由于它各项都是二次的，可使用以下方法对其求解：

其中，$Q\left( {\boldsymbol{\varGamma }} \right) = {\boldsymbol{\varGamma \varPhi}}$。

第2个子问题与标准核范数最小化问题类似，可将求解核范数最小化的迭代奇异阈值方法应用到式(10)中

其中，${{\boldsymbol{u}}_i},{{\boldsymbol{v}}_i},{\sigma _i}$分别是${{\boldsymbol{\varGamma}} _{L + 1}}$的奇异向量和奇异值，${Y}_{{\lambda }_{1}/{\beta }_{1}}$为奇异值阈值函数，阈值函数的定义如下所示

第3个子问题需要用到多维阈值的方法，具体计算如下：

通过不断迭代求解上述3个子问题，最终得到式(8)的解。算法流程如图1所示。

图  1  算法流程图

Figure  1.  Algorithm flowchart

下载: 全尺寸图片 幻灯片

4.   实验及结果分析

4.1   实验数据

为了说明基于稀疏和低秩结构的层析SAR三维成像算法的有效性，本文使用了仿真数据和机载实测数据进行实验，机载实验包含两组实测数据，一组为由德国宇航局(DLR)的机载F-SAR系统提供的HV极化通道的L波段数据^[20]，另一组为中科院空天院机载阵列干涉SAR系统获取的HH极化通道的Ku波段数据^[21]。其中L波段数据总航过数为9, Ku波段数据总通道数为12。

首先，仿真场景模拟了高程向0 m到90 m范围的森林区域，在该场景中，本文对不同航过数、稀疏度进行了模拟。图2(a)模拟了两个散射体，一个表示地面，另一个表示树冠散射中心，这两个散射体分别位于高程向9 m处和60 m处，对应的后向散射系数相等。图2(b)模拟了3个散射体，其中9 m处的表示地面，50 m和79 m处表示树冠散射中心，3个散射体的后向散射系数一致。仿真场景的相位服从$\left( {0,2\pi } \right)$的瑞利分布。仿真和L波段的F-SAR系统的实测数据所采用的雷达系统参数一致。雷达中心频率为1.325 GHz，方位向分辨率为0.4 m，距离向分辨率为1.5 m，实验分别采用9次航过数据和6次航过数据进行层析SAR成像。文献[20]给出了观测数据中的残余相位补偿文件和由基线间隔和中心斜距所计算的波数结果。

图  2  高程向归一化能量分布

Figure  2.  Normalized energy distribution in the elevation direction

下载: 全尺寸图片 幻灯片

其次，L波段F-SAR实测数据的采集区域为德国Trockenborn的部分区域，图3所示为该区域的光学图像，主要由森林和平地构成，非常适合验证本文方法在森林区域层析SAR成像的成像效果。Ku波段的实测数据为峨眉区域的城市数据，图4(a)为该区域的光学图像，图4(b)为该区域的SAR图像，该雷达系统的参数如表1所示。

图  3  F-SAR观测场景的光学图像

Figure  3.  Optical image of F-SAR measured scene

下载: 全尺寸图片 幻灯片

图  4  峨眉区域的图像

Figure  4.  Image of Emei area

下载: 全尺寸图片 幻灯片

表  1  Ku波段雷达系统参数

Table  1.  Ku-band radar system parameters

参数	数值
载波频率(GHz) 通道数N 成像场景海拔(m) 图像最近斜距(m) 航线海拔(m) 距离向像素尺寸(m) 方位向像素尺寸(m) 基线长度(m)	14.5 12
	420 1861 2157 0.1362 0.1051 2

下载: 导出CSV 

| 显示表格

4.2   仿真实验结果

仿真实验中，本文使用3种不同类型的方法对仿真数据进行重建，这3种方法分别是经典谱估计算法CAPON^[7,8]、压缩感知^[12-14]方法和基于稀疏和低秩结构的层析SAR成像方法。通过这些方法的仿真实验结果对比，以体现本文方法在层析SAR成像中的有效性。

首先在全航过(航过数为9)的情况下，分别采样以上方法对仿真数据进行重建。图5给出了航过数为9的情况下这3种方法的重建结果，可以看出CS和本文方法均能够准确恢复散射体的高程向位置，对归一化能量的值估计也大致准确，而CAPON算法对归一化能量值的估计存在误差，而且分辨率比另外两个方法要低。

图  5  航过数为9时各算法获得的高程向后向散射能量分布与仿真结果的对比

Figure  5.  Comparison of the simulation result of backscattered energy distribution of the elevation obtained by each algorithm with 9 interferograms

下载: 全尺寸图片 幻灯片

接下来将航过数减少为6，同样采样以上方法对仿真数据进行重建。图6给出了6条航过数的情况下这3种方法的重建结果，可以看出本文方法在减少航过数的情况下依然能够获得与航过数为9的情况下几乎一致的结果；而CS虽然可以恢复散射体的高程向位置，但归一化能量的估计出现误差，并且出现了虚假目标；CAPON算法恢复的散射体位置发生了偏移，无法准确估计目标的高程向位置。

图  6  航过数为6时各算法获得的高程向后向散射能量分布与仿真结果的对比

Figure  6.  Comparison of the simulation result of backscattered energy distribution of the elevation obtained by each algorithm with 6 interferograms

下载: 全尺寸图片 幻灯片

通过对图5和图6的4个图进行分析，可以发现本文方法在不同航过数和不同稀疏度条件下的结果均能够将多个散射体分离并进行重建，且对归一化能量的估计较好。

最后，本文针对两个散射体的场景分析噪声对航过数为9和航过数为6的仿真结果的影响，考察不同信噪比各种成像方法性能的对比。本文通过比较先验后向散射系数与重建的后向散射系数的均方误差(Mean Squared Error, MSE)说明算法的有效性，其具体计算公式如下：

其中，${y_i}$表示先验的后向散射系数真值，${{\hat y_i}}$为各方法重建后的后向散射系数。

实验结果如图7和图8所示。通过图7可知航过数为9的情况下，CAPON和其他方法相比误差较大，且不够稳定，而CS和本文方法随着信噪比的变化均有逐步减小的趋势，但本文方法拥有更低的均方误差，成像结果更加稳定。当航过数减少为6时，从图8可知本文方法的重建结果在不同信噪比下的均方误差均优于CAPON方法和CS方法。

图  7  航过数为9时不同信噪比情况下各算法的均方误差

Figure  7.  Mean square error of each algorithm at different SNR values with 9 interferograms

下载: 全尺寸图片 幻灯片

图  8  航过数为6时不同信噪比情况下各算法的均方误差

Figure  8.  Mean square error of each algorithm at different SNR values with 6 interferograms

下载: 全尺寸图片 幻灯片

4.3   机载实验结果

在机载实验中，本文仍然采用CAPON, CS和基于稀疏和低秩结构的层析SAR成像方法。首先对F-SAR系统森林数据进行层析成像处理。通过比较3种不同方法对实测森林区域数据的成像结果，由此验证本文方法在层析SAR的森林区域成像的性能优势。

在全航过数(航过数为9)的情况下，使用前文提到的3种方法对F-SAR的实测森林数据进行层析SAR成像，具体成像结果如图9所示。通过对比3种方法的成像结果，可以看出3种方法均能够将实测森林数据中的地面和森林区域分开，但CAPON方法与另外两种方法相比分辨率较低。

图  9  航过数为9时各算法获得的森林区域重建结果

Figure  9.  Reconstruction results of the forest area by each algorithm with 9 interferograms

下载: 全尺寸图片 幻灯片

随机抽取6次航过数据，依然使用以上3种方法在降低航过数的情况下对实测森林数据进行层析SAR成像。图10为航过数为6的情况下不同算法的成像结果。通过对比图10中各方法的成像结果，可以看出本文方法在减少航过数的情况下依然可得到高精度的成像结果，CAPON方法则出现了比较明显的分辨率下降，而CS对地面的成像效果出现地面不完整的现象。

图  10  航过数为6时各算法获得的森林区域重建结果

Figure  10.  Reconstruction results of the forest area by each algorithm with 6 interferograms

下载: 全尺寸图片 幻灯片

其次，我们将上面所说的3种方法应用于中科院空天院阵列干涉SAR系统获取的峨眉区域数据，对其进行层析成像处理。通过比较3种方法在峨眉区域的层析成像结果，说明本文方法在城市区域成像的有效性。

在全通道数(通道数为12)的情况下，各方法的成像结果如图11所示，可以看出CAPON方法对峨眉区域数据的建筑物进行层析成像结果较差，难以形成建筑外形，而CS和本文方法均能够区分出地面与建筑，且层析成像结果比较清晰。

图  11  通道数为12时各算法的峨眉区域重建结果

Figure  11.  Reconstruction results of Emei area by each algorithm with 12 channels

下载: 全尺寸图片 幻灯片

图12为通道数为8的情况下不同算法的成像结果。通过与图12对比可以看出CAPON方法对峨眉区域进行层析SAR成像的结果仍然较差，无法分辨建筑外形，CS虽然能够区分地面与建筑，但其成像结果中建筑顶部出现较多伪影，而本文方法在减少通道数的情况下依然能够区分峨眉区域的地面和建筑，并且基本不会出现伪影现象。

图  12  通道数为8时各算法的峨眉区域重建结果

Figure  12.  Reconstruction results of Emei area by each algorithm with 8 channels

下载: 全尺寸图片 幻灯片

本文对3种方法的成像结果进行距离向上的切片，并把3种方法的切片结果置于一张图上进行比较，如图13所示。

图  13  各方法成像结果的距离向切片对比

Figure  13.  Comparison of range slices of imaging results of each method

下载: 全尺寸图片 幻灯片

对图13进行分析，可以发现CAPON的成像结果中无法识别出建筑物；虽然CS方法能够对建筑和地面进行成像，但其在高度150 m附近出现了比较明显的伪影现象，如图中黑圈所示；而本文方法不仅对建筑和地面进行了高精度成像，并且基本没有出现伪影现象。

表2是各方法对峨眉数据进行成像的时间对比，实验运行平台：CPU为i9-10900X，内存为64GB，使用MATLAB软件计算。不难看出CAPON耗时非常少，本文方法耗时约为CS的1.67倍，基本在同一计算量级。

表  2  各方法耗时对比

Table  2.  Time comparison of each method

方法	耗时(s)
CAPON	61.844378
CS	43386.651744
本文方法	72769.757387

下载: 导出CSV 

| 显示表格

5.   结束语

本文提出了一种基于稀疏和低秩结构的层析SAR三维成像方法。本文从目标区域的特征入手，引入KL变换对相邻方位-距离单元高程向的低秩结构特性进行表征，构建稀疏和低秩结构相结合的层析SAR成像模型，采用ADMM方法推导层析SAR重构算法，获得高质量的三维重建结果。仿真与实测数据的实验结果表明，相比于传统的谱估计方法和稀疏信号处理方法，本文方法可以有效实现地面、树冠散射中心和建筑物成像，且在较少航过数和通道数时，保证目标区域三维重构精度。