太赫兹科学与电子信息学报  2020, Vol. 18 Issue (6): 1035-1039     DOI: 10.11805/TKYDA2019263
基于截断sinc插值的毫米波全息成像算法优化    [PDF全文]
邵文浩, 朱莉, 刘婕, 邹丽蓉     
南京理工大学 电子工程与光电技术学院,江苏 南京 210094
摘要: 针对毫米波系统在二维横断面成像时,空间频率间隔上非均匀性影响成像精确度问题,提出基于截断正弦基数(sinc)函数的插值优化算法。算法取sinc函数为卷积核,对数据进行插值,相较细胞元插值法,得到的散射插值幅值更接近实际值,成像精确度更高。在此基础上,由于sinc插值需要无穷项求和,所需计算的数据量巨大,因此对插值求和项进行截断处理,对不同截断点数所成横断面像作比较,进一步优化成像质量。仿真结果验证了算法的可行性和精确性。
关键词: 毫米波    全息成像    二维横断面像    sinc插值    
Optimization of millimeter wave holographic imaging algorithm based on truncated sinc interpolation
SHAO Wenhao, ZHU Li, LIU Jie, ZOU Lirong     
School of Electronic and Optical Engineering, Nanjing University of Science & Technology, Nanjing Jiangsu 210094, China
Abstract: Aiming at solving the problem that the non-uniformity of spatial frequency affects the accuracy in the process of two-dimensional cross-section imaging, an interpolation optimization algorithm based on truncated sine cardinal(sinc) function is proposed. The algorithm takes the sinc function as the convolution kernel to interpolate the data. Compared with the cell element interpolation method, the obtained scattering interpolation amplitude is closer to the actual value, and the imaging accuracy is higher. On this basis, sinc interpolation requires infinite summation and the amount of data to be calculated is huge, therefore, the interpolation summation should be truncated. The cross-sectional images formed by different truncation points are compared and processed by windowing to further optimize image quality. The simulation results verify the feasibility and accuracy of the optimization algorithm.
Keywords: millimeter wave    holographic imaging    two-dimensional cross-sectional image    sinc interpolation    

毫米波被广泛应用在成像系统中。一方面,毫米波波长较短,和低频段的微波相比,分辨力更高;另一方面,毫米波兼具红外和可见光的成像特性,有一定强度的穿透力,能够做到对隐匿物体的探测,且能够进行全时段的工作[1]。因此毫米波全息成像技术有着很广阔的研究前景,特别是在近场成像方面,例如生物医学上对患者的检查、安检上对隐匿违禁物品的检测、军事上在沙尘暴或雾霾等极端天气情况下对近场目标的探测等[23]。针对毫米波近场成像系统的使用环境,需要寻求一种精确度高、速度快的成像算法。基于波数域的对回波信号进行处理的成像算法,相较于基于空间域的算法,能够较好地抑制旁瓣拖影和像点发散的情况,提升成像质量,因此成为目前研究的重点[45]。本文首先介绍了基于波数域的毫米波二维横断面成像原理,然后针对算法实现中的关键步骤,即匹配滤波后非均匀数据的插值问题,提出了基于sinc函数的插值算法。使用该算法得到的散射插值幅度最接近实际值,因此理论上具有极高的成像精确度。由于sinc插值需要无穷项求和,计算数据量巨大,因此在具体实现中,需要对插值求和项进行截断处理,分别采用4点、8点、16点截断对数据进行插值处理,进一步优化成像质量。最后对优化算法进行仿真分析,验证了该算法的可行性和精确性。

1 近场毫米波二维横断面成像原理

被测目标的散射强弱等效到某一平面上形成二维平面像,横断面像反映的是被测目标沿着纵向和横向水平方向所形成平面上的散射分布。在坐标系中,横断面像是坐标平面xoz内的像,因此成像系统实现横断面像需要具备沿xz两个方向的分辨能力[67]

