CN104898118A - 一种基于稀疏频点的三维全息成像的重建方法 - Google Patents

Abstract

本发明提供一种基于稀疏频点的微波、毫米波、太赫兹三维全息成像的重建方法，属于图像处理技术领域。该方法通过对微波、毫米波、太赫兹雷达原始稀疏频点及其个数的选择，稀疏频点数据重建满频点数据，满频点数据处理得到的微波、毫米波、太赫兹三维全息成像结果。本发明方法能有效消除距离模糊现象，降低对雷达系统采样率和储存深度的要求，适合于实时应用成像中。

Description

一种基于稀疏频点的三维全息成像的重建方法

技术领域

本发明属于图像处理技术领域，尤其涉及一种基于稀疏频点的三维全息成像的重建方法，如微波、毫米波、太赫兹的三维全息成像。

背景技术

近年来，利用微波、毫米波、太赫兹波作为探测手段，获取探测目标的散射信息的微波、毫米波、太赫兹三维全息成像技术越来越受到人们的关注。微波、毫米波、太赫兹三维全息成像是在两个正交方向形成扫描孔径，并在距离向上发射宽带信号来形成三维成像几何，通过处理回波数据得到目标的三维复图像。微波、毫米波、太赫兹波可以穿透很多非极性材料，如皮革、塑料、泡沫等。通过对被检测物体进行高分辨率穿透成像，可以获得更准确的目标信息，大大降低误警率。微波、毫米波、太赫兹波与X射线相比不会因为光致电离而破坏被检测的物质，对人体来说是基本无害，具有较好的安全性。因此微波、毫米波、太赫兹三维全息成像在公共区域的安全检测、无损检测、雷达目标成像诊断等领域有广泛的应用前景。

随着人们对全息成像分辨率的要求越来越高，在对距离分布范围比较广的目标进行高分辨率全息成像时，为了避免距离模糊现象，需增加频域的采样点数。频域采样点数的增加会提高对雷达系统的数据采样率要求，同时雷达回波数据量变大，对雷达系统的储存深度也有了更高要求。这些指标要求的提高会增加雷达系统硬件的成本，同时高频点数的数据采集也增加了扫描时间，这对实时成像领域的应用是非常不利的。为了在不增加硬件成本的情况下消除距离模糊，可以通过一定算法来进行基于稀疏频点数据的无距离模糊的三维全息成像。一般来说，微波、毫米波、太赫兹波的穿透能力是有限的，其对于金属目标和人体皮肤来说是无法穿透的，因此雷达目标的后向散射在距离域表现出了稀疏特性。压缩感知理论表明，对于这种稀疏可压缩的信号，可以通过构建满足一定条件的观测矩阵，从较少的观测数据中恢复原信号。基于此理论，在雷达系统数据采集时，在频域进行对应于构造的观测矩阵的非均匀采样，利用压缩感知理论的稀疏信号恢复方法从采样得到的稀疏频点数据恢复出满频点的回波数据，从而消除了距离模糊现象。在实现本发明的过程中，发现上述现有技术存在如下技术缺陷：

1.在利用压缩感知理论进行稀疏频点数据恢复出满频点的回波数据的过程中，采样频点的选择需要满足对应于构造的观测矩阵的要求，如果获得的采样频点不满足观测矩阵的要求，会造成压缩感知理论的方法失效，因此对实际数据的采样获取造成很大的不便。

2.从采样得到的稀疏频点数据恢复出满频点的回波数据的压缩感知稀疏信号恢复方法，主要包括以正交匹配追踪算法为代表的贪婪追踪算法和以两步迭代收缩为代表的凸优化算法。这些稀疏信号恢复方法是通过多次迭代求解来找到对原始信号的最优逼近，如果迭代运算的次数过少，则会造成得到的信号结果偏离最优值，如果迭代运算的次数过多，则会造成运算量增加，处理时间过长，从而导致算法的运算效率很低，不利于其在实时成像中的应用。

并且迭代算法对噪声的敏感度非常高，噪声对信号的影响比较大，使得信号形式发生改变，从而求得的结果往往是某个局部最优解，而不是全局最优解，最终导致运算结果不准确，恢复图像散焦。

发明内容

