面向高分辨SAR成像的正则参数自学习

刘 伟 焦卫东 廖仙华 杨 磊

（中国民航大学天津市智能信号与图像处理重点实验室，天津 300300）

1 引言

合成孔径雷达（Synthetic Aperture Radar，SAR）是一种全天时、全天候的高分辨率微波成像雷达［1］，与常规雷达类似，通过接收发射信号的目标回波进行一系列处理后实现目标探测。与常规雷达不同的是，SAR 是雷达通过运动形成虚拟合成孔径，通过相干回波处理实现高分辨成像。

针对SAR 成像过程中的特征提取与增强问题，传统傅里叶成像方法受限于奈奎斯特采样定理，成像主瓣受到带宽限制，高旁瓣导致目标特征难以提取，进而无法精确恢复目标散射特征。2004 年，由Donoho 与Candes 等人提出的压缩感知（Compressed Sensing，CS）［2］理论系统性地解决了这一问题。CS利用信号的稀疏特性，在远小于奈奎斯特采样率的条件下，用较少的观测值去恢复原始信号，可以实现由低维观测信号到高维目标信号的精确恢复。常见的压缩感知算法在面对高维SAR 信号时，运算复杂度将会大大升高，导致运算效率较低。2011年，Boyd等人完善了交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）［3］，将一个比较复杂的优化问题分解为几个较为简单的凸的子问题，同时结合近端算法，有效降低了运算规模。相比较其他优化算法，ADMM最突出的优势在于它的可分解性和收敛速率，降低了问题的求解复杂度，实现了算法的快速收敛。

在ADMM 框架下，为了便于高精度地恢复信号散射特征，我们通常引入先验假设，其中最典型的莫过于Lasso（Least Absolute Shrinkage and Selection Op⁃erator，Lasso）模型，ADMM框架将Lasso模型分解为平滑的最小二乘项和非平滑的ℓ1范数正则项，ℓ1范数用来调节待成像目标的稀疏度。因此，针对稀疏特征恢复［4］Lasso ADMM 具有良好的鲁棒性。但是，针对特征增强问题，不止需要稀疏特征的恢复，还有平滑特征和结构特征等等，所以基于ℓ1范数的ADMM恢复信号特征精度有限。因此，正则项的选取成为了特征增强的关键。2010 年，J.Friedman 提出稀疏组Lasso（Sparse Group Lasso，SG-Lasso）［5］模型，该模型在Lasso模型基础上加入了组ℓ1/ℓF混合范数，在恢复稀疏特征的前提下也保留了图像的结构特征。相比Lasso ADMM，SG-Lasso有效提升了SAR成像性能。

SG-Lasso 模型相比Lasso 模型来说，正则项中加入了组ℓ1/ℓF混合范数，减小了稀疏恢复过度所带来的负面影响，恢复了成像目标更加精细的结构特征，但是正则项参数的配比会影响信号的散射特征恢复，并且正则项的数目越多，正则项参数配比就越复杂。相比于单正则项，多正则项的参数设定更显得格外困难。假设单正则项的参数设定是在一维直线上确定一个参数点，那么设定多正则项参数是在二维平面乃至高维空间上确定一条线甚至一个面，设置和调整难度呈指数倍增长，这就导致多正则项多参数调节问题成为复杂的“炼丹”，单纯的手动调节在多参数下是难以实现的，常常需要依靠经验值去调节，并且时常伴随着难以收敛的问题存在。因此，针对单正则项参数乃至多正则项参数调节问题，本文提出一种参数自学习算法，将贝叶斯方法与ADMM将结合，不仅可以完成单正则项的参数评估，而且还可以实现多正则项多参数的自适应协同优化。

