1.   引言

最近，基于超材料孔径(Metamaterial Aperture, MA)的计算微波成像(Computational Microwave Imaging, CMI)系统以其使用便携的硬件设备进行信息检索而受到广泛关注^[1-8]。现有的微波成像体制主要包括合成孔径微波成像和实孔径微波成像。合成孔径微波成像利用目标场景与测量系统的相对运动形成一个电大尺寸虚孔径以获得高方位向分辨率，例如合成孔径雷达^[9]。实孔径成像系统是利用多个收发通道构成的真实物理孔径获得方位向分辨率，例如MIMO (Multiple Input Multiple Output)^[10]。不同于传统的实孔径微波成像方案利用阵列天线生成窄波束对场景进行扫描，MA-CMI可视为压缩感知微波成像，通过控制MA生成多种随机辐射模式构建测量矩阵对场景进行多次测量，其中测量矩阵的每一行都表示一个有效辐射模式^[4]。MA-CMI的最大优势在于它仅需要一对收发通道即可，大大降低了系统的硬件成本。根据压缩感知理论，当测量矩阵各列的相关性足够小时，近似保持约束等距条件，可从压缩测量数据中准确恢复稀疏信号^[11]。现有的MA主要分为两类：频率分集超材料孔径(Frequency Diversity Metamaterial Apertures, FDMAs)^[2,5-8]和动态超材料孔径(Dynamic Metamaterial Apertures, DMAs)^[3,4]。FDMAs可以产生高度独立的随频率变化的随机辐射模式，然而这需要FDMA支持较大的工作带宽以获得足够多的有效测量模式。DMAs通过独立控制超材料单元上加载的开关器件的通断状态，进行随机孔径编码获得大量随机辐射模式。基于DMAs的CMI系统不依赖于频率变化，因此其可以工作在单频点下，相比于FDMAs大大节省了频谱资源。

近年来，众多国内外学者对基于MA-CMI系统的成像算法展开了研究。由于MA-CMI可视为压缩感知成像，稀疏信号重构(Sparse Signal Reconstruction, SSR)算法，如正交匹配追踪(Orthogonal Matching Pursuit, OMP)^[5]、两步迭代收缩/阈值(Two-step Iterative Shrinkage/Thresholding, TwIST)算法^[6]或稀疏贝叶斯学习(Sparse Bayesian Learning, SBL)^[7,8]等被用于MA-CMI系统。SSR方法需要将场景区域离散化并假设目标散射中心位于离散网格点上。然而，在实际应用中这一假设通常不能被满足，场景离散化必然造成网格失配使得成像质量变差。针对这一问题，国内外已有许多基于SBL的离网格(Off-grid)误差校正方法。例如，对目标的转向矢量在其邻域网格上进行泰勒近似展开，在此基础上，利用稀疏贝叶斯推理^[12,13]、块稀疏贝叶斯推理^[14]或多项式求根^[15]等方法校正离网格误差。最近，文献[16]提出了一种基于信号子空间匹配的Off-grid波达方向(Direction Of Arrival, DOA)估计方法，文献[17]将该方法应用于MIMO雷达微波关联成像中。然而，在上述所有离网格误差校正的SBL方法中，都要求转向矢量可以用数学公式描述。对于MA-CMI系统生成的随机辐射模式无法用数学公式描述，且MA-CMI是一个单快拍SSR问题。使得上述所有方法都不适用于MA-CMI系统。

本文针对MA-CMI系统测量过程存在离网格误差的情况提出一种基于稀疏贝叶斯学习的Sinc插值离网格成像方法(Off-Grid based on Sinc Interpolation Sparse Bayesian Learning, OGSISBL)方法。首先，第2节针对MA-CMI系统的测量矩阵无法用数学公式表示的问题，提出新的离网格观测模型，其对不在初始划分网格上的散射体的回波数据用2维(Two-Dimensional, 2D)Sinc函数插值逼近在网格上的MA的实测辐射模式。此外，由于Sinc函数的非线性性质，对其执行1阶泰勒展开得到其线性近似式，使得算法的推导易于处理；第3节基于上述观测模型开发了一种具有离网格误差校正能力的成像算法。提出的成像算法基于SBL并在期望最大化(Expectation Maximization, EM)算法框架下迭代消除离网格误差；第4节通过仿真实验对本文提出的成像方法进行验证；第5节为结论。

2.   基于超材料孔径计算微波成像系统的观测模型

如第1节所述，MA-CMI系统通过生成多种随机辐射模式探测目标场景，同时接收天线收集场景内目标的回波。将目标场景在2D角域划分均匀网格，网格点坐标矩阵D可表示为

