1 引言

基于传统合成孔径雷达(Synthetic Aperture Radar, SAR) 2维成像，干涉合成孔径雷达(Interferometric SAR, InSAR)通过在高度向配置基线，利用多个天线相位中心在高度上的差异能够获取目标场景多视角观测，具备对观测场景3维测绘的能力^[1–3]。相应地，InSAR在数字高程测绘(Digital Elevation Model, DEM)生成、场景高度估计和地形变化检测等方面已经得到了广泛应用^[4,5]。在InSAR处理中，获取高精度的干涉相位是其中一项关键技术。然而，干涉相位不可避免的存在噪声是InSAR固有问题，这是因为通道之间视角差异引入SAR复图像之间的去相干问题，对应干涉相位噪声^[6]。在幅度方面，相干斑噪声是SAR图像典型性特征，这是由SAR成像的相干积累特性所决定的。直观地讲，相干斑噪声造成相同散射介质的物体在SAR图像幅度表现出“椒盐”噪声特性，从而影响目标的分类和识别。因此，InSAR干涉图像幅度和相位同时存在噪声问题。为了满足后续图像应用的要求，InSAR图像的相干斑抑制和相位噪声滤波处理是必要过程。

近十几年来，许多学者较为系统的研究了干涉相位噪声滤波和相干斑抑制算法^[7–10]。在现有大多数研究工作中，干涉相位噪声建模为加性噪声模型^[11,12]，而相干斑则基于乘性噪声模型假设。绝大多数算法都是基于场景目标局部平稳假设，即同一散射介质目标相邻散射系数具有平稳特性，并且目标高度上连续变化(除了边缘区域)。近几年来，越来越多的学者研究利用结合幅度和相位提高两者噪声抑制的性能。文献[13,14]提出了一种联合幅度和干涉相位的稀疏正则化算法，通过联合构造幅相差分(Total Variation, TV)字典进行联合稀疏约束，以提升在城区场景的降噪能力。在文献[15]中，作者的前期工作研究了利用自适应表征字典进行InSAR图像干涉相位和幅度的联合稀疏表征，在稀疏成像中能够实现幅相滤波的特征化增强。文献[16]利用InSAR干涉处理，构造稀疏字典进行InSAR图像稀疏化表征，实现在降数据率采样下的图像重构。

在稀疏信号处理中，贝叶斯算法作为一种典型性的算法理论，相比其他算法，具有更高的稳健性和精度。相应的，稀疏贝叶斯算法在雷达成像中得到了较为广泛的研究和应用。进一步，结构化稀疏贝叶斯成像被提出用于提升雷达成像性能。文献[17]研究了联合稀疏模型在稀疏ISAR(Inverse SAR, ISAR)成像中的应用。文献[18]研究了结构化统计模型在ISAR成像中的应用。在文献[19]中，作者的前期工作研究了结构化稀疏贝叶斯模型在ISAR成像中的应用。文献[17–19]都是利用ISAR图像中散射目标相邻结构相关性以提升成像性能，其区别在于统计建模和算法求解的差异。对于InSAR而言，通道之间具有相干性，利用其相干性，通过结构化贝叶斯方法联合通道图像处理，能够进一步提升干涉相位滤波和相干斑抑制的性能，对应本文的主要研究内容。

本文提出了一种基于贝叶斯联合通道稀疏的InSAR成像算法，通过贝叶斯结构化稀疏化表征，提升干涉相位噪声滤波和相干斑抑制性能。首先，基于贝叶斯准则建立InSAR成像模型，通过多层级统计建模，将InSAR复图像联合幅度和干涉相位，建立为结构化稀疏统计模型。然后，利用最大期望(EM, Expectation Maximization)算法进行InSAR图像重构和多层级统计参数估计。相比较文献[15]的方法，本文方法由于利用贝叶斯结构化稀疏模型，能够进一步提升InSAR图像幅度和相位噪声滤波的能力。最后，通过实验分析进一步验证本文算法的有效性。

2 InSAR信号模型