为解决上述问题，本发明提供一种基于稀疏频点的微波、毫米波、太赫兹三维全息成像的重建方法，能快速精确恢复出满频点的回波数据，消除距离模糊现象，降低对雷达系统采样率和储存深度的要求，处理过程简单快速，适用于实时应用成像中。

本发明的基于稀疏频点的微波、毫米波、太赫兹三维全息成像的重建方法，适用于微波、毫米波、太赫兹波的雷达发射信号，其包括：

S1：确定雷达发射信号进行采样的稀疏频点个数N，所述稀疏频点个数N满足：

$N\≥\max \{(B+1+\frac{2}{B})r,\frac{2B_f\mathrm{\ΔR}}{c}\}---\left(1\right)$

其中，B为设定中间量，且B＝lnM-ln(1-β)，M为待恢复的频点个数，即满频点的个数；β为最大可忍受的峰值旁瓣比r的置信系数；B_f为雷达发射信号的带宽，c为电磁波在自由空间的传播速度，△R为在雷达天线有效波束角的照射范围内的重建目标表面在距离向最大分布范围长度；

S2：如果雷达发射信号属于线性调频连续波信号则先补偿掉该发射方式导致的剩余视频相位，消除剩余视频相位后得到和步进频连续波信号相同的回波信号后，再进行以下处理；如果雷达发射信号属于步进频连续波信号，则直接雷达发射信号进行以下处理：

在带宽内选取均匀随机分布的N个频点f₁,f₂,…,f_N，定义稀疏频点矩阵f_sparse＝(f₁,f₂,…,f_N)；

采样得到N个稀疏频点的回波数据S₁(x,y,k_sparse)，建立坐标系XYZ；其中定义X、Y为扫描方向，Z为雷达照射方向，(x,y,0)为在坐标系XYZ中的天线相位中心位置，k_sparse为f_sparse对应的自由空间的波数；

利用式(2)获得回波信号S₁(x,y,k_sparse)的转置

$S_1^T(x,y,k_{\mathrm{sparse}})=[s_1(x,y,k_1),s_1(x,y,k_2),...,s_1(x,y,k_N)]---\left(2\right)$

其中， $s_1(x,y,k_m)=\∫\∫\∫\mathrm{\σ}(x^{\′},y^{\′},z^{\′}){\mathrm{\α}}_{\mathrm{\θ}}\left[{\tan }^{-1}\frac{\mathrm{\ρ}(x^{\′},y^{\′})}{z^{\′}}\right]\exp (-2j{k}_{m}R){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}{\mathrm{dz}}^{\′}---\left(3\right)$

m＝1,2…N，σ(x′,y′,z′)为在坐标系XYZ中位于(x′,y′,z′)处目标的散射系数，a_θ(·)为雷达波束角的窗函数， $R=\sqrt{{(x-x^{\′})}^2+{(y-y^{\′})}^2+{(z-z^{\′})}^2};$

S3：对S2所得与设定的参考信号利用式(4)求相关：

$\mathrm{Corr}\left(R_j^{\′}\right)=S_1^T(x,y,k_{\mathrm{sparse}})S_2^{*}(x,y,k_{\mathrm{sparse}},R_j^{\′})---\left(4\right)$

其中， $S_2^T(x,y,k_{\mathrm{sparse}},R_j^{\′})=[\exp (-2j{k}_1{R}_j^{\′}),\exp (-2j{k}_2{R}_j^{\′}),...,\exp (-2j{k}_N{R}_j^{\′})],$ $R_j^{\′}=\frac{\mathrm{jc}}{2B_f},j=\mathrm{0,1},...,M-1;$ 为的共轭矩阵；

将Corr(R′_j)结果中的最大值max[Corr(R′_j)]对应的距离R′_j作为在雷达天线有效波束角的照射范围内的目标表面在距离向估计值，记为距离R′_max；

S4：对S2所得回波数据进行频谱搬移，利用式(5)结合S3获得的距离R′_max将其频谱搬移至基带，得到

式中 $S_3^T\left(k_{\mathrm{sparse}}\right)=[\exp \left(2j{k}_1{R}_{\max }^{\′}\right),\exp \left(2j{k}_2{R}_{\max }^{\′}\right),...,\exp \left(2j{k}_N{R}_{\max }^{\′}\right)],$ ο表示为矩阵的Hadamard积；

S5：对S4所得结果利用式(6)对波数进行线性插值得到