其中，ϕ和θ分别表示方位角和俯仰角，下标表示网格点位置。假设MA作为发射天线可以使用的有效辐射模式数为M，则采集在网格矩阵D上的MA与接收天线的第m种辐射模式可以分别表示为${\boldsymbol{E}}_m^{\text{E}}\left( {\boldsymbol{D}} \right)$和${\boldsymbol{E}}_m^{\text{R}}\left( {\boldsymbol{D}} \right)$。根据1阶玻恩(Bron)近似，测量矩阵与辐射模式间的关系可写为${{\boldsymbol{E}}_m} = {\boldsymbol{E}}_m^{\text{E}}\left( {\boldsymbol{D}} \right) \odot {\boldsymbol{E}}_m^{\text{R}}\left( {\boldsymbol{D}} \right)$，其中符号$\odot$表示Hadamard积。因此MA-CMI的第m次观测的测量模型可表示为

其中，${{\boldsymbol{h}}_m} = {\rm{vec}}\left( {{{\boldsymbol{E}}_m}} \right)$，vec(·)表示对矩阵矢量化算子，上标T表示转置操作，x代表目标场景的散射系数分布，w_m是第m次测量时的观测噪声。由此可知，MA-CMI的测量模型可表示为

其中，y=[y₁ y₂ ··· y_M]^T, H=[h₁ h₂ ··· h_M]^T, w=[w₁ w₂ ··· w_M]^T。

根据SSR理论，假设采样网格十分密集且目标落在采样网格点上，则可以从测量数据y中精确恢复x。然而，实际成像过程是无论网格多密集目标散射中心的位置总是不会落在这些采样网格点上，造成成像质量降低。为了解决这一问题，需要在观测模型中考虑离网格误差，并设计成像算法对其进行自适应估计和校正。根据衍射极限定理^[18]，由于MA的尺寸是有限的，可以将其辐射模式看作2维有限带宽信号，最大带宽为等效孔径。根据奈奎斯特(Nyquist)定理，在均匀采样网格上对测量模式进行2维Sinc函数插值，可以得到任意位置处辐射模式的幅值和相位。因此第m种测量模式在位置(${\tilde \phi _{{q_2},{q_1}}},{\tilde \theta _{{q_2},{q_1}}}$)处的值可以表示为

其中，${q_{\text{1}}} = {\text{1,}}2, \cdots, {P_{\text{1}}}， {q_{\text{2}}} = {\text{ 1,}} 2,\cdots, {P_{\text{2}}}，m{\text{=1,}} 2,\cdots {\text{,}}M$, ${\rm{sinc}}\left( x \right) = {{\sin\left( {{\pi}x} \right)}/{{\pi}x}}$, ${E_{m,{p_{\text{1}}},{p_{\text{2}}}}}$代表MA的第m种辐射模式在网格点(${\phi _{{p_1}}},{\theta _{{p_2}}}$)处的值。N=P₁×P₂表示划分网格点数量。由式(4)可知，任意网格处的第m种测量矩阵可写为

其中，${{\boldsymbol{B}}_\theta } \triangleq \left[ {{\boldsymbol{S}}_\theta ^1{\text{; }}{\boldsymbol{S}}_\theta ^2{\text{; }} \cdots ;{\boldsymbol{S }}_\theta ^{{P_1}}} \right] \in {\mathbb{R}^{N \times {P_2}}}$，${{\boldsymbol{B}}_\phi } \triangleq \left[ {{\boldsymbol{S}}_\phi ^1{\text{; }}{\boldsymbol{S}}_\phi ^2{\text{; }} \cdots ;{\boldsymbol{S }}_\phi ^{{P_2}}} \right] \in {\mathbb{R}^{N \times {P_1}}}$，K是一个行交换矩阵，作用是交换插值矩阵${{\boldsymbol{B}}_\theta } \otimes {{\boldsymbol{B}}_\phi }$的行顺序使得其前N行对应插值后的N个网格点的插值矩阵，[·]^1:N是取矩阵前N行矢量的操作算子，vecM(·)表示将列向量转化为矩阵的算子

S_ϕ和S_θ分别表示方位和俯仰维的插值矩阵，它们的上标对应坐标(i,j)处的网格点。