成像模型见图 1,多个天线阵元构成一维扫描阵列,阵列长度为${L_x}$,任意时刻只有一个天线单元处于工作状态,该天线单元工作时,发射步进频率信号,起始频率为${f_0}$,频率步进为$\Delta f$,共有M个频率点,最大频率为${f_{\max }}$。依靠直线阵列实现方位像成像,依靠步进频率信号实现距离像成像,二者结合可得二维横断面像[8]

Fig.1 Model of cross-sectional imaging system 图 1 横断面成像系统模型

设目标在以$(x = 0, z = 0)$点为中心的平面区域内具有电场分布${\mathit{\boldsymbol{E}}_1}(x, z;x = 0, z = 0)$,区域平面重合于坐标面xoz,在以$(x = 0, z = {R_0})$点为中心的区域具有电场分布${\mathit{\boldsymbol{E}}_2}(x, z;x = 0, z = {R_0})$,区域平面重合于坐标面xoz。目标区域电场分布${\mathit{\boldsymbol{E}}_1}(x, z;x = 0, z = 0)$等价于待求的目标像函数[9]

采用步进频信号一维横向采样的测量方式,得到的散射数据可以认为是不同波数$k$时,电场分布${\mathit{\boldsymbol{E}}_2}(x, z;x = 0, z = {R_0})$$z = {R_0}$的值[10],写成$\mathit{\boldsymbol{E}}(x, k;x = 0, z = {R_0})$

首先对$\mathit{\boldsymbol{E}}(x, k;x = 0, z = {R_0})$沿变量𝑥作傅里叶逆变换,将回波数据变换到波数域:

$ F({k}_{x}, k;x=0, z={R}_{0})={\displaystyle \int {E}_{\rm{2}}(x, k;x=0, z={R}_{0}){{\rm{e}}^{{\rm{j}}{\mathit{k}_\mathit{x}}\mathit{x}}}{\rm{d}}{\mathit{k}_\mathit{x}}} $ (1)

构建波数域相位补偿因子,即线性色散空间滤波器,其数学表达式为:

$ H({k_x}, {k_z}) = rect\left( {\frac{{\left| {{k_x}} \right|}}{{\sqrt {k_x^2 + k_z^2} }}} \right){{\rm{e}}^{ - {\rm{j}}{k_z}{R_0}}} $ (2)

利用相位补偿因子对测量平面波谱进行匹配滤波,得到目标平面波谱:

$ F({k_x}, k;x = 0, z = 0) = F({k_x}, k;x = 0, z = {R_0})H({k_x}, {k_z}) $ (3)

上式得到的目标平面波谱$F({k_x}, k;x = 0, z = 0)$以非正交的波数$k$和空间频率${k_x}$为自变量,需将波数$k$代换为${k_z}$。波数$k$与空间频率${k_x}$, ${k_z}$满足如下关系:

$ k_x^2 + k_z^2 = {k^2} $ (4)

经代换,得到以空间频率${k_x}$${k_z}$为自变量的目标平面波谱$F({k_x}, {k_z};x = 0, z = 0)$。代换后,在${k_x} - k$坐标下均匀分布的回波数据变换至${k_x} - {k_z}$坐标下是非均匀的,因而代换过程必须借助相应的插值算法实现。利用sinc插值,从非均匀分布的${k_x} - k$域回波数据得到均匀分布的${k_x} - {k_z}$域数据。

对插值完的$F({k_x}, {k_z};x = 0, z = 0)$进行二维傅里叶变换,得到目标区域电场分布${{\boldsymbol{E}}_1}(x, z;x = 0, z = 0)$,即:

$ {{\boldsymbol{E}}_{\rm{1}}}(x, z;x = 0, z = 0) = \frac{1}{{4{{\rm{ \mathsf{ π} }}^2}}}\iint {F({k_x}, {k_z};x = 0, z = 0){{\rm{e}}^{ - {\rm{j}}({k_x}x + {k_z}z)}}{\rm{d}}{k_x}{\rm{d}}{k_z}} $ (5)

