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

S A R 图像点目标仿真报告徐一凡 1 SAR 原理简介合成孔径雷达(Synthetic Aperture Radar ，简称SAR)是一种高分辨率成像雷达技术。

它利用脉冲压缩技术获得高的距离向分辨率，利用合成孔径原理获得高的方位向分辨率，从而获得大面积高分辨率雷达图像。

SAR 回波信号经距离向脉冲压缩后，雷达的距离分辨率由雷达发射信号带宽决定：2r rC B ρ=，式中r ρ表示雷达的距离分辨率，r B 表示雷达发射信号带宽，C 表示光速。

同为(PT x=,0z =；）， (;)PT R s r = = (2) (;)R s r 就表示任意时刻s 时，目标与雷达的斜距。

一般情况下，0v s s r -<<，于是通过傅里叶技术展开，可将(2)式可近似写为：220(;)()2v R s r r s s r =≈+- (3) 可见，斜距是s r 和的函数，不同的目标，r 也不一样，但当目标距SAR 较远时，在观测带内，可近似认为r 不变，即0r R =。

图2：空间几何关系 (a)正视图 (b)侧视图图2(a)中，Lsar 表示合成孔径长度，它和合成孔径时间Tsar 的关系是Lsar vTsar =。

(b)中，θ∆为雷达天线半功率点波束角，θ为波束轴线与Z 轴的夹角，即波束视角，min R 为近距点距离，max R 为远距点距离，W 为测绘带宽度，它们的关系为：2min (R H tg θθ∆=⋅-） 式中，rect()表示矩形信号，r K 为距离向的chirp 信号调频率，c f 为载频。

雷达回波信号由发射信号波形，天线方向图，斜距，目标RCS ，环境等因素共同决定，若不考虑环境因素，则单点目标雷达回波信号可写成式(6)所示：()()r n n s t wp t n PRT στ∞=-∞=-⋅-∑ (6) 其中，σ表示点目标的雷达散射截面，w 表示点目标天线方向图双向幅度加权，n τ表示载机发射第n 个脉冲时，电磁波再次回到载机时的延时2*(;)n R s r Cτ=，带入式(6)中得： 22(;)/()(exp[(2(;)/)]4 exp[-(;)]exp[2()]r n r r c n t n PRT R s r C s t w rect T j K t n PRT R s r C jR s r j f t n PRT σπππτλ∞=-∞-⋅-=⋅⋅⋅-⋅-⋅⋅-⋅-∑ (7) 式(7)就是单点目标回波信号模型，其中，2exp[(2(;)/)]r j K t n PRT R s r C π-⋅-是）） 为了进行方位向的压缩，方位向的回波数据必须在同一条直线上，也就是说必须校正距离徙动(,)R s r ∆。

由式（10）可知，不同的最近距离r 对应着不同的(,)R s r ∆，因此在时域处理距离徙动会非常麻烦。

因此，对方位向进行傅里叶变换，对距离向不进行变换，得到新的域。

由于方位向的频率即为多普勒频率，所以这个新的域也称为距离多普勒域。

将斜距R 写成多普勒fa 的函数，即(,)a R f r 。

众所周知，对最近距离为r 的点目标P ，回波多普勒a f 是倾斜角θ的函数，即2sin a Vf θλ=，斜距(,)/cos a R f r r θ=，于是22(,)/cos / /1 ()8a a R f r r r r r rf Vθλ===≈+ (11)所以距离多普勒域中的我距离徙动为(,)a R f r ∆=221()8a rf Vλ，可发现它不随慢时间变换，同一最短距离r 对应着相同大小的距离徙动。

因此在距离多普勒域对一个距离徙动校正就是对一组具有相同最短距离的点目标的距离徙动校正，这样可以节省运算量。

为了对距离徙动进行校正，需要得到距离徙动单元，即距离徙动体现在存储单元中的移5 点目标成像matlab 仿真5.1距离多普勒算法距离多普勒算法（RDA ）是在1976年至1978年为民用星载SAR 提出的，它兼顾了成熟、简单、高效和精确等因素，至今仍是使用最广泛的成像算法。

它通过距离和方位上的频域操作，到达了高效的模块化处理要求，同时又具有了一维操作的简便性。

图7示意了RDA 的处理流程。

这里主要讨论小倾斜角及短孔径下的基本RDA 处理框图。

1.当数据处在方位时域时，可通过快速卷积进行距离压缩。

也就是说，距离FFT 后随即进行距离向匹配滤波，再利用距离IFFT 完成距离压缩。

回波信号为：0020(,)[2()/]()exp{-4()/}exp{(-2()/)}r a c r s t s A w t R s c w s s j f R s c j K t R s c ππ=--⨯（14）距离向压缩后的信号为： 000(,){(,)()}[2()/]()exp{4()/}rc t t t r a c s t s IFFT S f s H f A t R s c w s s j f R s c ρπ==--- （15）20(){}exp{}exp{2}||t f f H f rect j j ft K T Kππ=-- （16） 2.通过方位FFT 将数据变换至距离多普勒域，多普勒中心频率估计以及大部分后续操作） ）3200004 (2/)()exp{}s s az s r a s sc f R A p t R c W f f j cπ=--- （20） 5.最后通过方位IFFT 将数据变换回时域，得到压缩后的复图像。

复原后的图像为：30000(,){(,)}(-2/)()4 exp{-}exp{2}c ac s s r a s s t s IFFT S t f A p t R c p s f R j j f s c ππ==⨯ （21）图8 距离多普勒算法流程图5.2 Chirp Scaling算法距离多普勒算法具有诸多优点，但是距离多普勒算法有两点不足：首先，当用较长的核函数提高距离徙动校正（RCMC）精度时，运算量较大；其次，二次距离压缩（SRC）对方位频率的依赖性问题较难解决，从而限制了其对某些大斜视角和长孔径SAR的处理精度。

Chirp Scaling算法避免了RCMC中的插值操作，通过对Chirp信号进行频率调制，实现了对该信号的尺度变换或平移。

图8显示了Chirp Scaling算法处理流程。

这里主要讨论小倾斜角及短孔径下的基本CSA 处理框图。

主要步骤包括四次FFT和三次相位相乘。

1.通过方位向FFT将数据变换到距离多普勒域。

2.通过相位相乘实现Chirp Scaling操作，使所有目标的距离徙动轨迹一致化。

这是第一步相位相乘。

用以改变线调频率尺度的Chirp Scaling二次相位函数为：中6.1 使用最近邻点插值的距离多普勒算法仿真结果本文首先对5个点目标的回波信号进行了仿真，5个点目标构成了矩形的4个顶点和中心，其坐标分别如下，格式为（方位向,距离向,后向反射系数）：0 9750 1100 9750 150 10000 10 10250 1100 10250 1图9的上图是距离向压缩后的图像，从图中可以看到5条回波信号（其中有几条部分重合，但仍能看出来）目标回波信号存在明显的距离徙动，需要进行校正。

图9的下图是通过最近邻点插值法校正后的图像，可以看出图像基本被校正为直线。

图9距离向压缩后最近邻点插值的结果图10为进行方位向压缩后形成的图像，可以明显看出5个点目标，并且5个点目标构成了矩形的四个顶点及其中心。

图10 通过最近邻点插值生成的点目标图像6.2 使用最近邻点插值的距离多普勒算法仿真结果图11上图为通过距离压缩后的图像，图11的下图为通过sinc插值法校正后的图像。

图11 距离向压缩后sinc插值的结果图12为进行方位向压缩后形成的图像，可以明显看出5个点目标，并且5个点目标构成了矩形的四个顶点及其中心。

图12 通过sinc插值生成的点目标图像本文讨论了距离多普勒算法和Chirp Scaling算法，其中距离多普勒算法考虑了最近邻点插值和sinc插值两种插值方法。

距离多普勒算法兼顾了成熟、简单、高效和精确等因素，至今仍被广泛使用，但是距离多普勒算法有两点不足：首先，当用较长的核函数提高距离徙动校正（RCMC）精度时，运算量较大；其次，二次距离压缩（SRC）对方位频率的依赖性问题较难解决，从而限制了其对某些大斜视角和长孔径SAR的处理精度。

最近邻点插值的优点是速度快，该插值的运行时间为2.267137秒，缺点是不够精确；sinc插值的优点是精确，该方法的运行时间为29.148728秒，缺点是速度慢；Chirp Scaling 算法避免了插值运算，提高了速度，运行时间为0.323327秒，但是其算法较为复杂。

%%================================================================%%文件名：NearSAR.m%%作者：徐一凡%%功能：合成孔径雷达距离多普勒算法点目标成像%%================================================================clear;clc;close all;%%================================================================%%常数定义C=3e8; %光速%%雷达参数Fc=1e9; %载频1GHzlambda=C/Fc; %波长的图理解Nslow=2^nextpow2(Nslow); %nextpow2为最靠近2的幂次函数，这里为fft变换做准备sn=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢时间域的时间矩阵PRT=(Xmax-Xmin+Lsar)/V/Nslow; %由于Nslow改变了，所以相应的一些参数也需要更新，周期减小了PRF=1/PRT;ds=PRT;%%快时间域参数设置Tr=5e-6; %脉冲持续时间5usBr=30e6; %chirp频率调制带宽为30MHzKr=Br/Tr; %chirp调频率Fsr=2*Br; %快时域采样频率，为3倍的带宽dt=1/Fsr; %快时域采样间隔Rmin=sqrt((Yc-Y0)^2+H^2);Rmax=sqrt((Yc+Y0)^2+H^2+(Lsar/2)^2);Nfast=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;K=Ntarget; %目标数目N=Nslow; %慢时域的采样数M=Nfast; %快时域的采样数T=Ptarget; %目标矩阵Srnm=zeros(N,M); %生成零矩阵存储回波信号for k=1:1:K %总共K个目标sigma=T(k,3); %得到目标的反射系数Dslow=sn*V-T(k,1); %方位向距离，投影到方位向的距离 R=sqrt(Dslow.^2+T(k,2)^2+H^2); %实际距离矩阵tau=2*R/C; %回波相对于发射波的延时Dfast=ones(N,1)*tm-tau'*ones(1,M); %(t-tau)，其实就是时间矩阵，ones(N,1)和ones(1,M)都是为了将其扩展为矩阵phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M));%相位，公式参见P.96Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*one s(1,M));%由于是多个目标反射的回波，所以此处进行叠加end%%================================================================%%距离-多普勒算法开始%%距离向压缩tic;：Sa_RD(n,m)=Sa_RD(n,M/2);elseif delta_RMC>=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;%%作者：徐一凡%%功能：合成孔径雷达距离多普勒算法点目标成像%%================================================================ clear;clc;close all;%%================================================================ %%常数定义C=3e8; %光速%%雷达参数Fc=1e9; %载频1GHzlambda=C/Fc; %波长%%目标区域参数Xmin=0; %目标区域方位向范围[Xmin,Xmax]Xmax=50;Yc=10000; %成像区域中线Y0=500; %目标区域距离向范围[Yc-Y0,Yc+Y0]%成像宽度为2*Y0%%轨道参数V=100; %SAR的运动速度100 m/sH=5000; %高度 5000 mR0=sqrt(Yc^2+H^2); %最短距离%%天线参数Kr=Br/Tr; %chirp调频率Fsr=2*Br; %快时域采样频率，为3倍的带宽dt=1/Fsr; %快时域采样间隔Rmin=sqrt((Yc-Y0)^2+H^2);Rmax=sqrt((Yc+Y0)^2+H^2+(Lsar/2)^2);Nfast=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,反射系数sigma]Ptarget=[Xmin,Yc-50*DY,1 %点目标位置，这里设置了5个点目标，构成一个矩形以及矩形的中心Xmin+50*DX,Yc-50*DY,1Xmin+25*DX,Yc,1Xmin,Yc+50*DY,1ones(N,1)和ones(1,M)都是为了将其扩展为矩阵phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M));%相位，公式参见P.96Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*one s(1,M));%由于是多个目标反射的回波，所以此处进行叠加end%%================================================================%%距离-多普勒算法开始%%距离向压缩tic;tr=tm-2*Rmin/C;Refr=exp(j*pi*Kr*tr.^2).*(0<tr&tr<Tr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr))));Gr=abs(Sr);%%开始进行距离弯曲补偿正侧视没有距离走动项主要是因为斜距的变化引起回波包络的徙动%%补偿方法：最近邻域插值法，具体为：先变换到距离多普勒域，分别对单个像素点计算出距离徙动量，得到距离徙动量与距离分辨率的比值，%%该比值可能为小数，按照四舍五入的方法近似为整数，而后在该像素点上减去徙动量%%方位向做fft处理再在频域做距离弯曲补偿Sa_RD = ftx(Sr); % 方位向FFT 变为距离多普域进行距离弯曲校正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(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]);xlabel('距离向'),ylabel('方位向'),title('距离向压缩，未校正距离徙动的图像'),kx=linspace(-1/dx/2,1/dx/2,Nfast);%kx域序列%%方位向参数cross-range:y domainTa=300;%时宽300m,合成孔径长度Ba=1;%带宽1(1/m)Ka=Fc/Xc;%调频斜率 Ka=Ba/Ta=Fc/XcNslow=1024;%为了快速运算Y0=200;y=linspace(-Y0,Y0,Nslow);%y域序列:-Y0~Y0dy=2*Y0/Nslow;ky=linspace(-1/dy/2,1/dy/2,Nslow);%ky域序列%%目标几何关系target geometry%x坐标,y坐标,复后向散射系数Ptar=[Xc,0,1+0jXc+50,-50,1+0jXc+50,50,1+0jXc-50,-50,1+0jXc-50,50,1+0j];disp('Position of targets');disp(Ptar)%%生成SAR正交解调后的回波数据Srnm=zeros(Nfast,Nslow);N=size(Ptar,1);%目标个数s1=ifty(scs_xky);%为显示存储数据scs_kxky=fftshift(fft(fftshift(scs_xky)));%距离向FFT（步骤3）srmc_kxky=scs_kxky.*exp(j*pi*(kx.^2'*ones(1,Nslow))./(1+Cs)./Ks...+j*2*pi*Xc*Cs.*(kx'*ones(1,Nslow)));%距离迁移校正&距离向匹配滤波（步骤4）srmc_xky=fftshift(ifft(fftshift(srmc_kxky)));%距离向IFFT（步骤5）f_xky=srmc_xky.*exp(-j*pi*Ks.*Cs.*(1+Cs).*((x-Xc).^2'*ones(1,Nslow))...-j*2*pi*phi0);%消除误差函数，方位向匹配滤波（步骤6）f_xy=fftshift(ifft(fftshift(f_xky).')).';%方位向IFFT（步骤7）%%CSA:7步结束toc;%%为了显示CSA算法的一致RCMC，将s1进行距离向压缩显示p0_x=exp(j*pi*Kr*(x-Xc).^2).*(abs(x-Xc)<Tr/2);%距离向LFM信号p0_kx=fftshift(fft(fftshift(p0_x)));p0_y=exp(-j*pi*Ka*y.^2).*(abs(y)<Ta/2);%方位向LFM信号p0_ky=fftshift(fft(fftshift(p0_y)));s_kxy=fftshift(fft(fftshift(s1)));%距离向FFTsxc_kxy=s_kxy.*(conj(p0_kx).'*ones(1,Nslow));sxc_kxky=fftshift(fft(fftshift(sxc_kxy).')).';%距离压缩后的2D频域信号sxc_xy=fftshift(ifft(fftshift(sxc_kxy)));%距离压缩后的信号sxc_xky=fftshift(fft(fftshift(sxc_xy).')).';%距离压缩后，距离-多普勒域imagesc(255-abs(f_xy));xlabel('方位向'),ylabel('距离向'),title('生成的点目标');。