$S_5^T(x,y,k)={\mathrm{interp}}_{k_{\mathrm{sparse}}}\left[S_4^T(x,y,k_{\mathrm{sparse}})\right]---\left(6\right)$

式中k＝(k′₁,k′₂,…,k′_M)为满频点的各频率对应的自由空间的波数，其中k′_i＝2π[f_c-B_f/2+(i-1)B_f/(M-1)]/c,i＝1,2,…,M，f_c为雷达发射信号的中心频率，M为满频点的个数；interp_X[Y(X)]为插值函数；

S6：对S5所得结果进行频谱搬移，结合S3获得的距离R′_max利用式(7)将其频谱的中心频率搬移至其实际所处的频率，得到

式中 $S_7^T\left(k\right)=[\exp \left(2{\mathrm{jk}}_1{R}_{\max }^{\′}\right),\exp \left(2{\mathrm{jk}}_2{R}_{\max }^{\′}\right),...,\exp \left(2{\mathrm{jk}}_M{R}_{\max }^{\′}\right)];$

S7：利用相位偏移算法对S6所得结果进行图像处理，得到恢复图像。

有益效果：

该方法通过对微波、毫米波、太赫兹雷达原始稀疏频点及其个数的选择，重建满频点的回波数据，满频点数据处理得到的微波、毫米波、太赫兹三维全息成像结果。在高运算效率的情况下，可以快速的处理稀疏频点的微波、毫米波、太赫兹雷达三维数据进行高分辨率快速成像，有效的消除了距离模糊现象，降低对雷达系统采样率和储存深度的要求。

1.步骤2中频点选取均匀随机分布，不存在限制条件，实际数据的采样获取简单可靠。

2.本方法不存在迭代过程，而是通过步骤3的方式直接找到距离最优估计值R′_max。

附图说明

图1为本发明的坐标系XYZ示意图；

图2为本发明的成像目标的光学图像；

图3为本发明的200频点数据处理结果示意图；

图4为本发明的目标IECAS 200频点数据距离向恢复结果示意图；

图5为本发明的目标剪刀200频点数据距离向恢复结果示意图；

图6为本发明的目标四个三角形200频点数据距离向恢复结果示意图；

图7为本发明的20频点数据处理结果示意图；

图8为本发明的目标IECAS 20频点数据距离向恢复结果示意图；

图9为本发明的目标剪刀20频点数据距离向恢复结果示意图；

图10为本发明的目标四个三角形20频点数据距离向恢复结果示意图。

具体实施方式

本发明的目的是提出一种用于稀疏频点的微波、毫米波、太赫兹三维全息成像中的重建方法，该方法通过对微波、毫米波、太赫兹雷达原始稀疏频点及其个数的选择，目标距离预估，稀疏频点回波信号的频谱搬移，稀疏频点数据重建满频点数据，满频点信号的频谱搬移得到满频点的回波数据。

为了达到上述目的，本发明的技术解决方案是：

步骤S1：确定雷达发射信号的采样频点个数，雷达发射信号包括：微波、毫米波、太赫兹的。根据采样稀疏频点个数与峰值旁瓣比的估计关系，利用统计学原理可得采样稀疏频点个数N需满足：

$N\≥\max \{(B+1+\frac{2}{B})r,\frac{2B_f\mathrm{\ΔR}}{c}\}---\left(1\right)$

其中B为设定中间量，且B＝lnM-ln(1-β)，M为待恢复的频点个数，β为最大可忍受的峰值旁瓣比r的置信系数；r为最大可忍受的峰值旁瓣比，B_f为发射信号带宽，c为电磁波在自由空间的传播速度，△R为在雷达天线有效波束角的照射范围内的重建目标表面在距离向最大分布范围长度。

步骤S2：如果雷达发射信号为线性调频连续波信号，则需要先补偿掉该发射方式导致的剩余视频相位，消除剩余视频相位后得到和步进频连续波信号相同的的回波信号后，再进行如下处理；如果雷达发射信号为步进频连续波信号，直接进行以下处理：

在带宽内选取均匀随机分布的N个频点f₁,f₂,…,f_N，定义稀疏频点矩阵f_sparse＝(f₁,f₂,…,f_N)；

