CN106772370A - 基于多发多收阵列的太赫兹近场稀疏成像方法 - Google Patents

Abstract

本发明涉及一种基于多发多收阵列的太赫兹近场稀疏成像方法，经过成像场景设置、获取太赫兹回波、回波离散化、通过Akaike信息量准则估计出散射点数目、最后通过稀疏重构算法估计目标点坐标，实现多发多收阵列的太赫兹近场稀疏成像。所提方法保留了目标信息在路径延迟中的非线性，利用观测矩阵列原子序号与目标场景之间的一一映射关系，基于稀疏重构技术的贪婪思想实现目标重构。其中，为了改善经典贪婪算法对目标稀疏度（即散射点个数）的依赖性，所提算法在每次迭代选择原子序列时，利用Akaike信息量准则，基于余量对目标的稀疏度进行自适应估计，从而进一步提高了目标重构算法的普适性和实用性。

Description

基于多发多收阵列的太赫兹近场稀疏成像方法

技术领域

本发明涉及一种成像技术，特别涉及一种基于多发多收阵列的太赫兹近场稀疏成像方法。

背景技术

太赫兹电磁波具有能量低、相干测量等特性，这些特性使得太赫兹成像在很多领域具有应用价值；并且太赫兹波有很强的穿透能力，而且具有方向性，因此在一些环境下太赫兹成像具有很多优势。

众所周知，在基于合成孔径的微波成像理论中，距离分辨率和奈奎斯特采样率这两个量共同决定了成像系统的复杂度。距离向分辨率上限由发射信号的带宽决定，而根据奈奎斯特采样定律，这要求系统的采样频率至少要两倍的信号带宽。这代表着提高成像的距离分辨率，信号带宽和系统采样频率都必须相应提高，随之带来的是增加的系统设计复杂度和海量的待处理数据。同时，受到当下实际模拟/数字电子器件工艺发展水平的限制，增加信号带宽和系统采样频率都使得成像系统的设计与实现变得困难重重。因此，根据目标在高频电磁区域中可以由少数有限个强散射点表达这一目标具有稀疏先验的事实，引入稀疏信号处理的稀疏成像技术为突破传统微波成像理论的限制和简化系统复杂度带来了希望。

同时由于近场时，目标信息存在于路径延迟的非线性项中，传统基于傅里叶变换的重构方法都需要对路径延迟项进行目标空间位置的线性表达近似，如经典的单站系统中球面波的平面波叠加近似或驻定相位法分解等。尽管如此，这些方法对多站系统很不适用，继续沿用类似方法对多站的路径延迟进行线性近似，引入的近似误差大，造成重构性能的下降。此外，对于多站系统，若对路径延迟采用一阶泰勒级数展开后再基于傅里叶变换进行目标散射点重构则由于忽略了延迟高阶项的影响，目标被限定在很小的范围内，也不是近场时多站系统路径延迟处理方法的优化选择。所以基于稀疏重构技术的贪婪思想，保留目标信息在路径延迟中的非线性，成为了太赫兹多发多收阵列近场成像的合理选择。

发明内容

本发明是针对传统基于傅里叶变换的目标重构算法在太赫兹近场多站成像中存在近似误差增加和目标重构范围缩小的问题，提出了一种基于多发多收阵列的太赫兹近场稀疏成像方法，利用一种稀疏度自适应估计的稀疏重构算法解决传统基于傅里叶变换的目标重构算法在太赫兹近场多站成像中存在近似误差增加和目标重构范围缩小的问题，以及低样本率下的太赫兹成像低分辨率的问题，实现线性紧凑型太赫兹多发多收阵列近场稀疏成像。

本发明的技术方案为：一种基于多发多收阵列的太赫兹近场稀疏成像方法，具体包括如下步骤：

1)成像场景设置：太赫兹多发多收阵列的M个发射机和N个接收机分布在二维平面的同一基线上，同时发射机、接收机和目标平面三者都位于同一平面XY，以成像场景中心为直角坐标原点O建立坐标系XOY，则第m个发射机和第n个接收机所在位置分别表示为和设目标散射点的直角坐标为r＝(x,y)，散射系数记为σ(r)，发射机和接收机天线组成的天线阵列基线中心到目标平面中心的距离为R₀；

