CN110426704B - 一种用于稀疏阵列的全变差快速成像算法 - Google Patents

Abstract

本发明涉及一种用于稀疏阵列的全变差快速成像算法，属于压缩感知成像技术领域。本发明提出一种单步Bi‑CGSTAB改进分割Bregman方法。分裂Bregman方法是一种基于交替方向乘法器ADMM的方法，用于解决TV正规化问题非常有效。首先，本发明改进了分裂Bregman方法的求解顺序，使得求解更高效；另一方面，引入Bi‑CGSTAB算法求解最优条件子问题，使算法快速收敛，并且在任何条件下，单步的Bi‑CGSTAB足以保证成像效果。分裂Bregman方法能够在几十次迭代之内以10^‑3的容差收敛到最终解。该方法可用于单输入多输出SIMO稀疏阵列系统的图像重建，也可不失一般性地扩展到MIMO阵列中。

Description

一种用于稀疏阵列的全变差快速成像算法

技术领域

本发明属于压缩感知成像技术领域，涉及一种用于稀疏阵列的全变差快速成像算法。

背景技术

毫米波成像广泛应用于人体安检系统中。从空间采样方面上，系统可分为单站扫描、线性阵列和多站(也称为多输入多输出，MIMO)三种模式。目前商用毫米波成像仪大多采用线性阵列扫描方式，在实时性和系统成本之间取得了较好的平衡。采用实孔径(RA)技术沿线性阵列排布方向获得分辨率，并将合成孔径(SA)技术应用于机械扫描方向。

为了进一步降低成像阵列的系统成本和复杂度，通常采用稀疏阵列技术来抑制阵元数目。这种技术的思想是在保证成像栅瓣和背景噪声水平可接受的前提下，尽量提高阵列稀疏度。即使如此，稀疏阵列仍然容易受到噪声影响而至今难以直接应用于成像系统中。目前比较常用的成像方法基于匹配滤波器原理，即获得自由空间中格林函数的伴随算子，如RA。该方法不需要迭代求解或矩阵求逆，因此即使不适用快速傅里叶变换(FFT)加速也能快速求解。

另一种成像质量更高但效率较低的方法是求解正则化逆散射问题，从压缩感知角度讲，范数能够很好地利用目标的点状特征。另外，全变差(TV)正则化是图像处理最常用和最成功的方法之一。该方法具有边缘恢复特点，对成片(二维)或成块(三维)目标的图像重建具有实用价值，并能够有效地抑制光栅瓣和背景噪声。然而，由于不可微性，该方法的计算量很大。目前提出的TV求解算法存在收敛速度慢，计算量大，内存占用率高的问题，因此很难应用于实时成像系统中。

发明内容

有鉴于此，本发明的目的在于提供一种用于稀疏阵列的全变差快速成像算法。

为达到上述目的，本发明提供如下技术方案：

一种用于稀疏阵列的全变差快速成像算法，首先，对分裂Bregman方法的求解顺序进行改进；其次，引入Bi-CGSTAB算法求解最优条件子问题，使算法快速收敛，并且在任何条件下，单步Bi-CGSTAB改进分割Bregman成像方法足以保证成像效果；最后，分裂Bregman方法在几十次迭代之内以10^-3的容差收敛到最终解；

在单输入多输出SIMO线性阵列中，接收回波表示为

S(x_r；k)＝∫∫σ(x，y)×exp[-jk(R_t+R_r)]dxdy (1)

其中，σ(x，y)是目标的反射率函数，k是波数；

设SIMO线性阵列沿X轴方向排布，散射场S(x_r；k)是接收阵元在(x_r，0)处接收的回波，S(x_r；k)视为关于x_r和k的二维矩阵；R_t和R_r分别是发射阵元和接收阵元与散射点(x，y)的单程距离；

将回波数据S(x_r；k)变形为单列数据

在成像域对反射率函数σ(x，y)进行网格剖分，并写成单列数据的形式/>

对于公式(1)，将积分算子表示为/>

公式(1)表示为矩阵乘法形式

Au＝f (4)

为提高图像质量，将TV正则化引入公式(4)的求解，即