依据上述原理分析,波数域相位补偿横断面成像算法的具体步骤如下[11]:1)对回波采样数据作方位维傅里叶逆变换,将其变换至波数域;2)根据天线阵列方位构建波数域相位补偿因子;3)利用相位补偿因子对回波信号作匹配滤波;4)对滤波后的数据进行sinc插值,实现数据在${k_x} - {k_z}$域的均匀分布;5)对${k_x} - {k_z}$域数据作二维快速傅里叶变换,得到横断面像。

基于波数域相位补偿的横断面全息成像算法流程如图 2所示。

Fig.2 Flow of cross-sectional imaging 图 2 二维横断面成像流程图
2 sinc插值算法原理

波谱矩阵$F({k_x}, k;x = 0, z = 0)$数据的分布情况见图 3。中间呈倒梯形分布的“●”数据代表成像的有效数据,倒梯形之外的数据为无效数据(隐失波)。有效数据代表的平面波可以传播能量,参与目标图像的形成,在${k_x} - {k_z}$坐标下均匀分布。根据色散关系式$k_x^2 + k_z^2 = {k^2}$,将${k_x} - k$坐标下数据变换到$ {k_x} - {k_z}$坐标下,数据的分布情况见图 4。经色散关系变换后,在${k_x} - k$坐标下均匀分布的数据在${k_x} - {k_z} $坐标下是非均匀的,呈扇形分布,${k_z}$的两端,即最大值与最小值端为双曲线,且随着${k_z}$的增大,数据的分布变得稠密[12]。sinc插值的目的,是使在${k_x} - {k_z}$坐标下呈扇形分布(分布在弧线上)的数据均匀分布,分布情况见图 5。对于任一${k_x}$,在不同波数$k$下,代换后的波数分量${k_z}$分布在不同弧线与${k_{x1}}$直线交点上,用“▲”表示,分布情况见图 6

Fig.3 Spatial frequency distribution of data in kx−k 图 3 kx−k坐标下数据空间频率分布
Fig.4 Spatial frequency distribution of nonuniform data in kx−kz 图 4 kx−kz坐标下非均匀数据空间频率分布
Fig.5 Spatial frequency distribution of uniform data in kx−kz 图 5 kx−kz坐标下均匀数据空间频率分布
Fig.6 Spatial frequency distribution of data under a single kx 图 6 单个kx下数据空间频率分布

插值的目的是要根据非均匀分布的波谱数据$F({k_x}, {k_z})$插值得出均匀分布的$F'({k_x}, {k_z})$。sinc函数插值算法是一种内插法,通常情况下函数本身表达式是未知的,利用已知位置的函数值通过内插的方法,得到任意位置处函数值[13]

根据Nyquist采样定理,对于带限信号$f(t)$,若采样频率${f_{\rm{s}}}$大于其带宽的两倍,则可用一系列离散时刻的采样值恢复原始信号,即:

$ f(t) = \sum\limits_{n = - \infty }^\infty {f\left( {\frac{n}{{{f_{\rm{s}}}}}} \right)\frac{{\sin \left[ {{\rm{ \mathsf{ π} }}{f_{\rm{s}}}\left( {t - \frac{n}{{{f_{\rm{s}}}}}} \right)} \right]}}{{{\rm{ \mathsf{ π} }}{f_{\rm{s}}}\left( {t - \frac{n}{{{f_{\rm{s}}}}}} \right)}}} = \sum\limits_{n = - \infty }^\infty {f\left( {\frac{n}{{{f_{\rm{s}}}}}} \right)} \operatorname{sinc} \left[ {{\rm{ \mathsf{ π} }}{f_{\rm{s}}}\left( {t - \frac{n}{{{f_{\rm{s}}}}}} \right)} \right] $ (6)

