CN102854507B - 一种基于gpu后向投影双站合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开了一种基于GPU后向投影双站合成孔径雷达成像方法，它通过优化了传统快速因式分解后向投影方法的分级参数的方法，与传统快速因式分解后向投影方法相比，提升了成像处理效率；另一方面，它通过将成像任务分离成若干个互不相关的子任务群，以便GPU对子任务群进行并行处理，并通过新颖的并行策略再次提升了处理效率。相较现有的双站SAR成像方法，本发明既保持了时域方法内存开销小、易进行运动补偿等优点，实现了大数据处理量的双站SAR快速成像，处理效率也大幅提升，同时还能够获得高质量的成像结果。

Description

一种基于GPU后向投影双站合成孔径雷达成像方法

技术领域

本技术发明属于雷达技术领域，它特别涉及了双站合成孔径雷达成像技术领域。

背景技术

双站合成孔径雷达（Bistatic synthetic aperture radar，简写为BSAR）是指收发天线分置于两个不同平台的SAR系统。与传统的单站SAR相比，双站SAR具有隐蔽性好，安全性高，抗干扰能力强的优点，并且能够实现一些单站SAR所无法实现的特殊应用模式，如前视成像，前视成像是一种非常具有应用价值的成像模式，可应用于导航、恶劣天气下的飞机导航及着陆等方面。鉴于双站SAR的多种优势，对双站SAR技术的研究具有重要意义。

快速因式分解后向投影方法（Fast Factorized Backprojection，简称为FFBP）是一种运算量低、适用于任意双站构型的SAR成像方法，它是后向投影方法（Backprojection，简称为BP）的一种快速实现。该方法首先对合成孔径进行子孔径分解，由于子孔径在方位向具有较低的分辨率，于是可以减少成像区域在方位向的采样点数，达到降低运算量的目的。详细请见文献Ulander,L.,Hellsten,H.,and Stenstrom,G.Synthetic aperture radar processing using fast factorizedback-projection.IEEE Transactions on Aerospace and Electronic Systems,39,3(July2003),760—776。然而在实时导航等一些实际应用中，快速因式分解后向投影方法还不能满足所要求的高实时性处理。

近年来，图形处理器（简称GPU）被引入SAR图像处理领域，由于GPU具有并行化处理数据的能力，而且硬件成本较低，所以可以用来提高SAR数据处理效率，如文献Ponce,O.,P.Prats,M.Rodriguez-Cassola,R.Scheiber,and A.Reigber,"Processing of circular SAR trajectories with fast factorizedback-projection,"Proc.IGARSS,3692-3695,Vancouver,Canada,2011.中便提到利用GPU进行数据处理，但是该文献并没有指出其具体实现过程。

发明内容

本发明的目的是针对现有双站SAR成像方法中存在的处理效率低的问题，提出改进快速因式分解后向投影方法与GPU相结合的一种基于GPU后向投影双站合成孔径雷达成像方法，它一方面优化了传统快速因式分解后向投影方法的分级参数，与传统快速因式分解后向投影方法相比，本发明提升了成像处理效率；另一方面，将成像任务分离成若干个互不相关的子任务群，以便GPU对子任务群进行并行处理，并通过新颖的并行策略再次提升了处理效率。相较现有的双站SAR成像方法，本发明成像处理效率大幅提升，能够满足双基地前视合成孔径雷达（Bistatic forward synthetic aperture radar，简写为BFSAR）的高实时性、高分辨要求。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、距离史、距离门

距离史是指接收机和发射机的天线相位中心到场景中散射点的距离之和。

距离门是指对应距离史的回波数据在整个回波数据中的位置。

定义2、雷达成像空间

雷达成像空间是指将场景空间中的散射点投影到距离向－方位向的二维空间坐标系，该空间由合成孔径雷达成像空间中的两个具有一定角度的坐标基确定。本发明中用以下数学关系表示成像空间M：

$M=\left\{\stackrel{\‾}{P}(v,u)\right|\stackrel{\‾}{P}(v,u)=v\·{\stackrel{\‾}{\mathrm{\ζ}}}_v+u\·{\stackrel{\‾}{\mathrm{\ζ}}}_u\};u,v\∈R$

其中和

表示构成成像空间M的坐标基，分别表示距离向和方位向，

为成像空间中的待观测点向量，u,v分别表示该点的距离和方位坐标,R表示实数。

定义3、合成孔径雷达成像场景参考点

合成孔径雷达成像场景参考点是指合成孔径雷达成像空间中的某个散射点，作为分析和处理场景中其他散射点的参照。

定义4、合成孔径雷达标准距离压缩方法

合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射参数，采用以下公式生成参考信号，并采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。

$f\left(t\right)=\exp (j\·\mathrm{\π}\·\frac{B}{T_r}\·t^2)t\∈[-\frac{T_r}{2},\frac{T_r}{2}]$

其中，j为虚数单位（即-1开根），f(t)为距离压缩参考函数，B为雷达发射基带信号的信号带宽，T_r为雷达发射信号脉冲宽度，t为时间变量，取值范围从

到

详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。

定义5、天线相位中心

天线相位中心是指雷达天线向外辐射信号的中心，本发明中天线相位中心指雷达天线的轨迹位置。

定义6、合成孔径与PRF时刻

合成孔径是指对于成像场景中的一个散射点从进入雷达波束照射范围至离开雷达波束照射范围的这段时间内，雷达波束中心所走过的长度。

雷达平台飞过一个合成孔径所需要的时间称为慢时间，雷达系统以一定的重复周期T_r发射接收脉冲，因此慢时间可以表示为一个以重复周期T_r为步长的离散化时间变量，其中每一个离散时间变量值为一个PRF时刻。

详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义7、标准辛格插值

标准辛格插值方法是指对于一个带限信号，在满足采样定理的情况下，采用卷积核为sinc的函数h(x)，h(x)的长度即窗长为W。

$h\left(x\right)=\sin c\left(x\right)=\frac{\sin \left(\mathrm{\πx}\right)}{\mathrm{\πx}}$

进行对已离散的信号g_d(i)插值，得到插值后所要的信号

$g\left(x\right)=\underset{i}{\mathrm{\Σ}}{g}_d\left(i\right)\sin c(x-i)$

详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。

定义8、子孔径与子图像

子孔径是指截取合成孔径部分长度得到的新合成孔径。

子图像是指对子孔径进行成像处理得到的雷达图像。

定义9、范数

设X是数域C上线性空间，称||·||为X上的范数(norm)，若它满足：1.正定性：||X||≥0，且||X||=0<=>X=0；2.齐次性：||aX||=|a|||X||；3.次可加性(三角不等式)：||X+Y||≤||X||+||Y||。若X=[x₁,x₂,…,x_N]^T为N×1维离散信号，向量XLP范数为 ${\left|\right|X\left|\right|}_P={\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\right|x^P\left|\right)}^{1/P},$ L1范数为 ${\left|\right|X\left|\right|}_1=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|x\right|,$ L2范数为 ${\left|\right|X\left|\right|}_2={\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|x\right|}^2\right)}^{\frac{1}{2}}.$