2)获取太赫兹回波，对回波进行去载波、滤波处理；

3)回波离散化：将获取的太赫兹回波根据成像场景离散化，构建观测矩阵，构建的观测矩阵中保留含有目标信息在延迟路径中的非线性，完整保留路径延迟项中的目标信息；

4)通过Akaike信息量准则估计出散射点数目；

5)通过贪婪算法估计目标点坐标，同时在每次迭代选择原子序列时，利用Akaike信息量准则，基于余量对目标的稀疏度进行自适应估计，实现多发多收阵列的太赫兹近场稀疏成像。

所述步骤2)具体步骤包括如下：

第m个发射机发射的信号：

其中u_m(t)是第m个发射机发射波形的复包络，其中m＝1,...,M，其带宽为B_m，f_m为第m个发射机的发射载频，发射频率分集信号，此时，任意两个相邻发射载频之间的差值足够大，使发射信号在频谱上占有分离不相交的频带，实现频域正交；

反射回波被第n个接收机接收的回波为：

其中n＝1,...,N，σ(x,y)是目标散射系数，τ_m,n(x,y)是第m个发射机到位于(x,y)的目标散射点再到第n个接收机的路径延迟；

该回波经去载波处理得到：

对回波利用发射信号频率正交性实现通道分离输出的第(n,m)个通道信号为：

并通过处理器得到：

其中表示卷积操作，上标*表示取共轭操作；

对做傅里叶变换，其频域形式为：

其中F{·}表示傅里叶变换；U_m(f)＝F{u_m(t)}，F{u_m(t)}是复包络u_m(t)的频谱，

设置低通匹配滤波器H_lp,m(f)为：

则经过频域滤波输出为：

所述步骤3)具体步骤包括如下：

a:将获取的太赫兹回波根据成像场景离散化：

将目标平面根据分辨率划分网格，对路径延迟项τ_m,n(x,y)进行展开：

其中(x_l,y_l)表示第l个强散射点的直角坐标，假设目标成像区间内共有L个待测的散射点，第m个发射天线到第l个散射点的距离为第n个接收天线到第l个散射的的距离为c为光速；

将获取的太赫兹回波根据目标平面划分的网格离散化处理得到：

其中表示第(n,m)个通道在第k个频率采样点f_k处的回波，表示第l个强散射点与第m个发射机和第n个接收机的距离，σ_l表示第l个强散射点的散射系数；

b:构建观测矩阵：

收集第(n,m)个通道的K个采样点构成向量为：

其中上标T为转置操作，

则收集MN个通道的回波可组成矩阵：

其中为MNK×L维的观测矩阵，σ为目标场景对应的L×1维向量，实际目标强散射点的个数为L₀，满足L₀＜＜L。

所述步骤4)具体步骤包括如下：

得到接收回波的协方差矩阵R＝S^r×(S^r)^H，

其中，上标H表示共轭转置；

Akaike信息量准则估计算法：

为似然函数

其中，λ_i为R的特征值，得到的散射点数目为0,1,...,MNK-1之间使Akaike信息量准则取最小值的估计值

所述步骤5)选择贪婪算法，并每次迭代时基于余量对目标的稀疏度进行自适应估计，以观测矩阵A和接收回波S^r作为输入，目标向量的估计值作为输出，目标点坐标进行估计。

本发明的有益效果在于：本发明基于多发多收阵列的太赫兹近场稀疏成像方法，所提算法保留了目标信息在路径延迟中的非线性，利用观测矩阵列原子序号与目标场景之间的一一映射关系，基于稀疏重构技术的贪婪思想实现目标重构。其中，为了改善经典贪婪算法对目标稀疏度(即散射点个数)的依赖性，所提算法在每次迭代选择原子序列时，利用Akaike信息量准则，基于余量对目标的稀疏度进行自适应估计，从而进一步提高了目标重构算法的普适性和实用性。

附图说明

图1为本发明基于多发多收阵列的太赫兹近场稀疏成像方法流程图；

图2为采用基于稀疏重构算法的多发多收阵列太赫兹近场成像方法的二维成像模型图；

图3为本发明基于多发多收阵列的太赫兹近场稀疏成像结果图；

