sar合成孔径雷达图像点目标仿真报告附matlab代码

1、SAR图像点目标仿真报告徐一凡1SAR原理简介合成孔径雷达(SyntheticApertureRadar，简称SAR)是一种高分辨率成像雷达技术。它利用脉冲压缩技术获得高的距离向分辨率，利用合成孔径原理获得高的方位向分辨率，从而获得大面积高分辨率雷达图像。SAR回波信号经距离向脉冲压缩后，雷达的距离分辨率由雷达发射信号带宽决定：CP二，式中P表示雷达的距离分辨率，B表示雷达发射信号带宽，C表示光速。同r2Brrr样，SAR回波信号经方位向合成孔径后，雷达的方位分辨率由雷达方位向的多谱勒带宽决定：P二?，式中P表示雷达的方位分辨率，B表示雷达方位向多谱勒带宽，v表示aBaaaa方位向SAR平台速

2、度。在小斜视角的情况下，方位分辨率近似表示为p二D，其中D为a2方位向合成孔径的长度。2SAR的几何关系雷达位置和波束在地面覆盖区域的简单几何模型如图1所示。此次仿真考虑的是正侧视的条带式仿真，也就是说倾斜角为零，SAR波束中心和SAR平台运动方向垂直的情况。图1雷达数据获取的几何关系建立坐标系XYZ如图2所示，其中XOY平面为地平面；SAR平台距地平面高H,以速度V沿X轴正向匀速飞行；P点为SAR平台的位置矢量，设其坐标为(x,y,z);T点为目标的位置矢量，设其坐标为(兀丁,打,zj；由几何关系，目标与SAR平台的斜距为：R二PT=J(x-x)2+(y-y)2+(z-z)2ITtt(1)由

3、图可知：y=°，z=H,Zf=0;令x=vs-，其中v为平台速度,s为慢时间变量(slowtime),假设需=vs，其中s表示SAR平台的x坐标为需的时刻；再令厂=込2+晋，r表示目标与SAR的垂直斜距，重写(1)式为：PT=R(s;r)=r2+v2-(s-s)2(2)R(s;r)就表示任意时刻s时，目标与雷达的斜距。一般情况下，v|s-s0|<<r，于是通过傅里叶技术展开，可将(2)式可近似写为：v2R(s;r)=r2+v2-(s一s)2qr+(s一s)2(3)丫02r0可见，斜距是sWr的函数，不同的目标，r也不一样，但当目标距SAR较远时，在观测带内，可近似认为r不

4、变，即r=R0。图2:空间几何关系(a)正视图(b)侧视图图2(a)中，Lsa厂表示合成孔径长度，它和合成孔径时间Tsa厂的关系是Lsar=vTsar.(b)中，&为雷达天线半功率点波束角，°为波束轴线与Z轴的夹角，即波束视角，Rmin为近距点距离，Rmax为远距点距离，w为测绘带宽度，它们的关系为：Rmin二H-tg(0-Rmax二H-tg(0+辱；)(4)W二Rmax-Rmin3SAR的回波信号模型SAR在运动中以一定的周期(1/PR巧发射和接收信号，具体过程如图3所示。发射机以工/的时间发射啁啾脉冲，然后切换天线开关接收回波信号。图3雷达发射脉冲串的时序当雷达不处于发射

5、状态时,它接收3反射回波。发射和接收回波的时间序列如图4所示。在机载情况下，每个回波可以在脉冲发射间隔内直接接收到。但是在星载情况下，由于距离过大，某个脉冲的回波要经过610个脉冲间隔才能接收到。这里仿真为了方便，默认为机载情况。图4脉冲雷达的发射与接收周期假设T为chirp信号持续时间，下标r表示距离向；PRF为重复频率，PRT为重复周期，r等于1/PRF。接收序列中，工二2*号s;r)表示发射第i个脉冲时，目标回波相对于发射nC序列的延时。雷达的发射序列数学表达式为式(5):s(t)=遠p(t-n*PRT)n=s(5)P(t)=rect()ej冗Krt2ej2kftr式中，WCtH表示矩形

6、信号，Kr为距离向的chirp信号调频率，f为载频。雷达回波信号由发射信号波形，天线方向图，斜距，目标RCS，环境等因素共同决定，若不考虑环境因素，则单点目标雷达回波信号可写成式(6)所示：s(t)=£Qwp(t-n-PRT-T)(6)rnn其中，Q表示点目标的雷达散射截面，w表示点目标天线方向图双向幅度加权，T表n示载机发射第n个脉冲时，电磁波再次回到载机时的延时T二2*R('；厂),带入式中得:nC8/nPRT2R(s;r)/Cs(t)二w-rect()-rTn=8rexpj兀K(tn-PRT2R(s;r)/C)2-r4兀exp-jR(s;r)-expj2兀f(tn-PR

7、TT)人cn式就是单点目标回波信号模型，其中，expjK(tn-PRT2R(s;r)/C)2是r4兀chirp分量，它决定距离向分辨率；exp-jR(s;r)为多普勒分量，它决定方位向分辨率。入对于任意一个脉冲，回波信号可表示为式(8)所示：s(t,s)=Aw(t2R(s;r)/C)w(ss)xexpj4兀fR(s;r)/C-r0rac'.0(8)expj兀K(t2R(s;r)/C)2r我们知道，由于R(s;r)随慢时间s的变化而变化，所以计算机记录到的回波数据存储形式如图5所示：图5目标照射时间内，单个点目标回波能量在信号处理器的二维存储器中的轨迹4距离徙动及校正根据图2可知，在倾斜

8、角为零或很小的时候，目标与雷达的瞬时距离为R(s;r)，根据几何关系可知，R(s;丫)=丫2+v2-(ss0)2，根据泰勒级数展开可得:R(s;r)=Jr2+v2-(ss)2沁r+(ss)2(9)v02r0q由式(9)可知，不同慢时间对应着不同的R(s;厂)，并且是一个双曲线形式或者近似为一个二次形式。如图5所示,同一目标的回波存储在计算机里不在同一直线上，存在距离徙动。从而定义距离徙动量：v2AR(s,r)=(ss)2(10)2r0为了进行方位向的压缩，方位向的回波数据必须在同一条直线上，也就是说必须校正距离徙动AR(s,r)。由式(10)可知，不同的最近距离r对应着不同的AR(s,r)，因

9、此在时域处理距离徙动会非常麻烦。因此，对方位向进行傅里叶变换，对距离向不进行变换，得到新的域。由于方位向的频率即为多普勒频率，所以这个新的域也称为距离多普勒域。将斜距R写成多普勒fa的函数，即R(f,r)。众所周知，对最近距离为r的点目标P,a2V回波多普勒f是倾斜角0的函数，即/=vsin0，斜距R(f,r)=r/cos9，于是aa九aR(f,r)=r/cos9=r/、；1一sin292(11)a1(九)£r+8(V)2fa1a所以距离多普勒域中的我距离徙动为AR(f,r)=-(-)2rf2，可发现它不随慢时间变换，a8Va同一最短距离r对应着相同大小的距离徙动。因此在距离多普勒域

10、对一个距离徙动校正就是对一组具有相同最短距离的点目标的距离徙动校正，这样可以节省运算量。为了对距离徙动进行校正，需要得到距离徙动单元，即距离徙动体现在存储单元中的移动数值，距离徙动单元可以表示为AR(f,r)/p，这个值通常为一个分数，由于存储单元ar都是离散的，所以不同通过在存储单元简单的移动得到准确的值。为了得到准确的徙动校正值，通常需要进行插值运算。本仿真采用了两种插值方法最近邻点插值和sinc插值，下面分别进行介绍。最近邻点插值法的优点是简单而快速，缺点是不够精确。AR(f,r)/p=N+n，其中N为整数部分，arn为小数部分，整数部分徙动可以直接通过平移消除，对于小数部分则通过四舍五

11、入的方法变为0或者1,这样就可以得到较为精确的插值。Sinc插值原理如下：在基带信号下，卷积核是sinc函数sin(兀x)h(x)=sinc(x)=(12)兀x，、_-'插值信号为g(x)=工g(i)sinc(x-i)(13)di即为所有输入样本的加权平均。可通过频域来理解，如图6所示，采样信号g(i)的频谱G(f)等于以采样率重复的信dd号频谱。为了重建信号g(x)，只需要一个周期频谱(如基带周期)，因此需要理想矩形低通滤波器在频域中提取基带频谱(如图6)所示。已知该理想滤波器在时域中是sinc函数。由于频域相乘相当于时域卷积，故插值可以通过与sinc核的卷积来实现。图6理想低通滤波

12、器怎样对采样信号进行插值5点目标成像matlab仿真5.1距离多普勒算法距离多普勒算法(RDA)是在1976年至1978年为民用星载SAR提出的，它兼顾了成熟、简单、高效和精确等因素，至今仍是使用最广泛的成像算法。它通过距离和方位上的频域操作，到达了高效的模块化处理要求，同时又具有了一维操作的简便性。图7示意了RDA的处理流程。这里主要讨论小倾斜角及短孔径下的基本RDA处理框图。1.当数据处在方位时域时，可通过快速卷积进行距离压缩。也就是说，距离FFT后随即进行距离向匹配滤波，再利用距离IFFT完成距离压缩。回波信号为：14)s(t,s)二Awt-2R(s)/cw(ss)00racxexp-j

13、4kfR(s)/cexpj兀K(t-2R(s)/c)20r距离向压缩后的信号为：s(t,s)=IFFTS(f,s)H(f)t0tt二Apt2R(s)/cw(ss)expj4kfR(s)/c0rac0rc15)ff2h(打=rectfKTTexp-加丘exp-j加如16)2通过方位FFT将数据变换至距离多普勒域，多普勒中心频率估计以及大部分后续操作都在该域进行。方位向傅里叶变换后信号为：S(t,f)=FFTs(t,s)1ssrc2R(f)=Apt“W(ff)0rca4kfRf2Xexp-jiexpj兀rcKa3.在距离多普勒域进行随距离时间及方位频率变化的RCMC，该域中同一距离上的一组目标轨迹