定义10、图形处理器（简称：GPU）与计算统一设备架构（简称：CUDA）

图形处理器GPU（Graphic Processing Unit）是指一种高度并行化的多核处理器，其特点是能够利用大量的处理单元进行并行计算，大大提高运算效率。

计算统一设备架构CUDA(Compute Unified Device Architecture)是指由NVIDIA公式提出的将GPU作为数据并行计算设备的软硬件体系，使得开发人员在不需要图形学相关知识的情况下，使用类C语言实现通用计算。

详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。

定义11、GPU内核函数（kernel）、线程块（block）和线程（thread）

由主机端（Host）调用在设备端（Device）运行的函数称为内核函数（kernel），一个内核函数是整个程序中的一个可以被并行执行的步骤。内核函数以线程网格（Grid）的形式组织，线程网格由若干个线程块（block）组成，而每个线程块又由若干个线程（thread）组成。

线程块（block）是内核函数执行并行计算的第一层并行单元，各线程块间无法通信，没有执行顺序。

线程（thread）是内核函数执行并行计算的第二层并行单元，同一线程块中的线程可以共享数据，没有执行顺序。

详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。

定义12、线程块索引与线程索引

线程块索引是指每个线程块在线程网格中的位置，是CUDA的内建变量，线程块的第一维索引为blockIdx.x，线程块的第二维索引为blockIdx.y，线程块的第三维索引恒为1。

线程索引是指每个线程在线程块中的位置，是CUDA的内建变量，线程的第一维索引为threadIdx.x，线程的第二维索引为threadIdx.y，线程的第三维索引为threadIdx.z。

详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。

定义13、GPU全局存储器

GPU全局存储器（global memory）是GPU片外的存储器，可读写，作用范围是整个程序，整个线程网格中的任意线程都能读写全局存储器的任意位置，存储空间大，生命周期为整个程序。

详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。

定义14、GPU内存分配方法和GPU数据拷贝方法

GPU内存分配方法：采用CUDA函数库中的函数cudaMalloc()实现在全局存储器中分配存储空间。

GPU数据拷贝方法：采用CUDA函数库中的函数cudaMemcpy()实现计算机内存与GPU全局存储器之间的数据复制。

详见文献NVIDIA CUDA C Programming Guide Version4.03/6/2011和“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。

定义15、取整函数

取整函数是指将任一实数映射到相近的整数的一类函数，本发明中用到三个取整函数，分别为ceil、floor和round。

ceil(x)为向上取整函数，表示取不超过x的最大一个整数。

floor(x)为向下取整函数，表示取不小于x的最小一个整数。

round(x)为近似取整函数，表示按四舍五入原则对x取整。

本发明提出一种基于GPU后向投影双站合成孔径雷达成像方法，它包括以下步骤：

步骤1：初始化参数

初始化成像系统参数包括：雷达系统工作的信号波长，记做λ，雷达系统发射信号带宽，记做B，雷达系统发射脉冲时宽，记做T_r，雷达系统采样频率，记做F_s，雷达系统脉冲重复频率，记做PRF，雷达发射平台速度矢量，记做雷达接收平台速度矢量，记做

雷达发射平台初始位置矢量，记做

雷达接收平台初始位置矢量，记做

场景参考点位置矢量，记做

雷达系统距离向采样点数，记做N_r，雷达系统方位向采样点数，记做N_a，光的传播速度，记做C，插值长度，记做W₀；

初始化场景参数：子图像距离向像素点间隔，记做d_r，第一级子图像方位向像素点间隔，记做

子图像距离向总的像素点数，记做s_r，第一级子图像方位向总的像素点数，记做

子孔径长度，记做b，即子孔径内有b个方位向采样点数；

雷达原始回波数据，记做

是一个二维的数据矩阵，雷达原始回波数据

数据矩阵每列数据是方位向回波信号的采样，每行的数据是逐个单脉冲距离向回波信号的采样；

线程块的第一维维度，记做N_blockx，线程块中线程的第一维索引，记做threadIdx.x，threadIdx.x=0,...,N_blockx，线程块的第二维维度，记做N_blocky，线程块中线程的第二维索引，记做threadIdx.y，threadIdx.y=0,...,N_blocky；初始化成像系统参数均由雷达系统提供，均为已知。

步骤2：计算天线相位中心历史位置

采用公式（1）计算得到雷达发射平台第n个PRF时刻的天线相位中心矢量 $\stackrel{\&OverBar;}{P_t}\left(n\right),$

$\stackrel{\&OverBar;}{P_t}\left(n\right)=\stackrel{\&OverBar;}{P_{t0}}+\stackrel{\&OverBar;}{V_t}\&CenterDot;n/\mathrm{PRF}---\left(1\right)$

公式（1）中，

是步骤1提供的雷达发射平台初始位置矢量，

是步骤1提供的雷达发射平台速度矢量，PRF是步骤1提供的雷达系统脉冲重复频率，n=1,...,N_a，n表示第n个PRF时刻，N_a是步骤1提供的雷达系统方位向采样点数；

采用公式（2）计算得到雷达接收平台第n个PRF时刻的天线相位中心矢量 $\stackrel{\&OverBar;}{P_r}\left(n\right),$

$\stackrel{\&OverBar;}{P_r}\left(n\right)=\stackrel{\&OverBar;}{P_{r0}}+\stackrel{\&OverBar;}{V_r}\&CenterDot;n/\mathrm{PRF}---\left(2\right)$

公式（2）中，

是步骤1提供的雷达接收平台初始位置矢量，

是步骤1提供的雷达接收平台速度矢量，PRF是步骤1提供的雷达系统脉冲重复频率，n=1,...,N_a，n表示第n个雷达系统脉冲重复频率PRF时刻，N_a是步骤1提供的雷达系统方位向采样点数；发射平台所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做P_{t_size}，接收平台所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做P_{r_size}；

步骤3：雷达原始回波数据距离压缩

对步骤1中的雷达原始回波数据采用传统的合成孔径雷达标准距离压缩方法进行压缩，得到距离压缩后的数据

占用的存储空间大小记做S_{RC_size}；

步骤4：为图形处理器GPU准备数据

步骤4.1：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P_{t_size}的存储空间，记做

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P_{r_size}的存储空间，记做

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为S_{RC_size}的存储空间，记做

步骤4.2：采用传统的GPU数据拷贝方法，将步骤2中得到雷达发射平台天线相位中心矢量

复制到步骤4.1得到的图形处理器GPU的存储空间

中；

采用传统的GPU数据拷贝方法，将步骤2中得到雷达接收平台天线相位中心矢量

复制到步骤4.1得到的图形处理器GPU的存储空间

中，n=1,...,N_a，N_a为雷达系统方位向采样点数；

采用传统的GPU数据拷贝方法，将步骤3中得到的距离压缩后的数据

复制到步骤4.1得到的图形处理器GPU的存储空间

中；

步骤5：第一级子图像生成

步骤5.1：采用公式（3）计算得到第一级子图像中像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)的位置矢量