图4为采用一阶泰勒级数展开后再基于傅里叶变换算法的多发多收阵列太赫兹近场成像结果图。

具体实施方式

如图1所示本发明基于多发多收阵列的太赫兹近场稀疏成像方法流程图，具体实施方案如下：

步骤A：成像场景设置；

采用线性紧凑型太赫兹多发多收阵列。图2所示为二维成像模型图，线性紧凑型太赫兹多发多收阵列的M个发射机和N个接收机分布在二维平面的同一基线上，同时发射机、接收机和目标平面三者都位于同一平面XY。以成像场景中心为直角坐标原点o建立坐标系XOY，则第m个发射机和第n个接收机所在位置分别表示为和设目标散射点的直角坐标为r＝(x,y)，散射系数记为σ(r)，天线阵列(发射机和接收机天线组成的阵列)基线中心到目标平面中心的距离为R₀。

步骤B：获取太赫兹回波；

第m个发射机发射的信号：

其中u_m(t)是第m个发射机发射波形的复包络，其中m＝1,...,M，其带宽为B_m，f_m为第m个发射机的发射载频，发射频率分集信号，此时，任意两个相邻发射载频之间的差值足够大，使发射信号在频谱上占有分离不相交的频带，实现频域正交。

反射回波被第n个接收机接收的回波为：

其中n＝1,...,N，σ(x,y)是目标散射系数，τ_m,n(x,y)是第m个发射机到位于(x,y)的目标散射点再到第n个接收机的路径延迟。

该回波经去载波处理得到：

对回波利用发射信号频率正交性实现通道分离输出的第(n,m)个通道信号为：

并通过相关处理器得到：

其中表示卷积操作，上标*表示取共轭操作。

对做傅里叶变换，其频域形式为：

其中F{·}表示傅里叶变换；U_m(f)＝F{u_m(t)}，F{u_m(t)}是复包络u_m(t)的频谱，

设置低通匹配滤波器H_lp,m(f)为：

则经过频域滤波输出为：

步骤C：回波离散化；

子步骤C1：将获取的太赫兹回波根据成像场景离散化：

将目标平面根据分辨率划分网格，对路径延迟项τ_m,n(x,y)进行展开：

其中(x_l,y_l)表示第l个强散射点(假设目标成像区间内共有L个待测的散射点)的直角坐标，第m个发射天线到第l个散射点的距离为第n个接收天线到第l个散射的距离为c为光速。

将获取的太赫兹回波根据目标平面划分的网格离散化处理得到：

其中表示第(n,m)个通道在第k个频率采样点f_k处的回波，表示第l个强散射点与第m个发射机和第n个接收机的距离，σ_l表示第l个强散射点的散射系数。

子步骤C2：构建观测矩阵：

收集第(n,m)个通道的K个采样点构成向量为：

其中上标T为转置操作，

则收集MN个通道的回波可组成矩阵：

其中为MNK×L维的观测矩阵，σ为目标场景对应的L×1维向量，实际目标强散射点的个数为L₀，满足L₀＜＜L。

在构建观测矩阵时保留了含有目标信息的延迟路径非线性项，未做任何近似及分解，能更完整保留相位非线性项中的目标信息。

步骤D：通过Akaike信息量准则估计出散射点数目；

子步骤D1：得到接收回波的协方差矩阵R

R＝S^r×(S^r)^H

其中，上标H表示共轭转置。

子步骤D2：Akaike信息量准则(Akaike information criterion、简称AIC)估计算法：

为似然函数

其中，λ_i为R的特征值，得到的散射点数目为0,1,...,MNK-1之间使Akaike信息量准则取最小值的估计值

步骤E：通过稀疏重构算法估计出目标点坐标，实现太赫兹近场成像；

现有的稀疏重构算法主要可以分为三类：凸优化方法，贪婪算法和组合算法。这些稀疏重构算法都可以应用于本发明中，以实现多发多收阵列的太赫兹近场稀疏成像。在本发明中借鉴贪婪算法的基本思想，在每次迭代中利用通过Akaike信息量准则估计余量的稀疏度(即散射点个数)，再根据此稀疏度更新搜索步长，提出一种稀疏度自适应估计的稀疏重构算法，以实现较好的稀疏成像效果。

