CN102788979B

CN102788979B - 一种基于后向投影InSAR成像配准的GPU实现方法

Info

Publication number: CN102788979B
Application number: CN201210252090.4A
Authority: CN
Inventors: 张晓玲; 韦顺军; 付涛
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2012-07-20
Filing date: 2012-07-20
Publication date: 2014-04-09
Anticipated expiration: 2032-07-20
Also published as: CN102788979A

Abstract

本发明公开了一种基于后向投影InSAR成像配准的GPU实现方法，它是利用后向投影算法将不同天线测量数据投影到同一个成像空间中进行成像，在成像过程中同时完成了InSAR不同天线的高保相成像和图像配准，实现了InSAR成像配准一体化，并结合GPU并行化处理技术实现后向投影算法的快速处理。相对于传统InSAR干涉相位生成，本发明具有保相性能好、成像配准一体化、并行化处理提高算法效率等优点，实现了高精度InSAR干涉相位的快速生成。本发明可以应用于合成孔径雷达成像，地球遥感等领域。

Description

一种基于后向投影InSAR成像配准的GPU实现方法

技术领域

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

背景技术

合成孔径雷达(SAR)是一种高分辨率的微波成像***。合成孔径雷达依靠雷达和目标之间的相对运动来合成虚拟阵列以获得方位向高分辨率，利用大带宽信号实现距离向高分辨率，可以对照射场景进行二维成像。合成孔径雷达干涉测量（InSAR）是一般SAR功能的延伸和扩展，利用两个或者多个位置不同的天线观测同一个目标场景，根据目标到不同天线的斜距差获得测量数据的干涉相位，再通过平台与地面观测场景的几何关系反演出地面场景的数字高程信息的技术。由于具有全天时、全天候的特点，InSAR已经成为当前提取大面积地表三维图像和地形高程变化信息的一项重要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。

高质量的干涉相位是InSAR获取高精度地形数字高程模型(DEM)的基础。如何获取高质量干涉相位是InSAR高程反演的关键技术，如文献‘An efficientsystem for creating synthetic InSAR images from simulations’,Xiaoru Yuan等，Pureand Applied Geophysics，2008年第165卷第3-4期。在InSAR数据处理过程中，干涉相位生成需要完成不同天线SAR数据二维成像和图像配准两个关键步骤。因此，二维成像算法的相位保持性能和图像配准精度直接决定了InSAR干涉相位的质量。目前，对于InSAR的干涉相位生成，现有处理方法是将数据二维成像和SAR图像配准两个处理过程分离进行，需先对InSAR数据进行二维成像处理，然后再通过SAR图像配准获得干涉相位。现有InSAR数据成像主要采用频域方法，如距离多普勒(RD)算法、变尺度(CS)算法、距离徙动(RM)算法等，对不同天线测量数据进行二维成像处理。虽然RD、CS和RM等频域成像算法计算复杂度小，能够实现大数据量的快速成像，但这些成像算法需要对距离压缩后的回波数据进行距离单元走动校正（RCMC），由于RCMC处理只利用目标到天线的近似斜距而不是真实斜距，造成成像结果的相位与真实相位存在一定的误差，降低了SAR图像相位保持精度。而且，在InSAR数据成像后，现有方法需要利用图像配准方法，如最大相关系数法、最大频谱法等，对不同天线SAR图像进行配准处理，将不同SAR图像中的相同散射点搬移到同一个像素。但是，利用配准方法进行InSAR图像配准往往不能实现完全精确配准，还增加了干涉相位生成处理的运算复杂度。

后向投影(BP)算法是一种基于时域相干处理的成像算法，其基本思想是通过计算成像区域内每一采样点到孔径长度内SAR天线平台之间的双程时延，然后将对应的时域回波信号进行相干累加，从而恢复出每个散射点的目标函数。与RD、CS和RM等频域处理方法不同，BP算法是一个信号相参积累的过程，不需要对数据进行距离徙动校正，直接利用场景目标到天线的精确斜距信息，所以可提高成像精度和保相精度。但是，相对于现有频域成像方法，BP算法运算量非常大，目前主要应用于SAR高精度二维成像，在InSAR数据成像还没有看到利用BP算法进行成像的报道。

GPU（Graphics Processing Unit，图形处理器）是一种高度并行化的多核处理器，它的特点是能够利用大量的处理单元进行并行计算，而CUDA(ComputeUnified Device Architecture，计算统一设备架构)则是由NVIDIA公式于2006年提出的利用GPU实现通用计算的类C编程模型，其自身带有特定的函数库（如空间分配、数据传输函数、CUFFT函数库等），见文献NVIDIA.CUDA CUFFTlibrary PG-05327-040_V01,February,2011和NVIDIA CUDA C ProgrammingGuide Version4.03/6/2011。

在InSAR成像中，由于不同天线观测同一个目标场景，不同天线具有相同的成像空间，因此可以利用后向投影技术将不同天线的测量数据投影到同一个成像空间进行成像，从而在数据成像的同时实现不同天线SAR图像的精确配准。由于后向投影算法是逐点累加成像，不同点之间的成像处理是可以并行的，所以利用CUDA语言以及基于GPU技术来实现算法的并行化，可以大大提高后向投影算法的运算效率。

发明内容

为了提高InSAR成像数据的相位保持性能，实现高质量的InSAR干涉相位快速生成，并减少不同天线图像配准过程对InSAR高程反演精度的影响，本发明提出了一种基于合成孔径雷达后向投影算法的InSAR成像配准一体化的GPU实现方法，即一种基于后向投影InSAR成像配准的GPU实现方法。该方法利用后向投影算法将不同天线测量数据投影到同一个成像空间中进行成像，在成像过程中同时完成了InSAR不同天线的高保相成像和图像配准，实现了InSAR成像配准一体化，并结合GPU并行化处理技术实现后向投影算法的快速处理。

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

