一种基于GPU的快速频域后向投影三维成像方法技术

本发明专利技术公开了一种基于GPU的快速频域后向投影三维成像方法，该方法基于GPU的并行化处理能力，运用三维频域后向投影算法，通过利用图形处理器GPU采用一种单程序多数据的指令模式来对其丰富的硬件资源进行调度的能力，使图形处理器GPU中的一个线程负责一个频谱网格点，在各个线程内并行计算当前频域网格点在所有阵列向的后向投影值并累加，最后通过三维傅里叶逆变换得到三维SAR图像，从而大幅度提升了后向投影算法的运行效率。与现有技术相比，本发明专利技术具有较高的运算效率实现了线阵SAR的三维成像的特点，适用于合成孔径雷达成像、地球遥感等领域。

【技术实现步骤摘要】
一种基于GPU的快速频域后向投影三维成像方法
：本技术专利技术属于雷达
，它特别涉及了合成孔径雷达(SAR)成像

技术介绍
：由于具有全天时、全天候和大场景观测等优势，合成孔径雷达(SAR)已成为当今大面积地形测绘的一项重要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。线阵SAR(LASAR)是传统二维SAR成像技术目标维数空间分辨能力的扩展，可以获得观测目标场景的三维雷达成像，能够更加精细地描述观测场景中目标的几何和散射特征，较传统二维SAR提高了雷达系统的目标特征提取和目标识别能力，成为了近年来SAR成像技术的热点研究课题(详见参考文献“张清娟,李道京,李烈辰.连续场景的稀疏阵列SAR侧视三维成像研究.电子与信息学报,2013,(5):1097-1102.”)。线阵SAR成像系统的基本原理是通过线阵天线的运动合成一个大的虚拟二维面阵天线，获得面阵平面内的二维高分辨，再结合脉冲压缩技术获得雷达视线方向高分辨率，从而实现对观测目标场景的三维成像。后向投影(BP)算法是一种精确的SAR时域成像算法，它首先将合成孔径雷达原始数据沿距离向进行距离压缩(脉冲压缩)，然后通过选择不同慢时间观测空间中任意像素点在距离压缩后SAR数据空间中的数据，补偿方位向多普勒相位，并进行相干积累，最终获得各像素点散射系数的成像算法。由于在精确已知天线相位中心(AntennaPhaseCenter,APC)的前提下，BP算法可以有效补偿运动误差，因而已被广泛应用，详见“师君.双基地SAR与线阵SAR原理及成像技术研究[D].电子科技大学博士论文.2009”。图形处理器(GraphicsProcessingUnit,GPU)是一种专为执行复杂的数学和几何运算而设计的微处理器，常用于图形运算，是显卡的核心部件。GPU非常适用于有大规模浮点运算的场景，有强大的并发计算能力，与中央处理器(CentralProcessingUnit,简称CPU)相比，计算效率有成百上千倍的提升，而在SAR成像过程中涉及大量的复数运算，因此非常适合用GPU来加速。频域BP成像算法将距离压缩回波的频谱投影到图像的波数域以得到图像的频谱，再将图像的频谱通过傅立叶逆变换到图像域得到三维SAR图像。频域BP算法不但具有时域BP算法成像精度高的优点，而且具有更高的成像效率。虽然频域BP算法与时域BP算法相比有较大的效率提升，但在场景较大的情况下仍需要耗费较多时间，因为随着场景的增大，成像时需要划分的网格点数呈指数级增加，而成像过程中需要逐个网格点计算后向投影值，计算量极为巨大，限制了其实际应用。但频域BP算法具有高度并行化的结构，算法中各个网格点后向投影值的计算是相互独立的，将其在GPU平台上进行并行加速后，计算效率有极大的提升。
技术实现思路
：本专利技术提供了一种基于GPU的快速频域后向投影三维成像方法，该方法基于GPU的并行化处理能力，运用三维频域后向投影算法，让一个线程负责一个频谱网格点，在各个线程内并行计算当前网格点在所有阵列向的后向投影值并累加，再经傅里叶逆变换得到三维SAR图像，从而极大地提升了后向投影算法的运行效率。为了方便描述本专利技术的内容，首先作以下术语定义：定义1、标准SINC函数插值方法插值利用采样后的离散信号重建连续时间信号，是采样的逆问题。采用公式进行SINC函数插值，其中，gd(i)为采样信号，g(x)为SINC函数插值后的信号，L为插值核的点数，sinc(·)为SINC函数，∑(·)表示求和运算符号，详细内容可参考文献“LanG.CummingFrank,HWong.合成孔径雷达成像:算法与实现[M].电子工业出版社,2012”。定义2、标准傅里叶变换和逆变换傅里叶变换是SAR处理中的一个重要工具，对于长度为N的离散信号g(n)，傅里叶变换可以写为傅里叶逆变换可以写为其中，g(n)为原始离散信号，N为原始离散信号g(n)的长度，G(k)为原始离散信号g(n)经过傅里叶变换后的信号，∑(·)表示求和运算符号，exp(·)表示e指数运算符号，j为虚数单位，π为圆周率。详细内容可参考文献“LanG.CummingFrank,HWong.合成孔径雷达成像:算法与实现[M].电子工业出版社,2012”。定义3、线阵合成孔径雷达(LASAR)线阵合成孔径雷达成像是将线性阵列天线固定于载荷运动平台上并与平台运动方向垂直，结合运动平台的运动以合成二维平面阵列以实现阵列平面二维成像，再利用雷达波束向回波延时实现距离一维成像，从而实现观测目标三维成像的一种合成孔径雷达技术，详见文献“王斌,王彦平,洪文,等.线阵SAR三维成像分辨率分析[J].计算机仿真,2011,28(3):282-286”。定义4、标准合成孔径雷达距离压缩方法标准合成孔径雷达距离压缩方法是指利用合成孔径雷达系统的发射信号参数，生成距离压缩参考信号，并采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程，详见文献“保铮,邢孟道,王彤.雷达成像技术[M].电子工业出版社,2005”。定义5、线阵合成孔径雷达后向投影成像算法线阵合成孔径雷达的后向投影成像算法，简称BP成像算法。BP成像算法首先利用雷达平台的轨迹信息求出雷达平台的与场景像素点的距离历史，然后通过距离历史找出回波数据中对应的复数数据，对回波数据进行相位补偿后进行相干累加，从而得到该像素点的复图像值。详见“师君，双基地SAR与线阵SAR原理及成像技术[D].电子科技大学博士论文.2009”。定义6、图形处理器(GraphicsProcessingUnit,简称GPU)图形处理器是一种专为执行复杂的数学和几何运算而设计的微处理器，常用于图形运算，是显卡的核心部件。它非常适用于有大规模浮点运算的场景，有强大的并发计算能力，与中央处理器(CentralProcessingUnit,简称CPU)相比，计算效率有成百上千倍的提升。定义7、Host端和Device端图形处理器GPU的应用程序需要CPU端和GPU端协作运行，由于GPU一般作为显卡设备安装在计算机中，所以常称GPU端为设备端(Device端)，CPU端为主机端(Host端)。定义8、图形处理器的变量blockDim、gridDim、threadIdx和blockIdx图形处理器GPU采用一种单程序多数据的指令模式来对其丰富的硬件资源进行调度，其并发执行的最小单元是thread，固定数量的thread组成一个block，一个图形处理器中又包含多个block。图形处理器产品提供了通用的并行计算架构CUDA，CUDA库中定义变量blockIdx来标识一个block的索引，内置变量gridDim表示一个GPU中总的block数目；定义变量threadIdx标识一个block中线程的索引，内置变量blockDim表示一个block中总的线程数目，其值由用户设定，一般取2的整数次幂。上述变量均有三个维度，GPU中的所有block和一个block中的线程均组织成三维形式。定义9、线阵合成孔径雷达的慢时刻与快时刻线阵合成孔径雷达运动平台飞过一个方位向合成孔径长度所需要的时间称为慢时间，雷达系统以一定时间长度的重复周期发射接收脉冲，因此慢时间可以表示本文档来自技高网...