以下算法中，p指代外循环迭代次数，z指代内循环迭代次数；为目标向量的估计值；S^r为接收回波；A为观测矩阵；J_p为第p次外循环迭代更新的支撑集，是第p次迭代时对目标场景估计出的散射点坐标序号的集合；r_p为第p次迭代后的余量；C_p为第p次迭代的支撑集候选集；G_p为第p次迭代的索引集；w_p为第p次迭代的搜索步长；r_p，z和J_p，z分别为第p次外循环下第z次内循环得到的候补余量和候补支撑集，当满足一定条件时赋值给余量r_p和支撑集J_p；为第p次迭代时候补余量r_p，z中散射点数目估计值；I_p，z为第p次外循环下第z次内循环的搜索步长；表示向上取整；|·|表示取绝对值；||a||₂表示a的2-范数；表示广义逆矩阵；A_C表示取矩阵A中下标为集合C中元素确定的列向量组成的子矩阵；∪表示取并集；表示空集。

具体流程为：

输入：观测矩阵A，接收回波S^r

输出：目标向量的估计值

子步骤E1：根据步骤D中估计出的散射点数目得到步长的初始值，令外循环初始搜索步长

子步骤E2：初始化向量余量r₀＝S^r，支撑集索引集候补支撑集令外循环迭代次数p＝1，内循环迭代次数z＝1。

子步骤E3：更新第p次外循环下第z次内循环的搜索步长为I_p,z＝w_p-1；计算|Ar_p-1|，将结果中模值最大的I_p，z个元素值所对应的下标存入索引集G_p。

子步骤E4：令支撑集候选集C_p＝J_p-1∪G_p；计算将结果中模值最大的I_p，z个元素值所对应的下标存入候补支撑集J_p，z；计算候补余量

子步骤E5：利用Akaike信息量准则估计出余量r_p，z中的散射点数目，估计值为步长更新

子步骤E6：如果满足迭代终止条件，则结束迭代，输出结果；否则，转至子步骤E7。

其中迭代终止条件为：当余量r_p，z的2-范数小于阈值时迭代终止，即||r_p,z||₂≤ε，其中在无噪声情况下ε＝0，有噪声情况下ε可选择为噪声的能量值。

子步骤E7：如果||r_p,z||₂＞||r_p-1||₂，则需要进行内循环，更新内循环迭代次数z＝z+1，支撑集的大小I_p,z＝zw_p，使得支撑集的大小越来越接近未知的稀疏度(即散射点个数L₀)，转至子步骤E3；否则，更新支撑集J_p＝J_p,z、余量r_p＝r_p,z和外循环迭代次数p＝p+1，内循环迭代次数重置z＝1，转至子步骤E3。

图3是在本实施例基于稀疏重构算法的多发多收阵列太赫兹近场成像方法成像结果图，图4是基于泰勒展开与傅里叶变换算法的多发多收阵列太赫兹近场成像结果图。其中：载波频率180GHz；采用7GHz带宽；发射机接收机所在基线与目标平面距离为5米；并根据分辨率划分目标平面网格；设置4个发射机，10个接收机均分布在二维平面的同一基线上，同时发射机、接收机和目标三者都位于同一平面，每相邻两个发射机间的距离均为0.1667米，每相邻两个接收机间的距离均为0.0167米。图3和4中白色圆点指示真实的散射点位置；4个发射机和10个接收机共组成40组通道，图3中示出在40组通道中随机选取8组通道，采用稀疏度自适应估计的稀疏重构算法得出的太赫兹近场成像结果图；图4中示出40组通道全部选取，采用一阶泰勒级数展开后再基于傅里叶变换算法得出的太赫兹近场成像结果图。

比较上述两维点目标仿真实验的成像结果，知此稀疏度自适应估计的稀疏重构算法在稀疏采样情况下，仍能较好的还原目标信息。并且本发明中利用的稀疏重构算法能够较好的解决多发多收阵列在近场目标信息位于延迟路径非线性项的问题，没有出现利用经典单站系统中球面波的平面波叠加近似或驻定相位法分解等引起的误差较大的问题，也没有出现采用一阶泰勒级数展开后再基于傅里叶变换引起的目标被限定在很小范围内的问题，所以本发明提出的算法具有普适性和实用性。

