CN110764088B - 一种超分辨率驻点扫描实时成像算法 - Google Patents

Abstract

本发明提供的是一种超分辨率驻点扫描实时成像算法。步骤(1)：天线在空间平面上对目标等间距扫描采样，并在频域记录回波数据；步骤(2)：采用距离偏移算法计算反射率函数的三维图像，在计算初始图像过程中、对插值后的回波进行三维逆快速傅里叶变换求取改进相干因子；步骤(3)：利用步骤(2)得到的改进相干因子校正初始图像。本发明的方法在波数域重新定义了相干因子的表达形式，在计算反射率函数非相干功率时，只需要一次三维快速傅里叶变换即可完成，因此在保证成像质量的前提下，将计算量从原来的2倍距离偏移算法降低到1.5倍，更适于实时性要求高的成像系统。

Description

一种超分辨率驻点扫描实时成像算法

技术领域

本发明涉及的是一种合成孔径成像方法，具体地说是一种基于改进相干因子的超分辨率驻点扫描实时成像算法。

背景技术

在毫米和太赫兹波成像中，广泛应用了合成孔径(synthetic aperture,SA)技术。除了适用于传统的雷达系统外,它还适用于材料的无损检测、安全监视和人体安检。无损检测是利用探测波对介质材料的内部结构进行精确的重建。在人体安检成像系统中，由于毫米波能够穿透服装布料和一些包装材料，因此可以利用毫米波对人体体表以毫米级空间分辨率进行重建。

为了提高毫米波成像的图像质量，一种可行的方法是利用散射公式建立了更为精确的数学模型。例如采用地震成像中常用的基尔霍夫偏移成像，相对于SA具有更高的成像质量，旁瓣电平被明显压低。将基尔霍夫偏移算法结合距离偏移算法，有利于提高算法的速度，同时具有更高的精度。将Rayleigh-Sommerfeld衍射公式应用于成像算法中，能够补偿传播损耗，提高图像质量。在声源描述方面，利用磁偶极矩和表面等效定理，精确地给出孔径场。另一类方法精确地测量点扩散函数标定系统性能从而提高分辨力。因此，只要选择适当的校准体，理论上就可以定量地重建目标反射率。

超分辨率(super resolution,SR)成像是一种新兴的改进型成像技术，能够有效地提高成像质量。SR最初的概念是指将一组降采样、低分辨率的图像序列结合成一幅高分辨率图像的过程。此外，从单个低分辨率输入图像生成高分辨率图像也属于SR的范畴。但是在本发明中，术语SR特指SA算法框架中有助于提高图像质量的方法，具体地说，在点扩展函数(PSF)方面，具有更尖锐的半最大值全宽度(FWHM)，更低的最大旁瓣和背景噪声。最近提出的一种半参数的、依赖于数据的功率谱密度估计器，称为迭代自适应方法(IAA)，它显示出优于传统估计器(如最小方差无失真响应估计器)的性能。进一步，通过探索矩阵Toeplitz结构的特征，在协方差矩阵求逆计算方面IAA显著加快。此外，有人提出引入奇异值分解方法，提出了一种快速的二维IAA算法。近年来，IAA已被广泛应用于SR成像。

此外，另一种重要的SR方法将成像问题转化为正则化优化问题，结合不同的正则化项，能够有效提高图像质量。如采用l₁范数正则化有利于描述点状目标，采用全变差方法有利于锐化连续目标体的边缘，从而提高成像质量。

总结起来，上述方法能够在不同程度上提高图像质量，但这些方法大多存在计算量和内存需求过大的问题，因此很少能直接应用于实时性成像系统中。

针对以上提出的超分辨成像方法的计算量高，内存需求大的技术缺陷，申请人之前层提出了一种“基干波数域相干因子的高速超分辨率驻点扫描成像方法”(专利号：ZL2017 1 0538119.8)，本发明申请是在此基础上提出的一种改进方案。