${\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})=\stackrel{\&OverBar;}{P_0}+(r^{\left(1\right)}-1)\&CenterDot;d_r\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a^{\left(1\right)}-1)\&CenterDot;d_a^{\left(1\right)}\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v---\left(3\right)$

公式（3）中，

是步骤1提供的场景参考点位置矢量，

和

表示构成雷达成像空间M的坐标基，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,s_r，s_r为子图像距离向总的像素点数，a⁽¹⁾表示像素点位于第一级子图像方位向的第a⁽¹⁾个位置，

为第一级子图像方位向总的像素点数，d_r是子图像距离向像素点间隔，

是第一级子图像方位向像素点间隔，第一级子图像所有像素点的位置矢量占用的存储空间大小记做P⁽¹⁾ _size；

步骤5.2：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P⁽¹⁾ _size的存储空间，记做

采用传统的GPU数据拷贝方法，将步骤5.1中提供的所有像素点的位置矢量

复制到图形处理器GPU的存储空间

中，a⁽¹⁾表示像素点位于第一级子图像方位向的第a⁽¹⁾个位置，为第一级子图像方位向总的像素点数，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,s_r，s_r为第一级子图像距离向总的像素点数；

步骤5.3：由步骤1提供的线程块第一维维度N_blockx和子图像距离向总的像素点数s_r，采用公式（4）得到线程网格的第一维维度

$N_{\mathrm{gridx}}^{\left(1\right)}=\mathrm{floor}((s_r+N_{\mathrm{blockx}}-1)/N_{\mathrm{blockx}})---\left(4\right)$

由步骤1提供的线程块第二维维度N_blocky和第一级子图像方位向总的像素点数

采用公式（5）得到线程网格的第二维维度

$N_{\mathrm{gridy}}^{\left(1\right)}=\mathrm{floor}((s_a^{\left(1\right)}+N_{\mathrm{blocky}}-1)/N_{\mathrm{blocky}})---\left(5\right)$

公式（5）中floor为向下取整函数，从而得到线程网格中线程块的第一维索引blockIdx.x⁽¹⁾，以及线程网格中线程块的第二维索引blockIdx.y⁽¹⁾， ${\mathrm{blockIdx}.y}^{\left(1\right)}=0,...,N_{\mathrm{gridy}}^{\left(1\right)};$

步骤5.4：将子图像像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)与线程网格中的线程一一对应，具体方法是：

对于像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)，将其与索引标识为blockIdx.x⁽¹⁾，blockIdx.y⁽¹⁾，threadIdx.x，threadIdx.y的线程对应，它们的对应关系为公式（6）和公式（7）：

a⁽¹⁾=threadIdx.x+blockIdx.x⁽¹⁾·N_blockx+1   （6）

r⁽¹⁾=threadIdx.y+blockIdx.y⁽¹⁾·N_blocky+1   （7）

公式（6）和公式（7）中，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，

为第一级子图像方位向总的像素点数，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,s_r，s_r为第一级子图像距离向总的像素点数，threadIdx.x为每个线程块中线程的第一维索引，threadIdx.x=0,...,N_blockx，N_blockx为每个线程块的第一维维度，threadIdx.y为每个线程块中线程的第二维索引，threadIdx.y=0,...,N_blocky，N_blocky为每个线程块的第二维维度，blockIdx.x⁽¹⁾是步骤5.3中得到的线程网格中线程块的第一维索引，blockIdx.y⁽¹⁾是步骤5.3中得到的线程网格中线程块的第二维索引；

步骤5.5：利用步骤2中提供的公式（1）

n=1+b·(j⁽¹⁾-1),+...,+b·j⁽¹⁾   （8）

得到发射平台的位置矢量

利用步骤2中提供的公式（2）

和公式（8），得到接收平台的位置矢量

其中，

是雷达发射平台初始位置矢量，

是雷达发射平台速度矢量，PRF是雷达系统脉冲重复频率，

是雷达接收平台初始位置矢量，是雷达接收平台速度矢量，j⁽¹⁾是子图像的索引，b为子孔径长度，n表示第n个PRF时刻；

再利用步骤5.1中得到的

采用公式（9）得到双基地雷达距离史R⁽¹⁾(n)，

$R^{\left(1\right)}\left(n\right)={\left|\right|\stackrel{\&OverBar;}{P_t}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})\left|\right|}_2---\left(9\right)$

公式（9）中，||·||₂为L₂范数，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，

为第一级子图像方位向总的像素点数，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,s_r，s_r为第一级子图像距离向总的像素点数；

步骤5.6：利用步骤5.5得到的双基地雷达距离史R⁽¹⁾(n)，采用公式（10）得到双基地雷达距离史R⁽¹⁾(n)对应的距离门位置ID(n)，

ID(n)=R⁽¹⁾(n)·F_s/C   （10）

公式（10）中，F_s是雷达系统采样频率，C是光的传播速度，在步骤3得到的距离压缩后的回波数据矩阵

中取出位于ID(n)-mm列的数据B(n)，其中，mm=-W₀/2,...,W₀/2，是一个整数序列，W₀为插值长度，采用标准辛格插值方法对取出的数据B(n)进行插值，得到插值重采样后的数据C(n)，n表示第n个PRF时刻，n=1+b·(j⁽¹⁾-1),...,b·j⁽¹⁾，j⁽¹⁾是子图像的索引，b是子孔径长度；

定义补偿相位因子为

将重采样后的数据C(n)与补偿相位因子K(n)相乘，得到相位补偿后的数据A(n)，其中j为虚数单位（即-1开根），R⁽¹⁾(n)为步骤5.5中得到的距离史，λ为雷达系统工作的信号波长，n表示第n个PRF时刻，n=1+b·(j⁽¹⁾-1),+...,+b·j⁽¹⁾，j⁽¹⁾是子图像的索引，b是子孔径长度；

n=1+b·(j⁽¹⁾-1),+...,+b·j⁽¹⁾，由于第n个雷达系统脉冲重复频率PRF时刻n是一个整数序列，因此得到的A(n)也是一个数据序列，

将数据序列A(n)中的所有数据相加，得到像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)的像素值矩阵

其中n表示第n个雷达系统脉冲重复频率PRF时刻，j⁽¹⁾是子图像的索引，b是子孔径长度，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，

为第一级子图像方位向总的像素点数，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,s_r，s_r为第一级子图像距离向总的像素点数；

步骤6：子图像融合

步骤6.1：采用公式（11）计算得到第二级子图像方位向像素点间隔

$d_a^{\left(2\right)}=d_a^{\left(1\right)}/b---\left(11\right)$

采用公式（12）计算得到第二级子图像方位向总的像素点数

$s_a^{\left(2\right)}=s_a^{\left(1\right)}\&CenterDot;b---\left(12\right)$

公式（11）、公式（12）中，

为第一级子图像方位向像素点间隔，b为子孔径长度，

为第一级子图像方位向总的像素点数；

采用公式（13）计算得到第二级子图像中像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)的位置矢量 ${\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)}),$

${\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})=\stackrel{\&OverBar;}{P_0}+(r^{\left(2\right)}-1)\&CenterDot;d_r\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a^{\left(2\right)}-1)\&CenterDot;d_a^{\left(2\right)}\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v---\left(13\right)$