图1所示InSAR成像几何，在SAR成像基础上，InSAR通过在高度向增加通道数，即通道 $1, 2,·\!·\!· ,L$ (本文主要考虑两个通道情况，此时 $L{\rm{ = }}2$ )，可实现对目标场景不同高度视角的观测。InSAR与SAR 2维成像几何相同，平台飞行方向为方位向，对应图1中X轴，与X轴垂直的波束指向为距离向，对应图1中的Y轴。如图1所示，在方位时间 $t$ 时刻，通道 $l$ 与目标P的距离可以表示为 ${R_l}\left( t \right)$ 。假设已经利用SAR算法进行成像处理，此时数据在距离时域和方位多普勒域，第 $l$ ( $l = 1,2, ·\!·\!· ,L$ )个通道的信号可以表示为

式中， ${{{s}}_l}$ 为第l个通道距离压缩方位频域数据， ${{{a}}_l}$ 为第 $l$ 个通道的SAR复图像， ${{{T}}_l}$ 表示方位傅立叶变换矩阵，对应SAR投影算子(更为复杂的，若 ${{{s}}_l}$ 为回波数据， ${{{T}}_l}$ 对应距离和方位回波调制的SAR投影算子)， ${{\overset{\frown} {{n}}} _l}$ 对应第 $l$ 个通道的噪声，包含雷达系统噪声。此时，若直接对 ${{{s}}_l}$ 进行方位傅里叶变换即可获得第 $l$ 个通道的SAR图像。由前面引言中的分析可知，InSAR图像存在干涉相位噪声和幅度相干斑。此处，本文考虑在成像处理过程中，进行特征化增强处理以实现相位滤波，这种方式区别于一般的SAR图像后续处理方法^[15]。为了进行噪声抑制，本文将 ${{{a}}_l}$ 分解成两个部分：待成像重构的“干净无噪声”的SAR图像 ${{{b}}_l}$ 和噪声项 ${\breve{{{n}}} _l}$ 。此处，可以认为 ${{{b}}_l}$ 是降噪处理后的结果，其幅度和相位均是噪声滤波后的。相应的，式(1)可以改写成

式中， ${{{n}}_l}{\rm{ = }}{{{T}}_l}{\breve{{{n}}} _l}{\rm{ + }}{{\overset{\frown} {{n}}} _l}$ 表示总的噪声分量，包含InSAR幅相噪声。基于式(2)，本文的主要目的是从数据 ${{{s}}_l}$ 直接重构 ${{{b}}_l}$ ，实现干涉相位噪声和相干斑抑制。

3 所提方法 3.1 成像模型

在第2节中，式(2)实质是一个求逆以及降噪问题。在文献[15]中，作者的前期工作已经研究了利用自适应表征字典实现InSAR图像稀疏化表征，其稀疏成像模型为

式中， ${\left\| \cdot \right\|_1}$ 和 ${\left\| \cdot \right\|_2}$ 分别为向量1-范数和2-范数， ${\sigma _b}$ 和 $\sigma _{{n_l}}^2$ 为统计建模中图像 ${{{b}}_l}$ 和噪声 ${{{n}}_l}$ 的统计参数^[15]， ${\varPhi _l}$ 为自适应表征字典，以稀疏化表征第 $l$ 个通道的SAR图像。式(3)中， ${\varPhi _l} = \psi {{{P}}_l}$ 由两部分组成，其中 $\psi $ 为复数小波基， ${{{P}}_l}$ 为相位矩阵。对于两个通道而言， ${{{b}}_1} = \rm{diag} \left[ {\exp \left( {{\rm j} \cdot \arg \left( {{{{b}}_2}} \right)} \right)} \right]$ 和 ${{{P}}_2} = \rm{diag} \cdot\left[\exp \left( - {\rm j} \right.\right.$ $\left.\left.\cdot \arg \left( {{{{b}}_1}} \right) \right) \right]$ 。显然可见， ${\varPhi _l}$ 的目的是通过补偿另外一个通道的相位，进行干涉相位处理，然后利用小波变换对复图像，包含幅度和干涉相位，进行联合稀疏表征，这样可以克服单幅SAR图像相位随机的特性，同时，干涉相位处理便于后续的相位噪声滤波。其中，选择小波基作为稀疏字典的依据是利用InSAR图像幅度和相位的局部相关性，小波变换可实现图像压缩，更为详细的分析可参考文献[15]。通过求解式(3)可实现稀疏特征化成像，达到相位和幅度噪声滤波的效果。在实际中，InSAR图像之间是相干的，利用联合通道图像处理可利用此相干信息，提升稀疏特征化成像的性能，即