令${\boldsymbol{S}} = {\left[ {{\boldsymbol{K}}\left( {{{\boldsymbol{B}}_\theta } \otimes {{\boldsymbol{B}}_\phi }} \right)} \right]^{1:N}} \in {\mathbb{R}^{N \times N}}$，矢量化${\tilde {\boldsymbol{E}}_m}$可得插值更新后的第m种辐射模式为${\tilde {\boldsymbol{h}}_m}{\text{ = vec}}\left( {{{\tilde {\boldsymbol{E}}}_m}} \right) = {\boldsymbol{S}}{{\boldsymbol{h}}_m}$，由此可得插值更新后的测量矩阵为$\tilde {\boldsymbol{H}} = {\boldsymbol{H}}{{\boldsymbol{S}}^{\text{T}}}$，基于上述推导，MA-CMI系统的测量模型可写为

由于式(7)中包含非线性插值矩阵S，直接使用此模型校正离网格误差很困难。为此，本文使用S的线性泰勒展开近似式。令$\left( {{{\bar \phi }_{{q_2},{q_{\text{1}}}}},{{\bar \theta }_{{q_2},{q_{\text{1}}}}}} \right)$是距离位于$\left( {{{\tilde \phi }_{{q_2},{q_{\text{1}}}}},{{\tilde \theta }_{{q_2},{q_{\text{1}}}}}} \right)$处目标最近的网格点，则此散射体的离网格误差为

然后，S_ϕ和S_θ分别对$\Delta {\tilde \phi _{{q_2},{q_{\text{1}}}}}$和$\Delta {\tilde \theta _{{q_2},{q_1}}}$的1阶导数可表示为

其中，$\zeta \left( {{{\tilde \phi }_{i{\text{,}}{q_1}}} - {\phi _{{p_1}}}} \right)$是${\rm{sinc}}\left( {{{\tilde \phi }_{i{\text{,}}{q_1}}} - {\phi _{{p_1}}}} \right)$对$\Delta {\tilde \phi _{i{\text{,}}{q_1}}}$求1阶导数，具体过程如式(10)所示

同理，$\zeta \left( {{{\tilde \theta }_{{q_{\text{2}}}{\text{,}}j}} - {\theta _{{p_{\text{2}}}}}} \right)$是${\rm{sinc}}\left( {{{\tilde \theta }_{{q_{\text{2}}}{\text{,}}j}} - {\theta _{{p_{\text{2}}}}}} \right)$对$\Delta {\tilde \theta _{{q_{\text{2}}}{\text{,}}j}}$的1阶导数为

由克罗内克积(Kronecker product)求导法则可知，插值矩阵S的1阶泰勒展开近似式可写为