采用公式（13）中，

是场景参考点位置矢量，

和

是定义2中定义的成像空间的坐标基，d_r是子图像距离向像素点间隔，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,s_r，s_r为子图像距离向总的像素点数，a⁽²⁾表示像素点位于第二级子图像方位向的第a⁽²⁾个位置，

为第二级子图像方位向总的像素点数，第二级子图像所有像素点的位置矢量占用的存储空间大小记做P⁽²⁾ _size；

步骤6.2：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P⁽²⁾ _size的存储空间，记做

采用传统的GPU数据拷贝方法，将步骤6.1中得到的所有像素点的位置矢量

复制到图形处理器GPU的存储空间

中，

为第二级子图像方位向总的像素点数，r⁽²⁾=1,...,s_r，s_r为子图像距离向总的像素点数；

步骤6.3：由步骤1提供的线程块第一维维度N_blockx和子图像距离向总的像素点数s_r，采用公式（14）得到线程网格的第一维维度

$N_{\mathrm{gridx}}^{\left(2\right)}=\mathrm{floor}((s_r+N_{\mathrm{blockx}}-1)/N_{\mathrm{blockx}})---\left(14\right)$

由步骤1提供的线程块第二维维度N_blocky和步骤6.1中提供的第二级子图像方位向总的像素点数

采用公式（15）得到线程网格的第二维维度

$N_{\mathrm{gridy}}^{\left(2\right)}=\mathrm{floor}((s_a^{\left(2\right)}+N_{\mathrm{blocky}}-1)/N_{\mathrm{blocky}})---\left(15\right)$

公式（14）、公式（15）中，floor为向下取整函数，从而得到线程网格中线程块的第一维索引blockIdx.x⁽²⁾，

以及线程网格中线程块的第二维索引blockIdx.y⁽²⁾， ${\mathrm{blockIdx}.y}^{\left(2\right)}=0,...,N_{\mathrm{gridy}}^{\left(2\right)};$

步骤6.4：将子图像像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)与线程网格中的线程一一对应，具体方法是：对于像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)，将其与索引标识为blockIdx.x⁽²⁾，blockIdx.y⁽²⁾，threadIdx.x，threadIdx.y的线程对应，它们的对应关系为公式（16）和公式（17）：

a⁽²⁾=threadIdx.x+blockIdx.x⁽²⁾·N_blockx+1   （16）

r⁽²⁾=threadIdx.y+blockIdx.y⁽²⁾·N_blocky+1   （17）

公式（16）和公式（17）中，a⁽²⁾表示像素点位于子图像方位向的第a⁽²⁾个位置，

为第二级子图像方位向总的像素点数，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,s_r，s_r为第二级子图像距离向总的像素点数，threadIdx.x为每个线程块中线程的第一维索引，threadIdx.x=0,...,N_blockx，N_blockx为每个线程块的第一维维度，threadIdx.y为每个线程块中线程的第二维索引，threadIdx.y=0,...,N_blocky，N_blocky为每个线程块的第二维维度，blockIdx.x⁽²⁾是步骤6.3中得到的线程网格中线程块的第一维索引，blockIdx.y⁽²⁾是步骤6.3中得到的线程网格中线程块的第二维索引；

步骤6.5：将n=round(b/2)+(j⁽¹⁾-1)·b代入步骤2中提供的公式（1）中，

$\stackrel{\&OverBar;}{P_t}\left(n\right)=\stackrel{\&OverBar;}{P_{t0}}+\stackrel{\&OverBar;}{V_t}\&CenterDot;n/\mathrm{PRF}---\left(1\right)$

得到发射平台的位置矢量

将n=round(b/2)+(j⁽¹⁾-1)·b代入步骤2中提供的公式（2）中，

$\stackrel{\&OverBar;}{P_r}\left(n\right)=\stackrel{\&OverBar;}{P_{r0}}+\stackrel{\&OverBar;}{V_r}\&CenterDot;n/\mathrm{PRF}---\left(2\right)$

得到接收平台的位置矢量

其中，

是雷达发射平台初始位置矢量，

是雷达发射平台速度矢量，PRF是雷达系统脉冲重复频率，

是雷达接收平台初始位置矢量，

是雷达接收平台速度矢量，round为近似取整函数，j⁽¹⁾步骤5得到的

的索引，b为子孔径长度，n表示第n个PRF时刻，再利用步骤6.1中得到的