式中， $\left\| \cdot \right\|_{2,1}^{}$ 为矩阵的21-范数，相比于式(3)，式(4)能够利用通道之间的联合稀疏特性，这在稀疏信号处理中经常被用到。由于式(4)的成像模型，缺乏统计上普适的表示意义，本文拟从贝叶斯角度，通过建立联合多通道结构化稀疏模型，构建联合特征化成像算法。

通过贝叶斯建模，将 ${{{b}}_l}$ 建立为联合统计模型，其概率密度函数可以表示为

式中， ${{b}} = \left[ {\begin{array}{*{20}{c}} {{{{b}}_1}}& {{{{b}}_2}}&·\!·\!· &{{{{b}}_L}} \end{array}} \right]$ , ${{{Σ}}_b} \in {\mathbb{C}^{M \cdot N \cdot L \times }}^{M \cdot N \cdot L}$ 为 $\bar {{B}}$ 的协方差矩阵，用于表示 ${{{Φ}}_l}{{{b}}_l}$ 之间的相关性， $M$ 和 $N$ 分别为距离和方位采样数， $\det \left( {{{{Σ}}_b}} \right)$ 为 ${{{Σ}}_b}$ 的行列式。其中， ${{{Σ}}_b}$ 可以进一步表示为

式中， ${{{C}}_{mn}} \in {\mathbb{C}^{L \times L}}$ 为图像 ${{{Φ}}_l}{{b}}$ 第 $\left( {m,n} \right)$ 个像素的协方差矩阵。为了对通道之间相关性进行建模， ${{{C}}_{mn}}$ 表示为 ${{{C}}_{mn}} = {\gamma _{mn}} \cdot {{{R}}_{mn}}$ ，其中 ${\gamma _{mn}}$ 为一个非负的参数， ${{{R}}_{mn}}$ 表示相关矩阵。简单起见，令 ${{{R}}_{mn}} = {{R}}$ ，并且 ${{R}}$ 可以表示为

由式(6)和式(7)可知，参数 ${\gamma _{mn}}$ 用来控制通道像素之间的相关性大小， ${\gamma _{mn}}$ 越小表示不同通道像素之间相干性越差，越不相关。相应地，本文称式(5)—式(7)为结构化稀疏模型。在似然函数统计建模中，将 ${{{n}}_l}$ 建模为每个分量独立同分布的复高斯噪声模型，那么联合通道的似然函数可以表示为

式中， $\alpha $ 为噪声方差的倒数。

基于式(5)和式(8)，由贝叶斯准则求得最大后验概率密度函数

其中，

至此，基于贝叶斯建模，本文已经建立了联合通道稀疏成像模型。假设给定一组估计参数 ${{Θ}} \!=\! \left\{ {\alpha ,{{{Σ}}_b}} \right\}$ , ${{{b}}_l}$ 的最大后验概率估计为 ${\hat {{b}}_l}{\rm{ = }}{{{μ}}_l}$ 。在实际中，由于InSAR数据量通常较大，因此式(10)中的矩阵求逆运算量巨大，为了避免矩阵求逆， ${{{b}}_l}$ 可以利用求解线性方程求得

式(11)的问题可以通过共轭梯度算法进行求解，参考文献[20]，利用句柄操作直接对2维图像 ${{{s}}_l}$ 或者 ${{{b}}_l}$ 进行距离和方位分维度处理，而不拉成向量操作，能有效提高算法的实际可行性。贝叶斯方法利用EM算法将问题求解分解为图像 ${{{b}}_l}$ 重构和参数 ${{Θ}}$ 的估计两个部分，两个部分交替迭代，构成一个整体。

3.2 算法求解

为了得到参数 $\Theta $ 的估计值，通过可以利用最大似然估计方法求得，其中2次型估计算子[17]可以表示为