【技术保护点】
一种基于GPU的快速频域后向投影三维成像方法，其特征是它包括以下步骤：步骤1、初始化线阵SAR系统参数和回波数据：初始化线阵SAR系统参数包括：平台速度矢量，记做V；线阵天线各阵元初始位置矢量，记做Pn(0)，其中n为第n个线阵天线的阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数；线阵天线长度，记做L；相邻天线阵元之间的间距，记做d；雷达平台高度，记做R0；雷达工作中心频率，记做fc；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做Br；雷达发射信号脉冲宽度，记做TP；雷达发射信号调频斜率，记做fdr；雷达接收系统的采样频率，记做fs；雷达发射系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记做PRI；天线在方位向的有效孔径长度，记做Da；光在空气中的传播速度，记做C；距离向快时刻，记做t，t＝1,2,…,T，T为距离向快时刻总数；方位向慢时刻，记做l，l＝1,2,…,I，I为方位向慢时刻总数；上述参数均为线阵SAR系统标准参数，其中线阵天线的阵元总数N，线阵天线长度L，相邻天线阵元之间的间距d，雷达平台高度R0，雷达中心频率fc，雷达载频波长λ，雷达发射基带信号的信号带宽Br，雷达发射信号脉冲宽度TP，雷达发射信号调频斜率fdr，雷达接收系统的采样频率fs，雷达系统的脉冲重复频率PRF，雷达系统的脉冲重复时间PRI，天线在方位向的有效孔径长度Da在线阵SAR系统设计和观测过程中已经确定；平台速度矢量V及线阵天线各阵元初始位置矢量Pn(0)在线阵SAR观测方案设计中已经确定；根据线阵SAR成像系统方案和观测方案，线阵SAR成像方法需要的初始化系统参数均为已知；在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据，记做s(t,l,n)，t＝1,2,…,T，l＝1,2,…,I，n＝1,2,…,N，其中，T为步骤1中初始化得到的距离向快时刻总数，I为步骤1中初始化得到的方位向慢时刻总数，N为步骤1中初始化得到的线阵天线的阵元总数，在仿真线阵SAR成像过程中，原始回波数据s(t,l,n)可根据线阵SAR成像系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在线阵SAR数据进行成像之前，原始回波数据s(t,l,n)均已知；步骤2、初始化线阵SAR的观测场景目标空间参数：初始化线阵SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标作为线阵SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分成大小相等的立体单元格，单元格在水平横向、水平纵向和高度向的边长分别记做dx、dy和dz，单元格在水平横向、水平纵向和高度向的总数分别记做Mx、My和Mz；观测场景目标空间Ω中的单元格总数，记做M，观测场景目标空间Ω在水平横向、水平纵向和高度向的宽度分别记做Wx、Wy和Wz；观测场景目标空间Ω中第m个单元格的位置矢量，记做Pm，m表示观测场景目标空间Ω中第m个单元格，m为正整数，m＝1,2,…,M,M为观测场景目标空间Ω中的单元格总数；将观测场景目标空间Ω中所有单元格的目标散射系数按位置顺序排列组成M行1列的向量，记做散射系数向量α，散射系数向量α在线阵SAR三维成像观测仿真过程中就已经确定，散射系数向量α中的第m个元素，记做αm；观测场景目标空间Ω在线阵SAR成像方案设计中已经确定；步骤3、原始回波数据预处理：步骤3.1、计算阵元位置矩阵：采用公式K＝I·N计算得到阵列向慢时刻总数，记做K，其中，I为步骤1中初始化得到的方位向慢时刻总数，N为步骤1中初始化得到的线阵天线的阵元总数；定义阵列向慢时刻k，记做k＝1,2,…,K；根据步骤1中初始化得到的平台速度矢量V、线阵天线各阵元初始位置矢量Pn(0)和雷达发射系统的脉冲重复频率PRF，采用公式Pn(l)＝Pn(0)+V·l/PRF，n＝1,2,…,N，l＝1,2,…,I，计算得到第n个线阵天线阵元在第l个方位向慢时刻的位置矢量，记做Pn(l)，其中，N为步骤1中初始化得到的线阵天线阵元总数，I为步骤1中初始化得到的方位向慢时刻总数；将所有阵元的所有方位向的位置矢量Pn(l)按先阵元位置后方位向的顺序排列组成3行K列的矩阵，记做阵元位置矩阵P；步骤3.2、对线阵SAR原始回波数据进行距离压缩：将步骤1中初始化得到的原始回波数据s(t,l,n)按先阵元位置后方位向的顺序排列组成回波信号矩阵S，回波信号矩阵S由K行T列组成，其中，K为步骤3.1中计算得到的阵列向慢时刻总数，T为步骤1中初始化得到的距离向快时刻总数；采用标准合成孔径雷达距离压缩方法对回波信号矩阵S进行距离向脉冲压缩，得到距离压缩后的回波数据，记做E；步骤3.3、对距离压缩后的回波数据进行傅里叶变换：采用标准傅...