本文针对SAR 成像推导其回波信号模型，在该模型下针对单正则项参数和多正则项参数问题依次进行参数自学习。正则项参数的自适应估计提出了贝叶斯边缘估计（Marginal Estimation Bayes，MEB）算法，针对目标后验采样问题引入了马尔科夫蒙特卡洛（Markov Chain Monte Carlo，MCMC）采样，但是常规的MCMC 只能采样标准的后验概率密度函数，无法对ℓ1/ℓF进行采样，因此本文利用Moreau-Yoshida未经调整的朗之万算法（Moreau-Yoshida Unadjusted Langevin Algorithm，MYULA）［6］进行采样，最后将边缘似然估计与梯度投影法相结合，实现参数自学习。MEB 算法根据先验特征和似然分布得到先验概率和似然函数，通过二者相乘求积分得到边缘似然函数，但因边缘似然函数存在两项多维积分，求解解析解实在太过复杂，所以通过MYULA和梯度投影法在其对数梯度方向上迭代最优值，最终迭代求出参数最优值，实现参数自学习。本文实验选取点目标仿真与实测数据，使用MEB 算法进行正则项参数自学习，将最终所得结论与参数遍历最优值作对比，通过热力图去验证所提算法的恢复性能，并且与其他自学习算法进行比对，从而验证MEB算法的有效性与实用性。

2 回波信号模型

在机载SAR 成像过程中，地面成像场景保持静止，雷达随飞机按照理想航迹保持运动，通过雷达相对于目标的运动形成虚拟合成孔径，雷达在不同的位置进行回波接收并进行相干积累，进而实现高分辨SAR 成像。图1 为SAR 模式运行的几何模型，载机沿预定航线以速度v水平飞行，载机与地面静止目标场景的中心OS的参考距离矢量为R0，其中q0(t)为参考通道位置矢量，Rs为地面目标位置矢量，t表示方位向时间。由矢量几何关系式可知，在t时刻天线参考通道（Antenna Phase Center，APC）与运动目标的距离R(t)为

在SAR 成像几何模型中，雷达随载机保持运动，地面目标与场景处于静止状态，设目标与场景中心距离为Rp，载机与地面场景中心的参考斜距为|R0-q0(t)|，此时将目标斜距在参考斜距|R0-q0(t)|处一阶泰勒展开

其中Rp(xp，yp)为静止散射点P与目标场景中心OS的位置矢量，O(t)为泰勒展开的二阶项及更高阶项，其影响在远场假设下可以忽略不计。由式（2）可得场景下的累积SAR回波为

其中ky=为PFA 距离向插值，kx=为PFA 方位向插值。式（4）进行变量替换kx=-vt′k0/R0，其中t′表示插值后的方位向时间，与kx呈线性正比关系。将变量替换后的式子通过傅里叶变换可得SAR数据域回波为

由式（5）可以看出，回波数据是由地面上P个散射点的回波信号叠加而来，sin c 函数代表距离向，表示第P个散射点的距离向包络，指数函数代表方位向，表示线性相位项，λ为发射信号波长。

综上所述，SAR 模式的回波信号模型为线性回归模型如

其中Y代表经过预处理后的回波数据，在式（5）中表示为，通常为一个N×M的矩阵，其中N为方位向，M为距离向，X代表待恢复的信号，其矩阵元素对应方位向和距离向位置，其中代表距离向包络，exp(-jkxxp) 代表方位向线性相位。W代表加性噪声。A为方位向傅里叶字典，可以实现方位向多普勒成像。在SAR 模式下，A为傅里叶字典如

其中[•]T代表矩阵的转置，fd(•)表示方位向上的多普勒频率。

3 典型稀疏SAR成像问题模型

求解式（6）回波模型是典型的逆问题求解，但由于加性噪声的存在或现实存在的降采样需求使问题变得较为复杂，可利用压缩感知进行求解，需要引入目标X本身的散射特征先验知识，创建目标X的先验模型尤为重要。压缩感知有两类框架，一种是凸优化学习框架，另一种是贝叶斯统计学习框架。在凸优化学习框架中，目标先验信息通过惩罚函数引入，ℓ1范数可以针对目标稀疏先验进行建模，ℓ2范数可以针对目标平滑先验进行建模等等，这些范数的引入可以有效解决欠定的目标恢复问题。本文通过优化学习建立正则化框架去求解式（6），使其具备良好的适应性。在凸优化学习中，正则化线性回归［8］是最普遍的，其数学模型为