发明内容

本发明的目的在于提供一种成像分辨率高，内存占用量少，更适用于高实时性成像系统的超分辨率驻点扫描实时成像算法。

本发明的目的是这样实现的：

步骤(1)：天线在空间平面上对目标等间距扫描采样，并在频域记录回波数据；

步骤(2)：采用距离偏移算法计算反射率函数的三维图像，在计算初始图像过程中、对插值后的回波进行三维逆快速傅里叶变换求取改进相干因子；

步骤(3)：利用步骤(2)得到的改进相干因子校正初始图像。

本发明的超分辨率驻点扫描实时成像算法还可以包括：

1.所述步骤(1)具体包括：驻点扫描成像系统采用收发天线同置方式，在采样平面上，沿两个方向等间隔均匀采样，采样间隔满足空间Nyquist采样定理，在每一个采样位置上记录频域回波的宽带响应数据；

所述步骤(2)具体包括：

步骤(2.1)：对回波数据在空间域维度上做二维快速傅里叶变换，得到波数域回波信号；

步骤(2.2)：忽略沿传播路径上的信号衰减，根据已知的场景中心与采样平面距离，对步骤(2.1)得到的三维波数域回波信号进行场景中心补偿，使得目标的真实位置与成像域对齐；

步骤(2.3)：对步骤(2.2)得到的三维回波信号，根据色散关系式，在波数域维度上将回波数据插值到均匀分布的距离向波数域内，插值后的数据具有在三个维度上均匀采样的性质，这一步也称为重采样；

步骤(2.4)：对波数域重采样的回波数据进行三维逆快速傅里叶变换，得到反射率函数的初始图像；

步骤(2.5)：根据反射率函数非相干功率在波数域的计算公式，对步骤(2.3)得到的回波数据求模平方，并求三维逆快速傅里叶变换；最后计算改进相干因子，定义为反射率函数相干功率与非相干功率的比；

所述步骤(3)具体包括：

利用步骤(2.5)计算得到的改进相干因子校正(2.4)得到的反射率函数初始图像。

2.步骤(1)中，设采样平面在直角坐标系内与平面XY重合，设天线的坐标为(x′，y′，0)，则接收的回波信号s(x′，y′；k)表示为

其中，R是从散射点到收发天线的单程距离；D(x，y，z)是散射点所在的目标域，σ(x，y，z)表示待重建的反射率函数；k是波数的幅值，c是电磁波在自由空间中的传播速度。

3.步骤(2.1)中，对回波信号s(x′，y′，0；k)的两个空间维度x′和y′做二维快速傅里叶变换，得到波数域回波信号

s(k_x′，k_y′，k)＝∫∫∫_D(x，y，z)σ(x，y，z)exp(-ik_x′x-ik_y′y-ik_zz)dxdydz

其中，k_z是与距离向z方向对应的波数，设z＞0。

4.步骤(2.2)中，将目标域中心位置与目标距离对齐，设目标与天线的实际距离为H_c，则

s_c(k_x′，k_y′，k)＝s(k_x′，k_y′，k)exp(ik_zH_c)。

5.步骤(2.3)中，反射率函数的重建公式为

将采样点重采样至均匀分布的三维波数域(k_x′，k_y′，k_z)；重采样后回波表示为s_c(k_x′，k_y′，k_z)。

6.步骤(2.4)中，对重采样数据s_c(k_x′，k_y′，k_z)进行三维逆快速傅里叶变换重建反射率初始图像，表示为

7.步骤(2.5)中，所述的反射率函数非相干功率的波数域包括如下三种：

i.

ii.

iii.

8.步骤(3)中，所述的反射率函数初始图像表示为：

σ_m(x，y，z)＝σ(x，y，z)ICF(x，y，z)。

