一种多频近场毫米波稀疏重建图像方法

本发明涉及图像重建方法，具体涉及一种多频近场毫米波稀疏重建图像方法。

背景技术：

对于工作在单频的近场毫米波成像系统而言，其成像质量往往会受到噪声、重建伪影、系统硬件限制等多方面的因素影响，为保证相对稳定的图像重建质量，可选择工作在多频的近场毫米波成像系统，对各频率下的重建图像进行平均处理，生成质量更为稳定的图像。对于单频近场毫米波稀疏成像系统，由l1范数和tv算子所构成的混合稀疏函数已经可以表征图像的稀疏性，使算法能够从低欠采样率下所获取的稀疏观测数据中重建图像。而对于多频近场毫米波稀疏成像系统，由于被测对象对不同频率的毫米波后向散射程度不同，重建图像所体现的稀疏程度也会有差异，此时需要考虑结合更适合度量图像稀疏性的混合稀疏函数。

技术实现要素：

针对现有技术中的上述不足，本发明提供的一种多频近场毫米波稀疏重建图像方法解决了如何有效从多频稀疏观测数据中重建图像的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种多频近场毫米波稀疏重建图像方法，包括以下步骤：

s1、通过多频近场毫米波器件扫描待测区域，得到多频稀疏观测序列；

s2、根据多频稀疏观测序列，通过基于基尼加权范数结构张量全变分融合算子的混合成像算法，得到重建而成的待测区域的扫描图像。

进一步地，所述步骤s2包括以下步骤：

s21、在计算机系统中设定迭代标志i、图像重建序列、第一一级缓存序列、中间因子、第一二级缓存序列和第二二级缓存序列的初值；

s22、根据当前的第一一级缓存序列、中间因子、第一二级缓存序列和第二二级缓存序列，基于基尼加权范数结构张量全变分融合算子的混合成像算法的迭代方程组迭代图像重建序列，并将迭代标志i自加1；

s23、判断基尼加权范数结构张量全变分融合算子是否大于中断容限若是，则跳转至步骤s24，若否，跳转至步骤s25；

s24、判断迭代标志i是否大于迭代上限，若是，则跳转至步骤s25，若否，则跳转至步骤s22；

s25、保存当前的图像重建序列，得到重建而成的待测区域的扫描图像。

进一步地，所述步骤s21中迭代标志i的初值设定为1，图像重建序列、第一一级缓存序列、中间因子、第一二级缓存序列和第二二级缓存序列的初值按以下各式设定：

g⁽ⁱ⁾＝m^#s

y⁽ⁱ⁾＝g⁽ⁱ⁾

γ⁽ⁱ⁾＝1

其中，g⁽ⁱ⁾为第i次迭代的图像重建序列，m^#为矩阵m的逆矩阵，矩阵m为多频近场毫米波器件的传输矩阵，s为多频稀疏观测序列，y⁽ⁱ⁾为第i次迭代的第一一级缓存序列，γ⁽ⁱ⁾为第i次迭代的中间因子，为第i次迭代的第一二级缓存序列，为第i次迭代的第二二级缓存序列，jk为加权块雅克比算子，w为加权矩阵。

进一步地，所述步骤s22中基于基尼加权范数结构张量全变分融合算子的混合成像算法的迭代方程组包括以下各式：

其中，g⁽ⁱ⁾为第i次迭代的图像重建序列，∏c(·)为实空间c的投影，b⁽ⁱ⁾为第一迭代暂存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，τ为极小化lipschitz常数，*为卷积，w为加权矩阵，jk为加权块雅克比算子，为的投影，为l∞-sq范数单位球空间，为的投影，为基尼指数gin单位球空间，y⁽ⁱ⁺¹⁾为第i+1次迭代的第一一级缓存序列，γ⁽ⁱ⁾为第i次迭代的中间因子，为第i次迭代的第一二级缓存序列，为第i次迭代的第二二级缓存序列。

进一步地，所述第一迭代暂存序列b⁽ⁱ⁾的计算公式为：

其中，y⁽ⁱ⁾为第i次迭代的第一一级缓存序列，τ为极小化lipschitz常数，m^#为矩阵m的逆矩阵，矩阵m为多频近场毫米波器件的传输矩阵，s为多频稀疏观测序列。

进一步地，所述的投影的计算公式为：

其中，为的投影，为第i次迭代的第一二级缓存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，τ为极小化lipschitz常数，ρ1为第一moreau包络近似性能平衡参数，l1为第一对偶目标子函数梯度lipschitz常数，jk为加权块雅克比算子，∏c(·)为实空间c的投影，为第二迭代暂存序列。

进一步地，所述第二迭代暂存序列的计算公式为：

其中，τ为极小化lipschitz常数，b⁽ⁱ⁾为第一迭代暂存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，ρ1为第一moreau包络近似性能平衡参数，g⁽ⁱ⁾为第i次迭代的图像重建序列，jk为加权块雅克比算子，为第i次迭代的第一二级缓存序列。

