CN108802726A - 基于图形处理器gpu的合成孔径雷达成像方法 - Google Patents

基于图形处理器gpu的合成孔径雷达成像方法 Download PDF

Info

Publication number
CN108802726A
CN108802726A CN201810315379.3A CN201810315379A CN108802726A CN 108802726 A CN108802726 A CN 108802726A CN 201810315379 A CN201810315379 A CN 201810315379A CN 108802726 A CN108802726 A CN 108802726A
Authority
CN
China
Prior art keywords
pulse
orientation
frequency
data
synthetic aperture
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
CN201810315379.3A
Other languages
English (en)
Other versions
CN108802726B (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.)
Xidian University
Original Assignee
Xidian University
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 Xidian University filed Critical Xidian University
Publication of CN108802726A publication Critical patent/CN108802726A/zh
Application granted granted Critical
Publication of CN108802726B publication Critical patent/CN108802726B/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
    • 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)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于图形处理器GPU的合成孔径雷达成像方法，克服了现有技术中脉冲压缩和徙动校正调用cuFFT库函数进行傅里叶变换需要提前创建傅里叶变换计划、占用辅助显存和多次读写显存的缺点，实现步骤为：(1)传输合成孔径雷达基带数据；(2)进行距离脉冲压缩；(3)进行距离徙动校正；(4)进行方位脉冲压缩；(5)进行成像。本发明减少了合成孔径雷达成像算法占用的辅助显存，通过将中间结果存放在共享内存中，避免多次读写显存，提高合成孔径雷达成像算法速度。

Description