14、相互重合。RCMC将距离徙动曲线拉直到与方位频率轴平行的方向。这里可以采用最近邻点插值法或者sine插值法，具体插值方法见前面。假设RCMC插值是精确的，信号变为：(17)2RS(t,f)=Ap(t)W(ff)20rcassc4kfR】rxexp-j00expjkc(18)4通过每一距离门上的频域匹配滤波实现方位压缩。匚Ka为进行方位压缩，将RCMC后的S(t,f)乘以频域匹配滤波器H(f)。azsf2H(f)二expjfrazsKa21)(19)S(t,f)=S(t,f)H(f)3s2sazsssc=Ap(t2R/c)W(ff)expjf0Ro)0r0asscc5最后通过方位IFFT将数据变

15、换回时域，得到压缩后的复图像。复原后的图像为:s(t,s)=IFFTS(t,f)acs3s二Ap(t-2R/c)p(s)0r0a4kfRxexp-jooexpj2兀fscsc图8距离多普勒算法流程图5.2 ChirpScaling算法距离多普勒算法具有诸多优点，但是距离多普勒算法有两点不足：首先，当用较长的核函数提高距离徙动校正（RCMC）精度时，运算量较大；其次，二次距离压缩（SRC）对方位频率的依赖性问题较难解决，从而限制了其对某些大斜视角和长孔径SAR的处理精度。ChirpScaling算法避免了RCMC中的插值操作，通过对Chirp信号进行频率调制，实现了对该信号的尺度变换或平移。图8