针对已有技术中的超分辨成像方法存在的计算量高、内存需求大的技术缺陷，本发明在之前提出的一种SR驻点扫描成像方法基础上(“基干波数域相干因子的高速超分辨率驻点扫描成像方法”，专利号：ZL 2017 10538119.8)，提出了基于改进相干因子的高速超分辨率驻点扫描实时成像方法及系统。该方法在波数域重新定义了相干因子的表达形式，在计算反射率函数非相干功率时，只需要一次三维快速傅里叶变换即可完成，因此在保证算法成像质量的前提下，将计算量从原来的2倍距离偏移算法降低到1.5倍，更适于实时性要求高的成像系统。

与现有技术相比，本发明的有益效果是：

1.将改进相干因子与RMA方法结合，显著提高了成像分辨率，并且对旁瓣和基底噪声具有明显的压制作用，获得更高的成像质量；

2.相对于本发明的前身方法——基于相干因子的RMA，即CF-RMA，本方法ICF-RMA具有更低的计算量和更少的内存占用量，更适用于高实时性成像系统。

附图说明

图1是驻点扫描成像示意图，其中：1为用于辐射和接收电磁信号的单站天线；

图2是基于改进相干因子的超分辨率驻点扫描实时成像算法流程图；

图3是采用RMA，CF-RMA和ICF-RMA三种方法得到的点扩展函数；

图4是仿真中采用的金属扇形目标；

图5是采用ICF-RMA的仿真成像结果；

图6是采用RMA的仿真成像结果；

图7是采用CF-RMA的仿真成像结果；

图8是成像实验中使用的金属目标；

图9是采用ICF-RMA的实验成像结果；

图10是采用RMA的实验成像结果；

图11是采用CF-RMA的实验成像结果。

具体实施方式

下面举例对本发明做更详细的描述。

结合图2，本发明的基于改进相干因子的超分辨率驻点扫描实时成像算法，包括如下步骤：

步骤(1)：天线在空间平面上对目标等间距扫描采样，并在频域记录回波数据；

步骤(2)：采用距离偏移算法(range migration algorithm，RMA)计算反射率函数的三维图像；所述步骤(2)在计算初始图像过程中，对插值后的回波进行三维逆快速傅里叶变换(IFFT)求取改进相干因子；

步骤(3)：利用步骤(2)得到的改进相干因子校正初始图像。

所述步骤(1)的步骤为：驻点扫描成像系统采用收发天线同置方式，在采样平面上，沿两个方向等间隔均匀采样，采样间隔需满足空间Nyquist采样定理，在每一个采样位置上记录频域回波的宽带响应数据。

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

步骤(2.1)：对回波数据在空间域维度上做二维快速傅里叶变换(FFT)，得到波数域回波信号；

步骤(2.2)：场景中心补偿：忽略沿传播路径上的信号衰减，根据已知的场景中心与采样平面距离，对步骤(2.1)得到的三维波数域回波信号进行场景中心补偿，使得目标的真实位置与成像域对齐；

步骤(2.3)：Stolt插值：对步骤(2.2)得到的三维回波信号，根据色散关系式，在波数域维度上将回波数据插值到均匀分布的距离向波数域内，插值后的数据具有在三个维度上均匀采样的性质，因此这一步也称为重采样；

步骤(2.4)：对波数域重采样的回波数据进行三维IFFT，得到反射率函数的初始图像；

步骤(2.5)：计算改进的相干因子：根据本发明提出的反射率函数非相干功率在波数域的计算公式，对步骤(2.3)得到的回波数据求模平方，并求三维IFFT；最后计算改进相干因子，定义为反射率函数相干功率与非相干功率的比。

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

利用步骤(2.5)计算得到的改进相干因子校正(2.4)得到的反射率函数初始图像。

驻点扫描成像是一种单站雷达成像方法，广泛应用于无损检测、医学成像等。