插值需要无穷项求和,具体实现时,在有限的计算资源下,需要对插值求和项进行截断处理。

对于P点的sinc插值[14],具体插值公式为:

$ F'({k_x}, {k_z}) = F'({k_x}, {k_j} + \Delta {k_{{k_z}}}) = \sum\limits_{n = - \frac{P}{2}}^{\frac{P}{2} - 1} {F({k_x}, {k_j} + ndk){l_n}} $ (7)
$ {k_j} + \Delta {k_{{k_z}}} = \sqrt {k_x^2 + k_z^2} $ (8)
$ {l_n} = \operatorname{sinc} \left[ {{\rm{ \mathsf{ π} }}( - n + \Delta {k_{{k_z}}})} \right], {\rm{ }}n = - \frac{P}{2}, \cdots , \frac{P}{2} - 1 $ (9)

式中:P为截断点数;$ F^{\prime}\left(k_{x}, k_{z}\right)$为待求的在${k_x} - {k_z}$坐标下均匀分布的波谱数据,$F'({k_x}, {k_j} + \Delta {k_{{k_z}}})$为均匀波谱数据$F({k_x}, {k_z})$根据色散公式代换所对应的非均匀波谱数据,$\Delta {k_{{k_z}}}$${k_x} - k$坐标下原始均匀数据$F({k_x}, k)$与非均匀数据$F'({k_x}, {k_j})$的差值;${k_j}$为宽带信号对应的第j个波数;$F({k_x}, {k_j} + ndk)$为均匀波谱数据$F({k_x}, {k_z})$在截断区间内的P个等间隔的取值,dk为空间频率间隔。

匹配滤波后,将得到的均匀分布的$F({k_x}, k)$利用色散关系式代换,得到不等间隔的空间频率分量${k_x}$${k_z}$。利用空间频率分量${k_z}$上的最大值和最小值以及数据个数,得到等间隔的${k_x}$${k_z}$,插值的目的就是要求得均匀的$F'({k_x}, {k_z})$

${k_x}$${k_z}$利用色散关系式$k_x^2 + k_z^2 = {k^2}$代换可得到在${k_x} - k$坐标下非均匀分布的$F'({k_x}, k)$$F'({k_x}, k)$等价于待求的$F'({k_x}, {k_z})$。为了得到$F'({k_x}, {k_z})$的数据,对均匀的$F({k_x}, k)$用sinc函数作为卷积核进行卷积操作,得到任意所求对应点处的$F'({k_x}, k)$,等价于求得$F'({k_x}, k)$

插值的流程图如图 7所示。由于${k_j}({k_x}, {k_{zi}})$$\Delta {k_{{k_{zi}}}}({k_x}, {k_{zi}})$${l_n}$都随着${k_x}$, ${k_{zi}}$的变化而变化,因此对于每一个待求的$F'({k_x}, {k_{zi}})$,都必须重新计算以上3个参数。

Fig.7 Flow of sinc interpolation 图 7 sinc插值流程图
3 二维横断面成像仿真分析

使用Matlab进行二维横断面像的仿真实验分析[1516]。实验条件为:起始工作频率${f_0} = 30{\rm{ GHz}}$,步进频率$\Delta f = 10{\rm{ MHz}}$,频率步进个数$M = 1\;024$,扫描点步进长度$\Delta d = 5{\rm{ mm}}$,天线阵列扫描长度${L_x} = 2{\rm{ m}}$,扫描点与目标中心距离${R_0} = 2{\rm{ m}}$。分别在xoz坐标平面内以(0, 0)为圆心,0.2 m为半径的圆上放置12个散射强度为1的目标小球,构成圆环图案。