16、显示了ChirpScaling算法处理流程。这里主要讨论小倾斜角及短孔径下的基本CSA处理框图。主要步骤包括四次FFT和三次相位相乘。1通过方位向FFT将数据变换到距离多普勒域。2. 通过相位相乘实现ChirpScaling操作，使所有目标的距离徙动轨迹一致化。这是第一步相位相乘。用以改变线调频率尺度的ChirpScaling二次相位函数为：H（t,f;R）=exp加丫（f;R）a（f）（t-2Rf；R）2（22）1asaBac3. 通过距离向FFT将数据变到二维频域。4. 通过与参考函数进行相位相乘，同时完成距离压缩、SRC和一致RCMC。这是第二步相位相乘。用于距离压缩，距离徙动校正的相位

17、函数写为：H(f,f;R)=expjn2ras1丫(f；r)1+af)aBaf2rXexpj4兀Rsa(f)fcr(23)5. 通过距离向IFFT将数据变回到距离多普勒域。6通过与随距离变化的匹配滤波器进行相位相乘，实现方位压缩。此外，由于步骤2中的ChirpScaling操作，相位相乘中还需要附加一项相位校正。这是第三步相位相乘。补偿由ChirpScaling引起的剩余相位函数是：H（t,f;R）二expj2Rf2-f2expj0（f;R）（24）2raBVBaMaAaB7.最后通过方位向IFFT将数据变回到二维时域，即SAR图像域。图8ChirpScaling算法流程图简而言之，R-D算法