Claims (5)

1.一种基于多发多收阵列的太赫兹近场稀疏成像方法，其特征在于，具体包括如下步骤：

1)成像场景设置：太赫兹多发多收阵列的M个发射机和N个接收机分布在二维平面的同一基线上，同时发射机、接收机和目标平面三者都位于同一平面XY，以成像场景中心为直角坐标原点O建立坐标系XOY，则第m个发射机和第n个接收机所在位置分别表示为和设目标散射点的直角坐标为r＝(x,y)，散射系数记为σ(r)，发射机和接收机天线组成的天线阵列基线中心到目标平面中心的距离为R₀；

2)获取太赫兹回波，对回波进行去载波、滤波处理；

3)回波离散化：将获取的太赫兹回波根据成像场景离散化，构建观测矩阵，构建的观测矩阵中保留含有目标信息在延迟路径中的非线性，完整保留路径延迟项中的目标信息；

4)通过Akaike信息量准则估计出散射点数目；

5)通过贪婪算法估计目标点坐标，同时在每次迭代选择原子序列时，利用Akaike信息量准则，基于余量对目标的稀疏度进行自适应估计，实现多发多收阵列的太赫兹近场稀疏成像。

2.根据权利要求1所述基于多发多收阵列的太赫兹近场稀疏成像方法，其特征在于，所述步骤2)具体步骤包括如下：

第m个发射机发射的信号：

其中u_m(t)是第m个发射机发射波形的复包络，其中m＝1,...,M，其带宽为B_m，f_m为第m个发射机的发射载频，发射频率分集信号，此时，任意两个相邻发射载频之间的差值足够大，使发射信号在频谱上占有分离不相交的频带，实现频域正交；

反射回波被第n个接收机接收的回波为：

$s_n^r\left(t\right)=\underset{m=1}{\overset{M}{\Σ}}\underset{\left(x,y\right)}{\∫\∫}\mathrm{\σ}(x,y)u_m(t-{\mathrm{\τ}}_{m,n}(x,y\left)\right)e^{j2{\mathrm{\πf}}_m(t-{\mathrm{\τ}}_{m,n}(x,y\left)\right)}dxdy$

其中n＝1,...,N，σ(x,y)是目标散射系数，τ_m,n(x,y)是第m个发射机到位于(x,y)的目标散射点再到第n个接收机的路径延迟；

该回波经去载波处理得到：

$s_n^r\left(t\right)=s_n^r\left(t\right)\·e^{-j2{\mathrm{\πf}}_{m}t}=\underset{m=1}{\overset{M}{\Σ}}\underset{\left(x,y\right)}{\∫\∫}\mathrm{\σ}(x,y)u_m(t-{\mathrm{\τ}}_{m,n}(x,y\left)\right)e^{-j2{\mathrm{\πf}}_m{\mathrm{\τ}}_{m,n}(x,y)}dxdy$

对回波利用发射信号频率正交性实现通道分离输出的第(n,m)个通道信号为：

${\tilde{s}}_{m,n}^r\left(t\right)=\underset{(x,y)}{\∫\∫}\mathrm{\σ}(x,y)u_m(t-{\mathrm{\τ}}_{m,n}(x,y\left)\right)e^{-j2{\mathrm{\πf}}_m{\mathrm{\τ}}_{m,n}(x,y)}dxdy$

并通过处理器得到：

$s_{m,n}^r\left(t\right)={\tilde{s}}_{m,n}^r\left(t\right)\⊗u_m^{*}\left(t\right)=\underset{\left(x,y\right)}{\∫\∫}\mathrm{\σ}(x,y)u_m(t-{\mathrm{\τ}}_{m,n}(x,y\left)\right)e^{-j2{\mathrm{\πf}}_m{\mathrm{\τ}}_{m,n}(x,y)}dxdy\⊗u_m^{*}\left(t\right)$

其中表示卷积操作，上标*表示取共轭操作；

对做傅里叶变换，其频域形式为：

$F\left\{s_{m,n}^r\left(t\right)\right\}=\underset{\left(x,y\right)}{\∫\∫}\mathrm{\σ}(x,y)U_m\left(f\right)U_m^{*}\left(f\right)e^{-j2\mathrm{\π}(f_m+f){\mathrm{\τ}}_{m,n}(x,y)}dxdy$