定义1、数字高程模型(DEM)

数字高程模型(Digital Elevation Model，DEM)，是指利用一组有序数值阵列形式表示地表或地面高程的一种实体地面模型。本发明中DEM表示成一系列地面点的平面坐标X、Y和高程坐标Z组成的数据阵列。对于一个地面区域D，地形DEM表示为

DEM＝{D_i|(x_i,y_i,z_i),i∈D}

其中(x_i,y_i)是第i个地面像素点对应的平面坐标，z_i是对应的高程坐标。

定义2、合成孔径雷达干涉测量(InSAR)

合成孔径雷达干涉测量(InSAR)指利用两个或者两个以上的SAR数据中的相位信息进行相干处理，结合雷达参数和雷达几何位置信息反演地表三维及其变化信息的遥感技术，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义3、雷达成像空间

雷达成像空间是指将场景空间中的散射点投影到距离向－方位向的二维空间坐标系，该空间由合成孔径雷达成像空间中的两个相互正交的坐标基确定。目前典型合成孔径雷达的成像空间包括距离向－方位向投影空间。本发明中用以下数学关系表示成像空间M：

其中

和

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

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

定义4、距离史与距离门

距离史是指不同天线相位中心到场景中散射点的距离组成的序列。距离门是指对应距离史的回波数据在整个回波数据中的位置，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

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

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

定义6、图像配准

图像配准是指将不同时间、不同传感器(成像设备)或不同条件下(天候、照度、摄像位置和角度等)对同一个地区获取的两幅或多幅图像进行地理坐标的匹配，使不同图像中的相同目标或者特征点位于图像中同一个位置的过程。本发明中的图像配准特指将不同天线SAR图像空间中对应的相同目标进行匹配的过程。

定义7、合成孔径雷达后向投影算法

后向投影算法是基于匹配滤波原理的合成孔径雷达成像算法，其主要通过相干累加实现合成孔径雷达数据的聚焦成像。详细内容可参考文献：“Research on Anovel fast back projection algorithm for strip map bistatic SAR imaging”，HuangYulin等。

定义8、范数

设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范数为

{| | X | |}_{P} = {(Σ_{i = 1}^{N} | x^{P} |)}^{1 / P},

L1范数为

{| | X | |}_{1} = Σ_{i = 1}^{N} | x |,

L2范数为

{| | X | |}_{2} = {(Σ_{i = 1}^{N} {| x |}^{2})}^{\frac{1}{2}} .

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

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

f (t) = \exp (j \cdot π \cdot \frac{B}{T_{r}} \cdot t^{2}) t &Element; [- \frac{T_{r}}{2}, \frac{T_{r}}{2}]

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

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

定义10、合成孔径雷达原始回波仿真方法

合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定***参数条件下具有SAR回波信号特性的原始信号的方法，详细内容可参考文献：“InSAR回波信号与***仿真研究”，张剑琦，哈尔滨工业大学硕士论文。

定义11、天线相位中心

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

定义12、辛格插值（sinc插值）方法

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

h (x) = \sin c (x) = \frac{\sin (πx)}{πx}

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

g (x) = \underset{i}{Σ} g_{d} (i) \sin c (x - i)

详见文献“合成孔径雷达成像----算法与实现”，Frank H.Wong等编著，电子工业出版社出版。

定义13、取整函数

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

定义18、GPU全局存储器

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

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

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

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

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

本发明提供一种基于后向投影InSAR成像配准的GPU实现方法，该方法的实现步骤如下：

步骤1：初始化InSAR成像***参数

InSAR成像空间由InSAR成像空间中的两个相互正交的坐标基确定，定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基，记做

该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基

垂直的单位向量作为InSAR成像空间的第二个坐标基，记做

该坐标基方向为距离向。

InSAR雷达平台包含两组天线，即主天线和副天线，两组天线之间的距离为基线长，记做B_l，主天线发射脉冲信号，经过T_d时间的延迟，主天线和副天线同时接收回波延迟信号。雷达平台主天线接收的回波数据，记做

雷达平台副天线接收的回波数据，记做

其中和

均为二维矩阵，第一维均对应方位向，第二维均对应距离向，即二维矩阵和

的行存储的是方位向数据，二维矩阵

和

的列存储的是距离向数据。

初始化InSAR成像***参数包括：雷达***工作的信号波长，记做λ，雷达平台主天线发射信号带宽，记做B，雷达平台主天线发射脉冲时宽，记做T_r，雷达平台接收***采样频率，记做F_s,雷达***脉冲重复频率，记做prf，雷达平台主天线合成孔径长度，记做L_m，雷达平台副天线合成孔径长度，记做L_s，雷达平台速度矢量，记做

雷达平台主天线初始位置矢量，记做

雷达平台副天线初始位置矢量，记做

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

雷达***距离向采样点数，记做N_r，雷达***方位向采样点数，记做N_a，场景距离向散射点间隔，记做d_r，场景方位向散射点间隔，记做d_a，场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，记做R_mc，场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，记做R_sc，场景参考点在主天线回波数据和副天线回波数据中的距离门相同，距离门位置记做I_c，线程网格中线程块的第一维索引，记做blockIdx.x，线程网格中线程块的第二维索引，记做blockIdx.y，每个线程块中线程的第一维索引，记做threadIdx.x，每个线程块中线程的第二维索引，记做threadIdx.y，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量。

采用公式

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

采用公式

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

n表示第n个PRF时刻，主天线所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做

副天线所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做

采用公式

\overset{&OverBar;}{P} (a, r) = {\overset{&OverBar;}{P}}_{c} + (r - 1) \cdot d_{r} \cdot {\overset{&OverBar;}{ζ}}_{u} + (a - 1) \cdot d_{a} \cdot {\overset{&OverBar;}{ζ}}_{v}