采样得到这N个频点的回波数据S₁(x,y,k_sparse)，如图1所示，建立坐标系XYZ，其中定义X、Y为扫描方向，Z为雷达照射方向，k_sparse＝(k₁,k₂,…,k_N)为各频点f₁,f₂,…,f_N对应的自由空间的波数，利用式(2)获得回波信号的转置可表示为：

$S_1^T(x,y,k_{\mathrm{sparse}})=[s_1(x,y,k_1),s_1(x,y,k_2),...,s_1(x,y,k_N)]---\left(2\right)$

(2)式中

$s_1(x,y,k_m)=\∫\∫\∫\mathrm{\σ}(x^{\′},y^{\′},z^{\′}){\mathrm{\α}}_{\mathrm{\θ}}\left[{\tan }^{-1}\frac{\mathrm{\ρ}(x^{\′},y^{\′})}{z^{\′}}\right]\exp (-2j{k}_{m}R){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}{\mathrm{dz}}^{\′}---\left(3\right)$

式中m＝1,2…N，σ(x′,y′,z′)为位于(x′,y′,z′)处目标的散射系数，a_θ(·)为雷达波束角的窗函数， $\mathrm{\ρ}(x^{\′},y^{\′})=\sqrt{{(x-x^{\′})}^2+{(y-y^{\′})}^2},R=\sqrt{{(x-x^{\′})}^2+{(y-y^{\′})}^2+{(z-z^{\′})}^2}.$

步骤S3：为了得到在雷达天线有效波束角范围内目标表面在距离向的估计值，对步骤S2所得回波数据求相关得到目标表面的距离向分布。对步骤S2所得回波数据与参考信号 $S_2^T(x,y,k_{\mathrm{sparse}},R_j^{\′})=[\exp (-2j{k}_1{R}_j^{\′}),\exp (-2j{k}_2{R}_j^{\′}),...,\exp (-2j{k}_N{R}_j^{\′})],$ $R_j^{\′}=\frac{\mathrm{jc}}{2B_f},j=\mathrm{0,1},...,M-1$ 求相关：

$\mathrm{Corr}\left(R_j^{\′}\right)=S_1^T(x,y,k_{\mathrm{sparse}})S_2^{*}(x,y,k_{\mathrm{sparse}},R_j^{\′})---\left(4\right)$

Corr(R′_j)结果中的最大值max[Corr(R′_j)]对应的距离R′_max为在雷达天线有效波束角的照射范围内的目标表面在距离向估计值。。

步骤S4：根据奈奎斯特采样定理，为使采样得到的原始稀疏频点的回波数据能正确恢复满频点数据，需对非基带的原始稀疏频点的回波数据搬移至基带。对原始稀疏频点的回波数据进行频谱搬移，将其频谱搬移至基带。

式中 $S_3^T\left(k_{\mathrm{sparse}}\right)=[\exp \left(2j{k}_1{R}_{\max }^{\′}\right),\exp \left(2j{k}_2{R}_{\max }^{\′}\right),...,\exp \left(2j{k}_N{R}_{\max }^{\′}\right)],$ ο表示为矩阵的Hadamard积。

步骤S5：对满足奈奎斯特采样定理的基带信号进行满频点数据重建。对步骤S4所得结果对波数进行线性插值得：

$S_5^T(x,y,k)={\mathrm{interp}}_{k_{\mathrm{sparse}}}\left[S_4^T(x,y,k_{\mathrm{sparse}})\right]---\left(6\right)$

式中k＝(k′₁,k′₂,…,k′_M)为满频点的各频率对应的自由空间的波数，其中k′_i＝2π[f_c-B_f/2+(i-1)B_f/(M-1)]/c,i＝1,2,…,M，f_c为雷达发射信号的中心频率，M为满频点的个数。

步骤S6：对步骤S5所得结果进行频谱搬移，将其频谱的中心频率搬移至其实际所处的频率。

式中 $S_7^T\left(k\right)=[\exp \left(2{\mathrm{jk}}_1{R}_{\max }^{\′}\right),\exp \left(2{\mathrm{jk}}_2{R}_{\max }^{\′}\right),...,\exp \left(2{\mathrm{jk}}_M{R}_{\max }^{\′}\right)].$

步骤S7：利用相位偏移算法对步骤S6所得结果进行处理，得到恢复图像。

此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：