采用公式 $R^{\left(2\right)}\left(n\right)={\left|\right|\stackrel{\&OverBar;}{P_t}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2$ 得到双基地雷达距离史R⁽²⁾(n)，采用公式 $\mathrm{RR}={\left|\right|\stackrel{\&OverBar;}{P_t}\left(m\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(m\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2$ 得到参考距离RR，其中，m=round(N_a/2)，是合成孔径的中心位置，||·||₂为L₂范数，N_a是雷达系统方位向采样点数，n表示第n个PRF时刻，a⁽²⁾表示像素点位于子图像方位向的第a⁽²⁾个位置，

为第二级子图像方位向总的像素点数，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,s_r，s_r为第二级子图像距离向总的像素点数；

步骤6.6：从步骤5得到的

中取出位于第aa行，第rr列的数据

其中，aa=round(a⁽²⁾/b)，表示像素点位于

的第aa行，rr=round(r⁽²⁾/b)，表示像素点位于

的第rr列，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，为第一级子图像方位向总的像素点数，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,s_r，s_r为第一级子图像距离向总的像素点数，round为近似取整函数，j⁽¹⁾是

的索引；

定义补偿相位因子 $K^{\left(2\right)}\left(n\right)=\exp \left\{j2\mathrm{\&pi;}\frac{R^{\left(2\right)}\left(n\right)-\mathrm{RR}}{\mathrm{\&lambda;}}\right\},$ 将取出的数据

与补偿相位因子K⁽²⁾(n)相乘，得到相位补偿后的数据序列A⁽²⁾(n)，其中j为虚数单位（即-1开根），R⁽²⁾(n)为步骤6.5中得到的距离史，RR是步骤6.5中得到的参考距离，λ为雷达系统工作的信号波长，n表示第n个PRF时刻；

将得到的相位补偿后的数据序列A⁽²⁾(n)中的所有数据相加，得到像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)的像素值矩阵

其中n表示第n个PRF时刻，a⁽²⁾表示像素点位于子图像方位向的第a⁽²⁾个位置，

为第二级子图像方位向总的像素点数，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,s_r，s_r为第二级子图像距离向总的像素点数；

步骤7：数据输出

采用传统的GPU数据拷贝方法，将步骤6.6中得到的成像后的数据从图形处理器GPU中导出，其中，

为第二级子图像方位向总的像素点数，r⁽²⁾=1,...,s_r，s_r为第二级子图像距离向总的像素点数。

本发明的创新点在于针对当前的双站SAR成像方法所具有的处理效率低下，实时性差的缺点，提出了一种新的双站SAR成像方案，此方法将快速因式分解后向投影方法与GPU并行实现相结合，实现了双站SAR快速成像。

本发明的优点在于利用了较易获得的硬件，实现了大数据处理量的双站SAR快速成像，与现有的方法相比，该方法既保持了时域方法内存开销小，易进行运动补偿等优点，而且处理效率也大幅提升，同时还能够获得高质量的成像结果。

附图说明

图1为本发明具体实施方式采用的双站合成孔径雷达飞行几何关系图

其中Pt为发射平台在坐标系中的飞行轨迹；Pr为接收平台在坐标系中的飞行轨迹；P表示场景中的待测点；X、Y、Z表示场景坐标轴；a_t表示发射平台雷达波束指向单位向量；a_r表示接收平台雷达波束指向单位向量；Vt表示发射机的速度矢量；Vr表示接收机的速度矢量。

图2是步骤5的流程示意图

图3为步骤6的流程示意图

图4为本发明的流程示意图

图5是本发明所使用的仿真参数列表

具体实施方式

本发明主要采用仿真实验的方法进行验证该方案的可行性，所有步骤、结论都在MATLAB7.0，Visual Studio2010和GeForce GT240显卡上验证正确。具体实施步骤如下：

步骤1：初始化参数

参数详见附图5；

步骤2：计算天线相位中心历史位置

采用公式

计算得到雷达发射平台第n个PRF时刻的天线相位中心矢量采用公式计算得到雷达接收平台第n个PRF时刻的天线相位中心矢量

n表示第n个PRF时刻，发射平台所有1024个PRF时刻的天线相位中心矢量占用的存储空间大小，记做P_{t_size}=3·1024·8byte，接收平台所有1024个PRF时刻的天线相位中心矢量占用的存储空间大小，记做P_{r_size}=3·1024·8byte；

步骤3：雷达原始回波数据距离压缩

对步骤1中的雷达原始回波数据

采用传统的合成孔径雷达标准距离压缩方法进行压缩，得到距离压缩后的数据

占用的存储空间大小记做S_{RC_size}=1024·2048·8·2byte；

步骤4：为图形处理器GPU准备数据

步骤4.1：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P_{t_size}的存储空间，记做

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P_{r_size}的存储空间，记做

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为S_{RC_size}的存储空间，记做

步骤4.2：采用传统的GPU数据拷贝方法，将步骤2中得到雷达发射平台所有1024个PRF时刻的天线相位中心矢量

复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤2中得到雷达接收平台所有1024个PRF时刻的天线相位中心矢量

复制到图形处理器GPU的存储空间

中，n=1,...,1024；采用传统的GPU数据拷贝方法，将步骤3中得到的距离压缩后的数据

复制到图形处理器GPU的存储空间

中；

步骤5：第一级子图像生成

步骤5.1：采用公式 ${\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})=\stackrel{\&OverBar;}{P_0}+(r^{\left(1\right)}-1)\&CenterDot;0.5\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a^{\left(1\right)}-1)\&CenterDot;6.4\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v$ 计算得到第一级子图像中像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)的位置矢量

r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,100，a⁽¹⁾表示像素点位于第一级子图像方位向的第a⁽¹⁾个位置，a⁽¹⁾=1,...,10，第一级子图像所有像素点的位置矢量占用的存储空间大小记做P⁽¹⁾ _size=10·100·8·2byte；

步骤5.2：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P⁽¹⁾ _size的存储空间，记做

采用传统的GPU数据拷贝方法，将步骤5.1中得到的所有像素点的位置矢量复制到图形处理器GPU的存储空间中，a⁽¹⁾=1,...,10，r⁽¹⁾=1,...,100；

步骤5.3：由步骤1提供的线程块第一维维度N_blockx=8和子图像距离向总的像素点数s_r=100，采用公式 $N_{\mathrm{gridx}}^{\left(1\right)}=\mathrm{floor}((s_r+N_{\mathrm{blockx}}-1)/N_{\mathrm{blockx}})$ 得到线程网格的第一维维度

由步骤1提供的线程块第二维维度N_blocky=8和第一级子图像方位向总的像素点数

采用公式 $N_{\mathrm{gridy}}^{\left(1\right)}=\mathrm{floor}((s_a^{\left(1\right)}+N_{\mathrm{blocky}}-1)/N_{\mathrm{blocky}})$ 得到线程网格的第二维维度 $N_{\mathrm{gridy}}^{\left(1\right)}=2,$ 其中floor为向下取整函数，从而得到线程网格中线程块的第一维索引blockIdx.x⁽¹⁾，blockIdx.x⁽¹⁾=0,...,12，以及线程网格中线程块的第二维索引blockIdx.y⁽¹⁾，blockIdx.y⁽¹⁾=0,...,1；

步骤5.4：将子图像像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)与线程网格中的线程一一对应，方法是：对于像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)，将其与索引标识为blockIdx.x⁽¹⁾，blockIdx.y⁽¹⁾，threadIdx.x，threadIdx.y的线程对应，他们的对应关系为a⁽¹⁾=threadIdx.x+blockIdx.x⁽¹⁾·8+1，r⁽¹⁾=threadIdx.y+blockIdx.y⁽¹⁾·8+1，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，a⁽¹⁾=1,...,10，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,100，threadIdx.x为每个线程块中线程的第一维索引，threadIdx.x=0,...,7，threadIdx.y为每个线程块中线程的第二维索引，threadIdx.y=0,...,7；

步骤5.5：利用步骤2中提供的公式代入n=1+32·(j⁽¹⁾-1),...,32·j⁽¹⁾，得到发射平台的位置矢量

利用步骤2中提供的公式

代入n=1+32·(j⁽¹⁾-1),...,32·j⁽¹⁾，得到接收平台的位置矢量

j⁽¹⁾是子图像的索引j⁽¹⁾=1,...,32，再利用步骤5.1中得到的

采用公式 $R^{\left(1\right)}\left(n\right)={\left|\right|\stackrel{\&OverBar;}{P_t}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})\left|\right|}_2$ 得到双基地雷达距离史R⁽¹⁾(n)，||·||₂为L₂范数，n表示第n个PRF时刻，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，a⁽¹⁾=1,...,10，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,100；

步骤5.6：利用步骤5.5得到的双基地雷达距离史R⁽¹⁾(n)，采用公式ID(n)=R⁽¹⁾(n)·390/300得到双基地雷达距离史R⁽¹⁾(n)对应的距离门位置ID(n)，在步骤3得到的距离压缩后的回波数据矩阵中取出位于ID(n)-mm列的数据B(n)，其中，mm=-8,...,8，是一个整数序列，采用标准辛格插值方法对取出的数据B(n)进行插值，得到插值重采样后的数据C(n)，n表示第n个PRF时刻，n=1+32·(j⁽¹⁾-1),...,32·j⁽¹⁾，j⁽¹⁾是子图像的索引；

将重采样后的数据C(n)与补偿相位因子K(n)相乘，得到相位补偿后的数据序列A(n)，其中j为虚数单位（即-1开根），R⁽¹⁾(n)为步骤5.5中得到的距离史，n表示第n个PRF时刻，n=1+32·(j⁽¹⁾-1),...,32·j⁽¹⁾，j⁽¹⁾是子图像的索引；

将子孔径内的32个PRF时刻的相位补偿后的数据A(n)相加，n=1+32·(j⁽¹⁾-1),...,32·j⁽¹⁾，得到像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)的像素值矩阵

其中n表示第n个PRF时刻，j⁽¹⁾是子图像的索引，a⁽¹⁾表示像素点位于子图像方位向的第a⁽¹⁾个位置，a⁽¹⁾=1,...,10，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾=1,...,100；

步骤六：子图像融合

步骤6.1：采用公式计算得到第二级子图像方位向像素点间隔

采用公式

计算得到第二级子图像方位向总的像素点数

为第一级子图像方位向像素点间隔，为第一级子图像方位向总的散射点数，采用公式 ${\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})=\stackrel{\&OverBar;}{P_0}+(r^{\left(2\right)}-1)\&CenterDot;0.5\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a^{\left(2\right)}-1)\&CenterDot;0.2\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v$ 计算得到第二级子图像中像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)的位置矢量r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,100，a⁽²⁾表示像素点位于第二级子图像方位向的第a⁽²⁾个位置，a⁽²⁾=1,...,320，第二级子图像所有像素点的位置矢量占用的存储空间大小记做P⁽²⁾ _size=320·100·8·2byte；

步骤6.2：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P⁽²⁾ _size的存储空间，记做

采用传统的GPU数据拷贝方法，将步骤6.1中得到的所有像素点的位置矢量

复制到图形处理器GPU的存储空间

中，a⁽²⁾=1,...,320，r⁽²⁾=1,...,100；

步骤6.3：由步骤1中的N_blockx=8和s_r=100，采用公式 $N_{\mathrm{gridx}}^{\left(2\right)}=\mathrm{floor}((s_r+N_{\mathrm{blockx}}-1)/N_{\mathrm{blockx}})$ 得到线程网格的第一维维度 $N_{\mathrm{gridx}}^{\left(2\right)}=13,$ 由步骤1中的N_blocky=8和 $s_a^{\left(2\right)}=320,$ 采用公式 $N_{\mathrm{gridy}}^{\left(2\right)}=\mathrm{floor}((s_a^{\left(2\right)}+N_{\mathrm{blocky}}-1)/N_{\mathrm{blocky}})$ 得到线程网格的第二维维度

其中floor为向下取整函数，从而得到线程网格中线程块的第一维索引blockIdx.x⁽²⁾，blockIdx.x⁽²⁾=0,...,12，以及线程网格中线程块的第二维索引blockIdx.y⁽²⁾，blockIdx.y⁽²⁾=0,...,39；

步骤6.4：将子图像像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)与线程网格中的线程一一对应，方法是：对于像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)，将其与索引标识为blockIdx.x⁽²⁾，blockIdx.y⁽²⁾，threadIdx.x，threadIdx.y的线程对应，他们的对应关系为a⁽²⁾=threadIdx.x+blockIdx.x⁽²⁾·8+1，r⁽²⁾=threadIdx.y+blockIdx.y⁽²⁾·8+1，a⁽²⁾表示像素点位于子图像方位向的第a⁽²⁾个位置，a⁽²⁾=1,...,320，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,100，threadIdx.x为每个线程块中线程的第一维索引，threadIdx.x=0,...,7，threadIdx.y为每个线程块中线程的第二维索引，threadIdx.y=0,...,7；

步骤6.5：利用步骤2中提供的公式

代入n=round(32/2)+(j⁽¹⁾-1)·32，得到发射平台的位置矢量

利用步骤2中提供的公式

代入n=round(32/2)+(j⁽¹⁾-1)·32，得到接收平台的位置矢量

j⁽¹⁾是第一级子图像

的索引，j⁽¹⁾=1,...,32，采用公式 $R^{\left(2\right)}\left(n\right)={\left|\right|\stackrel{\&OverBar;}{P_t}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2$ 得到双站距离史R⁽²⁾(n)，采用公式 $\mathrm{RR}={\left|\right|\stackrel{\&OverBar;}{P_t}\left(m\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(m\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2$ 得到参考距离RR，m是合成孔径中心位置，m=512，round为近似取整函数，||·||₂为L₂范数，n表示第n个PRF时刻，a⁽²⁾表示像素点位于子图像方位向的第a⁽²⁾个位置，a⁽²⁾=1,...,320，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,100；

步骤6.6：采用公式aa=round(a⁽²⁾/32)，rr=round(r⁽²⁾/32)从第一级子图像

中取出第aa行，第rr列的数据

其中，a⁽¹⁾=1,...,10，r⁽¹⁾=1,...,100，round为近似取整函数，j⁽ⁱ⁾是第一级子图像

的索引，j⁽¹⁾=1,...,32；

将取出的数据

与补偿相位因子K⁽²⁾(n)相乘， $K^{\left(2\right)}\left(n\right)=\exp \left\{j2\mathrm{\&pi;}\frac{R^{\left(2\right)}\left(n\right)-\mathrm{RR}}{0.03}\right\},$ 得到相位补偿后的数据序列A⁽²⁾(n)，其中j为虚数单位（即-1开根），R⁽²⁾(n)为步骤6.5中得到的距离史，n表示第n个PRF时刻，n=round(32/2)+(j⁽¹⁾-1)·32，j⁽¹⁾是第一级子图像的索引；

将得到的所有相位补偿后的数据A⁽²⁾(n)相加，得到像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)的像素值矩阵

其中n表示第n个PRF时刻，a⁽²⁾表示像素点位于子图像方位向的第a⁽²⁾个位置，a⁽²⁾=1,...,320，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾=1,...,100；

步骤七：数据输出

采用传统的GPU数据拷贝方法，将步骤6.6中得到的成像后的数据a⁽²⁾=1,...,320，r⁽²⁾=1,...,100，从图形处理器GPU中导出。

Claims (1)

1.一种基于GPU后向投影双站合成孔径雷达成像方法，其特征是它包括以下步骤：

步骤1：初始化参数

初始化成像系统参数包括：雷达系统工作的信号波长，记做λ，雷达系统发射信号带宽，记做B，雷达系统发射脉冲时宽，记做T_r，雷达系统采样频率，记做F_s，雷达系统脉冲重复频率，记做PRF，雷达发射平台速度矢量，记做

雷达接收平台速度矢量，记做

雷达发射平台初始位置矢量，记做

雷达接收平台初始位置矢量，记做

场景参考点位置矢量，记做雷达系统距离向采样点数，记做N_r，雷达系统方位向采样点数，记做N_a，光的传播速度，记做C，插值长度，记做W₀；

初始化场景参数：子图像距离向像素点间隔，记做d_r，第一级子图像方位向像素点间隔，记做

子图像距离向总的像素点数，记做s_r，第一级子图像方位向总的像素点数，记做

子孔径长度，记做b，即子孔径内有b个方位向采样点数；

雷达原始回波数据，记做

是一个二维的数据矩阵，雷达原始回波数据

数据矩阵每列数据是方位向回波信号的采样，每行的数据是逐个单脉冲距离向回波信号的采样；

线程块的第一维维度，记做N_blockx，线程块中线程的第一维索引，记做threadIdx.x，threadIdx.x＝0,...,N_blockx，线程块的第二维维度，记做N_blocky，线程块中线程的第二维索引，记做threadIdx.y，threadIdx.y＝0,...,N_blocky；初始化成像系统参数均由雷达系统提供，均为已知；

步骤2：计算天线相位中心历史位置

采用公式（1）计算得到雷达发射平台第n个PRF时刻的天线相位中心矢量

$\stackrel{\&OverBar;}{P_t}\left(n\right)=\stackrel{\&OverBar;}{P_{t0}}+\stackrel{\&OverBar;}{V_t}\&CenterDot;n/\mathrm{PRF}---\left(1\right)$

公式（1）中，

是步骤1提供的雷达发射平台初始位置矢量，是步骤1提供的雷达发射平台速度矢量，PRF是步骤1提供的雷达系统脉冲重复频率，n＝1,...,N_a，n表示第n个PRF时刻，N_a是步骤1提供的雷达系统方位向采样点数；

采用公式（2）计算得到雷达接收平台第n个PRF时刻的天线相位中心矢量

$\stackrel{\&OverBar;}{P_r}\left(n\right)=\stackrel{\&OverBar;}{P_{r0}}+\stackrel{\&OverBar;}{V_r}\&CenterDot;n/\mathrm{PRF}---\left(2\right)$

公式（2）中，

是步骤1提供的雷达接收平台初始位置矢量，是步骤1提供的雷达接收平台速度矢量；发射平台所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做P_{t_size}，接收平台所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做P_{r_size}；

步骤3：雷达原始回波数据距离压缩

对步骤1中的雷达原始回波数据

采用传统的合成孔径雷达标准距离压缩方法进行压缩，得到距离压缩后的数据

占用的存储空间大小记做S_{RC_size}；

步骤4：为图形处理器GPU准备数据

步骤4.1：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P_{t_size}的存储空间，记做

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P_{r_size}的存储空间，记做

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为S_{RC_size}的存储空间，记做

步骤4.2：采用传统的GPU数据拷贝方法，将步骤2中得到雷达发射平台天线相位中心矢量

复制到步骤4.1得到的图形处理器GPU的存储空间

中；

采用传统的GPU数据拷贝方法，将步骤2中得到雷达接收平台天线相位中心矢量

复制到步骤4.1得到的图形处理器GPU的存储空间

中；

采用传统的GPU数据拷贝方法，将步骤3中得到的距离压缩后的数据

复制到步骤4.1得到的图形处理器GPU的存储空间

中；

步骤5：第一级子图像生成

步骤5.1：采用公式（3）计算得到第一级子图像中像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)的位置矢量

${\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})=\stackrel{\&OverBar;}{P_0}+(r^{\left(1\right)}-1)\&CenterDot;d_r\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a^{\left(1\right)}-1)\&CenterDot;d_a^{\left(1\right)}\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v---\left(3\right)$

公式（3）中，是步骤1提供的场景参考点位置矢量，

和

表示构成雷达成像空间M的坐标基，r⁽¹⁾表示像素点位于子图像距离向的第r⁽¹⁾个位置，r⁽¹⁾＝1,...,s_r，s_r为第一级子图像距离向总的像素点数，a⁽¹⁾表示像素点位于第一级子图像方位向的第a⁽¹⁾个位置，为第一级子图像方位向总的像素点数，d_r是子图像距离向像素点间隔，

是第一级子图像方位向像素点间隔，第一级子图像所有像素点的位置矢量占用的存储空间大小记做P⁽¹⁾ _size；

步骤5.2：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P⁽¹⁾ _size的存储空间，记做

采用传统的GPU数据拷贝方法，将步骤5.1中提供的所有像素点的位置矢量复制到图形处理器GPU的存储空间

中；

步骤5.3：由步骤1提供的线程块第一维维度N_blockx和子图像距离向总的像素点数s_r，采用公式（4）得到线程网格的第一维维度

$N_{\mathrm{gridx}}^{\left(1\right)}=\mathrm{floor}((s_r+N_{\mathrm{blockx}}-1)/N_{\mathrm{blockx}})---\left(4\right)$

由步骤1提供的线程块第二维维度N_blocky和第一级子图像方位向总的像素点数

采用公式（5）得到线程网格的第二维维度

$N_{\mathrm{gridy}}^{\left(1\right)}=\mathrm{floor}((s_a^{\left(1\right)}+N_{\mathrm{blocky}}-1)/N_{\mathrm{blocky}})---\left(5\right)$

公式（5）中floor为向下取整函数，从而得到线程网格中线程块的第一维索引blockIdx.x⁽¹⁾，以及线程网格中线程块的第二维索引blockIdx.y⁽¹⁾， ${\mathrm{blockIdx}.y}^{\left(1\right)}=0,...,N_{\mathrm{gridy}}^{\left(1\right)};$

步骤5.4：将子图像像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)与线程网格中的线程一一对应，具体方法是：

对于像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)，将其与索引标识为blockIdx.x⁽¹⁾，blockIdx.y⁽¹⁾，threadIdx.x，threadIdx.y的线程对应，它们的对应关系为公式（6）和公式（7）：

a⁽¹⁾＝thredIdx.x+blockIdx.x⁽¹⁾·N_blockx+1  (6)

r⁽¹⁾＝threadIdx.y+blockIdx.y⁽¹⁾·N_blocky+1         （7）

公式（6）和公式（7）中，threadIdx.x为每个线程块中线程的第一维索引，threadIdx.x＝0,...,N_blockx，N_blockx为每个线程块的第一维维度，threadIdx.y为每个线程块中线程的第二维索引，threadIdx.y＝0,...,N_blocky，N_blocky为每个线程块的第二维维度，blockIdx.x⁽¹⁾是步骤5.3中得到的线程网格中线程块的第一维索引，blockIdx.y⁽¹⁾是步骤5.3中得到的线程网格中线程块的第二维索引；

步骤5.5：利用步骤2中提供的公式（1）

n＝1+b·(j⁽¹⁾-1),+...,+b·j⁽¹⁾        （8）

得到雷达发射平台第n个PRF时刻的天线相位中心矢量

利用步骤2中提供的公式（2）

和公式（8），得到雷达接收平台第n个PRF时刻的天线相位中心矢量

其中，j⁽¹⁾是子图像的索引，b为子孔径长度；

再利用步骤5.1中得到的

采用公式（9）得到双基地雷达距离史R⁽¹⁾(n)，

$R^{\left(1\right)}\left(n\right)={\left|\right|\stackrel{\&OverBar;}{P_t}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(1\right)}(a^{\left(1\right)},r^{\left(1\right)})\left|\right|}_2---\left(9\right)$

公式（9）中，||·||₂为L₂范数；

步骤5.6：利用步骤5.5得到的双基地雷达距离史R⁽¹⁾(n)，采用公式（10）得到双基地雷达距离史R⁽¹⁾(n)对应的距离门位置ID(n)，

ID(n)＝R⁽¹⁾(n)·F_s/C         （10）

公式（10）中，F_s是雷达系统采样频率，C是光的传播速度，在步骤3得到的距离压缩后的回波数据矩阵

中取出位于ID(n)-mm列的数据B(n)，其中，mm＝-W₀/2,...,W₀/2，是一个整数序列，W₀为插值长度，采用标准辛格插值方法对取出的数据B(n)进行插值，得到插值重采样后的数据C(n)；

定义补偿相位因子为

，将重采样后的数据C(n)与补偿相位因子K(n)相乘，得到相位补偿后的数据A(n)，其中j为虚数单位，即-1开根，R⁽¹⁾(n)为步骤5.5中得到的距离史，λ为雷达系统工作的信号波长；

将数据序列A(n)中的所有数据相加，得到像素点P⁽¹⁾(a⁽¹⁾,r⁽¹⁾)的像素值矩阵

步骤6：子图像融合

步骤6.1：采用公式（11）计算得到第二级子图像方位向像素点间隔

$d_a^{\left(2\right)}=d_a^{\left(1\right)}/b---\left(11\right)$

采用公式（12）计算得到第二级子图像方位向总的像素点数

$s_a^{\left(2\right)}=s_a^{\left(1\right)}\&CenterDot;b---\left(12\right)$

采用公式（13）计算得到第二级子图像中像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)的位置矢量

${\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})=\stackrel{\&OverBar;}{P_0}+(r^{\left(2\right)}-1)\&CenterDot;d_r\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a^{\left(2\right)}-1)\&CenterDot;d_a^{\left(2\right)}\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v---\left(13\right)$

公式（13）中，

和

是成像空间的坐标基，d_r是子图像距离向像素点间隔，r⁽²⁾表示像素点位于子图像距离向的第r⁽²⁾个位置，r⁽²⁾＝1,...,s_r，s_r为第二级子图像距离向总的像素点数，a⁽²⁾表示像素点位于第二级子图像方位向的第a⁽²⁾个位置，为第二级子图像方位向总的像素点数，第二级子图像所有像素点的位置矢量占用的存储空间大小记做P⁽²⁾ _size；

步骤6.2：采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为P⁽²⁾ _size的存储空间，记做

采用传统的GPU数据拷贝方法，将步骤6.1中得到的所有像素点的位置矢量

复制到图形处理器GPU的存储空间

中；

步骤6.3：由步骤1提供的线程块第一维维度N_blockx和子图像距离向总的像素点数s_r，采用公式（14）得到线程网格的第一维维度

$N_{\mathrm{gridx}}^{\left(2\right)}=\mathrm{floor}((s_r+N_{\mathrm{blockx}}-1)/N_{\mathrm{blockx}})---\left(14\right)$

由步骤1提供的线程块第二维维度N_blocky和步骤6.1中提供的第二级子图像方位向总的像素点数

采用公式（15）得到线程网格的第二维维度

$N_{\mathrm{gridy}}^{\left(2\right)}=\mathrm{floor}((s_a^{\left(2\right)}+N_{\mathrm{blocky}}-1)/N_{\mathrm{blocky}})---\left(15\right)$

从而得到线程网格中线程块的第一维索引blockIdx.x⁽²⁾，

以及线程网格中线程块的第二维索引blockIdx.y⁽²⁾， ${\mathrm{blockIdx}.y}^{\left(2\right)}=0,...,N_{\mathrm{gridy}}^{\left(2\right)};$

步骤6.4：将子图像像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)与线程网格中的线程一一对应，具体方法是：对于像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)，将其与索引标识为blockIdx.x⁽²⁾，blockIdx.y⁽²⁾，threadIdx.x，threadIdx.y的线程对应，它们的对应关系为公式（16）和公式（17）：

a⁽²⁾＝thredIdx.x+blockIdx.x⁽²⁾·N_blockx+1   （16）

r⁽²⁾＝threadIdx.y+blockIdx.y⁽²⁾·N_blocky+1       （17）

公式（16）和公式（17）中，blockIdx.x⁽²⁾是步骤6.3中得到的线程网格中线程块的第一维索引，blockIdx.y⁽²⁾是步骤6.3中得到的线程网格中线程块的第二维索引；

步骤6.5：将n＝round(b/2)+(j⁽¹⁾-1)·b代入步骤2中提供的公式（1）中，

$\stackrel{\&OverBar;}{P_t}\left(n\right)=\stackrel{\&OverBar;}{P_{t0}}+\stackrel{\&OverBar;}{V_t}\&CenterDot;n/\mathrm{PRF}---\left(1\right)$

得到雷达发射平台第n个PRF时刻的天线相位中心矢量

将n＝round(b/2)+(j⁽¹⁾-1)·b代入步骤2中提供的公式（2）中，

$\stackrel{\&OverBar;}{P_r}\left(n\right)=\stackrel{\&OverBar;}{P_{r0}}+\stackrel{\&OverBar;}{V_r}\&CenterDot;n/\mathrm{PRF}---\left(2\right)$

得到雷达接收平台第n个PRF时刻的天线相位中心矢量

其中，round为近似取整函数，j⁽¹⁾为步骤5得到的

的索引，再利用步骤6.1中得到的

，采用公式 $R^{\left(2\right)}\left(n\right)={\left|\right|{\stackrel{\&OverBar;}{P}}_t\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2+{\left|\right|{\stackrel{\&OverBar;}{P}}_r\left(n\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2$ 得到双基地雷达距离史R⁽²⁾(n)，采用公式 $\mathrm{RR}={\left|\right|\stackrel{\&OverBar;}{P_t}\left(m\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2+{\left|\right|\stackrel{\&OverBar;}{P_r}\left(m\right)-{\stackrel{\&OverBar;}{P}}^{\left(2\right)}(a^{\left(2\right)},r^{\left(2\right)})\left|\right|}_2$ 得到参考距离RR，其中，m＝round(N_a/2)，是合成孔径的中心位置；

步骤6.6：从步骤5得到的

中取出位于第aa行，第rr列的数据

其中，aa＝round(a⁽²⁾/b)，表示像素点位于

的第aa行，rr＝round(r⁽²⁾/b)，表示像素点位于

的第rr列；

定义补偿相位因子 $K^{\left(2\right)}\left(n\right)=\exp \left\{j2\mathrm{\&pi;}\frac{R^{\left(2\right)}\left(n\right)-\mathrm{RR}}{\mathrm{\&lambda;}}\right\}$ ，将取出的数据

与补偿相位因子K⁽²⁾(n)相乘得到相位补偿后的数据序列A⁽²⁾(n)，其中R⁽²⁾(n)为步骤6.5中得到的距离史，RR是步骤6.5中得到的参考距离；

将得到的相位补偿后的数据序列A⁽²⁾(n)中的所有数据相加，得到像素点P⁽²⁾(a⁽²⁾,r⁽²⁾)的像素值矩阵

步骤7：数据输出

采用传统的GPU数据拷贝方法，将步骤6.6中得到的成像后的数据

从图形处理器GPU中导出。

Images