其中，${\dot {\boldsymbol{B}}_\phi } \triangleq \left[ {\dot {\boldsymbol{S}}_\phi ^1{\text{; }}\dot {\boldsymbol{S}}_\phi ^2{\text{; }} \cdots ;{\text{ }}\dot {\boldsymbol{S}}_\phi ^{{P_2}}} \right] \in {\mathbb{R}^{N \times {P_1}}}$, ${\dot {\boldsymbol{B}}_\theta } \triangleq \left[ {\dot {\boldsymbol{S}}_\theta ^1{\text{; }}\dot {\boldsymbol{S}}_\theta ^2{\text{; }} \cdots ;{\text{ }}\dot {\boldsymbol{S}}_\theta ^{{P_1}}} \right] \in {\mathbb{R}^{N \times {P_2}}}$, diagV(·)表示将矢量张成对角矩阵，${\boldsymbol{g}} = {\rm{vec} }\left( {\left[ {\begin{array}{*{20}{c}} {\Delta {{\tilde \phi }_{1,1}}}&{\Delta {{\tilde \phi }_{1,2}}}& \cdots &{\Delta {{\tilde \phi }_{1,{P_1}}}} \\ {\Delta {{\tilde \phi }_{2,1}}}&{\Delta {{\tilde \phi }_{2,2}}}& \cdots &{\Delta {{\tilde \phi }_{2,{P_1}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\Delta {{\tilde \phi }_{{P_2},1}}}&{\Delta {{\tilde \phi }_{{P_2},2}}}& \cdots &{\Delta {{\tilde \phi }_{{P_2},{P_1}}}} \end{array}} \right]} \right) \in {\mathbb{R}^{N \times 1}}$和${\boldsymbol{t}} = {\rm{vec}} \left( {\left[ {\begin{array}{*{20}{c}} {\Delta {{\tilde \theta }_{1,1}}}&{\Delta {{\tilde \theta }_{1,2}}}& \cdots & {\Delta {{\tilde \theta }_{1,{P_1}}}} \\ {\Delta {{\tilde \theta }_{2,1}}}&{\Delta {{\tilde \theta }_{2,2}}}& \cdots &{\Delta {{\tilde \theta }_{2,{P_1}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\Delta {{\tilde \theta }_{{P_2},1}}}&{\Delta {{\tilde \theta }_{{P_2},2}}}& \cdots &{\Delta {{\tilde \theta }_{{P_2},{P_1}}}} \end{array}} \right]} \right) \in {\mathbb{R}^{N \times 1}}$。另外，由Sinc函数的定义可知${\left. {\boldsymbol{S}} \right|_{{\boldsymbol{g}} = {{ {\textit{0}}}},{\boldsymbol{t}} = {{ {\textit{0}}}}}} = {{\boldsymbol{I}}_N}$。将式(12)替换式(7)中的插值矩阵S，可得最终的观测模型为

其中，${\boldsymbol{\varPhi}} = {\boldsymbol{H}} \left( {{{\boldsymbol{I}}_N} + \tilde {\boldsymbol{S}}_\phi ^{\text{T}}{\rm{diagV}}({\boldsymbol{g}}) + \tilde {\boldsymbol{S}}_\theta ^{\text{T}}{\rm{diagV}}\left( t \right)} \right)$，${\tilde {\boldsymbol{S}}_\phi } = {\left[ {{\boldsymbol{K}}\left( {{{\left. {{{\boldsymbol{B}}_\theta }} \right|}_{{\boldsymbol{t}} = { {\textit{0}}}}}} \right) \otimes \left( {{{\left. {{{\dot {\boldsymbol{B}}}_\phi }} \right|}_{{\boldsymbol{g}} = { {\textit{0}}}}}} \right)} \right]^{1:N}}$和${\tilde {\boldsymbol{S}}_\theta }{\text{ = }}\left[ {\boldsymbol{K}}\left( {{{\left. {{{\dot {\boldsymbol{B}}}_\theta }} \right|}_{{\boldsymbol{t}} = { {\textit{0}}}}}} \right) \otimes \left( {{{\left. {{{\boldsymbol{B}}_\phi }} \right|}_{{\boldsymbol{g}} = { {\textit{0}}}}}} \right) \right]^{1:N}$。基于此模型的成像算法在第3节给出。

3.   2D离网格稀疏贝叶斯学习

3.1   稀疏贝叶斯模型

已知散射体在场景中是稀疏的，本文采用层次贝叶斯模型^[19]建模场景中散射体的分布。其中第1层，假设待重构场景x服从零均值复高斯分布，并且x的第n行元素的精度为${\delta _n}$，则x的概率密度函数(Probability Density Function, PDF)如式(14)

其中，${\boldsymbol{A}} = {\rm{diag}}\left( {{\boldsymbol{\delta}} } \right)$, ${{\boldsymbol{\delta}} } = {\left[ {\begin{array}{*{20}{c}} {{\delta _1}}&{{\delta _2}}& \cdots &{{\delta _N}} \end{array}} \right]^{\text{T}}}$, CN$(\cdot)$表示复高斯PDF。层次贝叶斯模型的第2层，假设精度矢量${{\boldsymbol{\delta }}}$的每个元素${\delta _n}$, $n = 1,2, \cdots ,N$是独立同分布的且服从伽马分布(Gamma distribution)，其PDF如式(15)所示

其中，a是伽马分布的形状参数，b是率参数，$\varGamma (\cdot)$表示伽马函数。选择伽马分布做先验是因为它是高斯分布的共轭先验^[19]。另外，假设观测噪声w是零均值复高斯白噪声，由此可得观测矢量y的PDF为

其中精度参数$\beta$服从伽马分布

3.2   稀疏贝叶斯学习(SBL)

根据第2节得到的观测模型。对MA-CMI系统测量结果进行离网格误差校正可描述为如式(18)的最大后验(MAximum Posteriori, MAP)推断问题

其中，${{\boldsymbol{\varTheta}} } = \left\{ {{\boldsymbol{x}},\beta ,{{\boldsymbol{\delta}} }} \right\}$是隐藏变量，${{\boldsymbol{\varOmega}} } = \left\{ {{\boldsymbol{g}},{\boldsymbol{t}}} \right\}$是待估计参数。式(18)的后验分布如式(19)所示

其中，式(19)的分子项$p\left( {\left. {\boldsymbol{y}} \right|{\boldsymbol{x}},{{\beta}} ;{\boldsymbol{g}},{\boldsymbol{t}}} \right)$, $p\left( {\left. {\boldsymbol{x}} \right|{\boldsymbol{\delta}} } \right)$, $p\left( {\boldsymbol{\delta}} \right)$和$p\left( \beta \right)$分别由式(16)、式(14)、式(15)和式(17)给出，分母项可表示为$p\left( {\boldsymbol{y}} \right) = \displaystyle\iiint {p\left( {{\boldsymbol{x}},\beta ,{\boldsymbol{\delta}} ,{\boldsymbol{y}};{\boldsymbol{g}},{\boldsymbol{t}}} \right)} {\text{d}}{\boldsymbol{x}}{\text{d}}{\boldsymbol{\delta}} {\text{d}}\beta$。由于p(y)没有闭型表达式，使得不能直接使用EM算法进行迭代求解。这种情况下，采用变分贝叶斯期望最大化(Variational Bayesian Expectation Maximization, VBEM)方法^[19]。

3.2.1   变分贝叶斯-期望阶段

对于VBEM算法的第1步是计算后验PDF$p\left( {\boldsymbol{x}},\beta ,{{\boldsymbol{\delta}} } \left| {{\boldsymbol{y}};{\boldsymbol{g}},{\boldsymbol{t}}} \right. \right)$，同时计算隐藏变量${{\boldsymbol{\varTheta}} } = \left\{ {{\boldsymbol{x}},\beta ,{{\boldsymbol{\delta}} }} \right\}$的MAP。根据VBEM算法理论，假设各个隐藏变量${\boldsymbol{x}},\beta ,{{\boldsymbol{\delta}} }$是相互独立的，则后验PDF可近似写成如式(20)的形式

最小化式(20)中$\tilde q\left( {{\boldsymbol{x}},{\boldsymbol{\delta}} ,\beta } \right)$与$p\left( {\left. {{\boldsymbol{x}},{\boldsymbol{\delta}} ,\beta } \right|{\boldsymbol{y}};{\boldsymbol{g}},{\boldsymbol{t}}} \right)$的KL(Kullback-Leibler)散度得到后验PDF的最佳近似式^[19]。

根据此准则，可得$\tilde q\left( {\boldsymbol{x}} \right)$服从多元复高斯分布^[19]，即

其中，均值矢量${\boldsymbol{\mu}} {\text{ = }}\beta {\boldsymbol{\varSigma}} {{\boldsymbol{\varPhi}} ^{\text{H}}}{\boldsymbol{y}}$，方差矩阵${\boldsymbol{\varSigma}} = {\left[ {\beta {{\boldsymbol{\varPhi}} ^{\text{H}}}{\boldsymbol{\varPhi}} {\text{ + }}{\boldsymbol{A}}} \right]^{ - 1}}$。$\tilde q\left( {\boldsymbol{\delta}} \right)$服从伽马分布可写为^[19]

其中，形状参数$\tilde a = a + 1$和尺度参数${\tilde b_n} = b + \left\langle {x_n^2} \right\rangle$。$\tilde q\left( \beta \right)$服从伽马分布^[19]

其中，形状参数$\tilde c = c + 1$和尺寸参数$\tilde d = d + \left\langle {\left\| {{\boldsymbol{y}} - {\boldsymbol{\varPhi}} {\boldsymbol{x}}} \right\|_2^2} \right\rangle$。

由式(21)可得x的最大后验估计为

由式(22)可得${\delta _n}$的最大后验估计为

其中，算子(·)_(n,n)表示提取矩阵第n行，第n列元素。由式(23)可得$\beta$的最大后验估计为

3.2.2   变分贝叶斯-最大化阶段

根据VBEM算法理论^[19]，目标函数定义为对数似然函数的期望，即

其中，$\langle \cdot\rangle$表示求期望的算子。根据3.1节所示的贝叶斯公式，目标函数式(27)重新写为

其中，C表示常数。由式(13)可知，利用最大似然估计g和t，即Q(g,t)分别对其求偏导。

首先计算g的估计值，忽略Q(g,t)中与g^H无关的项，简化式(28)后可得Q(g,t)和$\partial {{Q\left( {{\boldsymbol{g}},{\boldsymbol{t}}} \right)} \mathord{\left/ {\vphantom {{Q\left( {g,t} \right)} {\partial {g^{\text{H}}}}}} \right. } {\partial {{\boldsymbol{g}}^{\text{H}}}}}$分别为

其中，diagM(·)表示取矩阵对角线元素构成列向量，(·)*为共轭算子，(·)^H是共轭转置算子。

对方程式$\partial {{Q\left( {{\boldsymbol{g}},{\boldsymbol{t}}} \right)} \mathord{\left/ {\vphantom {{Q\left( {g,t} \right)} {\partial {{\boldsymbol{g}}^{\text{H}}}}}} \right. } {\partial {{\boldsymbol{g}}^{\text{H}}}}} = { {\textit{0}}}$求解即可得到g的估计值$\hat {\boldsymbol{g}}$

其中，${{\boldsymbol{G}}_1} = \mathcal{R}\left\{ {\left( {{{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{H}}} \right) \odot \left( {{{\tilde {\boldsymbol{S}}}_\phi }\left( {{{\boldsymbol{\mu}} }{{{\boldsymbol{\mu}} }^{\text{H}}} + {{\boldsymbol{\varSigma}} }} \right)\tilde {\boldsymbol{S}}_\phi ^{\text{T}}} \right)} \right\}$, ${{{\boldsymbol{\gamma}} }_1} = \mathcal{R}\left\{ {{\rm{diagV}}\left( {{{\tilde {\boldsymbol{S}}}_\phi }{{{\boldsymbol{\mu}} }^*}} \right){{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{y}}} \right\}$, ${{{\boldsymbol{\eta}} }_1} = \mathcal{R}\left\{ {\rm{diagM}} \left[ {{{\tilde {\boldsymbol{S}}}_\phi }\left( {{{\boldsymbol{\mu}} }{{{\boldsymbol{\mu }}}^{\text{H}}} + {{\boldsymbol{\varSigma}} }} \right){{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{H}}} \right] \right\}$, ${{{\boldsymbol{\zeta}} }_1} = \mathcal{R}\left\{ {{\rm{diagM}}\left[ {{{\tilde {\boldsymbol{S}}}_\phi }\left( {{{\boldsymbol{\mu}} }{{{\boldsymbol{\mu}} }^{\text{H}}} + {{\boldsymbol{\varSigma}} }} \right)\tilde {\boldsymbol{S}}_\theta ^{\text{T}}{\rm{diagV}}\left( {\boldsymbol{t}} \right){{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{H}}} \right]} \right\}$，算子$\mathcal{R}\{\cdot\}$表示取实部操作。同理，求解${{\partial Q\left( {\boldsymbol{t}} \right)} \mathord{\left/ {\vphantom {{\partial Q\left( t \right)} {\partial {t^{\text{H}}}}}} \right. } {\partial {{\boldsymbol{t}}^{\text{H}}}}} = 0$即可得到t的估计值$\hat {\boldsymbol{t}}$。由于g和t在Q(g,t)中具有对偶的关系，这里直接给出${{\partial Q\left( {\boldsymbol{t}} \right)} \mathord{\left/ {\vphantom {{\partial Q\left( t \right)} {\partial {t^{\text{H}}}}}} \right. } {\partial {{\boldsymbol{t}}^{\text{H}}}}} = 0$的解，如式(31)

其中，${{\boldsymbol{G}}_2} = \mathcal{R}\left\{ {\left( {{{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{H}}} \right) \odot \left( {{{\tilde {\boldsymbol{S}}}_\theta }\left( {{\boldsymbol{\mu}} {{\boldsymbol{\mu}} ^{\text{H}}} + {\boldsymbol{\varSigma}} } \right)\tilde {\boldsymbol{S}}_\theta ^{\text{T}}} \right)} \right\}$, ${{\boldsymbol{\gamma}} _2} = \mathcal{R}\left\{ {{\rm{diagV}}\left( {{{\tilde {\boldsymbol{S}}}_\theta }{{\boldsymbol{\mu}} ^*}} \right){{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{y}}} \right\}$, ${{\boldsymbol{\eta}} _2} = \mathcal{R}\Bigr\{ {\rm{diagM}} \Bigr[ {{\tilde {\boldsymbol{S}}}_\theta }\left( {\boldsymbol{\mu}} {{\boldsymbol{\mu}} ^{\text{H}}} + {\boldsymbol{\varSigma}} \right){{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{H}} \Bigr] \Bigr\}$, ${{\boldsymbol{\zeta}} _2} = \mathcal{R}\left\{ {{\rm{diagM}}\left[ {{{\tilde {\boldsymbol{S}}}_\theta }\left( {{\boldsymbol{\mu}} {{\boldsymbol{\mu }}^{\text{H}}} + {\boldsymbol{\varSigma}} } \right)\tilde {\boldsymbol{S}}_\phi ^{\text{T}}{{\rm{diagV}}} \left( {\boldsymbol{g}} \right){{\boldsymbol{H}}^{\text{H}}}{\boldsymbol{H}}} \right]} \right\}$。

算法迭代终止判据由式(32)定义

或算法达到最大迭代次数，其中i表示当前迭代次数，$\tau$表示容差。本文设置$\tau = 1{0^{ - 4}}$，最大迭代次数为500。

本节最后，给出本文所提OGSISBL算法的流程图(如图1所示)并做简要说明：输入测量回波数据y和初始测量矩阵H，在初始化参数$\beta ,{\boldsymbol{\delta}} ,{\boldsymbol{g}},{\boldsymbol{t}}$和超参数集合$\left\{ {a,b,c,d} \right\}$，并由式(9)和式(12)分别计算定实矩阵$\{ {\dot {\boldsymbol{S}}_\phi },{\dot {\boldsymbol{S}}_\theta },{\boldsymbol{S}}\}$后进入SBL迭代计算。首先，将当前各参数值代入式(21)计算更新${\boldsymbol{\mu }}$和${\boldsymbol{\varSigma}}$；其次，分别由式(25)和式(26)更新参数${\boldsymbol{ \delta}} ,\beta$；最后，离网格误差矢量估计值$\widehat{{\boldsymbol{g}}}\text{ }和\text{ }\widehat{{\boldsymbol{t}}}$分别由式(30)和式(31)更新；重复这个迭代更新过程，直到满足收敛准则(参见式(32))或达到最大迭代次数。最终，输出重构场景$\hat {\boldsymbol{x}}$和更新后的网格坐标矢量$\hat {\tilde {\boldsymbol{\phi}}}$和$\hat {\tilde {\boldsymbol{\theta}}}$。

图 1  本文OGSISBL算法流程图

下载: 全尺寸图片 幻灯片

4.   实验仿真与分析

本节通过几个数值模拟实验来验证所提方法的性能。本节的所有成像实验使用的测量矩阵H是通过对一款FDMA-CMI系统在微波暗室环境下采集的辐射方向图数据进行处理得到的。FDMA-CMI系统由FDMA作为发射天线，喇叭天线作为接收天线组成(如图2所示)，此系统的详细设计在本文之前的工作中已给出^[20]。图3为34.4 GHz, 36.65 GHz和38.15 GHz频点激励下的FDMA辐射模式图。最终，本文采集得到104种有效辐射模式。文献[21]作为在网(On-grid)模型算法的代表与本文所提算法进行比较，本文中称该算法为On-grid SBL。为了揭示本文所提算法对场景稀疏目标的重构性能，采用均方根误差(Root Mean Square Error, RMSE)作为目标位置的重构性能度量，定义为

图 2  FDMA-CMI系统

下载: 全尺寸图片 幻灯片

图 3  FDMA在不同频点的归一化辐射模式

下载: 全尺寸图片 幻灯片

其中，${\hat {\tilde {\boldsymbol{\phi}}} ^{\left( i \right)}}$和${\hat {\tilde {\boldsymbol{\theta}}} ^{\left( i \right)}}$分别是第i次蒙特卡罗实验得到的散射体方位角和俯仰角的估计值，${\boldsymbol{\phi}}$和θ是散射体在角域上的真实位置，K是目标场景中的真实散射体数，N_C为蒙特卡罗实验次数。

实验1　假设一个位于空域xoy平面的目标场景，其中$x \in [ - 860{\text{ }}{\rm{mm}}{\text{, 860 }}{\rm{mm}}]$, $y \in [ - 860{\text{ }}{\rm{mm}}, {\text{ 860 }}{\rm{mm}}]$，FDMA-CMI距离场景500 mm。根据空域与角域的转换关系可知场景在角域的范围为ϕ∈[–60°, 60°]和θ∈[–60°,60°]。将目标场景在角域均匀离散化为31×31个网格点，即俯仰和方位维的网格间距r均为4°。信噪比(Signal to Noise Ratio, SNR)设置为10 dB。为充分证明所提算法对稀疏目标场景的重构能力，本文仿真了3种不同场景，分别是：五边形目标、五角星形目标和雪花形目标。图4显示了本文所提算法与On-grid SBL算法^[21]对3种不同场景的重构结果。图4(a1)、图4(b1)和图4(c1)分别为加入随机离网格误差作为扰动后的3种目标场景。其中用红色十字标记散射体的真实位置，而黄色网格标记了散射体所在的初始网格；图4(a4)、图4(b4)和图4(c4)分别为3种场景的真实散射体在空域的归一化幅度。图4(a2)和图4(a5)、图4(b2)和图4(b5)、图4(c2)和图4(c5)分别为3种目标场景的On-grid SBL算法重构结果。可以看出On-grid SBL算法没有更新网格，即该算法不能处理离网格误差，且其成像结果中存在很多高幅度假峰。其中最高假峰的幅度约为0.43，几乎为重构真实目标散射体最低幅度的66%，并且随着场景中目标数量的增加，高幅度假峰的数量也在增加，这在实际应用中是不允许的。本文所提算法对3种稀疏目标场景的重构结果分别如图4(a3)和图4(a6)、图4(b3)和图4(b6)、图4(c3)和图4(c6)所示。由图4(a6)、图4(b6)和图4(c6)可以看出所提方法的重构结果中尽管仍存在假峰，但最高假峰幅度约为0.15，是重构真实目标散射体最低幅度的20%，成像质量得到了很大提高。为了定量分析本文所提算法对离网格误差的校正效果，表1给出了3种稀疏目标场景200次蒙特卡罗实验不同SNR下消除离网格误差的能力。由表1可知，本文所提算法在SNR=20 dB时，3种目标场景的平均离网格误差均可以被减小到初始误差的52.6%左右；即使在SNR=10 dB的情况下，对于场景1的平均离网格误差也可以减小约初始网格的66.0%左右，但随着场景中稀疏目标个数增加算法校正离网格误差的能力略有下降。当SNR=5 dB时，对于3种目标场景，算法均无消除离网格误差的能力，这是由于此时噪声能量可以比拟目标能量。值得注意的是对于场景3，结合图4(c2)和图4(c5)可以看出，On-grid SBL算法失效而本文提出的算法仍然保持很好的重构结果(参见图4(c3)和图4(c6))。

图 4  3种不同稀疏目标场景，2种算法对于目标位置和归一化幅度的重构结果

下载: 全尺寸图片 幻灯片

表 1  比较不同信噪比下本文所提算法对散射体位置的重构效果(3种不同场景)(mm)

设置(场景1)	SNR=20 dB		SNR=15 dB		SNR=10 dB		SNR=5 dB
平均离网格误差	运行前	运行后	运行前	运行后	运行前	运行后	运行前	运行后
	11.00	5.79	11.00	6.13	11.00	7.26	11.00	14.78
设置(场景2)	SNR=20 dB		SNR=15 dB		SNR=10 dB		SNR=5 dB
平均离网格误差	运行前	运行后	运行前	运行后	运行前	运行后	运行前	运行后
	5.17	2.76	5.17	3.01	5.17	3.94	5.17	9.02
设置(场景3)	SNR=20 dB		SNR=15 dB		SNR=10 dB		SNR=5 dB
平均离网格误差	运行前	运行后	运行前	运行后	运行前	运行后	运行前	运行后
	5.38	2.86	5.38	3.18	5.38	4.17	5.38	11.09

下载: 导出CSV 

| 显示表格

实验2　本实验研究了2种算法对于不同离散网格间隔r随SNR的RMSE变化曲线以及对于不同SNR随离散网格间隔r的算法平均运行时间变化曲线。本实验采用五角星形场景(如图4场景2)，图5(a)为2种算法200次蒙特卡罗实验不同离散网格间隔下RMSE随SNR的变化曲线。从中可以看出2种算法的RMSE均随SNR的增大而减小，且r越小重构性能越好；在SNR为0 dB和5 dB时，无论r取多少，2种算法重构失败；当r一定时，本文所提算法性能始终优于On-grid SBL；当选取相对较粗的网格时，On-grid SBL算法随着SNR提升重构效果并没有明显改善，相比之下，本文提出算法的重构效果大大提升，这是由于本文方法考虑了测量过程中存在离网格误差，并在算法中对其进行了校正，且SNR越高，网格间隔越小，2种算法的重构效果差异越显著。图5(b)为2种算法在不同离散网格间隔r和SNR下的平均运行时间。从中可以看出对于不同SNR，2种算法得平均运行时间均随r的增大快速下降，这是由于当目标场景范围确定时，网格间隔大小直接决定了划分的网格数量，进而决定了计算数据的维度；对于不同网格间隔，On-grid SBL的运行时间始终低于本文所提算法，这是因为On-grid SBL不考虑离网格误差使得其计算量大大降低，但也造成其重构效果很差；当r一定时，2种算法的平均运行时间均对SNR变化不敏感，这说明了网格间隔是影响算法计算效率的主要因素。由图5可知本文所提方法在重构性能与网格间隔间存在反比关系，使得本文方法在实际应用中需要对其进行折中考虑以选取合适的网格间隔。

图 5  单快拍、不同信噪比和网格间隔的2种算法的重构性能

下载: 全尺寸图片 幻灯片

5.   结束语

本文首先针对MA-CMI系统的测量矩阵是一个伪随机矩阵，无法直接进行数学操作的情况，提出了一种基于2维Sinc插值函数的离网格观测模型。在稀疏贝叶斯学习框架下，基于此观测模型对离网格误差进行校正，获得较好的成像效果，并通过数值仿真验证了本算法的有效性。针对本成像算法的快速算法有待于进一步的研究。