资源描述：
《基于后向投影的sar成像算法与gpu加速研究》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。
中图分类号：TN957论文编号：102870415-S094学科分类号：081001硕士学位论文基于后向投影的SAR成像算法与GPU加速研究研究生姓名班阳阳学科、专业通信与信息系统研究方向雷达成像指导教师张劲东副教授南京航空航天大学研究生院电子信息工程学院二О一四年十二月万方数据 万方数据 NanjingUniversityofAeronauticsandAstronauticsTheGraduateSchoolCollegeofElectronicandInformationEngineeringSARImagingAlgorithmandGPUAccelerationResearchBasedonBackProjectionAThesisinCommunicationandInformationSystembyBanYangyangAdvisedbyAssociateProf.ZhangJindongSubmittedinPartialFulfillmentoftheRequirementsfortheDegreeofMasterofEngineeringDecember,2014万方数据 万方数据 承诺书本人声明所呈交的硕士学位论文是本人在导师指导下进行的研究工作及取得的研究成果。除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得南京航空航天大学或其他教育机构的学位或证书而使用过的材料。本人授权南京航空航天大学可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。（保密的学位论文在解密后适用本承诺书）作者签名：日期：万方数据 南京航空航天大学硕士学位论文摘要合成孔径雷达（SyntheticApertureRadar，SAR）是一种全天时、全天候的微波成像系统，高分辨率的特点使它在军用和民用领域有着不可替代的作用。随着合成孔径雷达成像技术的发展，各种高分辨率成像算法应运而生。然而高分辨率带来的巨大计算量成为某些成像算法实际应用的瓶颈，其中最为典型的就是后向投影（BackProjection,BP）算法。BP成像算法是一种时域成像算法，与传统SAR成像算法相比，该成像算法原理简单，并且在原理上不存在任何理论近似，能够实现高分辨率SAR成像。因此，更具有实际的研究价值。基于以上背景，本文的主要工作如下：（1）分析了BP成像算法以及快速BP（FastBackProjection，FBP）成像算法的成像模型、实现原理和计算量。针对BP成像算法计算量巨大的特点，本文在BP算法的基础上实现了一种FBP成像算法。实验结果证明该FBP成像算法与BP成像算法成像质量相当，并且FBP算法在一定程度上从算法层面降低了计算量。（2）研究了BP成像算法在非理想航迹下的运动误差以及运动补偿技术。分析了BP成像算法误差的来源，建立运动误差模型，分别提出了基于对比度最优准则的自聚焦方法以及基于划分子孔径和对比度准则的BP成像自聚焦方法，实测数据成像结果验证了这两种运动补偿算法可分别有效应用于短孔径和长孔径成像处理。（3）介绍了基于计算统一设备架构（ComputeUnifiedDeviceArchitecture，简称CUDA）环境下的图形处理器（GraphicProcessingUnit，GPU）编程技术，分析了BP成像算法的内在并行性，提出了一种适合GPU加速实现的BP成像算法加速方案；针对SAR处理数据量较大以及GPU显存受限的问题，在此方案的基础上进一步提出基于流技术的GPU优化方案。实测数据处理结果为优化后比优化前平均成像速度提升约78.8%，表明该方案的有效性和可行性。关键词：合成孔径雷达，后向投影，运动补偿与自聚焦，图形处理器，计算统一设备架构I万方数据 基于后向投影的SAR成像算法与GPU加速研究ABSTRACTSyntheticapertureradar(SAR)isahighresolutionmicrowaveimagingsystemwithallweather,dayandnightcapability.Thecharacteristicofhighresolutionmakesitplayanirreplaceableroleinthemilitaryandcivilfields.WiththedevelopmentofSARimagingtechnology,allkindsofhighresolutionimagingalgorithmsarise.However,tremendouscomputationamountcausedbyhighresolutionbecomesthebottleneckofsomeimagingalgorithmsinpracticalapplication.Andamongthem,BPalgorithmisatypicaltime-domainimagingalgorithm.ComparedwiththetraditionalSARimagingalgorithms,itissimpleinprincipleandhasnotheoreticalapproximation.Therefore,BPalgorithmhasbetterimagingresultandmorepracticalstudyvalue.Basedontheabovebackgrounds,themainworkofthispaperareasfollows:(1)WeanalyzetheBPandFBPimagingalgorithms,includingtheirimagingmodels,realizationprinciplesandcomputationamountanalysis.TheBPalgorithmcan’tbeusedwidelybecauseoftremendouscomputation.ThispaperstudiesaFBPimagingalgorithmbasedonBPalgorithm.AndtheexperimentalresultshowsthattheimagingqualitiesofFBPalgorithmandBPalgorithmarealmostthesame.AndFBPalgorithmdecreasesthecomputationamounttoacertaindegreeatthealgorithmlevel.(2)WealsosudythemotionerrorandmotioncompensationmethodofBPalgorithminthecaseofnonidealtrajectory.TheerrorsourceofBPalgorithmisananlyzedandthemotionerrormodelisestablised.ByIntroducingthemotioncompensationtechnologyofBPalgorithm,weproposetwophasecompensationalgorithmsbasedoncontrastmaximizationcriterionandsubaperture-contrastmaximizationcriterion,respectively.Andthefinalimagingresultsofthemeasureddataverifythatthetwomotioncompensationalgorithmscanbeeffectivelyappliedtorespectivelyshortapertureandlongaperture.(3)WeintroducetheGPUprogrammingtechnologybasedonCUDA,andAnalyzetheinherentparallelismofBPalgorithm.AschemesuitableforGPUaccelerationisalsoproposed.ComparedwiththeSARrawdataamount,GPUmemoryisnotenough.ThenGPUaccelerationoptimizationisproposedbasedonstreamtechnology.Theresultofthemeasureddatashowstheeffectivenessandfeasibilityofthescheme,andtheavergespeedofoptimizedalgorithmimproves78.8percent.Keywords:SyntheticApertureRadar(SAR),BackProjection(BP),MotionCompensationAndAutofocus,GraphicProcessingUnit(GPU),ComputeUnifiedDeviceArchitecture(CUDA)II万方数据 南京航空航天大学硕士学位论文目录第一章绪论............................................................................................................................................11.1研究背景及意义........................................................................................................................11.2合成孔径雷达发展历史............................................................................................................21.3后向投影成像算法发展............................................................................................................31.4SAR运动补偿技术的研究........................................................................................................31.5GPU发展历史............................................................................................................................51.6GPU在SAR成像中的应用......................................................................................................81.7本文主要工作............................................................................................................................9第二章BP成像算法和FBP成像算法..............................................................................................102.1信号模型..................................................................................................................................102.1.1聚束SAR空间模型.......................................................................................................102.1.2回波信号分析.................................................................................................................112.2BP成像算法.............................................................................................................................122.2.1BP算法原理分析............................................................................................................122.2.2BP运算量分析................................................................................................................132.3FBP成像算法...........................................................................................................................142.3.1成像模型.........................................................................................................................142.3.2FBP成像原理.................................................................................................................142.3.3FBP算法流程.................................................................................................................152.3.4FBP算法运算量.............................................................................................................162.4数据仿真..................................................................................................................................172.5本章小结..................................................................................................................................20第三章BP算法运动误差分析以及自聚焦方法研究........................................................................213.1运动误差分析..........................................................................................................................213.1.1运动误差模型.................................................................................................................213.1.2相位误差对图像的影响分析.........................................................................................233.2BP成像算法自聚焦研究.........................................................................................................233.2.1基于对比度最优化准则的自聚焦算法.........................................................................243.2.2基于子孔径划分与对比度准则的BP成像自聚焦方法..............................................263.2.3实测数据处理结果.........................................................................................................28III万方数据 基于后向投影的SAR成像算法与GPU加速研究3.3本章小结..................................................................................................................................30第四章GPU及CUDA介绍..............................................................................................................314.1GPU并行处理技术..................................................................................................................314.2编程模型..................................................................................................................................324.2.1主机与设备.....................................................................................................................324.2.2内核函数.........................................................................................................................344.2.3线程结构.........................................................................................................................354.3存储器模型..............................................................................................................................364.4CUDA程序优化.......................................................................................................................384.4.1存储器访存优化.............................................................................................................394.4.2异步执行.........................................................................................................................404.4.3指令优化.........................................................................................................................414.5本章小结..................................................................................................................................42第五章基于GPU平台的BP成像算法实现....................................................................................435.1BP成像算法的流程分析.........................................................................................................435.1.1距离压缩.........................................................................................................................435.1.2反投影.............................................................................................................................435.2BP成像算法的并行实现.........................................................................................................445.2.1算法设计流程.................................................................................................................455.2.2距离压缩的并行实现.....................................................................................................465.2.3后向投影计算的并行实现.............................................................................................465.3BP成像算法的优化设计.........................................................................................................475.3.1计算精度与速度优化.....................................................................................................475.3.2基于流技术的结构优化.................................................................................................485.3.3反投影计算结构优化.....................................................................................................485.4实测数据结果处理分析..........................................................................................................495.5本章小结..................................................................................................................................52第六章总结与展望...............................................................................................................................53参考文献.................................................................................................................................................55致谢........................................................................................................................................................59在学期间的研究成果及发表的学术论文..............................................................................................60IV万方数据 南京航空航天大学硕士学位论文图表清单图1.1CPU与GPU组成对比图...........................................................................................................6图1.2NVIDIATesla系列GPU图.......................................................................................................7图2.1聚束SAR空间模型图..............................................................................................................11图2.2聚束SAR后向投影示意图......................................................................................................12图2.3FBP成像模型图.......................................................................................................................14图2.4极坐标系下误差图....................................................................................................................15图2.5第i步处理过程.........................................................................................................................16图2.6不同n值下的OPT比值.........................................................................................................17图2.7BP算法点目标仿真图..............................................................................................................18图2.8FBP算法点目标仿真图...........................................................................................................19图2.9BP和FBP实测数据成像.........................................................................................................19图3.1非理想航迹聚束SAR成像空间几何模型图..........................................................................21图3.2极坐标形式BP成像模型.........................................................................................................23图3.3对比度最优化自聚焦算法实现流程图....................................................................................26图3.4分子孔径自聚焦示意图............................................................................................................27图3.5划分子孔径结合对比度最优化自聚焦处理流程图................................................................27图3.6自聚焦前BP算法图像.............................................................................................................28图3.7自聚焦前BP算法图像放大图.................................................................................................28图3.8自聚焦前BP算法图像局部放大图.........................................................................................29图3.9自聚焦后BP算法图像.............................................................................................................29图3.10自聚焦后BP算法图像放大图...............................................................................................30图3.11自聚焦后BP算法图像局部放大图.......................................................................................30图4.1CUDA编程模型.......................................................................................................................33图4.2加法在GPU实现原理图..........................................................................................................34图4.3线程结构模型...........................................................................................................................35图4.4warp执行模型图.......................................................................................................................36图4.5存储器层次结构图....................................................................................................................37图4.6kernel执行方式比较图.............................................................................................................41图5.1后向投影算法的GPU并行计算分配图..................................................................................44V万方数据 基于后向投影的SAR成像算法与GPU加速研究图5.2利用CUDA实现基于CPU+GPU的后向投影成像算法流程图...........................................45图5.3多流异步执行流程图................................................................................................................48图5.4优化之前执行流程图................................................................................................................48图5.5优化之后执行流程图................................................................................................................49图5.6CPU和GPU目标仿真图.........................................................................................................50图5.7线性插值点目标仿真图............................................................................................................50图5.8脉冲压缩耗时图........................................................................................................................50图5.9反投影耗时图............................................................................................................................51图5.10实测数据成像图......................................................................................................................51表2.1SAR仿真数据成像参数...........................................................................................................18表4.1存储器比较................................................................................................................................38表5.1CPU与GPU型号及主频对比.................................................................................................49表5.2SAR实测数据成像参数...........................................................................................................49表5.3基于CPU/GPU下不同插值方法成像时间对比（混合精度）（单位：毫秒）...................51表5.4CPU与GPU成像质量比较.....................................................................................................52VI万方数据 南京航空航天大学硕士学位论文缩略词缩写英文全称注释BPBackProjectin后向投影SARSyntheticApertureRadar合成孔径雷达CUDAComputeUnifiedDeviceArchitecture计算统一设备架构CPUCentralProcessingUnit中央处理器GPUGraphicsProcessingUnit图形处理器CBPConvolutionBackProjection卷积反投影CATComputer-aidedTomography计算机辅助层析FBPFastBackProjection快速反投影FFBPFastFactorizedBackProjection快速分解反投影INSInertialNavigationSystem惯性导航系统GPSGlobalPositioningSystem全球定位系统RCMRangeCellMigration距离单元徙动MDMapDrift子视图偏移法QPEQuadraticPhaseError二次相位误差MAMMulti-ApertureMapDrift多子孔径子视图偏移法PDPhaseDifference相位差分SPESub-AperturePhaseError子孔径相位误差ROPERankOnePhaseEstimate秩一相位估计PGAPhaseGradientAutofocus相位梯度自聚焦MLEMaximumLikelihoodEstimate最大似然估计CRLBCramer-RaoLowerBound克拉美-劳下界PASTProjectionApproximationSubspaceTracking投影近似子空间跟踪RLSRecursiveLeastSquare递归最小二乘GEGeometryEngine几何引擎ALUArithmeticAndLogicUnit逻辑计算单元DRAMDynamicRandomAccessMemory动态随即存储器RDARangeDopplerAlgorithm距离多普勒算法PFAPolarFormatAlgorithm极坐标格式算法VII万方数据 基于后向投影的SAR成像算法与GPU加速研究CSChirpScaling调频算法OPTOperationTime运算时间SIMDSingleInstruction,MultipleData单指令多数据SIMTSingleInstruction,MultipleThread单指令多线程SMStreamingMultiprocessor流多处理器SPStreamingProcessor流处理器FFTFastFourierTransform快速傅立叶变换IFFTInverseFastFourierTransform逆快速傅立叶变换VIII万方数据 南京航空航天大学硕士学位论文第一章绪论1.1研究背景及意义[1][2]合成孔径的基本概念起源于二十世纪五十年代初期，美国古德依尔公司的CarlWiley第一个提出可以通过频率分析来改善雷达方位向分辨率，这是合成孔径雷达（SyntheticAperture[3][4][5][6]Radar，SAR）成像原理发展的最初阶段。SAR的出现极大地丰富了原有的雷达概念，是现代遥感应用领域的一项重要技术，使得雷达可以对感兴趣的目标（如海面舰船、地面坦克、建筑物等）进行高分辨成像与识别，而不再仅仅将目标视作一个点来分析处理，因而能够为人[7][8][9][10]们提供更多有价值的信息。SAR是一种全天时、全天候的微波成像雷达，因此SAR成像技术的发展为基于光学传感[11]器的成像技术提供了非常有用的补充。另外，SAR强透射性、远距离、高分辨率的特点使它在军用和民用领域有着不可替代的作用，对国防技术现代化、国民经济建设和科技的发展具有重大的意义。因此，SAR成像技术也越来越受到国际上先进技术国家的重视，是竞争激烈、发展迅速的技术领域。[12]合成孔径雷达采用合成孔径和脉冲压缩技术对目标进行二维高分辨率成像。平台与目标之间进行相对运动可以得到方位向的高分辨率，发射大带宽线性调频信号，然后对接收到的回波信号进行脉冲压缩处理可以得到距离向的高分辨率。SAR成像算法都是假设雷达平台沿着理想航迹运动，但雷达平台在实际运动中会由于发动机震动、气流等外界因素的影响偏离理想航迹而导致运动误差的形成，最终影响SAR成像效果。因此，进行有效的相位补偿是获得SAR高分辨率成像的关键所在。为了获得高分辨率SAR成像，SAR处理算法复杂度越来越大，需要处理的回波数据量也越来越大，传统的SAR处理模式已经无法满足成像处理的需求。尤其是对于BP成像算法，该算法是一种时域成像算法，其在原理上不存在任何理论近似，能够实现高分辨率成像。但是由于该算法是逐点成像，使得该算法运算量巨大，实现效率非常低。因此，国内外学者一方面致力于快速算法实现问题的研究，力图在算法层面降低运算量；另一方面致力于高速并行器件编程模式的研究，通过硬件加速来提高数据处理的效率。传统的CPU（CentralProcessingUnit，中央处理器）单核和多核（最多8~16个核）运算已经无法满足如此高度复杂的运算需求，采用处理器数目更多的硬件进行并行加速处理已成为一种趋势。图形处理器（GraphicProcessingUnit，GPU）技术的迅速发展为大规模数据的计算与分析提供了一种新的技术手段。GPU强大的处理能力、高度并行性及高显存带宽为高效的SAR成像算法提供了具有发展前景的新型运算平台，为一定程度上提高BP成像算法的实现效率提供了可能。1万方数据 基于后向投影的SAR成像算法与GPU加速研究1.2合成孔径雷达发展历史二十世纪八十年代后，美国有多家公司开始专门研究和制造可以作为武器系统的高性能[13]SAR系统，并将其陆续应用于军事领域成功地装备部队，如洛克希德-马丁公司研制的GSAR系统作用距离可达150公里，成像区宽度可达40公里，在条带模式下分辨率可达3米，在聚束模式下分辨率可达1米。诺思罗普-格鲁曼公司研制的TESAR系统在条带模式下分辨率可达1米，在聚束模式下分辨率可达0.3米，且系统的重量减轻至30.8kg。还有由美国Michigan大学ERIM实验室于1987年为海军研制的P-3SAR系统，在1994年ERIM又在P-3SAR上加装了[14]UHF频段超宽带系统，用来进行穿透植被的研究。P-3SAR可以在X、C、L和UHF四个波段工作，为全极化系统，具有多种成像模式，并可采用偏置相位中心天线（DPCA）技术进行动目标检测。美国国防部高级研究计划署（ARPA）利用该系统研究了叶簇穿透（FOPEN）和地表穿透（GPEN）雷达和干涉雷达检测和识别隐藏目标。其它一些国家也先后开展了机载SAR技术的研究，在1996年EUSAR会议上共有11个国[15][16][17][18]家介绍了18种机载SAR系统，具有代表性的机载SAR系统主要有德国的E-SAR，[19][20]瑞典的Carabas等。德国宇航院（DLR）为了测试新的信号处理技术而研制的E-SAR是一个比较独特的全极化雷达系统，它的搭载平台是DonierDO228型飞机，其机构设计简单，利用小天线、宽波束消除飞行扰动的影响，保证波束能稳定地照射地面区域，并将惯导天线直接固[21]定在机身上。信号处理过程中采用了很多新的技术，例如运动补偿时利用的是反射率偏移法[22]与惯导和全球定位系统（IMU/GPS）所提供数据结合的方法，精度很高，使得该系统的分辨率已达到2m以下。瑞典国防研究机构研制的CARABASII合成孔径雷达系统工作在VHF/UHF波段，也有穿透叶簇和穿透地表的能力。采用步进频率调制波形，在多次飞行试验中成功地对丛林中和地表下的目标成像并能够实时成像，其分辨率达2.5米。我国于1976年组建了国内第一个以微波成像雷达为研究方向的研究室，即中国科学院电子学研究所，开始系统地开展合成孔径雷达的研究，很多高校也先后开始SAR的研究。1979年9月17日，中科院电子所研制的第一台机载SAR样机在陕西武功机场试飞成功，为我国获取了第一张SAR图像。该系统距离分辨率仅为180m，这是因为在距离向并未采用脉冲压缩技术；而在方位向采用了合成孔径技术，因而方位分辨率达到30m，这标志着我国从此成为世界上为数不多的掌握这一核心技术的国家之一。1980年第二台采用了距离脉压技术的机载SAR系统试飞成功，距离和方位分辨率都提高到了15m。接下来单通道SAR系统，多通道、多极化SAR系统又陆续于1983年与1987年分别研制成功。随后，1994年研制成功的机载实时成像处理器意味着我国机载SAR处理成像技术又迈进了一大步。随后中科院电子所又成功研制了我国第一部L波段雷达系统（L-SAR），它具有实时成像的能力，分辨率为3m×3m。2000年10月，由2万方数据 南京航空航天大学硕士学位论文中国航空607研究所与南京航空航天大学联合研制的某型机载SAR试飞成功，该系统工作在X波段，采用全数字化处理方式，在国内首次得到3m×3m的实时机载SAR图像，并首次实现了机载前斜视成像；在对数据的进一步处理中，又在国内首次得到方位分辨率优于0.5m的SAR[23][24]图像。1.3后向投影成像算法发展在高分辨率SAR中，获得高的方位分辨率需要很长的合成孔径，由于各种因素的影响，飞机有很大抖动，很难保证沿着理想的航迹飞行，此时对运动误差的补偿成为成像算法中需要考虑的重要问题。对有较大运动误差的大转角超高分辨率SAR成像，最容易想到的是基于时域处理的成像算法，最经典也是最精确的一种时域成像方法是后向投影算法（BackProjection，简称BP）。后向投影成像算法源于计算机辅助层析（Computer-aidedTomography，简称CAT）成[25]像技术。D.C.Munson等人在1983年从算法原理上揭示了CAT技术和聚束SAR成像之间的相似关系，据此使得人们可以将己经相对成熟的CAT成像技术应用于聚束SAR成像的研究中。简单说来，BP算法是一种逐点成像的算法，是一个点对点的图像重建过程。在实际中，雷达发射的是球面波，那么散射点回波信号在距离压缩后的徙动轨迹是弯曲的，且不同距离散射点轨迹的弯曲程度不一样，因而不同散射点需要进行不同的聚焦处理。而BP逐点成像的特性恰好能满足这个要求，它可以通过计算每个像素到每个天线位置的距离，沿每个散射点的轨迹对其[26][27]进行时域相干叠加实现高分辨率成像。也就是说，在每个脉冲对应的天线位置，都逐点计算各个像素到此天线位置的距离，然后通过对距离压缩后回波数据进行插值的方法得出该脉冲[28]对各个像素点所做出的不同贡献。由于传统BP算法计算量巨大，给实际应用带来了难题，因此对快速BP成像算法的研究也成为了当今SAR成像领域的热点。国内外学者对快速BP算法的研究有显著的成果：1999年，AliF.Yegulalp提出了一种基于BP算法的FBP（FastBackProjection）算法；2003年，LARSM.H.ULANDER提出了一种基于BP算法的FFBP（FastFactorizedBackProjection）算法；2013年，KyraMoon，DavidG.Long提出了一种适用于条带SAR的FBP算法。在国内：2005年，余洪涛等提出了超宽带SAR成像的快速BP算法；2011年，李杨寰等提出了一种任意孔径的快速后向投影成像算法；2013年，李浩林，张磊等提出一种改进的快速分解后向投影SAR成像算法。就国内外目前的研究现状来看，快速BP算法大致分为两个研究方向的思路：一类是划分子孔径，另一类是划分子图像，这两种方法提高效率的基础是利用BP算法中的空间或时间采样的冗余。快速BP成像算法在算法层面降低了运算量，减少了BP成像的时间，使得对BP成像算法的研究进入一个新的阶段。1.4SAR运动补偿技术的研究[29]目前，SAR的运动补偿主要有两种方法：一是基于仪器测量；二是基于信号处理。前者3万方数据 基于后向投影的SAR成像算法与GPU加速研究利用装在SAR载体上的INS（InertialNavigationSystem，惯性导航系统）和GPS（GlobalPositioningSystem，全球定位系统）测得的载体运动参数获得SAR载体的运动轨迹从而实现运动补偿。1982年，J.N.Damoulakis等人通过模拟实验分析出结合惯导测量单元（IMU）和主惯导系统（INS）的导航信息应用于运动补偿，可获得良好的效果。而到了上世纪九十年代，由于全球定位系统（GPS）的定位精度不断得到提高，因此国外关于机载SAR运动补偿技术的[30][31]研究多是基于GPS/INS组合导航系统所获得的雷达平台的准确位置来实现。但是，实际使用时特别是在进行高分辨成像处理的情况下，仪器测得的运动参数往往仍不能满足要求。基于信号处理的运动补偿技术被广泛的应用到SAR成像处理过程中。这一思想源自于合成孔径雷达通过正常飞行录取的回波数据中就已经包含了运动不正常的情况这一事实。本文试图探讨的自聚焦技术就是信号处理里的一个重要研究内容，它能够有效的补偿仪器难以检测的SAR载体快速扰动的影响。因此，自聚焦技术是SAR成像处理中不可或缺的一个流程，它研究的是经运动补偿后残留RCM（RangeCellMigration，距离单元徙动）及方位相位畸变这两个问题，是运动测量单元（如IMU，GPS等）的补充。[32]1972年Hughes公司提出的MD（MapDrift，子视图偏移法）是最早被提出且成功处理了SAR散焦图像，它仅仅估计QPE（QuadraticPhaseError，二次相位误差）系数。MD算法是将全孔径的相位历史域分成两个不重叠的子孔径，然后分别对图像作相关处理获得相对偏移量，利用该相对偏移量与QPE系数成正比的关系获得QPE从而进行补偿。该方法的优点是计算量小，无需迭代，估计精度也高；缺点是它是基于参数的估计方法且只能估计QPE，无法估计高次相位误差。对MD算法作改进一般可以估计不超过5阶的相位误差，称为MAM（Multi-ApertureMapDrift，多子孔径子视图偏移法）。MAM算法将全孔径分成若干个子孔径，利用多个子图像[33]对的相对偏移量可以获得高阶的相位误差。一般地，将全孔径分成N个子孔径则可估计最高阶为N的相位误差。但是，当子孔径数增多时，脉冲响应的分辨率相应地会降低，因此MAM一般能估计不超过阶的相位误差。至1980年代末期，提出了PD（PhaseDifference，相位差分）[34]自聚焦算法。它是将全孔径划分为不重叠且等长的子孔径，利用SPE（Sub-AperturePhaseError，子孔径相位误差）的Fourier变换估计相位误差。PD算法不需要迭代，计算量小于MD[35]算法。1992年提出了ROPE（RankOnePhaseEstimate，秩一相位估计）算法。该算法利用相位的一次差分，通过方位向与距离向的不断交替迭代估计相位误差函数的一次差分和Doppler频率。尽管ROPE算法可以估计任意阶的相位误差，但是距离单元中仅存在一个强散射点的假设限制了实际的使用。[36]1989年提出了PGA（PhaseGradientAutofocus，相位梯度自聚焦）算法，是一种基于模型的估计算法，能够估计任意阶的相位误差，且具有较好的鲁棒性，是相位误差的MLE（MaximumLikelihoodEstimate，最大似然估计），理论上能够达到CRLB（Cramer-RaoLower4万方数据 南京航空航天大学硕士学位论文Bound，克拉美-劳下界）。因此直到今天，PGA算法仍然是使用最广泛的一种自聚焦算法。PGA算法主要分为四个步骤：循环移位；加窗；相位误差估计；迭代校正。1993年提出了基于特征[37]分解的自聚焦算法。该方法利用所有的脉冲信息直接估计相位误差，与PGA算法不同的是它直接估计相位误差，并且估计精度非常高。特征向量法需要估计距离单元样本的协方差矩阵以及对其特征分解，因此巨大的计算量也限制了它的应用。但是相位误差就是特征向量法中协方差矩阵最大特征值对应的特征向量，而这属于主成分分析的问题，可以通过其他方法计算得到。因此，为克服特征分解法的计算量大的问题并充分利用该算法估计精度高的优势，文献[38]提出了基于PAST（ProjectionApproximationSubspaceTracking，投影近似子空间跟踪）的自聚焦算法。它是通过RLS（RecursiveLeastSquare，递归最小二乘）的方法来逼近信号子空间。基于PAST算法的自聚焦算法估计能够有效的减小特征分解法计算量巨大的难题，其算法复杂度与估计性能均与PGA算法相当。国内关于运动补偿的研究起步相对来说要晚一些，不过也在这方面做了许多的工作。文献[39]在‘运7’机载SAR实测数据处理的过程中，运用瞬时多普勒调频率来估计方位向速度的变化量，从而对SAR图像的聚焦效果有明显的提升。文献[40]和[41]在无导航数据的无人机实测数据的处理过程中，在传统的瞬时多普勒调频率的运动补偿方法的基础上进行改善和提高，成功的首次实现了1m*1m的高分辨力的无人机SAR图像。虽然，我国在此领域取得了许多令人满意的结果，但是与国外的先进水平还是存在不少的差距，这些差距正是需要我们不断的去追赶。各种方法按照其性能、实现的复杂程度具有各自的优点和局限性，但它们都建立在这一假设的基础上，即各距离单元上的相位误差近似相同（空不变），并且除PD方法外，均以迭代方式把同一校正矢量用于图像内的所有散射点。上述自聚焦算法只能够适应一定程度上的运动误差，当运动误差过大，平台扰动剧烈的时候，频域成像算法就不如时域成像算法的运动补偿性能。尤其是当雷达平台驾驶员进行有意机动，形成大机动曲线航迹时，很多频域成像算法将失效，此时会考虑时域成像算法来解决这一问题。1.5GPU发展历史上个世纪60年代至70年代，对于图形图像的处理处于起步发展阶段，受到人员研发水平、研发条件与硬件设备等各方面的限制约束，作为计算机系统的对图形数据的一种处理方式，图形处理的规范性制定只局限于软件开发过程中的一些开发准则制定。到了后期，对于图形处理的技术慢慢开始发展起来，不断有新型的技术革新产品出现，如几何引擎(GeometryEngine，GE)的出现，它的出现是图形处理系统发展中的一个里程碑，标志着数字图像处理进入以图形处理器为主要产品的时代。5万方数据 基于后向投影的SAR成像算法与GPU加速研究几何引擎的重要组成部分之一是一个四位的浮点向量计算器，具体由存储器定制码元实现不同的运算，可对图形图像进行变换、滤波、去噪等处理。一个完整的三维图形图像处理系统可以由数量为12或以上的GE子模块组成，从而达到对图形的高度精确运算效果。以GE技术为核心，针对SGI图形图像处理，JamesClark创建了SGI多媒体公司，并拥有多项著名的革新产品，如GL、OpenGL等，为工业图形的处理器生产奠定了行业基础，并为后续的产品发展提供了技术引导。随着工业化的生产需求与技术力量的壮大，GE技术及图形图像的处理技术迅猛发展，与此同时，图形图像的发展慢慢从CPU为核心转向了GPU的方向。图1.1CPU与GPU组成对比图zkq20160222由图1.1可知，CPU主要由控制单元（Control）、缓存（Cache）、逻辑计算单元（ALU以及动态随即存储器（DRAM），CPU的大部分资源都被Cache占据，而GPU的大部分资源则主要是ALU，所以GPU的计算能力更大更强。对于图形数据的处理流程一般分为两步：首先将图像分割成各个小的元素，对于每个图像元素进行几何空间变换、元素光照推算等处理；然后对于图形进行绘制工作，主要是将图形光栅化，以便进行亮度调整、噪声去除、纹理分析、高通或低通图形滤波、图像合成等处理。上世纪末，电子巨头NVIDIA开始介入PC的3D图像视觉市场，并开发出了GeForece256产品，该图像处理系统首次运用了三维的几何变换与光照处理，可以呈现更加高清细致的图像。随着CPU被GPU产品的逐渐取代，市场需求与专家学者的关注性研究都推进了GPU产品的发展更新，平均以6个月左右的时间为周期进行产品更迭。在发展的历程中，先后出现了SGI公司生产的几何引擎系统、NVIDIA公司生产的GeForce系统等里程碑式产品，产品的芯片尺寸也大幅度的缩小，甚至可达到90nm的微小级别，更加利于集成式系统与嵌入式产品的开发，从而使得同样单位的硬件面积能够容纳更加多的逻辑电路，如今的集成强度相对于发展之初已经有了上千倍的提升。与此同时，为了便于电路更新与应用系统升级维护，不少的可汇编图形处理系统得到了提出与开发，代表生产商有ATI与NVIDIA。图1.2为NVIDIA生产的TESLA系列显卡。6万方数据 南京航空航天大学硕士学位论文图1.2NVIDIATesla系列GPU图以GPU为核心的图像图形处理系统在研究与发展中不断的改善其性能，壮大了处理系统的附带功能与可扩展体系结构，并已经有六代的产品革新阶段。在一代发展阶段中，代表的产品有TNT2(NVIDIA)、Voodoo3(3dfx)等，它们的处理架构核心都基本相似，利用过了三角形体系的硬件架构，速度上相较于CPU的产品有较大的提升，能同时进行两个像素的纹理提取或处理。虽然一代的GPU产品优于CPU产品，但是其没有引入T&L架构引擎，所以其在速度上还有冗余时延。zkq20160222第二个阶段的发展中，不少厂家将图形的处理设备加以改进，将T&L做成独立的架构体系，从传统的CPU系统中分离，使得处理速度上又有了飞速的发展，一些OpenGL、Geforce产品都将可支持设备顶点快速变换作为核心技术加以推广。这个阶段中的GPU产品相较于第一阶段的产品有了速度上的进一步提升，但缺乏可汇编的功能，不利于拓展性大型系统的搭建与维护。相比较于二代产品，三代图形处理架构首次将整体架构设计为可编程模式，使得工程师可以对系统进行顶点级操作，并在日后的需求更改或系统运营过程中定制切合实际需求的处理指令与产品结构，从而可以对硬件处理的运算流程进行流处理器式的数据采集与分析。同时，该项技术革新在通用的计算需求中体现了核心的使用价值。虽然三代产品拥有可汇编能力，对于系统中的某些模块进行了参数的可调整性设计，但并未形成完整的可编程体系，还有待开发与改善。第四代的产品是在第三代的产品基础上进行了可汇编性的进一步完善改进，相比较于上一代的处理系统，四代图形处理系统的优点有：像素或顶点级别的操作参数配置选择更加多元化；开发指令增加至上千条，利于汇编中的命令传送；图像纹理处理模式灵活易变；数据查找速度加快；支持浮点数据采集处理，处理数组对象无限制，对纹理数据的运算更加精确等等。第四代的GPU产品在通用型计算性问题中的处理明显优于前几代，不论是在处理的速度上，还是在7万方数据 基于后向投影的SAR成像算法与GPU加速研究处理对象的多元化和精确度上，都是前几代产品没有办法比拟的。同时，四代产品中引入了新的高级语言进行开发，如Cg等，更加利于开发者的使用，而一些难度较大的像素处理也可交由如OpenGL等扩展产品实现。第四代的图像处理产品已经有了较完善的应用功能和处理性能，第五阶段的发展中，将硬件电路的逻辑单元集成度进一步提升，如NVIDIAGeForce6800，其内部可以包含亿数量级的晶体管器件与近二十条的管道体系。并行设计的单元架构使得处理系统的速度又有了飞跃式的增加，并让集成系统的体积更小、实现更加灵活、功能更加完善，除了常规的像素纹理等元素的处理分析，第五代系统还在编程架构中添加了目标的动态分析、函数的存储调用等功能，大大简化了使用者的开发难度与开发成本，便于产品的推广与发展。在第五代产品的基础上，六代产品的技术革新点在于系统架构上的优化与性能提升，如图像展现系统界面的生成传送效率的提高；数字图像处理过程中可处理对象的复杂度进一步增加；加强了图像数据运算中的精确度，提升了图像质量，改善成品的清晰度与可视度。六代中代表产品有NVIDIA公司研发的GeForce7800系列，其核心引擎基于GPU架构，具体型号为CineFX4.0。通过新型架构的引进，该产品在图像处理中的质量与速度上都优于之前产品，可实现基于顶点的对像素着色处理、工业标准的图像检测等。并且在原有五代产品的优势上，7800系列的产品进一步改善了像素色调设定的参数配置，并结合经典的三角处理体系，使得六代产品在传统的图像处理运算中有着更加出色的表现，与此同时，在新的三维视觉建模中，六代图像处理zkq20160222软件也能达到开发应用的需求，完全保障了图像的色彩亮度饱和度上的低失真展现，给人类的视觉更加完美的感官体验，也能满足专业的工业化需求应用。GeForce8800产品是在7800系列上的改善，其中集成的晶体单元进一步增加，制造工艺维持了之前的90nm数量级别体系，系统频率大大提升，在速度和性能上都有着质的飞跃。1.6GPU在SAR成像中的应用近年来GPU技术的迅速发展为大规模数据的计算与分析提供了一种新的技术手段。2007年6月，NVIDIA推出了CUDA。将GPU作为数据并行计算设备，CUDA与GPU结合构成软[42]硬件体系。CUDA不需要借助图形学API（ApplicationProgrammingInterface，应用编程接口），而是采用开发人员熟悉的类C语言，因此不需要学习新的语言，只需要掌握GPU架构方面的知识和并行算法。与前期基于GPU的通用计算的开发方式相比，CUDA的出现改变了原来GPU端代码晦涩难懂的境况，引发了GPU通用计算的革命，加速了GPU的普及。与CPU相比，GPU明显具有成本低、性能高的特点。近年来，尤其是Fermi与Kepler架[43]构的出现，使得GPU通用计算的性能有了飞速的提升。随着GPU计算性能的提升，在科学计算领域得到充分的利用，在此基础上，各种复杂的SAR成像处理算法在GPU上得到加速实8万方数据 南京航空航天大学硕士学位论文[44][45][46]现。文献[47]对高分辨率星载SAR成像中ECSA（ExtendedChirpScaling）算法提出了一种加速方案，并在NVIDIATeslaC2075上面实现了6~8倍的有效加速；中科院电子所针对距离多普勒（RangeDoppler,RD）算法、CS（ChirpScaling）算法及k算法提出了不同的加速[48]方案，并在NVIDIAK20C上面实现了每秒约36M采样点的实时处理。文献[49]对复杂场景的SAR成像进行实时仿真，把计算机的相关技术应用到SAR图像仿真中去，充分利用GPU的可编程图形管线，实现了三维复杂场景的建立和实时阴影的生成。1.7本文主要工作本文主要研究了BP成像算法的原理以及快速BP成像算法的实现，另外还研究了BP成像算法的运动补偿技术和BP成像算法的GPU加速实现。本文各章节的主要内容如下：第一章：阐述论文的研究背景及意义，重点介绍了合成孔径雷达的发展史、BP算法、运动补偿技术以及GPU的发展状况。第二章：首先阐述了BP成像算法的信号模型以及成像原理，分析BP成像算法的计算量，并在此基础上提出一种快速BP实现算法。实验结果证明该快速算法的有效性。第三章：本章主要介绍了BP成像算法在非理想航迹下的运动误差问题。首先分析了BP成像算法误差的来源，建立运动误差模型；在此基础上介绍了BP成像算法的运动补偿技术，并提出了基于对比度最优准则和划分子孔径zkq的相位补偿20160222算法，通过实测数据成像结果分析可以看出该补偿算法的有效性。第四章：主要介绍了基于CUDA环境下的GPU并行编程技术，包括GPU编程模型、存储单元结构、优化技巧。第五章：首先分析基于GPU处理的BP算法内在的并行性，提出了一种针对BP算法的GPU加速实现方案；另外，由于SAR处理数据量较大，GPU显存受限，通过进一步研究GPU设计中的多流异步执行技术、数据传输模式和计算速度与精度，提出了一种GPU优化方案。该方案充分利用BP算法内在的并行性，最大程度地利用GPU设备的计算资源，使得成像速度与GPU非优化方案的实现相比有了近一倍的提升，大大提升了BP成像算法的计算效率。实测数据处理结果表明该方案的有效性和可行性。第六章：本章主要对论文的内容进行总结归纳，并对该文的不足之处提出进一步的研究思路，另外对BP算法的发展以及GPU优化实现的研究进行了展望。9万方数据 基于后向投影的SAR成像算法与GPU加速研究第二章BP成像算法和FBP成像算法随着合成孔径雷达成像技术的发展，各种高分辨率成像算法应运而生。目前，SAR成像处理算法主要分为频域成像算法和时域成像算法两类。频域成像算法主要包括距离—多普勒[50][51]（RangeDoppler，RD）算法、极坐标格式算法（PolarFormatAlgorithm，PFA）、CS（Chirp[52]Scaling）算法、-k算法等。时域成像算法主要包括后向投影算法（BackProjection，简称[53]BP）及基于以BP算法的快速BP算法等。本章将重点介绍时域成像算法—BP成像算法和快速BP成像算法的原理及仿真实现。2.1信号模型合成孔径雷达按照不同的照射角度可以分为正侧视、斜视、大斜视等，依据雷达波束的照射方式可以分为聚束式（Spotlight）SAR、条带式（Stripmap）SAR、扫描式（Scan）SAR等。BP成像算法属于时域成像算法，它通过计算雷达天线相位中心与成像目标点之间的双程延时进行相位补偿，因此该算法适用于各种SAR成像模式。本文主要分析研究正侧视下的聚束式SAR成像。下面主要分析BP成像算法的信号模型。zkq201602222.1.1聚束SAR空间模型聚束式SAR成像是一种常见的、具有超高分辨率的成像模式，其特点是在成像处理时间内雷达波束始终跟踪指向特定的成像区域。与其他SAR成像模式相比，聚束SAR的方位向分辨率不再由其合成天线的长度决定，而是由成像时间内雷达转动的观察角度决定。为了分析简单，我们假设雷达平台沿着平行于X轴的方向，在高度固定为H的Z轴上做匀速直线飞行。假设成像场景中心为坐标原点O，R为雷达飞行航迹与成像场景中心的垂直距离，载机理想速度为v，载机理想飞行高度为H，载机在Y轴上的理想坐标为y，则天线相位中心a的瞬时坐标为xyH,,，点Px,y,z为成像场景中任意一点的坐标，R为雷达天线相位中aaa心与原点的瞬时斜距。如图2.1所示为聚束SAR信号模型图。10万方数据 南京航空航天大学硕士学位论文Z飞行航迹XRHXRa·xyz,,PO'yOaY图2.1聚束SAR空间模型图2.1.2回波信号分析雷达发射线性调频信号，其表达式为2sArectexpjkexpj2f(2.1)cTr其中为快时间变量，A为信号的幅度值，T为脉冲宽度，k为信号线性调频斜率，f为发射rczkq20160222信号载频，BkTr为信号带宽，rect为矩形窗函数，其具体定义为：TrTr1,2rect=(2.2)TTrr0,2令2P=rectexpjkexpj2f(2.3)cTr则发射信号可以表示为sAP(2.4)假设雷达平台在X轴的位置为u，则接收到的回波信号可以表示为：2Runsrn,uP(2.5)nc其中为场景中第n个像素点相对于雷达在u位置时的反射系数，Ru为第n个像素点与nn11万方数据 基于后向投影的SAR成像算法与GPU加速研究2Run雷达天线相位中心的距离，是相应的双程时延。c对雷达回波信号解调去载频为：2Run2c24RunfRucnsrn,urectexpjkexpj(2.6)nTrcc2.2BP成像算法BP算法源于计算机层析成像技术，由McCorkle首先引用到冲激体制SAR成像中。BP成像算法是一种基于时域处理的成像算法，是一个点对点的图像重建过程，其基本思想是计算各方位时刻雷达平台的位置与成像区域内每个像素点的双程延时，再找出不同方位时刻对应的回波信号进行相干累加，最后恢复出每个像素的目标函数。在实际中，雷达发射的是球面波，那么散射点回波信号在距离压缩后的徙动轨迹是弯曲的，且不同距离散射点轨迹的弯曲程度不一样，因而不同散射点需要进行不同的聚焦处理。而BP逐点成像的特性恰好能满足这个要求，它可以通过计算成像网格中每个像素点到每个天线位置的距离，然后沿每个散射点的轨迹对其进行时域相干叠加实现高分辨率成像。2.2.1BP算法原理分析BP成像算法是比较经典的时域雷达成像算法，该算法在原理上不存在任何近似，能够实现高分辨率成像。如图2.2所示为对成像网格区域进行BP成像示意图：成像网格大小为MN，M为网格方位向点数，N为网格距离向点数，v为飞机在飞行方向上的速度。·距离像素点·N·················································方·····················位·························像·························素····················点·········Mv·成像区域图2.2聚束SAR后向投影示意图对回波信号做匹配滤波完成距离向脉冲压缩处理，如下：*s,,usuP(2.7)outr12万方数据 南京航空航天大学硕士学位论文*式中P是P时间反转后的共轭，为快时间卷积。通过时域相干累加得到成像网格点最终图像数据：fxy,suu,expj2fudu(2.8)mnuoutmncmn222xmnyu其中，xy,为成像网格中第mn,个网格点，u为雷达天线位于mnmnc0,u，信号到成像网格点xymn,的双程时延。对于成像网格中的像素点，首先计算每个像素点到当前方位时刻雷达天线相位中心的瞬时斜距，以此计算出雷达回波的双程时延，再通过插值从脉冲压缩后的数据中求出各像素点对应的回波值，并将此值沿着方位时刻进行相干叠加，最后得出像素点回波总能量值，即得出该点的图像。后向投影成像算法的具体实现步骤如下：步骤1：读入SAR回波数据；步骤2：距离向脉冲压缩；步骤3：计算网格点与雷达天线的斜距，并进行插值操作得到子图像；步骤4：相位补偿；步骤5：重复步骤3和4，子图像数据沿方位时刻相干叠加。步骤6：保存叠加数据，即成像结果。2.2.2BP运算量分析BP成像算法主要包括距离向压缩和反投影两个过程，距离压缩过程主要包括复数乘法、复数加法以及复数除法；反投影过程主要包括复数乘法和复数加法。下面计算主要步骤所需的计算量，假设回波原始数据大小为RL（R为距离向大小；L为方位向大小），成像网格大小为MN。RR（1）距离向FFT：LRlog次复数乘法，LRlog次复数加法。2222（2）参考系数点乘运算：RL次复数乘法，RL复数除法。RR（3）距离向IFFT运算：LRlog次复数乘法，LRlog次复数加法。2222（4）反投影运算：LMN次复数乘法运算，LMN复数加法。假设回波数据方位向大小L、成像网格方位向大小M和距离向大小N相等，即LMN==，3反投影中的插值运算次数就高达N，这个计算量是非常巨大的。那么由以上分析可知，BP成13万方数据 基于后向投影的SAR成像算法与GPU加速研究像算法的运算量巨大，而主要运算量集中在反投影操作上。2.3FBP成像算法BP算法具有很好的鲁棒性，但由于算法运算量巨大，限制了该算法在实际SAR信号处理[54]中的应用。假如一个L点长度的合成孔径天线，对像素点MN的区域进行成像，需要计算每个像素点到每个方位时刻天线位置的距离；然后对回波信号进行插值，完成像素矩阵的相干叠加。由算法流程可知，仅相干叠加的次数就有LMN次，这是一个非常大的运算量。由于初始阶段子孔径的角域带宽很小，相应的角域采样率可以很低，较少的角域采样点数就能得到“粗”的角分辨率图像，即没有必要将每个子孔径对应的数据后向投影到成像网格的每个像素点。因此，下面介绍一种经典快速BP算法，实验结果验证了该算法的有效性。2.3.1成像模型由上一节分析可知，对于合成孔径天线的每个天线位置，并不需要把得到的回波数据都反投到所有的图像像素点上。该算法具体原理如图2.3所示，不同的天线位置对应一组不同的同心圆，且在每个圆弧上，后向投影的回波数据是相同的。对于A、B处相邻的天线位置，在如图所示的小扇形波束中，可以近似认为它们有相同的一组同心圆，即具有一组相同的后向投影回波数据。而这组相同的后向投影回波数据，可以认为是与扇形波束的中心距离线（即图中的[55]长虚线）相关的，这是FBP算法成像最基本的原理。当相邻的天线位置变大时，从直观上看，扇形波束宽度必须减小，以使相邻天线所对应的圆弧达到满足条件的要求。图2.3FBP成像模型图2.3.2FBP成像原理在图2.3中，设两个相邻孔径天线位于A、B处，且AB长度为d。各位置到像素点的距离如图2.4所示。14万方数据 南京航空航天大学硕士学位论文P’PrR()rR()··dAB图2.4极坐标系下误差图当天线位置处于A点，且其视角为时，所得到的P点的投影数据是没有误差的。当天线转过角度时，所得到的位于P’处的投影数据偏离了前一投影圆弧，存在误差。根据余弦定理，可以得到：22Rrd2rdcos(2.9)和22Rrd2rdcos(2.10)则得到距离误差：22RRcoscosRRR(2.11)R12cosd/rdr/2当像素点与孔径距离远远大于孔径天线长度，且min1,tan时，上式可以简化为dR(2.12)2综上分析，可以得出以下结论：只要保持天线孔径长度d和扇形波束角宽度的乘积不变，距离误差就可以控制在一定的范围。这就是FBP算法实现的原理。2.3.3FBP算法流程由上节分析可知，FBP算法的处理过程可以看为子孔径叠加以及新波束形成的过程，且新波束的个数随着子孔径长度的增加而增加，同时其波束宽度变窄。我们以L点长度的孔径天线，对像素点为MN的区域成像来阐述FBP成像算法流程。为了提高运算效率，此时坐标系选用极坐标系。假设L点长度的合成孔径天线可以表示为Lll...l(2.13)12k其中，l表示第i步处理过程中，子孔径叠加后其个数减少的倍数。i15万方数据 基于后向投影的SAR成像算法与GPU加速研究为简单起见，可设LMN，ll...ln，则在第i步处理过程中，新子孔径12k的长度变为原子孔长度的n倍，但其数目减少了n倍；同时，每个新孔径位置对应的波束数与原孔径位置对应的波束数相比，其数目增加了n倍，但宽度变为原来波束的1n。由于扇形波束的后向投影数据与其中心距离线相关，因此，可以利用新波束的中心距离线来计算后向投影数据。以叠加后某个新孔径位置对应的一个新波束的中心距离线为例，其计算步骤说明如下：首先对这条中心距离线进行M点采样，将这M个采样点对新孔径位置进行后向投影数据的计算，得到一组后向投影数据。然后将这组后向投影数据与原来n个旧波束的中心距离线上的M个采样点的后向投影数据一一对应进行叠加，得到这条中心距离线最终的后向投影数据。由此可见，每个新波束的后向散射数据与原来n个旧波束后向散射数据是相关的。此外，为了保证得到最后所需要的直角坐标下的MN个像素点，还必须对新波束中心距离线上M个后向投影数据进行倍的插值。图2.5第i步处理过程图2.5说明了第i步处理过程，此时n2。在图2.5中，子孔径的数目减少为原来的一半，而波束的数目增加为原来的二倍，宽度变为原来的二分之一。而新波束的后向投影数据则是通过对合成前相邻两个子孔径所对应的两套波束的后向投影数据进行一系列计算后获得的。按照上述处理完成以后的操作，直至完成最后一步（即第k步）的操作。在完成最后一步操作后，可以将与MN个像素点相符的极坐标点保留下来，将与MN个像素点不相符的极坐标点去掉，这样就可以得到像素点为MN的整个区域图像。2.3.4FBP算法运算量由图2.5分析来说明FBP的总运算量。首先考虑第一步处理情况。在此过程中，原始的子孔径叠加后变为Ln个，每个新的子孔径新形成的波束数为n个，而每个新波束的形成需要原来n个旧波束的数据，波束的采样率为M，为了得到最终所需要的像素点，需要对其进行插值。设插值因子为，则其运算量为2LnnnMnN(2.14)2再考虑第二步的情况。在此过程中，子孔径叠加后其个数变为Ln个，每个新的子孔径2新形成的波束数为n个，而每个新的波束的形成需要原来n个旧波束的数据，波束的采样率仍16万方数据 南京航空航天大学硕士学位论文为M，插值因子仍为，则其运算量为222LnnnMnN(2.15)2由此我们看到每一步的运算量均为nN，直到最后所有的波束形成以获得最终所需要的像素点的高精度图像，此时所有的原始子孔径叠加为一个长的孔径。这时，总的运算量为2nNlogN。n一般地，对于SAR信号处理，8倍插值时就能满足成像精度的要求，即8。由上一节3分析可知，原始BP算法插值的运算量为N。由于运算量和运算时间成正比，因此可以得到FBP算法和BP算法所用运算时间（OPT,OperationTime）的比值：tNBP(2.16)t8lognNFBPn由上式可知，对于定值N，当ne（e为自然对数）时，上式可得到最大比值。但由于n必须为整数，因此上式最大比值出现在n3处。由于BP算法所用的时间只和N有关，因此对于给定的N，当n3时，FBP算法所用的时间是最少的，即此时速度最快。从上式还可以看出，当n2时，若N128，则BP算法比FBP算法快；若N128，则FBP算法比BP算法快，且N越大，FBP算法相对于BP算法的速度就越快。如图2.6所示为n23568、、、、时BP算法和FBP算法所用时间比值图，由图可知当n3的时候，OPT值最大。对于任意n值，随着N的变大，OPT值也变大。10090n=2n=380n=5n=670n=860比值50OPT40302010000.20.40.60.811.21.41.61.82N取值大小4x10图2.6不同n值下的OPT比值2.4数据仿真为了验证BP成像算法和FBP成像算法的成像效果，接下来用点目标回波进行成像，表2.1为本文SAR点目标实验数据的参数。另外，MN4096，网格间隔大小为0.2m×0.2m，17万方数据 基于后向投影的SAR成像算法与GPU加速研究FBP级数n2。表2.1SAR仿真数据成像参数参数名称数值调频斜率78956GHz/s信号带宽1.16GHz距离采样率1.5GHz脉冲重复频率PRF2200Hz波长0.03m雷达运动速度154m/s雷达高度5700m场景中心到孔径的垂直斜距10000m（a）BP算法点目标成像（b）点目标方位向剖面图（c）BP算法点目标等高线图图2.7BP算法点目标仿真图由图2.7可以看出，对理想点目标进行后向投影成像，方位向和距离向都具有较好的成像18万方数据 南京航空航天大学硕士学位论文效果，另外方位向的峰值旁瓣比也达到了-13.3dB。（a）FBP算法点目标成像（b）点目标方位向剖面图（c）FBP算法点目标等高线图图2.8FBP算法点目标仿真图(a)BP成像图(b)FBP成像图图2.9BP和FBP实测数据成像19万方数据 基于后向投影的SAR成像算法与GPU加速研究由图2.8可以看出，对理想点目标进行后向投影成像，方位向和距离向都具有较好的成像效果，另外方位向的峰值旁瓣比也达到了-13dB。另外，图2.9为实测数据仿真图，成像大小为40964096，由仿真结果可知，FBP成像效果相比于BP成像效果略差，但整体效果相当。BP成像耗时80.52分钟，FBP耗时18.25分钟，FBP比BP成像速度提升约5倍。FBP算法虽然运算量减小，但是旁瓣变高，说明FBP算法运算量的减小是有一定代价的。2.5本章小结本章节首先介绍了BP成像算法的信号模型，并在此基础上阐述了BP算法的成像原理以及BP算法实现步骤。另外，通过分析BP成像算法的计算量，得出BP算法计算量巨大的结论。另外，本章节还主要介绍了一种经典FBP算法，包括FBP成像算法原理和FBP算法实现流程，通过FBP算法的实现以在算法层面降低计算量以提高成像效率。通过对比两种算法的计算量以及成像仿真结果，证明了FBP算法的有效性。20万方数据 南京航空航天大学硕士学位论文第三章BP算法运动误差分析以及自聚焦方法研究3.1运动误差分析由上一节内容分析可知，雷达平台在沿着理想航迹运动时，通过BP成像算法对点目标进行成像，方位向和距离向都可以达到很好的聚焦效果。但是在实际的成像过程中，由于雷达平台受大气湍流及载机机械振动的影响往往偏离理想航迹，无法精确获得各方位时刻雷达天线与成像区域各像素点的斜距，从而在信号处理的过程中引入相位误差，进而影响最终的成像效果，因此有必要对后向投影成像算法成出的图像进行运动误差的补偿。3.1.1运动误差模型Z实际航迹X理想航迹RaHXR·(x,y,z)YPO’yaO图3.1非理想航迹聚束SAR成像空间几何模型图如图3.1所示，在没有航迹运动误差的情况下，假设SAR系统波束中心始终指向场景中心O点，载机沿着理想航迹（X轴）方向飞行，设慢时间为t，快时间为τ，T为合成孔径的时间，a载机理想速度为v，载机理想飞行高度为H，载机在Y轴上的理想坐标为y，则天线相位中心a的瞬时坐标为xyH,,，其中xvt，场景中任意目标P点xyz,,。理想航迹情况下天aaa线相位中心到场景中目标P点的距离为：22Rxxr,,xxr(3.1)Paa22式中，ryyHz，但在实际情况下，受大气湍流及振动的影响载机往往偏离a理想航迹飞行，如图3.1中虚线所示航迹。假设在实际情况下，天线相位中心的瞬时坐标为21万方数据 基于后向投影的SAR成像算法与GPU加速研究''xyHaa,,，则天线相位中心到场景中目标P点的距离为：''2'2Rxxr,,xxrRxxr,,R(3.2)Paa0a22'''式中ryyHz，R为理想情况下与实际情况下，天线相位中心到场景中目a标P点的距离之差。由于天线相位中心的位置随着慢时间t沿着X轴变化，对于场景中的目标''点P，其Rxxr,,与Rxxr,,也是随着天线相位中心的X轴坐标变化而变化，所以RPaPa对于给定的场景中一点P而言可以理解为是x的函数，即Rx。由2.1.1节可知，目标点Paa的回波信号为：'2RP2''c24RPfRcPsr,xaPrectexpjkexpj(3.3)Trcc'其中为P点相对于雷达天线在x处的散射系数，R()为天线相位中心与P点的距离。对PaP雷达回波信号解调去载频，做距离向脉冲压缩以后得到雷达回波信号可以表示为：''2RP4fcRPRxaS,xaPsincTkrexpjdxdy(3.4)ccBP算法是将各方位时刻的雷达回波数据反投影至成像区域进行成像，式（3.4）经过反投影后可以写成图像域形式：Smexy,mexpjm(3.5)式中01mN，m为方位向脉冲索引，x、y分别为图像在X方向（方位向）和Y方向（距a离向）索引；()m为整个成像场景各散射点对第m个脉冲的反射系数；m为Rx引ea起的相位误0,1,...,N1，m为第m个脉冲的相位误差，即第m个eeeeae脉冲所反投的图像空不变相位误差。由各方位时刻回波数据形成的层析图相干叠加，得到最后的SAR图像为：Na1Ixy,Sxy,expjm(3.6)mem0在式（3.6）中，最后一项指数项是为了消除由非理想航迹所带来的相位误差。但通常在惯道不精确的实际情况中，偏离理想航迹所引起的Rx是未知的，所以由运动误差所带来的a相位误差在BP算法中是无法消除的，因此需要在BP算法SAR成像过后进行自聚焦处理，消除由运动误差所带来的相位误差。22万方数据 南京航空航天大学硕士学位论文3.1.2相位误差对图像的影响分析由以上分析可知，运动误差会导致各阶相位误差的存在，然而不同阶次的相位误差对图像质量的影响也有所不同。在频域成像算法中，由于存在傅立叶变换，一阶相位误差会引起图像在方位向上整体偏移，高阶相位误差引起图像散焦。在BP成像算法中，在极坐标下对场景成像，如图3.2所示，X轴为孔径位置与航迹方向，O点为场景中心点，把成像场景区域划定为以孔径中心为原点的极坐标网格（扇形区域）。场景中任意点P的极坐标为r,，根据多普勒定义fv2sin，将成像区域坐标映射变换为r,sin，则成像区域就是距离—多普勒域，d在方位向与数据域形成傅立叶变换对，根据傅立叶变换性质可知，一阶相位误差会引起图像在方位向的整体偏移。在单孔径成像时，图像方位向偏移对图像处理不会产生很大问题；但在分子孔径进程成像时，图像偏移会使同一目标点在每个子孔径图像中位置不一样，导致子孔径图像融合成全孔径图像会出现各子孔径图像不能完全对准。在实际成像过程中，雷达平台的抖动导致无法获取精确的位置信息，因而无法进行精确的相位补偿，而导致图像散焦。这就是高阶误差所带来的影响。所以，在分子孔径成像时，需要消除由一阶相位误差带来的图像方位向的偏移以及高阶误差导致的散焦问题。YP(,)r·OrO’X图3.2极坐标形式BP成像模型3.2BP成像算法自聚焦研究自聚焦是一种根据回波数据对相位误差进行估计的算法，在SAR成像处理中占据着非常重要的地位。现有的自聚焦算法主要分为两类：第1类是基于图像特显点的自聚焦算法，第2类是基于图像整体质量的估计算法。对于基于图像特显点的方法，其代表算法为相位梯度自聚焦算法（PhaseGradientAutofocus，PGA）。该算法利用多个强点目标的相位历程来对相位误差进行估计。由于其不基于任何相位误差模型，因此PGA算法对高阶和低阶误差均有很好的补偿能力。后一类基于图像质量的估计算法主要将图像的整体信息作为评价准则，通过调节相位的23万方数据 基于后向投影的SAR成像算法与GPU加速研究估计值使图像质量达到最优来确定补偿的相位。常用的评价准则包括：最小熵，最大对比度，最大锐度等。根据误差模型的不同，该算法又可分为基于参数和不基于参数的两类，分别用于估计低频和高频误差。相对于PGA算法，由于该类算法不基于孤立强点，所以其在低信噪比和低对比度的条件下更具有鲁棒性，能够取得更好的聚焦效果。以上提到的各种自聚焦方法都是基于频域成像算法提出的。频域成像算法认为距离压缩后的信号与最终图像之间存在着傅立叶变换对的关系。因此，传统自聚焦算法的实现基于以下两个条件：（1）包络误差对图像的影响可以忽略；（2）图像与距离相位历史数据具有傅立叶变换对关系。自聚焦就是对存在于距离压缩信号域的数据相位误差进行估计，从而实现图像的重新聚焦。相对于频域算法来说，后向投影算法作为一种时域成像算法，其在成像精度、对复杂轨迹和地表起伏的适应性上有着天然的优势。由于在时域算法中不存在明显的傅立叶变换关系，因此无法将传统自聚焦方法直接应用于BP成像数据。为了解决这一问题，文献[56]和文献[57]对BP成像中的信号性质进行了分析。经研究发现，当图像在极坐标或伪极坐标格式下的时候，傅立叶变换对的关系是近似满足的。这样就可以针对（伪）极坐标下的图像设计相应的自聚焦算法。以上提到的算法均基于傅立叶变换对的关系，而对于不满足该条件下的自聚焦方法却缺乏足够的研究。为了弥补上述不足，本文提出了一种针对BP聚束SAR图像，不依赖于傅立叶变换对关系的时域自聚焦算法。该算法根据反投影原理并根据图像锐度值准则实现了在直角坐标系下的BP图像自聚焦。3.2.1基于对比度最优化准则的自聚焦算法对于BP算法成像，一阶相位误差主要引起图像在方位向的偏移，而自聚焦只能消除高阶相位误差引起的图像散焦。对比度最优化自聚焦算法是以图像对比度最优化为评判标准，判断是否完全补偿掉相位误差。下面来推导与BP算法相结合情况的对比度最优化自聚焦算法。设0,1,...,N1为相位误差补偿向量，根据式（3.5）、式（3.6），则相位a补偿后的SAR图像与单脉冲反投的图像的关系为：Na1Ixy,Smxy,expjm(3.7)m0当时，相位误差正好完全消除，图像完全聚焦。采用图像对比度定义：e2222EIxy,,EIxyIxy,C(3.8)22EIxy,,EIxy2式中，表示方差，E表示均值，Ixy,表示图像中第xy,点像素的强度。对于给定的SAR图像，E计算出来是一常数。因此上式可以进一步化简，得到对比度更简单的等24万方数据 南京航空航天大学硕士学位论文价式子，如下：Nx1Ny14CeqIxy,(3.9)xy00式中，Nx、Ny分布是图像在X方向（方位向）与Y方向（距离向）的像素单元总数。当SAR图像聚焦效果越好，图像的对比度越大；如果图像完全聚焦，则对比度达到最大。而图像聚焦过程与补尝的相位密切相关，相位补偿精度越高，图像的聚焦效果越好。而对比度是相位误差补偿向量的函数，所以，可以转化为对对比度求最大值问题。因此，对上式求()m求导得：2CINx1Ny1Ixy,eq22,Ixy(3.10)mmxy002*因为Ixy,IxyIxy,,，式中“*”表示求复共轭，所以2Ixy,Ixy,*2ReIxy,mm(3.11)2ImI*xyS,xy,expjmm将3.11式带入3.10式化简得：CINx1Ny1eqjm2*4ImeSmxy,Ixy,Ixy,(3.12)mxy00为了求得对比度最大时的相位误差补偿函数，需要使上式左边等于0，所以可以令Im[]=0，求得：*jmAme(3.13)Am式中：Nx1Ny12*AmSmxy,Ixy,Ixy,(3.14)xy00算法实现步骤：（1）输入单个脉冲所反投的复图像Sxy,，然后初试化相位补偿向量，求出Ixy,，计m算对比度。（2）估计相位误差，将各脉冲反投的图像域乘以相位误差。（3）计算图像对比度，算出对比度增加的比例，与设定的终止门限T比较，判断是否进入0下一次迭代。如图3.3中流程图所示。终止门限T越小，则迭代次数越多，计算量越大，025万方数据 基于后向投影的SAR成像算法与GPU加速研究所得到的相位补偿向量越精确；反之，T越大，则迭代次数越少，计算量越少，所得0到的相位补偿向量就不够精确，所以要权衡计算量与所要的精确度来确定T，一般取0-310。（4）总的相位误差补偿向量等于每次迭代所得误差向量的乘积，把最终得到的相位误差补偿向量乘到原每个脉冲所得的图像域上去，再把所有脉冲所成的图像相干叠加，便得到了最终的SAR图像。输入各脉冲图像域初始化迭代次数ll=l+1jm()计算elll1j()mS(,)xyS(,)xyem1ml计算对比度lCI()Nll1CI()CI()T0CI()Ylj()mj()meelii0图3.3对比度最优化自聚焦算法实现流程图3.2.2基于子孔径划分与对比度准则的BP成像自聚焦方法为了获得超高分辨率的SAR图像，将增加合成孔径的长度。由于雷达平台并非沿着理想航迹运动，如果进行全孔径成像将会导致误差积累，使成像效果恶化。因此，大孔径情况下，一般先将全孔径划分为若干个子孔径，将子孔径分别BP成像、对比度最优化自聚焦处理，然后再将子孔径图像相干融合。针对运动误差所导致的相位误差，在子孔径分别成像时，一阶相位误差主要引起图像在方26万方数据 南京航空航天大学硕士学位论文位向的整体偏移而不是散焦，基于对比度准则的自聚焦无法消除一阶相位误差，只能消除高阶相位误差。对子孔径图像进行对比度自聚焦处理只是补偿掉子孔径内的高阶相位误差，但各子孔径图像都有不同程度的方位向偏移，所以子孔径之间仍然存在相位误差。子孔径1子孔径2子孔径3子孔径4X图3.4分子孔径自聚焦示意图划分子孔径成像时，如图3.4所示，把全孔径分为四个子孔径，中间细虚线带有箭头指向的X轴为理想孔径，弯曲的虚线是实际孔径，四个子孔径分别成像作对比度最优化自聚焦处理后，此时可以认为子孔径内已经消除了运动误差，子孔径可以认为如图3.4中粗黑线所示，是与理想航迹方向相同的直线。但是，由于四个子孔径的图像分别做对比度自聚焦，所以子孔径与子孔径之间会存在相位误差，如图3.4中四段粗黑线所示，并且子孔径图像间存在方位向偏移。因此，在做完所有子孔径BP算法成像、对比度自聚焦处理之后，需要子孔径图像配准消除子孔径一阶相位误差，然后子孔径之间对比度自聚焦处理。为了在子孔径图像融合是更好的消除一阶相位误差，子孔径BP算法成像采用极坐标形式，其整个成像流程图如图3.5所示。划分子孔径子孔径图像配准子孔径极坐标成像所有子孔径图像间对比度最优化自聚焦对比度最优化自聚焦子孔径图像相干融合成全孔径SAR图像N子孔径是否处理完Y极坐标图像变为直角坐标图像各子孔径图像投影到以孔径中心为原点的极坐标网格上图3.5划分子孔径结合对比度最优化自聚焦处理流程图27万方数据 基于后向投影的SAR成像算法与GPU加速研究3.2.3实测数据处理结果本文中验证的实测数据来自某型机载雷达所录取，通过实测数据的验证，可以有效的证明了最优对比度自聚焦算法与BP算法的结合，有效的提高了SAR图像的质量。该雷达工作在X波段下，以聚束模式采集数据，雷达回波参数如表2.1所示。图3.6自聚焦前BP算法图像图3.7自聚焦前BP算法图像放大图28万方数据 南京航空航天大学硕士学位论文图3.8自聚焦前BP算法图像局部放大图图3.9自聚焦后BP算法图像29万方数据 基于后向投影的SAR成像算法与GPU加速研究图3.10自聚焦后BP算法图像放大图图3.11自聚焦后BP算法图像局部放大图图3.7、图3.10均为全孔径16K成像，成像网格大小为为8K*8K，方位向和距离向网格间隔均为0.1米。图3.7是对比度最优化自聚焦前的SAR图像，图3.10是对比度最优化自聚焦后的SAR图像，图3.8、图3.9和图3.11、图3.12分别为分别为图3.7、图3.10相同部分的局部放大和再次放大图，可以从以上六幅SAR图像比较得出，采用了对比度最优化自聚焦后的图像质量明显好于聚焦之前的图像，验证了对比度最优化自聚焦算法与后向反投影算法相结合的有效性与可行性。3.3本章小结本章节提出了一种适用于后向反投影成像算法的对比度最优化自聚焦方法，该方法结合BP算法成像特点，利用图像对比度最优化来估计其相位误差；并将该方法与多子孔径BP成像进行结合，消除相位误差，提高图像分辨率。通过实测数据验证和比较，证明了该方法估计图像高阶相位误差的精确性，以及能与多子孔径BP成像很好的结合，对BP算法处理后的图像质量有很大的提高，表明该自聚焦方法的有效行和可行性。30万方数据 南京航空航天大学硕士学位论文第四章GPU及CUDA介绍并行性是指在同一时刻或者同一时间间隔内完成两种或者两种以上性质相同或者不相同的工作，只要在时间上互相重叠，都存在并行性。所谓并行处理技术就是分析问题中存在的并行性进行，利用多个处理单元同时执行多个任务，并最终解决问题。按并行化程度的不同，我们通常将并行处理技术细分为以下四种：比特级并行处理、指令级并行处理、数据级并行处理以及线程级并行处理。而GPU技术的高速发展大大提升了线程级并行运算的性能。4.1GPU并行处理技术随着SAR成像技术的飞速发展，人们期望SAR具有更加多样化的工作模式，更高的分辨率，更好的成像效果。这需要增加成像算法的复杂度，但是人们又不希望耗费太长的运行时间。因此，探索新的运算平台，同时研究与之相适应的更高效的SAR成像算法是大势所趋。GPU的产生为高效SAR成像算法提供了一个新的具有前景的运算平台。GPU通用计算性能高同时成本又低，与传统的基于CPU的成像算法相比，基于GPU的经过CUDA优化处理实现的SAR成像算法效率大大提升。关于GPU编程，最初开发人员采用汇编语言来实现，但是这很难掌握。随着GPU的高速发展以及CUDA的产生，这一问题得到了很好的解决。CUDA对图形硬件单元和API进行了封装，然后在类似于CPU的编程环境下对GPU进行编程开发。CUDA程序优化的目的是在误差允许范围内尽快完成计算任务。CUDA编程要达到良好的效果主要考虑计算量、精度和延迟三个方面。计算量有绝对和相对两种方式，对于计算量太小的程序，使用CUDA来进行处理不是很合算。从绝对计算量来说，如果要优化的程序执行时间短，并且调用次数不多时，使用CUDA优化无法将运算效率提高很多。在计算量非常小的情况下，使用CUDA在GPU上进行处理的执行时间的减少量相比于访存和数据传输的延迟时间来说可以忽略不计。我们知道GPU的显存带宽远大于CPU，通过PCI-E连接，所以IO带宽决定了输入和输出的吞吐量。当程序计算密度很低，大部分时间用于IO消耗而不是执行计算时，PCI-E的带宽就会影响整个程序的性能。相对计算量是指并行运行部分占整个程序的时间，当整个程序中大部分执行时间用于串行时，总执行时间并未降低很多。双精度可以帮助我们在大数据计算时获得值得信赖的结果。为了提高CUDA程序精度，常在重要环节采用双精度。目前CUDA还不可以独自为某个处理核心分配任务。CUDA处理过程中，首先将缓冲数据传输到GPU，然后在GPU上进行并行计算，最后将处理后数据传输回CPU，比起在CPU进行串行处理，这个过程需要耗费更长的时间。31万方数据 基于后向投影的SAR成像算法与GPU加速研究4.2编程模型在本章节主要介绍了CUDA的编程模型，包括主机、设备的概念，接着阐述了内核函数的概念，给出了GPU设备端的线程组织模型，并介绍了CUDA编程模型在硬件上的映射。4.2.1主机与设备在CUDA编程模型下，CPU和GPU各自执行不同的任务，CPU通常被称为主机（Host），GPU则通常被称为设备（Device），在一个系统中一般会有一个Host和一个或者多个Device。在CUDA并行模型下，CPU与GPU协同工作，CPU主要负责串行运算，而GPU则主要负责程序的并行运算部分。CPU和GPU各自拥有相互独立的存储地址空间：主机端内存和设备端显存，操作显存需要调用CUDAAPI中的存储器管理函数，包括显存开辟、初始化和释放操作，以及在主机端和设备端之间进行数据传输等。将待设计的程序分为串行部分和并行部分，串行部分交给CPU处理，并行部分交给GPU处理。通常将运行在GPU端的函数称为内核函数，即kernel，一个kernel函数不是一个完整的程序，而是整个CUDA程序中的一个在GPU端并行执行的部分。在整个程序设计的过程中，首先要将任务分解，提取程序中需要并行计算的部分，然后按照CUDA编程规则将待并行部分在GPU端的kernel函数中进行编程实现。串行部分程序只是负责整个程序的串行执行，主要包括数据的传输控制、API函数的调用控制以及kernel函数的控制等等。在整个程序运算过程中，CPU和GPU各司其职，CPU主要负责整个程序中的串行执行部分，GPU则主要负责整个程序中的并行执行部分。一个CPU可以配置多个GPU，GPU通过PCI-E总线与CPU连接，CPU用于计算逻辑性较强的任务，GPU用于计算规模大、分支简单、密度高的并行任务。CUDA的核心是kernel函数，它是并行执行的。在kernel函数启动前，CPU串行代码进行数据准备和设备初始化的工作，然后GPU完成核函数中的并行计算。而kernel函数间的串行计算也依靠CPU串行代码完成。CUDA编程模型如图4.1所示。32万方数据 南京航空航天大学硕士学位论文C程序串行执行thread0thread1thread2thread3thread4thread5thread6thread7串行代码Grid0Block0Block1Block2Block3并行kernelBlock4Block5Block6Block7a„+++++++++串行代码b„c„Grid1Block0Block1Block2Block3并行kernelBlock4Block5Block6Block7图4.1CUDA编程模型由于所有线程执行同一个kernel函数，kernel函数在编写过程中需要能够区分不同的线程实例。CUDA线程管理采用层级结构，自顶向下依次为：线程格（Grid）、线程块（Block）和线程（thread）；Grid与GPU相对应，即一个GPU只有一个Grid；Block与GPU中的SM单元想对应，同一个Block中的线程在同一个SM单元上执行，但同一个SM单元可以串行执行多个Block；thread与SM中的CUDA核相对应，一个线程对应一个CUDA核，同一个CUDA核上可以通过时分方式执行多个thread。33万方数据 基于后向投影的SAR成像算法与GPU加速研究4.2.2内核函数上一节提到，运行在GPU端的部分称为内核函数，即kernel。内核函数必须通过__global__函数类型限定符定义，并且只能在主机端代码中调用。在调用时，必须声明内核函数的执行参数，例如//kernel定义__global__voidAdd(flaot*A,float*B,float*C){inttid=threadIdx.x;C[tid]=A[tid]+B[tid];}intmain(){//kernel调用Add<<<1,N>>>(A,B,C);}分析上一段代码，Add<<<1,N>>>(A,B,C)完成了对kernel函数Add的调用，其中<<<>>>运算符中是kernel函数的执行参数，用于说明kernel函数中的线程数量以及线程的组织形式。1和N是函数的执行参数，表示kernel的Grid中只有一个block，每个block中有N个thread；而A，B，C则是函数执行所需的参数。在设备端运行的线程之间是并行执行的，每个线程都要按照指令的顺序串行执行一次kernel函数。每个线程都有自己唯一个blockID和threadID，blockID和threadID是kernel中的内建变量，不需要由程序员定义，是设备中专用寄存器提供的。所以，内建变量是只读的，且只能在GPU端的kernel函数中使用。如上面例子所示，N个线程中的每个线程都独立执行了一次Add()函数，将对应的两个数进行求和运算，N个线程一共完成了N对元素的求和运算。如图4.2所示，向量A和向量B对应元素相加并将结果对应存储到向量C中。图4.2加法在GPU实现原理图34万方数据 南京航空航天大学硕士学位论文4.2.3线程结构如图4.3所示，kernel函数内部的并行结构分为两层：Grid中的线程块（block）间并行和block中的线程（thread）间并行。block间无法通信，而同一个block内的thread可以相互合作，通过同步执行来协调对存储器的访问。Grid0SPSPSPSPSPSP…Block0…SPSPBlock1SMBlock2SMSM图4.3线程结构模型CUDA将计算任务映射为大量并行thread，然后通过硬件动态调度和执行这些线程。一个block中的所有thread都必须位于同一个处理器核心中，因此处理器核心的存储器资源将制约每个block的thread数量。执行kernel的thread总数等于每个block的thread数乘以block的数量，这些block组织为线程块网格。一个Grid的block数量和由处理的数据总数决定，一般来说，数据总数远远超过处理器数量。CUDA中使用索引变量threadIdx和blockIdx来标识线程，构成一维、二维、或三维线程块。GPU计算单元主要由流多处理器（简称SM）构成，其中又包含8个标量流处理器（简称SP）。SM中的8个SP使用同一块共享存储器。同一block中的thread放在同一个SM中发射，共享数据，其中每一个thread被发射到每一个SP上执行。同一时刻一个SM中可以有多个活动block在等待执行，这样就会造成访存延迟。为了提高资源利用率，考虑在一个SM中发射多个block，这样一个block就可以利用另一个block进行访存或者同步的间隙执行运算。实际运行中，block被分割为warp，它是一条32路SIMD（SingleInstruction，MultipleData，单指令多数据)指令，是SM调度的基本单位，在CUDA编程模型中并不存在，但不能忽略其影响。实际上，访存基本单元是half-warp，它可以由warp中的前半部分线程组成，也可以由后半部分线程组成。warp中包含32条线程是因为每发射一条warp指令，SM中的8个SP将执行4遍。warp执行过程如图4.4所示。35万方数据 基于后向投影的SAR成像算法与GPU加速研究图4.4warp执行模型图[58]CUDA采用了SIMT（SingleInstruction,MultipleThread，单指令多线程）执行模型，在SIMT中，无法知道执行数据的宽度，硬件根据不同的执行宽度而变化，方便了编程。SIMT中每个线程都是私有的，线程之间只能通过共享存储器或者同步进行通信。4.3存储器模型图4.5给出了存储器层次结构图，thread在执行时会访问不同的存储器空间，grid中所有的thread都可以访问同一块全局存储器，每一个block拥有一块共享存储器，每一个thread都拥有自己的私有寄存器和局部存储器，所有thread可以访问常数存储器和纹理存储器。下面将详细介绍这6种存储器的基本使用方法。36万方数据 南京航空航天大学硕士学位论文图4.5存储器层次结构图寄存器（register）的基本单元是寄存器文件，大小为32bit。register是一种高速缓存器，数据传输速度很快。由于要平分给每个线程，所以每个线程只拥有有限的register。当register使用完后，数据就存储于局部存储器（localmemory）中，localmemory访问速度较慢，可以用来存储线程的输入变量和中间变量。全局存储器（globalmemory）存在于显存中，拥有很高的带宽，但由于其没有缓存导致了较高的访存延迟。globalmemory可以被CPU、GPU以及grid中的thread读写访问。要想提高全局存储器带宽必须遵守合并访问要求，并避免分区冲突。类似于寄存器，共享存储器（sharedmemory）也是一种高速可读写存储器，利用sharedmemory可以最快实现线程间通信。常数存储器（constantmemory）是一种只读存储器，同样位于显存。其占用空间小、含缓存的特性适合于存储需要频繁访问的情况。由于是只读的，所以不存在缓存一致性问题。纹理存储器（texturememory）是一种只读存储器，通过缓存加速访问。texturememory能够申请很大的存储空间，常用于图像处理和查表计算，还可以用于对大数据的随机访问或非对齐访问中。表4.1对以上存储进行了比较：37万方数据 基于后向投影的SAR成像算法与GPU加速研究表4.1存储器比较存储器位置拥有缓存否访问权限registerGPU片内(SM内)GPU片内的高速存储器device可读/写sharedmemoryGPU片内(SM内)GPU片内的高速存储器device可读/写texturememory板载显存有device可读，host可读/写constantmemory板载显存有device可读，host可读/写localmemory板载显存无device可读/写globalmemory板载显存无device可读/写，host可读/写另外，在Fermi架构中，每个SM除了拥有专用的纹理内存外，还拥有64KB容量的片上缓存，这部分缓存可配置为16KB的一级缓存（L1）+48KB的共享缓存，或者48KB的一级缓存（L1）+16KB的共享缓存。这种划分方式是完全动态执行的，一个时钟周期之后可自动根据任务需要即时切换而不需要程序主动干预。不管哪种配置，当寄存器变量不够用的时候，L1缓存可以充当临时寄存器的作用，而以前的GPU架构就会把溢出的寄存器变量存放到DRAM（在这里应该是localmemory，即本地存储器（局部存储器），虽然是SM独有，但是没有缓存）中，增加了延迟。对于图形渲染来说，重复或者固定的数据比较多，因此一般是划分48KB为共享内存，当然剩下的16KB一级缓存也不是完全没用，它可以充当寄存器溢出的缓冲区，让寄存器能够实现不俗的提升性能。而在并行计算中，一级缓存与共享缓存同样重要，它们可以让同一个线程块中的线程能够互相协作，从而促进了片上数据广泛的重复利用并减少片外的通信量。共享存储器是许多高性能CUDA应用程序成为可能的重要促成因素。Fermi架构除了给每个SM增加了可配置的L1缓存，还给全局存储器增加了统一的L2缓存。Fermi架构拥有768KB的统一的L2缓存，可以供存、取和纹理操作使用。L2与所有的SM相通，为全局存储器提供高效的、快速的数据传输。4.4CUDA程序优化CUDA程序优化的最终目的是以最短的时间，在允许的误差范围内完成给定的计算任务。在这里，“最短的时间”是指整个程序的运行时间，更侧重计算的吞吐量，而不是单个数据的延迟。因此，一个优秀的CUDA程序，需要不断地优化以提高程序的综合性能。38万方数据 南京航空航天大学硕士学位论文4.4.1存储器访存优化通常，内存访问的带宽远远小于处理器的计算能力，在4.3节中，简单介绍了各种存储器的基本情况，本节会对程序设计和优化中遇到的内存问题做一些详细的讨论。CUDA程序中每个warp只有两种指令：访存和运算。运算指令只有在访存结束，数据准备好以后才能被执行。因此，为了使编写的运算程序具有很高的性能，我们从一开始就要规划好对存储器的访问，同时在开发的过程中每时每刻注意存储器访问对整个程序性能的影响。支持CUDA的显卡在有限的面积上实现了可观的显存带宽，但为此也付出了一定的代价：只有满足合并访问条件同时又有足够的运算指令隐藏访存延迟时，才能使实际带宽提高。以运算为主的CUDA程序中，我们应该尽量避免让存储器访问限制程序性能；而在以存储器访问为主的CUDA程序中，我们应该尽量增大程序的可用带宽。全局存储器访存优化在CUDA的存储器结构中，全局存储器的访存延迟是最大的，同时其存储器容量也是最大的，通过cudaMalloc函数分配的指针来指向其存储器空间。最大的访存延迟和存储器容量使得全局存储器的使用严重影响着CUDA程序的性能。CUDA提供了一种称为合并访问的机制来提高访问全局存储器的效率。所谓合并访问就是尽量要求相邻的线程访问相邻的地址空间，在不同计算能力的设备上，其具体要求还是有所差别的：在计算能力为1.0和1.1的设备上，合并访问要求第k个线程必须对齐访问第k个存储地址，不可以交叉；在计算能力为1.2以上的设备上，合并访问没有这个要求。合并访问的带宽会随着数据类型的不同而不同，正常情况下，使用4字节的数据类型性能较好，但是并不绝对，在某些情况下，使用８字节或16字节的数据性能可能会更好，尤其是在满足1.0和1.1设备的合并访问条件下。在设计的时候要考虑到全局存储器的合并访问，并以此为要求来设计，比如在必要的时候对数据的组织方式进行适当的改变。共享存储器访存优化共享存储器位于GPU片内，比起local/globalmemory来其速度快了很多。在不发生bankconflict的情况下，sharedmemory的延迟几乎只有local或globalmemory的百分之一，访问速度与register相当。对于计算能力1.x的设备，sharedmemory被均匀划分为16个存储器模块（bank），以此来提高并行访问时获得的带宽。不同的存储器模块可以同时工作，互不干扰，同时访问16个bank上的16个地址可以使有效带宽变为原来的16倍。每个bank的固定宽度是32bit，在每个时钟周期内每个bank提供32bit的带宽。warp中前half-warp与half-warp间不会发生bankconflict，39万方数据 基于后向投影的SAR成像算法与GPU加速研究只有处于同一half-warp内的线程才可能发生bankconflict，当一个bank中出现多个被half-warp请求访问的地址时，就产生了bankconflict。硬件将存在conflict的访存请求划分成几次不存在conflict的独立请求，然后串行完成，此时的有效带宽会降低，降低的大小就是独立请求的个数。这里还存在一种特殊情况，当一个half-warp中的所有thread对同一地址请求访问时，只需要一次就可以响应所有thread的请求。常数存储器访存优化常数存储器主要用于存放指令中的常数。对于当前硬件来说，如果一个half-warp的线程访问常数存储器中相同的一个数据，可能只需要一个周期就可以获得这个数据。实际使用constantmemory时速度一般还是低于立即数或者sharedmemory，但还是明显高于将数据存放在显存中的情况。如果有因为尺寸太大无法存放在寄存器或者sharedmemory中的查找表一类的常数数组，可以将其放在常数存储器中获得一定的加速。纹理存储器访存优化纹理存储器能够从过缓存利用数据的局部性，提高效率。它的主要用途是用于存放图像和查找表。使用texture的好处有：1)不用严格遵守合并访问条件，也能获得很高带宽。2)对于随机访问，如果要访问的数据不是很多，效率也不会特别差。3)可以使用线性滤波和自动类型转换等功能调用硬件的不可编程计算资源，而不必占用可编程计算单元。4.4.2异步执行异步执行是指CPU进行的操作（内核启动或者异步存储器拷贝数据）从主机端启动后，在GPU真正完成这些操作之前，在主机端就可以得到这些函数的返回值，CPU线程就可以继续进行下一步操作。也就是说，CPU端的API函数和内核启动的结束，和GPU端真正完成这些操作是异步的。通过异步函数，CPU可以在GPU端进行运算或者数据传输的同时进行其他操作，更加有效地利用系统中的计算资源。内核启动和显存内的数据拷贝（DevicetoDevice）总是异步的，而内存和显存之间的数据拷贝函数则有异步和同步两个版本。例如，下面的代码中进行了一次同步的数据拷贝，只有确实完成CPU和GPU间数据传输后，cudaMemcpy()函数才会返回，此时CPU才能开始执行后面的cpuFunction()函数。cudaMemcpy(a_d,a_h,size,cudaMemcpyHostToDevice);cpuFunction();而调用异步的cudaMemcpyAsync()函数时，CPU线程通过存储器管理API函数启动了一次40万方数据 南京航空航天大学硕士学位论文数据传输，在GPU传输的同时，cudaMemcpyAsync()函数就已经返回了。此时，cpuFunction()函数实质上是与GPU和CPU间的数据传输同时进行的。cudaMemcpyAsync(a_d,a_h,size,cudaMemcpyHostToDevice);cpuFunction();对主机端来说，内核启动是异步的。在下面的代码中，cpuFunction()也会在kernel函数在GPU上运行结束前就开始执行：kernel<<>>(a_d);cpuFunction();如果主机端需要使用内核函数计算得到的结果，应该加入同步保证操作的顺序一致性。kernel<<>>(a_d);cudaThreadSynchronize();cpuFunction();属于同一个流中的内核启动函数总是同步的，使用多个流和异步来实现GPU和CPU同时进行运算。创建一系列流，每个流都进行异步传输和内核执行。此时，一个流的传输和另一个流的执行可以同时进行，因此GPU就能够持续地进行计算，而不必等待数据，如图4.6所示。数据拷贝内核执行(a)kernel串行执行数据拷贝内核执行(b)kernel并行执行图4.6kernel执行方式比较图从图4.6可知，不使用异步执行时的程序运行时间为：执行时间+传输时间，而使用流和异步后，时间将为：执行时间+传输时间/流的数量。此时，CPU和GPU间的数据传输时间被有效的隐藏了，程序的性能得到了改善。4.4.3指令优化SM处理一条warp指令，首先为warp里的每个线程读取指令操作数，执行指令，最后为warp里的每个线程写入计算结果。因此，有效的指令吞吐量不仅取决于名义上的指令吞吐量，还取决于内存延迟和带宽。指令优化主要包括算数指令优化、控制流指令优化、访存指令优化、同步指令优化等。因此，为了提高程序运行性能，一般要采取以下手段增大指令吞吐量：41万方数据 基于后向投影的SAR成像算法与GPU加速研究1)避免使用低吞吐量指令。2)对每种类型的存储器进行优化，有效利用带宽。3)允许线程调度单元尽量用多的数学计算来覆盖访存延迟，这就需要有高的算数密度（或者说，对于每一个内存操作都有大量的算数操作来覆盖），同时每个超多处理器有很多的活动线程。吞吐量是指每个多处理器在一个时钟周期下执行的操作数目。对于大小为32的warp，一条指令由32个操作构成。因此，如果记T为每个时钟下的操作数目，那么指令吞吐量就是每32/T个时钟周期一条指令。所有的吞吐量都是针对一个多处理器而言的。所以，要计算整个设备的吞吐量需要诚意设备的多处理器个数。对程序中需要多次运行的代码进行指令级优化可以有效地提高指令流的吞吐量。4.5本章小结本章介绍了关于GPU的基本概念，包括编程模型、存储器模型以及CUDA程序优化。关于编程模型，主要介绍了主机和设备，GPU端的线程组织结构以及硬件映射。关于存储器层次结构，主要介绍了6种常用的存储器：寄存器、局部存储器、共享存储器、全局存储器、常数存储器、纹理存储器。CUDA程序优化，主要介绍了全局存储器、共享存储器、常数存储器、纹理存储器的访存优化以及指令优化等。GPU为高效SAR成像算法提供了具有发展前景的新型运算平台。与CPU相比，GPU通用计算具有成本低，性能高的特点，这为应对SAR成像信号处理领域新的挑战提供了具有前景的研究方向。42万方数据 南京航空航天大学硕士学位论文第五章基于GPU平台的BP成像算法实现由第一、二章分析可知，BP成像算法与传统的频域成像算法相比原理简单，因在成像过程中不存在任何近似，因此成像效果更好，但是其计算量大的缺点成为SAR成像处理领域的瓶颈。随着高性能硬件技术的发展，尤其是GPU并行技术的出现，使得BP算法再次成为SAR工作者研究的热点。本章从数据处理入手，针对BP成像算法的GPU并行化进行研究分析，并在此基础上提出一种基于GPU的BP成像算法快速实现方案；并通过研究GPU的优化方法，将该方案进一步优化处理，获得了更高的加速性能。5.1BP成像算法的流程分析BP算法主要包括距离向压缩和反投影两个部分，反投影又包括插值和相位补偿。下面对距离压缩和反投影进行详细的分析。5.1.1距离压缩该算法的距离压缩部分与传统SAR成像处理算法的距离压缩原理相同，都是将回波信号与参考信号进行卷积，实现匹配滤波。在某一方位时刻，距离压缩过程可以表示如下：sn()sn()hn()(5.1)r式中，sn()为该方位时刻的回波信号，hn()为参考信号。在信号处理中，时域的卷积等效于频域的相乘，因此在SAR处理中将sn()和hn()先通过FFT（FastFourierTransform，快速傅立叶变换）到频域，二者相乘之后再IFFT（InverseFastFourierTransform，逆快速傅立叶变换）变回到时域即可完成卷积运算。在SAR信号处理时，距离向FFT运算过程在各个脉冲之间是不相关的，因此距离压缩处理的时候脉冲之间的运算可以并行操作。5.1.2反投影BP成像算法假设发射波是冲激球面，通过回波在时域的相干叠加实现高分辨率成像。由2.2.1节可知，慢时间域相干叠加的成像结果可表示为fxm,ymusoutmnuu,expj2fcmnudu(5.2)其中s,,usuP(5.3)outr*为脉冲压缩的结果，式中P是P时间反转后的共轭，为快时间卷积，x,y为mn43万方数据 基于后向投影的SAR成像算法与GPU加速研究222yxumn成像网格中第mn,个网格点，u为雷达天线位于0,u时，天线相mnc位中心到场景x,y的双程延时。mn因此，后向投影算法的过程是一个点对点的图像重建过程，在反投影过程中需要数量巨大的插值操作，导致运算量巨大。以NN×的回波数据为例（其中N和N分别为距离向和方arar位向点数），成像网格方位向和距离向点数分别为N和N，假设N、N和N都为N，那iairaiair3么后向投影中插值的计算量就为约为ON。除了插值运算本身之外，斜距计算和相位补偿等操作在CPU计算中也是非常耗时的，这些因素制约了该算法在CPU平台上的应用和发展。············Grid0SPSP···············································SPSP·············32×SP············…············……····················Block0…SPSP····················Block1SM········Block2SMSM成像网格软件硬件图5.1后向投影算法的GPU并行计算分配图对后向投影算法的分析可知，不同网格点的反向投影操作过程完全相同并相互独立，这使得后向投影算法在GPU上并行实现成为可能。我们可以通过为每个成像网格点分配一个独立工作的线程，每个线程都单独执行斜距计算、插值、相位补偿等操作，最终完成对每个像素点的反投影计算。后向投影算法的GPU并行计算分配示意图如图5.1所示，对于NN的成像网iair格，一共分配NN个并行线程，其相应的并行度为NN。iairiair[59]值得注意的是，在GPU中所谓的“并行”执行并不严谨，实际情况下是“并发”执行的。由于硬件资源的限制，在理想条件下，并行运算所能获得的加速比等同于并行度，但由于并行处理所涉及的GPU资源分配及调度等原因，在实际情况下加速比要低于并行度。尽管如此，并行度为NN的反投影操作经过并行处理后依然能获得相当高的加速比。iair5.2BP成像算法的并行实现在实际处理中，要根据GPU的架构和编程模型，充分利用GPU的硬件资源来对成像算法进行优化。BP成像算法主要包括距离压缩以及反投影两部分，经过上一部分的分析可知，充分利用GPU的高速并行能力可以大大提高该算法的计算效率。实际CUDA程序优化时，我们需要解决很多问题。首先需要对算法进行划分，将串行部分交给CPU来处理，并行部分交给GPU来处理。然后利用kernel函数来处理并行部分，如果在编程过程中需要的线程操作较少，可以采用while(tid
