《基于gpu的机载sar成像算法的设计与实现》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。
巧參觀於处圍硕±学位论文胃１基于ＧＰＵ的机载ＳＡＲ成像算法的设计与实现Ｉ作挑Ｓ麵指导教师姓名、职称邢孟道教授＿申请学位类别工学極圭Ｉ 西安电子科技大学学位论文独创性（或创新性）声明秉承学校严谨的学风和优良的科学道德，本人声明所呈交的论文是我个人巧导师指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别加ｙ？标注和致谢中所罗列的内容Ｗ外，论文中不包含其他人己经发表或撰写过的研巧成果；也不包含为获得西安电子科技大学或其它教育机构的学位或证书而使用过的材料一。与我同Ｘ作的同事对本研巧所做的任何贡献均已在论文中作了明确的说明并表示了谢意。一学位论文若有不实之处，本人承担切法律责任。＇＇‘ｌ＇本人豁名；Ｉ：日期＾皮ｈ蘇＇辟Ｊ西安电子科技大学关于论文便用授权的说明］本人完令了解西安电子科技大学有关保留和使用学位论文的规定，目：Ｉ研巧生在校攻读学位期间论文工作的知识产权属于西安电子科技大学。学校有权保留送巧论文的复印件，允许查阅、借阅论文；学校可Ｗ公布论文的全部或部分内容，允许采用影印、缩印或其它复制手段保存论文。同时本人保证，结合学位论文研究成果完成的论、文发明专利等成果，署名单位为西安电子科技大学。保密的学位论文在年解密后适本书。＿用授权签：：本人名导师签名．心‘如至苗１主兴＞护日期：Ｈ期： 学校代码10701学号1302121129分类号TN95密级公开西安电子科技大学硕士学位论文基于GPU的机载SAR成像算法的设计与实现作者姓名：张锐一级学科：信息与通信工程二级学科：信号与信息处理学位类别：工学硕士指导教师姓名、职称：邢孟道教授学院：电子工程学院提交日期：2015年12月 AirborneSARImagingAlgorithmDesignandImplementationBasedonGPUAthesissubmittedtoXIDIANUNIVERSITYinpartialfulfillmentoftherequirementsforthedegreeofMasterinSignalandInformationProcessingByZhangRuiSupervisor:XingMengdaoProfessorDecember2015 摘要摘要随着合成孔径雷达的不断发展，其在国防工业和民用领域的应用越来越广泛。为了更加充分地发挥SAR全天时、全天候、远距离、高分辨率、广域观测的特点，实际应用中，对SAR的分辨率和测绘带的要求不断提高。机载宽波束高分辨率SAR是实现SAR高分辨率和宽测绘带的一种方式。但是，机载宽波束高分辨率SAR的合成孔径时间长，数据量和计算量较大，在传统的基于CPU的工作站或者服务器上进行数据处理会耗费比较长的时间，无法满足实时性要求。通用计算图形处理器（GPGPU）的出现，GPGPU为宽波束高分辨率SAR数据的实时处理提供了一种可行的途径，通用计算GPU专为加速计算密集型的应用而设计，如今已经发展成具有高度并行、多线程、多核心、超大带宽、具有数百个计算单元的高性能处理平台。本文首先介绍SAR成像的基本原理和通用计算GPU的硬件架构，然后借助NVIDIA的统一计算架构（CUDA）,采用通用计算GPU的对机载宽波束高分辨率SAR成像算法进行实现，并针对GPU的特点对算法的实现进行优化加速，本文主要研究内容总结如下：1、针对传统RMA算法在处理机载宽波束SAR回波数据时频谱利用率不高的问题，对传统的RMA算法进行了改进，并采用一种对数据分段的方案使用GPU对算法进行了实现和优化。该方案解决了在GPU显存不足以容纳一景SAR数据时，数据处理环节与内存/显存间数据传输环节的并行化问题。并针对算法和GPU硬件的特点，使用共享内存优化技术和开启多个CUDA流等方法对算法进行了优化，充分利用了GPU设备的计算资源。使用NVIDIAK40c显卡对算法进行了测试，测试结果表明相对于传统的CPU，GPU能实现14倍左右的加速比。2、由于机载宽波束高分辨率SAR的合成孔径时间长，载机飞行轨迹不稳定，会导致SAR图像散焦和几何失真，需要对回波数据进行运动补偿，传统的基于全孔径的成像处理算法难以满足实时成像的要求，本文采用一种基于子孔径的机载宽波束高分辨率SAR成像处理方法，并用GPU并行处理技术对该成像处理方法的耗时部分进行优化加速，采用并行归约算法和合并内存访问的方法提高了算法执行效率，然后使用MATLAB和CUDAC混合编程技术对算法的各个耗时模块进行实现并封装成MEX文件，使得MATLAB可以直接调用，从而实现快速成像的要求。最后使用NVIDIAK40c显卡对程序的性能进行了测试，实验结果表明该算法的实现方案能对实测数据进行良好的成像并实现9倍左右的加速比，提高了算法的执行效率。关键词：合成孔径雷达（SAR），实时，通用计算图形处理器（GPGPU），统一计算架构（CUDA）I ABSTRACTABSTRACTWiththedevelopmentofsyntheticapertureradar(SAR),ithasbeenwidelyusedinthenationaldefenceandcivilindustry.SinceSARiscapableofobservingearthsurfacedayandnight,all-weather,longrange,highresolution,widearea,higherresolutionandwiderswatharerequiredinupdatingSARsystem.WideazimuthbeamhighresolutionairborneSARisawaytoachievehighresolutionandwideswath.However,inWideazimuthbeamhighresolutionSAR,thesyntheticaperturetimeisverylongandtheamountofdataandcomputationisverylarge,dataprocessingonatraditionalCPU-basedworkstationorserverisrathertime-consuming,makingreal-timeprocessingofSARdataimpossible.WiththeemergenceofGPGPU,GPGPUprovidesafeasibleapproachtothereal-timeprocessingofhighresolutionSARdata.GPGPUisdesignedforacceleratingcomputing-intensiveapplications.Ithasnowdevelopedintoahighlyparallel,multi-threaded,multi-coreandlarge-bandwidthplatformwithhundredsofhigh-performanceprocessingunit.FirstthisthesisintroducesthebasicprinciplesofSARalgorithmsandthehardwarearchitectureofGPGPU,thenbasedonComputeUnifiedDeviceArchitecture(CUDA)technology,itimplementsthewideazimuthbeamhighresolutionairborneSARalgorithmwithGPUandoptimizethealgorithmforcharacteristicsofGPU.Themaincontentofthisthesisissummarizedasfollows:1.AsthelowspectrumutilizationoftraditionalRMAalgorithm,anewplanforanimprovedRMAalgorithmoperatedonaGPUisproposed.ThenewproposalmakesitpossibleforthedataprocessingprocedureandtheCPU/GPUdataexchangetoexecuteconcurrently,especiallywhenthesizeofSARdataexceedsthetotalGPUglobalmemorysize.InordertoexploitthetotalcomputationalresourcesofGPU,theGPUtechnologyofsharedmemoryandCUDAstreamisusedtooptimizetheplan.IthasbeenshownbyanexperimentonanNVIDIAK40cthattheproposedplanacceleratesSARdataprocessingbyabout13times.2.ForawideazimuthbeamhighresolutionairborneSARsystem,thenon-idealmotionofradarplatformvelocitiesresultinseriousblurringandgeometricdistortion.Themotioncompensationneedtobeperformed.Thetraditionalmethodbasedonfullaperturedatacannotmeetreal-timeprocessingrequirement.Inthisthesis,anewmethodbasedonsubapertureoperatedonaGPUisproposed.Parallelreducealgorithmandconsolidationmemoryaccesstechnologyisperformedtoimprovetheefficiencyofthealgorithm.ThenIII 西安电子科技大学硕士学位论文useMATLABandCUDACmixedprogrammingalgorithmtoimplementthevarioustime-consumingmodulesandpackagedintoMEXfile,sothatMATLABcanbecalleddirectly.IthasbeenshownbyanexperimentonanNVIDIAK40cthattheproposedsolutionacceleratesSARdataprocessingbyabout9times.Keywords:SAR,real-time,GPGPU,CUDAIV 插图索引插图索引图1.1HSA异构系统架构..................................................................................................................3图2.1LFM信号的时域波形和幅频特性..........................................................................................8图2.2LFM信号的匹配滤波..............................................................................................................8图2.3对LFM信号匹配滤波后的波形.............................................................................................9图2.4条带SAR的成像几何模型...................................................................................................11图3.1并行处理的关键技术.............................................................................................................13图3.2早期GPU的图像渲染管线...................................................................................................16图3.3G80显卡的标量处理器结构.................................................................................................17图3.4TeslaK40显卡硬件结构........................................................................................................18图3.5KeplerG110架构...................................................................................................................19图3.6SMX的组成结构图...............................................................................................................20图3.7SMX内置的Warp调度器.....................................................................................................21图3.8Kepler架构的存储层次.........................................................................................................22图3.9CUDA程序在GPU中的执行流程.......................................................................................23图3.10CUDA的线程层次图..........................................................................................................24图3.11CUDA的存储器层次图.......................................................................................................25图3.12CUDA异构编程模型..........................................................................................................26图4.1正侧视SAR成像几何模型...................................................................................................28图4.2基频回波的二维波数谱在KK平面的支撑区...............................................................29xr图4.3基频回波的二维波数谱在KK平面的支撑区...............................................................29xy图4.4K取值范围的选取..............................................................................................................30y图4.5改进的RMA算法的实现流程..............................................................................................32图4.6线程块大小对CUDA中warp占用数量的影响..................................................................35图4.7使用3个流时的任务执行时间线.........................................................................................37图4.8分配和释放一次显存的方案.................................................................................................38图5.1子孔径划分示意图................................................................................................................45图5.2算法流程图............................................................................................................................47图5.3MATLAB数据类型...............................................................................................................48图5.4mexFunction()函数的执行流程.............................................................................................52图5.5并行归约算法示意图............................................................................................................53图5.6实测数据处理结果.................................................................................................................54V 表格索引表格索引表4.1运行环境和仿真参数........................................................................................39表4.2两种方案耗时对比............................................................................................39表5.1mexFunction函数的输入参数...........................................................................49表5.2程序运行环境和雷达参数.................................................................................54表5.3两种方案耗时对比............................................................................................55VII 符号对照表符号对照表符号符号名称t距离向快时间tm方位向慢时间c光速（波速）波长Tp脉冲时宽调频率m多普勒调频率fc中心频率fr距离向频率fa方位向频率（多普勒频率）B带宽t0系统的时延fa多普勒带宽TBP时宽带宽积时间分辨率a方位向分辨率R斜距R合成孔径中心到目标的距离0R最短斜距Bv载机速度BW波束方位向宽度Ta合成孔径时间D天线有效尺寸PRF脉冲重复频率波束斜视角rect矩形窗函数冲击响应函数FFT傅里叶变换IFFT逆傅里叶变换K方位波数xIX 西安电子科技大学硕士学位论文K距离波数RDimx_线程网格x轴维度Dim_y线程网格y轴维度()a距离向窗函数r()a方位向窗函数aX 缩略语对照表缩略语对照表缩略语英文全称中文对照SARSyntheticApertureRadar合成孔径雷达R-DRange-Doppler距离-多普勒RCMCRangeCellMigrationCorrection距离单元徙动校正CSChirpScaling线频调变标RMARangeMigrationAlgorithm距离徙动算法PFAPolarFormatAlgorithm极坐标格式算法GPUGraphicsProcessorUnit图形处理器GPGPUGeneral-PurposeGraphicalProcessorUnit通用图形处理单元CUDAComputeUnifiedDeviceArchitecture统一计算设备架构HSAHeterogeneousSystemArchitecture异构系统架构hUMAheterogeneousUniformMemoryAccess异构统一寻址APUAcceleratedProcessingUnit加速处理单元CPUCentralProcessorUnit中央处理器SMStreamingMultiprocessors流处理器簇SISDSingleInstructionSingleData单指令单数据SIMDSingleInstructionMultipleData单指令多数据MISDMultipleInstructionSingleData多指令单数据MIMDMultipleInstructionMultipleData多指令多数据T<ransforming&Lighting顶点转换和光照处理SPStreamingProcessor流处理器GPCGraphicsProcessingCluster图形处理簇SFUSpecialFunctionUnit特殊函数单元PRFPulseRepetitionFrequencies脉冲重复频率LFMLineFrequencyModulation线性调频POSPPrincipleOfStationaryPhase驻定相位定理TBPTimeBandwidthProduct时宽带宽积PGAPhaseGradientAutofocus相位梯度自聚焦FBPFastBackProjection快速后向投影XIII 目录目录摘要........................................................................................................................................IABSTRACT........................................................................................................................III插图索引..............................................................................................................................V表格索引............................................................................................................................VII符号对照表.........................................................................................................................IX缩略语对照表..................................................................................................................XIII第一章绪论......................................................................................................................11.1SAR成像算法研究现状.....................................................................................11.2GPU通用计算的发展及现状.............................................................................21.3课题研究背景......................................................................................................31.4论文主要内容及安排..........................................................................................4第二章SAR成像基础理论.............................................................................................72.1脉冲压缩技术......................................................................................................72.2合成孔径雷达两维高分辨率原理......................................................................92.2.1合成孔径雷达距离向高分辨率原理.....................................................102.2.2合成孔径雷达方位向高分辨率原理.....................................................112.3本章小结............................................................................................................12第三章GPU硬件架构和编程模型...............................................................................133.1并行处理技术....................................................................................................133.2GPU硬件架构...................................................................................................153.2.1GPU硬件架构发展历程........................................................................153.2.2NvidiaTeslaK40显卡硬件架构............................................................183.3CUDA编程模型...............................................................................................223.4本章小结............................................................................................................26第四章一种改进的RMA算法的GPU实现和优化...................................................274.1传统的RMA算法.............................................................................................274.2一种改进的适用于宽波束大场景的RMA算法.............................................304.3算法的GPU实现和优化方案..........................................................................314.3.1总体方案设计.........................................................................................314.3.2数据并行划分.........................................................................................324.3.3GPU存储优化........................................................................................34XV 西安电子科技大学硕士学位论文4.3.4程序实现步骤..........................................................................................384.4实验结果及分析...............................................................................................394.5本章小结...........................................................................................................40第五章基于GPU和MATLAB的SAR成像算法实现及优化.................................415.1一种机载宽波束子孔径SAR成像算法..........................................................415.1.1存在线性误差时时聚焦结果的推导.....................................................415.1.2线性误差项的估计.................................................................................455.2MATLABCUDAC混合编程技术..................................................................475.2.1MATLABCUDAC混合编程技术简介................................................475.2.2MATLAB调用CUDAC程序的实现方法...........................................485.3基于GPU和MATLAB的成像算法实现和优化...........................................505.3.1总体方案设计及数据划分.....................................................................505.3.2GPU的存储优化....................................................................................525.3.3实验结果及分析-...................................................................................535.4本章总结...........................................................................................................55第六章总结与展望...........................................................................................................576.1论文总结...........................................................................................................576.2工作展望...........................................................................................................57参考文献.............................................................................................................................59致谢.....................................................................................................................................61作者简介.............................................................................................................................63XVI 第一章绪论第一章绪论1.1SAR成像算法研究现状SAR是一种通过重建目标的后向散射特性对目标进行成像的雷达。它可以在低可见度的天气下获取高分辨率的图像，同时对植被和墙体具有高可穿透性。目前，SAR广泛应用于民用和军用领域，在民用领域，它广泛地应用于资源勘测、自然灾害的监测、农作物的评估等领域；在军用领域，也广泛地运用于军事目标的侦查、敌情的收集等方面[1]。SAR成像就是从回波信号中准确地重构目标的后向散射特性，SAR的距离向高分辨率是通过发射宽带的线性调频信号，对回波信号进行脉冲压缩处理得到的。方位向的高分辨率则是通过雷达与目标的相对运动来构成合成阵列，得到一系列顺序阵元信号，再进行相干叠加得到的。SAR成像的难点在于，由于距离徙动的存在，点目标的系统响应在距离向和方位向的二维平面呈现为曲线，使基于匹配滤波的成像的实际计算复杂化[2]。因此，成像算法的主要工作是对目标的回波响应进行距离徙动校正，将响应的二维耦合进行解耦，从而使目标精确聚焦成像。目前，比较成熟的成像算法主要有距离多普勒（RD）算法、CS（ChirpScaling）算法和距离徙动算法（RMA）等。距离多普勒算法是在1978年为高效处理SEASATSAR数据而提出的，至今仍广泛使用[3]。RD算法被称为距离多普勒的原因是它在距离时域方位频域中对回波的距离徙动进行校正，它的基本思想是利用距离和方位在一定条件下的可去耦性，在两个一维操作之间使用距离徙动校正，将二维成像处理转换为两个一维脉冲压缩处理过程。其中，距离压缩只需在距离频域对SAR回波信号做匹配滤波，而方位压缩则利用了同一距离单元上不同目标的徙动曲线在距离时间-方位频率域上近似一致的特点，在距离多普勒域中先进行徙动校正，然后方位匹配滤波，从而实现目标的方位分辨率[4]。CS算法的基本原理在距离多普勒域对信号乘以一个CS因子，将原始信号距离向的中心频率做小的频移，使在不同距离单元上的点目标具有相同的徙动曲线，然后将信号变换到二维频域，在二维频域完成距离徙动校正，最终通过距离压缩和方位压缩来实现图像的聚焦[5]。距离徙动算法源于地震信号处理，它基本处理步骤是先将回波信号变换到二维波数域，然后乘以一个参考函数，使参考距离处的目标完全聚焦，非参考距离处的目标部分聚焦，再通过Stolt插值完成非参考距离出目标的聚焦。RMA算法在算法推导过程中没有使用任何数学近似，它是原理上最优的成像算法。但是，RMA算法中存在1 西安电子科技大学硕士学位论文插值操作，会使算法实现时计算量比较大。另外，如果插值不准确的话，对成像质量也会产生影响。1.2GPU通用计算的发展及现状真正意义上能对图像数据进行处理的GPU是NVIDIA在1999年首先研制成功的，当时的GPU主要是为CPU分担一些固定的图形计算工作，功能比较单一，而且不能进行编程，是一种功能固定的专用处理器。2002年，NVIDIA在GPU中加入了着色器，使得GPU第一次具有了可编程性，此时的GPU可以图形处理接口来进行一些图形处理之外的一些其它计算工作，但是实现过程非常复杂，而且编程限制比较多。2007年，NVIDIA公司发布了G80系列显卡，该款显卡支持微软提出的统一计算架构，这是GPU发展史上的一款革命性显卡，它标识这GPU通用计算时代的到来。同时，NVIDIA为GPU增加了一个易用的编程接口，也就是所谓的统一计算架构（ComputeUnifiedDeviceArchitecture，CUDA）。这使得开发人员无须学习复杂的着色语言或者图形处理原语，就能使用GPU进行通用计算。2009年，NVIDIA推出了基于“费米”架构的GPU,该系列的计算能力大大提升，特别是对双精度浮点数的支持方面大大加强。2012年，NVIDIA又推出了的“开普勒”架构的GPU,该系列在软硬件的结合方面得到了加强，性能也相对“费米”架构提升了3~4倍。目前，在通用计算方面最强劲的GPUteslaK80就是基于“开普勒”架构的。如今，GPU在通用计算中的诸多工程领域都得到了应用，如云计算、人脸识别技术，搜索引擎技术等等，解决了大量实际的工程问题，也推动了神经网络、深度学习等一些学科的发展。另外，GPU也广泛应用于一些基础学科的研究领域，如物理学，气象学等等。GPU在超级计算机领域也越来越受到重视，大量的超级计算机系统开始采用GPU作为集群节点中的计算单元，连续三次摘得全球超级计算机TOP500桂冠的“天河二号”计算机就大量采用了GPU[6]。AMD公司在收购了ATI公司后，也推出了很多高性能的用于通用计算的显卡，不同与NVIDIA，AMD公司推出的GPU的编程模型是基于开源的OpenCL的。由于AMD公司在CPU领域也拥有强大的研发实力，因此，近年来，AMD公司提出了一种新的体系结构HSA（HeterogeneousSystemArchitecture，异构系统架构），并于2015年推出了一种新的产品APU（AcceleratedProcessingUnit）。APU是基于AMD的新技术—异构统一寻址技术（heterogeneousUniformMemoryAccess，hUMA），将CPU和独立显卡中的GPU集成到一个芯片上，通过hUMA技术，CPU和GPU可以使用2 第一章绪论相同的方式对主机内存进行寻址，即CPU和GPU可以共享主机内存，而不需要额外在显卡中配置大量的显存，这样在使用GPU进行计算时就不必再将数据在CPU端和GPU端来回复制，节省了大量的存储带宽，可以大大提高GPU计算的性能。hUMA(MEMORY)GPUCPUPARALLELWORKLOADSSERIALWORKLOADS图1.1HSA异构系统架构目前，在使用GPU进行通用计算时，由于CPU和GPU各自拥有自己的存储空间，使得数据需要来回在CPU端和GPU端复制多次，耗费大量时间，而且由于GPU显存的大小无法做的很大，使得编程的复杂性提高，程序的性能下降。一旦HAS技术趋于成熟，上述这些问题都将不复存在，必将把GPU通用计算的性能提升到一个新的高度，甚至引发一场GPU体系架构的革命。1.3课题研究背景随着SAR成像处理技术的发展，如今，SAR成像的分辨率越来越高，成像的测绘带大小也越来越大，使得SAR回波的数据量也越来越大，另一方面，人们对SAR成像处理的实时性的要求却越来越高。要处理这些数据量大、分辨率要求高的SAR回波数据，需要有好的能够满足实时性要求的成像算法，同时，也需要好的硬件实现方式。对于机载宽波束高分辨率SAR，合成孔径时间长，方位距离耦合严重，数据量大，传统的全孔径处理方法难以满足实时性要求，需要对算法进行改进。因此，本文采用一种新的成像方案，先将全孔径数据划分为多个子孔径，然后对子孔径数据进行运动补偿并使用一种改进的RMA算法对其成像，再将其子孔径成像结果投影到全孔径的成像网格中，得到全孔径图像，从而能够适应实时性要求。在机载宽波束高分辨率SAR成像处理过程中，大的数据量和计算量为SAR信号处理带来了新的挑战。针对这些挑战，必须要有高效率的信号处理器。而用于通用计算的GPU中有大量的运算核，专为处理数据量大、计算量密集的应用而设计，非常适合对SAR回波数据进行处理。相对与传统的运算平台，充分利用GPU中的运算资源可以获得十几倍甚至几十倍的性能提升。它为解决高分辨率SAR信号处理中面临的计算量大的难题提供了一条可行的道路。机载宽波束高分辨率SAR的数据处理过程中，有大量的向量相乘、向量转置、3 西安电子科技大学硕士学位论文插值和FFT等操作，这些操作如果用传统的方法在基于CPU的工作站上进行处理，由于CPU串行处理的特性，必然会耗费大量的时间[7]。同时，随着近年来CPU的制程逐渐达到极限，摩尔定律开始失效，CPU计算性能的提高越来越困难。而在实时、高清3D图形的巨大市场需求的驱动下，GPU却在通用计算领域取得了飞速的发展，在计算能力和内存带宽方面都远远超过了CPU。CPU和GPU之间在浮点计算能力上的差异的原因是GPU中有大量数据处理的单元，能同时开启大量的线程，线程在GPU内部以线程块为单位进行任务调度，由于GPU中有大量的寄存器使得GPU中的上下文切换是轻量级的。同时，区别与CPU，GPU的处理效率是随着线程数的增加而增加的，所以，GPU非常适合处理计算量大的问题。而在高分辨率宽幅SAR成像中，大量耗时且基本的操作都可以用并行处理技术以较高的效率将其并行化。因此，可以采用GPU利用NVIDIA的CUDA技术对机载SAR进行实时成像处理。目前，GPU的编程模型还是基于CPU-GPU的异构模式。因此，在实际的处理中，CPU主要负责读写数据和处理算法的整体框架，而对于大量计算密集型的计算模块则交由GPU处理，通过相关的并行处理技术和编程技巧，对程序进行优化，充分利用GPU的计算资源和带宽，来满足实时成像的要求。1.4论文主要内容及安排机载宽波束高分辨率SAR能获得宽测绘带高分辨率的SAR图像，在军事和民用方面均有广泛的应用前景。本文主要研究了GPU在机载宽波束高分辨率SAR数据处理过程中的应用。论文的主要内容及安排如下：第一章是绪论部分，主要介绍了SAR成像算法和GPU通用计算的研究现状以及课题的研究背景。第二章介绍了合成孔径雷达相关的基础理论。第三章主要介绍了并行计算基础理论和GPU的硬件架构，并重点对NVIDIATeslaK40显卡的硬件架构和CUDA编程模型进行了讨论。在第四章中，针对机载宽波束高分辨率SAR成像过程中，由于方位波束较宽，传统RMA算法二维频谱利用率不高的问题，对传统RMA算法进行了改进，并针对实测数据处理时数据量和计算量大的问题，采用GPU对算法进行优化加速。在第五章中，针对在算法开发过程中，算法设计方案需要经常修改的问题。为了能在算法开发过程中使用GPU对成像算法进行加速，本文采用MATLAB和CUDAC混合编程技术，对机载宽波束高分辨率SAR中涉及的比较耗时的模块使用GPU进行加速，并编写成MEX文件的形式供MATLAB调用，提高了算法的执行效率，也方便了在算法过程中对这些模块进行灵活调用。4 第一章绪论第六章是总结与展望，主要是对本文所做工作的总结，并对后续需要做的工作进行展望。5 西安电子科技大学硕士学位论文6 第二章SAR成像基础理论第二章SAR成像基础理论2.1脉冲压缩技术雷达的发射信号的脉冲宽度决定雷达的作用距离，信号的带宽决定雷达的距离分辨率，二者都是雷达的重要性能，但是，对于常规的脉冲雷达，这两者却不能兼顾[8]。因为常规的雷达发射的脉冲信号的载频是单一的，这种信号的时宽和带宽的乘积为一个定值（约为1），不能同时增大时宽和带宽，即不能同时改善雷达的作用距离和距离分辨率，所以需要考虑其他的具有大时宽带宽积的更复杂的信号形式。针对上述问题，可以采用脉冲压缩技术来解决这两者之间的矛盾。使用脉冲压缩技术的雷达发射的信号，必须要有大的时宽带宽积，这种信号不再是简单的单一载频脉冲，而是经过了某种调制[9]。线性调频信号就是一种符合上述条件的信号，目前，大多数脉冲雷达都采用了这种信号。LFM信号的数学表达式为：2tjft2(c)st()rect()e2(2-1)Tpt这里,f表示雷达载波的频率，rect()是矩形信号，其表达式如下：cTptt1,1rect()T(2-2)pTp0，其它B表示调频斜率，根据信号的调频斜率，可以得到瞬时频率的表达式为：TpfKtT(/2tT/2)。cpp式(2-1)可以还写成另一种形式：st()Ste()jf2ct(2-3)式中，tjt2St()rect()eTp(2-4)St()是信号st()的复包络。由式(2-3)可以知道，st()是St()频移后得到的，两者的幅度谱相同。对式(2-4)用MATLAB进行仿真，并画出信号在时域的波形和在频域的幅度谱，如图2.1所示。7 西安电子科技大学硕士学位论文线性调频信号的实部0.50-0.5-1-5-4-3-2-1012345Timeinus线性调频信号的幅频特性40302010-30-20-100102030FrequencyinMHz图2.1LFM信号的时域波形和幅频特性压缩网络是脉冲压缩体制的另一个重要的组成部分，本质上，一个理想的压缩网络是一个匹配滤波器，信号st()的匹配滤波器的时域脉冲响应为：*ht()stt()(2-5)0t是系统的时延，用于保证系统物理可实现[10]。在进行理论推导时，可以取t＝000，此时，式(2-5)可以写为：*ht()s()t(2-6)匹配滤波图2.2LFM信号的匹配滤波如图2.2,st()经过系统ht()得输出信号st():otsinT(1)pTtstT()prect()ejf2ct(2-7)0pTt2Tpp式(2-7)是一个单频信号，由此可知，雷达发射的线性调频信号经过匹配滤波处理8 第二章SAR成像基础理论后，得到的是一个单频信号。当tT时，包络近似为辛克（sinc）函数。pttStTSaTtrect()()()TSaBtrect()()(2-8)0pppTT2ppsinc函数的脉冲宽度定义为sinc函数3db点间的宽度，如图2.4所示，当Bt21时，t，此时，脉冲宽度的大小为：2B112(2-9)2BB脉冲压缩比为：TpDBT(2-10)p1因此，使用匹配滤波器对信号进行滤波后，信号的脉冲宽度变为。以下是用Bmatlab仿真的图2.2的过程：线性调频信号经过匹配滤波器后的波形0仿真结果sinc函数-20,dB幅度-40-15-10-5051015TimeinsecB线性调频信号经过匹配滤波器后的波形（放大）0-4-13.4,dB幅度-3-2-1-0.500.5123TimeinsecB图2.3对LFM信号匹配滤波后的波形2.2合成孔径雷达两维高分辨率原理由于合成孔径雷达的分辨单元很小，所以它可以对雷达波束照射范围内的目标进行成像。在雷达成像的过程中，将场景中雷达载机飞行的方向成为方位向，与雷达载9 西安电子科技大学硕士学位论文机飞行方向垂直的方向则称为距离向[11]。SAR成像的关键是提高雷达的距离向和方位向两个维度的分辨率，当分辨率足够高时，雷达的分辨单元就会远小于目标的尺寸，从而就能够对目标成像。本章分两节，分别介绍SAR在距离维获得高分辨率的原理和在方位维获得高分辨率的原理。2.2.1合成孔径雷达距离向高分辨率原理SAR系统发射的雷达波形通常是LFM信号，当LFM信号的时间带宽积远大于1时，对雷达回波信号进行脉压处理后，雷达的距离分辨单元就反比与信号的带宽，这样，通过增大雷达发射信号的带宽，就能提高雷达的距离向的分辨率。本节主要根据上节阐述的脉冲压缩理论，分析SAR是如何提高距离向分辨率的。假设雷达发射的LFM信号表示为：2tjf2(ctt)st()rect()e2Tp(2-11)这里，T是发射的LFM信号的脉冲宽度，f是雷达载波的中心频率，是调频pc斜率。则信号的带宽为BT，时间带宽积为DBT。利用傅立叶变换，将式(2-11)pp变换到频域，就可以得到信号的幅度谱和相位谱。由2.1节中的分析可以知道，若满足BT1，则信号经过匹配滤波器完成匹配滤波处理后，信号的包络变为如下p的形式：sin[B(tt)]0T(2-12)pBtt()0由式(2-12)可以看出，信号的包络的表达式是sinc函数的形式。sinc函数的主瓣是两个3dB点之间的部分。通过上节的分析可以得到，脉压后信号的脉冲宽度为112，由此得到，式(2-11)中的LFM信号的脉冲压缩比的表达式为2BBTpDBT。可以看出，LFM信号经过脉压后的脉冲压缩比等于LFM信号的时间p带宽积。邻近目标的名义分辨率为：1cc(2-13)r22B由上式可以看出，SAR系统在距离向的分辨率和雷达发射信号的带宽成反比，雷达发射信号的带宽越宽，距离向的分辨率越高[12]。因此，一般的SAR系统发射的LFM信号的带宽都很大，从而能够得到高的距离分辨率。10 第二章SAR成像基础理论2.2.2合成孔径雷达方位向高分辨率原理tm0vtmR(,)tRRm00图2.4条带SAR的成像几何模型由上节可知，合成孔径雷达距离向的高分辨率是通过发射大时间带宽积的线性调频信号，在接收时采用脉冲压缩技术获得的；方位向的高分辨率是通过雷达与目标之间的相对运动获得的[13]。可以证明，雷达的方位向回波信号近似为一个线性调频信号，这样，我们也可以通过脉冲压缩的方法实现方位向的高分辨率。对于条带SAR的情况，它的成像几何模型如图2.4所示。AB为一个合成孔径长度，t0时刻设在合成孔径中心，载机的速度为v，合成孔径中心到目标的距离为mR，在t时刻，载机与目标的距离为R(,)tR。雷达的斜视角为。由图2.4中成像0mm0几何模型的几何关系，使用余弦定理可以计算得到：22RtR(,)R()2()sinvtRvt(2-14)mm000m对于一般的SAR有Rvt，所以可以将R(,)tR在t0处展开，并忽略3次0mm0m及以上项，得(2-14)式的近似式为：2cos22R(,)tRRvtsin()vtRftt(2-15)mm00md0cmmm22R20222v2cvos其中，fsin为多普勒中心频率，为多普勒调频率。dcmR0SAR回波信号的多普勒项为4jR(,)tRm0stRe(,)(2-16)am0将(2-15)式代入(2-16)式得：11 西安电子科技大学硕士学位论文24stR(,)exp(2jftjtjR)(2-17)am00dcmmm由式(2-17)可以看出，方位向信号具有LFM信号的形式，根据上节的分析，只需对上式进行脉压处理就可以得到方位向高分辨率。根据图2.4中的成像几何模型，可以计算得到雷达的合成孔径时间的近似表达式为RR00T(2-18)scosDcos由式(2-17)可以计算得到方位向信号的调频率为：222cvos(2-19)mR0将式(2-18)和式(2-19)相乘，就可以得到方位向信号的带宽：2vfTcos(2-20)dmsD根据上节的分析，由(2-20)式可以计算得到方位信号进行脉压处理后的时宽：1D(2-21)0fv2cosd则点目标的横向分辨率为：Dv(2-22)a02cos由式(2-22)可以看出，在条带SAR系统中，方位向分辨率与雷达天线孔径成正比，而与天线斜视角的余弦值成反比。但是，实际中，当斜视角过大时式(2-22)的关系将不能得到满足，而天线孔径的大小也不可能做的很大，所以方位分辨率的大小受到很多因素的限制[14]。2.3本章小结SAR系统发射的波形是LFM信号，本章基于LFM信号对SAR的二维分辨率进行了分析。首先，通过对匹配滤波过程的推导，介绍了脉冲压缩技术，并对脉压的过程进行了仿真，之后，分两节分别推导了SAR距离维和方位维分辨率的表达式，阐述了SAR二维高分辨率的原理。12 第三章GPU硬件架构和编程模型第三章GPU硬件架构和编程模型3.1并行处理技术并行处理技术的提出最早可以追溯到上世纪50年代。并行处理技术从形式上划分可以分为时间并行和空间并行，时间并行主要是采用流水线技术，通过让多个不同的处理过程轮流重叠使用同一部件的不同部分，使得不同处理过程在时间上错开，从而使硬件的处理速度加快。空间并行主要是增加运算部件的数量，通过多个运算部件的协同工作来共同解决问题。并行处理技术的实现需要三个方面的技术支撑：并行体系结构、并行编程模型以及并行优化算法[15]。其中，以并行优化算法为核心，以并行体系结构为载体，以并行编程模型为描述方式，三种技术相互联系又相互制约组成一个并行处理系统。在一个并行处理系统中，并行体系结构是并行处理的运行平台，并行编程模型为并行优化算法在并行体系结构上的实现提供接口，并行优化算法是问题解法在并行编程模型基础上的具体算法实现。图3.1展示了并行处理系统中的三种关键技术的层次关系图。应用并行优化算法并行编程模型平台（系统软件和硬件）并行体系结构图3.1并行处理的关键技术并行体系结构的种类多达几十种，分类方式也多种多样。其中，最通用的分类方式是弗林分类法。弗林分类法的分类标准是计算机的运行机制，根据不同计算机的运行机制，它把计算机系统的体系结构分为以下四种类型：SISD(SingleInstructionStreamSingleDataStream，单指令单据)、SIMD(SingleInstructionStreamMultipleDataStream，单指令多数据)、MISD(MultipleInstructionstreamSingleDatastream，多指令单数据)和MIMD(MultipleInstructionstreamMultipleDatastream，多指令多数据)[16]。在单指令单数据（SISD）系统中，一个指令流处理只处理一个数据流，它是最普通的顺序处理串行机。在单指令多数据（SIMD）系统中，有一个控制单元（controlunit）和多个处理单元（processingunit），每个处理单元都接收控制单元发出的相同指令，并针对不同的数据集进行处理。多指令多数据（MISD）系统的提出主要是为了分类的完整性。目前，还没有属于13 西安电子科技大学硕士学位论文[17]此类型的系统。多指令多数据（MIMD）系统拥有多个处理单元，每个处理单元都有自己的控制单元和操作的数据，各个处理单元之间通过共享内存或者互联网络的方式进行通信。因[18]此，MIMD系统中的各个处理单元可以执行不同的指令。广义上来说，前面的三种体系架构都可以对应到MIMD架构上来。目前，绝大多数的并行处理系统都属于MIMD类型的系统。并行编程模型是对并行体系结构的抽象，主要是为并行程序的开发提供并行编程框架和应用程序接口等并行编程环境。并行编程模型依赖于并行体系结构，而好的并行编程模型能够充分地挖掘出硬件的性能，它是硬件和软件之间沟通的桥梁。由于并行体系结构有很多种，所以并行编程模型的种类也有多种。根据这些编程模型的特征，可以将并行编程模型分为以下几类：数据并行编程模型、共享存储编程模型、消息传递编程模型和混合编程模型。数据并行编程模型是通过将数据划分为不同的部分显式地分配到各个进程中，各个进程在已分配的数据集上执行操作，进程之间无须显式地通信。数据并行编程模型的抽象级别比较高，编程比较简单，而且一大类科学计算问题也都可以用数据并行编程模型很好的解决，但是数据并行编程模型在编程的灵活性方面欠佳，仅适用于数据并行的情况[19]。共享存储编程模型的所有线程共享全局地址空间，多个线程之间通过全局地址空间进行通信，所有的数据也都保存在全局地址空间上，所以不需要像数据并行编程模型一样进行数据的显式分配。共享存储编程模型的多个线程之间采用异步并行，需要进行同步操作来协同各个线程之间的工作。共享存储编程模型的实例有很多，像常用的JavaThread、OpenMp都属于这一类编程模型。消息传递编程模型是通过消息传递来实现线程同步、数据共享等一系列的操作，一般在实现时都会提高一套运行时库来供程序员调用，通过调用运行时库来对数据划分和任务分配等各个方面进行控制。消息传递编程模型的数据划分和任务分配都需要程序员来完成，使得消息传递编程模型的灵活性更高，编程的复杂性也增加了，像著名的MPI就属于消息传递性编程模型。高性能计算正在向着异构计算的方向发展，而不同的体系结构也必然需要不同的编程模型来提供支持，在这种情况下，融合多种并行编程模型来对异构计算提高支持就成为一种好的选择，因此混合编程模型就呼之欲出了。混合编程模型可是是以上几种编程模型的组合，它能充分地发挥各个编程模型的优势，来提高并行计算的性能。并行算法的优化可以从任务并行、数据并行和循环并行三个方面来实现。任务并行就是将一个任务划分成多个子任务，且子任务能够并行执行，通过分配各个子任务到各个子线程来实现并行优化。数据并行就是将对数据进行划分，通过分配数据到各个子线程，使得各个子线程在各自的数据集上进行操作来实现并行优化。循环并行就是将一个主循环分成多个可以并行执行的子循环，再将各个子循环独立执行，从而对14 第三章GPU硬件架构和编程模型程序进行优化。在并行计算中，性能的瓶颈往往不是计算资源的稀缺，而是访问内存的延时和通信的开销。因此，在具体的算法优化中，往往集中于减小或隐藏访问内存的延时和通信延时。可以通过将具体计算和访问内存或通信分开进行，从而减小部分访问内存的延时及通信的开销。3.2GPU硬件架构3.2.1GPU硬件架构发展历程1999年，英伟达公司发布了型号为Geforce256的显卡，正式提出了GPU（GraphicsProcessingUnit）的概念。GPU在诞生之初主要是作为协处理器为图形显示进行加速。GPU的硬件架构在最初借鉴了CPU设计中的流水线技术，CPU的流水线技术是将CPU的指令执行分为若干步，每一步使用专门的硬件电路进行处理，这样就能像工厂的流水线一样，流水线中的每一级都能连续不断地执行，从而使不同指令的不同部分可以并行执行，提高指令的执行效率。图像的渲染就是给定三维场景的描述以及观测点的位置和方向，通过计算得到显示器中显示的二维光栅化图像的过程。图像渲染一般也分为顶点生成、顶点转换、光照、裁剪等几个不同的阶段。因此早期的GPU硬件架构主要由两部分组成：负责顶点转换和光照处理（Transforming&Lighting，T&L）的硬件处理单元和像素渲染流水线。T&L硬件处理单元是一个可以进行浮点型运算的功能固定的硬件处理单元。像素渲染流水线是根据像素渲染的几个不同的阶段，分别用不同的硬件电路进行专门处理，形成了一条图像渲染的流水线，这就是传统的固定渲染管线。传统的固定渲染管线一般分为三个部分，即像素着色单元（PixelShaderUnit），主要完成顶点的像素处理；纹理贴图单元（TextureMapUnit），主要完成纹理渲染；光栅化处理单元（RasterOperationsUnits），将三维场景光栅化为可以在显示器上显示的二维图像数据。描述三维场景的矢量数据输入GPU后，经过GPU的T&L处理和图像渲染管线处理后就输出了场景的二维观测图像。传统的图像渲染管线各个部分的功能是固定的，不具备可编程性，虽然对图像渲染进行了加速，但灵活性欠佳，而且图像渲染的效果仍不理想。15 西安电子科技大学硕士学位论文像素着色单元（PixelShaderUnit）纹理贴图单元（TextureMapUnit）光栅化处理单元（RasterOperationsUnits）图3.2早期GPU的图像渲染管线2001年，NVIDIA在它的新一代GPU中引入了一项新的重大革新—Shader（着色器），用来替代传统的功能固定的图像渲染管线。着色器（Shader）分为顶点着色器（VertexShader）和像素着色器（PixelShader）两种。顶点着色器主要负责顶点矩阵转换和光照处理等一系列操作，像素着色器主要负责图像的纹理渲染和三维场景的光栅化。它们都是采用了流水线技术的可编程的硬件逻辑单元。这样GPU就从固定渲染管线架构迁移到了可编程的流水线架构，GPU第一次具有了可编程性。由于顶点数据和像素数据都是四维矢量，所以着色器都被设计成四个通道的SIMD型的处理单元，即四元组结构。这样顶点着色器和像素着色器就能够在一个时钟周期内完成四次算术逻辑运算，处理一个完整的顶点几何转换和像素渲染。着色器的引入使得图像渲染的灵活性和效率大大提高，但是此时的着色器的编程性仍然有限，受到硬件的限制很大。第二代GPU一直采用这种顶点着色器和像素着色器分离的架构，并在此基础上不断改进，着色器的性能和可编程性不断提高。在ShaderModel2.0中，顶点着色器和像素着色器都开始支持双精度浮点运算，这为GPU在通用计算方面的应用奠定了基础。同时，顶点着色器中引入了流程控制，支持的指令数也大大增加，使得编程的灵活性得到加强。最终GPU的硬件架构中取消了早期GPU中固定的T&L处理单元，使得图像渲染完全由着色器完成。在一款GPU的架构设计完成时，它的顶点着色器和像素着色器的比例就会根据实际的计算需求而确定下来。早期的着色器尽管在处理矢量图像数据方面的效率很高，但是，在着色器在处理标量数据或者2D、3D数据时，会有大量计算资源闲置，使得计算效率降低。为了解16 第三章GPU硬件架构和编程模型决这一问题，混合型架构被采用到GPU的设计中，混合型架构中不仅有4D指令执行单元，还会搭配2D、3D、甚至标量指令执行单元，这样使得GPU计算资源闲置的状况得到缓解。尽管如此，计算资源闲置的问题仍然没有彻底解决，着色器的流水线中没有分支预测机制，因此，在遇到分支时，着色器只是简单地让某些处理模块闲置，使得硬件的计算资源没有充分利用，降低了计算效率，这样就使得GPU很难充分地发挥它的计算性能。由于不同图像的渲染对顶点渲染操作和像素渲染操作的需求是不一样的。但是在GPU的顶点着色器和像素着色器分离的架构中，这两种着色器的比例是固定的，这样就造成了图像渲染的过程中计算资源配置的不合理，大量的计算资源被闲置。为了解决这个问题，微软在Directx10中提出了统一渲染架构，将顶点着色器和像素着色器统一成一种相同的处理单元，这种处理单元既能进行顶点渲染，还能进行像素渲染，甚至还能进行几何渲染。这样就最大限度的减少了计算资源的闲置，还大大增加了编程的灵活性，也为GPU进行通用计算铺平了道路。2006年，NVIDIA推出了第一款基于统一渲染架构的GPUG80，并于2007年推出了CUDA（ComputeUnifiedDeviceArchitecture，统一设备计算架构），引发了GPU通用计算领域的一次革命[20]。HostSetup/Rstr/ZCullInputAssemblerVtxThreadIssueGeomThreadIssuePixelThreadIssueProcessorThreadSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPTTTTTTTTFFFFFFFFL1L1L1L1L1L1L1L1L2L2L2L2L2L2FBFBFBFBFBFB图3.3G80显卡的标量处理器结构G80显卡引入了很多革命性的技术，它也是世界上第一款支持微软Directx10的显卡。图3.3是G80的硬件架构，可以看出G80GPU主要由四部分组成，即几何处理单元，8组流处理器簇（SM，StreamMultiprocessers）,6组纹理寻址单元和6组光栅化处理单元。这四部分通过高速交叉总线相连。G80显卡首次将执行4D指令的矢量处理单元都改为了执行1D指令的标量处理单元，每个这样的标量处理单元称为SP17 西安电子科技大学硕士学位论文（StreamProcesser）。每个SP有自己的指令发射单元，每16个SP又组成一个SM，同一个SM中的SP可以通过寄存器文件或者共享存储空间进行通信。同时SM中还集成了一些专用的处理模块（SPU），执行一些特殊操作。这样无论是1D、2D、3D还是4D指令都能够通过SP高效率的执行，而不造成部分资源的闲置。NVIDIA在G80之后的GPU中基本都遵循了这种基本架构，但在SM的处理单元数量和性能方面都有很大的提高，在G100之后还引入了GPC（GraphicsProcessingCluster）的概念，GPC包含有几个SM，同时还有自己的几何处理模块和光栅化处理模块，这样就使得GPU的渲染速度进一步提高。3.2.2NvidiaTeslaK40显卡硬件架构NvidiaTeslaK40显卡是一款基于GK110BGPU的，主要用于工作站和服务器的高性能计算加速卡。如图3.4所示，是TeslaK40显卡的硬件结构。TeslaK40显卡拥有12GB的GDDR5显存，显存的主频达到3GHz。BIOS中包含大小为2M的ROM。K40所使用的GPU的型号是GK110B，这款GPU是基于NVIDIA的Kepler架构的。12GBGDDR524Pieces256M×16384bGPUFBVBIOSPower2MbitGK110BPEXSupplyROMGPIOGen3×1612V/3V312V12VAuxPCIExpressEdgeConnectorPower图3.4TeslaK40显卡硬件结构Kepler架构主要是针对高性能计算（HPC）领域的应用进行加速，支持双精度浮点数。一个完整的KeplerG110GPU包括15个SMX单元和6个64位的存储器管理单元。如图3.5所示，是KeplerG110的硬件架构。18 第三章GPU硬件架构和编程模型PCIExpress3.0HostInterfaceGigaThreadEngineMemoryControllerMemoryControllerSMXSMXSMXSMXSMXSMXSMXSMXL2CacheMemoryControllerMemoryControllerSMXSMXSMXSMXSMXSMXSMX图3.5KeplerG110架构KeplerG110采用了一种新的SM处理器架构SMX。如图3.6所示，是SMX的组成架构图。每一个KeplerG110SMX单元拥有192个单精度CUDA核，每一个CUDA核都完整地拥有一个浮点运算流水线和一个整数和逻辑运算单元。Kepler架构继承了Fermi结构对单精度和双精度浮点运算的良好支持。SMX中还包含有32个SFU（SpecialFunctionUnit），用于一些特殊函数的快速近似计算。19 西安电子科技大学硕士学位论文InstructionCacheWarpSchedulerWarpSchedulerWarpSchedulerWarpSchedulerdispatchdispatchdispatchdispatchdispatchdispatchdispatchdispatchRegisterFilecorecorecoreDPUnitcorecorecoreDPUnitLD/STSFUcorecorecoreDPUnitcorecorecoreDPUnitLD/STSFUcorecorecoreDPUnitcorecorecoreDPUnitLD/STSFUcorecorecoreDPUnitcorecorecoreDPUnitLD/STSFUcorecorecoreDPUnitcorecorecoreDPUnitLD/STSFUInterconnectNetwork64KBSharedMemory/L1Cache48KBRead-OnlyDataCacheTexTexTexTexTexTexTexTexTexTexTexTexTexTex图3.6SMX的组成结构图SMX调度线程时以32个线程为一组进行调度，这32个线程就称为一个warp。在KeplerGK110中，每一个SMX拥有4个warp调度器和8个指令发射器，支持4个warp同时并行执行。SMX内置了warp调度器，内置的warp调度器每次选择4个warp和2个独立的可以分发的指令进行调度。如图3.7所示，是SMX内置的warp调度器示意图。20 第三章GPU硬件架构和编程模型WarpSchedulerInstructionDispatchUnitInstructionDispatchUnitWarp8instruction11Warp8instruction12Warp2instruction42Warp2instruction43Warp4instruction95Warp4instruction96......Warp8instruction13Warp8instruction14timeWarp2instruction43Warp2instruction44Warp4instruction97Warp4instruction98图3.7SMX内置的Warp调度器Kepler架构的存储层次组织如图3.8所示。Kepler架构的存储管理单元支持内存存取的虚拟统一寻址，即无论是CPU的内存还是显存，都可以直接访问，这样就无需程序员显式地在主机端和设备端进行内存的拷贝。在KeplerG110中，每个SMX中拥有64KB的L1缓存，这64KB的缓存除了可以配置为16KB的共享内存和48KB的缓存或者配置为48KB的共享内存和16KB的缓存外，还可以配置为32KB的共享内存和32KB的缓存。另外，SMX中还拥有48KB的缓存由于存放函数计算时的只读数据。除了SMX中的L1缓存外，KeplerGK110中还拥有容量为1536KB专用的L2缓存，比Fermi架构的L2缓存容量提高了一倍。21 西安电子科技大学硕士学位论文ThreadSharedL1Read-OnlyMemoryCacheDataCacheL2CacheDRAM图3.8Kepler架构的存储层次Kepler架构最大的改进之一是增加了动态并行的机制（DynamicParallelism），即允许在GPU线程的内部再开启新的线程，而无需再转到CPU端，这样就提高了程序执行的效率，也大大增加了编程的灵活性。Kepler架构的另一个新特性是支持Hyper-Q技术，即允许多个CPU同时使用同一个GPU,这样就减少了CPU的等待时间，提高了执行效率和可编程性。特别是对在集群中使用的GPU的执行效率和可编程性的提高更加明显。另外，Tesla架构还支持双复制引擎，使得在CPU端内存向GPU端内存复制数据时，GPU端内存可以同时向CPU端内存复制数据，这样使得主机端和设备端的数据交换速度得到提高。3.3CUDA编程模型CUDA（ComputeUnifiedDeviceArchitecture,统一计算架构）是2007年NVIDIA发布的基于NVIDIA的GPU的一个通用并行计算平台和编程模型。主要用于大量复杂计算的GPU并行加速和优化。CUDA是一个软件环境，它使得编程者无须掌握复杂的图像渲染技术，仅仅使用C、C++、Fortran等高级语言即可进行GPU通用计算的编程，实现程序的并行加速。CUDA设计的目的主要是支持不同的语言和应用编程接口，提高GPU编程的通用性。CUDA的核心是三个方面的抽象：线程层次抽象、存储器抽象和同步屏障的抽象[21]。这些抽象的编程模型向程序员提供内嵌在粗粒度的数据并行和任务并行中的细粒度的数据并行和线程并行。它指导程序员将一个问题地分为几个子问题，这些子问题可以并行地在不同的线程块中执行，而每个子问题本身又可以通过多个线程的并行协作得到解决，最终整个问题得到解决。22 第三章GPU硬件架构和编程模型GPU里面内置了一个流处理器的阵列，一个有大量GPU线程的程序会被划分到线程块中，所有的线程块在GPU的调度下在流处理器上依次、独立地执行，因此，GPU中的流处理器越多，程序执行的时间就越少，如图3.9所示。每一个线程块可以在任何时间以任意的顺序调度到GPU中任何一个可用的流多处理器中，这样一个被编译的CUDA程序就可以在拥有不同数量的流多处理器的GPU中执行，而只有运行时系统需要知道实际的流多处理器的数量。这样CUDA就成为一个可扩展的编程模型，它使得GPU的架构可以通过增加流多处理器的数量和使用不同存储器划分来提高GPU的性能。MultithreadCUDAProgramBlock0Block1Block2Block3Block4Block5Block6Block7GPUwith2SMsGPUwith4SMsSM0SM1SM0SM1SM2SM3Block0Block1Block0Block1Block2Block3Block2Block3Block4Block5Block6Block7Block4Block5Block6Block7图3.9CUDA程序在GPU中的执行流程在CUDA中，当在主机端调用一个kernel时，kernel中的程序将在GPU中多个不同的CUDA线程中执行。这些线程组织成一个个线程块，线程块可以是一维、二维和三维的，线程块中的线程都有自己的线程ID，同一个线程块中的线程可以通过共享内存进行通信[22]。一个线程块中可以容纳的线程是有限的，如TeslaK40显卡的一个线程块中最多可以容纳2048个线程，而很多应用中又要开启大量的线程，因此CUDA又把线程块组织成Grid。Grid也可以是一维、二维或者三维的，Grid中的每个线程块也都有自己的线程块ID。各个线程块之间都是独立执行的，不能进行通信。这是由于在GPU内部，是以线程块的形式对任务进行调度的，而在进行任务调度时23 西安电子科技大学硕士学位论文各个线程块之间执行的先后次序无法得到保证，因此各个线程块中的任务不能有任何联系。在一个Gird中，一个线程的ID与它在线程块中的ID和它所在的线程块在Grid中的ID相关，同时，还和线程块和Grid的大小有关。可以通过CUDA内建的threadIdx和blockIdx两个变量通过计算得到每个线程在Grid中的线程ID。如图3.10所示，是CUDA的线程层次结构图。GridBlock(0,0)Block(1,0)Block(2,0)Block(0,1)Block(1,1)Block(2,1)Block(1,1)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)Thread(0,0)图3.10CUDA的线程层次图当CUDA程序运行时，它可能会从不同的存储器获取数据，这些存储器的层次结构图如图3.11所示。CUDA中的每个线程都有自己的localMemory，不同线程之间不能共享localMemory。每个线程块都有自己的共享内存（SharedMemory），共享内存在一个线程块内对所有的线程都可见，所以同一个线程块中的线程可以通过共享内存进行通信。常量内存（ConstantMemory）和纹理内存（TextureMemory）是两种特殊的只读内存，他们对所有的线程都可见。全局内存（GlobalMemory）是一种对所有的线程都可见的可读可写内存。当在线程中定义localMemory时，CUDA会优先分配寄存器来充当localMemory，寄存器的存取速度和运算器的速度相同，所以它的存取速度是最快的。当寄存器不够时，CUDA会使用GlobalMemory来充当localMemory，此时用来充当localMemory的GlobalMemory只对该线程可见，由于Global24 第三章GPU硬件架构和编程模型Memory是片外存储器，所以，它的存取速度较慢，此时对本地变量存取的速度将会大大降低。共享内存是片内存储器，所以它的存取速度比寄存器慢，但比GlobalMemory快很多。常量内存和纹理内存是一种特殊的片外存储器，当GPU对常量内存和纹理内存中的数据进行存取时，GPU会对存取的数据的以及相关的数据进行缓存，因此它的速度介于共享内存和全局内存之间。ThreadPer-threadlocalmemoryPer-blockThreadBlockSharedMemoryGrid0Block(0,0)Block(1,0)Block(2,0)Block(0,1)Block(1,1)Block(2,1)Grid1GlobalMemoryBlock(0,0)Block(1,0)Block(0,1)Block(1,1)Block(0,2)Block(1,2)图3.11CUDA的存储器层次图CUDA编程模型假设CUDA线程在设备端（GPU端）执行，而程序的剩下部分在主机端（如CPU端）执行，CUDA编程模型还假设设备端和主机端都拥有各自相互独立的存储空间，分别叫做设备端存储器和主机端存储器，因此CUDA编程模型是一种异构编程模型，如图3.12所示。在一个CUDA程序中，程序员要通过CUDA运行时显式地管理包括本地内存、共享内存、常量内存、纹理内存和全局内存在内的设备端内存以及主机端内存，还需要管理数据在主机端和设备端的相互拷贝。由于，主机端的线程和设备端线程是单独执行的，所以需要有同步机制来对线程进行控制，默认情况下，主机端和设备端的执行是同步的，主机端和设备端是否同步可以通过CUDA运行时来控制。另外，由于同一个线程块中线程可以通过共享内存进行通信，而线程块中的线程的执行顺序是不可知的，这样，就会存在线程间的竞争。为了解决这个问题，CUDA运行时提供了线程间的同步机制，还提供了原子操作，这些机制都25 西安电子科技大学硕士学位论文保证了线程之间能够有序地执行，而避免竞争。SerialCodeHostDeviceGrid0ParallelkernelBlock(0,0)Block(1,0)Block(2,0)Block(0,1)Block(1,1)Block(2,1)SerialCodeHostGrid1Block(0,0)Block(1,0)ParallelkernelBlock(0,1)Block(1,1)Block(0,2)Block(1,2)图3.12CUDA异构编程模型3.4本章小结本章首先讲述了并行计算相关的基础知识，然后介绍了GPU硬件架构的历史变迁。针对本文的设计方案中使用的NVIDIATelsaK40型号的显卡，对该显卡的硬件架构和性能进行的说明，并对CUDA编程模型进行了介绍，为下文中GPU程序的设计和优化奠定基础。26 第四章基于GPU和MATLAB的机载SAR成像算法实现及第四章一种改进的RMA算法的GPU实现和优化RMA（RangeMigrationAlgorithm）是一种比较成熟的SAR成像频域算法，在宽测绘带、高分辨率SAR成像方面应用比较广泛。RMA算法基于散射点模型，在公式的推导过程中没有进行任何的近似，对SAR成像中大斜视角的情况也有很好的适应性。因此，RMA算法在原理上是最优的成像算法。但是，RMA算法在实现的过程中要通过Stolt插值来完成距离徙动校正，插值的精度直接影响成像质量，但是插值过程的计算量巨大，对于宽测绘带、高分辨率SAR成像，插值过程往往会成为实时成像的瓶颈。另外，在正侧视宽波束的情况下，由于回波信号的波数谱弯曲地比较厉害，使用传统的RMA算法在插值过程中的频谱利用率较低，使得图像的分辨率降低。本章针对传统RMA算法存在的问题，对传统RMA算法进行改进，使算法适用于正侧视宽波束的情况，并采用GPU对改进的RMA算法进行实现和优化，提高算法的处理速度，使算法满足实时性要求。4.1传统的RMA算法RMA算法是在波数域采用插值的方法完成对SAR信号的距离徙动校正，下面从点目标的基本回波信号开始，在波数域给出RMA算法的推导。如图2.4所示，载机沿X轴飞行，以合成孔径中心为坐标原点，载机飞到A点处的X轴坐标为vt,0，m场景中任一点P到载机的飞行航线的垂直距离为R，则点目标所在的坐标为BPXR,B。假设载机沿X轴匀速飞行，则对于场景中的点目标P，与载机的瞬时斜距可以表示为：22RtR;XvtR(4-1)aBaB则雷达接收到的点目标P的回波，经过去载频后得到的基频回波可以表示为：2;RtRaBStech,;tRataBrattaa0c22;RtRaBexpjt(4-2)c4expjRtRaB;其中，t为距离向快时间，t为方位向慢时间，a()为距离向窗函数，a()为方ara位向窗函数，为发射信号调频率，t为波束中心指向P点时的方位时刻，为载波0波长。27 西安电子科技大学硕士学位论文xy图4.1正侧视SAR成像几何模型对基频回波信号进行距离向傅立叶变换，将信号变换到距离频域，则基频回波信号在距离频域可以表示为：SftAfatt()(,=-)()echrarraa0æöf2´-expççjpr÷÷ççg÷÷÷èøæöç4p÷´-expçççèøjR()tRffaBr;()+c÷÷÷c(4-3)其中，其中fr为距离频率，c为光速，fc为雷达载频。对基频回波信号进行脉压处理后，(4-3)式可以表示为：SftAfatt(,)=-()()echrarraa0æöç4p÷(4-4)´-expçççèøjRtRf()aBr;()+fc÷÷÷c采用驻相点原理将(4-4)式变换到方位频域，可以得到基频回波信号的二维频谱，忽略幅度项，则基频回波信号在二维频域的表达式为：22fffSffR,;exp4jRrcaechraBBcv2Xexpjf2av(4-5)令Kf4fc，Kf2v分别表示距离波数和方位波数，将其代入式rrcxa(4-5)，则式(4-5)可以表示为：SKK1()rx,e=-xp()jRKKBr22--xjXKx(4-6)28 第四章基于GPU和MATLAB的机载SAR成像算法实现及其中KKrx,的取值范围可以近似为：4(fB2)4(fB2)cBWcBWKsin,sinxcc22(4-7)44Kf(2Bf),(2B)rcccc其中，B为雷达发射信号的频带宽度，为波数宽度。结合式(4-7)可以得到式BW(4-6)的波数谱支撑区，如图4.2所示。对基频回波进行数字采样后得到一系列离散采样点数据，其波数谱在KK平面内是一组等间隔的点。在波数平面，距离波数Kxrr222可以分解为K和K两个分量，即KKK，则xyrxy22KKK(4-8)yrx将式(4-8)代入式(4-6)可得SKK1()(rx,e=--xpjRKjXKByx)(4-9)则二维波数谱在KK平面的支撑区如图4.3所示。xyKRKx图4.2基频回波的二维波数谱在KK平面的支撑区xrKyKx图4.3基频回波的二维波数谱在KK平面的支撑区xy由于，K和K在支撑区内是均匀的，则K是不均匀的。因此，不能直接对(4-9)rxx29 西安电子科技大学硕士学位论文式直接进行二维傅立叶变换使图像聚焦，需要在支撑区内将K在支撑区内均匀分段，y根据式(4-8)出相应的K的值，再通过插值得到相应点的波数谱的值。K在ry4KK,之间均匀取值，一般KKfB，yyminmaxyrminmincc24KKyrminmaxfcB2为固定值。对于不同Kx出的值，均将Ky在KKyymin,maxc之间均匀取值，再通过Stolt插值求得相应的波数谱。这样，插值后的数据在K和Kxy维都具有均匀的分布，即插值后的二维波数谱在KK平面具有矩形网格的形式，xy所以直接对插值后的信号进行二维FFT处理，或者两个一维的FFT处理，得到聚焦的图像。4.2一种改进的适用于宽波束大场景的RMA算法在正侧视波束较宽情况下，图4.3对应的扇形区域弯曲的很厉害。若还是采用传统的RMA算法，插值后的(KKxy,)取值范围还是矩形。若矩形选区小了，对于扇形区域的左右边缘部分，部分支撑区内的点会取不到，这样波数谱利用率不高，会在距离维损失部分分辨率；若矩形选大，完全覆盖扇形区域，会导致插值点数变多，影响插值效率。因此，需要对传统的RMA算法进行改进。Ky22KKKyr00xKx图4.4K取值范围的选取y由于K取值范围的固定，使得插值后(KKxy,)所在的矩形区域无法满足要求，y因此，可以对K的取值范围进行改进。如图4.4所示，可以以K为中心，在K两yy0y0224fc侧对K进行均匀取值，这里KKK,K。yyr00xr0c对于基频回波信号，对于每一个K取值单元，求出回波信号二维波数谱支撑区x中K的最大值和最小值。K的最大值为KKK，这里yyyrmaxmaxx30 第四章基于GPU和MATLAB的机载SAR成像算法实现及4fcB24fcB2K。K的最小值为KKK，这里K。rmaxyyrminminxrmincc通过K、K和K计算出K取值范围的大小。令ymaxyminy0yKKmaxyymax00mK,KyKyin，这里是取绝对值符号，则KyKKyy0p，其中K在K,K之间均匀取值。这样，对于每个K取值单元处，可以得到：ypx22KKKKK(4-10)yyy0prx则22KKKK(4-11)ry0ypx根据式(4-11)可以计算出每个K取值对应的K取值，再通过sinc插值就可以计yr算出各个K取值点处的波数谱的取值。y经过以上的处理后，式(4-9)变为SKK,expjRKKjXK(4-12)10rxByypx对式(4-12)做逆傅里叶变换，变换到距离时域，可以得到2RB22SRK,etxpjRKKjXK(4-13)20BxBrxxc这样，基频回波的距离向就得到聚焦，对式(4-13)中的剩余相位项进行补偿，补偿函数为：22HRK,expjRKK(4-14)BxBr0x补偿后的回波信号的形式如下：2RBSRK,etxpjXK(4-15)2Bxxc对式(4-15)进行逆傅里叶变换就可实现方位聚焦。4.3算法的GPU实现和优化方案4.3.1总体方案设计根据上一节的讨论，我们可以知道上述算法对RMA算法的改进主要集中在Stolt插值。改进的RMA算法的实现流程图如图4.5所示。GPU程序的设计必须首先根据数据和算法的特点，对数据和任务进行并行划分，明确数据的分配策略和以及确定算法中哪些部分需要并行，哪些部分必须串行，得到最优的程序总体设计方案。31 西安电子科技大学硕士学位论文读取数据距离向FFT转置方位FFT参考函数相乘/Stolt插值/转置距离向IFFT/转置方位向IFFTSAR图像图4.5改进的RMA算法的实现流程由图4.5可知，除了读写数据的操作，整个处理过程中所有的操作都有一定的并行性。因此，在数据量较小的情况下，可以考虑在主机端（CPU端）将数据读取完成后，直接将数据拷贝到设备端（GPU端）进行处理，在完成整个处理过程完成后，再将数据复制回主机端。但是，这种主机端只负责数据读取，设备端负责成像处理整个过程的方案需要一次性在GPU中开辟大量的临时空间。对于数据量小的情况，显卡的显存可以满足整个处理过程中的内存开销，但是对小数据量的数据进行处理往往不能获得很高的加速比。而且在实际的高分辨率的SAR实测数据的处理中，数据量往往都很大，显卡的显存不足以一次性读入所有的数据。例如，TeslaK40显卡的显存是12G，而这12G并不能全部都用来存储读入的数据，必须留出大量的存储空间用来存放临时数据，实际可以使用的只有几个G，而对于宽幅高分辨率SAR的数据量往往有几十个G，甚至上百G，因此显卡的显存难以满足实际数据处理的需要。另外，如果一次性读入数据，将无法使用CUDA中流的机制对数据处理进行加速，这样数据处理的效率将会降低。考虑到以上因素，本文采用将需要处理的回波数据进行分块处理的策略，设备端每次从主机端读入一块数据，处理完毕后复制会主机端，CPU主要负责整个处理流程的控制和数据存取的工作，GPU主要数据处理，同时采用CUDA流技术，减少内存拷贝带来的时延。4.3.2数据并行划分在实际SAR回波信号的获取过程中，SAR信号都是以矢量的形式进行存储的，SAR信号的处理往往也都是以一个矢量为基本的处理单元，且有大量随机读取数据的操作。根据SAR信号的这些特点，本文采用线性存储结构对读入的回波数据进行32 第四章基于GPU和MATLAB的机载SAR成像算法实现及存储。SAR信号是二维信号，在将其线性存储的过程中，需要以某一维为主进行数据的重排，然后再将其分段写入设备端。由图26可知，在整个算法的处理流程中有大量的矩阵转置操作，在采用线性存储的情况下，矩阵转置操作需要对数据进行重排，转置操作存在大量的并行性，在数据量不大，为数据开辟新的存储空间的情况下，可以使用GPU对矩阵转置操作进行加速，此时整个转置操作都在设备端完成。对于数据量比较大的情况，由于数据无法一次性读入设备端，使得GPU无法对整个矩阵完成转置操作，只能对矩阵的部分数据进行转置，而采用GPU对矩阵的部分数据进行转置并将数据拷贝的CPU端后，仍需要对这些数据进行重排，计算量仍然很大。因此，传统的做法是将矩阵转置的操作在CPU端完成，而在数据量大的情况下，CPU对矩阵转置操作是极其低效的。针对大数据量数据转置的传统做法存在的问题，本文根据分块矩阵的思想，采用一种改进的方法，在不增加内存拷贝次数的基础上，在GPU端完成矩阵转置的操作。本文采用按列线性存储的方式对矩阵数据进行存储，因此，每次转置就相当于对整个线性序列进行重排。假设矩阵M为需要转置的大型矩阵，由于GPU一次无法处理整个矩阵中的全部数据，因此，需要将M矩阵分块，将分块数据复制到GPU端完成转置处理后再拷贝的CPU端。为了使矩阵M的各个子块在GPU端完成转置后，能够在CPU端直接复制到相应的存储单元完成相应的转置操作，而无需再进行其他操作，这里采用按行分块的方式将矩阵M分块。设矩阵M按行分块后可表示成如下形式：ABM(4-16)C...则根据分块矩阵的性质，矩阵M转置后可表示为：TTTMABC...(4-17)这样，经过GPU转置后的各个分块矩阵就可以在按列存储的方式下，无须CPU端的其他操作，直接顺序排列在主机的线性内存中，完成整个矩阵的转置操作。另外，这种以分块的方式在GPU中进行矩阵转置的方案，可以在GPU中完成矩阵转置处理后直接进行后续的处理，处理完毕后再将完成转置并处理后的数据复制到主机端拼接成完整的矩阵，没有增加额外的内存复制操作。由于数据采用线性存储的方式，且在对矩阵进行存储时采用按列排列的方式，所以，在对原始矩阵按行分块并复制到GPU端时，如果采用常规的方法将会存在困难，所以，在实际实现中将采用CUDA中的2D线性内存，使用cudaMemcpy2D()函数来完成数据在主机端和设备端的复制。在CUDA中2D线性内存仍然采用一维存储的方式，但在进行数据拷贝时可以根据pitch参数确定数据所在矩阵的列数，根据width33 西安电子科技大学硕士学位论文参数确定实际需要拷贝的数据长度。这样，就解决了在按列线性存储的矩阵中按行选取数据的问题。在将数据分块复制到GPU端后，就需要考虑GPU中线程网格和线程块划分问题。在GPU的并行处理中，有两个维度的并行：线程网格的粗粒度并行和线程块的细粒度并行。在并行计算中，粒度是并行处理中计算量和通信量的比值，是衡量并行计算性能的一个重要标准。在线程网格中，大量的线程块集中于处理各自的计算任务，各个线程块之间无法进行通信，GPU中的调度器以线程块为单位对各个线程块进行调度，根据GPU中的线程块调度机制，GPU中每个CUDA核心至少要有两个线程块才能使每个CUDA核心中的计算计算资源充分利用，而GPU中有大量的CUDA核心，例如，TealaK40中有2880个CUDA核心，因此，线程网格中的线程块的划分不宜过少。线程网格的各维的维度可以采用下列的公式来计算：Datax__BlockDimx1Dimx_BlockDimx_(4-18)Data_yBlockDim_y1Dim_yBlockDim_y其中，Dimx_表示线程网格x轴的维度，Dim_y表示线程网格y轴的维度，Datax_是输入数据的x轴的维度，Data_y是输入数据的y轴的维度，BlockDimx_是线程块的x轴的维度，BlockDim_y是线程块的y轴的维度，式(4-18)中之所以要减1是为了使数据的边界处也得到处理。从上式可以看出，在线程块大小确定的情况下可以通过调整数据分块大小来改变线程网格的大小，使之与相应的设备相适应，完成性能调优。但是，在线程块大小确定的情况下，实际的线程网格大小还受制与GPU的显存大小，因此，必须结合GPU的显存大小才能确定线程网格的大小，具体的计算方法在下文给出。4.3.3GPU存储优化在NVIDIA的GPU中，一个CUDA核中的一条指令一次能广播给一个warp中的所有线程，warp的大小是有NVIDIA确定的，目前，在NVIDIA所有的产品中，warp的大小都为32。因此，一个线程块中的线程大小最好是32的倍数，这样能提高CUDA核的指令执行效率。如图4.6所示，是每一个CUDA核中的warp占用数量随线程块大小的变化曲线图。从图中可以看出在线程块大小为256、512和1024时，warp的占用数量都达到了最大，但是，在TeslaK40中，一个线程块中至多有1024个线程，因此，如果线程块大小为1024，则一个CUDA核中至多有一个线程块，而在一个CUDA核中至少要有两个线程块进行任务调度才能是CUDA中的计算资源不至于闲置，所以，一般取线程块大小为256或者512，具体选哪个还要根据内存使用34 第四章基于GPU和MATLAB的机载SAR成像算法实现及情况来具体确定。线程块的x轴的维度最好是16的倍数，这样，在对全局内存的访问时内存能够对齐，使得线程块中的线程对全局内存的多次访问可以合并，从而提高内存的存取速度[23]。图4.6线程块大小对CUDA中warp占用数量的影响原始矩阵的子块矩阵复制到GPU端后，需要在GPU端对子块矩阵进行转置操作。在GPU端进行矩阵转置的基本思想是在显存中开辟与子块矩阵相同大小的存储空间，建立一个零矩阵，然后开启大量线程根据转置后矩阵和原矩阵中元素位置的对应关系，将新建立的零矩阵的各个位置对应的数据填入。但是，这种从全局内存读取数据到各个线程中执行，然后直接写回全局内存的实现方式仍然非常低效。因为，影响GPU程序执行效率的一个重要因素是存储带宽，GPU的存储带宽的大小远远小于GPU运算的速度，例如，TeslaK40的峰值双精度浮点性能可以达到1.66Tflops，而它的存储器带宽却只有288GB/sec，两者差了将近两个数量级。因此，GPU对全局内存的访问是非常低效的，必须对它进行优化。由上文可知，划分的线程块在x维的大小是16的倍数，而在一个half-warp中所有线程，如果这些线程GPU在对全局内存进行访问时，访问的是内存首地址对齐的段，即内存首地址是32的倍数的连续的一段数据，GPU就可以将这16次对内存的访问合并为一次访问，从而大大减少访问内存的时间，可以根据GPU的这一特点，对程序进行性能调优。在上述GPU实现矩阵转置的方案中，当同一个线程块中的线程对全局内存中的子块矩阵数据进行读取时，由于满足合并读取的条件，因此同一个half-warp中的16线程对全局内存的16次访问可以合并为一次访问。但是，当每个线程处理完后向全局内存中写入数据时，由于同一个段的数据分别位于不同的half-warp中，因此不能满足全局内存合并方位的条件，这样一个half-warp对全局内存的一次访问将分成1635 西安电子科技大学硕士学位论文次独立的访问，大大降低访存的效率。在本文的方案中，我们采用共享内存的方法对该问题进行了解决。具体做法是，在每一个线程块中开辟一个与读取数据同等大小的共享内存，当一个线程块中的线程将数据处理完成后，并不直接将数据写回全局内存，而是写入共享内存中，然后使用CUDA中的__Synchronize()函数对所有线程进行同步，这样当数据在写入全局内存时就能满足合并访问的条件，将同一个half-warp中的16次写入操作合并为一次写入操作，提高了程序的执行效率。由于在GPU中对子块矩阵进行转置时要使用共享内存，因此在确定线程块中的线程大小时将要考虑共享内存的大小。由第三章可知，GPU中的SMX中有64KB的L1缓存，可以配置为48KB的共享内存和16KB的缓存，也可以配置为16KB的共享内存和48KB的缓存。这里，我们将它配置为48KB的共享内存和16KB的缓存，这样同一个线程块中可以有48KB6K个线程，远大于1024个线程的限制，结合之8前的分析，可以确定本方案中线程块中的线程个数确定为512，线程块x轴的维度为16，y轴的维度为32。在数据量巨大的情况下，将数据分段复制到GPU中处理的方法可以解决显卡中显存不够的问题，但是随之带来了大量的主机端和设备端内存拷贝操作。因为，CPU端和GPU端的数据传输效率远远低于全局内存的速度，所以大量的内存复制操作也将影响程序的性能[24]。本文在在对数据进行分段的前提下，采用CUDA中的流技术对上述问题进行解决。CUDA中的流可以看出GPU中的一个有序的操作序列，提交到每一个CUDA流中的任务，在满足一定条件的情况下，流都会按照任务提交的顺序对任务逐个执行，不同的流之间存在一定的依赖关系。不同流之间的核函数不能同时执行，但是，在一个流执行核函数时，其他的流可以执行主机和设备端的内存复制操作。要使用CUDA中的流技术，GPU必须要支持设备重叠功能，这些GPU中拥有一个执行核函数的引擎，同时还拥有一个或者多个执行主机和设备端内存拷贝的引擎，两个引擎各自拥有自己的任务执行队列。因此，在这些GPU中执行核函数的同时可以执行内存复制的操作。在将回波数据进行分块后，就可以将数据提交给各个不同的流进行处理，每个流中必须开辟自己独立的存储空间。由于使用流可以同时进行和核函数执行、从CPU端向GPU端复制数据和从GPU端向CPU端复制数据，因此，开启三个流后，将能够保证在数据读取完毕后，所有流在任何时刻都有任务在执行，本文中将采用三个流对回波数据进行处理，如图4.7所示。要采用流技术，必须要使用cudaMallocHost()函数在主机端为数据分配页锁定主机内存，使得在整个程序的执行期间，数据不会被交换出内存，而使用malloc()函数分配的主机内存是可分页的，操作系统会在内存紧36 第四章基于GPU和MATLAB的机载SAR成像算法实现及张的时候将其交换出内存，而当它需要被使用而换回内存时，它的存储地址可能已经改变，所以在进行数据拷贝时可能会出错[25]。在不使用流技术的情况下使用页锁定主机内存也能够提高程序的执行效率，因为GPU在主机端和设备端进行内存拷贝时，是现将可分页内存中的数据拷贝到页锁定内存中，再从页锁定内存中把数据拷贝到设备端，而直接使用页锁定主机内存就省去了上述的第一个步骤，提高了执行速度。第0个流第1个流第2个流CPU端拷贝数据到GPU端CPU端拷贝数据到数据处理时间GPU端GPU端拷贝到CPUCPU端拷贝数据到数据处理端GPU端GPU端拷贝数据到…数据处理CPU端图4.7使用3个流时的任务执行时间线在成像的过程中，在使用大量的FFT和IFFT操作，这些操作通过CUDACUFFT函数库中的相关函数来完成。在使用CUFFT函数库进行FFT操作时，需要使用cufftPlan1D函数来进行来创建一个执行FFT的任务，而cufftPlan1D函数在创建任务时会开辟与需要处理的数据同等大小的临时空间以及一些配置数据，因此，在进行数据划分和内存分配时必须为FFT的操作预留存储空间。另外，由于开启了流操作，所以在执行FFT操作时必须要用cufftSetStream()函数将每一个FFT的plan与对应的CUDA流关联起来。在GPU处理数据的过程中，显存的分配和释放会占用一定的时间。在本方案中，由于GPU要处理大量的数据块，如果每处理一块数据就分配和释放一次内存，将会耗费大量的时间，影响系统的实时性。因此，在本方案的3个CUDA流中，对于每一个CUDA流，只在建立CUDA流时分配一次显存，在整个处理过程中循环使用分配的内存，直到整个处理流程结束后才释放显存，如图4.8所示。这样，就可以避免频繁地分配和释放显存。37 西安电子科技大学硕士学位论文任务序号s=0显存分配参数读取成像处理任务存在s=s+1任务结束显存释放终止图4.8分配和释放一次显存的方案4.3.4程序实现步骤在进行成像处理之前，需要先计算出数据分块的具体方案。设原始数据的大小为M，距离向采样点数为N，方位向采用点数为N，每个分块子矩阵的大小为M，raD则在主机端需要分配M=M的存储空间。在设备端，由于要开启三个CUDA流，因H此，每个CUDA流必须都独立的分配存储空间，在每个流中，除了要为每个子矩阵块分配存储空间外，由于转置操作和FFT操作需要额外的存储空间，所以，在每个流中还必须开辟大小为每个子矩阵块大小2倍的临时存储空间。本文采用的显卡型号是TeslaK40系列，其显存大小是12G。根据以上分析可以得到：91M2GB(4-19)D即M1.5GB，因此，单个子矩阵块的数据量大小必须小于1.5GB。所以，原D始数据所分的块数为MMD1K(4-20)MDM按方位向取分块数据时，分块数据的距离向脉冲点数为D，按距离向取分NaM块数据时，分块数据的方位向脉冲点数为D。Nr基于以上分析，我们可以知道，要用GPU实现改进的RMA算法，并使得程序能够处理大数据量实测数据，需要以下几个步骤：（1）在主机端使用cudaMallocHost()函数为原始回波数据分配页锁定主机内38 第四章基于GPU和MATLAB的机载SAR成像算法实现及存。（2）确定分块子矩阵的大小，根据以上计算子矩阵的数据量必须小于1.5G，在本方案的设计中，取M512M。D（3）在GPU中开启3个CUDA流，按距离向将数据分块，每一个分块数据的数据量大小为M。依次将各个分块数据分配到3个CUDA流中对子块中的距离向数D据进行FFT,然后执行距离向脉冲压缩处理，处理完毕后将数据拷贝回主机端。（4）按方位向将数据分块，分块大小为M，使用cudaMemcpy2D()函数将分块D后数据拷贝的3个流中依次进行转置、方位FFT、参考函数相乘、插值等一系列操作，处理完成后将数据拷贝回主机端。（5）按距离向将数据分块，分块大小为M，拷贝数据到对应的CUDA流，依D次执行转置和距离向IFFT操作，处理完成后将数据从设备端拷贝到主机端。（6）按方位向将数据分块，分块大小为M，拷贝数据到对应的CUDA流，依D次执行转置和方位向IFFT操作，处理完成后将数据从设备端拷贝到主机端，完成成像过程。4.4实验结果及分析为了测试本章中算法实现的性能，本节采用L波段的一段机载SAR回波的仿真数据，分别使用CPU实现算法的方案和上文实现的CPU+GPU实现算法的方案对该段数据进行处理，得到处理后的图像，并测量两者的性能进行比较。测试时，对数据进行处理的的运行环境和仿真参数如下：表4.1运行环境和仿真参数CPU酷睿i7-4790（一颗）GPUNVIDIATeslaK40c（一部）主机内存128GB操作系统MicrosoftWindows7（x64）SAR仿真数据字节数3GBsinc插值点数16测试时，CPU程序同时开启16个线程，测试共进行3次，测试结果如下表所示：表4.2两种方案耗时对比CPU方案耗时CPU+GPU方案耗时加速比319.728s21.75s14.70321.41s22.43s14.33321.72s22.71s14.17从测试结果中，使用GPU程序进行处理时，处理速度明显加快，三次测试得到39 西安电子科技大学硕士学位论文的GPU加速比分别是14.70、14.33和14.17，即达到了14倍左右的加速比。4.5本章小结本章针对机载宽幅SAR成像中合成孔径时间长，二维波数谱弯曲度大的问题，通过对传统RMA算法的Stolt插值方法进行了改进，解决了回波二维波数谱的利用率低的问题。使用GPU编程对算法进行了实现，考虑到实测数据的数据量大的问题，对程序的数据并行划分、线程划分和存储等方面进行了优化，并使用一组仿真数据对程序性能进行了测试，测试结果表明，相对于CPU程序，GPU程序能实现14倍左右的加速。40 第五章基于GPU和MATLAB的SAR成像算法实现及优化第五章基于GPU和MATLAB的SAR成像算法实现及优化对于机载SAR而言，由于气流的扰动和载机速度的不稳定性，雷达载机会偏离理想的运动轨迹，录取的回波数据存在运动误差，使得SAR图像会出现严重的散焦和几何失真[26]。对于宽波束的SAR成像系统而言，由于合成孔径时间长，这个问题更加严重，因此，在成像处理之前，需要对回波数据进行运动补偿。相位梯度自聚焦算法（PhaseGradientAutofocus，PGA）能对运动误差进行很好的估计，但是，一般情况下，相位梯度自聚焦算法都是基于全孔径信号对运动误差进行估计，而机载宽波束高分辨率SAR成像的数据量很大，如果基于全孔径信号对运动误差进行估计，就无法满足快速实时成像的要求。本章采用一种基于子孔径的机载宽波束高分辨率SAR成像处理方法，并采用GPU并行处理技术对该算法的耗时部分进行优化加速，然后使用MATLAB和C混合编程技术对算法进行实现，使其可以实现快速成像的要求，最后对程序的性能进行了测试，实验结果表明该算法的实现方案能对实测数据进行良好的成像并实现9倍左右的加速比。5.1一种机载宽波束子孔径SAR成像算法5.1.1存在线性误差时时聚焦结果的推导根据图4.1中的几何模型，雷达接收到的点目标P的回波，考虑运动误差，并忽略幅度项，经过去载频处理后得到的基频回波，可以表示为：22;RtRaBRtaStech,;tRaBexpjtc(5-1)4expjRtaB;RRta其中，t为距离向快时间，t为方位向慢时间，为发射信号调频率，为载波a波长。将式(5-1)变换到距离频域，并进行脉冲压缩处理后，可得：4ffcrStech,;tRaBexpjRtRRta;Ba(5-2)c其中22RtR;RvtX(5-3)aBBa23Rtaatatat(5-4)aa012a3a这里，RtR;表示载机在理想的飞行轨迹下，雷达与点目标P的瞬时斜距。aBRta是载机在非理想轨迹下的斜距误差，这里我们用多项式对斜距误差进行近似，41 西安电子科技大学硕士学位论文多项式的系数为a(i0,1,2,)，v表示载机的理想飞行速度，f表示距离频率，c表ir示光速。对于子孔径数据，相位梯度自聚焦（PGA）算法能对运动误差中的二次项和二次项以上的高次项进行较好的估计。但是，对于运动误差中的常数项和线性项却不能得到良好的估计，因为PGA是根据SAR图像中的一些特显点散焦的情况进行运动误差的估计的。因此，式(5-2)中的二次项误差和二次项以上的高次误差可以得到良好的补偿，还剩余常数项误差和线性项误差。式(5-2)经过高次项补偿后，可以表示为4ffcr22Sftra,expjRvtXBac(5-5)4ffcrexpja01atac如果直接将常规的成像算法直接应用于式(5-5)，由于常数项误差a和线性项误差0a的存在，将会使点目标P出现散焦。因此，需要计算出常数项误差a和线性项误10差a的值，从而对运动误差进行补偿，得到点目标精确的聚焦位置。将式(5-5)变换到1距离波数域可以得到22axSKx,expjKRxXjKa1(5-6)10rrBrv4fcrf其中，xvt为载机的方位向坐标，K为距离波数。然后，将式(5-6)arc由方位空域变换到方位波数域，得到信号的二维波数谱22SKK,expjRKKKjXKKjaK(5-7)rxBrxx00xx0raK1rK2fa其中，K，为方位波数。这里，采用第四章中的改进的x0vxvRMA算法对式(5-7)进行处理，并推导出带有常数项和线性项运动误差的聚焦结果。将式(5-7)与一个参考函数相乘，以匹配场景中心点目标的回波相位，参考函数为22HKK,expjRKKjRK(5-8)0rxsrxsr将式(5-8)与式(5-7)相乘，并使用第四章的改进的Stolt插值的方法对式(5-7)进行Stolt插值，插值表达式可以表示为2222KKKKK(5-9)rxrxy0插值完成后，信号的表达式变为42 第五章基于GPU和MATLAB的SAR成像算法实现及优化221a1KKKK222rxyx0vSKK1yx,expjRBa2122222KKKKKKxr0xyxxva2(5-10)1222expjKKKKKXxr0xyxv2222expjKKKKaRrxyx00s22expjRKKKsrxy0由于，K的二次项和二次项以上的高次项都很小，所以在式(5-10)可以忽略它。y这样，式(5-10)可以变为aa22112SKK20yx,expjRAKBr2KKKjKKxr0xxr0Xvv22expjKrs00aRjRsKr0KxaKK222221rx0AKKKrxx0vKr0expjRKBy(5-11)22a12AK2KKKrxrx00vKK22arx01jXKyKvr0exp22KKrx0jaR0sKysjRKyKr02A1a1这里，。式(5-11)中K的第一项表示点目标的实际聚焦位置与精vy确聚焦位置在距离向的偏移量。得到式(5-11)后，对式(5-11)的距离向做IFFT将信号变换到距离空域，变换后的信号表达式为43 西安电子科技大学硕士学位论文aKK222221rx0AKKKrxx0vKr0rRB22a12SK2r,xAKrxrx002KKKv2222KKaKKrx001rxXaR0ssRKvKrr00aa22112expjRAKKKKjKK2X(5-12)Brxrxx00vvr022expjKrs00aRjRsKr0Kx这里，是距离脉压后得到的距离包络。由式(5-12)可以看出，距离包络是方位波数的函数，即距离包络中仍然存在距离单元徙动（RCM）。但是，由于这个误差很小，所以在本章的推导中将它忽略。这样，式(5-12)表示的距离脉压后的信号可以写为a2212jRAKKKK2Brxrx00SKr,rRexpv2xeqnajKK1X(5-13)xr0v22expjKrs00aRjRsKr0Kx其中2aa11RRX1a(5-14)eqnB0vvins式(5-13)中，有一个二次相位项需要补偿。对于距离单元R处的点目标，相位eqn补偿函数可以构造为22HexpjRRKKjRRK(5-15)Mateqnsr00xeqnsr将式(5-13)与式(5-15)相乘，并将相乘结果的相位项展开成K的泰勒级数，则式x(5-13)可以近似表示为SrK20,exerRqnexpjRKqnrK2(5-16)xexpjXKjRRxeqn2Kr0其中RaB1XXAv(5-17)RRB3A44 第五章基于GPU和MATLAB的SAR成像算法实现及优化为了避免方位混叠，还需要对式(5-16)进行去斜操作，去斜函数为2KxHjexpR(5-18)Rampeqn2Kr0去斜操作完成后，将信号变换到方位空域，则信号可以表示为Kr02srRexpjRKexpjxX(5-19)30eqneqnr2R在式(5-19)中，有一个二次相位项需要补偿。这里，构造一个deramp函数，对信号进行deramp处理，使信号聚焦。构造的deramp函数为Kr02Hjexpx(5-20)deramp2Reqn将式(5-20)与式(5-19)相乘，并将相乘的结果变换到方位波数域，进行方位脉压后，信号可以表示为4XKr02sr(r,K)RKexpjRKjX(5-21)4xeqnxeqnr0RR2从式(5-21)可以看出，信号已经在二维平面聚焦了，而且点目标聚焦的位置和相位由R、X和R三个参数决定。eqn5.1.2线性误差项的估计由上述推导可以知道，运动误差使得点目标成像的位置偏离了理想的位置，但是我们能够点目标的位置和运动误差的关系，如果能对运动误差进行良好的估计，这样就使得我们可以连续地对每个子孔径数据进行快速成像，得到聚焦良好的子图像。因此，对于子孔径成像来说，对斜距误差中的常数项a和线性项a的估计是非常关键的。01下面讨论斜距误差中的系数a和a的估计。01在对子孔径数据进行处理的过程中，可以通过多次PGA估计获取每一子孔径数据的二次及高次运动误差。而在将全孔径数据划分为多个子孔径时，各个子孔径之间有部分数据重叠，如下图所示，阴影部分为子孔径间重叠区域。图5.1子孔径划分示意图图5.1为子孔径数据划分示意图。通过PGA估计得到的二次或者高次运动误差45 西安电子科技大学硕士学位论文比较准确，对于重叠部分区域，前后两次子孔径估计结果只存在线性运动误差的差别。在子孔径k中，重叠部分采用PGA估计得到的相位误差可以表示为：kktaatt(5-22)ka01aakk其中，t为二次及高次项运动误差对应的相位误差，a，a为对应常数项和a01一次项系数。在第k1个子孔径中重叠部分估计得到的相位误差为：kk11taatt(5-23)ka101aa假设全孔径信号的相位误差为t，t对时间求导可以减少常数项的影响，aa''用t表示，在全孔径内t应该是连续的，因此，可以通过重叠部分估计出一aa'次项，对t的积分结果与t相比只差一个常数。对t和t求导，可aakaka1以得到''ktat(5-24)ka1a'1k'tat(5-25)ka11a两式相减，可以得到'''kk1tttaa(5-26)akak11a1'以子孔径k为基准，只需在式(5-25)的基础上减去t，相位求导后的数据在a子孔径k和子孔径k1上是连续的，对每一重叠部分进行处理，就可以得到相位求导式在全孔径时间内的连续表达式''tat(5-27)aa1''对tt进行积分有aa'''ttaadtaat1attaa(5-28)'通过线性拟合可以得到积分项的线性项at1aat，式(5-29)减去线性项，就可以得到二次及高次项线性误差t。再从PGA估计得到相位误差中减去t就aakkkk可以得到aat，通过拟合就可以得到a和a。01a01对子孔径数据进行运动补偿，并得到聚焦的子图像后，受制与子孔径信号的方位向带宽，子图的方位分辨率没有全孔径图像的分辨率高。在将子图拼接成全孔径图像时，可以借鉴快速后向投影算法（FastBackProjection，FBP）的思想，对于大场景中的每一点，在子图像中找对应的聚焦点目标，对子孔径信号补偿相位后进行相干叠加，随着子孔径信号不断地相干叠加，网格图像的方位分辨率越来越高，最终得到高分辨的大场景图像。整个算法的设计流程如图5.2所示。46 第五章基于GPU和MATLAB的SAR成像算法实现及优化aat01a图5.2算法流程图5.2MATLABCUDAC混合编程技术5.2.1MATLABCUDAC混合编程技术简介在SAR成像领域，MATLAB是进行SAR成像算法仿真和实现的重要编程工具。MATLAB以矩阵为基本的数据处理单位，拥有大量与信号处理相关的函数库和强大的绘图功能，对于SAR成像算法的开发和仿真非常方便，能大大提高算法开发的效率，因此，SAR成像算法的开发领域，MATLAB的使用非常广泛。但是，由于MATLAB是解释型的脚本语言，使得MATLAB程序的执行效率非常低下。尤其是在MATLAB程序中有大量密集的计算和循环时，由于每条MATLAB语言的执行都要进行解释，就会使MATLAB程序执行的性能大大下降，甚至会比C版本的程序慢几倍。而GPU中有大量用于计算的运算单元，非常擅长处理计算密集型的问题，当任务有良好的并行度时，对任务进行并行划分后，采用GPU进行加速往往能够得到很好的性能。对于MATLAB程序中的计算密集的部分可以经过并行划分后直接交付GPU来处理，而程序中的循环部分也可以通过寻找循环中的并行性，将循环分解后再用GPU进行处理，MATLAB的主程序则主要进行简单的流程控制以及一些数据交互的操作。因此，采用MATLAB和GPU混合编程技术对SAR成像算法进行仿真和实现，可以将GPU的高性能和MATLAB的易用性结合起来，进一步提高成像算法的开发效率。在市场上，已经有了一些比较好的MATLAB调用GPU的开发工具包，如Jacket47 西安电子科技大学硕士学位论文工具包，但是这些工具包使用不够灵活，不能对GPU中的内存划分和内存分配进行控制，而这两个方面真是对GPU程序加速性能影响最大的两个方面，因此，这些工具包的处理效率并不高，另外，当要处理的数据量较大时，这些工具包不能解决显卡内存不够的问题，使得它的应用范围也受到限制。所以，本章采用GPU和MATLAB混合编程技术，直接对程序进行实现。GPU的编程语言有多种，如NVIDIA的CUDAC和开源的OpenCL，本文所采用的GPU开发语言是NVIDIA的CUDAC。由于CUDAC是C语言的扩展，与C语言完全兼容，因此，MATLAB和CUDAC混合编程时，可以利用MATLAB提供的C/C++编程接口进行操作。MATLAB和C/C++混合编程的方式有多种，如MATLAB调用C/C++程序，C/C++程序调用MATLAB程序，C/C++程序调用MATLAB引擎等，本文主要用到了MATLAB调用CUDAC程序，所以这里重点对MATLAB调用C/C++程序技术进行介绍。5.2.2MATLAB调用CUDAC程序的实现方法要实现MATALB调用C/C++程序，首先必须要解决数据交互问题。由于MATLAB和C/C++都有自己的数据类型，而MATLAB的各种数据类型又都有自己的底层数据结构，因此，就要进行MATLAB和C/C++之间数据类型的转换。MATLAB中的数据都是以矩阵的形式进行组织的，从广义上来说，MATLAB其实只有阵列这一种数据类型，阵列中的数据类型各不相同就形成了不同类型的阵列[27]，如图5.3所示是MATLAB所有的阵列类型。在MATLAB和C/C++的混合编程接口中定义了一种mxArray类型，当MATLAB调用C/C++程序时，先将要处理的所有数据都转换转换成mxArray类型，C/C++程序在执行的过程中调用MATLAB提供的一系列API函数将MATLAB输入的mxArray类型数据转换成对应的C/C++数据类型，这样就实现了MATLAB和C/C++程序的交互。操作mxArray类型的API函数可以参见MATLAB的开发文档。阵列逻辑型字符型数值型元组类型结构体类型Java类函数句柄整型单精度型双精度型图5.3MATLAB数据类型MATLAB调用C/C++程序是通过MEX文件的形式来实现的。当编写完C/C++程序后，采用MATLAB的编译器将C/C++程序编译成MEX文件，这样，MATLAB48 第五章基于GPU和MATLAB的SAR成像算法实现及优化就可以通过调用MEX文件的形式对C/C++程序进行调用。在Windows平台下，MEX文件是用动态链接库（DLL）技术来实现的，它本质上是一系列动态链接库文件的集合。实际上，早期的MATLAB就是直接采用动态链接库文件来实现MATLAB调用C/C++程序的。在MATLAB中使用mexfilename.c命令就能够将C/C++程序编译成可以被MATLAB调用的MEX文件，由于MEX文件有固定的结构，所以C/C++程序的编写也必须遵循指定的结构。在MEX文件中，有一个mexFunction()函数，当MATLAB调用MEX文件后，程序将从mexFunction()的入口处开始执行，直至程序执行结束后返回到调用MEX文件处，继续执行后面的MATLAB程序[28]。mexFunction()函数的返回值为空，并拥有四个固定的参数，如表5.1所示。mexFunction()函数就相当于MATLAB程序和C/C++程序交互的接口，在编写用于MATLAB调用的C/C++程序时，每个C/C++程序中都必须有一个mexFunction()函数，通过mexFunction()函数可以实现C/C++程序和MATLAB的数据交换，在mexFunction函数的内部可以编写程序具体的执行逻辑。当MATLAB调用mexFunction()函数时，会从输入参数数组中读取输入数据进行处理，处理完毕后将结果写入输出参数数组，MATLAB程序从输出参数数组中得到数据处理的结果。表5.1mexFunction函数的输入参数参数含义intnlhs输出参数个数mxArray*plhs输出参数的mxArray数组intnrhs输入参数个数mxArray*prhs输入参数的mxArray数组根据以上规则编写完C/C++程序后，还需要使用MATLAB的编译器对程序进行编译，将C/C++程序编译成MEX文件。使用MATLAB编译器之前需要使用mex–setup指令对MATLAB编译器进行配置，配置完成后就可以使用mcc命令对C/C++程序进行编译，在mcc命令中可以指定编译选项，包括头文件和函数库等[29]。CUDAC是C语言的扩展，因此，可以使用上述方法编写mexFunction()函数，并将程序中的C语言部分编译成MEX文件供MATLAB调用。但是，由于程序中有GPU端的代码，这些代码无法使用CPU端的编译器进行编译，必须使用NVIDIA提供的NVCCGPU编译器进行编译，因此，对于程序中的GPU端代码，需要在编译前指定GPU编译器对它进行编译。编译器的指定必须在编译前通过配置文件进行配置，配置完成后，还需要在GPU编译器的编译选项中对需要使用的头文件和函数库进行指定。完成以上操作后，就可以使用编译器对整个CUDAC程序进行编译，这是CPU端的代码将由普通的C语言编译器编译，而GPU端代码则由指定的GPU编译器编译。将程序编译49 西安电子科技大学硕士学位论文成MEX文件后，把MEX文件拷贝到当前目录下，之后就可以使用MATLAB对它进行调用了。在调用MEX时，可以输入MEX文件的名称（不包括扩展名）直接调用。MEX文件的名称最好保证不能与MATLAB函数库已有的函数同名，如果MEX文件与这些函数同名，则MATLAB在调用该函数时会优先调用MEX文件，而不是MATLAB函数库中原来的函数。5.3基于GPU和MATLAB的成像算法实现和优化5.3.1总体方案设计及数据划分采用MATLAB和CUDAC对5.1节的成像算法进行实现的基本思想是：成像算法的主体部分由MATLAB实现，MATLAB主要执行算法中一些计算量不大的操作和以及一些简单的流程控制，而对于算法中计算量大和循环次数多的部分,MATLAB通过调用MEX文件，将其交付给GPU执行。根据图31所示的算法流程图，可以确定出算法中计算量比较大以及循环次数多的操作有以下几个：最小二乘拟合，脉冲压缩，多普勒中心估计，PGA、FFTSHIFT、改进的RMA算法以及大量的矩阵与参考向量对应点相乘的操作。由于在成像算法开发的过程中，成像的方案可能会经常变动，因此，可以将以上这些操作都编写成MEX文件，使其模块化，方便在成像算法的开发中调用以及灵活配置。为了降低MATLAB通过MEX文件调用GPU的复杂性，根据本方案所使用的GPU的具体情况，本方案中的所有GPU处理子模块将根据实际处理的数据量的大小自动计算线程网格和线程块的划分方案。本方案中，线程块大小固定取为1632，线程网格大小根据数据量由以下公式计算D15DimGrid_xx16(5-30)D31DimGrid_yy32其中，DimGrid_x为线程块x轴的维度，DimGrid_y为线程块y轴的维度，Dx为输入数据x轴的维度，D为输入数据y轴的维度。y由于GPU的显存有限，因此当数据量较大时，需要对数据分块。分块方案分两种情况讨论，第一种情况是GPU处理模块中没有FFT操作，第二种是GPU处理模块中有FFT操作。对于GPU处理模块中没有FFT操作的情况，设输入数据的数据量大小为M。D由于所有的GPU子模块中都要临时存储原始数据，所以必须要开辟与数据存储空间相同大小的存储空间，另外，本方案使用的TeslaK40显卡的显存大小为12GB，由此可得21M2GB，即M6GB，由于GPU处理存储输入数据外还要存储配置数DD50 第五章基于GPU和MATLAB的SAR成像算法实现及优化据、中间变量等一些其它的数据，必须为这些数据留出存储空间，因此，综合考虑，本方案取M5GB为是否分块的临界值，当数据量小于5GB时，由于GPU能一次mid处理完所有数据，所以不对数据分块；当数据量大于5GB时，必须对数据分块，设子块大小为M，且由于要多次执行内存复制操作，所以还需要开启3个CUDA流，s这样61M2GB，即M2GB，本方案取M1GB。sss对于GPU处理模块中没有FFT操作的情况，这里仍然设输入的数据量大小为M，子块大小为M。由于GPU执行FFT操作需要额外的存储空间，所以有Ds31M2GB，即M4GB，考虑到中间变量等其他存储空间的需要，这里取DDM3GB作为是否分块的临界点。当数据量小于3GB时，GPU能一次处理完所有mid数据，所以，不对数据进行分块；当数据量大于3GB时，GPU需要对数据进行分块，并开启3个CUDA流，此时91M2GB，即M1.5GB，本方案取M512MB。sss在上述GPU处理子模块中，最小二乘拟合，脉冲压缩，多普勒中心估计，PGA、FFTSHIFT和矩阵与参考向量对应点相乘这些操作都数据没有FFT操作，属于第一种情况；而改进的RMA算法的实现属于第二种情况。在使用MATLAB和CUDAC混合编程对算法进行开发时，数据在CPU端和GPU端的交互是一个关键问题。MATLAB在对复数进行存储时，会将复数的实部和虚部分开存储，因此，一个复数矩阵在实际存储时会分成实数部矩阵和虚数部矩阵两个部分。而在GPU中对回波数据进行处理时必须要定义一个复数结构体，以复数的形式进行运算。在复数结构体中，数据的实部和虚部是存储在连续的存储单元的。因此，GPU端在处理数据之前必须对数据的存储结构进行转换，具体的转换方法，要根据数据是否分块来确定。当GPU不需要对输入数据进行分块时，可以直接通过MATLAB提供的操作mxArray的API函数获取输入数据实数部矩阵和虚数部矩阵的首地址指针，并在GPU端开辟显存将数据拷贝到显存，在GPU的线程进行处理时，每个线程接收一个表示实数部实数变量和一个表示虚实部的实数变量，并在线程中定义一个复数变量，将实数部变量的值赋给复数变量的实部，将虚数变量赋给复数变量的虚部，这样就完成了MATLAB中的复数的存储结构到GPU中复数存储结构的转换，而且，避免了数据在主机端的多次拷贝。当GPU需要对输入数据进行分块时，必须开启CUDA流来掩盖内存复制所带来的时延。此时，主机端的数据必须在页锁定主机内存中，而通过MATLAB中的数据都存储在可分页的内存中，因此，在获取了输入数据实数部矩阵和虚数部矩阵的首地址指针后，必须在主机端在开辟与输入数据同等大小的页锁定主机内存，并将输入数据拷贝到页锁定主机内存中，最后再将数据复制到GPU端进行分块。根据以上的分析，在MEX文件中的mexFunction()函数的执行流程如图5.4所示。51 西安电子科技大学硕士学位论文获取数据位置Y数据量大于阀值吗？N分配页锁定主机内存数据复制CPU->GPU复制数据到页锁定主机内存数据处理将数据从页锁定主机内存复制到GPU端数据复制GPU->CPU数据处理复制数据GPU->CPU结束图5.4mexFunction()函数的执行流程5.3.2GPU的存储优化GPU的存储优化是GPU性能调优的主要内容。在本方案中，可以在两个方面对GPU进行存储优化：第一个方面是采用共享内存使GPU对全局内存进行合并访问，第二个方面是采用归约算法提高GPU的执行效率和程序的精度。在第四章我们已经介绍了共享内存在实现矩阵转置时对程序进行优化的过程，本章主要介绍共享内存通过合并访问全局内存对GPU实现FFTSHITF操作的加速。在整个成像的流程中存在大量的FFTSHITF操作，无论是在MATLAB执行的部分还是在GPU执行的部分都有大量的FFTSHITF操作，因此，在本方案中，FFTSHIFT操作不仅会编译成MEX文件，还会变成一个设备端函数，供GPU端的代码调用。FFTSHIFT操作主要是对向量前后的数据进行交换，因此，在执行FFTSHIFT操作的过程中，向量中的各个数据的位置都发生了变化。在用GPU实现FFTSHIFT操作时，需要在全局内存中开辟新的存储空间，用于存放处理后的数据。当执行52 第五章基于GPU和MATLAB的SAR成像算法实现及优化FFTSHIFT操作时，GPU需要从全局内存中读取数据，此时满足全局内存合并访问的条件，因此，可以把一个half-warp中16个线程对全局内存的16次访问合并为一次访问[30]。但是，GPU处理完毕后向全局内存中写入数据时，全局内存合并访问的条件不能得到满足，所以，必须在每一个线程块中开辟与线程块中输入数据同等大小的共享内存，在每个线程块处理完毕后将数据写如共享内存，而不是直接写入全局内存。等到所有线程都处理完后，再将共享内存中的数据写入到全局内存中，此时就可以满足合并访问的条件，使数据的写入性能大幅提高。在PGA、多普勒中心估计和最小二乘拟合时，有大量的求和、求均值操作，这些操作都涉及到浮点数的累加。在GPU中执行累加，可以采用并行归约算法进行优化。并行归约算法的操作过程如图5.所示，当线程块中的所有线程执行完毕需要进行累加时，将需要累加的数据分成等长的前后两段，两段中对应位置的数据在一个线程中执行相加操作，执行完后，需要累加的数据就减少了一半，如果这时数据个数一个大于1就重复之前分段相加的过程，直至需要累加的数据个数为1，此时该数就是所有数据的累加和。第0趟0123……303132333435……6263第1趟01…151617…31...第4趟0123第5趟01第6趟0图5.5并行归约算法示意图使用归约算法能大大地减少计算量，假设需要进行累加的数据个数为N，每次求和需要的时间为t。如果只用一个线程采用连续累加的方式进行求和，所需要的总s时间为Nt1；如果采用归约算法进行求和，需要的总时间为tNlog。当数据量ss2较大时，tNlog远远小于Nt1，因此归约算法使得计所需的时间大大减小。另s2s外，连续累加的方式对浮点数求和容易造成“大数吃小数”的问题，使得计算精度降低，而使用归约算法则能够大大减小这种情况发生的概率。5.3.3实验结果及分析-本节采用一段机载SARP波段实测数据对上述使用MATLAB和CUDAC实现的53 西安电子科技大学硕士学位论文算法的有效性进行验证，并测量程序执行的性能。分别用算法的上述实现方案和算法的MATLAB实现方案对测试数据进行处理，并计算两者的处理时间，得到了程序的加速比。程序的运行环境和雷达参数如下表所示。表5.2程序运行环境和雷达参数CPU酷睿i7-4790（一颗）GPUNVIDIATeslaK40c（一部）主机内存128GB操作系统MicrosoftWindows7（x64）实测数据字节数2.6GBsinc插值点数16平台速度120m/s信号带宽100MHz波长0.5m场景中心距雷达距离22km测试时，在CPU端仍然开启16个线程，并进行3次对比测试。用本文方案对实测数据处理后的成像结果如图5.6所示。图5.6实测数据处理结果MATLAB执行时间和混合编程方案执行时间的对比如表5.3所示。54 第五章基于GPU和MATLAB的SAR成像算法实现及优化表5.3两种方案耗时对比MATLAB执行时间混合编程方案执行时间加速比1097.21s118.15s9.291093.34s119.94s9.121094.51s121.37s9.02从表中可以看出，采用MATLAB和CUDA混合编程的模式，对算法实现了良好的加速，三次测试分别达到了9.29倍、9.12倍和9.02倍的加速比。5.4本章总结本章介绍了一种适用于机载宽幅SAR实时成像处理的算法设计方案，将回波数据划分成多个子孔径，在子孔径中对回波数据进行运动补偿和成像处理，借鉴FBP算法的思想，通过相干叠加的方式将子孔径图拼接成高分辨率的全孔径图。采用GPU和MATLAB混合编程的形式对算法设计方案进行了实现和优化，程序主体部分采用MATLAB编写，计算密集的部分使用MATLAB和CUDAC混合编程的方法调用GPU进行实现，并对程序进行了数据并行划分优化和存储优化。使用实测数据对程序性能进行的测试，测试结果表明GPU对程序实现了良好的加速。55 西安电子科技大学硕士学位论文56 第六章总结与展望第六章总结与展望6.1论文总结本文主要研究了一种适用于机载宽波束SAR成像的算法的GPU实现和优化。首先，本文阐述了SAR成像的基本原理，介绍了GPU架构的发展历程和CUDA编程模型，针对本文中设计方案所使用的TeslaK40型号的显卡，详细分析了TeslaK40硬件结构和存储器层次，为之后的成像算法的GPU实现和优化奠定基础。然后，针对机载宽波束SAR成像中合成孔径时间长、数据量大，传统RMA算法对回波二维波数谱利用率不高的问题，采用了一种改进的RMA算法，并使用GPU对算法进行了实现和优化。考虑到实际数据处理过程中数据量大的情况，在使用GPU实现算法时对数据进行了分块并开启了多个CUDA，并提出了一种在数据量超过GPU显存的情况下，在GPU端完成矩阵转置的新方法。使用GPU中的共享内存对程序进行了存储优化，并使用仿真的数据对程序的性能进行了测试，测试结果表明采用GPU对算法进行实现能达到14倍左右的加速比。最后，针对一种处理机载宽幅SAR实测数据的算法设计方案，采用MATLAB和CUDAC混合编程的方法对方案进行了实现和加速优化。将算法中比较耗时的部分使用GPU对其进行实现和加速，并封装成MEX文件的形式供MATLAB直接调用，这样，既对程序进行了加速，又方便了在算法开发过程中对各个模块灵活调用。使用一组实测数据，对程序的性能进行了测试，测试结果表明各个子模块都达到了9倍左右的加速比，整个程序也得到了良好的加速。6.2工作展望根据以上对本文工作的总结，对未来工作的展望主要有以下几个方面：（1）针对机载宽幅SAR实测数据处理中数据量和计算量大的问题，采用多GPU的设计方案对算法进行实现和加速，提高算法执行的实时性。（2）将SAR成像处理领域的各个常用子模块，采用GPU编程对其进行优化加速，并编写成MEX文件，形成一个SAR程序算法开发的工具包，供MATLAB调用，提高SAR成像算法开发的效率。57 西安电子科技大学硕士学位论文58 参考文献参考文献[1]朱柳.高分辨率SAR场景分布的重建方法研究[D].成都:电子科技大学,2013.[2]保铮,邢孟道,王彤.雷达成像技术[M].北京:电子工业出版社,2005.[3]LanG.Cumming,FrankH.Wong.合成孔径雷达——算法与实现[M].北京:电子工业出版社,2012.[4]龚伟.机载SAR成像处理及软件实现[D].北京:中国科学院电子学研究所,2006.[5]仇晓兰,丁赤飚,胡东辉.双站SAR成像处理技术[M].北京:科学出版社,2010.[6]李俊艳.超级计算与物理学研究[J].中国科技纵横,2014(7):243-244.[7]孟大地.基于NVIDIAGPU的机载SAR实时成像处理算法CUDA设计与实现[J].雷达学报,2013(12):482-483.[8]陆珺.近程脉冲雷达关键技术研究[D].南京:南京理工大学,2008.[9]刘齐发.基于现代数字信号处理理论提高SAR图像分辨率[D].武汉:武汉理工大学,2008.[10]王平.基于距离-多普勒算法的调频连续波SAR成像研究[D].南京:南京理工大学,2008.[11]李磊.合成孔径雷达成像及自聚焦算法研究[D].南京:南京理工大学,2009.[12]皮亦鸣.合成孔径雷达成像原理[M].成都:电子科技大学出版社,2007.[13]张娜.合成孔径雷达成像及成像质量的改善[D].西安:西安电子科技大学,2010.[14]刘寒艳,宋红军,程增菊.条带模式、聚束模式和滑动聚束模式的比较[J].中国科学院研究生院学报,2011(3):410-417.[15]汪连栋,杜静.成像雷达并行仿真优化技术[M].北京:电子工业出版社,2012.[16]王磊.基于多核架构的生物序列比对算法的设计与实现[D].重庆:重庆邮电大学,2011.[17]TimothyG.Mattson,BeverlyA.Sanders,BernaL.Massingill.并行编程模式[M].北京:机械工业出版社,2015.[18]李松涛.并行多处理器系统容错的研究与实现[D].成都:电子科技大学,2006.[19]王登平.基于PC集群的并行FP-Growth算法的研究与实现[D].西安:西安电子科技大学,2011.[20]孙成.基于GPU加速的电力系统并行潮流的研究[D].哈尔滨:哈尔滨工业大学,2012.[21]侯毅,沈彦男,王睿索,余宗超.基于GPU的数字影像的正射纠正技术的研究[J].现代测绘,2009(3):10-11.[22]商凯,胡艳.基于GPU的AES算法实现[J].电子技术,2011(5):9-11.[23]胡松,基于CUDA的三维非局部均值滤波并行算法设计[D].武汉:华中科技大学,2012.[24]ShaneCook.CUDAProgramming[M].北京:机械工业出版社,2014.[25]JasonSanders,EdwardKandrot.GPU高性能编程CUDA实战[M].北京:机械工业出版社,59 西安电子科技大学硕士学位论文2010.[26]张彬,李凉海.机载SAR运动补偿技术研究[J].遥测遥控,2007(1):130-134.[27]刘维.精通MATLAB和C/C++混合程序设计[M].北京:北京航空航天大学出版社,2008.[28]潘大夫,汪渤,周志强.MATLAB和C/C++混合编程技术研究[J].计算机工程与设计,2009(2):465-468.[29]陈聪,蒋丽.基于MATLAB和C/C++的仿真程序设计[J].科技风,2013(16):11-11.[30]刘来国.基于GPU的LARED-P加速技术的研究与实现[D].长沙:国防科技大学,2009.60 西安电子科技大学硕士学位论文致谢光阴似箭，日月如梭，不经意间，两年半的研究生生涯即将结束。回首两年半的研究生学习生涯，在各位师长、各位同学和各位朋友的帮助下，我一步步地成长，逐渐对SAR成像算法和GPU编程技术有了一定的研究。执笔落成之际，各位师长的谆谆教诲，各位同学和朋友的无私帮助仍然历历在目，在此，向他们致以最诚挚的感谢！首先，我要感谢我的导师邢孟道教授。从进入实验室开始，邢老师在科研上就给予了我很多的指导和帮助，我的每一点成绩的取得都离不开邢老师辛勤的培养。邢老师渊博的学识，对科学研究孜孜不倦的追求，一丝不苟的作风和严谨求实的治学态度使我受益匪浅，也激励着我在科研的路上不断前进。本论文也是在邢老师的悉心指导下完成的，从论文的选题到论文的撰写，无不渗透着刑老师的心血，值此论文完稿之际，谨对邢老师的辛勤培育以及谆谆教诲表示最衷心的感谢!这里，还要感谢孙光才老师一直以来对我的关心和指导。在我研究生生涯的各个阶段，孙老师都给予了我很大的帮助和支持。孙光才老师渊博的学识，对科学的不懈追求和忘我的工作精神，将是我在科研道路上的榜样。感谢杨军师兄一直以来对我的帮助和鼓励，杨军师兄不仅在科研上给予了我很大的帮助，在生活中，杨军师兄也教会了我很多道理。感谢实验室的张双喜、杨桃丽、吴玉峰、张佳佳、周芳、杨泽民、李学士、邵鹏、段佳、曾乐天、李震宇、孟自强、段佳、吴敏、左邵山、方午梅等各位师兄师姐对我的帮助和鼓励，在与各位师兄师姐的交流和讨论中让我获益良多。感谢景国彬、陈溅来、董琪、潘梦冠一直以来对我学习和生活中的帮助和支持，感谢张冲、胡贵彬、汪路峰、陆小鱼、劳丹涤、施佳刚、刘士杰、杜凡、刘瑞红、马倩等同学对我的帮助，还要感谢同宿舍的张云骥、肖恒、荆善灏对我生活和学习上的帮助，和你们度过的两年半的美好时光我将一直珍藏。感谢所有关心过、帮助过和支持过我的人！。61 西安电子科技大学硕士学位论文62 作者简介作者简介1.基本情况张锐，男，湖北蕲春人，1990年3月出生，西安电子科技大学电子工程学院信号与信息处理专业2013级硕士研究生。2.教育背景2009.08～20013.07西安电子科技大学，本科，专业：电子信息工程2013.08～西安电子科技大学，硕士研究生，专业：信号与信息处理3.攻读硕士学位期间的研究成果3.1参与科研项目及获奖[1]504所纵向项目,GEOSAR成像处理算法研究,2013.9-2014.9,项目通过验收,根据轨道特性参数建立逼真星地几何模型和精确的卫星和场景的运动模型。针对同步轨道SAR的特点，对同步轨道SAR进行偏航牵引分析。用MATLAB对以上模型进行仿真验证，并编写GUI界面。[2]xx所横向项目，xx对地探测技术，2014.5-2014.12，项目通过验收，针对算法中比较耗时的部分，采用GPU对算法进行加速。63 ．Ｉ爲參柄術專ＸｍｉＡＮＵＮＩＶＥＲＳＩＴＹ地址：西香市太白南路２号ｄｕ网址．ｘｉｄｉａｎ．ｅ．ｃｎ：ｗｗｗ－：Ｉ：！
此文档下载收益归作者所有