其中F{·}表示傅里叶变换；U_m(f)＝F{u_m(t)}，F{u_m(t)}是复包络u_m(t)的频谱，

设置低通匹配滤波器H_lp,m(f)为：

则经过频域滤波输出为：

$S_{m,n}^r\left(f\right)=F\left\{s_{m,n}^r\left(t\right)\right\}{H}_{lp,m}\left(f\right)=\underset{\left(x,y\right)}{\∫\∫}\mathrm{\σ}(x,y)e^{-j2\mathrm{\π}(f_m+f){\mathrm{\τ}}_{m,n}(x,y)}dxdy.$

3.根据权利要求2所述基于多发多收阵列的太赫兹近场稀疏成像方法，其特征在于，所述步骤3)具体步骤包括如下：

a:将获取的太赫兹回波根据成像场景离散化：

将目标平面根据分辨率划分网格，对路径延迟项τ_m,n(x,y)进行展开：

$\begin{array}{c}{\mathrm{\τ}}_{m,n}\left(x_l,y_l\right)=\frac{R_{Tx,m}^l\left(x_l,y_l\right)+R_{Rx,n}^l\left(x_l,y_l\right)}{c}\\ =\frac{\sqrt{{\left(x_m^t-x_l\right)}^2+{\left(y_m^t-y_l\right)}^2}}{c}+\frac{\sqrt{{\left(x_n^r-x_l\right)}^2+{\left(y_n^r-y_l\right)}^2}}{c}\end{array}$

其中(x_l,y_l)表示第l个强散射点的直角坐标，假设目标成像区间内共有L个待测的散射点，第m个发射天线到第l个散射点的距离为第n个接收天线到第l个散射的的距离为c为光速；

将获取的太赫兹回波根据目标平面划分的网格离散化处理得到：

$S_{m,n}^r\left(f_k\right)=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\σ}}_l{e}^{-j2\mathrm{\π}(f_m+f_k){\mathrm{\τ}}_{m,n}(x,y)}=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\σ}}_l{e}^{-j2\mathrm{\π}\frac{(f_m+f_k)}{c}{R}_{m,n}\left(l\right)}$

其中表示第(n,m)个通道在第k个频率采样点f_k处的回波，表示第l个强散射点与第m个发射机和第n个接收机的距离，σ_l表示第l个强散射点的散射系数；

b:构建观测矩阵：

收集第(n,m)个通道的K个采样点构成向量为：

$S_{m,n}^r={\[S_{m,n}^r\left(f_1\right)...S_{m,n}^r\left(f_K\right)\]}^T\stackrel{\mathrm{\Δ}}{=}{A}_{m,n}\mathrm{\σ}$

其中上标T为转置操作，

则收集MN个通道的回波可组成矩阵：

$S^r={\[{\left(S_{1,1}^r\right)}^T,...,{\left(S_{M,N}^r\right)}^T\]}^T\stackrel{\mathrm{\Δ}}{=}A\mathrm{\σ}$

其中为MNK×L维的观测矩阵，σ为目标场景对应的L×1维向量，实际目标强散射点的个数为L₀，满足L₀＜＜L。

4.根据权利要求3所述基于多发多收阵列的太赫兹近场稀疏成像方法，其特征在于，所述步骤4)具体步骤包括如下：

得到接收回波的协方差矩阵R＝S^r×(S^r)^H，

其中，上标H表示共轭转置；

Akaike信息量准则估计算法：

$AIC\left({\hat{L}}_0\right)=2K(MNK-{\hat{L}}_0)\ln \mathrm{\Λ}\left({\hat{L}}_0\right)+2{\hat{L}}_0(2MNK-{\hat{L}}_0)$

为似然函数

其中，λ_i为R的特征值，得到的散射点数目为0,1,...,MNK-1之间使Akaike信息量准则取最小值的估计值

5.根据权利要求4所述基于多发多收阵列的太赫兹近场稀疏成像方法，其特征在于，所述步骤5)选择贪婪算法，并每次迭代时基于余量对目标的稀疏度进行自适应估计，以观测矩阵A和接收回波S^r作为输入，目标向量的估计值作为输出，目标点坐标进行估计。