基于图形处理器GPU的合成孔径雷达成像方法
技术领域
本发明属于通信技术领域，更进一步涉及合成孔径雷达成像技术领域中的一种基于图形处理器GPU(Graphic Processing Unit)的距离多普勒成像方法。本发明可用于合成孔径雷达大规模回波数据的高速实时成像。
背景技术
实时处理基带数据是合成孔径雷达SAR(Synthetic Aperture Radar)信号处理系统的主要需求。随着合成孔径雷达成像技术向超高分辨率方向发展，对硬件处理能力提出了更高的要求。图形处理器GPU在浮点计算能力和存储带宽上相对传统的DSP(DigitalSignal Processor,DSP)有明显优势，容易进行并行程序实现，更适合对合成孔径雷达的大规模回波数据进行实时成像。
电子科技大学在其申请的专利文献“一种基于GPU的机载前视扫描雷达并行处理方法”(公开号：CN105137402A，申请号：2015105111842，申请日：2015年8月19日)中公开了一种基于GPU的机载前视扫描雷达并行处理方法。该方法的具体步骤是，首先利用NVIDIA公司CUDA运算平台中的cuFFT库来创建快速傅里叶变换(Fast Fourier Transformation,FFT)计划，并且综合利用CPU和GPU完成解卷积过程中的雷达天线方向图预处理；接着利用CPU读取雷达系统的原始数据并对其进行预处理并将预处理后的数据复制到GPU中；然后在GPU中分别并行完成距离向脉冲压缩和距离走动校正处理，并用GPU并行实现基于最大似然准则的解卷积成像算法；最后由CPU将最终的雷达成像结果保存。该方法存在的不足之处是，使用cuFFT库函数计算前必须先创建快速傅里叶变换计划，有较大时间消耗，cuFFT库函数计算时要占用与数据量大小相同的辅助显存，限制了算法能处理的回波数据大小和系统的实时处理能力。
西安电子科技大学在其申请的专利文献“基于GPU的弹载SAR前斜视成像方法”(公开号：CN104459693A，申请号：2014107194879，申请日：2014年12月1日)中公开了一种基于GPU的弹载SAR前斜视成像方法。该方法的步骤是，首先在CPU端设置合成孔径雷达的工作参数，将合成孔径雷达的工作参数存储于GPU中，将SAR原始回波数据以矩阵的形式存储于GPU中；在GPU端，对SAR原始回波数据矩阵作距离向处理，得出距离方位二维时域数据；在GPU端，将距离方位二维时域数据矩阵沿方位向进行补零扩展，得出大小为Na′×Nr的方位向补零扩展后的距离方位二维时域数据矩阵；在GPU端，对方位向补零扩展后的距离方位二维时域数据矩阵依次进行方位向处理，得出聚焦后的SAR图像数据。该方法存在的不足之处是，在SAR成像算法的距离向和方位向脉冲压缩过程中，使用cuFFT库函数需要将数据在图形处理器GPU和显存之间进行多次交换，脉冲压缩时间较长。
发明内容
本发明的目的在于针对上述已有技术的不足，提出一种基于图形处理器GPU的合成孔径雷达成像方法。
实现本发明目的的思路是，通过编写快速傅里叶变换算法，将脉冲压缩过程的中间计算结果存入共享内存中，实现高速的数据存取。本发明将矩阵转置操作在数据显存空间内进行，不占用辅助显存，可以处理更大规模的合成孔径雷达回波数据。
实现本发明目的具体步骤如下：
(1)传输合成孔径雷达基带数据：
(1a)选取与合成孔径雷达工作模式对应的方位采样点数和距离采样点数；
(1b)将实时采集的合成孔径雷达基带数据传输到图形处理器GPU的显存中；
(2)进行距离脉冲压缩：
(2a)在图形处理器GPU中分配与方位采样点数相同的线程块，每个线程块对应一个距离脉冲；
(2b)每个线程块从显存读取合成孔径雷达基带数据中对应的距离脉冲，利用斯托克汉姆Stockham快速傅里叶变换，将所读取方位的采样点从时域变换到距离频域；
(2c)对距离频域脉冲进行距离匹配滤波，将匹配滤波后的距离频域脉冲拷贝到与其对应的线程块的共享内存中；
(2d)每个线程块从共享内存读取匹配滤波后的距离频域脉冲，利用斯托克汉姆Stockham快速傅里叶逆变换，将所读取的距离频域脉冲变换到时域，得到压缩后的距离脉冲；
(2e)将压缩后的距离脉冲，拷贝到显存中；
(3)进行距离徙动校正：
(3a)对压缩后的距离脉冲，做二维傅里叶变换，得到二维频域脉冲；
(3b)对二维频域脉冲做距离徙动校正；
(3c)对校正后的二维频域脉冲，做二维傅里叶逆变换，得到校正后的距离脉冲；
(3d)将校正后的距离脉冲，拷贝到显存中；
(4)进行方位脉冲压缩：
(4a)在图形处理器GPU中分配与距离采样点数相同的线程块，每个线程块对应一个方位脉冲；
(4b)每个线程块从显存读取对应的方位脉冲，利用斯托克汉姆Stockham快速傅里叶变换，将方位脉冲从时域变换到方位频域；
(4c)对方位频域脉冲进行方位匹配滤波，将匹配滤波后的方位频域脉冲拷贝到与其对应的线程块的共享内存中；
(4d)每个线程块从共享内存读取匹配滤波后的方位频域脉冲，利用斯托克汉姆Stockham快速傅里叶逆变换，将方位频域脉冲变换到时域，得到压缩后的方位脉冲；
(4e)将压缩后的方位脉冲，拷贝到显存中；
(5)进行成像：
(5a)对压缩后的方位脉冲，进行方位多视操作；
(5b)对多视操作后的方位脉冲，进行灰度转换操作。
与现有技术相比，本发明具有以下优点：
第一，本发明在进行距离和方位脉冲压缩操作时，利用斯托克汉姆Stockham快速傅里叶变换进行频域匹配滤波，将中间结果序列放在共享内存中，克服了现有技术调用cuFFT库函数进行傅里叶变换需要提前创建傅里叶变换计划和多次读写显存的缺点，使得本发明的合成孔径雷达成像速度更快。
第二，本发明在进行距离徙动校正操作时，采用原址二维傅里叶变换操作，克服了现有技术中采用传统矩阵转置方法和调用cuFFT库函数进行时傅里叶变换时需要辅助显存的缺点，使得本发明在同等硬件条件下能处理更大规模的合成孔径雷达基带数据。
附图说明
图1是本发明的流程图；
图2是本发明显存中存放合成孔径雷达基带数据的示意图；
图3是本发明脉冲数据矩阵划分的示意图；
图4是本发明脉冲数据方阵子单元组划分的示意图。
具体实施方法
下面结合附图对本发明做进一步的描述。
本发明采用统一计算设备架构CUDA(Compute Unified Device Architecture)语言，可以在NVIDIA公司支持统一计算设备架构CUDA的图形处理器GPU上实现。
参照附图1，对本发明的具体步骤做进一步的描述。
步骤1，传输合成孔径雷达基带数据。
选取与合成孔径雷达工作模式对应的方位向采样点数Na和距离向采样点数Nr，Na与Nr的大小一般是2的次幂。通过cudaMemcpyAsync()接口函数，将实时采集到的基带数据从中央处理器CPU的锁页内存拷贝到图形处理器GPU的显存。基带数据存放形式如图2，图2中的横向箭头表示每个方位的采样点线性存放在显存中，纵向箭头表示各方位的采样数据按顺序存放在显存中，一个方格表示一个采样点。
步骤2，进行快速傅里叶变换。
对于Na个方位，每个方位Nr个采样点的基带数据，距离向快速傅里叶变换(FastFourier Transform,FFT)需要进行Na次长度为Nr点的计算，方位向FFT需要进行Na次长度为Nr点的计算，每次的计算由一个线程块完成。
本发明在CUDA环境下编写并行化的Stockham FFT计算函数。Stockham FFT算法与一般FFT算法相比，有数据顺序输入，顺序输出的特点，不需要对输入序列进行按位倒序操作，存储器的访问冲突较少。各级计算的数据读写是向量化的，更适合在GPU上用多核心同时执行。针对NVIDIA公司图形处理器GPU的硬件架构，本发明构造以Stockham基-4FFT蝶形操作作为基础计算单元的计算框架。对脉冲长度开4次方，得到基-4蝶形序列的总层数t。每个线程块内的线程，用t层Stockham基-4FFT蝶形操作完成序列长度为4t的FFT计算，其中每个线程块中的线程负责多个基-4蝶形操作的数据读写和计算，各层数据的读写在共享内存中进行。
Stockham基-4FFT蝶形单元的计算结果可由下式得到：
Y(r*i+k)＝(d1+d3)+(d2+d4)
其中，Y表示当前基-4蝶形序列，r*i+k、分别表示一个序列元素在当前基-4蝶形序列中的对应位置，r表示当前基-4蝶形序列的数据组总数，r＝4t-q，t表示基-4蝶形序列的总层数，q表示当前基-4蝶形序列所在的层数，q的取值范围为0到t-1，*表示相乘操作，i表示当前基-4蝶形操作的序号，i的取值范围为0到L表示当前基-4蝶形序列中一个数据组的数据总量，L＝4q，k表示线程块中的线程序号，d1、d2、d3、d4分别表示基-4蝶形操作的一个中间结果， X表示上一层基-4蝶形序列，4*r*i+k、4*r*i+r+k、4*r*i+2*r+k、4*r*i+3*r+k分别表示一个输入序列元素在上一层基-4蝶形序列中的对应位置，分别表示基-4蝶形操作的一个蝶形因子，j表示取复数的虚部符号。
对于4t点输入序列，在第1层，基-4蝶形序列一个数据组的数据总量L＝4，基-4蝶形序列的数据组总数r＝4t-1，每个数据组内有1个基-4蝶形操作单元，用r个线程从共享内存读入数据进行计算，线程号与数据组编号k对应，计算完成后，根据输出序列的索引，将数据放回共享内存；在第2层，基-4蝶形序列一个数据组的数据总量L＝42＝16，基-4蝶形序列的数据组总数r＝4t-2，每个数据组内有4个基-4蝶形操作单元，仍用r个线程从共享内存读入数据进行计算计算完成后，根据输出序列的索引，将数据放回共享内存；以此类推，直到最后一级蝶形计算，基-4蝶形序列一个数据组的数据总量L＝4t，共有1组数据，数据组内有4t-1个基-4蝶形操作单元，用r个线程从共享内存读入数据进行计算计算完成后，根据输出序列的索引，将数据放回共享内存。
对于t为奇数的2t点输入序列，可先用一次Cooley-Tokey基-2DIF-FFT操作，将序列分解为两个独立的4t'点序列，再使用所述的4t点Stockham FFT计算方法进行计算。
步骤3，进行距离匹配滤波。
在与步骤2相同的CUDA核函数内继续执行操作，每个线程块内的线程为其Stockham FFT结果计算对应的距离匹配滤波系数，距离匹配滤波系数公式如下：
其中，SRu表示第u个采样点的匹配滤波系数，u的取值范围为1到Z，Z表示采样数据长度，cos表示取余弦操作，fs表示采样频率，π表示圆周率，tp表示雷达脉冲宽度，B表示雷达发射信号带宽，j表示取复数的虚部符号，sin表示取正弦操作。
读取每个线程块共享内存中的距离频域脉冲，将距离频域脉冲与对应的距离匹配滤波系数相乘，得到匹配滤波后的距离频域脉冲。
步骤4，进行快速傅里叶逆变换。
本发明采用对傅里叶变换的输入输出取共轭的方法来实现傅里叶逆变换操作。在与步骤2、步骤3相同的CUDA核函数内继续执行操作，首先对距离匹配滤波数据进行取共轭操作。使用步骤2所述的StockhamFFT计算方法，进行Na次长度为Nr点的快速傅里叶变换。最后，对完成Stockham FFT计算的距离匹配滤波数据进行取共轭操作，得到基带数据的距离脉冲压缩结果。
步骤5，进行二维傅里叶变换。
本发明采用把两个维度的数据各进行一次傅里叶变换的方法，将采样数据变换到二维频域。为了最大化利用显存带宽，提高算法执行速度，需要保证输入数据连续地存放在显存中。所以在完成距离维数据的傅里叶变换后，先将数据矩阵进行转置，再对方位向数据进行傅里叶变换。传统的矩阵转置操作需要分配与数据量相等的存储空间，用于存放转置后的数据，针对SAR成像计算过程，本发明用一种原址矩阵转置方法，直接在数据的存储空间内进行矩阵转置。
第一步，利用步骤2所述的Stockham FFT计算方法，在距离向进行Na次长度为Nr点的快速傅里叶变换，将脉冲变换到距离频域。
第二步，对脉冲数据矩阵进行原址矩阵转置操作。
参照图3脉冲数据矩阵划分示意图，对脉冲数据矩阵的划分做进一步描述。图3中的每个方框代表一个脉冲数据方阵，斜线为方阵的对角线，箭头表示矩阵转置位置。将Na×Nr大小的数据矩阵划分为大小的数据方阵。
参照图4脉冲数据方阵子单元组划分示意图，对脉冲数据方阵子单元组的划分做进一步描述。图4中的方框表示一个脉冲数据方阵，每个带虚线的小方框表示一个32×32大小的方阵子单元，箭头表示数据转置位置，斜线表示方阵的对角线。把每个数据方阵划分为大小的方阵子单元。以方阵对角线为轴，把两个对称的方阵子单元划分为一个方阵子单元组。方阵子单元A和方阵子单元B是一个方阵子单元组，方阵子单元C就在方阵对角线上，自身作为一个子单元组。为每个子单元组，分配一个32线程的线程块，从显存读取子单元组内的数据，按行存放到共享内存，再按列从共享内存中读出，按行写入显存的对称位置，得到原址矩阵转置数据。需要注意的是，转置之后，原来连续的Na点数据现在分组存放在在个方阵中，后续的计算需要有跨度地读取转置后的数据，数据的读写在组内是连续的，所以能占满显存带宽。
第三步，利用步骤2所述的Stockham FFT计算方法，对原址矩阵转置后的脉冲数据矩阵进行Nr次长度为Na点的快速傅里叶变换，将脉冲从距离频域变换到二维频域，
第四步，将二维频域脉冲写入共享内存。
步骤6，进行距离徙动校正。
在与步骤5第四步相同的CUDA核函数内继续执行操作，每个线程块内的线程为其二维频域脉冲计算对应的距离徙动校正系数，距离徙动校正系数公式如下：
其中，QRu表示第u个采样点的徙动校正系数，Rs表示场景中心距离，fa表示回波的多普勒，faM表示回波的最大多普勒，fr表示距离频率；
读取每个线程块共享内存中的二维频域脉冲，将二维频域脉冲和对应的距离徙动校正系数相乘，得到校正后的二维频域脉冲。
步骤7，进行二维傅里叶逆变换。
本发明采用对二维傅里叶变换的输入输出序列取共轭的方法来实现二维傅里叶逆变换操作，将数据从二维频域变换到时域。首先对二维频域数据进行取共轭操作。再使用步骤5所述二维快速傅里叶变换方法计算输入序列。最后，对输出序列进行取共轭操作，得到二维傅里叶逆变换结果。
步骤8，进行快速傅里叶变换。
利用步骤2所述的StockhamFFT计算方法，在距离向进行Nr次长度为Na点的快速傅里叶变换，将采样数据变换到方位频域，结果写入共享内存。
步骤9，方位匹配滤波。
在与步骤8相同的CUDA核函数内继续执行操作，每个线程块内的线程为其StockhamFFT结果计算对应的方位匹配滤波系数，方位匹配滤波系数公式如下：
其中，SAu表示第u个采样点的匹配滤波系数，V表示雷达载机的速度，RB表示目标垂直距离；
读取共享内存中的方位频域脉冲，将方位频域脉冲与对应的方位匹配滤波系数相乘，得到匹配滤波后的方位频域脉冲。
步骤10，进行快速傅里叶逆变换。
本发明采用对傅里叶变换的输入输出取共轭的方法来实现傅里叶逆变换操作。首先对方位匹配滤波数据进行取共轭操作。再使用步骤2所述的StockhamFFT计算方法，进行Nr次长度为Na点的快速傅里叶变换。最后，对完成Stockham FFT计算的方位匹配滤波数据进行取共轭操作，得到基带数据的方位压缩结果。
步骤11，进行方位多视。
多视操作可以有效抑制图像的斑点噪声，提高辐射分辨率。大多数雷达系统中的固有方位分辨率高于距离分辨率，本发明选择视数为3的方位多视。将方位采样点数Na除以3，取整数部分，得到方位向成像点数。将方位向连续的3个方位压缩数据做相加操作，即可得到方位多视复图像数据。
分配Nr个线程，每个线程块分配32个线程，每个线程负责Na点的方位多视处理。各线程循环地读入3个方位压缩数据、累加、输出一个多视数据。多视操作完成后得到大小为复图像数据。其中，floor表示向下取整操作。
步骤12，进行灰度转换。
将方位多视复图像数据转换为人眼可以识别的灰度图像数据。
分配个线程块，每个线程块分配32个线程，每个线程块负责纵向数据块的处理。各线程按行读取方位多视复图像数据，求绝对值，放回显存，同时累加到私有求和变量，循环次完成整个图像的处理。
将Nr个累加得到的求和变量放到显存，用一个线程从显存中读取、累加，将累加结果除以图像大小得到实图像的平均值。
各线程循环读入实图像数据，除以平均值，再乘以灰度系数80进行归一化，将大于255的结果限制为255，得到8位灰度图像数据，写入显存。
本发明的效果可以通过以下仿真实验进一步说明。
1.仿真条件：
本发明的仿真条件如下表所示：
表1仿真条件表
操作系统 Windowsx64
编程环境 VS2013
GPU计算平台 CUDA7.5
GPU硬件 NVIDIAGTX970(Maxwell架构)
数据量 8192×4096单精度浮点复数
2.仿真内容：
选取与合成孔径雷达工作模式对应的方位向采样点数8192和距离向采样点数4096。通过cudaMemcpyAsync()接口函数，将实时采集到的合成孔径雷达基带数据从中央处理器CPU的锁页内存拷贝到图形处理器GPU的显存。
通过6层Stockham基-4FFT蝶形单元计算，可得到顺序的4096点FFT计算结果。使用8192个线程块计算8192个方位的基带数据，每个线程块设置256个线程，分配32768字节的共享内存用于Stockham基-4FFT计算的层间数据交换。每个线程分配16个浮点型复数寄存器，线程块内的线程同时计算4096点数据。第一层计算，每线程的连续四个寄存器的数据读取地址偏移按1024递增，在并行情况下256线程实际读取了4096点数据的4个1024数据段中的256个数据。所以之后的3组数据在读取数据时，索引要加上256，共4组数据。每次读取时，各线程是合并访问显存的，数据按256点向量化读取。每个线程计算4组基-4单元后，把数据放回。第二层计算读取时，基-4单元内部的数据读取以1024/4＝256为跨度，仍然有1024个基-4单元，256个线程中每个线程负责4个基-4单元的计算。第三层计算读取时，基-4单元内部的数据读取以256/4＝64为跨度。第四层计算读取时，基-4单元内部的数据读取以64/4＝16为跨度。第五层计算读取时，基-4单元内部的数据读取以16/4＝4为跨度。第六层计算读取时，基-4单元内部的数据读取以4/4＝1为跨度。第六层计算完成后，基带数据被转换到距离频域，将计算结果写入共享内存。
各线程计算与其距离频域数据对应的距离匹配滤波系数，从共享内存中读取FFT计算结果，将将距离频域数据与距离匹配滤波系数做相乘操作，得到距离匹配滤波结果，写入共享内存。
对距离匹配滤波数据进行取共轭操作。使用所述的4096点FFT计算方法，进行快速傅里叶变换。对完成FFT计算的距离匹配滤波数据进行取共轭操作，得到基带数据的距离脉冲压缩结果，写入显存。
二维快速傅里叶变换过程中使用的Stockham FFT计算函数与距离和方位脉冲压缩中的一致，除此之外还是使用了原址矩阵转置操作。
将8192×4096大小的数据矩阵分为2个4096×4096大小的数据方阵，4096×4096大小的数据方阵划分为128×128个32×32大小的方阵子单元，按方阵对角线轴对称的两个方阵子单元划分划分为一个子单元组。
在线程网络中分配128×128的二维线程块，每个线程块中设置32个线程，在共享内存中分配两块32×32大小的空间。32个线程连续读取方阵子单元数据，按行放入共享内存。按列将共享内存中的数据读出，按行写入当前方阵子单元在显存的对称子单元位置，实现原址矩阵转置。
在原址矩阵转置后进行方位向8192点FFT计算时，前后4096点输入序列在显存中的有4096×4096大小的地址偏移。再次进行原址矩阵转置操作后，数据存放方式和进行距离徙动校正前一致。
分配4096个线程块，每个线程块负责一个方位的8192点FFT计算，一个线程块内设置256个线程，与4096点FFT计算的区别是每个线程分配32个浮点复数寄存器，从而完整地存储8192点数据。将8192点输入序列先进行一次Cooley-Tokey基-2DIF-FFT，分解为2组4096点序列。再使用所述的4096点Stockham FFT算法进行计算。
各线程计算与其方位频域数据对应的方位匹配滤波系数，从共享内存中读取FFT计算结果，将将方位频域数据与方位匹配滤波系数做相乘操作，得到方位匹配滤波结果，写入共享内存。
对方位匹配滤波数据进行取共轭操作。使用所述的8192点FFT计算方法，进行快速傅里叶变换。对完成FFT计算的方位匹配滤波数据进行取共轭操作，得到基带数据的方位脉冲压缩结果，写入显存。
分配128个线程块，每个线程块32个线程，共4096个线程。
多视计算。将方位点数8192除3取整得到方位向成像点数2730，每个线程负责纵向8192点的方位多视处理。各线程循环地读入3个方位数据、累加、得到1个多视数据，写入显存，得到大小为4096×2730的复图像数据。
灰度转换。各线程负责对纵向2730点复数图像数据进行取绝对值操作，每次的结果写入显存，得到大小为4096×2730的实图像数据。同时，将每次的绝对值结果累加到线程私有变量SUM，最后得到4096个纵向累加数据，写入显存。用一个线程将4096个纵向累加数据进行求和，并处以图像大小4096×2730，得到图像平均值。各线程读入实图像数据，除以平均值，再乘以灰度系数80，最后把数值大小限制在0-255范围，强制转换为8bit整型数据，得到8位灰度图像数据，写入显存。
3.仿真结果分析：
本发明的仿真结果如表2、表3所示：
表2显存空间占用表
项目 传统GPU方法 本发明所述方法 占用减少
数据空间占用 256MB 256MB 0％
FFT计算空间占用 256MB 0 100％
矩阵转置空间占用 256MB 0 100％
其他暂存空间 0.78MB 0.78MB 0％
总计 768.78MB 256.78MB 66.6％
表3算法执行时间表
项目 传统GPU方法 本发明所述方法 加速比
创建cuFFT计划 345ms 0ms —
距离脉冲压缩 14.9ms 6.03ms 2.47
距离徙动校正 25.2ms 20.7ms 1.22
方位脉冲压缩 12.2ms 7.03ms 1.74
表2为显存空间占用表。与传统方法相比，本发明减少了传统方法调用cuFFT库函占用的辅助显存空间，减少了传统矩阵转置方法占用的转置空间。总计减少三分之二的显存占用，因此可以支持更大点数SAR成像算法。
表3为算法执行时间表。与传统方法相比，本发明减少了传统方法中调用cuFFT库函数创建cuFFT计划的时间。传统方法在调用cuFFT库函数之前，要把数据写入显存，后续的计算要等待库函数把数据写入显存。本发明所述方法通过实现Stockham FFT计算函数，将数据类型转化、距离向FFT、距离匹配滤波、距离向IFFT的计算过程整合到一起，中间结果通过共享内存实现快速存取，不需要通过显存存取数据，从而将计算时间由14.9ms减少到6.03ms，加速比达到2.47倍。除此之外距离徙动校正和方位脉冲压缩也分别有1.22倍和1.74倍的加速比。