所述步骤(1)的步骤为：采样天线采用收发同置方式，假设采样平面在直角坐标系内与平面XY重合，天线沿两个方向等间隔均匀采样，采样间隔需满足空间Nyquist采样定理。假设天线的坐标为(x′，y′，0)，则接收的回波信号s(x′，y′；k)表示为

其中，R是从散射点到收发天线的单程距离；D(x，y，z)是散射点所在的目标域，σ(x，y，z)表示待重建的反射率函数；k是波数的幅值，c是电磁波在自由空间中的传播速度。

所述步骤(2)首先采用RMA计算反射率函数的三维图像，并对插值后的回波进行三维逆快速傅里叶变换(IFFT)求取改进相干因子。

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

对回波信号s(x′，y′，0；k)的两个空间维度x′和y′做FFT，得到波数域回波信号

s(k_x′，k_y′，k)＝∫∫∫_D(x，y，2)σ(x，y，z)exp(-ik_x′x-ik_y′y-ik_zz)dxdydz (3)

其中，k_z是与距离向z方向对应的波数，推导中假设z＞0。

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

场景中心补偿：考虑到FFT计算中具有周期性，必须将目标域中心位置与目标距离对齐，防止重建结果在距离向出现反褶。假设目标与天线的实际距离为H_c，则

s_c(k_x′，k_y′，k)＝s(k_x′，k_y′，k)exp(ik_zH_c) (5)

所述步骤(2.3)的步骤为：由色散关系式(4)可见，距离向波数k_z与k，k_x′和k_y′存在非线性关系。考虑到一般情况下，k，k_x′和k_y′是均匀分布的，因此采样点在k_z上是非均匀的。

根据匹配滤波原理，反射率函数的重建公式为

由公式(6)，σ(x，y，z)的重建公式从形式上符合IFFT，为适应后续的三维IFFT运算，有必要将采样点重采样至均匀分布的三维波数域(k_x′，k_y′，k_z)。其波数域分布是不均匀的，如果要采用FFT加速计算，必须先完成采样点的插值，使得s_c(k_x′，k_y′，k)在三维波数域(k_x′，k_y′，k_z)中均匀排布。假设重采样后回波表示为sc₍k_x′，k_y′，k_z)。

所述步骤(2.4)的步骤为：对重采样数据s_c(k_x′，k_y′，k_z)进行三维IFFT即可重建反射率初始图像，表示为

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

计算改进相干因子。首先给出三种反射率函数非相干功率的波数域新定义：

iv.

v.

vi.

三种定义方法在计算量上稍有差别，这是三维IFFT的计算顺序带来的。三维IFFT的计算过程相当于每个维度单独计算IFFT，因此允许使用中间量来减少计算量。本发明以定义i为主进行说明和展示实施例。改进相干因子定义为反射率函数相干功率与非相干功率的比值，即

所述步骤(3)的步骤为：利用步骤(2.5)计算得到的改进相干因子校正(2.4)得到的反射率函数初始图像，即

σ_m(x，y，z)＝σ(x，y，z)ICF(x，y，z) (9)

驻点扫描成像系统的工作方式如图1所示，天线1用于同时发射和接收电磁信号，天线扫描构成平面与XY平面重合，图中的网格点为对应天线驻留采样的位置，天线的扫描间隔需满足空间Nyquist采样定律，目标的散射点以(x，y，z)表示。基于改进相干因子的超分辨率驻点扫描实时成像算法流程图如图2所示。

一：点扩展函数

实施例中，采用工作频率75-110GHz，天线的空间采样间隔为3.5mm，在各维度上的扫描孔径尺寸为444.5mm，点目标在空间中位于(0，0，500)毫米处，采用RMA，CF-RMA和ICF-RMA三种方法得到的点扩展函数如图3所示。进一步，为了定量比较三种方法点扩展函数的性能，分析曲线的三种参量，即半峰全宽(FWHM)，背景噪声电平和最大旁瓣幅值，其具体定义为：

