CN108957450A - 一种毫米波雷达gpu实时三维成像方法 - Google Patents

一种毫米波雷达gpu实时三维成像方法 Download PDF

Info

Publication number
CN108957450A
CN108957450A CN201810749676.9A CN201810749676A CN108957450A CN 108957450 A CN108957450 A CN 108957450A CN 201810749676 A CN201810749676 A CN 201810749676A CN 108957450 A CN108957450 A CN 108957450A
Authority
CN
China
Prior art keywords
wave
vertical
angle
dimensional
number domain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201810749676.9A
Other languages
English (en)
Other versions
CN108957450B (zh
Inventor
张秀锋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xi'an Hengfan Electronic Technology Co Ltd
Original Assignee
Xi'an Hengfan Electronic Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xi'an Hengfan Electronic Technology Co Ltd filed Critical Xi'an Hengfan Electronic Technology Co Ltd
Priority to CN201810749676.9A priority Critical patent/CN108957450B/zh
Publication of CN108957450A publication Critical patent/CN108957450A/zh
Application granted granted Critical
Publication of CN108957450B publication Critical patent/CN108957450B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9004—SAR image acquisition techniques
    • G01S13/9011—SAR image acquisition techniques with frequency domain processing of the SAR signals in azimuth
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/904—SAR modes

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明属于雷达成像领域，公开了一种毫米波雷达GPU实时三维成像方法，包括：在CPU端，获取雷达系统参数和三维场景回波数据；并对目标场景进行网格划分，获取三维目标场景的成像区域；在GPU端，从CPU端获取系统参数和回波数据，调用CUDA傅里叶变换接口函数实现回波数据在距离‑竖直‑角度三个维度的傅里叶变换；在GPU端设计线性插值内核函数，实现波数域从距离‑竖直‑角度到X‑Y‑Z采样均匀化；在GPU端，调用CUDA傅里叶变换接口函数实现插值后的数据在X‑Y‑Z三个维度的逆傅里叶变换，实现三维成像；本发明用于提高毫米波雷达三维实时成像效率，达到实时处理雷达回波数据的目的。

Description

一种毫米波雷达GPU实时三维成像方法
技术领域
本发明属于毫米波雷达成像技术领域，尤其涉及一种毫米波雷达GPU实时三维成像方法，具体为在CUDA架构中实现毫米波雷达三维成像的并行设计，用于提高成像效率，达到实时处理毫米波雷达三维成像的目的。
背景技术
结合旋转门合成孔径技术的毫米波雷达近场成像具有三维分辨率高特点，但数据维度高，算法需要方位、高度和距离纵深三维耦合处理，包含大量的矩阵运算和迭代运算，计算复杂度较高，采用传统串行计算方式进行三维成像聚焦处理，运行时间过长，难以达到实时成像需求；GPU(图形处理器)从专业图形处理具备多线程并行高性能浮点计算能力。CUDA是英伟达针对GPU平行计算专门设计并广泛应用编程模型，支持C/C++、FORTAN等多语言扩展，包含丰富函数库，比如 CUBLAS、CUFFT、CUDNN等。在CUDA架构中通过合理设计毫米波三维成像计算流程，可有效综合利用GPU大规模平行计算能力有效提高处理效率，达到毫米波雷达实时三维成像目的。
发明内容
针对上述现有技术的缺点，本发明的目的在于提供一种针对毫米波雷达GPU 实时三维成像方法，以解决毫米波雷达三维成像运算复杂度高、实时处理困难的问题。
本发明实现的技术思路为：在对旋转门工作模式下的毫米波雷达三维回波进行三维聚焦成像处理时，首先在CPU端按照串行的方式输入系统参数和三维场景回波数据，将系统参数和回波数据拷贝至GPU端；在GPU端进行角度维匹配滤波，主要包含矩阵FFT和IFT以及矩阵点乘操作；采用CUDA编程架构利用Stolt插值的方法在空间波数域进行非均匀采样向均匀采样的插值运算；在GPU端实现三维的IFT运算，重建直角坐标系下三维场景的散射系数分布。
为达到上述目的，本发明采用如下技术方案予以实现。
一种毫米波雷达GPU实时三维成像方法，包括如下步骤：
步骤1，在CPU端获取系统参数和三维场景回波信号，并对所述目标场景进行网格划分，将系统参数和三维场景回波数据拷贝至GPU端；
步骤2，在GPU端将所述三维场景回波数据变换到距离-竖直-角度三维波数域，进行角度维匹配滤波，并将匹配滤波后的三维波数域数据变换到距离-竖直二维波数域；
步骤3，将所述角度维匹配滤波之后的距离-竖直二维波数域的回波数据在空间波数域进行非均匀采样向均匀采样的Stolt插值运算，并进行三维逆傅里叶变换获得直角坐标系下三维场景的图像。
进一步地，在步骤1中，在CPU端获取系统参数和三维场景回波信号的方法如下：
在CPU端，采用串行方式读入系统参数和回波信号，回波信号具有如下形式：
其中，Kr＝2π(fc+fr)/c表示距离向波数，fc为雷达系统发射信号的载频，fr为线性调频信号的频率，z'表示雷达天线阵元的竖直(Z)方向坐标，θ表示雷达阵列的旋转角度；(x,y,z)表示场景目标散射点的坐标，σ(x,y,z)表示位于(x,y,z)散射点的散射系数，R表示雷达天线阵列的旋转半径，上式为球面波的表示形式，其平面波形式表示为：
其中，Kxy为Kr在XOY平面的分量，Kz'为Kr在竖直(Z)方向的分量，并存在关系式φ为积分式引入的积分变量，表示角度， Fσ(2Kxycosφ,2Kxysinφ,Kz')表示σ(x,y,z)在水平(X)-垂直(Y)-竖直(Z)三个维度的傅里叶变换结果，即在水平(X)-垂直(Y)-竖直(Z)三维的波数域表达。
进一步地，在步骤1中，对所述目标场景进行网格划分的方法如下：
对所述目标场景成像区域在空间波数域进行网格划分，表示如下：
Kx＝2π(-Xnum/2:Xnum/2-1)/(Xnum*Xres)
Ky＝2π(-Ynum/2:Ynum/2-1)/(Ynum*Yres)
Kz'＝2π(-Znum/2:Znum/2-1)/(Znum*Zres)
其中Kx为Kr在水平(X)方向的分量，Ky为Kr在垂直(Y)方向的分量，并满足关系Xnum、Ynum、Znum分别为X、Y、Z三个方向网格点数，Xres、Yres、 Zres分别为X、Y、Z三个方向的成像分辨率，(a:b)表示由a到b依次递增1的闭区间。
进一步地，步骤2具体方法为：
(2a)在GPU端进行初始化：分配显存空间，将从CPU端获取系统参数和三维场景回波数据，并将其复制到GPU显存上；
(2b)在GPU端调用cufft库函数，对所述目标场景回波信号在竖直(Z) 方向进行快速傅里叶变换，将回波信号变换到距离-竖直二维波数域，表达式如下：
改写为卷积的形式为：
S(Kxy,Kz',θ)＝g(θ-φ)*Fσ(2Kxycosθ,2Kxysinθ,Kz')
(2c)在GPU端，对距离-竖直二维波数域信号在角度频域进行角度维匹配滤波。
进一步地，步骤2c中，对距离-竖直二维波数域信号在角度频域进行角度维匹配滤波的方法具体为：
(2c1)在GPU端，调用cufft库函数，对所述距离-竖直二维波数域信号在角度(θ)方向进行快速傅里叶变换，将距离-竖直二维波数域信号变换到距离- 竖直-角度三维波数域，表达式如下：
S(Kxy,Kz',ξ)＝G(ξ,Kxy)Fσ(2Kxycosξ,2Kxysinξ,Kz')
其中ξ为θ的波数域形式，G(ξ,Kxy)为g(θ,Kxy)沿θ方向的傅里叶变换；
(2c2)在GPU端进行匹配滤波，设计复数乘法并行运算核函数，通过距离- 竖直-角度三维波数域信号与匹配滤波函数相乘完成角度维匹配滤波，匹配滤波信号为：
H(ξ,Kxy)＝G-1(ξ,Kxy)
角度维匹配滤波后的距离-竖直-角度三维波数域信号具有如下形式：
S(Kxy,Kz',ξ)＝Fσ(2Kxycosξ,2Kxysinξ,Kz')
(2c3)在GPU端调用cufft库函数，对所述角度维匹配滤波后的距离-竖直 -角度三维波数域信号在角度(θ)方向进行快速逆傅里叶变换，将角度维匹配滤波后的距离-竖直-角度三维波数域信号变换到距离-竖直二维波数域，表达式如下：
S(Kxy,Kz',θ)＝Fσ(2Kxycosθ,2Kxysinθ,Kz')。
进一步地，步骤3具体步骤为：
(3a)X、Y两个维度的波数域采样(2Kxycosθ,2Kxysinθ)是非均匀的，在变换到 X-Y-Z三维直角坐标系之前需要进行Stolt插值进行采样均匀化处理，在GPU端，采用步骤(1b)所述的(Kx,Ky)作为插值输入，(2Kxycosθ,2Kxysinθ)作为插值样本，对 Fσ(2Kxycosθ,2Kxysinθ,Kz')进行波数域采样均匀化插值，插值后的信号具有如下形式：
S(Kx,Ky,Kz')＝Fσ(Kx,Kx,Kz')
(3b)在GPU端，调用cufft库函数，对波数域采样均匀化后的信号S(Kx,Ky,Kz') 依次进行X、Y、Z三个方向的傅里叶变换(FFT)，如步骤(1a)所述，σ(x,y,z)与 Fσ(2Kxycosθ,2Kxysinθ,Kz')为三维傅里叶变换对，对角度滤波后的距离-竖直二维波数域信号在水平(X)-垂直(Y)-竖直(Z)进行三维逆傅里叶变换即可重建σ(x,y,z)。
进一步地，在步骤(2c)中，在GPU端实现距离-竖直二维波数域信号在角度维的匹配滤波，匹配滤波的实质是计算两个矩阵的哈达马积，通过GPU多线程并行技术，每个线程在一个时刻负责一个矩阵元素的乘法，通过同时调用多个线程来实现匹配滤波的快速化。
进一步地，在步骤(3a)中，在GPU端实现距离-竖直二维波数域到X-Y-Z 三个维度的均匀化插值，每个线程在一个时刻负责一个向量一个元素的插值运算，通过同时调用多个线程来提升插值效率，同时在寻找插值位置上，采用二分法降低插值算法的复杂度，提升插值算法效率。
本发明与现有技术相比具有如下优点：1)采用三维ωK算法实现毫米波雷达三维场景目标成像，算法运算量小，便于实时成像；2)基于GPU平台，并行实现毫米波雷达GPU实时三维成像处理，可以提高成像效率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的毫米波雷达GPU实时三维成像方法的并行过程示意图；
图2为本发明实施例提供的波数均匀化插值方式示意图；
图3为本发明实施例提供的基于CUDA编程的程序流程示意图；
图4为本发明实施例提供的矩阵转置内核设计方法示意图；
图5为本发明实施例提供的矩阵哈达马积内核设计方法示意图；
图6为本发明实施例提供的实测数据成像结果图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例提供一种毫米波雷达GPU实时三维成像方法，如图1所示，所述方法包括如下步骤：
步骤1，在CPU端，获取系统参数和三维场景回波信号，并对所述目标场景进行网格划分，将系统参数和回波数据拷贝至GPU端；
步骤1具体为：
(1a)在CPU端，采用串行方式读入系统参数和回波信号，回波信号具有如下形式：
其中，Kr＝2π(fc+fr)/c表示距离向波数，fc为雷达系统发射信号的载频，fr为线性调频信号的频率，z'表示雷达天线阵元的竖直(Z)方向坐标，θ表示雷达阵列的旋转角度；(x,y,z)表示场景目标散射点的坐标，σ(x,y,z)表示位于(x,y,z)散射点的散射系数，R表示雷达天线阵列的旋转半径。上式为球面波的表示形式，其平面波形式表示为：
其中，Kxy为Kr在XOY平面的分量，Kz'为Kr在竖直(Z)方向的分量，并存在关系式φ为积分式引入的积分变量，表示角度， Fσ(2Kxycosφ,2Kxysinφ,Kz')表示σ(x,y,z)在水平(X)-垂直(Y)-竖直(Z)三个维度的傅里叶变换结果，即在水平(X)-垂直(Y)-竖直(Z)三维的波数域表达。
(1b)对所述目标场景成像区域在空间波数域进行网格划分：
Kx＝2π(-Xnum/2:Xnum/2-1)/(Xnum*Xres)
Ky＝2π(-Ynum/2:Ynum/2-1)/(Ynum*Yres)
Kz'＝2π(-Znum/2:Znum/2-1)/(Znum*Zres)
其中Kx为Kr在水平(X)方向的分量，Ky为Kr在垂直(Y)方向的分量，并满足关系Xnum、Ynum、Znum分别为X、Y、Z三个方向网格点数，Xres、Yres、 Zres分别为X、Y、Z三个方向的成像分辨率，(a:b)表示由a到b依次递增1的闭区间。
步骤2，在GPU端，将所述三维场景回波数据变换到距离-竖直-角度三维波数域，进行角度维匹配滤波，并将匹配滤波后的三维波数域数据变换到距离-竖直二维波数域；
步骤2具体为：
(2a)在GPU端进行初始化：分配显存空间，将从所述CPU端获取所述雷达系统参数和所述回波信号，并将其复制到GPU显存上；
(2b)在GPU端，调用cufft库函数，对所述目标场景回波信号在竖直(Z) 方向进行快速傅里叶变换(FFT)，将回波信号变换到距离-竖直二维波数域，表达式如下：
改写为卷积的形式为：
S(Kxy,Kz',θ)＝g(θ-φ)*Fσ(2Kxycosθ,2Kxysinθ,Kz')
(2c)在GPU端，对距离-竖直二维波数域信号在角度频域进行角度维匹配滤波，具体包括：
(2c1)在GPU端，调用cufft库函数，对所述距离-竖直二维波数域信号在角度(θ)方向进行快速傅里叶变换(FFT)，将距离-竖直二维波数域信号变换到距离-竖直-角度三维波数域，表达式如下：
S(Kxy,Kz',ξ)＝G(ξ,Kxy)Fσ(2Kxycosξ,2Kxysinξ,Kz')
其中ξ为θ的波数域形式，G(ξ,Kxy)为g(θ,Kxy)沿θ方向的傅里叶变换；
(2c2)在GPU端进行匹配滤波，设计复数乘法并行运算核函数，通过距离- 竖直-角度三维波数域信号与匹配滤波函数相乘完成角度维匹配滤波，匹配滤波信号为：
H(ξ,Kxy)＝G-1(ξ,Kxy)
角度维匹配滤波后的距离-竖直-角度三维波数域信号具有如下形式：
S(Kxy,Kz',ξ)＝Fσ(2Kxycosξ,2Kxysinξ,Kz')
(2c3)在GPU端，调用cufft库函数，对所述角度维匹配滤波后的距离-竖直-角度三维波数域信号在角度(θ)方向进行快速逆傅里叶变换(IFT)，将角度维匹配滤波后的距离-竖直-角度三维波数域信号变换到距离-竖直二维波数域，表达式如下：
S(Kxy,Kz',θ)＝Fσ(2Kxycosθ,2Kxysinθ,Kz')
步骤3，在GPU端，将所述角度维匹配滤波之后的距离-竖直二维波数域回波数据在空间波数域进行非均匀采样向均匀采样的Stolt插值运算，并进行三维逆傅里叶变换获得直角坐标系下三维场景的图像。
步骤3具体为：
(3a)由步骤(2c3)所示，X、Y两个维度的波数域采样(2Kxycosθ,2Kxysinθ)是非均匀的，在变换到X-Y-Z三维直角坐标系之前需要进行Stolt插值进行采样均匀化处理，在GPU端，采用步骤(1b)所述的(Kx,Ky)作为插值输入，(2Kxycosθ,2Kxysinθ) 作为插值样本，对Fσ(2Kxycosθ,2Kxysinθ,Kz')进行波数域采样均匀化插值，插值过程如图2所示，插值后的信号具有如下形式：
S(Kx,Ky,Kz')＝Fσ(Kx,Kx,Kz')
(3b)在GPU端，调用cufft库函数，对波数域采样均匀化后的信号S(Kx,Ky,Kz') 依次进行X、Y、Z三个方向的傅里叶变换(FFT)，如步骤(1a)所述，σ(x,y,z)与 Fσ(2Kxycosθ,2Kxysinθ,Kz')为三维傅里叶变换对，对角度滤波后的距离-竖直二维波数域信号在水平(X)-垂直(Y)-竖直(Z)进行三维逆傅里叶变换即可重建σ(x,y,z)。
为提高成像效率，本发明利用GPU的并行特性，采用CUDA编程模型等技术手段实现三维ωK算法在GPU平台上并行实现。将成像算法并行化，首先需对三维ωK算法进行运算量分析，找出算法耗时步骤。经分析可知，该成像算法的运算量集中快速傅里叶变换和X-Y-Z三维波数域均匀化插值。
三维ωK算法包含大量的矩阵向量操作，尤其是对矩阵进行一维快速傅里叶变换(FFT)操作，而CUDA函数库提供了简单高效的常用函数，其中CUFFT库是一个用于多矩阵和向量进行一维和多维快速傅里叶变换操作的函数库，提供了丰富的快速傅里叶变换(FFT)类型接口，调用库函数时无需按照硬件特性设计复杂的内核函数就能获得很高的性能，因此本发明中以调用CUFFT库函数为主，设计编写CUDA内核函数为辅，进行毫米波雷达GPU实时三维成像的实现。
上述步骤2，3中，主要采用CUDA编程模型，结合CUFFT库函数实现了三维ωK算法(三维ωK算法的GPU实现流程图如图1所示)的并行编程，完成了毫米波雷达GPU实时三维成像的并行化。
CUDA编程模型的一般流程如图3所示，主体包含以下五个步骤：
1.GPU初始化，分配显存空间，
2.数据传输：将数据从主机内存复制到设备显存中
3.内核调用：调用CUDA提供的接口函数或自己设计的内核函数，实现算法并行处理
4.数据传输：将计算结果从设备拷贝回主机
5.释放显存空间，并销毁句柄指令
本发明中主要调用的CUFFT库函数有：cufftPlan1d(配置一维快速傅里叶变换句柄)，cufftExecC2C(执行浮点复数到浮点复数的傅里叶变换)， cufftDestory(销毁cufft句柄)。虽然使用CUFFT函数可以使计算性能达到及优，但是针对算法中表达式复杂的部分如插值部分则没有CUDA提供的接口，因此需对上述算法的部分步骤进行内核函数设计，内核函数设计的性能决定了并行化设计效率的提升值，内核函数设计如下：
矩阵转置并行实现：
由于CUDA中通常将数据组织成一维的形式，步骤2和步骤3中对回波信号进行竖直(Z)方向、角度(θ)方向傅里叶变换和傅里叶逆变换时，需要对矩阵进行转置操作。由于在矩阵转置操作过程中矩阵元素是相互独立的，而且执行的算法是相同的，因此这些操作通过设计并行核函数在GPU中实现，既可以发挥 GPU并行优势，同时可以避免数据在内存(CPU)和显存(GPU)之间的频繁传输。
本发明中，需要执行转置操作的矩阵规模不大，为了最大限度发挥GPU的多线程优势，可以在显存(GPU)上为转置结果分配与原矩阵相同的空间大小，通过CUDA核函数内部的线程索引来操作矩阵，如图4所示，在一次循环中，一个线程负责一个矩阵元素的位置转换，由于是多线程并行，因此在一次循环中就同时很多个矩阵元素的位置调换，一次循环结束时增加每个线程对应的元素位置偏移继续同时执行多个矩阵元素的位置调换，直到原矩阵所有元素位置调换完成，就完成了整个矩阵的转置操作。
哈达马积并行实现：
步骤(2c3)中距离-竖直-角度三维波数域信号的角度维匹配滤波主要通过计算矩阵间的哈达马积，即两个矩阵间对应元素的乘积，由于计算哈达马积时，矩阵元素之间互不干扰，因此非常适合采用GPU的多线程并行技术来提高运算效率，同时由于匹配滤波函数是确定的值，因此可以提前生成并存储起来，匹配滤波就简化为求矩阵间的哈达马积。
本发明中，对矩阵求哈达马积应尽可能使更多的线程进行并行计算，并行示意图如图5所示，通过CUDA核函数内部的线程索引来操作矩阵，在一次循环中，一个线程负责两个矩阵对应元素的乘积，由于是多线程并行，因此在一次循环中就同时计算了两个矩阵中很多对应元素的乘积，一次循环结束时增加每个线程对应的元素位置偏移继续同时执行矩阵多个元素的哈达马积，直到两个矩阵所有对应元素哈达马积计算完成，就完成了两个矩阵的哈达马积，即距离-竖直-角度三维波数域信号角度维匹配滤波。
线性插值并行实现：
步骤3中，对Fσ(2Kxycosθ,2Kxysinθ,Kz')进行波数域采样均匀化插值具体是通过线性插值的方式来实现，由于插值过程中每一个矩阵元素都要遍历(2Kxycosθ,2Kxysinθ) 寻找插值点的位置，采用CPU实现时程序循环层次过多，导致效率偏低，而采用 GPU实现，可以发挥GPU的多线程优势，使一个矩阵的多个元素同时执行插值算法，可以节约程序的运行时间。
从(2Kxycosθ,2Kxysinθ)到(Kx,Ky)的均匀化插值需要每一个待插值位置都遍历一次波数定义域，若采用普通的遍历方式，则需要在CUDA内核函数中频繁进行复杂的逻辑判断，会大大折损GPU的运算性能，由于是线性插值，为了提升效率，可以通过“二分法”在GPU上寻找插值点位置，相比CPU可以提升运算性能，线性插值的实现按照下式执行：
其中，S为插值的结果，X为输入数据在插值区间的位置Xmax和Xmin分别为插值区间位置的最大值和最小值，Smax和Smin分别为插值区间信号的最大值和最小值。
指令优化方式
由于整个三维ωK算法中需要频繁调用cufft库函数实现傅里叶变换和逆变换，而调用cufft库需要定义和销毁cufftHandle，由于初始化cufftHandle比较耗时，因此可以在算法执行前初始化cufftHandle，在每一次调用cufft库函数之后不要立即销毁cufftHandle，等到整个三维ωK成像算法执行完毕再销毁，可以大幅节约成像算法执行时间。
本发明的效果可以通过以下仿真实验说明：
1.仿真内容
为了验证本发明在CUDA架构下并行计算的优越性，通过一组仿真实验对本发明所需的时间与CPU串行所需的时间定量分析，仿真测试硬件平台参数如表1 所示，软件平台参数如表2所示：
表1硬件平台参数
CPU Intel(R)i7-7700K
内存 16GB
显卡 NVIDIA GeForce GTX1070
显卡显存 8GB
计算能力 6.1
表2软件平台参数
2.仿真结果及分析
为了测试本发明中GPU实时三维成像的效率和效果，通过一组实测数据通过本发明提到的方法进行三维实时成像，原始回波信号数据为双精度复数类型，数据个数为460*384*320，大小为862.5MB，由于本发明中成像算法主要的运算量集中与傅里叶变换和波数域插值两个部分，通过对比两个部分在CPU与GPU下的运算时间验证本发明所提方法的效率，通过成像结果验证本发明所提方法的正确性。
表3一维傅里叶变换，数据每组点数为460点，数据共384*320组，执行时间：
表4波数域均匀化插值，插值矩阵元素为双精度浮点型复数，插值前每组数据460个元素，插值后每组数据16个元素，数据共384*320组：
从表3、表4可看出，本发明采用的方法对毫米波实时三维成像速率有明显提升，成像算法中比较耗时的两个部分在GPU平台下实现，耗时均在毫秒量级，完全满足实时成像的要求。
1.通过调用cufft函数对算法进行并行实现，加快了三维ωK成像算法的处理速度；
综上，本发明能带来的有益效果是：
1.通过调用cufft函数对算法进行并行实现，加快了三维ωK成像算法的处理速度；
2.在内核函数设计上，通过优化指令流，使成像效率获得提升；
3.将耗时较长的波数域插值算法基于GPU平台并行化，实现毫米波雷达实时三维成像，进而缩短成像处理时间，提高成像效率。
本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括： ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (8)

1.一种毫米波雷达GPU实时三维成像方法，其特征在于，包括如下步骤：
步骤1，在CPU端获取系统参数和三维场景回波信号，并对所述目标场景进行网格划分，将系统参数和三维场景回波数据拷贝至GPU端；
步骤2，在GPU端将所述三维场景回波数据变换到距离-竖直-角度三维波数域，进行角度维匹配滤波，并将匹配滤波后的三维波数域数据变换到距离-竖直二维波数域；
步骤3，将所述角度维匹配滤波之后的距离-竖直二维波数域的回波数据在空间波数域进行非均匀采样向均匀采样的Stolt插值运算，并进行三维逆傅里叶变换获得直角坐标系下三维场景的图像。
2.根据权利要求1所述的方法，其特征在于，在步骤1中，在CPU端获取系统参数和三维场景回波信号的方法如下：
在CPU端，采用串行方式读入系统参数和回波信号，回波信号具有如下形式：
其中，Kr＝2π(fc+fr)/c表示距离向波数，fc为雷达系统发射信号的载频，fr为线性调频信号的频率，z'表示雷达天线阵元的竖直(Z)方向坐标，θ表示雷达阵列的旋转角度；(x,y,z)表示场景目标散射点的坐标，σ(x,y,z)表示位于(x,y,z)散射点的散射系数，R表示雷达天线阵列的旋转半径，上式为球面波的表示形式，其平面波形式表示为：
其中，Kxy为Kr在XOY平面的分量，Kz'为Kr在竖直(Z)方向的分量，并存在关系式φ为积分式引入的积分变量，表示角度，Fσ(2Kxycosφ,2Kxysinφ,Kz')表示σ(x,y,z)在水平(X)-垂直(Y)-竖直(Z)三个维度的傅里叶变换结果，即在水平(X)-垂直(Y)-竖直(Z)三维的波数域表达。
3.根据权利要求1所述的方法，其特征在于，在步骤1中，对所述目标场景进行网格划分的方法如下：
对所述目标场景成像区域在空间波数域进行网格划分，表示如下：
Kx＝2π(-Xnum/2:Xnum/2-1)/(Xnum*Xres)
Ky＝2π(-Ynum/2:Ynum/2-1)/(Ynum*Yres)
Kz'＝2π(-Znum/2:Znum/2-1)/(Znum*Zres)
其中Kx为Kr在水平(X)方向的分量，Ky为Kr在垂直(Y)方向的分量，并满足关系Xnum、Ynum、Znum分别为X、Y、Z三个方向网格点数，Xres、Yres、Zres分别为X、Y、Z三个方向的成像分辨率，(a:b)表示由a到b依次递增1的闭区间。
4.根据权利要求1所述的方法，其特征在于，步骤2具体方法为：
(2a)在GPU端进行初始化：分配显存空间，将从CPU端获取系统参数和三维场景回波数据，并将其复制到GPU显存上；
(2b)在GPU端调用cufft库函数，对所述目标场景回波信号在竖直(Z)方向进行快速傅里叶变换，将回波信号变换到距离-竖直二维波数域，表达式如下：
改写为卷积的形式为：
S(Kxy,Kz',θ)＝g(θ-φ)*Fσ(2Kxycosθ,2Kxysinθ,Kz')
(2c)在GPU端，对距离-竖直二维波数域信号在角度频域进行角度维匹配滤波。
5.根据权利要求4所述的方法，其特征在于，步骤2c中，对距离-竖直二维波数域信号在角度频域进行角度维匹配滤波的方法具体为：
(2c1)在GPU端，调用cufft库函数，对所述距离-竖直二维波数域信号在角度(θ)方向进行快速傅里叶变换，将距离-竖直二维波数域信号变换到距离-竖直-角度三维波数域，表达式如下：
S(Kxy,Kz',ξ)＝G(ξ,Kxy)Fσ(2Kxycosξ,2Kxysinξ,Kz')
其中ξ为θ的波数域形式，G(ξ,Kxy)为g(θ,Kxy)沿θ方向的傅里叶变换；
(2c2)在GPU端进行匹配滤波，设计复数乘法并行运算核函数，通过距离- 竖直-角度三维波数域信号与匹配滤波函数相乘完成角度维匹配滤波，匹配滤波信号为：
H(ξ,Kxy)＝G-1(ξ,Kxy)
角度维匹配滤波后的距离-竖直-角度三维波数域信号具有如下形式：
S(Kxy,Kz',ξ)＝Fσ(2Kxycosξ,2Kxysinξ,Kz')
(2c3)在GPU端调用cufft库函数，对所述角度维匹配滤波后的距离-竖直-角度三维波数域信号在角度(θ)方向进行快速逆傅里叶变换，将角度维匹配滤波后的距离-竖直-角度三维波数域信号变换到距离-竖直二维波数域，表达式如下：
S(Kxy,Kz',θ)＝Fσ(2Kxycosθ,2Kxysinθ,Kz')。
6.根据权利要求1所述的方法，其特征在于，步骤3具体步骤为：
(3a)X、Y两个维度的波数域采样(2Kxycosθ,2Kxysinθ)是非均匀的，在变换到X-Y-Z三维直角坐标系之前需要进行Stolt插值进行采样均匀化处理，在GPU端，采用步骤(1b)所述的(Kx,Ky)作为插值输入，(2Kxycosθ,2Kxysinθ)作为插值样本，对Fσ(2Kxycosθ,2Kxysinθ,Kz')进行波数域采样均匀化插值，插值后的信号具有如下形式：
S(Kx,Ky,Kz')＝Fσ(Kx,Kx,Kz')
(3b)在GPU端，调用cufft库函数，对波数域采样均匀化后的信号S(Kx,Ky,Kz')依次进行X、Y、Z三个方向的傅里叶变换(FFT)，如步骤(1a)所述，σ(x,y,z)与Fσ(2Kxycosθ,2Kxysinθ,Kz')为三维傅里叶变换对，对角度滤波后的距离-竖直二维波数域信号在水平(X)-垂直(Y)-竖直(Z)进行三维逆傅里叶变换即可重建σ(x,y,z)。
7.根据权利要求4所述的方法，其特征在于，在步骤(2c)中，在GPU端实现距离-竖直二维波数域信号在角度维的匹配滤波，匹配滤波的实质是计算两个矩阵的哈达马积，通过GPU多线程并行技术，每个线程在一个时刻负责一个矩阵元素的乘法，通过同时调用多个线程来实现匹配滤波的快速化。
8.根据权利要求4所述的方法，其特征在于，在步骤(3a)中，在GPU端实现距离-竖直二维波数域到X-Y-Z三个维度的均匀化插值，每个线程在一个时刻负责一个向量一个元素的插值运算，通过同时调用多个线程来提升插值效率，同时在寻找插值位置上，采用二分法降低插值算法的复杂度，提升插值算法效率。
CN201810749676.9A 2018-07-10 2018-07-10 一种毫米波雷达gpu实时三维成像方法 Active CN108957450B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810749676.9A CN108957450B (zh) 2018-07-10 2018-07-10 一种毫米波雷达gpu实时三维成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810749676.9A CN108957450B (zh) 2018-07-10 2018-07-10 一种毫米波雷达gpu实时三维成像方法

Publications (2)

Publication Number Publication Date
CN108957450A true CN108957450A (zh) 2018-12-07
CN108957450B CN108957450B (zh) 2022-12-09

Family

ID=64482465

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810749676.9A Active CN108957450B (zh) 2018-07-10 2018-07-10 一种毫米波雷达gpu实时三维成像方法

Country Status (1)

Country Link
CN (1) CN108957450B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110208753A (zh) * 2019-06-27 2019-09-06 电子科技大学 一种基于gpu的雷达目标回波信号获取方法
CN110632595A (zh) * 2019-09-24 2019-12-31 西南交通大学 一种主动毫米波成像方法及系统、储存介质、成像设备
CN110764088A (zh) * 2019-10-25 2020-02-07 哈尔滨工程大学 一种超分辨率驻点扫描实时成像算法
CN111289975A (zh) * 2020-01-21 2020-06-16 博微太赫兹信息科技有限公司 一种多gpu并行计算的快速成像处理系统
CN111650585A (zh) * 2020-08-04 2020-09-11 中国人民解放军国防科技大学 近场毫米波稀疏mimo扫描阵列全聚焦成像方法和装置
CN112261036A (zh) * 2020-10-20 2021-01-22 苏州矽典微智能科技有限公司 数据传出方法和装置
CN112540381A (zh) * 2020-11-17 2021-03-23 中国科学院西安光学精密机械研究所 非均匀快速傅里叶变换的非视域单进多出三维重建方法
CN113805174A (zh) * 2021-09-13 2021-12-17 博微太赫兹信息科技有限公司 一种基于gpu的圆周合成孔径雷达图像重建方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2120063A1 (en) * 2008-05-15 2009-11-18 The European Community, represented by the European Commission Radar-imaging of a scene in the far-field of a one-or two-dimensional radar array
CN103675810A (zh) * 2013-11-13 2014-03-26 中国科学院电子学研究所 穿墙雷达成像的方法
CN107064930A (zh) * 2017-03-29 2017-08-18 西安电子科技大学 基于gpu的雷达前视成像方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2120063A1 (en) * 2008-05-15 2009-11-18 The European Community, represented by the European Commission Radar-imaging of a scene in the far-field of a one-or two-dimensional radar array
CN103675810A (zh) * 2013-11-13 2014-03-26 中国科学院电子学研究所 穿墙雷达成像的方法
CN107064930A (zh) * 2017-03-29 2017-08-18 西安电子科技大学 基于gpu的雷达前视成像方法

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110208753A (zh) * 2019-06-27 2019-09-06 电子科技大学 一种基于gpu的雷达目标回波信号获取方法
CN110208753B (zh) * 2019-06-27 2023-04-25 电子科技大学 一种基于gpu的雷达目标回波信号获取方法
CN110632595A (zh) * 2019-09-24 2019-12-31 西南交通大学 一种主动毫米波成像方法及系统、储存介质、成像设备
CN110632595B (zh) * 2019-09-24 2022-11-01 西南交通大学 一种主动毫米波成像方法及系统、储存介质、成像设备
CN110764088A (zh) * 2019-10-25 2020-02-07 哈尔滨工程大学 一种超分辨率驻点扫描实时成像算法
CN110764088B (zh) * 2019-10-25 2023-10-27 哈尔滨工程大学 一种超分辨率驻点扫描实时成像算法
CN111289975B (zh) * 2020-01-21 2022-04-22 博微太赫兹信息科技有限公司 一种多gpu并行计算的快速成像处理系统
CN111289975A (zh) * 2020-01-21 2020-06-16 博微太赫兹信息科技有限公司 一种多gpu并行计算的快速成像处理系统
CN111650585A (zh) * 2020-08-04 2020-09-11 中国人民解放军国防科技大学 近场毫米波稀疏mimo扫描阵列全聚焦成像方法和装置
CN111650585B (zh) * 2020-08-04 2020-11-03 中国人民解放军国防科技大学 近场毫米波稀疏mimo扫描阵列全聚焦成像方法和装置
CN112261036A (zh) * 2020-10-20 2021-01-22 苏州矽典微智能科技有限公司 数据传出方法和装置
CN112261036B (zh) * 2020-10-20 2021-09-24 苏州矽典微智能科技有限公司 数据传出方法和装置
CN112540381B (zh) * 2020-11-17 2023-02-10 中国科学院西安光学精密机械研究所 非均匀快速傅里叶变换的非视域单进多出三维重建方法
CN112540381A (zh) * 2020-11-17 2021-03-23 中国科学院西安光学精密机械研究所 非均匀快速傅里叶变换的非视域单进多出三维重建方法
CN113805174A (zh) * 2021-09-13 2021-12-17 博微太赫兹信息科技有限公司 一种基于gpu的圆周合成孔径雷达图像重建方法

Also Published As

Publication number Publication date
CN108957450B (zh) 2022-12-09

Similar Documents

Publication Publication Date Title
CN108957450A (zh) 一种毫米波雷达gpu实时三维成像方法
CN104459666B (zh) 基于LabVIEW的弹载SAR回波仿真及成像方法
Shams et al. A survey of medical image registration on multicore and the GPU
Krone et al. Fast Visualization of Gaussian Density Surfaces for Molecular Dynamics and Particle System Trajectories.
CN103544682B (zh) 一种三维超声图像非局部均值滤波方法
JP3818369B2 (ja) 入力画像に最も類似の複数の３次元モデルから画像を選択する方法
CN104459693A (zh) 基于gpu的弹载sar前斜视成像方法
CN105549078B (zh) 不规则地震数据的五维插值处理方法及装置
CN102073982A (zh) 用gpu实现超大sar图像各向异性扩散滤波加速方法
CN109557540A (zh) 基于目标散射系数非负约束的全变差正则化关联成像方法
CN102663738A (zh) 三维图像配准方法及系统
CN107064930A (zh) 基于gpu的雷达前视成像方法
CN108919220A (zh) 基于嵌入式gpu的弹载sar前侧视成像方法
Khan et al. Implicit neural representations for medical imaging segmentation
US10564271B2 (en) Systems, methods and devices for highly-parallelized QUS-value determination for characterizing a specimen
Rahman et al. Parallel implementation of a spatio-temporal visual saliency model
CN106777536A (zh) 基于精细电磁仿真的电磁远场二、三维可视化处理方法
CN109871249A (zh) 一种远程桌面操作方法、装置、可读存储介质及终端设备
CN116563096B (zh) 用于图像配准的形变场的确定方法、装置以及电子设备
CN106291551B (zh) 一种基于gpu的并行结构isar距离对准方法
CN103927721A (zh) 基于gpu的运动目标边缘增强方法
CN110426688A (zh) 一种基于地形背景目标的sar回波模拟方法
Preethi et al. Gaussian filtering implementation and performance analysis on GPU
CN110176021A (zh) 结合亮度校正的显著性信息的水平集图像分割方法及系统
CN104483670B (zh) 基于gpu的合成孔径雷达回波仿真方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant