[发明专利]一种基于GPU的快速频域后向投影三维成像方法

权利要求书

1.一种基于GPU的快速频域后向投影三维成像方法，其特征是它包括以下步骤：

步骤1、初始化线阵SAR系统参数和回波数据：

初始化线阵SAR系统参数包括：平台速度矢量，记做V；线阵天线各阵元初始位置矢量，记做P_n(0)，其中n为第n个线阵天线的阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数；线阵天线长度，记做L；相邻天线阵元之间的间距，记做d；雷达平台高度，记做R₀；雷达工作中心频率，记做f_c；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做B_r；雷达发射信号脉冲宽度，记做T_P；雷达发射信号调频斜率，记做f_dr；雷达接收系统的采样频率，记做f_s；雷达发射系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记做PRI；天线在方位向的有效孔径长度，记做D_a；光在空气中的传播速度，记做C；距离向快时刻，记做t，t＝1,2,…,T，T为距离向快时刻总数；方位向慢时刻，记做l，l＝1,2,…,I，I为方位向慢时刻总数；上述参数均为线阵SAR系统标准参数，其中线阵天线的阵元总数N，线阵天线长度L，相邻天线阵元之间的间距d，雷达平台高度R₀，雷达中心频率f_c，雷达载频波长λ，雷达发射基带信号的信号带宽B_r，雷达发射信号脉冲宽度T_P，雷达发射信号调频斜率f_dr，雷达接收系统的采样频率f_s，雷达系统的脉冲重复频率PRF，雷达系统的脉冲重复时间PRI，天线在方位向的有效孔径长度D_a在线阵SAR系统设计和观测过程中已经确定；平台速度矢量V及线阵天线各阵元初始位置矢量P_n(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的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分成大小相等的立体单元格，单元格在水平横向、水平纵向和高度向的边长分别记做d_x、d_y和d_z，单元格在水平横向、水平纵向和高度向的总数分别记做M_x、M_y和M_z；观测场景目标空间Ω中的单元格总数，记做M，观测场景目标空间Ω在水平横向、水平纵向和高度向的宽度分别记做W_x、W_y和W_z；观测场景目标空间Ω中第m个单元格的位置矢量，记做P_m，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、线阵天线各阵元初始位置矢量P_n(0)和雷达发射系统的脉冲重复频率PRF，采用公式P_n(l)＝P_n(0)+V·l/PRF，n＝1,2,…,N，l＝1,2,…,I，计算得到第n个线阵天线阵元在第l个方位向慢时刻的位置矢量，记做P_n(l)，其中，N为步骤1中初始化得到的线阵天线阵元总数，I为步骤1中初始化得到的方位向慢时刻总数；将所有阵元的所有方位向的位置矢量P_n(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函数进行回波数据插值：

采用标准SINC函数插值方法对距离向傅里叶变换后的回波数据进行SINC插值，得到SINC插值后的回波数据，记做

步骤5、计算三维图像的频谱网格：

采用公式和计算得到在跨航向、方位向和距离向的空间频率间隔，分别记做△k_x、△k_y和△k_z，其中，W_x、W_y和W_z分别为步骤2中初始化得到的观测场景目标空间Ω在水平横向、水平纵向和高度向的宽度，π为圆周率；

采用公式和计算得到在跨航向、方位向和距离向的频谱网格点总数，分别记做N_x、N_y和N_z，其中，表示向上取整操作，d_x、d_y和d_z分别为步骤2中初始化得到的单元格在水平横向、水平纵向和高度向的边长，π为圆周率；

步骤6、计算图形处理器GPU中所需的block数量：

定义计算图形处理器GPU中的并行计算架构CUDA的内置变量blockDim在x、y和z方向上的值分别为BX、BY和BZ，采用公式和计算得到并行计算架构CUDA的内置变量gridDim在x、y和z方向上的值，分别记做gridDim.x、gridDim.y和gridDim.z，其中，N_x、N_y和N_z分别为步骤5中计算得到的在跨航向、方位向和距离向的频谱网格点总数；

步骤7、计算线程的索引：

采用公式idx＝threadIdx.x+blockIdx.x·BX、idy＝threadIdx.y+blockIdx.y·BY和idz＝threadIdx.z+blockIdx.z·BZ，计算得到线程的索引，分别记做idx、idy和idz，其中，threadIdx.x、threadIdx.y和threadIdx.z分别为并行计算架构CUDA的变量threadIdx在x、y和z方向上的值，其取值范围分别为0,1,…,BX、0,1,…,BY和0,1,…,BZ，BX、BY和BZ分别为步骤6中得到的并行计算架构CUDA的内置变量blockDim在x、y和z方向上的值，blockIdx.x、blockIdx.y和blockIdx.z分别为并行计算架构CUDA的变量blockIdx在x、y和z方向上的值，其取值范围分别为0,1,…,gridDim.x、0,1,…,gridDim.y和0,1,…,gridDim.z，gridDim.x、gridDim.y和gridDim.z分别为步骤6中计算得到的并行计算架构CUDA的内置变量gridDim在x、y和z方向上的值；

步骤8、信号频率离散化：

采用公式k_x＝(idx-1)·△k_x，k_y＝(idy-1)·△k_y和k_z＝2·π·f_c/C+(idz-1)·△k_z，计算得到信号频率在跨航向、方位向和距离向的波数域离散表示，分别记做k_x、k_y和k_z，其中，△k_x、△k_y和△k_z分别为步骤5中计算得到的在跨航向、方位向和距离向的空间频率间隔，idx、idy和idz分别为步骤7中计算得到的线程的索引，f_c为步骤1中初始化得到的雷达工作中心频率，C为步骤1中初始化得到的光在空气中的传播速度，π为圆周率；

步骤9、计算后向投影值的索引p：

采用公式k₁＝2π(f_c-f_s/2)/C，计算得到距离向起始频率，记做k₁，其中，f_c为步骤1中初始化得到的雷达工作中心频率，f_s为步骤1中初始化得到的雷达接收系统的采样频率，C为步骤1中初始化得到的光在空气中的传播速度，π为圆周率；

采用公式△k＝2π(f_s/T)/C，计算得到距离向频率间隔，记做△k，其中，f_s为步骤1中初始化得到的雷达接收系统的采样频率，T为步骤1中初始化得到的距离向快时刻总数，C为步骤1中初始化得到的光在空气中的传播速度，π为圆周率；

采用公式计算得到后向投影值的索引，记做p，其中，k_x、k_y和k_z分别为步骤8中计算得到的信号频率在跨航向、方位向和距离向的波数域离散表示；

步骤10、为后向投影值补偿相位因子并累加：

采用公式计算得到在跨航向、方位向和距离向需要补偿的相位因子，分别记做和其中，exp(·)表示e指数运算符号，j为虚数单位，k_x、k_y和k_z为步骤8中计算得到的信号频率在跨航向、方位向和距离向的波数域离散表示，k＝1,2,…,K为步骤3.1中得到的阵列向慢时刻，K为步骤3.1中计算得到的阵列向慢时刻总数，P(1,k)为步骤3.1中计算得到的阵元位置矩阵P第1行、第k列的元素，P(2,k)为步骤3.1中计算得到的阵元位置矩阵P第2行、第k列的元素，R₀为步骤1中初始化得到的雷达平台高度；

采用公式计算得到补偿相位因子后的后向投影值，记做其中，k＝1,2,…,K为步骤3.1中得到的阵列向慢时刻，K为步骤3.1中计算得到的阵列向慢时刻总数，p为步骤9中计算得到的后向投影值的索引，为步骤4中SINC插值后的回波数据的第k行、第p列元素，∑(·)表示求和运算符号；

步骤11、计算三维成像结果：

采用标准傅里叶逆变换对步骤10中计算得到的补偿相位因子后的后向投影值进行三维傅里叶逆变换，得到线阵SAR观测场景目标空间Ω的三维成像结果。