半峰全宽(FWHM)：靠近0dB处的两个负峰值被认为是主瓣的两端，主瓣幅值下降到3dB处对应的宽度是FWHM，该参数主要表示分辨率；

背景噪声电平：主瓣采样点之外，所有采样点的幅值平均值；

最大旁瓣幅值：主瓣外区域中采样点的最大幅值，其中最大旁瓣和背景噪声电平是图像动态范围和图像质量的指标。

表1点扩展函数的性能

表1对点扩展函数的三个参数性能进行计算，发现本发明提出的ICF-RMA相对于传统的RMA，分辨率提高了0.3mm，最大旁瓣和基底噪声电平均具有约10dB的压制。CF-RMA的结果比ICF-RMA稍好。

二：电磁仿真成像

该实施例主要用于验证和对比不同成像方法对连续目标体的成像性能。仿真中采用的工作频率与采样方式与实施例一是一致的。仿真中采用的目标体是扇形金属片，其尺寸如图4所示。电磁仿真采用矩量法，辐射源采用理想电偶极子。

分别采用ICF-RMA，RMA和CF-RMA三种方法对散射场进行目标图像重建，得到结果如图5-7所示。该图像是目标位置对应的三维切片图，即距离向500mm对应切片图。图像的动态范围去40dB。对比三幅图像结果能明显看出，RMA算法的成像结果最差，目标周围具有明显的旁瓣；而其他两种方法的结果图像质量相似。为进一步定量估计三种方法的成像质量，本实施例引用SSIM系数对三幅图像的结果进行定量估计，参考图像为图4目标生成的灰度图。得到SSIM系数结果如表2所示。理论上，SSIM系数越接近1，则目标图像与参考图像相似度越高，成像质量越好。由表2的结果，RMA的图像质量不及其他两种方法，这一点与我们的视觉结果是一致的。

表2 SSIM系数估计图像质量

三：实验成像结果

本实施例通过驻点扫描成像实验验证并比较三种成像方法的性能。成像实验采用；矢量网络分析仪(Agilent，N5247A)作为发射和相干接收信号，配合机械扫描台完成目标采样。目标为八叶金属片，如图8所示。表3给出了详细的实验参数。

表3实验参数

分别采用ICF-RMA，传统的RMA和CF-RMA三种方法的成像实验结果如图9-11所示。图像的动态范围是20dB，相比于仿真结果，动态范围存在恶化，这是由于采样仪器引入的噪声。由成像结果可以看出，ICF-RMA和CF-RMA两种方法对RMA的成像分辨率有明显提高，并且具有更干净的图像背景，证明两种超分辨方法都具有旁瓣基底噪声压制能力。

下面比较三种方法的计算效率。由于算法均采用复数运算，其中复数乘法运算占用绝大部分资源，因此采用复数乘法运算次数定量衡量算法的计算量。假设在X和Y方向上天线的空间采样点数分别为N_x和N_y，波数域采样点数为N_r。几种算法中各步骤的计算量列于表4中。

本实施例中提出采用改进相干因子与RMA计算量的比衡量由于改进相干因子的加入，ICF-RMA方法计算效率的下降程度：

其中，T_RMA和T_ICF分别表示RMA和改进相干因子的计算量。为了方便进一步讨论，下面简化为N_x＝N_y＝N_r＝N情况，得到近似比率

当N＝128时，由(11)计算得到α_est＝0.59，表面此时ICF的计算量大约相当于RMA计算量的0.59倍。但是对于CF-RMA，可以计算得到α_est＝1.03，即CF的计算量大致相当于一次额外的RMA计算。下面比较三种方法的程序实际运行时间，采用计算机的配置为：双核Intel(R)Xeon(R)CPU E5-2650，64GB RAM，运行时间列于表5中。从时间上看，ICF-RMA的运算量大约相当于RMA的1.54倍，相当于比率α＝0.54，该结果与理论估计结果一致性较好。而CF-RMA的计算量大约是RMA的两倍，由此在计算量上，本发明提出的ICF-RMA计算效率比其前身CF-RMA更高。