其中‖•‖1代表ℓ1范数，‖•‖F代表Frobenius范数，λ>0代表正则项参数。式（9）右边第一项为最小二乘项，也称保真项，用来表示恢复信号与原信号的相似度，第二项为正则项，也称惩罚项，用来控制恢复图像的稀疏度。在Lasso模型中，ℓ1范数正则项虽然可以很好的去除杂波和噪声，保留强散射特征，但是它也会影响弱散射特征的恢复，使得信号散射特征恢复精度有限。而SG-Lasso模型可以保留目标结构特征的前提下还能有效去除噪声与杂波，其数学模型为

其中Xi表示目标X的组结构，i表示分组数，范围从1 取到L。如果λ2=0，则表达式就与Lasso 模型一致；如果λ1=0，ℓ1范数正则项不存在，模型将变为组Lasso 模型，只能实现组间稀疏与组内元素平滑。式（10）右边第一项仍然为最小二乘项，第二项为ℓ1范数正则项，用来实现组结构内部元素的稀疏特征，而第三项正则项先用Frobenius范数实现组内元素平滑，然后求块结构的ℓ1范数，实现各组之间的稀疏表征，从而在恢复图像中既保证了有用信息不会被消除，又抑制了噪声与杂波。

综上所述，正则项引入了先验知识，实现了特征增强，所以正则项对信号恢复有着举足轻重的影响。正则项参数的大小决定着正则项的影响比重，但是传统的手动“炼丹”调参一直都是一个比较困难的问题，因为不同的成像方式、场景和噪声条件通常要设定不同比重的正则项，尤其面临多正则项多参数问题时，传统手动搜索的难度呈指数项增加，效率也会变得极为低下。

4 MEB参数自学习算法

正则化框架的建立可以较好地恢复目标的散射特征，而正则项参数的设定决定着先验模型和目标真实特征的匹配程度，进而影响信号的恢复效果。相比于传统的手动调参，MEB 算法在结合了梯度投影法后，对单正则项和多正则项的参数都可进行估计，并且针对非可微的正则项函数问题也引入了MYULA 进行解决，提升了ADMM 算法的恢复效率。首先根据贝叶斯原理推导贝叶斯边缘估计模型。

4.1 边缘估计贝叶斯建模

设连续凸函数f(X)代表式（9）或式（10）中最小二乘项，g(X)为式（9）或式（10）中正则项。根据传统贝叶斯理论，式（6）中Y的似然函数为

先验X的概率密度为

其中Q(λ)表示的是先验的归一化函数，表示为

在经典贝叶斯［9］方法下，正则项参数λ直接由观测数据Y来估计。例如在边缘似然估计中，可得到

当边缘似然函数取得最大值时，图像效果恢复最佳，λ也达到最佳值。其中边缘概率密度可由贝叶斯公式求得

由上述推导，虽然边缘概率密度可以由式（15）表示，但由式（13）及式（15）可知其涉及两个m 维积分，求解λ太过复杂。本文通过使用近似优化技术［10］，在成像逆问题中采用贝叶斯边缘估计算法实现正则化参数λ的自动选择。下面从单正则项和多正则项两方面进行参数自学习。

4.2 单正则项

对于单正则项而言，在Lasso 模型下，g(X)=λ‖vec(X)‖1为线性齐次函数，由于式（15）形式太过复杂难以求解，但针对其对数梯度

进行求解相对更为容易，其中l(X)=‖vec(X)‖1满足：对于任意X∈Rm，m 表示数据维度，l(tX)=tl(X)。因此可以推断出

对于所有的λ>0，将式（17）代入式（16）可得

采用随机近似梯度算法［11］，式（18）中的p(X|Y，λ)dX可以被一个马尔可夫采样器所代替，最后得到的梯度对数边缘似然为