式中， ${\left\langle \cdot \right\rangle _{q(x)}}$ 表示对 $q\left( x \right)$ 的期望， ${{{Θ}}^{\left( {\rm old} \right)}}$ 表示前一迭代过程的估计值。在式(12)中， $Q\left( {{Θ}} \right)$ 可以分解为两个独立部分，分别为：

分别代表对 ${{{Σ}}_b}$ 和 $\alpha $ 的最大估计值。

(1) ${\gamma _{mn}}$ 和 $r$ 估计由式(5)、式(10)和式(13)，可以得到

设定 $Q\left( {{{{Σ}}_b}} \right)$ 在 ${{{C}}_{mn}}$ 处的导数为零， ${{{C}}_{mn}}$ 的估计算子表示为

式(15)中， ${{{Σ}}_{{b\ '}}}^{\!\!\!\!mn}$ 和 ${{{μ}}_{{b\ '}}}^{\!\!\!\!\!\!mn}$ 对应2维SAR图像第 $\left( {m,n} \right)$ 个像素联合通道矢量的协方差和均值。那么， ${\gamma _{mn}}$ 和 $r$ 的估计算子为

为了避免过度拟合过程，通过实验分析和经验获得， $r$ 的取值不大于0.25。

(2) $\alpha $ 估计通过式(8)、式(9)和式(12)，可以得到

设定 $Q\left( \alpha \right)$ 在 $\alpha $ 处的导数为零， $\alpha $ 的估计算子表示为

利用以上参数估计算法，用于估计式(10)中 ${{{{μ}}}_l}$ 和 ${{{{Σ}}}_l}$ ，通过不断迭代，直至收敛，最终可以获得稀疏成像结果，其算法流程框图如表1所示。需要指出的是，本文的方法通过联合通道贝叶斯结构化稀疏建模，能够提升稀疏特征化结果，对应更好的干涉相位噪声滤波和相干斑抑制结果。

运算量分析：通过分析可知，本文算法的主要运算量在于式(10)和式(11)中 ${{{μ}}_l}$ 和 ${{{Σ}}_l}$ 的计算，即矩阵的求逆操作。为了避免矩阵求逆，前面已经分析，可以利用式(11)通过共轭梯度算法求解。对于式(11)而言，其运算复杂度主要依赖于 ${{{T}}_l}$ 和 ${{{{T}}}_l}^{\rm{H}} $ 操作，分别对应于SAR投影算子以及其逆过程操作，在本文中为方位傅立叶变换矩阵及其逆变换矩阵，其运算复杂度为 ${\rm O}\left( {M \cdot N \cdot {{\log }_2}N} \right)$ 。假设EM算法迭代次数为 ${N_{\rm EM}}$ ，在每次迭代中，求解式(11)的共轭梯度迭代次数为 ${N_{\rm CG}}$ ，那么本文算法的运算量为 ${\rm O}\left( {L \cdot {N_{\rm EM}} \cdot {N_{\rm CG}} \cdot M \cdot N \cdot {{\log }_2}N} \right)$ 。其中，运算量与通道数 $L$ 成正比。

4 实验分析

下面通过实验分析以验证本文算法的有效性，并通过与现有方法比较，揭示本文算法的优势。在本节中，将实验分为2个部分。第1个实验部分为仿真数据实验，第2个部分为实测数据实验，包含山区数据和城区数据。下面所有实验均在个人电脑进行运行，电脑配置为CPU 3.50 GHz双核处理器，编程软件为R2015b版本MATLAB。

在仿真数据实验中，理想情况下图像幅度和干涉相位值如图2(a)和图2(b)所示。其中，幅度值对应平地和山区两个区域，每个区域对应相同的散射系数。图2(c)和图2(d)对应存在相干斑和相位噪声的图像。为了抑制噪声，利用本文算法进行稀疏化特征成像处理，其成像结果如图2(e)和图2(f)所示。可见，幅度的相干斑噪声和干涉相位噪声都得到了很好的抑制。为了更好地评估本文算法的性能，此处本文利用小波变换相位滤波^[21]算法进行比较，其结果如图2(g)所示。通过比较图2(e)和图2(g)可知，本文算法在地形平坦区域和高度陡变区域具有更好的相位滤波效果。图2(h)和图2(i)进一步给出两种算法结果与真实相位的差异图。通过比较可知，本文算法的相位估计误差更小。同时，计算了两种算法相位估计的均方误差，本文算法和文献[21]算法结果分别为0.0082 rad²和0.0405 rad²，与图2(h)和图2(i)一致。在运算量方面，EM算法的迭代次数设置为20，算法可近似达到收敛，此时两个通道图像处理时间约为118 s。