18、是将徙动曲线逐一校正，CS算法是以某一徙动曲线为参考，在Doppler域内消除不同距离门的徙动曲线的差异，令这些曲线成为一组相互，平行”的曲线，然后在二维频率域内统一的去掉距离徙动。通俗一点就是，RD算法是将弯曲的信号一根根掰直，而CS算法是先把所有信号都掰得一样弯，然后再统一掰直。6仿真结果6.1使用最近邻点插值的距离多普勒算法仿真结果本文首先对5个点目标的回波信号进行了仿真，5个点目标构成了矩形的4个顶点和中心，其坐标分别如下，格式为（方位向,距离向,后向反射系数）：09750110097501501000010102501100102501图9的上图是距离向压缩后的图像，从图中可以看到5

19、条回波信号（其中有几条部分重合，但仍能看出来）目标回波信号存在明显的距离徙动，需要进行校正。图9的下图是通过最近邻点插值法校正后的图像，可以看出图像基本被校正为直线。图9距离向压缩后最近邻点插值的结果图10为进行方位向压缩后形成的图像，可以明显看出5个点目标，并且5个点目标构成了矩形的四个顶点及其中心。图10通过最近邻点插值生成的点目标图像6.2 使用最近邻点插值的距离多普勒算法仿真结果图11上图为通过距离压缩后的图像，图11的下图为通过sine插值法校正后的图像。图11距离向压缩后sine插值的结果图12为进行方位向压缩后形成的图像，可以明显看出5个点目标，并且5个点目标构成了矩形的四个顶点

20、及其中心。图12通过sine插值生成的点目标图像6.3 ChirpScaling算法仿真结果同样，在ChirpSealing中，对5个点目标的回波信号进行了仿真，5个点目标构成了矩形的4个顶点和中心，其坐标分别如下，格式为（方位向,距离向，后向反射系数）：1200011250-50112505011150-5011150501图13是仿真的雷达回波信号图图13仿真出来的SAR回波信号图14是经过第一次相位校正之后，通过距离向压缩后的距离时域-方位时域信号图（ChirpSealing算法的七个步骤中并不包含该信号，该信号是将步骤2之后的信号通过方位向傅里叶逆变换，再进行距离向压缩得到的，只为了验

21、证原理）。按照理论，该图中所有点的距离徙动都应该一样。从图中大致可看出，五个点的距离徙动是差不多的。图14ChirpSealing第2步之后、经过距离向压缩得到的图图15为步骤5之后，信号距离压缩，距离徙动校正之后的距离多普勒域中的信号图。图15距离徙动校正之后的图图16为步骤6之后，消除相位偏移的图。图16消除相位偏移的图图17为通过ChirpSealing算法生成的点目标图像图17通过ChirpSealing算法生成的点目标图像6.4几种算法比较本文讨论了距离多普勒算法和ChirpSealing算法，其中距离多普勒算法考虑了最近邻点插值和sine插值两种插值方法。距离多普勒算法兼顾了成熟、

22、简单、高效和精确等因素，至今仍被广泛使用，但是距离多普勒算法有两点不足：首先，当用较长的核函数提高距离徙动校正RCMC）精度时，运算量较大；其次，二次距离压缩（SRC）对方位频率的依赖性问题较难解决，从而限制了其对某些大斜视角和长孔径SAR的处理精度。最近邻点插值的优点是速度快，该插值的运行时间为2.267137秒，缺点是不够精确；sine插值的优点是精确，该方法的运行时间为29.148728秒，缺点是速度慢；ChirpSealing算法避免了插值运算，提高了速度，运行时间为0.323327秒，但是其算法较为复杂。%=%文件名：NearSAR.m%作者：徐一凡%功能：合成孔径雷达距离多普勒算法