计算得到场景散射点P(a,r)的位置矢量，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，场景所有散射点的位置矢量占用的存储空间大小记做为雷达平台主天线场景散射点P(a,r)成像后的数据，

初始化为0，a＝1,...,s_a，s_a为场景方位向总的散射点数，r＝1,...,s_r，s_r为场景距离向总的散射点数，雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小记做

为雷达平台副天线场景散射点P(a,r)成像后的数据，初始化为0，a＝1,...,s_a，s_a为场景方位向总的散射点数，r＝1,...,s_r，s_r为场景距离向总的散射点数，雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小记做

步骤2：InSAR原始回波数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据进行压缩，得到平台主天线距离压缩后数据，记做

的数据大小记做

采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据

进行压缩，得到平台副天线距离压缩后数据，记做

的数据大小记做

步骤3：为图形处理器GPU分配数据

采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为

的存储空间，记做

的存储空间，记做

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

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

。

采用传统的GPU数据拷贝方法，将步骤2中得到的平台主天线距离压缩后数据

从计算机内存复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤2中得到的平台副天线距离压缩后数据

从计算机内存复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤1中得到的主天线所有N_a个PRF时刻的天线相位中心矢量

n＝1,...,N_a，N_a为雷达***方位向采样点数，从计算机内存复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤1中得到的副天线所有N_a个PRF时刻的天线相位中心矢量

中；采用传统的GPU数据拷贝方法，将步骤1中得到的场景所有散射点的位置矢量a＝1,...,s_a，s_a为场景方位向总的散射点数，r＝1,...,s_r，s_r为场景距离向总的散射点数，从计算机内存复制到图形处理器GPU的存储空间

中。

步骤4：为图形处理器GPU内核函数划分线程网格

线程网格中的线程块以二维形式划分，线程块中的线程也以二维形式划分，线程网格第一维总的线程数等于场景方位向总的散射点数s_a，线程网格第二维总的线程数等于场景距离向总的散射点数s_r，即一个场景散射点对应一个线程，具体划分步骤是：

首先确定线程块中的线程数，方法是：Db_x为线程块第一维分配的线程数，Db_x设置为8的整数倍，同时不超过GPU线程硬件规定的范围；Db_y为线程块第二维分配的线程数，Db_y设置为8的整数倍，同时不超过GPU线程硬件规定的范围。

然后确定线程块的数目，方法是：Dg_x为线程网格第一维分配的线程块数目，Dg_x由场景方位向总的散射点数s_a和线程块第一维分配的线程数Db_x共同确定，计算公式为Dg_x＝floor((s_a+Db_x-1)/Db_x)；Dg_y为线程网格第二维分配的线程块数目，Dg_y由场景距离向总的散射点数s_r和线程块第二维分配的线程数Db_y共同确定，计算公式为Dg_y＝floor((s_r+Db_y-1)/Db_y)，其中floor为向下取整函数。

步骤5：平台主天线距离压缩后数据沿方位向成像

图形处理器GPU实现平台主天线距离压缩后数据

沿方位向成像的具体步骤是：