进一步地，所述的投影的计算公式为：

其中，为的投影，为第i次迭代的第二二级缓存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，τ为极小化lipschitz常数，ρ2为第二moreau包络近似性能平衡参数，l2为第二对偶目标子函数梯度lipschitz常数，w为加权矩阵，πc(·)为实空间c的投影，为第三迭代暂存序列。

进一步地，所述第三迭代暂存序列的计算公式为：

其中，τ为极小化lipschitz常数，b⁽ⁱ⁾为第一迭代暂存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，ρ2为第二moreau包络近似性能平衡参数，g⁽ⁱ⁾为第j次迭代的图像重建序列，w为加权矩阵，为第i次迭代的第二二级缓存序列。

进一步地，所述步骤s23中基尼加权范数结构张量全变分融合算子的计算公式为：

其中，g⁽ⁱ⁺¹⁾为第i+1次迭代的图像重建序列，g⁽ⁱ⁾为第i次迭代的图像重建序列，||·||2为2范数。

综上，本发明的有益效果为：考虑到毫米波在不同频率下对被测对象的后向散射程度的不同，重建图像的稀疏程度也具有差异。针对这个问题，本发明提出了一种多频近场毫米波稀疏重建图像方法，设计的混合成像算法可以有效表征多频近场毫米波情况下的稀疏性，与结合l1范数+tv算子混合稀疏函数的多频成像算法相比，本发明提出的基尼加权范数+stv(结构张量)算子混合稀疏函数的多频成像算法有着更好的成像质量。

附图说明

图1为一种多频近场毫米波稀疏重建图像方法的流程图；

图2(a)为实测实验被测对象；

图2(b)为36ghz～44ghz频率下多频近场毫米波全采样重建图像；

图3(a)为欠采样率为14％时，基尼加权l1范数+stv算子所重建图像；

图3(b)为欠采样率为21％时，基尼加权l1范数+stv算子所重建图像；

图3(c)为欠采样率为28％时，基尼加权l1范数+stv算子所重建图像；

图3(d)为欠采样率为14％时，基尼加权l1范数+tv算子所重建图像；

图3(e)为欠采样率为21％时，基尼加权l1范数+tv算子所重建图像；

图3(f)为欠采样率为28％时，基尼加权l1范数+tv算子所重建图像。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种多频近场毫米波稀疏重建图像方法，包括以下步骤：

s1、通过多频近场毫米波器件扫描待测区域，得到多频稀疏观测序列；

s2、根据多频稀疏观测序列，通过基于基尼加权范数结构张量全变分融合算子的混合成像算法，得到重建而成的待测区域的扫描图像。

所述步骤s2包括以下步骤：

s21、在计算机系统中设定迭代标志i、图像重建序列、第一一级缓存序列、中间因子、第一二级缓存序列和第二二级缓存序列的初值；

所述步骤s21中迭代标志i的初值设定为1，图像重建序列、第一一级缓存序列、中间因子、第一二级缓存序列和第二二级缓存序列的初值按以下各式设定：

g⁽ⁱ⁾＝m^#s

y⁽ⁱ⁾＝g⁽ⁱ⁾

γ⁽ⁱ⁾＝1

其中，g⁽ⁱ⁾为第i次迭代的图像重建序列，m^#为矩阵m的逆矩阵，矩阵m为多频近场毫米波器件的传输矩阵，s为多频稀疏观测序列，y⁽ⁱ⁾为第i次迭代的第一一级缓存序列，γ⁽ⁱ⁾为第i次迭代的中间因子，为第i次迭代的第一二级缓存序列，为第i次迭代的第二二级缓存序列，jk为加权块雅克比算子，w为加权矩阵。

s22、根据当前的第一一级缓存序列、中间因子、第一二级缓存序列和第二二级缓存序列，基于基尼加权范数结构张量全变分融合算子的混合成像算法的迭代方程组迭代图像重建序列，并将迭代标志i自加1；

步骤s22中基于基尼加权范数结构张量全变分融合算子的混合成像算法的迭代方程组包括以下各式：

其中，g⁽ⁱ⁾为第i次迭代的图像重建序列，πc(·)为实空间c的投影，b⁽ⁱ⁾为第一迭代暂存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，τ为极小化lipschitz常数，*为卷积，w为加权矩阵，jk为加权块雅克比算子，为的投影，为l∞-sq范数单位球空间，为的投影，为基尼指数gin单位球空间，y⁽ⁱ⁺¹⁾为第i+1次迭代的第一一级缓存序列，γ⁽ⁱ⁾为第i次迭代的中间因子，为第i次迭代的第一二级缓存序列，为第i次迭代的第二二级缓存序列。

第一迭代暂存序列b⁽ⁱ⁾的计算公式为：