根据Bregman方法，求解公式(5)等效于求解如下问题：

将图像重建问题转化为求解公式(6)的最小化问题。

可选的，所述单步Bi-CGSTAB改进分割Bregman成像方法，具体包括如下步骤：

S1：SIMO稀疏阵列数据采样；

在SIMO稀疏阵列中，发射阵元位于线性阵列的中心位置，接收阵元以任意方式排布于线性维度上；即在笛卡尔坐标系中，天线阵列与X轴重合；发射阵元将探测信号辐射到空间中，被目标体反射，由接收阵元接收回波数据，建立方程(1)并给出初始化条件；

S2：求解如下子问题

其中，i＝x和y；上角标(k)表示第k次迭代结果；

该子问题用收缩算子有效地求解，这是一个闭型解；

S3：进一步，采用单步Bi-CGSTAB算法求解如下子问题

S4：采用如下方程迭代更新变量

S5：将S2-S4构成迭代，直至达到收敛要求，即满足达到最大迭代次数或容差小于设定值。

可选的，所述S1的步骤为：

在SIMO阵列中，线性阵列沿X轴方向排布，散射场S(x_r；k)是接收阵元在(x_r，0)处接收的回波，k表示波数，并表示为

/>

其中，f是采样信号的频率分量，c是信号在空间中的传播速度；

初始化条件为u⁽⁰⁾＝A^Hf；

其中，上角标(0)表示初始化阶段；

所述S2的步骤为：

利用收缩算子求解最小化子问题(7)，得到

其中，i＝x和y；

和/>

分别表示X和Y维度上的方向导数；

所述S3的步骤为：

最小化公式(8)等效于求解如下最佳条件子问题

Ku^(k+1)＝rhs^(k) (11)

利用单步Bi-CGSTAB算法求解(11)。

本发明的有益效果在于：

1、将分裂Bregman方法改进为一种更合理的求解步骤，能够有效加快算法的收敛速度；

2、引入Bi-CGSTAB算法作为内部迭代求解最优条件子问题。Bi-CGSTAB算法具有收敛性速度快的特点，因此一步迭代Bi-CGSTAB足以达到分裂Bregman方法的收敛要求。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：

图1是SIMO阵列成像原理图；

图2是SIMO稀疏阵列采用RA的仿真成像结果；

图3是SIMO稀疏阵列采用本发明提出方法的仿真成像结果；

图4是SIMO稀疏阵列实验装置；

图5是SIMO稀疏阵列采用RA的实验成像结果；

图6是SIMO稀疏阵列采用本发明提出方法的实验成像结果；

图7是实验中容差与迭代次数的关系曲线。

图中：1、SIMO阵列的发射天线阵元，2、SIMO阵列的接收天线阵元，3、待成像目标，4、实验中的发射天线阵元，5、沿箭头方向的接收天线阵列，6、金属球目标。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

针对TV正则化问题求解计算量大的技术缺陷，本发明提出一种单步Bi-CGSTAB改进分割Bregman方法。分裂Bregman方法是一种基于交替方向乘法器(ADMM)的方法，用于解决TV正规化问题非常有效。首先，本发明改进了分裂Bregman方法的求解顺序，使得求解更高效；另一方面，引入Bi-CGSTAB算法求解最优条件子问题，使算法快速收敛，并且在任何条件下，单步的Bi-CGSTAB足以保证成像效果。分裂Bregman方法能够在几十次迭代之内以10^-3的容差收敛到最终解。该方法可用于单输入多输出(SIMO)稀疏阵列系统的图像重建，也可不失一般性地扩展到MIMO阵列中。

在SIMO阵列中，接收回波的可以表示为

S(x_r；k)＝∫∫σ(x，y)×exp[-jk(R_t+R_r)]dxdy (1)

其中，σ(x，y)是目标的反射率函数，k是波数；假设SIMO线性阵列沿X轴方向排布，散射场S(x_r；k)是接收阵元在(x_r，0)处接收的回波，因此S(x_r；k)可以视为关于x_r和k的二维矩阵。R_t和R_r分别是发射阵元和接收阵元与散射点(x，y)的单程距离。