步骤5.1：确定场景散射点与线程网格中线程的对应关系，方法是：对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数，blockIdx.x为线程网格中线程块的第一维索引，blockIdx.y为线程网格中线程块的第二维索引，threadIdx.x为线程块中线程的第一维索引，threadIdx.y为线程块中线程的第二维索引，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量，散射点P(a,r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockIdx.x和线程块中线程的第一维索引threadIdx.x的对应关系为a＝threadIdx.x+blockIdx.x·Db_x+1，散射点P(a,r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockIdx.y和线程块中线程的第二维索引threadIdx.y的对应关系为r＝threadIdx.y+blockIdx.y·Db_y+1，其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数。

步骤5.2：l_m为主天线半个合成孔径内的PRF时刻个数，计算公式为l_m＝round(L_md_a2)，其中round为近似取整函数，L_m为主天线合成孔径长度，d_a为场景方位向散射点间隔，L_m和d_a均由步骤1初始化得到；i_m为主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，i_m的取值满足条件：i_m＝a-l_m,...,a+l_m，若a-l_m＜1，则a-l_m＝1，若a+l_m＞N_a，则a+l_m＝N_a，其中N_a为雷达***方位向采样点数，即i_m只取雷达***采样范围内的主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。采用公式

计算得到雷达平台主天线i_m时刻处散射点P(a,r)的距离史

其中||□||₂为L2范数，

为步骤1中得到的雷达平台主天线第i_m个PRF时刻的天线相位中心矢量，

为步骤1中得到的场景散射点P(a,r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数。

步骤5.3：

为场景散射点P(a,r)的距离史

对应的距离门，计算公式为

其中round为近似取整函数，

为步骤5.2中得到的雷达平台主天线i_m时刻处散射点P(a,r)的距离史，R_mc为场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，I_c为场景参考点在主天线回波数据中的距离门，d_r为场景距离向散射点间隔，R_mc、I_c和d_r均由步骤1初始化得到，i_m为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数。从步骤2中得到的平台主天线距离压缩后数据

的第i_m行中取第个数据的前W₀个数据和后W₀个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据

其中W₀为标准辛格插值的半个窗长。

步骤5.4：根据补偿相位因子的计算公式

得到散射点P(a,r)在PRF时刻i_m应补偿的相位因子

其中j为虚数单位（即-1开根），为步骤5.2中得到的雷达平台主天线i_m时刻处散射点P(a,r)的距离史，λ为雷达***工作的信号波长，λ由步骤1初始化得到，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数。将步骤5.3中得到的插值重采样后的数据

与散射点P(a,r)在PRF时刻i_m应补偿的相位因子

相乘，得到相位补偿后的数据

把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据

相加，得到雷达平台主天线场景散射点P(a,r)成像后的数据

即

其中i_m为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。

步骤5.5：采用传统的GPU数据拷贝方法，将步骤5.4中得到的雷达平台主天线场景所有散射点成像后的数据a＝1,...,s_a，s_a为场景方位向总的散射点数，r＝1,...,s_r，s_r为场景距离向总的散射点数，从图形处理器GPU的

复制到计算机内存中，其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。

步骤6：平台副天线距离压缩后数据沿方位向成像

图形处理器GPU实现平台副天线距离压缩后数据沿方位向成像的具体步骤是：

步骤6.1：确定场景散射点与线程网格中线程的对应关系，方法是：对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数，blockIdx.x为线程网格中线程块的第一维索引，blockIdx.y为线程网格中线程块的第二维索引，threadIdx.x为线程块中线程的第一维索引，threadIdx.y为线程块中线程的第二维索引，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量，散射点P(a,r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockIdx.x和线程块中线程的第一维索引threadIdx.x的对应关系为a＝threadIdx.x+blockIdx.x·Db_x+1，散射点P(a,r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockIdx.y和线程块中线程的第二维索引threadIdx.y的对应关系为r＝threadIdx.y+blockIdx.y·Db_y+1，其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数。

步骤6.2：l_s为副天线半个合成孔径内的PRF时刻个数，计算公式为l_s＝round(L_sd_a2)，其中round为近似取整函数，L_s为副天线合成孔径长度，d_a为场景方位向散射点间隔，L_s和d_a均由步骤1初始化得到；i_s为副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，i_s的取值满足条件：i_s＝a-l_s,...,a+l_s，若a-l_s＜1，则a-l_s＝1，若a+l_s＞N_a，则a+l_s＝N_a，其中N_a为雷达***方位向采样点数，即i_s只取雷达***采样范围内的副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。采用公式

计算得到雷达平台副天线i_s时刻处散射点P(a,r)的距离史

采用公式

计算得到雷达平台主天线i_s时刻处散射点P(a,r)的距离史

其中||□||₂为L2范数，

为步骤1中得到的雷达平台副天线第i_s个PRF时刻的天线相位中心矢量，

为步骤1中得到的雷达平台主天线第i_s个PRF时刻的天线相位中心矢量，

步骤6.3：

为场景散射点P(a,r)的距离史

对应的距离门，计算公式为

{\overset{&OverBar;}{I}}_{s} (i_{s}; a, r) = ((I_{c} + round {\overset{&OverBar;}{R}}_{s} (i_{s}; a, r) - R_{sc}) / d_{r}),

其中round为近似取整函数，

为步骤6.2中得到的雷达平台副天线i_s时刻处散射点P(a,r)的距离史，R_sc为场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，I_c为场景参考点在主天线回波数据中的距离门，d_r为场景距离向散射点间隔，R_sc、I_c和d_r均由步骤1初始化得到，i_s为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数。从步骤2中得到的平台副天线距离压缩后数据

的第i_s行中取第

个数据的前W₀个数据和后W₀个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据

其中W₀为标准辛格插值的半个窗长。

步骤6.4：根据补偿相位因子的计算公式

{\overset{&OverBar;}{F}}_{s} (i_{s}; a, r) = \exp {j \cdot 4 \cdot π \cdot {\overset{&OverBar;}{R}}_{m} (i_{s}; a, r) / λ},

得到散射点P(a,r)在PRF时刻i_s应补偿的相位因子其中j为虚数单位（即-1开根），

为步骤6.2中得到的雷达平台主天线i_s时刻处散射点P(a,r)的距离史，λ为雷达***工作的信号波长，λ由步骤1初始化得到，a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数。将步骤6.3中得到的插值重采样后的数据

与散射点P(a,r)在PRF时刻i_s应补偿的相位因子相乘，得到相位补偿后的数据

把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据

相加，得到雷达平台副天线场景散射点P(a,r)成像后的数据即其中i_s为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。

步骤6.5：采用传统的GPU数据拷贝方法，将步骤6.4中得到的雷达平台副天线场景所有散射点成像后的数据a＝1,...,s_a，s_a为场景方位向总的散射点数，r＝1,...,s_r，s_r为场景距离向总的散射点数，从图形处理器GPU的

复制到计算机内存中，其中

为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。

步骤7：将步骤5中得到的雷达平台主天线场景所有散射点成像后的数据

和步骤6中得到的雷达平台副天线场景所有散射点成像后的数据

进行共轭相乘，得到InSAR成像的干涉相位

完成InSAR数据成像配准一体化处理，其中a表示散射点位于场景方位向的第a个位置，a＝1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r＝1,...,s_r，s_r为场景距离向总的散射点数。

本发明的创新点在于利用后向投影算法，将InSAR不同天线测量数据投影到同一个成像空间进行成像，实现InSAR数据成像和图像配准一体化，并根据后向投影算法逐点相干积累原理，利用GPU并行化技术实现算法快速处理。通过仿真实验验证，相对于传统InSAR干涉相位生成，该算法具有保相性能好、成像配准一体化、并行化处理提高算法效率等优点，实现了高精度InSAR干涉相位的快速生成。

本发明的优点在于在数据成像过程中同时完成InSAR高保相成像与配准，在数据成像时同时获得高质量干涉相位，不需要再进行图像配准处理，并且可以利用GPU并行化处理提高算法效率。本发明可以应用于合成孔径雷达成像，地球遥感等领域。

附图说明

图1为本发明所提供方法的流程示意框图；

图2为本发明具体实施方式采用的干涉合成孔径雷达平台飞行几何关系及仿真圆锥场景图；

其中，横坐标（X轴）为方位向，纵坐标（Y轴）为距离向，垂直坐标（Z轴）为高度向，B_l为雷达平台主天线和副天线之间的基线，基线的两个端点分别为雷达平台主天线和副天线，P(a,r)为场景中任一散射点，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512，

为雷达平台飞行速度矢量，

为雷达平台主天线第i_m个PRF时刻的天线相位中心矢量，

为雷达平台副天线第i_s个PRF时刻的天线相位中心矢量，

为雷达平台主天线i_m时刻处散射点P(a,r)的距离史，

为雷达平台副天线i_s时刻处散射点P(a,r)的距离史，i_m为雷达平台主天线距离散射点P(a,r)半个合成孔径内的PRF时刻，i_s为雷达平台副天线距离散射点P(a,r)半个合成孔径内的PRF时刻，s_a为仿真场景方位向总的散射点数，s_r为仿真场景距离向总的散射点数，d_a为仿真场景方位向散射点间隔，d_r为仿真场景距离向散射点间隔。

图3本发明具体实施方式采用的线程网格划分；

其中，Block为线程网格中线程块的二维划分，Thread为第一维索引为0、第二维索引为15的线程块的二维线程划分。

具体实施方式

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

步骤1：初始化InSAR成像***参数

选择与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基

{\overset{&OverBar;}{ζ}}_{v} = [\begin{matrix} 0 & 1 & 0 \end{matrix}],

该坐标基方向为方位向；选择InSAR成像空间的第二个坐标基

{\overset{&OverBar;}{ζ}}_{u} = [\begin{matrix} 1 & 0 & 0 \end{matrix}],

该坐标基方向为距离向。

初始化InSAR成像***参数包括：雷达***工作的信号波长λ＝0.03m，雷达平台主天线发射信号带宽B＝150MHz，雷达平台主天线发射脉冲时宽T_r＝1μs，雷达平台接收***采样频率F_s＝300MHz,雷达***脉冲重复频率prf＝500Hz，雷达***距离向采样点数N_r＝1024，雷达***方位向采样点数N_a＝2048，场景距离向散射点间隔d_r＝0.5m，场景方位向散射点间隔d_a＝0.3m，雷达平台主天线合成孔径长度L_m＝150m，雷达平台副天线合成孔径长度L_s＝150m，雷达平台速度矢量

\overset{&OverBar;}{V} = [\begin{matrix} 1 & 150 & 0 \end{matrix}],

速度的单位为m/s，雷达平台主天线初始位置矢量

{\overset{&OverBar;}{P}}_{m} (0) [\begin{matrix} 0 & 0 & 6000 \end{matrix}],

雷达平台副天线初始位置矢量

{\overset{&OverBar;}{P}}_{s} (0) = [\begin{matrix} 10 & 0 & 6000 \end{matrix}],

场景参考位置矢量

{\overset{&OverBar;}{P}}_{c} = [\begin{matrix} 8000 & 0 & 0 \end{matrix}],

位置的单位为m，场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值R_mc＝10000m，场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值R_sc＝9992，场景参考点在主天线回波数据和副天线回波数据中的距离门I_c＝1。

仿真场景为一立体圆锥，仿真场景方位向总的散射点数s_a＝512，仿真场景距离向总的散射点数s_r＝512，圆锥最小高度为0m，最大高度为50m。

采用传统的合成孔径雷达原始回波仿真方法，生成雷达平台主天线仿真回波数据

采用传统的合成孔径雷达原始回波仿真方法，生成雷达平台副天线仿真回波数据

采用公式

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

采用公式

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

n＝1,...,2048，主天线所有2048个PRF时刻的天线相位中心矢量占用的存储空间大小

副天线所有2048个PRF时刻的天线相位中心矢量占用的存储空间大小采用公式

\overset{&OverBar;}{P} (a, r) = {\overset{&OverBar;}{P}}_{c} + (r - 1) \cdot 0.5 \cdot {\overset{&OverBar;}{ζ}}_{u} + (a - 1) 0.3 \cdot {\overset{&OverBar;}{ζ}}_{v}

计算得到场景散射点P(a,r)的位置矢量，r表示散射点位于场景距离向的第r个位置，r＝1,...,512，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，场景所有散射点的位置矢量占用的存储空间大小为雷达平台主天线场景散射点P(a,r)成像后的数据，

初始化为0，a＝1,...,512，r＝1,...,512，雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小

为雷达平台副天线场景散射点P(a,r)成像后的数据，

初始化为0，a＝1,...,512，r＝1,...,512，雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小

步骤2：InSAR原始回波数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据进行压缩，得到平台主天线距离压缩后数据

的大小

进行压缩，得到平台副天线距离压缩后数据

的大小

步骤3：为图形处理器GPU分配数据

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

从计算机内存复制到图形处理器GPU的存储空间

从计算机内存复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤1中得到的主天线所有2048个PRF时刻的天线相位中心矢量

n＝1,...,2048，从计算机内存复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤1中得到的副天线所有2048个PRF时刻的天线相位中心矢量

n＝1,...,2048，从计算机内存复制到图形处理器GPU的存储空间

中；采用传统的GPU数据拷贝方法，将步骤1中得到的场景所有散射点的位置矢量

a＝1,...,512，r＝1,...,512，从计算机内存复制到图形处理器GPU的存储空间

中。

步骤4：为图形处理器GPU内核函数划分线程网格

线程块第一维分配的线程数Db_x＝16，线程块第二维分配的线程数Db_y＝16，线程网格第一维分配的线程块数目Dg_x＝32，线程网格第二维分配的线程块数目Dg_y＝32。

步骤5：平台主天线距离压缩后数据沿方位向成像

图形处理器GPU实现平台主天线距离压缩后数据

沿方位向成像的具体步骤是：

步骤5.1：确定场景散射点与线程网格中线程的对应关系，方法是：对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512，blockIdx.x为线程网格中线程块的第一维索引，blockIdx.y为线程网格中线程块的第二维索引，threadIdx.x为线程块中线程的第一维索引，threadIdx.y为线程块中线程的第二维索引，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量，散射点P(a,r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockIdx.x和线程块中线程的第一维索引threadIdx.x的对应关系为a＝threadIdx.x+blockIdx.x·16+1，散射点P(a,r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockIdx.y和线程块中线程的第二维索引threadIdx.y的对应关系为r＝threadIdx.y+blockIdx.y·16+1。

步骤5.2：l_m为主天线半个合成孔径内的PRF时刻个数，l_m＝250；i_m为主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，i_m的取值满足条件：i_m＝a-250,…,a+250，若a-250＜1，则a-l_m＝1，若a+250＞2048，则a+l_m＝2048。采用公式

计算得到雷达平台主天线i_m时刻处散射点P(a,r)的距离史

其中||□||₂为L2范数，为步骤1中得到的雷达平台主天线第i_m个PRF时刻的天线相位中心矢量，

为步骤1中得到的场景散射点P(a,r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a＝1,...,5，r表示散射点位于场景距离向的第r个位置，r＝1,...,512。

步骤5.3：

为场景散射点P(a,r)的距离史对应的距离门，计算公式为

)其中round为近似取整函数，

为步骤5.2中得到的雷达平台主天线i_m时刻处散射点P(a,r)的距离史，i_m为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512。从步骤2中得到的平台主天线距离压缩后数据

的第i_m行中取第个数据的前4个数据和后4个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据

步骤5.4：根据补偿相位因子的计算公式

{\overset{&OverBar;}{F}}_{m} (i_{m}; a, r) = \exp {j \cdot 4 \cdot π \cdot (i_{m}; a, r) / 0.03},

得到散射点P(a,r)在PRF时刻i_m应补偿的相位因子其中j为虚数单位（即-1开根），

为步骤5.2中得到的雷达平台主天线i_m时刻处散射点P(a,r)的距离史，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512。将步骤5.3中得到的插值重采样后的数据

与散射点P(a,r)在PRF时刻i_m应补偿的相位因子

相乘，得到相位补偿后的数据

相加，得到雷达平台主天线场景散射点P(a,r)成像后的数据

即

步骤5.5：采用传统的GPU数据拷贝方法，将步骤5.4中得到的雷达平台主天线场景所有散射点成像后的数据

从图形处理器GPU的复制到计算机内存中，其中

为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。

步骤6：平台副天线距离压缩后数据沿方位向成像

图形处理器GPU实现平台副天线距离压缩后数据

沿方位向成像的具体步骤是：

步骤6.1：确定场景散射点与线程网格中线程的对应关系，方法是：对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512，blockIdx.x为线程网格中线程块的第一维索引，blockIdx.y为线程网格中线程块的第二维索引，threadIdx.x为线程块中线程的第一维索引，threadIdx.y为线程块中线程的第二维索引，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量，散射点P(a,r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockIdx.x和线程块中线程的第一维索引threadIdx.x的对应关系为a＝threadIdx.x+blockIdx.x·16+1，散射点P(a,r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockIdx.y和线程块中线程的第二维索引threadIdx.y的对应关系为r＝threadIdx.y+blockIdx.y·16+1。

步骤6.2：l_s为副天线半个合成孔径内的PRF时刻个数，l_s＝250；i_s为副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，i_s的取值满足条件：i_s＝a-250,...,a+250，若a-250＜1，则a-250＝1，若a+250＞2048，则a+250＝2048。采用公式

计算得到雷达平台副天线i_s时刻处散射点P(a,r)的距离史

采用公式计算得到雷达平台主天线i_s时刻处散射点P(a,r)的距离史

其中||□||₂为L2范数，为步骤1中得到的雷达平台副天线第i_s个PRF时刻的天线相位中心矢量，

为步骤1中得到的场景散射点P(a,r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512。步骤6.3：为场景散射点P(a,r)的距离史对应的距离门，计算公式为

{\overset{&OverBar;}{I}}_{s} (i_{s}; a, r) = 1 + round (({\overset{&OverBar;}{R}}_{s} (i_{s}; a, r) - 9992) / 0.5),

其中round为近似取整函数，

为步骤6.2中得到的雷达平台副天线i_s时刻处散射点P(a,r)的距离史，i_s为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512。从步骤2中得到的平台副天线距离压缩后数据

的第i_s行中取第

个数据的前4个数据和后4个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据

步骤6.4：根据补偿相位因子的计算公式

{\overset{&OverBar;}{F}}_{s} (i_{s}; a, r) = \exp {j \cdot 4 \cdot π \cdot {\overset{&OverBar;}{R}}_{m} (i_{s}; a, r) / 0.03},

得到散射点P(a,r)在PRF时刻i_s应补偿的相位因子

其中j为虚数单位（即-1开根），

为步骤6.2中得到的雷达平台主天线i_s时刻处散射点P(a,r)的距离史，a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r表示散射点位于场景距离向的第r个位置，r＝1,...,512。将步骤6.3中得到的插值重采样后的数据

与散射点P(a,r)在PRF时刻i_s应补偿的相位因子

相乘，得到相位补偿后的数据

相加，得到雷达平台副天线场景散射点P(a,r)成像后的数据

即

其中i_s为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。

步骤6.5：采用传统的GPU数据拷贝方法，将步骤6.4中得到的雷达平台副天线场景所有散射点成像后的数据

从图形处理器GPU的

复制到计算机内存中，其中

为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。

进行共轭相乘，得到InSAR成像的干涉相位

完成InSAR数据成像配准一体化处理，其中a表示散射点位于场景方位向的第a个位置，a＝1,...,512，r＝1,...,512。

通过本发明具体实施方式的仿真及测试，本发明所提供的基于后向投影InSAR成像配准的GPU实现方法能够实现InSAR不同天线成像和图像配准同时完成，与现有的InSAR干涉相位生成方法相比，本发明克服现有方法数据成像和图像配准相互隔离的缺点，可直接从原始数据成像处理过程中获得配准后的干涉相位，减少数据处理流程并且提高了成像相位保持精度，为InSAR高质量干涉相位生成提供了一种新方法。同时，利用GPU并行化技术实现了算法的快速处理。

Claims

1.一种基于后向投影InSAR成像配准的GPU实现方法，其特征是它包括以下步骤：

步骤1：初始化InSAR成像***参数

InSAR成像空间由InSAR成像空间中的两个相互正交的坐标基确定，定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基，记做该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基

垂直的单位向量作为InSAR成像空间的第二个坐标基，记做

该坐标基方向为距离向；

InSAR雷达平台包含两组天线，即主天线和副天线，两组天线之间的距离为基线长，记做B_l，主天线发射脉冲信号，经过T_d时间的延迟，主天线和副天线同时接收回波延迟信号；雷达平台主天线接收的回波数据，记做

雷达平台副天线接收的回波数据，记做其中

和

均为二维矩阵，第一维均对应方位向，第二维均对应距离向，即二维矩阵

和的行存储的是方位向数据，二维矩阵

和

的列存储的是距离向数据；

雷达平台主天线初始位置矢量，记做

雷达平台副天线初始位置矢量，记做

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

雷达***距离向采样点数，记做N_r，雷达***方位向采样点数，记做N_a，场景距离向散射点间隔，记做d_r，场景方位向散射点间隔，记做d_a，场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，记做R_mc，场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，记做R_sc，场景参考点在主天线回波数据和副天线回波数据中的距离门相同，距离门位置记做I_c，线程网格中线程块的第一维索引，记做blockIdx.x，线程网格中线程块的第二维索引，记做blockIdx.y，每个线程块中线程的第一维索引，记做threadIdx.x，每个线程块中线程的第二维索引，记做threadIdx.y，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量；

采用公式

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

采用公式

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

n=1,...,N_a,n表示第n个PRF时刻，主天线所有N_a个PRF时刻的天线相位中心矢量占用的存储空间大小，记做

采用公式

\overset{&OverBar;}{P} (a, r) = {\overset{&OverBar;}{P}}_{c} + (r - 1) \cdot d_{r} \cdot {\overset{&OverBar;}{ζ}}_{u} + (a - 1) \cdot d_{a} \cdot {\overset{&OverBar;}{ζ}}_{v}

计算得到场景散射点P(a,r)的位置矢量，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，场景所有散射点的位置矢量占用的存储空间大小记做

为雷达平台主天线场景散射点P(a,r)成像后的数据，

初始化为0，a=1,...,s_a，s_a为场景方位向总的散射点数，r=1,...,s_r，s_r为场景距离向总的散射点数，雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小记做

为雷达平台副天线场景散射点P(a,r)成像后的数据，初始化为0，a=1,...,s_a，s_a为场景方位向总的散射点数，r=1,...,s_r，s_r为场景距离向总的散射点数，雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小记做

步骤2：InSAR原始回波数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据

进行压缩，得到平台主天线距离压缩后数据，记做

的数据大小记做

进行压缩，得到平台副天线距离压缩后数据，记做

的数据大小记做

步骤3：为图形处理器GPU分配数据

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

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

的存储空间，记做

从计算机内存复制到图形处理器GPU的存储空间

从计算机内存复制到图形处理器GPU的存储空间

n=1,...,N_a,N_a为雷达***方位向采样点数，从计算机内存复制到图形处理器GPU的存储空间

a=1,...,s_a,s_a为场景方位向总的散射点数，r=1,...,s_r，s_r为场景距离向总的散射点数，从计算机内存复制到图形处理器GPU的存储空间

中；

步骤4：为图形处理器GPU内核函数划分线程网格

首先确定线程块中的线程数，方法是：Db_x为线程块第一维分配的线程数，Db_x设置为8的整数倍，同时不超过GPU线程硬件规定的范围；Db_y为线程块第二维分配的线程数，Db_y设置为8的整数倍，同时不超过GPU线程硬件规定的范围；

然后确定线程块的数目，方法是：Dg_x为线程网格第一维分配的线程块数目，Dg_x由场景方位向总的散射点数s_a和线程块第一维分配的线程数Db_x共同确定，计算公式为Dg_x=floor((s_a+Db_x-1)/Db_x)；Dg_y为线程网格第二维分配的线程块数目，Dg_y由场景距离向总的散射点数s_r和线程块第二维分配的线程数Db_y共同确定，计算公式为Dg_y=floor((s_r+Db_y-1)/Db_y)，其中floor为向下取整函数；

步骤5：平台主天线距离压缩后数据沿方位向成像

图形处理器GPU实现平台主天线距离压缩后数据

沿方位向成像的具体步骤是：

步骤5.1：确定场景散射点与线程网格中线程的对应关系，方法是：对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，blockIdx.x为线程网格中线程块的第一维索引，blockIdx.y为线程网格中线程块的第二维索引，threadIdx.x为线程块中线程的第一维索引，threadIdx.y为线程块中线程的第二维索引，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量，散射点P(a,r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockIdx.x和线程块中线程的第一维索引threadIdx.x的对应关系为a=threadIdx.x+blockIdx.x·Db_x+1，散射点P(a,r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockIdx.y和线程块中线程的第二维索引threadIdx.y的对应关系为r=threadIdx.y+blockIdx.y·Db_y+1，其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数；

步骤5.2：l_m为主天线半个合成孔径内的PRF时刻个数，计算公式为l_m=round(L_m/d_a/2)，其中round为近似取整函数，L_m为主天线合成孔径长度，d_a为场景方位向散射点间隔，L_m和d_a均由步骤1初始化得到；i_m为主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，i_m的取值满足条件：i_m=a-l_m,...，a+l_m，若a-l_m<1，则a-l_m=1，若a+l_m>N_a，则a+l_m=N_a，其中N_a为雷达***方位向采样点数，即i_m只取雷达***采样范围内的主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻；采用公式

计算得到雷达平台主天线i_m时刻处散射点P(a,r)的距离史

其中

为L2范数，

为步骤1中得到的雷达平台主天线第i_m个PRF时刻的天线相位中心矢量，为步骤1中得到的场景散射点P(a,r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；

步骤5.3：

为场景散射点P(a,r)的距离史对应的距离门，计算公式为

{\overset{&OverBar;}{I}}_{m} (i_{m}; a, r) = I_{c} + round (({\overset{&OverBar;}{R}}_{m} (i_{m}; a, r) - R_{mc}) / d_{r}),

其中round为近似取整函数，

为步骤5.2中得到的雷达平台主天线i_m时刻处散射点P(a,r)的距离史，R_mc为场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，I_c为场景参考点在主天线回波数据中的距离门，d_r为场景距离向散射点间隔，R_mc、I_c和d_r均由步骤1初始化得到，i_m为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；从步骤2中得到的平台主天线距离压缩后数据的第i_m行中取第

其中W₀为标准辛格插值的半个窗长；

步骤5.4：根据补偿相位因子的计算公式

{\overset{&OverBar;}{F}}_{m} (i_{m}; a, r) = \exp {j \cdot 4 \cdot π \cdot {\overset{&OverBar;}{R}}_{m} (i_{m}; a, r) / λ},

得到散射点P(a,r)在PRF时刻i_m应补偿的相位因子

其中j为虚数单位（即-1开根），

为步骤5.2中得到的雷达平台主天线i_m时刻处散射点P(a,r)的距离史，λ为雷达***工作的信号波长，λ由步骤1初始化得到，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；将步骤5.3中得到的插值重采样后的数据

与散射点P(a,r)在PRF时刻i_m应补偿的相位因子相乘，得到相位补偿后的数据

相加，得到雷达平台主天线场景散射点P(a,r)成像后的数据

即

其中i_m为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻；

a=1,...,s_a,s_a为场景方位向总的散射点数，r=1,...,s_r，s_r为场景距离向总的散射点数，从图形处理器GPU的

复制到计算机内存中，其中

为步骤3中在图形处理器GPU的全局存储器上分配的存储空间；

步骤6：平台副天线距离压缩后数据沿方位向成像

步骤6.1：确定场景散射点与线程网格中线程的对应关系，方法是：对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，blockIdx.x为线程网格中线程块的第一维索引，blockIdx.y为线程网格中线程块的第二维索引，threadIdx.x为线程块中线程的第一维索引，threadIdx.y为线程块中线程的第二维索引，blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量，散射点P(a,r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockIdx.x和线程块中线程的第一维索引threadIdx.x的对应关系为a=threadIdx.x+blockIdx.x·Db_x+1，散射点P(a,r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockIdx.y和线程块中线程的第二维索引threadIdx.y的对应关系为r=threadIdx.y+blockIdx.y·Db_y+1，其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数；

步骤6.2：l_s为副天线半个合成孔径内的PRF时刻个数，计算公式为l_s=round(L_s/d_a/2)，其中round为近似取整函数，L_s为副天线合成孔径长度，d_a为场景方位向散射点间隔，L_s和d_a均由步骤1初始化得到；i_s为副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，i_s的取值满足条件：i_s=a-l_s,...,a+l_s，若a-l_s<1，则a-l_s=1，若a+l_s>N_a，则a+l_s=N_a，其中N_a为雷达***方位向采样点数，即i_s只取雷达***采样范围内的副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻；采用公式

计算得到雷达平台副天线i_s时刻处散射点P(a,r)的距离史采用公式

计算得到雷达平台主天线i_s时刻处散射点P(a,r)的距离史其中

为L2范数，

为步骤1中得到的场景散射点P(a,r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；

步骤6.3：

为场景散射点P(a,r)的距离史

对应的距离门，计算公式为

{\overset{&OverBar;}{I}}_{s} (i_{s}; a, r) = I_{c} + round (({\overset{&OverBar;}{R}}_{s} (i_{s}; a, r) - R_{sc}) / d_{r}),

其中round为近似取整函数，

为步骤6.2中得到的雷达平台副天线i_s时刻处散射点P(a,r)的距离史，R_sc为场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，I_c为场景参考点在主天线回波数据中的距离门，d_r为场景距离向散射点间隔，R_sc、I_c和d_r均由步骤1初始化得到，i_s为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；从步骤2中得到的平台副天线距离压缩后数据的第i_s行中取第

其中W₀为标准辛格插值的半个窗长；

步骤6.4：根据补偿相位因子的计算公式

{\overset{&OverBar;}{F}}_{s} (i_{s}; a, r) = \exp {j \cdot 4 \cdot π \cdot {\overset{&OverBar;}{R}}_{m} (i_{s}; a, r) / λ},

得到散射点P(a,r)在PRF时刻i_s应补偿的相位因子

其中j为虚数单位（即-1开根），

为步骤6.2中得到的雷达平台主天线i_s时刻处散射点P(a,r)的距离史，λ为雷达***工作的信号波长，λ由步骤1初始化得到，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；将步骤6.3中得到的插值重采样后的数据