CN102788979B - 一种基于后向投影InSAR成像配准的GPU实现方法 - Google Patents
一种基于后向投影InSAR成像配准的GPU实现方法 Download PDFInfo
- Publication number
- CN102788979B CN102788979B CN201210252090.4A CN201210252090A CN102788979B CN 102788979 B CN102788979 B CN 102788979B CN 201210252090 A CN201210252090 A CN 201210252090A CN 102788979 B CN102788979 B CN 102788979B
- Authority
- CN
- China
- Prior art keywords
- scene
- thread
- scattering point
- data
- scattering
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于后向投影InSAR成像配准的GPU实现方法,它是利用后向投影算法将不同天线测量数据投影到同一个成像空间中进行成像,在成像过程中同时完成了InSAR不同天线的高保相成像和图像配准,实现了InSAR成像配准一体化,并结合GPU并行化处理技术实现后向投影算法的快速处理。相对于传统InSAR干涉相位生成,本发明具有保相性能好、成像配准一体化、并行化处理提高算法效率等优点,实现了高精度InSAR干涉相位的快速生成。本发明可以应用于合成孔径雷达成像,地球遥感等领域。
Description
技术领域
本技术发明属于雷达技术领域,它特别涉及了合成孔径雷达(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={Di|(xi,yi,zi),i∈D}
其中(xi,yi)是第i个地面像素点对应的平面坐标,zi是对应的高程坐标。
定义2、合成孔径雷达干涉测量(InSAR)
合成孔径雷达干涉测量(InSAR)指利用两个或者两个以上的SAR数据中的相位信息进行相干处理,结合雷达参数和雷达几何位置信息反演地表三维及其变化信息的遥感技术,详见文献“合成孔径雷达成像原理”,皮亦鸣等编著,电子科技大学出版社出版。
定义3、雷达成像空间
雷达成像空间是指将场景空间中的散射点投影到距离向-方位向的二维空间坐标系,该空间由合成孔径雷达成像空间中的两个相互正交的坐标基确定。目前典型合成孔径雷达的成像空间包括距离向-方位向投影空间。本发明中用以下数学关系表示成像空间M:
定义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=[x1,x2,…,xN]T为N×1维离散信号,向量XLP范数为 L1范数为 L2范数为
定义9、合成孔径雷达标准距离压缩方法
合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射参数,采用以下公式生成参考信号,并采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。
其中,j为虚数单位(即-1开根),f(t)为距离压缩参考函数,B为雷达发射基带信号的信号带宽,Tr为雷达发射信号脉冲宽度,t为时间变量,取值范围从到详见文献“雷达成像技术”,保铮等编著,电子工业出版社出版。
定义10、合成孔径雷达原始回波仿真方法
合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定***参数条件下具有SAR回波信号特性的原始信号的方法,详细内容可参考文献:“InSAR回波信号与***仿真研究”,张剑琦,哈尔滨工业大学硕士论文。
定义11、天线相位中心
天线相位中心是指雷达天线向外辐射信号的中心,本发明中天线相位中心指雷达天线的轨迹位置。
定义12、辛格插值(sinc插值)方法
辛格插值方法是指对于一个带限信号,在满足采样定理的情况下,采用卷积核为sinc的函数h(x),h(x)的长度即窗长为W。
进行对已离散的信号gd(i)插值,得到插值后所要的信号
详见文献“合成孔径雷达成像----算法与实现”,Frank H.Wong等编著,电子工业出版社出版。
定义13、取整函数
取整函数是指将任一实数映射到相近的整数的一类函数,本发明中用到三个取整函数,分别为ceil、floor和round。
ceil(x)为向上取整函数,表示取不超过x的最大一个整数。
floor(x)为向下取整函数,表示取不小于x的最小一个整数。
round(x)为近似取整函数,表示按四舍五入原则对x取整。
定义14、合成孔径与PRF时刻
合成孔径是指对于测绘场景中的一个散射点从进入雷达波束照射范围至离开雷达波束照射范围的这段时间内,雷达波束中心所走过的长度。
雷达平台飞过一个合成孔径所需要的时间称为慢时间,雷达***以一定的重复周期Tr发射接收脉冲,因此慢时间可以表示为一个以重复周期Tr为步长的离散化时间变量,其中每一个离散时间变量值为一个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)是内核函数执行并行计算的第二层并行单元,同一线程块中的线程可以共享数据,没有执行顺序。
详见文献“GPU高性能运算之CUDA”,张舒等编著,中国水利水电出版社出版。
定义17、线程块索引与线程索引
线程块索引是指每个线程块在线程网格中的位置,是CUDA的内建变量,线程块的第一维索引为blockIdx.x,线程块的第二维索引为blockIdx.y,线程块的第三维索引恒为1。
线程索引是指每个线程在线程块中的位置,是CUDA的内建变量,线程的第一维索引为threadIdx.x,线程的第二维索引为threadIdx.y,线程的第三维索引为threadIdx.z。
详见文献“GPU高性能运算之CUDA”,张舒等编著,中国水利水电出版社出版。
定义18、GPU全局存储器
GPU全局存储器(global memory)是GPU片外的存储器,可读写,作用范围是整个程序,整个线程网格中的任意线程都能读写全局存储器的任意位置,存储空间大,生命周期为整个程序。
详见文献“GPU高性能运算之CUDA”,张舒等编著,中国水利水电出版社出版。
定义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雷达平台包含两组天线,即主天线和副天线,两组天线之间的距离为基线长,记做Bl,主天线发射脉冲信号,经过Td时间的延迟,主天线和副天线同时接收回波延迟信号。雷达平台主天线接收的回波数据,记做雷达平台副天线接收的回波数据,记做其中和均为二维矩阵,第一维均对应方位向,第二维均对应距离向,即二维矩阵和的行存储的是方位向数据,二维矩阵和的列存储的是距离向数据。
初始化InSAR成像***参数包括:雷达***工作的信号波长,记做λ,雷达平台主天线发射信号带宽,记做B,雷达平台主天线发射脉冲时宽,记做Tr,雷达平台接收***采样频率,记做Fs,雷达***脉冲重复频率,记做prf,雷达平台主天线合成孔径长度,记做Lm,雷达平台副天线合成孔径长度,记做Ls,雷达平台速度矢量,记做雷达平台主天线初始位置矢量,记做雷达平台副天线初始位置矢量,记做场景参考点位置矢量,记做雷达***距离向采样点数,记做Nr,雷达***方位向采样点数,记做Na,场景距离向散射点间隔,记做dr,场景方位向散射点间隔,记做da,场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值,记做Rmc,场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值,记做Rsc,场景参考点在主天线回波数据和副天线回波数据中的距离门相同,距离门位置记做Ic,线程网格中线程块的第一维索引,记做blockIdx.x,线程网格中线程块的第二维索引,记做blockIdx.y,每个线程块中线程的第一维索引,记做threadIdx.x,每个线程块中线程的第二维索引,记做threadIdx.y,blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量。
采用公式计算得到雷达平台主天线第n个PRF时刻的天线相位中心矢量采用公式计算得到雷达平台副天线第n个PRF时刻的天线相位中心矢量n表示第n个PRF时刻,主天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小,记做副天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小,记做采用公式 计算得到场景散射点P(a,r)的位置矢量,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,场景所有散射点的位置矢量占用的存储空间大小记做为雷达平台主天线场景散射点P(a,r)成像后的数据,初始化为0,a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小记做 为雷达平台副天线场景散射点P(a,r)成像后的数据,初始化为0,a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小记做
步骤2:InSAR原始回波数据进行距离压缩
步骤3:为图形处理器GPU分配数据
采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做。
采用传统的GPU数据拷贝方法,将步骤2中得到的平台主天线距离压缩后数据从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤2中得到的平台副天线距离压缩后数据从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤1中得到的主天线所有Na个PRF时刻的天线相位中心矢量n=1,...,Na,Na为雷达***方位向采样点数,从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤1中得到的副天线所有Na个PRF时刻的天线相位中心矢量n=1,...,Na,Na为雷达***方位向采样点数,从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤1中得到的场景所有散射点的位置矢量a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,从计算机内存复制到图形处理器GPU的存储空间中。
步骤4:为图形处理器GPU内核函数划分线程网格
线程网格中的线程块以二维形式划分,线程块中的线程也以二维形式划分,线程网格第一维总的线程数等于场景方位向总的散射点数sa,线程网格第二维总的线程数等于场景距离向总的散射点数sr,即一个场景散射点对应一个线程,具体划分步骤是:
首先确定线程块中的线程数,方法是:Db_x为线程块第一维分配的线程数,Db_x设置为8的整数倍,同时不超过GPU线程硬件规定的范围;Db_y为线程块第二维分配的线程数,Db_y设置为8的整数倍,同时不超过GPU线程硬件规定的范围。
然后确定线程块的数目,方法是:Dg_x为线程网格第一维分配的线程块数目,Dg_x由场景方位向总的散射点数sa和线程块第一维分配的线程数Db_x共同确定,计算公式为Dg_x=floor((sa+Db_x-1)/Db_x);Dg_y为线程网格第二维分配的线程块数目,Dg_y由场景距离向总的散射点数sr和线程块第二维分配的线程数Db_y共同确定,计算公式为Dg_y=floor((sr+Db_y-1)/Db_y),其中floor为向下取整函数。
步骤5:平台主天线距离压缩后数据沿方位向成像
步骤5.1:确定场景散射点与线程网格中线程的对应关系,方法是:对于场景散射点P(a,r),a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数,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:lm为主天线半个合成孔径内的PRF时刻个数,计算公式为lm=round(Lmda2),其中round为近似取整函数,Lm为主天线合成孔径长度,da为场景方位向散射点间隔,Lm和da均由步骤1初始化得到;im为主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,im的取值满足条件:im=a-lm,...,a+lm,若a-lm<1,则a-lm=1,若a+lm>Na,则a+lm=Na,其中Na为雷达***方位向采样点数,即im只取雷达***采样范围内的主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。采用公式计算得到雷达平台主天线im时刻处散射点P(a,r)的距离史其中||□||2为L2范数,为步骤1中得到的雷达平台主天线第im个PRF时刻的天线相位中心矢量,为步骤1中得到的场景散射点P(a,r)的位置矢量,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。
步骤5.3:为场景散射点P(a,r)的距离史对应的距离门,计算公式为其中round为近似取整函数,为步骤5.2中得到的雷达平台主天线im时刻处散射点P(a,r)的距离史,Rmc为场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值,Ic为场景参考点在主天线回波数据中的距离门,dr为场景距离向散射点间隔,Rmc、Ic和dr均由步骤1初始化得到,im为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。从步骤2中得到的平台主天线距离压缩后数据的第im行中取第个数据的前W0个数据和后W0个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据其中W0为标准辛格插值的半个窗长。
步骤5.4:根据补偿相位因子的计算公式得到散射点P(a,r)在PRF时刻im应补偿的相位因子其中j为虚数单位(即-1开根),为步骤5.2中得到的雷达平台主天线im时刻处散射点P(a,r)的距离史,λ为雷达***工作的信号波长,λ由步骤1初始化得到,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。将步骤5.3中得到的插值重采样后的数据与散射点P(a,r)在PRF时刻im应补偿的相位因子相乘,得到相位补偿后的数据把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加,得到雷达平台主天线场景散射点P(a,r)成像后的数据即其中im为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。
步骤5.5:采用传统的GPU数据拷贝方法,将步骤5.4中得到的雷达平台主天线场景所有散射点成像后的数据a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,从图形处理器GPU的复制到计算机内存中,其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。
步骤6:平台副天线距离压缩后数据沿方位向成像
图形处理器GPU实现平台副天线距离压缩后数据沿方位向成像的具体步骤是:
步骤6.1:确定场景散射点与线程网格中线程的对应关系,方法是:对于场景散射点P(a,r),a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数,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:ls为副天线半个合成孔径内的PRF时刻个数,计算公式为ls=round(Lsda2),其中round为近似取整函数,Ls为副天线合成孔径长度,da为场景方位向散射点间隔,Ls和da均由步骤1初始化得到;is为副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,is的取值满足条件:is=a-ls,...,a+ls,若a-ls<1,则a-ls=1,若a+ls>Na,则a+ls=Na,其中Na为雷达***方位向采样点数,即is只取雷达***采样范围内的副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。采用公式计算得到雷达平台副天线is时刻处散射点P(a,r)的距离史采用公式计算得到雷达平台主天线is时刻处散射点P(a,r)的距离史其中||□||2为L2范数,为步骤1中得到的雷达平台副天线第is个PRF时刻的天线相位中心矢量,为步骤1中得到的雷达平台主天线第is个PRF时刻的天线相位中心矢量,为步骤1中得到的场景散射点P(a,r)的位置矢量,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。
步骤6.3:为场景散射点P(a,r)的距离史对应的距离门,计算公式为 其中round为近似取整函数,为步骤6.2中得到的雷达平台副天线is时刻处散射点P(a,r)的距离史,Rsc为场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值,Ic为场景参考点在主天线回波数据中的距离门,dr为场景距离向散射点间隔,Rsc、Ic和dr均由步骤1初始化得到,is为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。从步骤2中得到的平台副天线距离压缩后数据的第is行中取第个数据的前W0个数据和后W0个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据其中W0为标准辛格插值的半个窗长。
步骤6.4:根据补偿相位因子的计算公式 得到散射点P(a,r)在PRF时刻is应补偿的相位因子其中j为虚数单位(即-1开根),为步骤6.2中得到的雷达平台主天线is时刻处散射点P(a,r)的距离史,λ为雷达***工作的信号波长,λ由步骤1初始化得到,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。将步骤6.3中得到的插值重采样后的数据与散射点P(a,r)在PRF时刻is应补偿的相位因子相乘,得到相位补偿后的数据把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加,得到雷达平台副天线场景散射点P(a,r)成像后的数据即其中is为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。
步骤6.5:采用传统的GPU数据拷贝方法,将步骤6.4中得到的雷达平台副天线场景所有散射点成像后的数据a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,从图形处理器GPU的复制到计算机内存中,其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。
步骤7:将步骤5中得到的雷达平台主天线场景所有散射点成像后的数据和步骤6中得到的雷达平台副天线场景所有散射点成像后的数据进行共轭相乘,得到InSAR成像的干涉相位完成InSAR数据成像配准一体化处理,其中a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数。
本发明的创新点在于利用后向投影算法,将InSAR不同天线测量数据投影到同一个成像空间进行成像,实现InSAR数据成像和图像配准一体化,并根据后向投影算法逐点相干积累原理,利用GPU并行化技术实现算法快速处理。通过仿真实验验证,相对于传统InSAR干涉相位生成,该算法具有保相性能好、成像配准一体化、并行化处理提高算法效率等优点,实现了高精度InSAR干涉相位的快速生成。
本发明的优点在于在数据成像过程中同时完成InSAR高保相成像与配准,在数据成像时同时获得高质量干涉相位,不需要再进行图像配准处理,并且可以利用GPU并行化处理提高算法效率。本发明可以应用于合成孔径雷达成像,地球遥感等领域。
附图说明
图1为本发明所提供方法的流程示意框图;
图2为本发明具体实施方式采用的干涉合成孔径雷达平台飞行几何关系及仿真圆锥场景图;
其中,横坐标(X轴)为方位向,纵坐标(Y轴)为距离向,垂直坐标(Z轴)为高度向,Bl为雷达平台主天线和副天线之间的基线,基线的两个端点分别为雷达平台主天线和副天线,P(a,r)为场景中任一散射点,a表示散射点位于场景方位向的第a个位置,a=1,...,512,r表示散射点位于场景距离向的第r个位置,r=1,...,512,为雷达平台飞行速度矢量,为雷达平台主天线第im个PRF时刻的天线相位中心矢量,为雷达平台副天线第is个PRF时刻的天线相位中心矢量,为雷达平台主天线im时刻处散射点P(a,r)的距离史,为雷达平台副天线is时刻处散射点P(a,r)的距离史,im为雷达平台主天线距离散射点P(a,r)半个合成孔径内的PRF时刻,is为雷达平台副天线距离散射点P(a,r)半个合成孔径内的PRF时刻,sa为仿真场景方位向总的散射点数,sr为仿真场景距离向总的散射点数,da为仿真场景方位向散射点间隔,dr为仿真场景距离向散射点间隔。
图3本发明具体实施方式采用的线程网格划分;
其中,Block为线程网格中线程块的二维划分,Thread为第一维索引为0、第二维索引为15的线程块的二维线程划分。
具体实施方式
本发明主要采用仿真实验的方法进行验证,所有步骤、结论都在VC++、MATLAB7.0、CUDA4.0和GeForce GTX560Ti显卡上验证正确。具体实施步骤如下:
步骤1:初始化InSAR成像***参数
选择与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基 该坐标基方向为方位向;选择InSAR成像空间的第二个坐标基 该坐标基方向为距离向。
初始化InSAR成像***参数包括:雷达***工作的信号波长λ=0.03m,雷达平台主天线发射信号带宽B=150MHz,雷达平台主天线发射脉冲时宽Tr=1μs,雷达平台接收***采样频率Fs=300MHz,雷达***脉冲重复频率prf=500Hz,雷达***距离向采样点数Nr=1024,雷达***方位向采样点数Na=2048,场景距离向散射点间隔dr=0.5m,场景方位向散射点间隔da=0.3m,雷达平台主天线合成孔径长度Lm=150m,雷达平台副天线合成孔径长度Ls=150m,雷达平台速度矢量 速度的单位为m/s,雷达平台主天线初始位置矢量 雷达平台副天线初始位置矢量 场景参考位置矢量 位置的单位为m,场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值Rmc=10000m,场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值Rsc=9992,场景参考点在主天线回波数据和副天线回波数据中的距离门Ic=1。
仿真场景为一立体圆锥,仿真场景方位向总的散射点数sa=512,仿真场景距离向总的散射点数sr=512,圆锥最小高度为0m,最大高度为50m。
采用公式计算得到雷达平台主天线第n个PRF时刻的天线相位中心矢量采用公式计算得到雷达平台副天线第n个PRF时刻的天线相位中心矢量n=1,...,2048,主天线所有2048个PRF时刻的天线相位中心矢量占用的存储空间大小副天线所有2048个PRF时刻的天线相位中心矢量占用的存储空间大小采用公式 计算得到场景散射点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内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做
采用传统的GPU数据拷贝方法,将步骤2中得到的平台主天线距离压缩后数据从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤2中得到的平台副天线距离压缩后数据从计算机内存复制到图形处理器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:平台主天线距离压缩后数据沿方位向成像
步骤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:lm为主天线半个合成孔径内的PRF时刻个数,lm=250;im为主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,im的取值满足条件:im=a-250,…,a+250,若a-250<1,则a-lm=1,若a+250>2048,则a+lm=2048。采用公式计算得到雷达平台主天线im时刻处散射点P(a,r)的距离史其中||□||2为L2范数,为步骤1中得到的雷达平台主天线第im个PRF时刻的天线相位中心矢量,为步骤1中得到的场景散射点P(a,r)的位置矢量,a表示散射点位于场景方位向的第a个位置,a=1,...,5,r表示散射点位于场景距离向的第r个位置,r=1,...,512。
步骤5.3:为场景散射点P(a,r)的距离史对应的距离门,计算公式为)其中round为近似取整函数,为步骤5.2中得到的雷达平台主天线im时刻处散射点P(a,r)的距离史,im为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,512,r表示散射点位于场景距离向的第r个位置,r=1,...,512。从步骤2中得到的平台主天线距离压缩后数据的第im行中取第个数据的前4个数据和后4个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据
步骤5.4:根据补偿相位因子的计算公式 得到散射点P(a,r)在PRF时刻im应补偿的相位因子其中j为虚数单位(即-1开根),为步骤5.2中得到的雷达平台主天线im时刻处散射点P(a,r)的距离史,a表示散射点位于场景方位向的第a个位置,a=1,...,512,r表示散射点位于场景距离向的第r个位置,r=1,...,512。将步骤5.3中得到的插值重采样后的数据与散射点P(a,r)在PRF时刻im应补偿的相位因子相乘,得到相位补偿后的数据把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加,得到雷达平台主天线场景散射点P(a,r)成像后的数据即其中im为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。
步骤5.5:采用传统的GPU数据拷贝方法,将步骤5.4中得到的雷达平台主天线场景所有散射点成像后的数据从图形处理器GPU的复制到计算机内存中,其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。
步骤6:平台副天线距离压缩后数据沿方位向成像
步骤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:ls为副天线半个合成孔径内的PRF时刻个数,ls=250;is为副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,is的取值满足条件:is=a-250,...,a+250,若a-250<1,则a-250=1,若a+250>2048,则a+250=2048。采用公式计算得到雷达平台副天线is时刻处散射点P(a,r)的距离史采用公式计算得到雷达平台主天线is时刻处散射点P(a,r)的距离史其中||□||2为L2范数,为步骤1中得到的雷达平台副天线第is个PRF时刻的天线相位中心矢量,为步骤1中得到的雷达平台主天线第is个PRF时刻的天线相位中心矢量,为步骤1中得到的场景散射点P(a,r)的位置矢量,a表示散射点位于场景方位向的第a个位置,a=1,...,512,r表示散射点位于场景距离向的第r个位置,r=1,...,512。步骤6.3:为场景散射点P(a,r)的距离史对应的距离门,计算公式为 其中round为近似取整函数,为步骤6.2中得到的雷达平台副天线is时刻处散射点P(a,r)的距离史,is为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,512,r表示散射点位于场景距离向的第r个位置,r=1,...,512。从步骤2中得到的平台副天线距离压缩后数据的第is行中取第个数据的前4个数据和后4个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据
步骤6.4:根据补偿相位因子的计算公式 得到散射点P(a,r)在PRF时刻is应补偿的相位因子其中j为虚数单位(即-1开根),为步骤6.2中得到的雷达平台主天线is时刻处散射点P(a,r)的距离史,a表示散射点位于场景方位向的第a个位置,a=1,...,512,r表示散射点位于场景距离向的第r个位置,r=1,...,512。将步骤6.3中得到的插值重采样后的数据与散射点P(a,r)在PRF时刻is应补偿的相位因子相乘,得到相位补偿后的数据把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加,得到雷达平台副天线场景散射点P(a,r)成像后的数据即其中is为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻。
步骤6.5:采用传统的GPU数据拷贝方法,将步骤6.4中得到的雷达平台副天线场景所有散射点成像后的数据从图形处理器GPU的复制到计算机内存中,其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间。
步骤7:将步骤5中得到的雷达平台主天线场景所有散射点成像后的数据和步骤6中得到的雷达平台副天线场景所有散射点成像后的数据进行共轭相乘,得到InSAR成像的干涉相位完成InSAR数据成像配准一体化处理,其中a表示散射点位于场景方位向的第a个位置,a=1,...,512,r=1,...,512。
通过本发明具体实施方式的仿真及测试,本发明所提供的基于后向投影InSAR成像配准的GPU实现方法能够实现InSAR不同天线成像和图像配准同时完成,与现有的InSAR干涉相位生成方法相比,本发明克服现有方法数据成像和图像配准相互隔离的缺点,可直接从原始数据成像处理过程中获得配准后的干涉相位,减少数据处理流程并且提高了成像相位保持精度,为InSAR高质量干涉相位生成提供了一种新方法。同时,利用GPU并行化技术实现了算法的快速处理。
Claims (1)
1.一种基于后向投影InSAR成像配准的GPU实现方法,其特征是它包括以下步骤:
步骤1:初始化InSAR成像***参数
InSAR成像空间由InSAR成像空间中的两个相互正交的坐标基确定,定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基,记做该坐标基方向为方位向;定义在地平面内,并与InSAR成像空间的第一个坐标基垂直的单位向量作为InSAR成像空间的第二个坐标基,记做该坐标基方向为距离向;
InSAR雷达平台包含两组天线,即主天线和副天线,两组天线之间的距离为基线长,记做Bl,主天线发射脉冲信号,经过Td时间的延迟,主天线和副天线同时接收回波延迟信号;雷达平台主天线接收的回波数据,记做雷达平台副天线接收的回波数据,记做其中和均为二维矩阵,第一维均对应方位向,第二维均对应距离向,即二维矩阵和的行存储的是方位向数据,二维矩阵和的列存储的是距离向数据;
初始化InSAR成像***参数包括:雷达***工作的信号波长,记做λ,雷达平台主天线发射信号带宽,记做B,雷达平台主天线发射脉冲时宽,记做Tr,雷达平台接收***采样频率,记做Fs,雷达***脉冲重复频率,记做prf,雷达平台主天线合成孔径长度,记做Lm,雷达平台副天线合成孔径长度,记做Ls,雷达平台速度矢量,记做雷达平台主天线初始位置矢量,记做雷达平台副天线初始位置矢量,记做场景参考点位置矢量,记做雷达***距离向采样点数,记做Nr,雷达***方位向采样点数,记做Na,场景距离向散射点间隔,记做dr,场景方位向散射点间隔,记做da,场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值,记做Rmc,场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值,记做Rsc,场景参考点在主天线回波数据和副天线回波数据中的距离门相同,距离门位置记做Ic,线程网格中线程块的第一维索引,记做blockIdx.x,线程网格中线程块的第二维索引,记做blockIdx.y,每个线程块中线程的第一维索引,记做threadIdx.x,每个线程块中线程的第二维索引,记做threadIdx.y,blockIdx.x、blockIdx.y、threadIdx.x和threadIdx.y均为计算统一设备架构CUDA的内建变量;
采用公式计算得到雷达平台主天线第n个PRF时刻的天线相位中心矢量采用公式计算得到雷达平台副天线第n个PRF时刻的天线相位中心矢量n=1,...,Na,n表示第n个PRF时刻,主天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小,记做副天线所有Na个PRF时刻的天线相位中心矢量占用的存储空间大小,记做采用公式 计算得到场景散射点P(a,r)的位置矢量,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,场景所有散射点的位置矢量占用的存储空间大小记做为雷达平台主天线场景散射点P(a,r)成像后的数据,初始化为0,a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,雷达平台主天线场景所有散射点成像后的数据占用的存储空间大小记做为雷达平台副天线场景散射点P(a,r)成像后的数据,初始化为0,a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,雷达平台副天线场景所有散射点成像后的数据占用的存储空间大小记做
步骤2:InSAR原始回波数据进行距离压缩
步骤3:为图形处理器GPU分配数据
采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做采用传统的GPU内存分配方法,在图形处理器GPU的全局存储器上分配大小为的存储空间,记做
采用传统的GPU数据拷贝方法,将步骤2中得到的平台主天线距离压缩后数据从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤2中得到的平台副天线距离压缩后数据从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤1中得到的主天线所有Na个PRF时刻的天线相位中心矢量n=1,...,Na,Na为雷达***方位向采样点数,从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤1中得到的副天线所有Na个PRF时刻的天线相位中心矢量n=1,...,Na,Na为雷达***方位向采样点数,从计算机内存复制到图形处理器GPU的存储空间中;采用传统的GPU数据拷贝方法,将步骤1中得到的场景所有散射点的位置矢量a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,从计算机内存复制到图形处理器GPU的存储空间中;
步骤4:为图形处理器GPU内核函数划分线程网格
线程网格中的线程块以二维形式划分,线程块中的线程也以二维形式划分,线程网格第一维总的线程数等于场景方位向总的散射点数sa,线程网格第二维总的线程数等于场景距离向总的散射点数sr,即一个场景散射点对应一个线程,具体划分步骤是:
首先确定线程块中的线程数,方法是:Db_x为线程块第一维分配的线程数,Db_x设置为8的整数倍,同时不超过GPU线程硬件规定的范围;Db_y为线程块第二维分配的线程数,Db_y设置为8的整数倍,同时不超过GPU线程硬件规定的范围;
然后确定线程块的数目,方法是:Dg_x为线程网格第一维分配的线程块数目,Dg_x由场景方位向总的散射点数sa和线程块第一维分配的线程数Db_x共同确定,计算公式为Dg_x=floor((sa+Db_x-1)/Db_x);Dg_y为线程网格第二维分配的线程块数目,Dg_y由场景距离向总的散射点数sr和线程块第二维分配的线程数Db_y共同确定,计算公式为Dg_y=floor((sr+Db_y-1)/Db_y),其中floor为向下取整函数;
步骤5:平台主天线距离压缩后数据沿方位向成像
步骤5.1:确定场景散射点与线程网格中线程的对应关系,方法是:对于场景散射点P(a,r),a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数,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:lm为主天线半个合成孔径内的PRF时刻个数,计算公式为lm=round(Lm/da/2),其中round为近似取整函数,Lm为主天线合成孔径长度,da为场景方位向散射点间隔,Lm和da均由步骤1初始化得到;im为主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,im的取值满足条件:im=a-lm,...,a+lm,若a-lm<1,则a-lm=1,若a+lm>Na,则a+lm=Na,其中Na为雷达***方位向采样点数,即im只取雷达***采样范围内的主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻;采用公式计算得到雷达平台主天线im时刻处散射点P(a,r)的距离史其中为L2范数,为步骤1中得到的雷达平台主天线第im个PRF时刻的天线相位中心矢量,为步骤1中得到的场景散射点P(a,r)的位置矢量,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数;
步骤5.3:为场景散射点P(a,r)的距离史对应的距离门,计算公式为 其中round为近似取整函数,为步骤5.2中得到的雷达平台主天线im时刻处散射点P(a,r)的距离史,Rmc为场景参考点到雷达平台主天线各PRF时刻天线相位中心距离的最小值,Ic为场景参考点在主天线回波数据中的距离门,dr为场景距离向散射点间隔,Rmc、Ic和dr均由步骤1初始化得到,im为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数;从步骤2中得到的平台主天线距离压缩后数据的第im行中取第个数据的前W0个数据和后W0个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据其中W0为标准辛格插值的半个窗长;
步骤5.4:根据补偿相位因子的计算公式 得到散射点P(a,r)在PRF时刻im应补偿的相位因子其中j为虚数单位(即-1开根),为步骤5.2中得到的雷达平台主天线im时刻处散射点P(a,r)的距离史,λ为雷达***工作的信号波长,λ由步骤1初始化得到,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数;将步骤5.3中得到的插值重采样后的数据与散射点P(a,r)在PRF时刻im应补偿的相位因子相乘,得到相位补偿后的数据把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加,得到雷达平台主天线场景散射点P(a,r)成像后的数据即其中im为步骤5.2中得到的雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻;
步骤5.5:采用传统的GPU数据拷贝方法,将步骤5.4中得到的雷达平台主天线场景所有散射点成像后的数据a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,从图形处理器GPU的复制到计算机内存中,其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间;
步骤6:平台副天线距离压缩后数据沿方位向成像
图形处理器GPU实现平台副天线距离压缩后数据沿方位向成像的具体步骤是:
步骤6.1:确定场景散射点与线程网格中线程的对应关系,方法是:对于场景散射点P(a,r),a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数,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:ls为副天线半个合成孔径内的PRF时刻个数,计算公式为ls=round(Ls/da/2),其中round为近似取整函数,Ls为副天线合成孔径长度,da为场景方位向散射点间隔,Ls和da均由步骤1初始化得到;is为副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,is的取值满足条件:is=a-ls,...,a+ls,若a-ls<1,则a-ls=1,若a+ls>Na,则a+ls=Na,其中Na为雷达***方位向采样点数,即is只取雷达***采样范围内的副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻;采用公式计算得到雷达平台副天线is时刻处散射点P(a,r)的距离史采用公式计算得到雷达平台主天线is时刻处散射点P(a,r)的距离史其中为L2范数,为步骤1中得到的雷达平台副天线第is个PRF时刻的天线相位中心矢量,为步骤1中得到的雷达平台主天线第is个PRF时刻的天线相位中心矢量,为步骤1中得到的场景散射点P(a,r)的位置矢量,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数;
步骤6.3:为场景散射点P(a,r)的距离史对应的距离门,计算公式为 其中round为近似取整函数,为步骤6.2中得到的雷达平台副天线is时刻处散射点P(a,r)的距离史,Rsc为场景参考点到雷达平台副天线各PRF时刻天线相位中心距离的最小值,Ic为场景参考点在主天线回波数据中的距离门,dr为场景距离向散射点间隔,Rsc、Ic和dr均由步骤1初始化得到,is为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数;从步骤2中得到的平台副天线距离压缩后数据的第is行中取第个数据的前W0个数据和后W0个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据其中W0为标准辛格插值的半个窗长;
步骤6.4:根据补偿相位因子的计算公式 得到散射点P(a,r)在PRF时刻is应补偿的相位因子其中j为虚数单位(即-1开根),为步骤6.2中得到的雷达平台主天线is时刻处散射点P(a,r)的距离史,λ为雷达***工作的信号波长,λ由步骤1初始化得到,a表示散射点位于场景方位向的第a个位置,a=1,...,sa,sa为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,sr,sr为场景距离向总的散射点数;将步骤6.3中得到的插值重采样后的数据与散射点P(a,r)在PRF时刻is应补偿的相位因子相乘,得到相位补偿后的数据把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有PRF时刻的相位补偿后的数据相加,得到雷达平台副天线场景散射点P(a,r)成像后的数据即其中is为步骤6.2中得到的雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的PRF时刻;
步骤6.5:采用传统的GPU数据拷贝方法,将步骤6.4中得到的雷达平台副天线场景所有散射点成像后的数据a=1,...,sa,sa为场景方位向总的散射点数,r=1,...,sr,sr为场景距离向总的散射点数,从图形处理器GPU的复制到计算机内存中,其中为步骤3中在图形处理器GPU的全局存储器上分配的存储空间;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210252090.4A CN102788979B (zh) | 2012-07-20 | 2012-07-20 | 一种基于后向投影InSAR成像配准的GPU实现方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210252090.4A CN102788979B (zh) | 2012-07-20 | 2012-07-20 | 一种基于后向投影InSAR成像配准的GPU实现方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102788979A CN102788979A (zh) | 2012-11-21 |
CN102788979B true CN102788979B (zh) | 2014-04-09 |
Family
ID=47154431
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201210252090.4A Expired - Fee Related CN102788979B (zh) | 2012-07-20 | 2012-07-20 | 一种基于后向投影InSAR成像配准的GPU实现方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102788979B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103809177A (zh) * | 2013-12-30 | 2014-05-21 | 南京大学 | 一种基于fpga的雷达成像并行化方法 |
CN103869315B (zh) * | 2014-03-18 | 2016-05-25 | 电子科技大学 | 临近空间圆周合成孔径雷达快速后向投影成像方法 |
CN104181510A (zh) * | 2014-09-05 | 2014-12-03 | 电子科技大学 | 一种并行雷达脉冲压缩方法 |
CN104375136A (zh) * | 2014-11-18 | 2015-02-25 | 无锡悟莘科技有限公司 | 一种主天线为反弓形天线的雷达*** |
CN106646664B (zh) | 2016-10-09 | 2018-10-23 | 华讯方舟科技有限公司 | 基于gpu的人体微波回波模拟方法及*** |
CN107180014B (zh) * | 2017-04-28 | 2018-10-23 | 华讯方舟科技有限公司 | 一种快速sinc插值方法及*** |
CN108008389B (zh) * | 2017-12-01 | 2019-12-10 | 电子科技大学 | 一种基于gpu的快速频域后向投影三维成像方法 |
CN108711335A (zh) * | 2018-06-15 | 2018-10-26 | 青岛擎鹰信息科技有限责任公司 | 一种实现分布式大场景雷达成像仿真方法及其*** |
CN113640797B (zh) * | 2021-08-09 | 2022-04-12 | 北京航空航天大学 | 一种用于参考条带模式InSAR的前斜视测高方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101441271A (zh) * | 2008-12-05 | 2009-05-27 | 航天恒星科技有限公司 | 基于gpu的sar实时成像处理设备 |
CN101937082A (zh) * | 2009-07-02 | 2011-01-05 | 北京理工大学 | 基于gpu众核平台的合成孔径雷达并行成像方法 |
-
2012
- 2012-07-20 CN CN201210252090.4A patent/CN102788979B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101441271A (zh) * | 2008-12-05 | 2009-05-27 | 航天恒星科技有限公司 | 基于gpu的sar实时成像处理设备 |
CN101937082A (zh) * | 2009-07-02 | 2011-01-05 | 北京理工大学 | 基于gpu众核平台的合成孔径雷达并行成像方法 |
Non-Patent Citations (6)
Title |
---|
A fast InSAR raw signal simulation based on GPGPU;Zhihua He 等;《International Conference on Graphic and Image Processing (ICGIP 2011)》;20110930;第8285卷;82855O-1-82855O-7 * |
Accelerating InSAR raw data simulation on GPU using CUDA;Zhang Fan等;《Geoscience and Remote Sensing Symposium (IGARSS), 2010 IEEE International》;20100730;2932-2935 * |
An efficient system for creating synthetic InSAR images from simulations;Xiaoru Yuan等;《Pure and Applied Geophysics》;20080430;第165卷(第3-4期);671-691 * |
Xiaoru Yuan等.An efficient system for creating synthetic InSAR images from simulations.《Pure and Applied Geophysics》.2008,第165卷(第3-4期),671-691. |
Zhang Fan等.Accelerating InSAR raw data simulation on GPU using CUDA.《Geoscience and Remote Sensing Symposium (IGARSS), 2010 IEEE International》.2010,2932-2935. |
Zhihua He 等.A fast InSAR raw signal simulation based on GPGPU.《International Conference on Graphic and Image Processing (ICGIP 2011)》.2011,第8285卷82855O-1-82855O-7. |
Also Published As
Publication number | Publication date |
---|---|
CN102788979A (zh) | 2012-11-21 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102788979B (zh) | 一种基于后向投影InSAR成像配准的GPU实现方法 | |
CN103713288B (zh) | 基于迭代最小化稀疏贝叶斯重构线阵sar成像方法 | |
CN103018740B (zh) | 一种基于曲面投影的InSAR成像方法 | |
CN105677942B (zh) | 一种重复轨道星载自然场景sar复图像数据快速仿真方法 | |
CN109471080B (zh) | 基于simulink的高速平台雷达回波信号模拟*** | |
CN104730520B (zh) | 基于子孔径合成的圆周sar后向投影自聚焦方法 | |
Benetazzo et al. | Stereo wave imaging from moving vessels: Practical use and applications | |
CN103576137B (zh) | 一种基于成像策略的多传感器多目标定位方法 | |
CN101876704B (zh) | 干涉合成孔径雷达三维陆地场景回波仿真方法 | |
CN102854507B (zh) | 一种基于gpu后向投影双站合成孔径雷达成像方法 | |
CN103698763A (zh) | 基于硬阈值omp的线阵sar稀疏成像方法 | |
CN104391295A (zh) | 一种图像熵最优的压缩传感sar稀疏自聚焦成像方法 | |
CN103941243A (zh) | 一种基于sar三维成像的自旋式飞行器测高方法 | |
CN102662171A (zh) | 一种sar层析三维成像方法 | |
CN103616682B (zh) | 一种基于曲面投影的多基线InSAR处理方法 | |
CN106680812A (zh) | 一种基于解析面元的微波关联成像仿真方法 | |
CN104698457A (zh) | 一种迭代曲面预测InSAR成像及高度估计方法 | |
CN114047511B (zh) | 一种基于csa算法的时变海面机载sar成像仿真方法 | |
CN105699969A (zh) | 基于广义高斯约束的最大后验估计角超分辨成像方法 | |
CN103207387A (zh) | 一种机载相控阵pd雷达杂波的快速模拟方法 | |
CN110646782A (zh) | 一种基于波形匹配的星载激光在轨指向检校方法 | |
CN103778633B (zh) | 确定数字高程模型单元网格遮挡的方法及装置 | |
CN106872977A (zh) | 一种基于分段弱正交匹配追踪的层析sar三维成像方法 | |
Willis et al. | Hardware-accelerated SAR simulation with NVIDIA-RTX technology | |
CN103076608A (zh) | 轮廓增强的聚束式合成孔径雷达成像方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20140409 Termination date: 20170720 |