【技术特征摘要】
1.一种基于GPU的快速频域后向投影三维成像方法，其特征是它包括以下步骤：步骤1、初始化线阵SAR系统参数和回波数据：初始化线阵SAR系统参数包括：平台速度矢量，记做V；线阵天线各阵元初始位置矢量，记做Pn(0)，其中n为第n个线阵天线的阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数；线阵天线长度，记做L；相邻天线阵元之间的间距，记做d；雷达平台高度，记做R0；雷达工作中心频率，记做fc；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做Br；雷达发射信号脉冲宽度，记做TP；雷达发射信号调频斜率，记做fdr；雷达接收系统的采样频率，记做fs；雷达发射系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记做PRI；天线在方位向的有效孔径长度，记做Da；光在空气中的传播速度，记做C；距离向快时刻，记做t，t＝1,2,…,T，T为距离向快时刻总数；方位向慢时刻，记做l，l＝1,2,…,I，I为方位向慢时刻总数；上述参数均为线阵SAR系统标准参数，其中线阵天线的阵元总数N，线阵天线长度L，相邻天线阵元之间的间距d，雷达平台高度R0，雷达中心频率fc，雷达载频波长λ，雷达发射基带信号的信号带宽Br，雷达发射信号脉冲宽度TP，雷达发射信号调频斜率fdr，雷达接收系统的采样频率fs，雷达系统的脉冲重复频率PRF，雷达系统的脉冲重复时间PRI，天线在方位向的有效孔径长度Da在线阵SAR系统设计和观测过程中已经确定；平台速度矢量V及线阵天线各阵元初始位置矢量Pn(0)在线阵SAR观测方案设计中已经确定；根据线阵SAR成像系统方案和观测方案，线阵SAR成像方法需要的初始化系统参数均为已知；在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据，记做s(t,l,n)，t＝1,2,…,T，l＝1,2,…,I，n＝1,2,…,N，其中，T为步骤1中初始化得到的距离向快时刻总数，I为步骤1中初始化得到的方位向慢时刻总数，N为步骤1中初始化得到的线阵天线的阵元总数，在仿真线阵SAR成像过程中，原始回波数据s(t,l,n)可根据线阵SAR成像系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在线阵SAR数据进行成像之前，原始回波数据s(t,l,n)均已知；步骤2、初始化线阵SAR的观测场景目标空间参数：初始化线阵SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标作为线阵SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分成大小相等的立体单元格，单元格在水平横向、水平纵向和高度向的边长分别记做dx、dy和dz，单元格在水平横向、水平纵向和高度向的总数分别记做Mx、My和Mz；观测场景目标空间Ω中的单元格总数，记做M，观测场景目标空间Ω在水平横向、水平纵向和高度向的宽度分别记做Wx、Wy和Wz；观测场景目标空间Ω中第m个单元格的位置矢量，记做Pm，m表示观测场景目标空间Ω中第m个单元格，m为正整数，m＝1,2,…,M,M为观测场景目标空间Ω中的单元格总数；将观测场景目标空间Ω中所有单元格的目标散射系数按位置顺序排列组成M行1列的向量，记做散射系数向量α，散射系数向量α在线阵SAR三维成像观测仿真过程中就已经确定，散射系数向量α中的第m个元素，记做αm；观测场景目标空间Ω在线阵SAR成像方案设计中已经确定；步骤3、原始回波数据预处理：步骤3.1、计算阵元位置矩阵：采用公式K＝I·N计算得到阵列向慢时刻总数，记做K，其中，I为步骤1中初始化得到的方位向慢时刻总数，N为步骤1中初始化得到的线阵天线的阵元总数；定义阵列向慢时刻k，记做k＝1,2,…,K；根据步骤1中初始化得到的平台速度矢量V、线阵天线各阵元初始位置矢量Pn(0)和雷达发射系统的脉冲重复频率PRF，采用公式Pn(l)＝Pn(0)+V·l/PRF，n＝1,2,…,N，l＝1,2,…,I，计算得到第n个线阵天线阵元在第l个方位向慢时刻的位置矢量，记做Pn(l)，其中，N为步骤1中初始化得到的线阵天线阵元总数，I为步骤1中初始化得到的方位向慢时刻总数；将所有阵元的所有方位向的位置矢量Pn(l)按先阵元位置后方位向的顺序排列组成3行K列的矩阵，记做阵元位置矩阵P；步骤3.2、对线阵SAR原始回波数据进行距离压缩：将步骤1中初始化得到的原始回波数据s(t,l,n)按先阵元位置后方位向的顺序排列组成回波信号矩阵S，回波信号矩阵S由K行T列组成，其中，K为步骤3.1中计算得到的阵列向慢时刻总数，T为步骤1中初始化得到的距离向快时刻总数；采用标准合成孔径雷达距离压缩方法对回波信号矩阵S进行距离向脉冲压缩，得到距离压缩后的回波数据，记做E；步骤3.3、对距离压缩后的回波数据进行傅里叶变换：采用标准傅里叶变换对距离压缩后的回波数据E进行距离向傅里叶变换，得到距离向傅里叶变换后的回波数据，记做步骤4、利用SINC函...

【专利技术属性】
技术研发人员：张晓玲，蒲黎明，师君，韦顺军，余鹏，
申请(专利权)人：电子科技大学，
类型：发明
国别省市：四川,51