(1)雷达发射信号为线性调频连续波信号，可对回波数据进行均匀随机分布的采样，得到N个频点的回波数据；

(2)对步骤S3中的求相关操作，可利用快速傅里叶变换进行计算；

(3)步骤S5的线性插值操作可以通过其他插值操作来实现，例如三阶样条插值，多项式插值法；

(4)本方法也可以用于两个正交方向采用合成孔径与实孔径相结合的成像体制，以利于快速的数据采集和成像处理。

下面通过实测数据对该算法进行了验证。逐点扫描三维全息成像雷达系统参数如表1所示：

表1

图2为待成像目标的光学图像。其中三个目标分别位于三个不同的距离平面上，这三个距离平面的间隔为25cm。图3为满频点200个频点数据恢复结果，可以看出不同距离面上的目标都得到了很好的聚焦。图4为目标IECAS通过200频点数据得到的距离向恢复结果，图5为目标剪刀通过200频点数据得到的距离向恢复结果，图6为目标四个三角形通过200频点数据得到的距离向恢复结果。通过这三幅图可以看出满频点的数据能够精确重建目标所在距离位置的分布，无距离模糊现象出现。图7为利用本方法对20频点数据处理得到的结果，不同距离面上的目标也都得到了很好的聚焦，图8、图9、图10为三个目标通过20随机频点数据得到的距离向恢复结果，与200频点数据得到的距离向恢复结果相吻合，同样无距离模糊现象出现。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (1)

1.一种基于稀疏频点的微波、毫米波、太赫兹三维全息成像的重建方法，适用于微波、毫米波、太赫兹波的雷达发射信号，其特征在于，包括：

S1：确定雷达发射信号进行采样的稀疏频点个数N，所述稀疏频点个数N满足：

其中，B为设定中间量，且B＝ln M-ln(1-β)，M为待恢复的频点个数，即满频点的个数；β为最大可忍受的峰值旁瓣比r的置信系数；B_f为雷达发射信号的带宽，c为电磁波在自由空间的传播速度，ΔR为在雷达天线有效波束角的照射范围内的重建目标表面在距离向最大分布范围长度；

S2：如果雷达发射信号属于线性调频连续波信号则先补偿掉该发射方式导致的剩余视频相位，消除剩余视频相位后得到和步进频连续波信号相同的回波信号后，再进行以下处理；如果雷达发射信号属于步进频连续波信号，则直接雷达发射信号进行以下处理：

在带宽内选取均匀随机分布的N个频点f₁，f₂，…，f_N，定义f_sparse＝(f₁，f₂，…，f_N)；

采样得到N个稀疏频点的回波数据S₁(x，y，k_sparse)，建立坐标系XYZ；其中定义X、Y为扫描方向，Z为雷达照射方向，(x，y，0)为在坐标系XYZ中的天线相位中心位置，k_sparse为f_sparse对应的自由空间的波数；

利用式(2)获得回波信号S₁(x，y，k_sparse)的转置

其中，

m＝1，2…N，σ(x′，y′，z′)为在坐标系XYZ中位于(x′，y′，z′)处目标的散射系数，a_θ(·)为雷达波束角的窗函数，

S3：对S2所得与设定的参考信号利用式(4)求相关：

其中， 为的共轭矩阵；

将Corr(R′_j)结果中的最大值max[Corr(R′_j)]对应的距离R′_j作为在雷达天线有效波束角的照射范围内的目标表面在距离向估计值，记为距离R′_max；

S4：对S2所得回波数据进行频谱搬移，利用式(5)结合S3获得的距离R′_max将其频谱搬移至基带，得到

式中о表示为矩阵的Hadamard积；

S5：对S4所得结果利用式(6)对波数进行线性插值得到 

式中k＝(k′₁，k′₂，…，k′_M)为满频点的各频率对应的自由空间的波数，其中k′_i＝2π[f_c-B_f/2+(i-1)B_f/(M-1)]/c，i＝1，2，…，M，f_c为雷达发射信号的中心频率，M为满频点的个数；interp_X[Y(X)]为插值函数；

S6：对S5所得结果进行频谱搬移，结合S3获得的距离R′_max利用式 (7)将其频谱的中心频率搬移至其实际所处的频率，得到

式中

S7：利用相位偏移算法对S6所得结果进行图像处理，得到恢复图像。