Claims (10)

1.一种基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，利用图形处理器GPU的共享内存快速存取脉冲压缩操作中的数据，具体步骤包括如下：
(1)传输合成孔径雷达基带数据：
(1a)选取与合成孔径雷达工作模式对应的方位采样点数和距离采样点数；
(1b)将实时采集的合成孔径雷达基带数据传输到图形处理器GPU的显存中；
(2)进行距离脉冲压缩：
(2a)在图形处理器GPU中分配与方位采样点数相同的线程块，每个线程块对应一个距离脉冲；
(2b)每个线程块从显存读取合成孔径雷达基带数据中对应的距离脉冲，利用斯托克汉姆Stockham快速傅里叶变换，将所读取方位的采样点从时域变换到距离频域；
(2c)对距离频域脉冲进行距离匹配滤波，将匹配滤波后的距离频域脉冲拷贝到与其对应的线程块的共享内存中；
(2d)每个线程块从共享内存读取匹配滤波后的距离频域脉冲，利用斯托克汉姆Stockham快速傅里叶逆变换，将所读取的距离频域脉冲变换到时域，得到压缩后的距离脉冲；
(2e)将压缩后的距离脉冲，拷贝到显存中；
(3)进行距离徙动校正：
(3a)对压缩后的距离脉冲，做二维傅里叶变换，得到二维频域脉冲；
(3b)对二维频域脉冲做距离徙动校正；
(3c)对校正后的二维频域脉冲，做二维傅里叶逆变换，得到校正后的距离脉冲；
(3d)将校正后的距离脉冲，拷贝到显存中；
(4)进行方位脉冲压缩：
(4a)在图形处理器GPU中分配与距离采样点数相同的线程块，每个线程块对应一个方位脉冲；
(4b)每个线程块从显存读取对应的方位脉冲，利用斯托克汉姆Stockham快速傅里叶变换，将方位脉冲从时域变换到方位频域；
(4c)对方位频域脉冲进行方位匹配滤波，将匹配滤波后的方位频域脉冲拷贝到与其对应的线程块的共享内存中；
(4d)每个线程块从共享内存读取匹配滤波后的方位频域脉冲，利用斯托克汉姆Stockham快速傅里叶逆变换，将方位频域脉冲变换到时域，得到压缩后的方位脉冲；
(4e)将压缩后的方位脉冲，拷贝到显存中；
(5)进行成像：
(5a)对压缩后的方位脉冲，进行方位多视操作；
(5b)对多视操作后的方位脉冲，进行灰度转换操作。
2.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(2b)、步骤(4b)中所述斯托克汉姆Stockham快速傅里叶变换的步骤如下：
第一步，对每个线程块对应的脉冲长度开4次方，得到斯托克汉姆Stockham快速傅里叶变换的基-4蝶形序列的总层数；
第二步，利用下式，计算斯托克汉姆Stockham快速傅里叶变换每一层基-4蝶形序列的元素：
Y(r*i+k)＝(d1+d3)+(d2+d4)
其中，Y表示当前基-4蝶形序列，r*i+k、分别表示一个序列元素在当前基-4蝶形序列中的对应位置，r表示当前基-4蝶形序列的数据组总数，r＝4t-q，t表示基-4蝶形序列的总层数，q表示当前基-4蝶形序列所在的层数，q的取值范围为0到t-1，*表示相乘操作，i表示当前基-4蝶形操作的序号，i的取值范围为0到L表示当前基-4蝶形序列中一个数据组的数据总量，L＝4q，k表示线程块中的线程序号，d1、d2、d3、d4分别表示基-4蝶形操作的一个中间结果， X表示上一层基-4蝶形序列，4*r*i+k、4*r*i+r+k、4*r*i+2*r+k、4*r*i+3*r+k分别表示一个输入序列元素在上一层基-4蝶形序列中的对应位置，分别表示基-4蝶形操作的一个蝶形因子，j表示取复数的虚部符号。
3.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(2c)中所述距离匹配滤波的步骤如下：
第一步，利用下式，计算距离匹配滤波的系数：
其中，SRu表示第u个采样点的匹配滤波系数，u的取值范围为1到Z，Z表示脉冲长度，cos表示取余弦操作，fs表示采样频率，π表示圆周率，tp表示合成孔径雷达的脉冲宽度，B表示雷达发射信号带宽，sin表示取正弦操作；
第二步，读取每个线程块共享内存中的距离频域脉冲，将距离频域脉冲与对应的距离匹配滤波系数相乘，得到匹配滤波后的距离频域脉冲。
4.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(2d)、步骤(4d)中所述的斯托克汉姆Stockham快速傅里叶逆变换的步骤如下：
第一步，对输入序列进行取共轭操作；
第二步，对取共轭的输入序列，使用权利要求2所述的斯托克汉姆Stockham快速傅里叶变换，得到傅里叶逆变换的中间结果序列；
第三步，对傅里叶逆变换的中间结果序列进行取共轭操作，得到斯托克汉姆Stockham快速傅里叶逆变换的输出序列。
5.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(3a)中所述的二维傅里叶变换的步骤如下：
第一步，对脉冲数据矩阵，使用权利要求2所述的斯托克汉姆Stockham快速傅里叶变换，将一个维度的脉冲变换到频域；
第二步，将Na×Nr大小的脉冲数据矩阵划分为个Nr×Nr大小的脉冲数据方阵，Na表示方位采样点数，Nr表示距离采样点数；
第三步，将Nr×Nr大小的脉冲数据方阵划分为32×32大小的方阵子单元；
第四步，将按方阵对角线对称的两个方阵子单元划分为一个方阵子单元组；
第五步，在图形处理器GPU中为每个方阵子单元组分配一个线程块，将显存中的每一个方阵子单元组，按行写入对应线程块的共享内存；
第六步，按列从共享内存中读取方阵子单元的脉冲数据，再按行写入另一个子单元在显存中的对应位置，得到原址矩阵转置后的脉冲数据矩阵；
第七步，对原址矩阵转置后的脉冲数据矩阵，使用权利要求2所述的斯托克汉姆Stockham快速傅里叶变换，将另一个维度的脉冲变换到频域，得到二维频域脉冲。
6.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(3b)中所述距离徙动校正的步骤如下：
第一步，利用下式，计算距离徙动校正的系数：
其中，QRu表示第u个采样点的距离徙动校正系数，Rs表示场景中心距离，fa表示回波的多普勒频率，faM表示回波的最大多普勒频率，fr表示距离频率；
第二步，读取每个线程块共享内存中的二维频域脉冲，将二维频域脉冲和对应的距离徙动校正系数相乘，得到校正后的二维频域脉冲。
7.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(3c)中所述二维傅里叶逆变换的步骤如下：
第一步，对输入序列进行取共轭操作；
第二步，对取共轭的输入序列，使用权利要求5所述的二维傅里叶变换，得到二维傅里叶逆变换的中间结果序列；
第三步，对二维傅里叶逆变换的中间结果序列进行取共轭操作，得到二维傅里叶逆变换的输出序列。
8.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(4c)中所述方位匹配滤波的步骤如下：
第一步，利用下式，计算方位匹配滤波的系数：
其中，SAu表示第u个采样点的匹配滤波系数，V表示雷达载机的速度，RB表示目标垂直距离，表示开平方根操作；
第二步，读取共享内存中的方位频域脉冲，将方位频域脉冲与对应的方位匹配滤波系数相乘，得到匹配滤波后的方位频域脉冲。
9.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(5a)中所述的方位多视操作的步骤如下：
第一步，将雷达基带数据的方位采样点数除以3的结果取整数部分，得到方位向成像点数；
第二步，对连续3个方位的采样数据做相加操作，得到方位多视复图像数据。
10.根据权利要求1所述的基于图形处理器GPU的合成孔径雷达成像方法，其特征在于，步骤(5b)中所述的灰度转换操作的步骤如下：
第一步，对方位多视复图像数据做取绝对值操作，得到实图像数据；
第二步，对实图像数据做取平均值操作；
第三步，对实图像数据做除以平均值操作；
第四步，将实图像数据归一化到0到255之间，得到8位整形灰度图像数据。
CN201810315379.3A 2017-12-29 2018-04-10 基于图形处理器gpu的合成孔径雷达成像方法 Active CN108802726B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2017114776613 2017-12-29
CN201711477661 2017-12-29