由式（19）可知，通过对数梯度边缘似然近似迭代λ的条件为后验X的采样序列必须已知或者可以求解，这可通过马尔科夫蒙特卡洛（Markov Chain Monte Carlo，MCMC）［12］采样解决。MCMC 中未经调整的朗之万算法（Unadjusted Langevin Algorithm，ULA）［13］只能解决函数可微问题，而MYULA算法可以求解凸且不平滑的高维逆问题。而式（9）中目标函数

即为MCMC中后验概率密度的势能函数。

对于Lasso模型而言，最小二乘项为可微项，而正则项为不可微项，Lasso模型的代价函数为凸函数，并且是非平滑函数。针对此类非可微正则项求解问题，传统的ULA方法并不能对非可微函数进行梯度求解，因此本文采用MYULA算法进行求解，其迭代公式为

其中φ为步长，ε为平滑参数，N为d维均值为0的高斯随机变量。针对式（20）中势能函数的非可微部分，通过其邻近算子求次梯度［14］，可得到

对于任意的φ>0，此时式（21）变为

将式（23）代入式（19），可以求出边缘似然函数的对数梯度，则对于任意的n∈N，由梯度投影法有

其中∏Ω(•) 为投影算子，表示在凸闭集Ω 上的投影，δ>0为迭代步长。通过迭代在对数梯度边缘似然函数方向∇λlogp(Y|λ)上搜索，最终可得到参数λ的估计值。

4.3 多正则项

对于多正则项而言，存在多个正则项参数。本文考虑在双层稀疏组Lasso 模型下，因为g(X)=是可分离齐次函数，由两个线性齐次函数构成。在这种情况下，对于所有的λ∈Ω，先验的归一化函数Q(λ)为

Ap=dm×dp，dm和dp为X的维度。而混合ℓ1/ℓF范数和ℓ1范数都属于线性齐次函数，因此

利用MYULA 对多正则项情况下势能函数类似于式（23）进行采样，得到X的采样值，再通过梯度投影法对λ1和λ2进行并行迭代求解。

5 算法流程

本文算法结合贝叶斯理论，利用MYULA 采样并最终利用梯度投影法进行迭代求解参数。首先确定X0以及参数的初值，设置迭代次数，然后根据式（23）利用MYULA采样目标后验X，将采样得到的X带入式（24）进行迭代，将前25 次结果舍去，最后将所求结果求和取平均，得到的值即为参数自学习所得到的值。具体步骤如表1所示。

表1 贝叶斯边缘估计算法流程Tab.1 Algorithm Flowchart of MEB

6 实验验证

为验证所提参数自适应算法的有效性和准确性，本文首先基于目标点仿真在Lasso和SG-Lasso两种正则化框架下利用所提算法进行成像对比，并对所求参数与遍历参考值进行误差分析，遍历参考值是根据网格搜索的方式求得，具体求解方式为将参数先划分在一个大致区域内，然后通过网格搜索方式进行参数成像对比，成像恢复最好的那个参数即为遍历参考值，其中误差=。同时利用美国Sandia 实验室公布的复图像SAR 数据，在这两种正则化框架下进行成像比对。针对迭代次数，由于算法进行了大量实验，得出MEB 算法基本在25 次之后就达到收敛，所以本文将25 次前所求正则参数进行舍去，并将25 次之后（一般N取300即可）的正则参数进行平均值求解，以达到正则参数求解的精确性。最后将仿真模型的蒙特卡洛数据用于热力图成像恢复实验进行性能评估，用以定量分析算法的参数恢复精度。

6.1 仿真数据实验