23、点目标成像%=clear;clc;closeall;%=%常数定义C=3e8;%雷达参数Fc=1e9;lambda=C/Fc;%目标区域参数Xmin=0;Xmax=50;Yc=10000;Y0=500;%光速%载频1GHz%波长%目标区域方位向范围Xmin,Xmax%成像区域中线%目标区域距离向范围Yc-YO,Yc+YO%成像宽度为2*Y0%轨道参数V=100;H=5000;R0=sqrt(Yc"2+2);%天线参数D=4;%SAR的运动速度100m/s%高度5000m%最短距离%方位向天线长度Lsar=lambda*RO/D;%SAR合成孔径长度，合成孔径雷达成像算法与实现P.1O

24、OTsar=Lsar/V;%慢时间域参数Ka=-2*V"2/lambda/R0;Ba二abs(Ka*Tsar);PRF二Ba;又为复频率，所以等于PRT=1/PRF;ds=PRT;%SAR照射时间%多普勒频域调频率P.93%多普勒频率调制带宽%脉冲重复频率，PRF其实为多普勒频率的采样率,%脉冲重复时间%慢时域的时间步长Nslow=ceil(Xmax-Xmin+Lsar)/V/ds);%慢时域的采样数，ceil为取整函数，结合P.76的图理解%nextpow2为最靠近2的幕次函数，这里为Nslow=2nextpow2(Nslow);fft变换做准备sn=linspace(Xmin-L

25、sar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢时间域的时间矩阵PRT=(Xmax-Xmin+Lsar)/V/Nslow;%由于Nslow改变了，所以相应的一些参数也需要更新，周期减小了PRF=1/PRT;ds=PRT;%快时间域参数设置%脉冲持续时间5usTr=5e-6;Br=30e6;%chirp频率调制带宽为30MHzKr二Br/Tr;%chirp调频率Fsr=2*Br;%快时域采样频率，为3倍的带宽dt=1/Fsr;%快时域采样间隔Rmin=sqrt(Yc-Y0厂2+2);Rmax=sqrt(Yc+Y0厂2+2+(Lsar/2厂2);Nfast=ceil(2*(Rm

26、ax-Rmin)/C/dt+Tr/dt);%快时域的采样数量Nfast=2"nextpow2(Nfast);%更新为2的幕次，方便进行fft变换tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);%快时域的离散时间矩阵dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;%更新间隔Fsr=1/dt;%分辨率参数设置%距离向分辨率%方位向分辨率''%点目标的数量%点目标位置，这里设置了5个点目标，构DY=C/2/Br;DX=D/2;%点目标参数设置Ntarget=5;%点目标格式x,y,反射系数sigmaPtarget=Xmin,Y

27、c-50*DY,l成一个矩形以及矩形的中心Xmin+50*DX,Yc-50*DY,1Xmin+25*DX,Yc,1Xmin,Yc+50*DY,1Xmin+50*DX,Yc+50*DY,1;disp('Parameters:')%参数显示%=%生成回波信号K=Ntarget;N=Nslow;M=Nfast;T=Ptarget;Srnm=zeros(N,M);fork=1:1:Ksigma=T(k,3);Dslow=sn*V-T(k,1);disp('SamplingRateinfast-timedomain');disp(Fsr/Br)disp('Samp

28、lingNumberinfast-timedomain');disp(Nfast)disp('SamplingRateinslow-timedomain');disp(PRF/Ba)disp('SamplingNumberinslow-timedomain');disp(Nslow)disp('RangeResolution');disp(DY)disp('Cross-rangeResolution');disp(DX)disp('SARintegrationlength');disp(Lsar)disp(

29、'Positionoftargets');disp(Ptarget)%目标数目%慢时域的采样数%快时域的采样数%目标矩阵%生成零矩阵存储回波信号%总共K个目标%得到目标的反射系数%方位向距离，投影到方位向的距离%实际距离矩阵R=sqrt(Dslow.2+T(k,2)2+H2);tau=2*R/C;%回波相对于发射波的延时Dfast=ones(N,l)*tm-tau'*ones(l,M);%(t-tau),其实就是时间矩阵，ones(N,l)和ones(1,M)都是为了将其扩展为矩阵phase二pi*Kr*Dfast."2-(4*pi/lambda)*(R

30、9;*ones(1,M);%相位，公式参见P.96Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*(abs(Dslow)<Lsar/2)'*ones(1,M);%由于是多个目标反射的回波，所以此处进行叠加end%=%距离-多普勒算法开始%距离向压缩tic;tr=tm-2*Rmin/C;Refr=exp(j*pi*Kr*tr."2).*(0tr&trTr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr);Gr二abs(Sr);%开始进行距离弯曲补偿正侧

31、视没有距离走动项主要是因为斜距的变化引起回波包络的徙动%补偿方法：最近邻域插值法，具体为：先变换到距离多普勒域，分别对单个像素点计算出距离徙动量，得到距离徙动量与距离分辨率的比值，%该比值可能为小数，按照四舍五入的方法近似为整数，而后在该像素点上减去徙动量%方位向做fft处理再在频域做距离弯曲补偿Sa_RD=ftx(Sr);%方位向FFT变为距离多普域进行距离弯曲校正%距离徙动运算，由于是正侧视，fdc=0,只需要进行距离弯曲补偿。Kp=1;%计算或者预设预滤波比h=waitbar(0,'最近邻域插值');%首先计算距离迁移量矩阵forn=1:N%总共有N个方位采样form=1

32、:M%每个方位采样上有M个距离采样delta_R=(1/8)*(lambda/V厂2*(R0+(m-M/2)*C/2/Fsr)*(n-N/2)*PRF/N厂2;%距离迁移量P.160；(R0+(m-M/2)*C/2/Fsr):每个距离向点m的R0更新；(n-N/2)*PRF/N：不同方位向的多普勒频率不一样RMC=2*delta_R*Fsr/C;%此处为delta_R/DY，距离徒动了几个距离单元delta_RMC=RMC-round(RMC);%距离徒动量的小数部分ifm+round(RMC)M%判断是否超出边界Sa_RD(n,m)=Sa_RD(n,M/2);elseifdelta_RMC=

33、0.5%五入Sa_RD(n,m)=Sa_RD(n,m+round(RMC)+1);else%四舍Sa_RD(n,m)=Sa_RD(n,m+round(RMC);endendendwaitbar(n/N)endclose(h)%=Sr_rmc=iftx(Sa_RD);%距离徙动校正后还原到时域Ga=abs(Sr_rmc);%方位向压缩ta=sn-Xmin/V;Refa二exp(j*pi*Ka*ta."2).*(abs(ta)Tsar/2);Sa=iftx(ftx(Sr_rmc).*(conj(ftx(Refa).'*ones(1,M);Gar=abs(Sa);toc;%=%绘图

34、colormap(gray);figure(l)subplot(211);row二tm*C/2-2008;col二sn*V-26;imagesc(row,col,255-Gr);%距离向压缩，未校正距离徙动的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel('距离向'),ylabel('方位向'),title('距离向压缩，未校正距离徙动的图像')，subplot(212);imagesc(row,col,255-Ga);%距离向压缩，校正距离徙动后的图像axis(Yc-Y0,Yc+Y0,Xmi

35、n-Lsar/2,Xmax+Lsar/2);xlabel('距离向'),ylabel('方位向'),title('距离向压缩，校正距离徙动后的图像')，figure(2)colormap(gray);imagesc(row,col,255-Gar);%方位向压缩后的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel('距离向'),ylabel('方位向'),title('方位向压缩后的图像'),%=%文件名:SincSAR.m%作者：徐一凡%功能:合

36、成孔径雷达距离多普勒算法点目标成像%=clear;clc;closeall;%=%常数定义C=3e8;%光速%雷达参数Fc=1e9;%载频1GHzlambda=C/Fc;%波长%目标区域参数Xmin=O;Xmax=5O;%目标区域方位向范围Xmin,XmaxYc=1OOOO;YO=5OO;%成像区域中线%目标区域距离向范围Yc-Y0,Yc+Y0%成像宽度为2*Y0%轨道参数V=1OO;H=5OOO;R0=sqrt(Yc"2+2);%天线参数D=4;%SAR的运动速度100m/s%高度5OOOm%最短距离%方位向天线长度Lsar=lambda*RO/D;%SAR合成孔径长度，合成孔径雷

37、达成像算法与实现P.1OOTsar=Lsar/V;%慢时间域参数Ka=-2*V"2/lambda/R0;Ba=abs(Ka*Tsar);PRF=Ba;又为复频率，所以等于PRT=1/PRF;ds=PRT;%SAR照射时间%多普勒频域调频率P.93%多普勒频率调制带宽%脉冲重复频率，PRF其实为多普勒频率的米样率，%脉冲重复时间%慢时域的时间步长Nslow二ceil(Xmax-Xmin+Lsar)/V/ds);%慢时域的采样数，ceil为取整函数，结合P.76的图理解Nslow=2nextpow2(Nslow);%nextpow2为最靠近2的幕次函数，这里为fft变换做准备sn=lin

38、space(Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢时间域的时间矩阵PRT=(Xmax-Xmin+Lsar)/V/Nslow;%由于Nslow改变了，所以相应的一些参数也需要更新，周期减小了PRF=1/PRT;ds=PRT;%快时间域参数设置Tr=5e-6;Br=30e6;%脉冲持续时间5us%chirp频率调制带宽为30MHzKr=Br/Tr;Fsr=2*Br;%chirp调频率%快时域采样频率，为3倍的带宽dt=1/Fsr;%快时域采样间隔Rmin=sqrt(Yc-YO)2+H2);Rmax=sqrt(Yc+Y0厂2+2+(Lsar/2厂2);Nfa

39、st=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%快时域的采样数量Nfast=2"nextpow2(Nfast);%更新为2的幕次，方便进行fft变换tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);%快时域的离散时间矩阵dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;%更新间隔Fsr=1/dt;%分辨率参数设置DY=C/2/Br;DX=D/2;%点目标参数设置Ntarget=5;%点目标格式x,y,反射系数sigmaPtarget=Xmin,Yc-50*DY,1成一个矩形以及矩形的中心%距离向分辨率%方位向分辨率%点

40、目标的数量%点目标位置，这里设置了5个点目标，构Xmin+50*DX,Yc-50*DY,1Xmin+25*DX,Yc,1Xmin,Yc+50*DY,1Xmin+50*DX,Yc+50*DY,l;disp('Parameters:')%参数显示disp('SamplingRateinfast-timedomain');disp(Fsr/Br)disp('SamplingNumberinfast-timedomain');disp(Nfast)disp('SamplingRateinslow-timedomain');disp(PRF

41、/Ba)disp('SamplingNumberinslow-timedomain');disp(Nslow)disp('RangeResolution');disp(DY)disp('Cross-rangeResolution');disp(DX)disp('SARintegrationlength');disp(Lsar)disp('Positionoftargets');disp(Ptarget)%=%生成回波信号K二Ntarget;N=Nslow;M=Nfast;T=Ptarget;Srnm=zeros(N

42、,M);fork=1:1:Ksigma=T(k,3);Dslow=sn*V-T(k,1);R=sqrt(Dslow."2+T(k,2厂2+2);tau=2*R/C;Dfast=ones(N,1)*tm-tau'*ones(1,M);%目标数目%慢时域的采样数%快时域的采样数%目标矩阵%生成零矩阵存储回波信号%总共K个目标%得到目标的反射系数%方位向距离，投影到方位向的距离%实际距离矩阵%回波相对于发射波的延时%(t-tau)，其实就是时间矩阵，ones(N,1)和ones(1,M)都是为了将其扩展为矩阵phase=pi*Kr*Dfast.2-(4*pi/lambda)*(R&

43、#39;*ones(1,M);%相位，公式参见P.96Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*(abs(Dslow)<Lsar/2)'*ones(1,M);%由于是多个目标反射的回波，所以此处进行叠加end%=%距离-多普勒算法开始%距离向压缩tic;tr=tm-2*Rmin/C;Refr二exp(j*pi*Kr*tr."2).*(0tr&trTr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr);Gr=abs(Sr);%开始进行距离弯曲补偿

44、正侧视没有距离走动项主要是因为斜距的变化引起回波包络的徙动%补偿方法：最近邻域插值法，具体为：先变换到距离多普勒域，分别对单个像素点计算出距离徙动量，得到距离徙动量与距离分辨率的比值，%该比值可能为小数，按照四舍五入的方法近似为整数，而后在该像素点上减去徙动量%方位向做fft处理再在频域做距离弯曲补偿Sa_RD=ftx(Sr);%方位向FFT变为距离多普域进行距离弯曲校正%距离徙动运算，由于是正侧视，fdc=O,只需要进行距离弯曲补偿。Kp=1;%计算或者预设预滤波比%第二种方法进行截断sinc插值进行距离徒动校正'h=waitbar(O,'Sinc插值')；P=4;%

45、4点sinc插值RMCmaxtix=zeros(N,M);forn=1:Nform=P:Mdelta_R=(l/8)*(lambda/V厂2*(R0+(m-M/2)*C/2/Fsr)*(n-N/2)*PRF/N厂2;%首先计算距离迁移量计算方法就是把斜距变换到距离多普勒域就知道了RMC=2*delta_R*Fsr/C;%距离徒动了几个距离单元delta_RMC=RMC-round(RMC);%距离徒动量的小数部分fori=-P/2:P/2-1ifm+RMC+i>M%判断是否超出边界RMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,M)*sinc(pi*(-i+R

46、MC);elseRMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,m+round(RMC)+i)*sinc(pi*(-i+delta_RMC);endendendwaitbar(n/N)endclose(h)%=Sr_rmc=iftx(RMCmaxtix);%距离徙动校正后还原到时域Ga=abs(Sr_rmc);%方位向压缩ta=sn-Xmin/V;Refa=exp(j*pi*Ka*ta."2).*(abs(ta)<Tsar/2);Sa=iftx(ftx(Sr_rmc).*(conj(ftx(Refa).'*ones(1,M);Gar=abs(

47、Sa);toc;%=%绘图colormap(gray);figure(1)subplot(211);row=tm*C/2-2008;col=sn*V-26;imagesc(row,col,255-Gr);%距离向压缩，未校正距离徙动的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabelC距离向，),ylabel(，方位向)，title('距离向压缩，未校正距离徙动的图像')，subplot(212);imagesc(row,col,255-Ga);%距离向压缩，校正距离徙动后的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsa

48、r/2,Xmax+Lsar/2);xlabelC距离向'),ylabel('方位向'),title('距离向压缩，校正距离徙动后的图像')，figure(2)colormap(gray);imagesc(row,col,255-Gar);%方位向压缩后的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabelC距离向'),ylabel('方位向'),title('方位向压缩后的图像'),%ChirpScaling算法%徐一凡clearall;clc;%距离向参数range

49、:xdomainTr=200;%时宽200mBr=1;%带宽1Kr=Br/Tr;%调频斜率Fc=4;%载频4Nfast=512;%为了快速运算Xc=1200;X0=150;%定义距离向范围x=Xc+linspace(-X0,X0,Nfast);%x域序列:Xc-X0Xc+X0dx=2*X0/Nfast;%定义步长kx=linspace(-1/dx/2,1/dx/2,Nfast);%kx域序列%方位向参数cross-range:ydomainTa=300;%时宽300m,合成孔径长度Ba=1;%带宽1(1/m)Ka=Fc/Xc;%调频斜率Ka=Ba/Ta=Fc/XcNslow=1024;%为了快

50、速运算Y0=200;y=linspace(-Y0,Y0,Nslow);%y域序列:-Y0Y0dy=2*Y0/Nslow;ky=linspace(-1/dy/2,1/dy/2,Nslow);%ky域序列%目标几何关系targetgeometry%x坐标,y坐标，复后向散射系数Ptar=Xc,0,1+0jXc+50,-50,1+0jXc+50,50,1+0jXc-50,-50,1+0jXc-50,50,1+0j;disp('Positionoftargets');disp(Ptar)%生成SAR正交解调后的回波数据Srnm=zeros(Nfast,Nslow);N=size(Pta

51、r,l);%目标个数h=waitbar(0,'SAR回波生成');fori=1:1:Nxn=Ptar(i,1);yn=Ptar(i,2);sigma=Ptar(i,3);%提取每个目标的信息X=x.'*ones(1,Nslow);%扩充为矩阵Y=ones(Nfast,1)*y;%扩充为矩阵DX=sqrt(xn"2+(Y-yn)."2);%中间变量phase二pi*Kr*(X-DX)."2-2*pi*Fc*DX;%回波相位Srnm二Srnm+sigma*exp(j*phase).*(abs(X-DX)Tr/2).*(abs(Y-yn)Ta/2);%回波累加waitbar(i/N)endclose(h)tic;%数据准备phi0=-x'*sqrt(Fc“2-ky.“2);phi1二-Fc*x'*(1./sqrt(Fc“2-ky.“2);phi2=1/2*x'*(ky2./(F"2-ky2).J.5);Cs=ones(Nfast,1)*(Fc./sqrt(Fc"2-