表4成像方法中各步骤的计算量

表5三种成像方法的实际运行时间

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (3)

1.一种超分辨率驻点扫描实时成像算法，其特征是：

步骤(1)：天线在空间平面上对目标等间距扫描采样，并在频域记录回波数据；

驻点扫描成像系统采用收发天线同置方式，在采样平面上，沿两个方向等间隔均匀采样，采样间隔满足空间Nyquist采样定理，在每一个采样位置上记录频域回波的宽带响应数据；

设采样平面在直角坐标系内与平面XY重合，设天线的坐标为(x′,y′,0)，则接收的回波信号s(x′,y′；k)表示为:

其中，R是从散射点到收发天线的单程距离；D(x,y,z)是散射点所在的目标域，σ(x,y,z)表示待重建的反射率函数；k是波数的幅值，c是电磁波在自由空间中的传播速度；

步骤(2)：采用距离偏移算法计算反射率函数的三维图像，在计算初始图像过程中、对插值后的回波进行三维逆快速傅里叶变换求取改进相干因子；

步骤(2.1)：对回波数据在空间域维度上做二维快速傅里叶变换，得到波数域回波信号；

对回波信号s(x′,y′,0；k)的两个空间维度x′和y′做二维快速傅里叶变换，得到波数域回波信号

s(k_x′,k_y′,k)＝∫∫∫_D(x,y,z)σ(x,y,z)exp(-uk_x′x-ik_y′y-ik_zz)dxdydz

其中，k_z是与距离向z方向对应的波数，设z>0；

步骤(2.2)：忽略沿传播路径上的信号衰减，根据已知的场景中心与采样平面距离，对步骤(2.1)得到的三维波数域回波信号进行场景中心补偿，使得目标的真实位置与成像域对齐；将目标域中心位置与目标距离对齐，设目标与天线的实际距离为H_c，则：

s_c(k_x′,k_y′,k)＝s(k_x′,k_y′,k)exp(ik_zH_c)；

步骤(2.3)：对步骤(2.2)得到的三维回波信号，根据色散关系式，在波数域维度上将回波数据插值到均匀分布的距离向波数域内，插值后的数据具有在三个维度上均匀采样的性质，这一步也称为重采样；

反射率函数的重建公式为：

将采样点重采样至均匀分布的三维波数域(k_x′,k_y′,k_z)；重采样后回波表示为s_c(k_x′,k_y′,k_z)；

步骤(2.4)：对波数域重采样的回波数据进行三维逆快速傅里叶变换，得到反射率函数的初始图像；

步骤(2.5)：根据反射率函数非相干功率在波数域的计算公式，对步骤(2.3)得到的回波数据求模平方，并求三维逆快速傅里叶变换；最后计算改进相干因子，定义为反射率函数相干功率与非相干功率的比；

所述的反射率函数非相干功率的波数域包括如下三种：

i.

ii.

iii.

步骤(3)：利用步骤(2)得到的改进相干因子校正初始图像；

利用步骤(2.5)计算得到的改进相干因子校正(2.4)得到的反射率函数初始图像。

2.根据权利要求1所述的超分辨率驻点扫描实时成像算法，其特征是步骤(2.4)中，对重采样数据s_c(k_x′,k_y′,k_z)进行三维逆快速傅里叶变换重建反射率初始图像，表示为：

3.根据权利要求1所述的超分辨率驻点扫描实时成像算法，其特征是：步骤(3)中，改进相干因子定义为反射率函数相干功率与非相干功率的比值，即

利用步骤(2.5)计算得到的改进相干因子校正(2.4)得到的反射率函数初始图像表示为：

σ_m(x,y,z)＝σ(x,y,z)ICF(x,y,z)。