首先仿真实验通过模拟“SAR”为成像场景，并在场景内随机分布散射点，模拟回波数据大小为256×512，其中仿真参数设置如下：中心频率为10 GHz，脉冲宽度为5 μs，发射信号带宽为150 MHz，脉冲重复频率为1000 Hz 以及作用距离为8 km。图2（a）表示无降采样无噪声下通过傅里叶变换得到的仿真参考恢复图像，图2（b）是降采样加噪声后通过傅里叶变化得到的成像结果，其中降采样率为0.7，信噪比为10 dB，图2（c）为Lasso 模型恢复结果，其中正则项参数由遍历获取，此时真值可认定为遍历结果，图2（d）为本文参数自学习成像结果。由上面分析可以看出，在降采样加噪声之后，通过成像恢复算法可看出本文算法估计的参数用于成像时噪声基本被抑制。由表2 所示，参数自学习数值为0.0038，与手动遍历参考值0.0035 相比，误差在10%以内。在效率方面，相比遍历参考值，算法复杂性得到了较大的降低。

表2 Lasso模型仿真数据参数误差分析Tab.2 Parameter Error Analysis of Lasso Model Simulation Data

另外，由图2（c）、（d）看出，在Lasso 框架下，图中黑色阴影基本被去除，代表信号噪声被抑制，但是仿真数据本身结构轮廓也没有完整恢复，代表一些弱散射点也被消除，只保留了信号的稀疏特征，结构特征无法完整保留。所以在此前提下，本文还在SG-Lasso 模型下利用算法对多正则项参数进行估计，通过实验来验证算法可行性。

本组仿真实验是通过随机模拟四个图形中的散射点而形成的，其大小为256×256，实验参数设置为：中心频率为10 GHz，脉冲宽度为10 μs，发射信号带宽为300 MHz，脉冲重复频率为1020 Hz。图3（a）是无降采样无噪声时傅里叶变换后的参考恢复图像，图3（b）是降采样加噪声后傅里叶变换的成像结果，降采样率为0.8，信噪比为10 dB，图3（c）表示遍历正则项参数的参考值在SG-Lasso 模型下的恢复图像，图3（d）为本文算法在SG-Lasso 模型下的恢复图像。

从图3 可以看出，恢复算法不仅保留了图像的稀疏特性，去除了（b）中的阴影部分，而且由于组ℓ1/ℓF范数的存在，图像的边缘轮廓与结构特征也被保留下来，恢复效果较好。而参数值的精确程度由表3 可知，参数自学习的估计值分别为0.008 和2.197，而在二维方向上手动遍历参考值为0.007和2.2，误差均不超过15%。在效率方面，参数自学习的复杂性相比遍历参考值相比有了明显的下降。

表3 SG-Lasso模型仿真数据参数误差分析Tab.3 Parameter Error Analysis of SG-Lasso Model Simulation Data

6.2 实测数据实验

本组实验采用美国Sandia 实验室的SAR 实测复数据实验图像，其大小为401×701。雷达作用于X 波段，发射信号带宽为1500 MHz，脉冲重复频率为2100 Hz，作用距离为10 km，距离向和方位向上的分辨率都为0.1 m。图4（a）表示未降采样未加噪的傅里叶变换结果，图4（b）为降采样加噪声的傅里叶变换结果，其中降采样率为0.7，信噪比为10 dB，图4（c）表示遍历正则项参数最优值在Lasso 模型下的恢复图像，图4（d）是经过本文参数算法在Lasso模型下的恢复图像。由图4 可以看出，降采样加噪声处理后，通过成像恢复算法可看出本文算法估计的参数用于成像时噪声基本被抑制。如表4 所示，参数自学习的估计值为0.184，与遍历最优值0.2相比，误差不超过10%。

表4 Lasso模型实测数据参数误差分析Tab.4 Parameter Error Analysis of Lasso Model Real Data

与仿真数据类似，Lasso模型对图像稀疏度恢复效果较好，但是图像中的有效信息也一并被稀疏掉了，导致了图像的结构特征丢失。因此本文也对SG-Lasso模型进行了参数自学习。