对数据仅进行匹配滤波,不进行插值,所成的二维横断面像如图 8所示。由图 8可知,基于波数域的二维横断面成像算法能够显示出目标像素点,但是不够清晰,这是由于各个目标点处的像素点波瓣展宽引起的,各个点的像之间又发生混叠,造成成像质量的下降,无法有效识别各个目标点[7]。因此引入sinc插值算法。对匹配滤波后的数据进行4点截断sinc插值,经仿真后所成的二维横断面像如图 9所示。从图 9可知,经过sinc插值算法处理后,与图 8相比,各个像素点发散的情况得到一定抑制,成像质量显著提升,验证了sinc插值算法的可行性。

Fig.8 Two-dimensional cross-sectional image using wavenumber domain phase compensation 图 8 波数域相位补偿二维横断面像
Fig.9 Two-dimensional cross-sectional image with 4-point truncated sinc interpolation 图 9 4点截断sinc插值二维横断面像

在此基础上,为了进一步提高成像质量,对匹配滤波数据分别进行8点、16点截断sinc插值处理,得到的二维横断面像见图 10图 11。经更高截断点数的sinc插值处理后,图 10图 11的各个像点进一步收敛,能量更加集中,成像质量有所提升,随着计算量的增大,运算时间也将增加。

Fig.10 Two-dimensional cross-sectional image with 8-point truncated sinc interpolation 图 10 8点截断sinc插值二维横断面像
Fig.11 Two-dimensional cross-sectional image with 16-point truncated sinc interpolation 图 11 16点截断sinc插值二维横断面像

表 1以(0.1, 0.173 2)处的散射点为例,比较使用插值算法前后及不同截断点数插值的仿真结果的精确度。表 1中的数据显示了使用插值算法前后像点在x, z方向的误差,sinc插值算法能有效降低成像的畸变,降低旁瓣发散的影响,提高成像质量。此外,采用不同截断点数的sinc插值成像的像点位置无明显差异(差值小于Matlab仿真结果图片分辨力而无法测得),可知更高截断点数的插值对成像畸变的改善甚小。此外,4点、8点、16点截断插值算法的计算量是成倍提升的,使得插值更加耗时,但更高截断点数的精确度提升却很小。在有限的计算资源下,8点截断sinc插值比较适合数据处理,且能满足成像精确度要求。

表 1 插值算法前后(0.1, 0.173 2)处散射点的成像精确度比较 Table 1 Comparison of imaging accuracy of scattering points (0.1, 0.173 2) before and after interpolation algorithm
4 结论

本文首先介绍了近场毫米波二维横断面成像原理,对成像基础算法加以说明,并进行仿真实验。针对基础算法仿真结果中存在的图像畸变、像点发散、混叠等问题,提出了截断sinc插值算法对匹配滤波数据进行处理,该算法的核心思想是在插值过程中采用截断sinc函数作为卷积核,以降低数据空间频率上的非均匀性对成像结果的影响。仿真结果验证了二维横断面成像基础算法的可行性,以及截断sinc插值算法在抑制图像畸变、像点发散,提升成像质量方面的有效性。后续的研究将在毫米波三维成像算法以及数据处理方向进行。

