1 引言 稀疏微波成像，是指将稀疏信号处理理论引入以合成孔径雷达(Synthetic Aperture Radar,SAR)为代表的微波成像中，将稀疏信号处理与微波成像理论两者相结合所形成的理论体系。稀疏微波成像的研究内容包括其基础理论、成像体制和实现方法等多方面内容。与传统SAR成像体制相比，稀疏微波成像不仅可降低系统复杂度，还能在目标分辨能力、模糊抑制、旁瓣抑制等方面提升系统的成像质量，具有广泛的应用前景，是下一代微波成像技术的重要发展方向[1, 2, 3, 4, 5, 6]。 机载平台是一个重要的稀疏微波成像系统部署平台。与传统的机载SAR系统类似，相对于星载平台，机载稀疏微波成像系统面临其独有的问题，其中最主要的，就是机载系统运动误差的处理。机载平台飞行轨迹受到运载平台性能、天气、机师操作等诸多因素的影响，很难保持严格的匀速直线运动轨迹，这将直接导致雷达回波信号与理想情况相比产生误差。运动误差一般表现为回波的相位误差，相位误差往往导致图像重建质量下降，出现散焦、位移等现象[7, 8]。随着相位误差的增大，甚至可能导致稀疏重建失败。在传统SAR成像中，相位误差的补偿主要有两种方法。一是利用运载平台的传感器，如惯性导航系统、GPS等，测量出载机平台的准确运动状态数据，从而计算出误差相位并进行相应的补偿。二是基于雷达回波数据的补偿，即自聚焦算法，直接从雷达回波数据中估算出相位误差并加以补偿[9, 10, 11]。常用的自聚焦算法有子孔径相关法(MD算法)、相位梯度自聚焦算法(PGA算法)等[9, 10, 11, 12]。虽然现有机载雷达系统的运动传感器性能越来越高，但自聚焦在高精度成像与提升现有系统性能等应用中仍有很大的价值[13]。 在稀疏微波成像系统的信号处理流程中，传统SAR成像的匹配滤波成像算法被取代，改用以${\ell _q}$正则化算法为代表的稀疏重建算法进行微波图像的成像与重建。与匹配滤波相比，稀疏重建算法具有低旁瓣、低采样率要求、可抑制模糊等诸多性能优势[1, 5, 14, 15, 16]，但由于稀疏重建算法的应用，针对匹配滤波发展出的传统自聚焦算法无法直接应用于稀疏微波成像系统中。现有针对稀疏微波成像自聚焦算法的研究集中在稀疏信号处理理论的框架下寻求对误差相位的估计与重建，其典型代表就是两步迭代算法[17, 18, 19, 20, 21, 22, 23]。该方法将稀疏微波成像视为一个非2次正则化框架下的一个优化问题，在迭代中同时建立误差项和稀疏重建场景两个迭代目标，分布交替迭代，在稀疏重建迭代的同时对误差相位进行纠正。该方法虽然在仿真数据中得到了成功的应用，但其本质上是基于迭代优化对相位误差进行求解，对误差模型不进行任何假设，运算量很大，难以在实际数据处理中得到应用。 另一方面，作为一种非线性算法，稀疏重建算法的运算量相比较于比配滤波算法亦较为巨大。在处理真实数据时，随着场景尺寸的增加，运算耗时急剧增加，使得稀疏重建算法在实际应用中受到限制。针对这个问题，文献[1, 24, 25]提出了基于回波模拟算子的快速稀疏微波成像算法。在该算法中，利用回波模拟算子及其逆算子的反复迭代实现正则化过程的加速。该算法实现原始数据稀疏重建时的2维解耦，可大大降低稀疏重建算法的运算复杂度，使得大场景稀疏重建成为可能。这也为我们提供了一个使稀疏微波成像自聚焦可行化思路，即将该快速算法的优势引入稀疏微波成像自聚焦中，实现自聚焦迭代的加速。 在本文中，我们提出一种结合传统自聚焦算法与回波模拟算子快速稀疏重建算法的稀疏微波成像自聚焦方法。回波模拟算子本质是上利用传统匹配滤波算法(如距离多普勒算法、Chrip Scaling算法等)及其逆过程代替传统稀疏重建中的观测矩阵以实现2维解耦与加速。因此，我们可以将传统自聚焦算法引入回波模拟算子的正向过程中，结合使得基于匹配滤波的自聚焦算法得以兼容于稀疏信号处理框架。本文以MD自聚焦算法为例，提出一种结合MD自聚焦算法与回波模拟算子的快速稀疏微波成像误差补偿算法。我们称之为MD-回波模拟算子自聚焦算法。 本文结构如下：首先在第2节建立稀疏微波成像系统的相位误差模型，并介绍现有基于两步迭代的稀疏微波成像自聚焦方法；第3节中，基于回波模拟算子的快速稀疏微波成像算法进行简单介绍；然后在第4节简单介绍MD算法，推导MD自聚焦算法与回波模拟算子的快速稀疏微波成像误差补偿算法；第5节将利用仿真和机载数据实验说明本文算法的有效性；最后将给出总结。 2 稀疏微波成像系统相位误差模型 2.1 机载稀疏微波成像相位误差模型 机载稀疏微波成像雷达的模型可以用式(1)表示[1]： ${y} = {Φ}{{x} + {N}}$ (1) 其中，y是回波信号时间采样后组成的列向量；x是观测场景的后向散射系数空间离散化后组成的列向量；${Φ} $称为系统的观测矩阵，其元素由雷达系统参数决定，为一卷积阵；N为回波中的加性噪声。相关向量与矩阵的具体定义在附录中给出。 按照稀疏信号处理理论，只要向量x稀疏，式(1)就可以进行稀疏重建，即从${Φ} $与y反求出x。不同于传统SAR的匹配滤波算法，稀疏微波成像中采用稀疏重建算法。常用的稀疏重建算法为${\ell _q}$正则化算法[26]： $ \mathop{\rm {min}}\limits _x\left\{ |\!| {{y} - {Φ} {x} |\!| _2^2 + \lambda |\!| {x} |\!| _q^q} \right\}$ (2) 其中，$0 < q \le 1$,$\lambda $是正则化系数，决定了重建过程中场景稀疏性的权重。 在实际机载系统中，系统不可避免面临平台非理想运动导致的运动误差，如图 1所示。 图 1 机载成像雷达运动误差示意图 Fig.1 Diagram of airborne SAR motion error 运动误差对于对波信号的影响主要体现在回波的相位误差[7]。在考虑相位误差后，式(1)所述观测方程中，每一个回波采样都将被附加一个误差相位： ${y} = {Θ} {Φ} {x} + {N}$ (3) 其中， ${Θ} = \left( {\begin{array}{*{20}{c}} {\exp ({\rm{j}}{\theta _1})}&{}&{}&{}\\ {}&{\exp ({\rm{j}}{\theta _2})}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{\exp ({\rm{j}}{\theta _M})} \end{array}} \right)$ (4) 而${\theta _1},{\theta _2},\cdots ,{\theta _M}$即为每个回波采样的误差相位，$M$为回波的采样总点数。 因为误差矩阵${Θ} $的存在，观测矩阵由${Φ} $变为${\tilde {Φ}} = {Θ} {Φ} $，若我们仍然使用${Φ} $作为观测矩阵进行稀疏重建，将导致观测矩阵与目标场景产生失配，导致成像性能下降，甚至稀疏重建失败。自聚焦的目的，就是对误差矩阵${Θ}$进行有效估计。 2.2 基于两步迭代的稀疏微波成像自聚焦算法 现有的稀疏微波成像自聚焦算法普遍基于两步迭代的方法[17]。该方法将稀疏微波成像的自聚焦问题视为一个多目标优化问题，将相位误差和稀疏场景同时视为优化项： $J({Θ} ,{x}) = \arg { \mathop {\rm {min}} \limits _{{Θ} ,{x}} }\left\{ |\!| {{y} - {Θ} {Φ} {x} |\!| _2^2 + \lambda |\!| {{x}_1}} \right\}$ (5) $J({Θ} ,{x})$就是优化的代价函数，表达为一个 ${\ell _1}$正则化的形式。$\lambda $是正则化参数，在这里表达了求解过程中稀疏性的权重。为最小化代价函数，传统的方法就是两步迭代法，即在每个迭代步骤中，首先固定相位误差项${Θ} $不动，更新场景的估计值x；再保持场景估计值x不动，更新相位误差项${Θ} $。因为每个迭代步骤中包含交替的两个子迭代，因此被称为两步迭代法。该算法的详细流程在下面给出。 基于两步迭代的稀疏微波成像自聚焦算法 Sparse microwave imaging autofocus algorithm based on two-step iteration 输入：回波信号y，观测矩阵${Φ} $ 输出：相位误差估计${{Θ} ^{(K)}}$，场景重建值${{x}^{(K)}}$ 1：设定迭代计数器$k = 1$；相位误差初值 ${{Θ} ^{(1)}} = {\bf{0}}$；场景 重建初值 ${{x}^{(1)}} = {\bf{0}}$ 2： ${{x}^{(k + 1)}} = \arg {{\rm{min}} _{x}}J({{Θ} ^{(k)}},{x})$ 3： ${{Θ} ^{(k + 1)}} = \arg {{\rm{min}} _{Θ} }J({Θ} ,{{x}^{(k + 1)}})$； 4：if 迭代不收敛 then 5： $k = k + 1$； 6： goto step2； 7：end if； 8：返回${{Θ} ^{(K)}},{{x}^{(K)}}$ 两步迭代算法完全通过数据对相位误差进行估计与补偿，是一种自聚焦算法。参考文献[17]用仿真数据说明了其有效性。但该算法基于正则化模型，由于每个迭代步骤需进行两次子迭代，运算量甚至还高于常用的基于${\ell _1}$正则化的稀疏重建算法。我们知道，稀疏重建算法的复杂度和观测矩阵的规模成平方关系，随着场景尺寸的增加，算法耗时会迅速增加。在实际应用中，问题规模往往会变的十分巨大，导致稀疏重建变得实际上不可能。例如，简单的估计表明，我们如果对尺寸为1024×8192的场景进行稀疏重建，观测矩阵的尺寸会达到107×107级别，运算消耗的内存空间可能达到数十TB级别，时间消耗在一般的工作站上会达到数月至半年。这当然是不可接受的。而基于两步迭代的自聚焦算法，其运算量至少是 ${\ell _1}$正则化算法的两倍。因此，在实际应用中，该方法基本不具有实施的可行性。 3 基于回波模拟算子的快速稀疏微波成像算法 基于回波模拟算子的快速稀疏微波成像算法对稀疏重建的过程进行加速。算法加速的主要思路就是将原始数据在方位向和距离向进行2维解耦，将大尺寸的观测矩阵按照方位向和距离向分解为两个小尺寸的观测矩阵[24]。 如式(2)所示的正则化问题一般通过阈值迭代(Iterative Soft-Thresholding,IST)方法求解，其迭代步骤为： ${x^{(k + 1)}} = f{\rm{(}}{x^{(k)}} + \mu {\Phi ^{\rm{H}}}\left( {y - \Phi {x^{(k)}}} \right),\lambda {\rm{)}}$ (6) 其中，$u$是迭代步长参数，${\{ \cdot \} ^{\rm{H}}}$代表共轭转置操作，$f( \cdot )$是软阈值操作： $ f(z,\lambda ) = \left\{ \!\! {\begin{array}{*{20}{c}} {(|z| - \lambda )\left( {\frac{z}{{|z|}}} \right)} {,\ |z| > \lambda }\\ \!\!\!\! \!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! 0,\quad {\rm{else}} \end{array}} \right.$ (7) 在式(6)所示的迭代步骤中，观测矩阵${Φ} $的尺寸正比于场景规模的平方。随着场景增大、采样点数增加，迭代中所需的乘加操作急剧增加。基于回波模拟算子的快速稀疏微波成像方法，将观测矩阵在方位向和距离向进行解耦，把大尺寸的矩阵分解为小尺寸的先行算子的组合，达到加速的效果。 回波模拟算子基于传统匹配滤波算法的操作算子进行构建，常见的匹配滤波算法如RD算法、CS算法等都可以用作2维解耦中的算子。我们以RD算法为例，用符号 ${P_{\rm{rd}}}\{ \cdot \} $表示RD算法的操作过程，称其为正向算子，则： $\hat {x} = {P_{\rm{rd}}}\{ {y}\} $ (8) 其中y是观测回波，$\hat {x}$是RD算法的成像结果——它可以视为真实场景x的一个有效近似。我们把${P_{{\rm{rd}}}}\{ \cdot \} $分解为3个线性算子的合成： ${P_{{\rm{rd}}}}\{ \cdot \} = {P_{\rm{a}}}\left\{ {C\left\{ {{P_{\rm{r}}}\{ \cdot \} } \right\}} \right\} $ (9) 其中 ${P_{\rm{r}}}\{ \cdot \} $表示距离压缩算子， $C\{ \cdot \} $代表距离徙动校正算子， ${P_{\rm{a}}}\{ \cdot \} $代表方位压缩算子。文献[24]证明了3个算子的线性型。而RD算子的逆过程，就称为回波模拟算子$T \,\{ \cdot \} $： $T\;\{ \cdot \} = P_{\rm{r}}^{\rm{H}}\left\{ {{C^{ - 1}}\left\{ {P_{\rm{a}}^{\rm{H}}\{ \cdot \} } \right\}} \right\} $ (10) 回波模拟算子可以视为观测矩阵${Φ} $的一种近似： ${y} = {Φ} {x} \hspace{12pt}$ (11) ${y} = T\;\{ \hat {x}\} $ (12) 这样，基于回波模拟算子的快速稀疏微波成像算法可用以下正则化形式表示： $\hat {x} = \arg \mathop {\min }\limits_{x} \left\{ |\!| {{y} - T \,\{ {x}\} |\!| _2^2 + \lambda |\!| {x}_q^q} \right\}$ (13) 其中$\lambda $是正则化参数。 该算法将大尺寸的观测矩阵${Φ} $近似地表达为3个小规模线性算子$P_{\rm{r}}^{\rm{H}}\{ \cdot \} $,${C^{ - 1}}\{ \cdot \} $和$P_{\rm{a}}^{\rm{H}}\{ \cdot \} $的合成，问题规模也得到了大大的降低。计算可以得到，该算法的时间复杂度由一般的IST算法的${\rm{O}}\left( {{n^2}} \right)$降至${\rm{O}}\left( {n\log \left( n \right)} \right)$，其中n代表场景规模。2维解耦的过程可简单用图 2说明。 图 2 基于回波模拟算子的快速稀疏微波成像算法示意图 Fig.2 Schematic diagram of accelerated sparse microwave imaging algorithm based on SAR raw data simulator 4 MD-回波模拟算子自聚焦算法 4.1 子孔径相关(MD)算法 Map-Drift(MD)算法[10]，又名子孔径相关法，是一种成熟的SAR成像自聚焦算法，本文不再加以赘述。MD算法基于有限阶的相位误差模型，对相位误差函数采用多项式拟合并估计误差系数，其往往需要经过多次迭代才能得到相位误差的完全补偿。 4.2 MD-回波模拟算子算法 由于MD算法是方位向的操作，因此，我们将单步的MD操作加入基于回波模拟算子的稀疏微波成像算法每一步迭代的方位向正向算子中。以基于RD算法的回波模拟算子为例，原基于回波模拟算子的稀疏微波成像算法的迭代过程为图 3所示，而MD-回波模拟算子算法的迭代过程则为图 4所示。 图 3 基于回波模拟算子的快速稀疏微波成像算法迭代步骤 Fig.3 Iteration step of accelerated sparse microwave imaging algorithm based on SAR raw data simulator 图 4 MD-回波模拟算子算法迭代步骤 Fig.4 Iteration step of MD-SAR raw data simulator algorithm 由于基于回波模拟算子的算法本身是正则化迭代算法，因此我们在每次迭代中加入的是单步MD操作。在每次迭代完成后，除需要判断正则化算法的收敛性，也需要同时判断MD过程的收敛性。 下面给出算法具体的迭代步骤。设输入为观测回波y，第$k - 1$步迭代的输出为${{x}^{(k - 1)}}$，那么本算法第$k$步迭代步骤如下： (1) 计算残差 IST算法在每次迭代中，首先需要计算图像的残差。记残差为$\hat {y}$： ${\hat {y}^{(k)}} = {y} - {{y}^{(k - 1)}}$ (14) 其中${{y}^{(k - 1)}}$是第$k - 1$步回波模拟算子输出的结果。 (2) RD正向算子操作 通过RD正向算子，从残差计算出本次迭代的输出。记RD正向算子输出为${x}_3^{(k)}$， $\begin{aligned}{x}_3^{(k)} & = {P_{{\rm{rd}}}}\{ {\hat {y}^{(k)}}\} \,= {P_{\rm{a}}}^{(k)}\left\{ {C\left\{ {{P_{\rm{r}}}\left\{ {{{\hat {y}}^{(k)}}} \right\}} \right\}} \right\} \\ & = \!{P_{\rm{a}}}^{(k)}\!\!\left\{\! {{x}_2^{(k)}} \!\right\}\! = {P_{\rm{a}}}^{(k)}\left\{ {C\left\{ {{x}_1^{(k)}} \right\}} \right\}\end{aligned}$ (15) 其中距离压缩算子为： ${P_{\rm{r}}}\{ \cdot \} = \exp \left( { - {\rm{j}}{{π}} /{K_{\rm{r}}}f_{\rm{r}}^2} \right)$ (16) 其输出为${x}_1^{(k)}$。RCMC算子是一个sinc插值核，难以用解析形式直接给出。其输出为${x}_2^{(k)}$；而这个插值可以表示为： ${x}_2^{(k)}({f_{\rm{a}}},\tau )\! = \!\!\sum\limits_\tau {{x}_1^{(k)}} ({\tilde f_{\rm{a}}},\tilde \tau ){\rm{sinc}}(\tilde \tau \! - \!(\tau \! + \!\Delta R({f_{\rm{a}}},\tau )))$ (17) 方位压缩算子为： $P_{\rm{a}}^{(k)}\left\{ \cdot \right\} = \exp \left( { - {\rm{j}}{{π}} /K_{\rm{a}}^{(k)}f_{\rm{a}}^{\rm{2}}} \right)$ (18) 其中${K_{\rm{r}}}$是距离向调频率，${f_{\rm{r}}}$是距离向信号变换到频域后的频率参数，${f_{\rm{a}}}$代表多普勒频率，$\tau $为快时间。 这里要强调的是，方位向调频率$K_{\rm{a}}^{(k)}$在每一次迭代中都是进行更新的；更新$K_{\rm{a}}^{(k)}$就是通过下面的单步MD操作。 (3) 单步MD操作 MD操作的目的是在每一次迭代中使用单步MD算法对方位向调频率进行一次更新。设${x}_2^{(k)}$为正向RCMC算子第$k$步的输出，将其在方位向划分子孔径，以两视为例，第1视和第2视分别记为： $\quad {x}_{2,1}^{(k)} = {x}_2^{(k)}({f_{\rm{a}}}:[- {B_{\rm{a}}}/2,0],\tau )$ (19) $\,{x}_{2,2}^{(k)} = {x}_2^{(k)}({f_{\rm{a}}}:[0,{B_{\rm{a}}}/2],\tau )$ (20) 其中${B_{\rm{a}}}$为方位向带宽。将前一次迭代中更新的方位压缩算子作用于两个子视并求取自相关，计算本次迭代中方位向调频率的偏差值： ${x}_{3,1}^{(k)} = P_{\rm{a}}^{(k - 1)} \Big\{ {x}_{2,1}^{(k)} \Big\}$ (21) ${x}_{3,2}^{(k)} = P_{\rm{a}}^{(k - 1)}\Big\{ {x}_{2,4}^{(k)}\Big\}$ (22) $\ {c^{(k)}} = {\rm{E}}\left[{{x}_{3,1}^{(k)},{x}_{3,4}^{(k)}} \right] \hspace{7pt}$ (23) $\quad \quad \Delta K_{\rm{a}}^{(k)} = \frac{{2K_{\rm{a}}^{(k - 1)} \ \ \Delta {x^{(k)}}}}{{vT}} \hspace{5pt}$ (24) 其中$\Delta {x^{(k)}}$是两个子孔径图像在方位向偏移量，由自相关${c^{(k)}}$求出；$T$是合成孔径时间，$v$是载机运动速度。 更新后方位向调频率为： $K_{\rm{a}}^{(k)} = K_{\rm{a}}^{(k - 1)} + \Delta K_{\rm{a}}^{(k)}$ (25) (4) 软阈值操作 软阈值操作是IST算法可以实现稀疏重建的关键步骤。通过软阈值操作，重建场景中的小目标被抑制，强目标被放大，从而实现稀疏约束。设阈值门限是$\lambda $，软阈值输出为： ${x^{(k)}} = f({\hat x^{(k)}},\lambda ) = \left\{ \matrix{ (\left| {{{\hat x}^{(k)}}} \right| - \lambda )\left( {{{{{\hat x}^{(k)}}} \over {\left| {{{\hat x}^{(k)}}} \ \ \right|}}} \right),\left| {{x^{(k)}}} \right| > \lambda \hfill \cr 0,\quad {\rm{else}} \hfill \cr} \right.$ (26) 其中 ${\hat {x}^{(k)}} = {{x}^{(k - 1)}} + \mu {x}_3^{(k)}$ (27) $u$是迭代步长参数。 (5) 迭代并进行回波模拟算子操作 判断输出${{x}^{(k)}}$是否收敛： $\!|\!|{{x}^{(k)}} - {{x}^{(k - 1)}|\!|} \lt {\epsilon_{x}}$ (28) 判断MD自聚焦是否收敛： $\Delta K_{\rm{a}}^{(k)} \lt {\epsilon_{K}}$ (29) 若未达到收敛目标，则计算其模拟回波，返回初始步骤重新开始迭代。 ${{y}^{(k)}} = T\left\{ {{{x}^{(k)}}} \right\}$ (30) 回波模拟算子$T \,\{ \cdot \} $定义如式(10)。 最终，MD-回波模拟算子算法可以用式(31)所示IST迭代求解框架表示： $\left. {\matrix{ {{x^{(k + 1)}} = f\left( {{x^{(k)}} + \mu P_{\rm{a}}^{(k + 1)}\left\{ {C\left\{ {{P_{\rm{r}}}\left\{ {y - T\left\{ {{x^{(k)}}} \right\}} \right\}} \right\}} \right\},\lambda } \right)} \hfill \cr {P_{\rm{a}}^{(k + 1)}\{ \cdot \} = \exp \left( { - {{{\rm{j}}\pi } \over {K_{\rm{a}}^{(k + 1)}}}f_{\rm{a}}^2} \right) = \exp \left( { - {{{\rm{j}}\pi } \over {K_{\rm{a}}^{(k)} + \Delta K_{\rm{a}}^{(k + 1)}}}f_{\rm{a}}^2} \right)} \hfill \cr {\Delta K_{\rm{a}}^{(k + 1)} = {{2K_{\rm{a}}^{(k)}\Delta {x^{(k + 1)}}} \over {vT}},\quad \Delta {x^{(k + 1)}} \sim {\rm{E}}\left[ {x_{3,1}^{(k + 1)},x_{3,2}^{(k + 1)}} \right]} \hfill \cr {x_{3,1}^{(k + 1)} = P_{\rm{a}}^{(k)}\left\{ {{f_{\rm{a}}}:[ - {B_{\rm{a}}}/2,0];C\left\{ {{P_{\rm{r}}}\left\{ {y - T\left\{ {{x^{(k)}}} \right\}} \right\}} \right\}} \right\}} \hfill \cr {x_{3,2}^{(k + 1)} = P_{\rm{a}}^{(k)}\left\{ {{f_{\rm{a}}}:[0,{B_{\rm{a}}}/2];C\left\{ {{P_{\rm{r}}}\left\{ {y - T\left\{ {{x^{(k)}}} \right\}} \right\}} \right\}} \right\}} \hfill \cr } } \right\}$ (31) 其中$\lambda $是正则化参数， $\mu $是迭代步长参数，$f( \cdot )$表示软阈值操作。 算法清单如下所示。 MD-回波模拟算子算法 MD-SAR raw data simulator algorithm 输入：带有相位误差的回波信号${{y}^{(1)}}$；距离压缩算子${P_{\rm{r}}}\{ \cdot \} $；初始带 误差的方位压缩算子$P_{\rm{a}}^{(0)}\{ \cdot \} $; RCMC算子$C\{ \cdot \} $；回波模拟 算子$T \,\{ \cdot \} \! =\! P_{\rm{r}}^{\rm{H}}\left\{ {{C^{ - 1}}\left\{ {P_{\rm{a}}^{\rm{H}}\{ \cdot \} } \right\}} \right\}$；迭代步长参数μ；阈值$\lambda $； 重建误差收敛门限值${\epsilon _x}$；多普勒调频率误差收敛门限值${\epsilon_K}$； 输出：重建的场景${{x} ^{(K)}}$ 1：初始化迭代计数器$k = 1$；初始化重建场景${{x} ^{(0)}} = {\bf{0}}$； 2：${\hat {y} ^{(k)}} = {y} - T\left\{ {{{x} ^{(k - 1)}}} \right\}$； 3：${x}_1^{(k)} = {P_{\rm{r}}}\{ {\hat {y}^{(k)}}\} $； 4：${x}_2^{(k)} = C\{ {x}_1^{(k)}\} $； 5：对${x}_2^{(k)}$划分子孔径，求取互相关，计算方位向调频率的误差 $\Delta K_{\rm{a}}^{(k)}$； 6：${x}_3^{(k)} = P_{\rm{a}}^{(k)}\{ {x}_2^{(k)}\} $； 7：${\hat x^{(k)}} = {x^{(k - 1)}} + \mu x_3^{(k)}$； 8：对${\hat {x}^{(k)}}$进行软阈值操作，得到${{x}^{(k)}} = f\left( {{{\hat {x}}^{(k)}},\lambda } \right)$； 9：if ${|\!|{x}^{(k)}} - {{x}^{(k - 1)}}|\!| < {\epsilon_x}$ and $\parallel \Delta \mathop k\nolimits_a^{(k)} \parallel < {\epsilon_K}$ 10： return ${{x}^{(k)}}$； 11：else 12： ${{y}^{(k)}} = T\left\{ {{{x}^{(k)}}} \right\}$； 13： $k = k + 1$； 14： 跳转到步骤2； 15：endif 算法流程图如图 5所示。 图 5 MD-回波模拟算子算法流程图 Fig.5 Flow chart of MD-SAR raw data simulator algorithm 5 仿真与实验 5.1 计算机仿真 本文采用点目标仿真数据验证MD-回波模拟算子算法的有效性，并和现有基于两步迭代的自聚焦算法进行对比。仿真参数如表 1所示。 表 1 仿真参数 Tab. 1 Simulation setup 实验场景中设置1个点目标，如图 6(a)所示，点目标的方位向切片在图 7(a)中给出。我们在方位向加入最大为0.15 ${{π}} $的2次相位误差，采用稀疏重建算法得到的含误差的重建结果如图 6(b)所示。特别地，其目标所在距离门的方位向切片在图 7(b)中给出。可以明显看到2次相位误差造成的散焦。在图 6(c)中，我们给出了采用本文所述算法进行自聚焦补偿后的稀疏重建结果，方位向切片在图 7(c)中给出。可以看到，本文所述算法较好地完成了2次相位误差的补偿。在图 6(d)和图 7(d)中，我们提供了两步迭代法的结果进行对比。 图 6 2维点目标仿真结果 Fig.6 Simulation of two-dimensional point target 图 7 2维点目标仿真方位向切片 Fig.7 Azimuth slice of simulation of two-dimensional point target 我们计算了点目标的峰值旁瓣比，在表 2中给出。本算法相比较于基于两步迭代的稀疏微波成像自聚焦算法，最大的优势在本算法基于回波模拟算子快速算法构建，运算效率较高。为说明本特点，我们亦比较本算法的运算耗时。在笔者使用的工作站上，运算1个1024×1024尺寸的场景，运算耗时也在表 2中给出。可以看到，因没有引入额外的迭代步骤，本算法对比于原基于回波模拟算子快速重建算法，并未引入太多额外耗时。相比于两步迭代算法优势巨大。 表 2 聚焦性能与运算耗时 Tab. 2 Focusing performance and time elapsed 5.2 机载数据实验 我们采用机载实验数据验证本文所述方法。机载实验参数如表 3所列[2]。我们利用图像熵值评估聚焦性能[27]： $S = - \sum {p_{i,j}}\ln {p_{i,j}}$ (32) 其中${p_{i,j}}$为图像中每点的能量密度。 表 3 机载实验参数 Tab. 3 Airborne experiment setup 实验结果在图 8中给出。图 8(a)为对存在相位误差的原始数据直接应用基于回波模拟算子的快速稀疏微波成像方法的成像结果；图 8(b)为应用本文自聚焦方法后的稀疏微波成像结果。可以看到，原数据直接成像存在较大的散焦现象。计算得到，此时图像熵值为16.91。而应用自聚焦方法后，成像质量有较大提升，基本实现聚焦，图像熵值也提升为12.93。 图 8 机载实验结果图 Fig.8 Airborne experiment result 6 结束语 本文首先介绍了稀疏微波成像模型和稀疏微波成像相位误差模型，简要介绍了现有基于两步迭代的稀疏微波成像自聚焦算法；然后介绍了基于回波模拟算子的快速稀疏微波成像算法；最后将MD自聚焦算法引入回波模拟算子框架中，构建了MD-回波模拟算子自聚焦算法。仿真和实验结果表明，本算法可以较好地实现稀疏微波成像中回波2次相位误差的补偿和自聚焦，收敛速度较快。 本文所述算法基于回波模拟算子和MD算法，因此只能处理2次相位误差。但在实际应用中，2次相位误差往往是对成像质量影响最大的一种误差形式。本算法利用了回波模拟算子进行2维解耦，聚焦过程未引入额外的迭代操作，在运算效率上具有较大优势，在实际机载数据的处理中具有广泛的应用前景。 附录 稀疏微波成像模型 设载机以速度$v$相对地面场景做匀速直线运动，根据SAR理论，对于一般的场景${Ω}$，雷达回波解调后的基带信号可以表示为： $\begin{aligned}s(t,\tau ) = & \int_{(a,b) \in {Ω} } \sigma (a,b)\exp \left\{ { - \frac{{4{\rm{j}} {π} {f_{\rm{c}}}R(t,a,b)}}{{\rm{c}}}} \right\}\\ & . \; p\left( {\tau - \frac{{2R(t,a,b)}}{{\rm{c}}}} \right){\rm{d}}a{\rm{d}}b + N\end{aligned}$ (A-1) 其中， $R(t,a,b) = \sqrt {{R_0}{{(a,b)}^2} + {v^2}{t^2}}$ (A-2) 是雷达天线相位中心到场景中的点$(a,b) \in {Ω} $的瞬时斜距；${R_0}(a,b)$是雷达天线相位中心到点$(a,b)$的近距；$\sigma (a,b)$是点$(a,b)$的后向散射系数；t是慢时间；$\tau $是快时间；c是光速；${f_{\rm{c}}}$是载波频率；$p(\tau )$是发射波形；N是回波中附加的加性噪声。 对上式进行时间和空间上的离散化，即可表示为一个卷积矩阵和场景后向散射系数相成的形式。我们记： $ \ \ a = {a_1},{a_2},\cdots ,{a_m},\cdots ,{a_M} \\ \; \ b = {b_1},{b_2},\cdots ,{b_n},\cdots ,{b_N} \\ \ \ \tau = {\tau _1},{\tau _2},\cdots ,{\tau _l},\cdots ,{\tau _L} \\ \ \ t = {t_1},{t_2},\cdots ,{t_k},\cdots ,{t_K}$ 这样，式(1)进行时间和空间采样后可表示为： $\begin{aligned}s[l,k] = & \mathop {\rm{\sum}} \limits_{{_{\Large{1 \le m \le M,1 \le n \le N}}}}\sigma ({a_m},{b_n}) \\ & . \; \exp \left\{ { - \frac{{4{\rm{j}}{π} {f_{\rm{c}}}R({t_k},{a_m},{b_n})}}{{\rm{c}}}} \right\} \\ & . \; p\left( {{\tau _l} - \frac{{2R({t_k},{a_m},{b_n})}}{{\rm{c}}}} \right) + N\end{aligned}$ (A-3) 并可写为矩阵形式： ${{y} = }{Φ}{{x} + {N}}$ (A-4) 其中， y = ${[s({t_1},{\tau _1})\,s({t_1},{\tau _2}) \,\cdots \,s({t_k},{\tau _l})\,\cdots \,s({t_K},{\tau _L})]^{\rm{T}}}$ (A-5) 是回波信号时间采样结果组成的向量；K和L分别是方位向和距离向的采样点数； x = ${[\sigma ({a_1},{b_1})\,\sigma ({a_1},{b_2}) \,\cdots \,\sigma ({a_m},{b_n}) \,\cdots\,({a_M},{b_N})]^{\rm{T}}}$ (A-6) 是场景后向散射系数空间采样后排成的向量； ${Φ} = \left( \!\! {\begin{array}{*{20}{c}} {{\phi _{1,1}}(1,1)} \ \cdots \ {{\phi _{M,N}}(1,1)}\\ {{\phi _{1,1}}(1,2)} \ \cdots \ {{\phi _{M,N}}(1,2)}\\ \vdots \quad \quad \quad \quad \quad \quad \vdots \\ {{\phi _{1,1}}(1,L)} \ \cdots \ {{\phi _{M,N}}(1,L)}\\ {{\phi _{1,1}}(2,1)} \ \cdots \ {{\phi _{M,N}}(2,1)}\\ \vdots \ \quad \quad \ddots \ \quad \quad \vdots \\ {{\phi _{1,1}}(K,L)} \ \cdots \ {{\phi _{M,N}}(K,L)} \end{array}} \!\! \right)$ (A-7) 为系统的观测矩阵，其中 $\begin{aligned}{\phi _{m,n}}(k,l) = & \exp \left\{ { - \frac{{4{\rm{j}}{π} {f_{\rm{c}}}R({t_k},{a_m},{b_n})}}{{\rm{c}}}} \right\}\\ & . \,p\left( {{\tau _l} - \frac{{2R({t_k},{a_m},{b_n})}}{{\rm{c}}}} \right)\end{aligned}$ (A-8)