专利名称：一种基于后向投影InSAR成像配准的GPU实现方法
技术领域：
本技术发明属于雷达技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。
背景技术：
合成孔径雷达(SAR)是一种高分辨率的微波成像系统。合成孔径雷达依靠雷达和目标之间的相对运动来合成虚拟阵列以获得方位向高分辨率，利用大带宽信号实现距离向高分辨率，可以对照射场景进行二维成像。合成孔径雷达干涉测量(InSAR)是一般SAR功能的延伸和扩展，利用两个或者多个位置不同的天线观测同一个目标场景，根据目标到不同天线的斜距差获得测量数据的干涉相位，再通过平台与地面观测场景的几何关系反演出地面场景的数字高程信息的技术。由于具有全天时、全天候的特点，InSAR已经成为当前提取大面积地表三维图像和地形高程变化信息的一项重要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。高质量的干涉相位是InSAR获取高精度地形数字高程模型(DEM)的基础。如何获取高质量干涉相位是InSAR高程反演的关键技术。在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 (Compute Unified DeviceArchitecture,计算统一设备架构)则是由NVIDIA公式于2006年提出的利用GPU实现通用计算的类C编程模型，其自身带有特定的函数库(如空间分配、数据传输函数、CUFFT函数库等)，见文献 NVIDIA. CUDA CUFFT library PG-05327-040_V01, February, 2011 和 NVIDIACUDA C Programming Guide Version 4.03/6/2011。在InSAR成像中，由于不同天线观测同一个目标场景，不同天线具有相同的成像空间，因此可以利用后向投影技术将不同天线的测量数据投影到同一个成像空间进行成像，从而在数据成像的同时实现不同天线SAR图像的精确配准。由于后向投影算法是逐点累加成像，不同点之间的成像处理是可以并行的，所以利用CUDA语言以及基于GPU技术来实现算法的并行化，可以大大提高后向投影算法的运算效率。
发明内容
为了提高InSAR成像数据的相位保持性能，实现高质量的InSAR干涉相位快速生成，并减少不同天线图像配准过程对InSAR高程反演精度的影响，本发明提出了一种基于合成孔径雷达后向投影算法的InSAR成像配准一体化的GPU实现方法，即一种基于后向投 同一个成像空间中进行成像，在成像过程中同时完成了 InSAR不同天线的高保相成像和图像配准，实现了 InSAR成像配准一体化，并结合GPU并行化处理技术实现后向投影算法的快速处理。为了方便描述本发明的内容，首先作以下术语定义定义I、数字高程模型(DEM)数字高程模型(Digital Elevation Model, DEM),是指利用一组有序数值阵列形式表示地表或地面高程的一种实体地面模型。本发明中DEM表示成一系列地面点的平面坐标X、Y和高程坐标Z组成的数据阵列。对于一个地面区域D，地形DEM表示为DEM=(Di) (xi； Yi, Zi), i e D}其中(Xi，Yi)是第i个地面像素点对应的平面坐标，Zi是对应的高程坐标。定义2、合成孔径雷达干涉测量(InSAR)合成孔径雷达干涉测量(InSAR)指利用两个或者两个以上的SAR数据中的相位信息进行相干处理，结合雷达参数和雷达几何位置信息反演地表三维及其变化信息的遥感技术，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。定义3、雷达成像空间雷达成像空间是指将场景空间中的散射点投影到距离向-方位向的二维空间坐标系，该空间由合成孔径雷达成像空间中的两个相互正交的坐标基确定。目前典型合成孔径雷达的成像空间包括距离向-方位向投影空间。本发明中用以下数学关系表示成像空间M M = {Ρ(ν,ιι) I P(v，w) = v· . +u.hveR其中孓和&表示构成成像空间M的相互正交的坐标基，分别表示距离向和方位向。P(v，w)为成像空间中的待观测点向量，U，V分别表示该点的距离和方位坐标，R表示实数。定义4、距离史与距离门距离史是指不同天线相位中心到场景中散射点的距离组成的序列。距离门是指对应距离史的回波数据在整个回波数据中的位置，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。定义5、合成孔径雷达成像场景参考点合成孔径雷达成像场景参考点是指合成孔径雷达成像空间中的某个散射点，作为分析和处理场景中其他散射点的参照。定义6、图像配准图像配准是指将不同时间、不同传感器(成像设备)或不同条件下(天候、照度、摄像位置和角度等)对同一个地区获取的两幅或多幅图像进行地理坐标的匹配，使不同图像中的相同目标或者特征点位于图像中同一个位置的过程。本发明中的图像配准特指将不同天线SAR图像空间中对应的相同目标进行匹配的过程。定义7、合成孔径雷达后向投影算法 后向投影算法是基于匹配滤波原理的合成孔径雷达成像算法，其主要通过相干累加实现合成孔径雷达数据的聚焦成像。详细内容可参考文献“Research on A novel fastback projection algorithm for strip map bistatic SAR imaging”，Huang Yulin 等。定义8、范数设X是数域C上线性空间，称Il · Il为X上的范数(norm)，若它满足1.正定性Il X Il 彡 0，且 Il X Il =0〈=>X=0 ；2.齐次性Il aX Il =Ia Il X Il ；3.次可加性(三角不等式)I|X+Y< Il X Il+ 11 Y II。若父二^^…七^为狀丨维离散信号’向量父⑶范数为
= ΣΜΓ，L1 范数为 14=ΣΗ，L2 范数为I ^II =i ZW2 T
V/=I J =ιV /=1 J定义9、合成孔径雷达标准距离压缩方法合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射参数，采用以下公式生成参考信号，并采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。/(/) = expO-a-~·/2) ε[-γ,γ]其中，j为虚数单位(即-I开根)，f(t)为距离压缩参考函数，B为雷达发射基带信号的信号带宽，I；为雷达发射信号脉冲宽度，t为时间变量，取值范围从到详见文献
“雷达成像技术”，保铮等编著，电子工业出版社出版。定义10、合成孔径雷达原始回波仿真方法合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定系统参数条件下具有SAR回波信号特性的原始信号的方法，详细内容可参考文献“InSAR回波信号与系统仿真研究”，张剑琦，哈尔滨工业大学硕士论文。定义11、天线相位中心天线相位中心是指雷达天线向外辐射信号的中心，本发明中天线相位中心指雷达天线的轨迹位置。定义12、辛格插值(sine插值)方法辛格插值方法是指对于一个带限信号，在满足采样定理的情况下，采用卷积核为sine的函数h(x)，h(x)的长度即窗长为W。「 π “、 ，、 sin(^x)/ (λ. j = sin c (χ)=-----
πχ进行对已离散的信号gd(i)插值，得到插值后所要的信号sr(x)=Z^(/)sinc(x-0
i详见文献“合成孔径雷达成像算法与实现”，Frank H. Wong等编著，电子工业出版社出版。
定义13、取整函数取整函数是指将任一实数映射到相近的整数的一类函数，本发明中用到三个取整函数，分别为ceil、floor和round。ceil(x)为向上取整函数，表示取不超过X的最大一个整数。floor (x)为向下取整函数,表示取不小于x的最小一个整数。round(x)为近似取整函数,表示按四舍五入原则对x取整。定义14、合成孔径与PRF时刻合成孔径是指对于测绘场景中的一个散射点从进入雷达波束照射范围至离开雷达波束照射范围的这段时间内，雷达波束中心所走过的长度。雷达平台飞过一个合成孔径所需要的时间称为慢时间，雷达系统以一定的重复周期I；发射接收脉冲，因此慢时间可以表示为一个以重复周期I；为步长的离散化时间变量，其中每一个离散时间变量值为一个PRF时刻。详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。定义15、图形处理器(GPU)与计算统一设备架构(CUDA)图形处理器GPlXGraphic 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)是内核函数执行并行计算的第二层并行单元，同一线程块中的线程可以共享数据，没有执行顺序。详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。定义17、线程块索引与线程索引线程块索引是指每个线程块在线程网格中的位置，是CUDA的内建变量，线程块的第一维索引为blockldx. X，线程块的第二维索引为blockldx. y，线程块的第三维索引恒为Io线程索引是指每个线程在线程块中的位置，是CUDA的内建变量，线程的第一维索引为threadldx. X,线程的第二维索引为threadldx. y,线程的第三维索引为threadldx. z。详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。定义18、GPU全局存储器GPU全局存储器(global memory)是GPU片外的存储器，可读写，作用范围是整个程序，整个线程网格中的任意线程都能读写全局存储器的任意位置，存储空间大，生命周期为整个程序。详见文献“GPU高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。
定义19、GPU内存分配方法和GPU数据拷贝方法GPU内存分配方法采用CUDA函数库中的函数cudaMallocO实现在全局存储器中分配存储空间。GPU数据拷贝方法采用CUDA函数库中的函数cudaMemcpyO实现计算机内存与GPU全局存储器之间的数据复制。详见文献NVIDIA CUDA C Programming Guide Version 4· 03/6/2011 和 “GPU 高性能运算之CUDA”，张舒等编著，中国水利水电出版社出版。本发明提供一种基于后向投影InSAR成像配准的GPU实现方法，该方法的实现步骤如下步骤I :初始化InSAR成像系统参数InSAR成像空间由InSAR成像空间中的两个相互正交的坐标基确定，定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基，记做孓，该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基'^垂直的单位向量作为InSAR成像空间的第二个坐标基，记做瓦，该坐标基方向为距离向。InSAR雷达平台包含两组天线，即主天线和副天线，两组天线之间的距离为基线长，记做B1,主天线发射脉冲信号，经过Td时间的延迟，主天线和副天线同时接收回波延迟信号。雷达平台主天线接收的回波数据，记做I，雷达平台副天线接收的回波数据，记做民，其中Ew和良均为二维矩阵，第一维均对应方位向，第二维均对应距离向，即二维矩阵S,,,和艮的行存储的是方位向数据，二维矩阵瓦和瓦的列存储的是距离向数据。初始化InSAR成像系统参数包括雷达系统工作的信号波长，记做λ，雷达平台主天线发射信号带宽，记做B，雷达平台主天线发射脉冲时宽，记做 ；，雷达平台接收系统采样频率，记做Fs，雷达系统脉冲重复频率，记做prf，雷达平台主天线合成孔径长度，记做Lm，雷达平台副天线合成孔径长度，记做Ls，雷达平台速度矢量，记做 ，雷达平台主天线初始位置矢量，记做1,(0),雷达平台副天线初始位置矢量，记做艮(0)，场景参考点位置矢量，记做瓦,雷达系统距离向采样点数，记做队，雷达系统方位向采样点数，记做Na，场景距离向散射点间隔，记做4，场景方位向散射点间隔，记做da，场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，记做Rm。，场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，记做Rs。，场景参考点在主天线回波数据和副天线回波数据中的距离门相同，距离门位置记做I。，线程网格中线程块的第一维索引，记做blockldx. X，线程网格中线程块的第二维索引，记做blockldx. y,每个线程块中线程的第一维索引，记做threadldx. x,每个线程块中线程的第二维索引，记做 threadldx. y, blockldx. x、blockIdx.Y、threadldx. x和threadldx. y均为计算统一设备架构CUDA的内建变量。采用公式P ( ) = Pw(O)+ . // /计算得到雷达平台主天线第η个PRF时刻的天线相位中心矢量5 ( )，采用公式瓦(η) =Ρι(0)+ ·η//计算得到雷达平台副天线第η个PRF时刻的天线相位中心矢量h(n)，n = Na，η表示第η个PRF时亥Ij,主天线所有Nf^PRF时刻的天线相位中心矢量占用的存储空间大小,记做，副天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小，记做Fs —size。采用公式F( ，r) = I + (r -1) · d + ( -1) · H计算得到场景散射点P (a，r)的位置矢量，r表示散射点位于场景距离向的第！■个位置，r=l, . . . , sr, sr为场景距离向总的散射点数，a表示散射点位于场景方位向的第a个位置，a=l，...，sa，sa为场景方位向总的散射点数，场景所有散射点的位置矢量占用的存储空间大小记做。B (o，r)为雷达平台主天线场景散射点P(a，r)成像后的数据,1,,,(0^)初始化为O, a=l, . . . , sa, sa为场景方位向总的散射点数，r=l, . . . , sr, Sr为场景距离向总的散射点数，雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小记做 ，B (ο,r)为雷达平台副天线场景散射点P (a, r)成像后的数据,5s(a，r)初始化为O,a=l, . . .，sa, sa为场景方位向总的散射点数，r=l, . . . , sr, Sr为场景距离向总的散射点数,雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小记做民。步骤2 =InSAR原始回波数据进行距离压缩采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据I进行压缩，得到平台主天线距离压缩后数据，记做, I的数据大小记做§ _SiZe ο采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据
I,进行压缩，得到平台副天线距离压缩后数据，记做民，民的数据大小记做民—俞e。步骤3 :为图形处理器GPU分配数据采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为Sm 的存储空间，记做民,_D :采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为的存储空间，记做;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为E,_ ze的存储空间，记做采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为E_s/Ze的存储空间，记做;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为《Te的存储空间，记做β ;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为的存储空间，记做Sjb —D ;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为民的存储空间，记做Si _D ο采用传统的GPU数据拷贝方法，将步骤2中得到的平台主天线距离压缩后数据豆 从计算机内存复制到图形处理器GPU的存储空间iw_D中；采用传统的GPU数据拷贝方法，将步骤2中得到的平台副天线距离压缩后数据艮从计算机内存复制到图形处理器GPU的存储空间艮_^>中；采用传统的GPU数据拷贝方法，将步骤I中得到的主天线所有Na个PRF时亥Ij的天线相位中心矢量，n=l，. . . , Na, Na为雷达系统方位向采样点数，从计算机内存复制到图形处理器GPU的存储空间瓦_0中；采用传统的GPU数据拷贝方法，将步骤I中得到的副天线所有Na个PRF时刻的天线相位中心矢量瓦00，n = I..., N1，Na为雷达系统方位向采样点数，从计算机内存复制到图形处理器GPU的存储空间瓦中；采用传统的GPU数据拷贝方法，将步骤I中得到的场景所有散射点的位置矢量^ = ，sa为场景方位向总的散射点数，r=l,. . . , sr, sr为场景距离向总的散射点数，从计算机内存复制到图形处理器GPU的存储空间歹—D中。步骤4 :为图形处理器GPU内核函数划分线程网格线程网格中的线程块以二维形式划分，线程块中的线程也以二维形式划分，线程网格第一维总的线程数等于场景方位向总的散射点数Sa，线程网格第二维总的线程数等于场景距离向总的散射点数S，即一个场景散射点对应一个线程，具体划分步骤是首先确定线程块中的线程数，方法是Db_x为线程块第一维分配的线程数，Db_x设置为8的整数倍，同时不超过GPU线程硬件规定的范围；Db_y为线程块第二维分配的线程数，Db_y设置为8的整数倍，同时不超过GPU线程硬件规定的范围。
然后确定线程块的数目，方法是Dg_x为线程网格第一维分配的线程块数目，Dg_X由场景方位向总的散射点数sa和线程块第一维分配的线程数Db_x*同确定，计算公式为Dg_x=f Ioor ((sr+Db_x-l) /Db_x) ；Dg_y为线程网格第二维分配的线程块数目，Dg_y由场景距离向总的散射点数\和线程块第二维分配的线程数Db_y*同确定，计算公式为Dg_y=fIoor ((sa+Db_y-l) /Db_y),其中 floor 为向下取整函数。步骤5 :平台主天线距离压缩后数据沿方位向成像图形处理器GPU实现平台主天线距离压缩后数据Ia沿方位向成像的具体步骤是步骤5. I :确定场景散射点与线程网格中线程的对应关系，方法是对于场景散射点P(a，r), a表示散射点位于场景方位向的第a个位置，a=l, . . .，sa, saS场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数，blockldx. X为线程网格中线程块的第一维索引，blockldx. y为线程网格中线程块的第二维索引，threadldx. X为线程块中线程的第一维索引，threadldx. y为线程块中线程的第二维索引，blockldx. X、blockldx. y、threadldx. x 和 threadldx. y 均为计算统一设备架构CUDA的内建变量，散射点P(a，r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockldx. X和线程块中线程的第一维索引threadldx. x的对应关系为a=threadldx. x+blockldx. x *Db_x+l,散射点P (a, r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockldx. y和线程块中线程的第二维索引threadldx. y的对应关系为r=threadldx. y+blockldx. y · Db_y+1,其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数。步骤5. 2 :lm为主天线半个合成孔径内的PRF时刻个数,计算公式为lm=round(Lm/da/2)，其中ixnmd为近似取整函数，Lm为主天线合成孔径长度，4为场景方位向散射点间隔，Lffl和da均由步骤I初始化得到；im为主天线距离散射点P(a，r)前后半个合成孔径内的 PRF 时刻，im 的取值满足条件im = a_lm，, a+lm，若 a_lm〈l，则 a_lm = 1，若 a+lm>Na，则a+lm=Na，其中Na为雷达系统方位向采样点数，即1只取雷达系统采样范围内的主天线距离散射点P (a, r)前后半个合成孔径内的PRF时刻。采用公式=
算得到雷达平台主天线乜时刻处散射点P(a，r)的距离史其中Il □ ||2为1^2范数，？,4/m)为步骤I中得到的雷达平台主天线第imfPRF时刻的天线相位中心矢量，？fer)为步骤I中得到的场景散射点P(a，r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=l,. . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数。步骤5. 3为场景散射点P (a, r)的距离史对应的距离门，计算公式为 = ，其中 round 为近似取整函数,为步骤
5.2中得到的雷达平台主天线im时刻处散射点P (a，r)的距离史，Rfflc为场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，I。为场景参考点在主天线回波数据中的距离门，dr为场景距离向散射点间隔，Rmc、I。和dr均由步骤I初始化得到，iffl为步骤5. 2中得到的雷达平台主天线距离散射点P (a，r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l，...,sr,sr为场景距离向总的散射点数。从步骤2中得到的平台主天线距离压缩后数据的第^行中取第^^^个数据的前…个数据和后…个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据Ar)，其中Wtl为标准辛格插值的半个窗长。 步骤5.4:根据补偿相位因子的计算公式瓦(4;0^)=找|>{_/.41^(4);￥)/4，得到散射点P (a，r)在PRF时刻im应补偿的相位因子，其中j为虚数单位(即-I开根),瓦,也;为步骤5. 2中得到的雷达平台主天线im时刻处散射点P (a，r)的距离史，λ为雷达系统工作的信号波长，λ由步骤I初始化得到，a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l，...,sr,sr为场景距离向总的散射点数。将步骤5. 3中得到的插值重采样后的数据与散射点P(a，r)在PRF时刻im应补偿的相位因子相乘，得到相位补偿后的数据I。把雷达平台主天线距离散射点P (a，r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加，得到雷达平台主天线场景散射点P(a，r)
成像后的数据s=其中im为步骤5. 2中得到的雷达平台主天
°m Ka^r)，hf}
线距离散射点P (a, r)前后半个合成孔径内的PRF时刻。步骤5. 5 :采用传统的GPU数据拷贝方法，将步骤5. 4中得到的雷达平台主天线场景所有散射点成像后的数据民,(^，a = l，...,sa，sa为场景方位向总的散射点数，r=l,. . . , sr, sr为场景距离向总的散射点数，从图形处理器GPU的复制到计算机内存中，其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。步骤6 :平台副天线距离压缩后数据沿方位向成像图形处理器GPU实现平台副天线距离压缩后数据民沿方位向成像的具体步骤是步骤6. I :确定场景散射点与线程网格中线程的对应关系，方法是对于场景散射点P(a，r), a表示散射点位于场景方位向的第a个位置，a=l, . . .，sa, saS场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数，blockldx. X为线程网格中线程块的第一维索引，blockldx. y为线程网格中线程块的第二维索引，threadldx. X为线程块中线程的第一维索引，threadldx. y为线程块中线程的第二维索引，blockldx. X、blockldx. y、threadldx. x 和 threadldx. y 均为计算统一设备架构CUDA的内建变量，散射点P(a，r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockldx. X和线程块中线程的第一维索引threadldx. x的对应关系为a=threadldx. x+blockldx. x *Db_x+l,散射点P (a, r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockldx. y和线程块中线程的第二维索引threadldx. y的对应关系为r=threadldx. y+blockldx. y · Db_y+1,其中Db_x为线程块第一维分配的线程数,Db_y为线程块第二维分配的线程数。步骤6. 2 :ls为副天线半个合成孔径内的PRF时刻个数,计算公式为ls=round(Ls/da/2)，其中round为近似取整函数，Ls为副天线合成孔径长度，da为场景方位向散射点间隔，匕和4均由步骤I初始化得到；is为副天线距离散射点P(a，r)前后半个合成孔径内的 PRF 时刻，is 的取值满足条件is=a-ls, · · ·，a+ls，若 a_ls〈l，则 a_ls=l，若 a+ls>Na，则a+ls=Na，其中Na为雷达系统方位向采样点数，即is只取雷达系统采样范围内的副天线距离散射点P (a，r)前后半个合成孔径内的PRF时刻。采用公式= 计算得到雷达平台副天线is时刻处散射点P(a，r)的距离史瓦(^r)，采用公式R" (I ; ，〃) = ||ρ (I) - ρ(α，〃)||2计算得到雷达平台主天线i s时刻处散射点P (a, r)的距离史R (K^,r)，其中Il □ Il 2为L2范数，E(g为步骤I中得到的雷达平台副天线第Is个PRF时刻的天线相位中心矢量,瓦(υ为步骤I中得到的雷达平台主天线第isfPRF时刻的天线相位中心矢量，为步骤I中得到的场景散射点P(a，r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数。步骤6. 3为场景散射点P(a, r)的距离史见对应的距离门，计算公式%\(is\a,r)=Ic +roundaR^i/.aj-y-R^/d,),其中 round 为近似取整函数,iM 's:a，r)为步骤 6. 2 中得到的雷达平台副天线is时刻处散射点P (a，r)的距离史，Rsc为场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，I。为场景参考点在主天线回波数据中的距离门，dr为场景距离向散射点间隔，Rsc、I。和dr均由步骤I初始化得到，is为步骤6. 2中得至IJ的雷达平台副天线距离散射点P (a，r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l，. . . , sr, sr为场景距离向总的散射点数。从步骤2中得到的平台副天线距离压缩后数据民的第^行中取第个数据的前Wtl个数据和后Wtl个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据AAr)，其中Wtl为标准辛格插值的半个窗长。步骤6. 4 :根据补偿相位因子的计算公式旯(4;0^) = 6\卩卜4‘[豆上;0，0/2}，得到散射点P (a，r)在PRF时刻is应补偿的相位因子(4; ,O，其中j为虚数单位(即_1开根)，瓦,(Uir)为步骤6. 2中得到的雷达平台主天线is时刻处散射点P(a，r)的距离史，λ为雷达系统工作的信号波长，λ由步骤I初始化得到，a表示散射点位于场景方位向的第a个位置，a=l,. . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l，...,sr,sr为场景距离向总的散射点数。将步骤6. 3中得到的插值重采样后的数据亡#;^与散射点P(a，r)在PRF时刻is应补偿的相位因子相乘，得到相位补偿后的数据 、(dr)。把雷达平台副天线距离散射点P (a，r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加，得到雷达平台副天线场景散射点P(a，r)成像后
的数据1>力=’其中is为步骤6. 2中得到的雷达平台副天线距离散
射点P (a, r)前后半个合成孔径内的PRF时刻。步骤6.5 :采用传统的GPU数据拷贝方法，将步骤6. 4中得到的雷达平台副天线场景所有散射点成像后的数据B α = 1，…,5a，Sa为场景方位向总的散射点数，r=l,. . . , sr, sr为场景距离向总的散射点数，从图形处理器GPU的R—D复制到计算机内存中，其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。步骤7 :将步骤5中得到的雷达平台主天线场景所有散射点成像后的数据δ ,⑷.r)和步骤6中得到的雷达平台副天线场景所有散射点成像后的数据5s( ,r)进行共轭相乘，得到InSAR成像的干涉相位i(o，r),完成InSAR数据成像配准一体化处理，其中a表示散射点位于 场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数。本发明的创新点在于利用后向投影算法，将InSAR不同天线测量数据投影到同一个成像空间进行成像，实现InSAR数据成像和图像配准一体化，并根据后向投影算法逐点相干积累原理，利用GPU并行化技术实现算法快速处理。通过仿真实验验证，相对于传统InSAR干涉相位生成，该算法具有保相性能好、成像配准一体化、并行化处理提高算法效率等优点，实现了高精度InSAR干涉相位的快速生成。本发明的优点在于在数据成像过程中同时完成InSAR高保相成像与配准，在数据成像时同时获得高质量干涉相位，不需要再进行图像配准处理，并且可以利用GPU并行化处理提高算法效率。本发明可以应用于合成孔径雷达成像，地球遥感等领域。
图I为本发明所提供方法的流程示意框图；图2为本发明具体实施方式
采用的干涉合成孔径雷达平台飞行几何关系及仿真圆锥场景图；其中，横坐标(X轴)为方位向，纵坐标(Y轴)为距离向，垂直坐标(Z轴)为高度向，B1为雷达平台主天线和副天线之间的基线，基线的两个端点分别为雷达平台主天线和副天线，P (a, r)为场景中任一散射点，a表示散射点位于场景方位向的第a个位置，a=l,512，r表示散射点位于场景距离向的第r个位置，r=l，512， 为雷达平台飞行速度矢量，艮(/, )为雷达平台主天线第im APRF时刻的天线相位中心矢量，^从)为雷达平台副天线第is个PRF时刻的天线相位中心矢量,为雷达平台主天线im时刻处散射点P (a, r)的距离史，瓦O; w)为雷达平台副天线is时刻处散射点P(a，r)的距离史，im为雷达平台主天线距离散射点P (a, r)半个合成孔径内的PRF时刻，is为雷达平台副天线距离散射点P (a, r)半个合成孔径内的PRF时刻，Sa为仿真场景方位向总的散射点数，Sr为仿真场景距离向总的散射点数，da为仿真场景方位向散射点间隔，dr为仿真场景距离向散射点间隔。图3本发明具体实施方式
采用的线程网格划分；其中，Block为线程网格中线程块的二维划分，Thread为第一维索引为O、第二维索引为15的线程块的二维线程划分。
具体实施例方式本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在VC++、MATLAB7. O、CUDA4. O和GeForce GTX 560Τ 显卡上验证正确。具体实施步骤如下步骤I :初始化InSAR成像系统参数选择与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基i=
，该坐标基方向为方位向；选择InSAR成像空间的第二个坐标基孓=[I O Oj，该坐标基方向为距离向。初始化InSAR成像系统参数包括雷达系统工作的信号波长λ = O. 03m,雷达平台主天线发射信号带宽B=150MHz,雷达平台主天线发射脉冲时宽T1=I μ s,雷达平台接收系统采样频率Fs=300MHz，雷达系统脉冲重复频率prf=500Hz，雷达系统距离向采样点数 Nr=1024,雷达系统方位向采样点数Na=2048，场景距离向散射点间隔d,=0. 5m,场景方位向散射点间隔4=0. 3m,雷达平台主天线合成孔径长度Lm = 150m,雷达平台副天线合成孔径长度Ls=150m，雷达平台速度矢量 =
，速度的单位为m/s，雷达平台主天线初始位置矢量系m(0) =
,雷达平台副天线初始位置矢量瓦(0) = [10 O 6000],场景参考位置矢量，记做瓦=[8000 O O],位置的单位为m，场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值Rme=IOOOOm,场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值Rs。= 9992，场景参考点在主天线回波数据和副天线回波数据中的距尚门Ic= I。仿真场景为一立体圆锥，仿真场景方位向总的散射点数sa=512，仿真场景距离向总的散射点数S,=512，圆锥最小高度为0m，最大高度为50m。采用传统的合成孔径雷达原始回波仿真方法，生成雷达平台主天线仿真回波数据K，采用传统的合成孔径雷达原始回波仿真方法，生成雷达平台副天线仿真回波数据民。采用公式Ρ ,(Η) = Ρ (0) + ν· /^/计算得到雷达平台主天线第η个PRF时刻的天线相位中心矢量5 ,(H)，采用公式？, (H) = E (O) + . H / prf计算得到雷达平台副天线第η个PRF时刻的天线相位中心矢量Pj(H)，n 2048，主天线所有2048个PRF时刻的天线相位中心矢量占用的存储空间大小Pm —Size = 3.2048.Kyies，副天线所有2048个PRF时刻的天线相位中心矢量占用的存储空间大小瓦—size = 3 ■ 204E-Shytes。采用公式P(a,r) = Pc +0-1)-0.5- +(α-1)·0.3· ；.计算得到场景散射点P (a, r)的位置矢量，r表不散射点位于场景距离向的第r个位置，r=l,. . .，512, a表不散射点位于场景方位向的第a个位置，a=l，...，512，场景所有散射点的位置矢量占用的存储空间大小Ρ_^ = 3·512·512· Wytes。艮紅r)为雷达平台主天线场景散射点P (a，r)成像后的数据，
初始化为O，a=l, 512，r=l,. . .，512，雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小5 _size = 512.512揚卿，B (W)为雷达平台副天线场景散射点P (a, r)成像后的数据，Ss(fl，r)初始化为O，a=l，512，r=l，. . .，512，雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小民_size = 512-512-Bbytes。步骤2 =InSAR原始回波数据进行距离压缩
采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据,进行压缩，得到平台主天线距离压缩后数据的大小_size = 1024-2048-Sfevtev。采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据氧进行压缩，得到平台副天线距离压缩后数据豆民的大小艮—1^ = 1024.2048.% ^步骤3 :为图形处理器GPU分配数据采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为s, _.ν/α的存储空间，记做-D;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为的存储空间，记做;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为5 , 的存储空间，记做瓦— :采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为的存储空间，记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为P_ iz的存储空间，记做5—β ;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小 为的存储空间，记做民,—￡>;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为民的存储空间，记做S, —D。采用传统的GPU数据拷贝方法，将步骤2中得到的平台主天线距离压缩后数据民t从计算机内存复制到图形处理器GPU的存储空间中；采用传统的GPU数据拷贝方法，将步骤2中得到的平台副天线距离压缩后数据瓦从计算机内存复制到图形处理器GPU的存储空间中；采用传统的GPU数据拷贝方法，将步骤I中得到的主天线所有2048个PRF时刻的天线相位中心矢量瓦》( )，n=l,. . .，2048，从计算机内存复制到图形处理器GPU的存储空间D中；采用传统的GPU数据拷贝方法，将步骤I中得到的副天线所有2048个PRF时刻的天线相位中心矢量瓦(》)，// = 1,...,2048 ,从计算机内存复制到图形处理器GPU的存储空间瓦―β中；采用传统的GPU数据拷贝方法，将步骤I中得到的场景所有散射点的位置矢量5(a,r), α = 1,...,512 , r = l，...，512 ,从计算机内存复制到图形处理器GPU的存储空间P_D 中。步骤4 :为图形处理器GPU内核函数划分线程网格线程块第一维分配的线程数Db_x=16，线程块第二维分配的线程数Db_y=16，线程网格第一维分配的线程块数目Dg_x=32，线程网格第二维分配的线程块数目Dg_y=32。步骤5 :平台主天线距离压缩后数据沿方位向成像图形处理器GPU实现平台主天线距离压缩后数据沿方位向成像的具体步骤是步骤5. I :确定场景散射点与线程网格中线程的对应关系，方法是对于场景散射点P(a，r), a表不散射点位于场景方位向的第a个位置，a=l, . . . , 512, r表不散射点位于场景距离向的第r个位置，r=l,512，blockldx. x为线程网格中线程块的第一维索引，blockldx. y为线程网格中线程块的第二维索引，threadldx. x为线程块中线程的第一维索引，threadldx. y 为线程块中线程的第二维索引，blockldx. x、blockldx. y、threadldx.X和threadldx. y均为计算统一设备架构CUDA的内建变量，散射点P(a，r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockldx. X和线程块中线程的第一维索引threadldx. x 的对应关系为 a=threadldx. x+blockldx. x · 16+1,散射点 P (a, r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockldx. y和线程块中线程的第二维索引 threadldx. y 的对应关系为 r=threadldx. y+blockldx. y · 16+1。步骤5. 2 : Im为主天线半个合成孔径内的PRF时刻个数，Iffl = 250 ； im为主天线距离散射点P (a, r)前后半个合成孔径内的PRF时刻，im的取值满足条件im=a_250，, a+250,若 a-250〈l，则 a-1^1，若 a+250>2048，则 a+1^2048。采用公式瓦,,0 况”)=|瓦,0;,,)-5( 计算得到雷达平台主天线乜时刻处散射点P(a，r)的距离史瓦其中Il □ ||2SL2范数，？.4U为步骤I中得到的雷达平台主天线第imfPRF时刻的天线相位中心矢量为步骤I中得到的场景散射点P(a，r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=l,. . .，512, r表示散射点位于场景距离向的第r个位置，r=l,. . .，512。步骤5. 3为场景散射点P (a, r)的距离史* (/ ,:a.r)对应的距离门，计算公式为;o,r) = l + TOMHi/((民,(ira;ij,r)-1000())/().5),其中 round 为近似取整函数,R,Ai ,;a.r)为步骤5. 2中得到的雷达平台主天线im时刻处散射点P (a，r)的距离史，im为步骤5. 2中得到的雷达平台主天线距离散射点P(a，r)前后半个合成孔径内的PRF时刻，a表示散射点位 于场景方位向的第a个位置，a=l,512，r表示散射点位于场景距离向的第r个位置，r=l，...，512。从步骤2中得到的平台主天线距离压缩后数据I的第乜行中取第I,
个数据的前4个数据和后4个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据A,,(i, ;a.r)。步骤5.4:根据补偿相位因子的计算公式Fm (/, ；a,r) = exp{j-4-π-Rm(/；,;a,r)/0.03},得到散射点 P (a, r)在 PRF 时刻 in 应补偿的相
位因子,其中j为虚数单位(即-I开根)，瓦dr)为步骤5. 2中得到的雷达平台主天线im时刻处散射点P(a，r)的距离史，a表示散射点位于场景方位向的第a个位置，a=l，...，512，r表示散射点位于场景距离向的第r个位置，r=l，. . .，512。将步骤5. 3中得到的插值重采样后的数据与散射点P(a，r)在PRF时刻im应补偿的相位因子相乘，得到相位补偿后的数据。把雷达平台主天线距离散射点P(a，r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据I； (4^，r)相加，得到雷达平台主
天线场景散射点P(a，r)成像后的数据s .即=’其中im为步骤5. 2
中得到的雷达平台主天线距离散射点P(a，r)前后半个合成孔径内的PRF时刻。步骤5. 5 :采用传统的GPU数据拷贝方法，将步骤5. 4中得到的雷达平台主天线场景所有散射点成像后的数据δ,,,(ο,Γ), a = 1,...,512 , r = 1,...,512 ,从图形处理器GPU的iw_￡ 复制到计算机内存中，其中i~_i 为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。步骤6 :平台副天线距离压缩后数据沿方位向成像图形处理器GPU实现平台副天线距离压缩后数据民沿方位向成像的具体步骤是步骤6. I :确定场景散射点与线程网格中线程的对应关系，方法是对于场景散射点P(a，r), a表不散射点位于场景方位向的第a个位置，a=l, . . . , 512, r表不散射点位于场景距离向的第r个位置，r=l,512，blockldx. x为线程网格中线程块的第一维索引，blockldx. y为线程网格中线程块的第二维索引，threadldx. x为线程块中线程的第一维索引，threadldx. y 为线程块中线程的第二维索引，blockldx. x、blockldx. y、threadldx.X和threadldx. y均为计算统一设备架构CUDA的内建变量，散射点P(a，r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockldx. X和线程块中线程的第一维索引threadldx. x 的对应关系为 a=threadldx. x+blockldx. x · 16+1,散射点 P (a, r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockldx. y和线程块中线程的第二维索引 threadldx. y 的对应关系为 r=threadldx. y+blockldx. y · 16+1。步骤6. 2 ls为副天线半个合成孔径内的PRF时刻个数，ls=250 ；is为副天线距离散射点P (a, r)前后半个合成孔径内的PRF时刻，is的取值满足条件is=a_250，a+250,若 a-250〈l，则 a-250=l，若 a+250>2048，则 a+250=2048。采用公式瓦ft;w) = ||瓦(O-计算得到雷达平台副天线is时刻处散射点P(a，r)的距离史瓦(/ ;ν)，采用公式i w0:;a,r) = | M(is)- (fl，r>l计算得到雷达平台主天线is时刻处散射点P(a，r)的距离史
其中Il □ Il 2SL2范数，瓦ft)为步骤I中得到的雷达平台副天线第is fPRF时刻的天线相位中心矢量，5 )力步骤I中得到的雷达平台主天线第^个PRF时刻的天线相位中心矢量为步骤I中得到的场景散射点P(a，r)的位置矢量，a表示散射点位 于场景方位向的第a个位置，a=l,512，r表示散射点位于场景距离向的第r个位置，r=l，…，512。步骤6. 3 :1>:;0!，/0为场景散射点？(&，r)的距离史瓦(z;;a,r)对应的距离门，计算公式 (/；;O,r) = I + round((Rt(/；a,r)-9992)/0.5)，其中 round 为近似取整函数,见( ::0,;·)为步骤 6. 2 中
得到的雷达平台副天线is时刻处散射点P (a，r)的距离史，is为步骤6. 2中得到的雷达平台副天线距离散射点P (a, r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=l,. . . , 512,r表不散射点位于场景距离向的第r个位置，r=l,. . . , 512。从步骤2中得到的平台副天线距离压缩后数据艮的第is行中取第个数据的前4个数据和后4个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据 c,,,(i,n;e，r) ο步骤6.4:根据补偿相位因子的计算公式艮(4;0^) = — [/.41反几;￥)/0.03}，得到散射点P(a，r)在PRF时刻is应补偿的相位因子Fs.，其中j为虚数单位(即_1开根为步骤6. 2中得到的雷达平台主天线is时刻处散射点P(a，r)的距离史，a表不散射点位于场景方位向的第a个位置，a=l, . . .，512, r表不散射点位于场景距离向的第r个位置，r=l,. . .，512。将步骤6. 3中得到的插值重采样后的数据与散射点P (a, r)在PRF时刻is应补偿的相位因子Fs 相乘，得到相位补偿后的数据冗
把雷达平台副天线距离散射点P(a，r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加，得到雷达平台副天线场景散射点P(a，r)成像后的数据5s(a，r),即
=其中is为步骤6. 2中得到的雷达平台副天线距离散射点P(a，r)前
h
后半个合成孔径内的PRF时刻。步骤6. 5 :采用传统的GPU数据拷贝方法，将步骤6. 4中得到的雷达平台副天线场景所有散射点成像后的数据Ss(w·)，a = 1,...,512, /^ = 1，...，512，从图形处理器6 "勺瓦—1)复制到计算机内存中，其中_D为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。
步骤7 :将步骤5中得到的雷达平台主天线场景所有散射点成像后的数据和步骤6中得到的雷达平台副天线场景所有散射点成像后的数据5s(a，r)进行共轭相乘，得到InSAR成像的干涉相位，完成InSAR数据成像配准一体化处理，其中a表示散射点位于场景方位向的第a个位置，a=l,. . .，512，r=l, . . .，512。通过本发明具体实施方式
的仿真及测试，本发明所提供的基于后向投影的InSAR成像配准一体化及GPU实现方法能够实现InSAR不同天线成像和图像配准同时完成，与现有的InSAR干涉相位生成方法相比，本发明克服现有方法数据成像和图像配准相互隔离的缺点，可直接从原始数据成像处理过程中获得配准后的干涉相位，减少数据处理流程并且 提高了成像相位保持精度，为InSAR高质量干涉相位生成提供了一种新方法。同时，利用GPU并行化技术实现了算法的快速处理。
权利要求
1.一种基于后向投影InSAR成像配准的GPU实现方法，其特征是它包括以下步骤 步骤I :初始化InSAR成像系统参数 InSAR成像空间由InSAR成像空间中的两个相互正交的坐标基确定，定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标棊,记做I，该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基 垂直的单位向量作为InSAR成像空间的第二个坐标基，记做孓，该坐标基方向为距离向； InSAR雷达平台包含两组天线，即主天线和副天线，两组天线之间的距离为基线长，记做B1，主天线发射脉冲信号，经过Td时间的延迟，主天线和副天线同时接收回波延迟信号；雷达平台主天线接收的回波数据，记做氣,，雷达平台副天线接收的回波数据，记做瓦，其中i, 和艮均为二维矩阵，第一维均对应方位向，第二维均对应距离向，即二维矩阵ijPlJ勺行存储的是方位向数据，二维矩阵1,,,和民的列存储的是距离向数据； 初始化InSAR成像系统参数包括雷达系统工作的信号波长，记做λ，雷达平台主天线发射信号带宽，记做B，雷达平台主天线发射脉冲时宽，记做 ；，雷达平台接收系统采样频率，记做Fs，雷达系统脉冲重复频率，记做prf，雷达平台主天线合成孔径长度，记做Lm，雷达平台副天线合成孔径长度，记做Ls，雷达平台速度矢量，记做ψ ,雷达平台主天线初始位置矢量，记做瓦(O)，雷达平台副天线初始位置矢量，记做瓦(O),场景参考点位置矢量，记做瓦,雷达系统距离向采样点数，记做队，雷达系统方位向采样点数，记做Na，场景距离向散射点间隔，记做4，场景方位向散射点间隔，记做da，场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，记做Rm。，场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，记做Rs。，场景参考点在主天线回波数据和副天线回波数据中的距离门相同，距离门位置记做I。，线程网格中线程块的第一维索引，记做blockldx. x,线程网格中线程块的第二维索引，记做blockldx. y,每个线程块中线程的第一维索引，记做threadldx. x,每个线程块中线程的第二维索引，记做 threadldx. y, blockldx. x、blockldx. y、threadldx. x和threadldx. y均为计算统一设备架构CUDA的内建变量； 采用公式P (η) = P, (0) + V-n/prf计算得到雷达平台主天线第η个PRF时刻的天线相位中心矢量 》( )，采用公式瓦(η)= 40)+ ·η/；μ/计算得到雷达平台副天线第η个PRF时刻的天线相位中心矢量瓦(《)，n = I,..., Na，η表示第η个PRF时刻,主天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小,记做—size，副天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小，记做瓦—size ;采用公式？(o，r) =瓦+(r-l)<- , +(α-1)·<- ；计算得到场景散射点P (a，r)的位置矢量，r表示散射点位于场景距离向的第！■个位置，r=l, . . . , sr, sr为场景距离向总的散射点数，a表示散射点位于场景方位向的第a个位置，a=l，...，sa，sa为场景方位向总的散射点数，场景所有散射点的位置矢量占用的存储空间大小记做5_.由￡;1,,,(^)为雷达平台主天线场景散射点？(&，1·)成像后的数据B (,,)初始化为O, a=l, . . . , sa, sa为场景方位向总的散射点数，r=l, . . . , sr, Sr为场景距离向总的散射点数，雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小记做B (V)为雷达平台副天线场景散射点P(a，r)成像后的数据，民初始化为0，a=l, . . .，sa, sa为场景方位向总的散射点数，r=l, . . . , sr, Sr为场景距离向总的散射点数,雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小记做民步骤2 : InSAR原始回波数据进行距离压缩 采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据I,,,进行压缩，得到平台主天线距离压缩后数据，记做，民^的数据大小记做； 采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据氧进行压缩，得到平台副天线距离压缩后数据，记做艮，艮的数据大小记做民JLl步骤3 :为图形处理器GPU分配数据 采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为I的存储空间，记做采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为Ss的存储空间，记做民_D ;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为瓦的存储空间，记做 采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大/J力？ —的存储空间，记做采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为？」細的存储空间，记做;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为的存储空间，记做Sm _D ;采用传统的GPU内存分配方法，在图形处理器GPU的全局存储器上分配大小为民的存储空间，记做Si _D ； 采用传统的GPU数据拷贝方法，将步骤2中得到的平台主天线距离压缩后数据5 从计算机内存复制到图形处理器GPU的存储空间中；采用传统的GPU数据拷贝方法，将步骤2中得到的平台副天线距离压缩后数据民从计算机内存复制到图形处理器GPU的存储空间艮_0中；采用传统的GPU数据拷贝方法，将步骤I中得到的主天线所有Na个PRF时刻的天线相位中心矢量，n=l,. . . , Na, Na为雷达系统方位向采样点数,从计算机内存复制到图形处理器GPU的存储空间 中；采用传统的GPU数据拷贝方法，将步骤I中得到的副天线所有NaAPRF时刻的天线相位中心矢量豕(》)， = 1,... N N1为雷达系统方位向采样点数，从计算机内存复制到图形处理器GPU的存储空间E 中；采用传统的GPU数据拷贝方法，将步骤I中得到的场景所有散射点的位置矢量巧^)，a = l,...5tse，sa为场景方位向总的散射点数，r=l,. . . , sr, sr为场景距离向总的散射点数，从计算机内存复制到图形处理器GPU的存储空间中； 步骤4 :为图形处理器GPU内核函数划分线程网格 线程网格中的线程块以二维形式划分，线程块中的线程也以二维形式划分，线程网格第一维总的线程数等于场景方位向总的散射点数Sa，线程网格第二维总的线程数等于场景距离向总的散射点数S，即一个场景散射点对应一个线程，具体划分步骤是 首先确定线程块中的线程数，方法是Db_x为线程块第一维分配的线程数，Db_x设置为8的整数倍，同时不超过GPU线程硬件规定的范围；Db_y为线程块第二维分配的线程数，Db_y设置为8的整数倍，同时不超过GPU线程硬件规定的范围； 然后确定线程块的数目，方法是Dg_x为线程网格第一维分配的线程块数目，Dg_x由场景方位向总的散射点数Sa和线程块第一维分配的线程数Db_x*同确定，计算公式为Dg_x=fIoor ((sr+Db_x-l)/Db_x) ;Dg_y为线程网格第二维分配的线程块数目，Dg_y由场景距离向总的散射点数\和线程块第二维分配的线程数Db_y共同确定，计算公式为Dg_y=fIoor ((sa+Db_y-l) /Db_y),其中 floor 为向下取整函数；步骤5 :平台主天线距离压缩后数据沿方位向成像 图形处理器GPU实现平台主天线距离压缩后数据沿方位向成像的具体步骤是 步骤5. I :确定场景散射点与线程网格中线程的对应关系，方法是对于场景散射点P (a, r), a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, Sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数，blockldx. X为线程网格中线程块的第一维索引，blockldx. y为线程网格中线程块的第二维索引，threadldx. X为线程块中线程的第一维索引，threadldx. y为线程块中线程的第二维索引，blockldx. X、blockldx. y、threadldx. x 和 threadldx. y 均为计算统一设备架构CUDA的内建变量，散射点P(a，r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockldx. X和线程块中线程的第一维索引threadldx. x的对应关系为a=threadldx. x+blockldx. x *Db_x+l,散射点P (a, r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockldx. y和线程块中线程的第二维索引threadldx. y的对应关系为r=threadldx. y+blockldx. y · Db_y+1,其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数； 步骤5. 2 lm为主天线半个合成孔径内的PRF时刻个数，计算公式为lm=round(Lm/da/2)，其中ixnmd为近似取整函数，Lm为主天线合成孔径长度，4为场景方位向散射点间隔，Ln^P 4均由步骤I初始化得到；im为主天线距离散射点P(a，r)前后半个合成孔径内的PRF 时刻，im 的取值满足条件im = a-lm，, a+lm，若 a-lm〈l，则 a_lm = 1，若 a+lm>Na，则a+lm=Na，其中Na为雷达系统方位向采样点数，S卩im只取雷达系统采样范围内的主天线距离散射点P(a，r)前后半个合成孔径内的PRF时刻；采用公式1(/,,,; ，「) = ||1(( )-再《，〃)||2计算得到雷达平台主天线im时刻处散射点P (a，r)的距离史I(^r)，其中Il □ ||2SL2范数，为步骤I中得到的雷达平台主天线第imfPRF时刻的天线相位中心矢量为步骤I中得到的场景散射点P (a，r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=l,. . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数； 步骤5. 3 ;fl，r)为场景散射点P (a, r)的距离史瓦 (dr)对应的距离门，计算公式为(L；= 4 + round((Rm(im;a,r)-Rnic)/dr)，其中 round 为近似取整函数,为步骤 5. 2 中得到的雷达平台主天线im时刻处散射点P (a，r)的距离史，Rfflc为场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值，I。为场景参考点在主天线回波数据中的距离 门，dr为场景距离向散射点间隔，Rmc、I。和dr均由步骤I初始化得到，iffl为步骤5. 2中得至IJ的雷达平台主天线距离散射点P (a，r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=l，sa，saS场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数；从步骤2中得到的平台主天线距离压缩后数据的第im行中取第T ,(/,,,; ，;·)个数据的前Wtl个数据和后Wtl个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据，其中Wtl为标准辛格插值的半个窗长； 步骤5.4:根据补偿相位因子的计算公式艮(4,凡『)^乂]3丨/.4.[11 (4; ^)/2}，得到散射点P (a, r)在PRF时刻im应补偿的相位因子艮，其中j为虚数单位卿-I开根)，瓦4dr)为步骤5. 2中得到的雷达平台主天线im时刻处散射点P(a，r)的距离史，λ为雷达系统工作的信号波长，λ由步骤I初始化得到，a表示散射点位于场景方位向的第a个位置，a=l,. . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数；将步骤5. 3中得到的插值重采样后的数据与散射点P (a, r)在PRF时刻im应补偿的相位因子見(4; V)相乘，得到相位补偿后的数据I;把雷达平台主天线距离散射点P (a，r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据t相加，得到雷达平台主天线场景散射点P(a，r)成像后的数据,ΒΡΒΜ(α,Γ) = ΣΤ (4；α,Γ) ’其中乜为步骤5. 2中得到的雷达平台主天线距离散射点P (a, r)前后半个合成孔径内的PRF时刻；步骤5. 5 :采用传统的GPU数据拷贝方法， 将步骤5. 4中得到的雷达平台主天线场景所有散射点成像后的数据民,( ·)，《 = 1，...，丨，sa为场景方位向总的散射点数，r=l,. . . , sr, sr为场景距离向总的散射点数，从图形处理器GPU的复制到计算机内存中，其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间； 步骤6 :平台副天线距离压缩后数据沿方位向成像 图形处理器GPU实现平台副天线距离压缩后数据民沿方位向成像的具体步骤是 步骤6. I :确定场景散射点与线程网格中线程的对应关系，方法是对于场景散射点P (a, r), a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, Sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数，blockldx. X为线程网格中线程块的第一维索引，blockldx. y为线程网格中线程块的第二维索引，threadldx. X为线程块中线程的第一维索引，threadldx. y为线程块中线程的第二维索引，blockldx. X、blockldx. y、threadldx. x 和 threadldx. y 均为计算统一设备架构CUDA的内建变量，散射点P(a，r)位于场景方位向的位置a与线程网格中线程块的第一维索引blockldx. X和线程块中线程的第一维索引threadldx. x的对应关系为a=threadldx. x+blockldx. x *Db_x+l,散射点P (a, r)位于场景距离向的位置r与线程网格中线程块的第二维索引blockldx. y和线程块中线程的第二维索引threadldx. y的对应关系为r=threadldx. y+blockldx. y · Db_y+1,其中Db_x为线程块第一维分配的线程数，Db_y为线程块第二维分配的线程数； 步骤6. 2 ls为副天线半个合成孔径内的PRF时刻个数，计算公式为ls=round(Ls/da/2)，其中round为近似取整函数，Ls为副天线合成孔径长度，da为场景方位向散射点间隔，匕和4均由步骤I初始化得到；is为副天线距离散射点P(a，r)前后半个合成孔径内的 PRF 时刻，is 的取值满足条件is=a-ls, · · ·，a+ls，若 a_ls〈l，则 a_ls=l，若 a+ls>Na，则a+ls=Na，其中Na为雷达系统方位向采样点数，即is只取雷达系统采样范围内的副天线距离散射点P(a，r)前后半个合成孔径内的PRF时刻；采用公式I(/s;rt，r) =计算得到雷达平台副天线is时刻处散射点P(a，r)的距离史瓦(^r)，采用公式R'"(I;α，〃) = ψ' (I)-Pfer)||2计算得到雷达平台主天线is时刻处散射点P (a, r)的距离史瓦其中Il □ ll2SL2范数，为步骤I中得到的雷达平台副天线第isfPRF时刻的天线相位中心矢量，为步骤I中得到的雷达平台主天线第isfPRF时刻的天线相位中心矢量，5(a，r)为步骤I中得到的场景散射点P(a，r)的位置矢量，a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数；步骤6. 3为场景散射点P (a, r)的距离史/对应的距离门，计算公式为Ii(Zi;a,r) =Ic+ round((Rj(is;a,r)-Rsc)/dT)，其中 round 为近似取整函数,见((:a,r)为步骤 6. 2 中得到的雷达平台副天线is时刻处散射点P (a，r)的距离史，Rsc为场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值，I。为场景参考点在主天线回波数据中的距离门，dr为场景距离向散射点间隔，Rs。、I。和火均由步骤I初始化得到，is为步骤6. 2中得到的雷达平台副天线距离散射点P (a，r)前后半个合成孔径内的PRF时刻，a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数；从步骤2中得到的平台副天线距离压缩后数据民的第is行中取第个数据的前Wtl个数据和后Wtl个数据，采用标准辛格插值方法对这组数据进行插值，得到插值重采样后的数据己((,,Ar)，其中Wtl为标准辛格插值的半个窗长； 步骤6. 4 :根据补偿相位因子的计算公式瓦=，得到散射点P(a，r)在PRF时刻is应补偿的相位因子，其中j为虚数单位(即_1开根)，I(^r)为步骤6. 2中得到的雷达平台主天线is时刻处散射点P(a，r)的距离史，λ为雷达系统工作的信号波长，λ由步骤I初始化得到，a表示散射点位于场景方位向的第a个位置，a=l,. . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数；将步骤6. 3中得到的插值重采样后的数据C ，;·)与散射点P(a，r)在PRF时刻is应补偿的相位因子相乘，得到相位补偿后的数据把雷达平台副天线距离散射点P(a，r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加，得到雷达平台副天线场景散射点P(a，r)成像后的数据g (ajr),BP s(a,r) = Tsa；fl,r) ’其中is为步骤6. 2中得到的雷达平台副天线距离散射点P (a, r)前后半个合成孔径内的PRF时刻； 步骤6. 5 :采用传统的GPU数据拷贝方法，将步骤6. 4中得到的雷达平台副天线场景所有散射点成像后的数据民(《，〃)， = ，sa为场景方位向总的散射点数，r=l,. . . , sr, sr为场景距离向总的散射点数，从图形处理器GPU的E复制到计算机内存中，其中E_￡ 为步骤3中在图形处理器GPU的全局存储器上分配的存储空间； 步骤7 :将步骤5中得到的雷达平台主天线场景所有散射点成像后的数据iw从r)和步骤6中得到的雷达平台副天线场景所有散射点成像后的数据民(《j进行共轭相乘，得到InSAR成像的干涉相位,完成InSAR数据成像配准一体化处理，其中a表示散射点位于场景方位向的第a个位置，a=l, . . . , sa, sa为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=l,. . . , sr, sr为场景距离向总的散射点数。
全文摘要
本发明公开了一种基于后向投影InSAR成像配准的GPU实现方法，它是利用后向投影算法将不同天线测量数据投影到同一个成像空间中进行成像，在成像过程中同时完成了InSAR不同天线的高保相成像和图像配准，实现了InSAR成像配准一体化，并结合GPU并行化处理技术实现后向投影算法的快速处理。相对于传统InSAR干涉相位生成，本发明具有保相性能好、成像配准一体化、并行化处理提高算法效率等优点，实现了高精度InSAR干涉相位的快速生成。本发明可以应用于合成孔径雷达成像，地球遥感等领域。
文档编号G01S13/90GK102788979SQ201210252090
公开日2012年11月21日 申请日期2012年7月20日 优先权日2012年7月20日
发明者付涛, 张晓玲, 韦顺军 申请人:电子科技大学