在实测数据实验中，我们利用RADARSAT-2重复航过干涉数据进行实验，每个通道为单视复(Single Look Complex, SLC)图像。两个通道的录取回波数据间隔时间为24 h，并被应用于DEM生成，可以用来高程估计。在第一个实验中，选取山区场景的部分数据。通道1幅度的图像如图3(a)所示，两个通道的原始干涉相位如图3(b)所示。明显可见，干涉相位存在噪声，SAR幅度图像存在相干斑。然后，利用稀疏特征化增强的方法，在稀疏成像中进行噪声滤波处理。其中，本文利用文献[15]中作者前期工作的算法作为对比。在稀疏化表征中，都选择双树小波变换作为稀疏表征字典。图3(d)–图3(i)分别为文献[15]和本文方法的成像结果。相比传统SAR成像，两种稀疏方法都能够较为有效的进行干涉相位噪声滤波和相干斑抑制。同时，我们给出了传统SAR成像、文献[15]和本文方法结果的相位“残点”，如图3(c)、图3(f)和图3(i)所示。可见，文献[15]和本文方法在山区地形陡变部分，能够保留一定的细节信息。相比于文献[15]，本文方法在地形平坦局域具有更好的滤波效果。需要指出的是，文献[15]和本文方法由于稀疏字典构造等局限性因素，其滤波效果具有一定的低通特性，这是后续工作需要改进的地方。

第2个实验部分，本文选取城区场景进行实验。文献[15]和本文的成像结果如图4所示。由图4可知，本文的方法在在建筑物区域具有更好的噪声抑制能力，同时在能够较好地避免相位的不连续性，相应地，干涉相位图像的突变区域对应幅度图像的突变部分，表现为观测场景建筑物的边缘区域。为了定性分析，本文通过计算图像的等效视数(Equivalent Number of Looks, ENL)评估相干斑抑制的性能，其定义为 $\rm{ENL} {\rm{ = }}\displaystyle\frac{{{\mu _{\rm enl}}^{\!\!\!\!\!\!\!\! 2}}}{{{\sigma _{enl}}^{\!\!\!\!\!\!\! 2}}}\ \ $ ，其中 ${\mu _{\rm enl}}$ 和 ${\sigma _{\rm enl}}^{\!\!\!\!\!\!\! 2}\ \ $ 分别为均匀场景区域强度图像的均值和方差。本文选取2个不同场景区域，以目标像素为中心，通过加窗(大小 $11 \times 11$ )选取相邻窗口像素计算ENL，其结果如表2所示。由表2可知，传统SAR成像的ENL结果近似为1，与使用的单视复图像一致。可见，本文的算法能够更加有效实现在同一散射介质区域的相干斑噪声抑制。相比于文献[15]的方法，本文算法可以取得更大的ENL值，对应较好的相干斑抑制性能，这得益于算法的结构化稀疏贝叶斯模型。为了评估算法对于独立散射点的分辨率保持能力，利用孤立散射点进行分辨率测量。传统SAR成像2维分辨率为2.67 m×4.13 m(理论值为1.92 m×3.27 m)。文献[15]和本文算法分辨率与传统成像比较接近，图5给出了方位脉冲响应形式，图中蓝色、红色和黑色分别对应传统成像、文献[15]和本文方法结果。

5 结论

本文针对InSAR稀疏成像，提出了一种通道联合的结构化稀疏贝叶斯成像算法。基于贝叶斯准则建立了联合通道结构化稀疏模型，并通过EM算法进行了图像重构和参数估计，相比于传统稀疏算法，本文算法利用了通道之间联合稀疏的相关特性，能够更加有效地提高InSAR成像算法中的噪声滤波特性，降低幅度噪声和相位噪声对成像结果的影响，实验结果证明了本文算法的有效性。