首先，将回波数据S(x_r；k)变形为单列数据

在成像域对反射率函数σ(x，y)进行网格剖分，并写成单列数据的形式/>

对于公式(1)，将积分算子表示为/>

由此公式(1)可以表示为矩阵乘法形式/>

Au＝f (4)

为了提高图像质量，将TV正则化引入公式(4)的求解，即

根据Bregman方法，求解公式(5)等效于求解如下问题：

因此，我们将图像重建问题转化为求解公式(6)的最小化问题。

本发明的目的就是为了解决上述问题，提供：

单步Bi-CGSTAB改进分割Bregman成像方法，包括如下步骤：

步骤(1)：SIMO稀疏阵列数据采样：在SIMO稀疏阵列中，发射阵元位于线性阵列的中心位置，接收阵元以任意方式排布于线性维度上。即在笛卡尔坐标系中，天线阵列与X轴重合。发射阵元将探测信号辐射到空间中，被目标体反射，由接收阵元接收回波数据，建立方程(1)并给出初始化条件；

步骤(2)：求解如下子问题

其中，i＝x和y。上角标(k)表示第k次迭代结果。该问题可以用收缩算子有效地求解，这是一个闭型解。

步骤(3)：进一步，采用单步Bi-CGSTAB算法求解如下子问题

步骤(4)：采用如下方程迭代更新变量

步骤(5)：将步骤(2)——步骤(4)构成迭代，直至达到收敛要求，即满足达到最大迭代次数或容差小于设定值。

所述步骤(1)的步骤为：

在SIMO阵列中，线性阵列沿X轴方向排布，散射场S(x_r；k)是接收阵元在(x_r，0)处接收的回波，k表示波数，并可表示为

其中，f是采样信号的频率分量，c是信号在空间中的传播速度。初始化条件为

u⁽⁰⁾＝A^Hf；

其中，上角标(0)表示初始化阶段。

所述步骤(2)的步骤为：

利用收缩算子求解最小化子问题(7)，得到

其中，i＝x和y。

和/>

分别表示X和Y维度上的方向导数。

所述步骤(3)的步骤为：

最小化问题(8)等效于求解如下最佳条件子问题

Ku^(k+1)＝rhs^(k) (11)

这一步在算法中占有主要计算量，因此有必要采用一种高效方法求解(11)。利用矩阵求逆或伪逆求解(11)是最直接的方法。然而，矩阵求逆计算量非常大，有时算子A作为隐函数是无法显式表达，或者表达成矩阵后占用内存过大，特别是在三维情况下。更合理的解决方案是使用迭代求解方法。本发明提出了单步Bi-CGSTAB算法来求解(11)，因此这种方法不需要迭代收敛性检查，在任何情况下，一步迭代就足够了。原因是，一方面是由于Bi-CGSTAB具有超快收敛速度；另一方面，分裂Bregman方法收敛于最终解，不需要严格求解步骤(3)中的子问题。因此，可以有效地解决TV正则化问题。

SIMO稀疏阵列成像系统原理图如图1所示。发射天线阵元(1)位于原点，向自由空间辐射探测信号与待成像目标(3)相互作用，反射信号由接收天线阵元(2)接收并记录，本实施例假设线性阵列与X轴重合。阵列中共包含82个接收天线，阵元间距为6mm，系统的工作频段为26–40GHz。系统的工作方式如下：82个接收天线中，随机选择阵元总数约1/3的接收阵元工作，即27个接收阵元工作，其他阵元不工作。这样就获得了在该频段的SIMO稀疏阵列。

实施例一：电磁仿真成像

该实施例通过电磁仿真验证本发明提出的成像算法的成像性能。电磁仿真采用矩量法，源采用理想电偶极子，在工作的接收阵元位置上计算目标的散射场。成像目标是直径为5cm的金属球。分别采用传统的RA和改进的单步Bi-CGSTAB分割Bregman方法(即本发明提出的方法)对图像进行重建，结果如图2和图3所示。图2具有明显的旁瓣和栅瓣和较高的背景噪声，这是由于采用了稀疏阵列引起的。而图3的图像重建质量明显更高，这是由于该方法具有更高的精度，而且TV正则化具有较强的压制基底噪声能力。