其中，y⁽ⁱ⁾为第i次迭代的第一一级缓存序列，τ为极小化lipschitz常数，m^#为矩阵m的逆矩阵，矩阵m为多频近场毫米波器件的传输矩阵，s为多频稀疏观测序列。

所述的投影的计算公式为：

其中，为的投影，为第i次迭代的第一二级缓存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，τ为极小化lipschitz常数，ρ1为第一moreau包络近似性能平衡参数，l1为第一对偶目标子函数梯度lipschitz常数，jk为加权块雅克比算子，πc(·)为实空间c的投影，为第二迭代暂存序列。

所述第二迭代暂存序列的计算公式为：

其中，τ为极小化lipschitz常数，b⁽ⁱ⁾为第一迭代暂存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，ρ1为第一moreau包络近似性能平衡参数，g⁽ⁱ⁾为第i次迭代的图像重建序列，jk为加权块雅克比算子，为第i次迭代的第一二级缓存序列。

所述的投影的计算公式为：

其中，为的投影，为第i次迭代的第二二级缓存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，τ为极小化lipschitz常数，ρ2为第二moreau包络近似性能平衡参数，l2为第二对偶目标子函数梯度lipschitz常数，w为加权矩阵，πc(·)为实空间c的投影，为第三迭代暂存序列。

所述第三迭代暂存序列的计算公式为：

其中，τ为极小化lipschitz常数，b⁽ⁱ⁾为第一迭代暂存序列，λ1为第一稀疏权衡参数，λ2为第二稀疏权衡参数，ρ2为第二moreau包络近似性能平衡参数，g⁽ⁱ⁾为第i次迭代的图像重建序列，w为加权矩阵，为第i次迭代的第二二级缓存序列。

s23、判断基尼加权范数结构张量全变分融合算子是否大于中断容限若是，则跳转至步骤s24，若否，跳转至步骤s25；

所述步骤s23中基尼加权范数结构张量全变分融合算子的计算公式为：

其中，g⁽ⁱ⁺¹⁾为第i+1次迭代的图像重建序列，g⁽ⁱ⁾为第i次迭代的图像重建序列，||·||2为2范数。

s24、判断迭代标志i是否大于迭代上限，若是，则跳转至步骤s25，若否，则跳转至步骤s22；

s25、保存当前的图像重建序列，得到重建而成的待测区域的扫描图像。

实验：

在多频近场毫米波稀疏成像实测实验中，成像系统通过工作在36ghz～44ghz频段下的天线探头对128mm×128mm采样平面下方约60mm处的被测对象以约2mm的网格步进进行扫描，扫频间隔为0.1ghz(nf＝61)。图2给出了实测被测对象以及运行在36ghz～44ghz工作频率的近场毫米波系统对被测对象进行数据采集后重建的全采样重建图像。其中，图2(a)的被测对象为金属剪刀，图2(b)为36ghz～44ghz频段下的多频近场毫米波全采样重建图像，多频全采样重建图像可用来作为图像质量评估标准ssim和psnr的参考图像。

在本实验中，基尼加权l1范数选择具有8阶消失矩的daubechies小波，stv算子选择标准差为0.5的3×3大小高斯卷积核。此外，算法中所给出的稀疏权衡参数λ1＝4×10^-4，λ2＝2×10^-4；moreau包络近似性能平衡参数ρ1＝1，ρ2＝1；优化极小化lipschitz常数τ＝8；对偶目标子函数梯度lipschitz常数：

中断容限

图3展示了在不同欠采样率(14％，21％，28％)下，结合不同混合稀疏函数的多频近场毫米波稀疏成像算法重建图像效果。其中，图3(a)到图3(c)为结合基尼加权l1范数+stv算子混合稀疏函数的多频成像算法所重建图像，图3(d)到图3(f)为结合l1范数+tv算子混合稀疏函数的多频成像算法所重建图像。通过比较图3(a)和图3(d)可知，在欠采样率为14％时，选用基尼加权l1范数+stv算子的多频成像算法可将剪刀轴心的螺丝图案重建出来，而选用l1范数+tv算子的多频成像算法却没有成功重建。而在欠采样率为21％和28％时，选用基尼加权l1范数+stv算子的多频成像算法对右下角剪刀手柄重建效果优于选用l1范数+tv算子的多频成像算法。

为了更进一步说明两种混合稀疏函数下成像效果的比较，表1给出了在中断容限下独立重复50次实验后，混合稀疏函数的多频近场毫米波稀疏成像算法，在不同欠采样率(14％，21％，28％)下算法收敛时各指标的平均效果与对比。通过表1的对比可知，在不同欠采样率时，结合所提出的基尼加权l1范数+stv算子混合稀疏函数的成像算法，有着比结合l1范数+tv算子混合稀疏函数的成像算法更好的图像重建能力，说明了所提出的混合稀疏函数在实测实验中的有效性。

表1