Publications (2)

Publication Number Publication Date
CN108802726A true CN108802726A (zh) 2018-11-13
CN108802726B CN108802726B (zh) 2020-04-14

Family

ID=64094890

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810315379.3A Active CN108802726B (zh) 2017-12-29 2018-04-10 基于图形处理器gpu的合成孔径雷达成像方法

Country Status (1)

Country Link
CN (1) CN108802726B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109581362A (zh) * 2018-11-28 2019-04-05 中国科学院国家空间科学中心 合成孔径雷达高度计在可变脉冲簇模式下的信号处理方法
CN109828271A (zh) * 2018-12-29 2019-05-31 北京行易道科技有限公司 Sar成像方法及装置，sar成像系统
CN110109115A (zh) * 2019-05-09 2019-08-09 西安电子科技大学 基于fpga和ddr3的sar快速成像装置及方法
CN113805174A (zh) * 2021-09-13 2021-12-17 博微太赫兹信息科技有限公司 一种基于gpu的圆周合成孔径雷达图像重建方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101441271A (zh) * 2008-12-05 2009-05-27 航天恒星科技有限公司 基于gpu的sar实时成像处理设备
CN102854507A (zh) * 2012-09-12 2013-01-02 电子科技大学 一种基于gpu后向投影双站合成孔径雷达成像方法
CN104459666A (zh) * 2014-12-01 2015-03-25 西安电子科技大学 基于LabVIEW的弹载SAR回波仿真及成像方法
CN105974415A (zh) * 2016-06-24 2016-09-28 西安电子科技大学 一种机载sar方位空变运动误差的高精度补偿方法
JP2017106799A (ja) * 2015-12-09 2017-06-15 株式会社東芝 合成開口レーダ装置及びそのレーダ信号処理方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101441271A (zh) * 2008-12-05 2009-05-27 航天恒星科技有限公司 基于gpu的sar实时成像处理设备
CN102854507A (zh) * 2012-09-12 2013-01-02 电子科技大学 一种基于gpu后向投影双站合成孔径雷达成像方法
CN104459666A (zh) * 2014-12-01 2015-03-25 西安电子科技大学 基于LabVIEW的弹载SAR回波仿真及成像方法
JP2017106799A (ja) * 2015-12-09 2017-06-15 株式会社東芝 合成開口レーダ装置及びそのレーダ信号処理方法
CN105974415A (zh) * 2016-06-24 2016-09-28 西安电子科技大学 一种机载sar方位空变运动误差的高精度补偿方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
DU P 等: ""From CUDA to OpenCL:Towards a performance portable solution for multi-platform GPU programming"", 《PARALLEL COMPUTING》 *
俞惊雷 等: ""一种基于GPU的高效合成孔径雷达信号处理器"", 《信息与电子工程》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109581362A (zh) * 2018-11-28 2019-04-05 中国科学院国家空间科学中心 合成孔径雷达高度计在可变脉冲簇模式下的信号处理方法
CN109828271A (zh) * 2018-12-29 2019-05-31 北京行易道科技有限公司 Sar成像方法及装置，sar成像系统
CN110109115A (zh) * 2019-05-09 2019-08-09 西安电子科技大学 基于fpga和ddr3的sar快速成像装置及方法
CN110109115B (zh) * 2019-05-09 2022-12-02 西安电子科技大学 基于fpga和ddr3的sar快速成像装置及方法
CN113805174A (zh) * 2021-09-13 2021-12-17 博微太赫兹信息科技有限公司 一种基于gpu的圆周合成孔径雷达图像重建方法