实施例二：成像实验结果

实验采用了与实施例一相同的稀疏阵列几何结构。将直径为5cm的金属球作为目标(6)位于距离阵列0.5m左右的位置，如图4所示。图4中，发射阵元(4)位于阵列中心，即坐标原点。接收阵列(5)沿箭头方向排布。分别采用传统的RA和单步Bi-CGSTAB改进分割Bregman方法(即本发明提出的方法)对图像进行重建，结果如图5和图6所示。由实验结果可见，本发明提出的成像方法在抑制栅瓣和背景噪声方面明显优于RA，与仿真结果吻合较好。在本实验中，容差与迭代次数的关系曲线如图7所示，只需要16次迭代就能够达到10^-3的相对容差。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (2)

1.一种用于稀疏阵列的全变差快速成像算法，其特征在于：首先，对分裂Bregman方法的求解顺序进行改进；其次，引入Bi-CGSTAB算法求解最优条件子问题，使算法快速收敛，并且在任何条件下，单步Bi-CGSTAB改进分割Bregman成像方法足以保证成像效果；最后，分裂Bregman方法在几十次迭代之内以10^-3的容差收敛到最终解；

在单输入多输出SIMO线性阵列中，接收回波表示为

S(x_r；k)＝∫∫σ(x，y)×exp[-jk(R_t+R_r)]dxdy (1)

其中，σ(x，y)是目标的反射率函数，k是波数；

设SIMO线性阵列沿X轴方向排布，散射场S(x_r；k)是接收阵元在(x_r，0)处接收的回波，S(x_r；k)视为关于x_r和k的二维矩阵；R_t和R_r分别是发射阵元和接收阵元与散射点(x，y)的单程距离；

将回波数据S(x_r；k)变形为单列数据

在成像域对反射率函数σ(x，y)进行网格剖分，并写成单列数据的形式

对于公式(1)，将积分算子表示为

公式(1)表示为矩阵乘法形式

Au＝f (4)

为提高图像质量，将TV正则化引入公式(4)的求解，即

根据Bregman方法，求解公式(5)等效于求解如下问题：

将图像重建问题转化为求解公式(6)的最小化问题；

所述单步Bi-CGSTAB改进分割Bregman成像方法，具体包括如下步骤：

S1：SIMO稀疏阵列数据采样；

在SIMO稀疏阵列中，发射阵元位于线性阵列的中心位置，接收阵元以任意方式排布于线性维度上；即在笛卡尔坐标系中，天线阵列与X轴重合；发射阵元将探测信号辐射到空间中，被目标体反射，由接收阵元接收回波数据，建立方程(1)并给出初始化条件；

S2：求解如下子问题

其中，i＝x和y；上角标(k)表示第k次迭代结果；

该子问题用收缩算子有效地求解，这是一个闭型解；

S3：进一步，采用单步Bi-CGSTAB算法求解如下子问题

S4：采用如下方程迭代更新变量

S5：将S2-S4构成迭代，直至达到收敛要求，即满足达到最大迭代次数或容差小于设定值。

2.根据权利要求1所述的一种用于稀疏阵列的全变差快速成像算法，其特征在于：所述S1的步骤为：

在SIMO阵列中，线性阵列沿X轴方向排布，散射场S(x_r；k)是接收阵元在(x_r,0)处接收的回波，k表示波数，并表示为

其中，f’是采样信号的频率分量，c是信号在空间中的传播速度；

初始化条件为u⁽⁰⁾＝A^Hf'；

其中，上角标(0)表示初始化阶段；

所述S2的步骤为：

利用收缩算子求解最小化子问题(7)，得到

其中，i＝x和y；

和

分别表示X和Y维度上的方向导数；

所述S3的步骤为：

最小化公式(8)等效于求解如下最佳条件子问题

Ku^(k+1)＝rhs^(k) (11)

利用单步Bi-CGSTAB算法求解(11)。

Images