本组实验同样采用的是美国Sandia 实验室的SAR 实测复数据实验图像，场景为马路，其大小为401×321。雷达作用于X 波段，发射信号带宽为1500 MHz，脉冲重复频率为2100 Hz，作用距离为10 km，距离向和方位向上的分辨率都为0.1 m。图5（a）是未降采样未加噪的傅里叶变换结果，图5（b）表示降采样加噪声的傅里叶变换结果，其中降采样率为0.7，信噪比为10 dB，图5（c）为经过遍历最优值正则项参数在SG-Lasso 模型下的恢复图像，图5（d）表示经过参数估计算法在SG-Lasso 模型下的恢复图像。

从图5 可以看出，恢复算法不仅保留了图像的稀疏特性，而且由于组ℓ1/ℓF范数的存在，图像的结构特征也被保留下来，恢复效果较好。而参数值的精确程度由表5 可知，参数自学习的估计值分别为199 和1.66，与最优值200 和1.6 相比，误差在4%之内。

表5 SG-Lasso模型实测数据参数误差分析Tab.5 Parameter Error Analysis of SG-Lasso Model Real Data

6.3 热力图恢复实验

本节实验采用相变热力图（Phase Transition Diagram，PTD）［15］进行定量评估，该方法通过100 组蒙特卡洛实验，对每一个点求100次结果的平均值，可以计算在不同参数下算法成像的恢复结果与参考恢复图像的相关度。相关度越接近1，恢复效果越好，图像颜色也越深。

针对单正则项，实验利用ℓ1范数的参数计算恢复图像与参考图像的NMSE（Normalized Mean Square Error，NMSE）

其中Xk、X分别表示恢复图像和参考恢复信号。NMSE越小，图像恢复度越好。

根据图6 可以看出，当λ=0.0035 时，NMSE 最小，恢复效果最佳，与算法所得参数0.0038对比，误差在10%以内。

针对多正则项，以两个正则项参数为横纵坐标，通过100组蒙特卡洛实验验证算法的成功概率，如图7 所示，颜色从浅到深渐变表示恢复图像效果从差到好，其中降采样率为0.8，信噪比为10 dB。根据图7 可以看出，λ2对图像恢复起主导作用。当λ1为0.007，λ2为2.2 时，图像恢复效果最佳，与算法所得参数0.008 与2.197 相比，误差在15%以内，体现了算法的准确性和实用性。

6.4 目标恢复对比实验

为验证算法的实用性，本文引入传统联合后验估计算法［16］进行对比，传统联合后验估计主要是利用目标X与参数λ的互相迭代求解最优值，其迭代公式如下

其中t为迭代次数，从1 开始；n为目标X的维度，当n较大时，α和β对参数影响不大，一般都取为1，该算法一般在5～10次后达到收敛。

仿真实验是通过随机模拟两个图形中的散射点而形成的，其大小为256×512，实验参数设置为：中心频率为5 GHz，脉冲宽度为5 μs，发射信号带宽为300 MHz，脉冲重复频率为1000 Hz。图8（a）是无降采样无噪声时通过傅里叶变换后的参考恢复图像，图8（b）是降采样加噪声后的傅里叶成像结果，降采样率为0.7，信噪比为10 dB，图8（c）表示运用本文算法的恢复图像，图8（d）为运用联合后验估计算法的恢复图像。

根据图8 可以看出，相比于传统联合后验估计算法的恢复结果，本文所提算法的恢复效果更佳，并且传统联合后验估计算法只能对单正则项参数进行估计，局限性较大，而MEB 可对多正则项多参数进行协同自学习，体现了算法的有效性与实用性。

7 结论

本文针对优化问题中常见的正则项参数“炼丹”问题，引入了参数自学习算法，实现了单特征乃至多特征的参数自学习。贝叶斯边缘估计算法通过与贝叶斯原理相结合，对边缘似然函数进行建模求解，并采用MYULA 采样解决了先验分布非可微问题，最终使用边缘似然估计求出参数最优值，解决了传统手动搜索的繁琐冗长问题，使ADMM 算法更加高效，最后通过仿真数据和实测数据验证了算法的有效性和实用性。