Also Published As

Publication number Publication date
CN108802726B (zh) 2020-04-14

Similar Documents

Publication Publication Date Title
CN108802726A (zh) 基于图形处理器gpu的合成孔径雷达成像方法
KR102353241B1 (ko) 가속화된 수학 엔진
CN109949255B (zh) 图像重建方法及设备
CN107506828A (zh) 计算装置和方法
CN108898554A (zh) 提高图像分辨率的方法及相关产品
CN109992743A (zh) 矩阵乘法器
KR0130772B1 (ko) 고속디지탈신호처리프로세서
CN107239824A (zh) 用于实现稀疏卷积神经网络加速器的装置和方法
CN111461311B (zh) 基于众核处理器的卷积神经网络运算加速方法及装置
CN111383741B (zh) 医学成像模型的建立方法、装置、设备及存储介质
CN107451097B (zh) 国产申威26010众核处理器上多维fft的高性能实现方法
CN105866774A (zh) 线性调频信号极坐标格式成像算法的fpga实现方法
CN112703511B (zh) 运算加速器和数据处理方法
CN109117455A (zh) 计算装置及方法
WO2022151779A1 (zh) 卷积运算的实现方法、数据处理方法及装置
CN114461978A (zh) 数据处理方法、装置、电子设备及可读存储介质
US6728742B1 (en) Data storage patterns for fast fourier transforms
CN107085827A (zh) 基于硬件平台实现的超分辨力图像复原方法
CN106204669A (zh) 一种基于gpu平台的并行图像压缩感知方法
Chan et al. High-throughput 64k-point FFT processor for THz imaging radar system
CN106291551A (zh) 一种基于gpu的并行结构isar距离对准方法
Jin et al. GPU-based parallel implementation of SAR imaging
US20220035890A1 (en) Time Domain Unrolling Sparse Matrix Multiplication System and Method
CN113344765A (zh) 一种频域天文图像目标检测方法及系统
Liu et al. Efficient large-scale 1D FFT vectorization on multi-core vector accelerator

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