参考文献
[1]
MOHAMMADIAN Nafiseh, FURXHI Orges, SHORT Robert, et al. SAR millimeter wave imaging systems[C]// Passive and Active Millimeter-Wave Imaging XXⅡ. Baltimore, Maryland, USA: SPIE, 2019: 109940A-1-13.
[2]
LI X, LI S Y, ZHAO G Q, et al. Multi-polarized millimeter-wave imaging for concealed weapon detection[C]// IEEE International Conference on Microwave and Millimeter Wave Technology. Beijing: IEEE, 2016: 892-894.
[3]
LEE D J, KIM S, SONG D R, et al. Photonic-assisted on-wafer millimeter-wave imaging system[C]// 41st International Conference on Infrared, Millimeter, and Terahertz Waves. Copenhagen, Denmark: [s.n.], 2016: 1-2.
[4]
保铮, 邢孟道, 王彤. 雷达成像技术[M]. 北京: 电子工业出版社, 2005. (BAO Zheng, XING Mengdao, WANG Tong. Radar imaging technology[M]. Beijing: Electronic Industry Press, 2005.)
[5]
APPLEBY Roger, ROBERTSON Duncan A, WIKNER David. Millimeter wave imaging: a historical review[C]// SPIE Defense + Security. Anaheim, California, USA: SPIE, 2017: 1018902-1-13.
[6]
QIAO L, WANG Y, ZHAO Z, et al. Exact reconstruction for near-field three-dimensional planar millimeter-wave holographic imaging[J]. Journal of Infrared, Millimeter, and Terahertz Waves, 2015, 36(12): 1221-1236.
[7]
王子野, 乔灵博, 王迎新. 高分辨力亚毫米波全息成像系统[J]. 太赫兹科学与电子信息学报, 2016, 14(6): 833-837. (WANG Ziye, QIAO Lingbo, WANG Yingxin. Wide-band three-dimensional submillimeter-wave holographic imaging system[J]. Journal of Terahertz Science and Electronic Information Technology, 2016, 14(6): 833-837.)
[8]
张铁男.毫米波三维全息成像图像重构算法研究[D].南京: 南京理工大学, 2016. (ZHANG Tienan. Research on image reconstruction algorithm of millimeter wave 3D holographic image[D]. Nanjing, China: Nanjing University of Science and Technology, 2016.)
[9]
史耀亮.毫米波雷达近程平面扫描散射成像算法研究[D].南京: 南京理工大学, 2015. (SHI Yaoliang. Research on short-range planar scanning scattering imaging algorithm for millimeter wave radar[D]. Nanjing, China: Nanjing University of Science and Technology, 2015.)
[10]
张麟兮, 李南京, 胡楚锋, 等. 雷达目标散射特性测试与成像诊断[M]. 北京: 中国宇航出版社, 2009: 145-220. (ZHANG Lingxi, LI Nanjing, HU Chufeng, et al. Radar target scattering characteristics test and imaging diagnosis[M]. Beijing: China Aerospace Publishing House, 2009: 145-220.)
[11]
KAN Y Z, ZHU Y F, FU Q. A direct integral imaging method for near-field 3-D imaging[C]// IEEE Progress in Electromagnetic Research Symposium. Shanghai, China: IEEE, 2016.
[12]
徐枫, 朱莉, 李小辉, 等. 基于Stolt算法的毫米波全息二维成像优化[J]. 电光与控制, 2019, 26(2): 62-65. (XU Feng, ZHU Li, LI Xiaohui, et al. Two-dimensional imaging optimization of millimeter wave holography based on Stolt algorithm[J]. Electro-Optic and Control, 2019, 26(2): 62-65.)
[13]
张贵平, 贾鑫, 尹灿斌, 等. SAR成像处理中插值算法的比较研究[J]. 国外电子测量技术, 2008, 27(8): 9-11. (ZHANG Guiping, JIA Xin, YIN Canbin, et al. Comparative study of interpolation algorithms in SAR imaging processing[J]. Foreign Electronic Measurement Technology, 2008, 27(8): 9-11.)
[14]
MAO N H, ZHANG X M. Novel technique for planar near-field scattering measurements[C]// IEEE Asia-pacific Microwave Conference. Hong Kong, China: IEEE, 1997: 897-900.
[15]
陈丽.合成孔径声呐成像算法研究及MATLAB仿真[D].杭州: 杭州电子科技大学, 2016. (CHEN Li. Synthetic aperture sonar imaging algorithm research and MATLAB simulation[D]. Hangzhou, China: Hangzhou University of Electronic Science and Technology, 2016.)
[16]
朱莉, 李兴国, 王本庆. 近程目标毫米波全息成像算法及仿真[J]. 电光与控制, 2011, 18(6): 37-40. (ZHU Li, LI Xingguo, WANG Benqing. MMW holographic imaging algorithm and simulation for short-distance target[J]. Electro-Optic and Control, 2011, 18(6): 37-40.)