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

徐一凡

1 SAR原理简介

合成孔径雷达(Synthetic Aperture Radar .简称SAR)是一种高分辨率成像雷达技术。它利用脉冲压缩技术获得高的距离向分辨率.利用合成孔径原理获得高的方位向分辨率.从而获得大面积高分辨率雷达图像。

SAR回波信号经距离向脉冲压缩后.雷达的距离分辨率由雷达发射信号带宽决定：

rC.式中r表示雷达的距离分辨率.Br表示雷达发射信号带宽.C表示光速。同2Br样.SAR回波信号经方位向合成孔径后.雷达的方位分辨率由雷达方位向的多谱勒带宽决定：

ava.式中a表示雷达的方位分辨率.Ba表示雷达方位向多谱勒带宽.va表示方位向BaD.其中D为方位向合成2SAR平台速度。在小斜视角的情况下.方位分辨率近似表示为a孔径的长度。

2 SAR的几何关系

雷达位置和波束在地面覆盖区域的简单几何模型如图1所示。此次仿真考虑的是正侧视的条带式仿真.也就是说倾斜角为零.SAR波束中心和SAR平台运动方向垂直的情况。

图1 雷达数据获取的几何关系

建立坐标系XYZ如图2所示.其中XOY平面为地平面；SAR平台距地平面高H.以速度V沿X轴正向匀速飞行；P点为SAR平台的位置矢量.设其坐标为(x,y,z)； T点为目标的位置矢量.设其坐标为(xT,yT,zT)；由几何关系.目标与SAR平台的斜距为：

RPT(xxT)2(yyT)2(zzT)2(1)

. .

y0,zH,zT0；vs.其中v为平台速度.s为慢时间变量由图可知：令x（slow time）.

假设xTvs.其中s表示SAR平台的x坐标为xT的时刻；再令r与SAR的垂直斜距.重写(1)式为：

PTR(s;r)H2yT2.r表示目标

r2v2(ss0)2 (2) R(s;r)就表示任意时刻s时.目标与雷达的斜距。一般情况下.vss0r.于是通过傅里

叶技术展开.可将(2)式可近似写为：

v2R(s;r)rv(ss0)r(ss0)2 (3)

2r222可见.斜距是s和r的函数.不同的目标.r也不一样.但当目标距SAR较远时.在观测带内.可近似认为r不变.即rR0。

图2：空间几何关系 (a)正视图 (b)侧视图

图2(a)中.Lsar表示合成孔径长度.它和合成孔径时间Tsar的关系是LsarvTsar。(b)中.为雷达天线半功率点波束角.为波束轴线与Z轴的夹角.即波束视角.Rmin为近距点距离.Rmax为远距点距离.W为测绘带宽度.它们的关系为：

RminHtg(2) RmaxHtg(2) (4)

WRmaxRmin. .

3 SAR的回波信号模型

SAR在运动中以一定的周期(1/PRF)发射和接收信号.具体过程如图3所示。发射机以

l的时间发射啁啾脉冲.然后切换天线开关接收回波信号。

图3 雷达发射脉冲串的时序 当雷达不处于发射状态时.它接收3反射回波。发射和接收回波的时间序列如图4所示。在机载情况下.每个回波可以在脉冲发射间隔内直接接收到。但是在星载情况下.由于距离过大.某个脉冲的回波要经过6~10个脉冲间隔才能接收到。这里仿真为了方便.默认为机载情况。

图4 脉冲雷达的发射与接收周期

假设Tr为chirp信号持续时间.下标r表示距离向；PRF为重复频率.PRT为重复周期.

等于1/PRF。接收序列中.n2*R(s;r)表示发射第i个脉冲时.目标回波相对于发射序C列的延时。雷达的发射序列数学表达式为式(5):

s(t)

np(tn*PRT)tjKrtj2fct)eeTr2 （5）

p(t)rect(

式中.rect()表示矩形信号.Kr为距离向的chirp信号调频率.fc为载频。

雷达回波信号由发射信号波形.天线方向图.斜距.目标RCS.环境等因素共同决定.若不考虑环境因素.则单点目标雷达回波信号可写成式(6)所示：

sr(t)

nwp(tnPRTn) (6)

其中.表示点目标的雷达散射截面.w表示点目标天线方向图双向幅度加权.n表示

载机发射第n个脉冲时.电磁波再次回到载机时的延时n2*R(s;r).带入式(6)中得： C. .

sr(t)nwrect(4tnPRT2R(s;r)/C)Tr (7)

exp[jKr(tnPRT2R(s;r)/C)2] exp[-j

R(s;r)]exp[j2fc(tnPRTn)]式(7)就是单点目标回波信号模型.其中.exp[jKr(tnPRT2R(s;r)/C)2]是

chirp分量.它决定距离向分辨率；exp[-j

4 R(s;r)]为多普勒分量.它决定方位向分辨率。

对于任意一个脉冲.回波信号可表示为式(8)所示：

sr(t,s)A0wr(2R(s;r)/C)wa(ssc)exp{j4f0R(s;r)/C} exp{jKr(2R(s;r)/C)}2（8）

我们知道.由于R(s;r)随慢时间s的变化而变化.所以计算机记录到的回波数据存储形

式如图5所示：

图5 目标照射时间内.单个点目标回波能量在信号处理器的二维存储器中的轨迹

4 距离徙动及校正

根据图2可知.在倾斜角为零或很小的时候.目标与雷达的瞬时距离为R(s;r).根据几

何关系可知.R(s;r)r2v2(ss0)2.根据泰勒级数展开可得：

222v2R(s;r)rv(ss0)r(ss0)2 (9)

2r. .

由式(9)可知.不同慢时间对应着不同的R(s;r).并且是一个双曲线形式或者近似为一

个二次形式。如图5所示.同一目标的回波存储在计算机里不在同一直线上.存在距离徙动。从而定义距离徙动量：

v2R(s,r)(ss0)2 （10）

2r为了进行方位向的压缩.方位向的回波数据必须在同一条直线上.也就是说必须校正距离徙动R(s,r)。由式（10）可知.不同的最近距离r对应着不同的R(s,r).因此在时域处理距离徙动会非常麻烦。因此.对方位向进行傅里叶变换.对距离向不进行变换.得到新的域。由于方位向的频率即为多普勒频率.所以这个新的域也称为距离多普勒域。

将斜距R写成多普勒fa的函数.即R(fa,r)。众所周知.对最近距离为r的点目标P.回波多普勒fa是倾斜角的函数.即fa2Vsin.斜距R(fa,r)r/cos.于是

R(fa,r)r/cosr/1sin2 r/1(fa2)2V (11)

1 r()2rfa28V所以距离多普勒域中的我距离徙动为R(fa,r)=()2rfa2.可发现它不随慢时间变换.同一最短距离r对应着相同大小的距离徙动。因此在距离多普勒域对一个距离徙动校正就是对一组具有相同最短距离的点目标的距离徙动校正.这样可以节省运算量。

为了对距离徙动进行校正.需要得到距离徙动单元.即距离徙动体现在存储单元中的移动数值.距离徙动单元可以表示为R(fa,r)/r.这个值通常为一个分数.由于存储单元都是离散的.所以不同通过在存储单元简单的移动得到准确的值。为了得到准确的徙动校正值.通常需要进行插值运算。

本仿真采用了两种插值方法最近邻点插值和sinc插值.下面分别进行介绍。最近邻点插值法的优点是简单而快速.缺点是不够精确。R(fa,r)/r=Nn.其中N为整数部分.n为小数部分.整数部分徙动可以直接通过平移消除.对于小数部分则通过四舍五入的方法变为0或者1.这样就可以得到较为精确的插值。

Sinc插值原理如下：在基带信号下.卷积核是sinc函数

18Vh(x)sinc(x)插值信号为

sin(x) x （12）

g(x)gd(i)sinc(xi)

i （13）

即为所有输入样本的加权平均。

可通过频域来理解.如图6所示.采样信号gd(i)的频谱Gd(f)等于以采样率重复的信

. .

号频谱。为了重建信号g(x).只需要一个周期频谱（如基带周期）.因此需要理想矩形低通滤波器在频域中提取基带频谱（如图6）所示。已知该理想滤波器在时域中是sinc函数。由于频域相乘相当于时域卷积.故插值可以通过与sinc核的卷积来实现。

图6 理想低通滤波器怎样对采样信号进行插值

5 点目标成像matlab仿真

5.1距离多普勒算法

距离多普勒算法（RDA）是在1976年至1978年为民用星载SAR提出的.它兼顾了成熟、简单、高效和精确等因素.至今仍是使用最广泛的成像算法。它通过距离和方位上的频域操作.到达了高效的模块化处理要求.同时又具有了一维操作的简便性。

图7示意了RDA的处理流程。这里主要讨论小倾斜角及短孔径下的基本RDA处理框图。 1.当数据处在方位时域时.可通过快速卷积进行距离压缩。也就是说.距离FFT后随即进行距离向匹配滤波.再利用距离IFFT完成距离压缩。回波信号为：

s0(t,s)A0wr[t2R(s)/c]wa(ssc) exp{-j4f0R(s)/c}exp{jKr(t-2R(s)/c)}距离向压缩后的信号为：

2（14）

src(t,s)IFFTt{S0(ft,s)H(ft)} A0r[t2R(s)/c]wa(ssc)exp{j4f0R(s)/c} （15）

ff2H(ft)rect{}exp{j}exp{j2ft0} （16）

|K|TK 2.通过方位FFT将数据变换至距离多普勒域.多普勒中心频率估计以及大部分后续操作

都在该域进行。方位向傅里叶变换后信号为：

S1(t,fs)FFTs{src(t,s)}2Rrd(fs)]Wa(fsfsc) （17） c4f0R0fs2 exp{-j}exp{j}cKa A0pr[t 3.在距离多普勒域进行随距离时间及方位频率变化的RCMC.该域中同一距离上的一组

目标轨迹相互重合。RCMC将距离徙动曲线拉直到与方位频率轴平行的方向。这里可以采用最近邻点插值法或者sinc插值法.具体插值方法见前面。假设RCMC插值是精确的.信号变为：

. .

2R0)Wa(fsfsc)c （18）

4f0R0fs2 exp{-j}exp{j}cKaS2(t,fs)A0pr(t

4.通过每一距离门上的频域匹配滤波实现方位压缩。为进行方位压缩.将RCMC后的

S2(t,fs)乘以频域匹配滤波器Haz(fs)。

fs2Haz(fs)exp{j}

KaS3(t,fs)S2(t,fs)Haz(fs) （19）

A0pr(t2R0/c)Wa(fsfsc)exp{j4f0R0 （20）

}c 5.最后通过方位IFFT将数据变换回时域.得到压缩后的复图像。复原后的图像为：

sac(t,s)IFFTs{S3(t,fs)} A0pr(t-2R0/c)pa(s) exp{-j4f0R0}exp{j2fscs}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二次相位函数为：

H1(t,fa;Rs)exp[j(fa;RB)a(fa)(t2R(fa;Rs)2)] c （22）

3.通过距离向FFT将数据变到二维频域。

4.通过与参考函数进行相位相乘.同时完成距离压缩、SRC和一致RCMC。这是第二步相位相乘。用于距离压缩.距离徙动校正的相位函数写为：

H2(fr,fa;Rs)exp[j

1fr2](fa;RB)[1a(fa)]4Rsa(fa) exp[j]fr]c （23）

5.通过距离向IFFT将数据变回到距离多普勒域。

6.通过与随距离变化的匹配滤波器进行相位相乘.实现方位压缩。此外.由于步骤2中的Chirp Scaling操作.相位相乘中还需要附加一项相位校正。这是第三步相位相乘。补偿由Chirp Scaling引起的剩余相位函数是：

H2(tr,fa;RB)exp[j

2RBVfaM2fa2]exp[j(fa;RB)]

（24）

7.最后通过方位向IFFT将数据变回到二维时域.即SAR图像域。

图8Chirp Scaling算法流程图

. .

简而言之.R-D算法是将徙动曲线逐一校正.CS算法是以某一徙动曲线为参考.在

Doppler域内消除不同距离门的徙动曲线的差异.令这些曲线成为一组相互“平行”的曲线.然后在二维频率域内统一的去掉距离徙动。通俗一点就是.RD算法是将弯曲的信号一根根掰直.而CS算法是先把所有信号都掰得一样弯.然后再统一掰直。

6 仿真结果

6.1 使用最近邻点插值的距离多普勒算法仿真结果

本文首先对5个点目标的回波信号进行了仿真.5个点目标构成了矩形的4个顶点和中心.其坐标分别如下.格式为（方位向,距离向,后向反射系数）： 0 9750 1 100 9750 1 50 10000 1 0 10250 1 100 10250 1

图9的上图是距离向压缩后的图像.从图中可以看到5条回波信号（其中有几条部分重合.但仍能看出来）目标回波信号存在明显的距离徙动.需要进行校正。图9的下图是通过最近邻点插值法校正后的图像.可以看出图像基本被校正为直线。

距离向压缩，未校正距离徙动的图像-400-200方位向02004000.950.960.970.980.991距离向1.011.021.031.041.05x 104距离向压缩，校正距离徙动后的图像-400-200方位向02004000.950.960.970.980.991距离向1.011.021.031.041.05x 104

图9距离向压缩后最近邻点插值的结果

图10为进行方位向压缩后形成的图像.可以明显看出5个点目标.并且5个点目标构成了矩形的四个顶点及其中心。

方位向压缩后的图像-400-300-200-100X: 9777 Y: -5.347Index: -5.635e+04RGB: 0.698, 0.698, 0.698X: 9777 Y: 95.5Index: -6.914e+04RGB: 0.619, 0.619, 0.619方位向0X: 9998 Y: 48.56Index: -9.081e+04RGB: 0.508, 0.508, 0.508X: 1.022e+04 Y: -0.1306Index: -6.263e+04RGB: 0.667, 0.667, 0.667X: 1.022e+04 Y: 97.24Index: -7.73e+04RGB: 0.587, 0.587, 0.5871002003004000.950.960.970.980.991距离向1.011.021.031.041.05x 104

图10 通过最近邻点插值生成的点目标图像

. .

6.2 使用最近邻点插值的距离多普勒算法仿真结果

图11上图为通过距离压缩后的图像.图11的下图为通过sinc插值法校正后的图像。

距离向压缩，未校正距离徙动的图像-400-200方位向02004000.950.960.970.980.9911.01距离向1.021.031.041.05x 104距离向压缩，校正距离徙动后的图像-400-200方位向02004000.950.960.970.980.9911.01距离向1.021.031.041.05x 104

图11 距离向压缩后sinc插值的结果

图12为进行方位向压缩后形成的图像.可以明显看出5个点目标.并且5个点目标构成了矩形的四个顶点及其中心。

方位向压缩后的图像-400-300-200-100X: 9777 Y: -1.869Index: -4.011e+04RGB: 0.46, 0.46, 0.46X: 9775 Y: 98.98Index: -3.731e+04RGB: 0.508, 0.508, 0.508方位向0X: 1e+04 Y: 48.56Index: -4.628e+04RGB: 0.381, 0.381, 0.381X: 1.022e+04 Y: -0.1306Index: -3.207e+04RGB: 0.571, 0.571, 0.571X: 1.022e+04 Y: 97.24Index: -4.445e+04RGB: 0.413, 0.413, 0.4131002003004000.950.960.970.980.991距离向1.011.021.031.041.05x 104

图12 通过sinc插值生成的点目标图像

6.3 Chirp Scaling算法仿真结果

同样.在Chirp Scaling中.对5个点目标的回波信号进行了仿真.5个点目标构成了矩形的4个顶点和中心.其坐标分别如下.格式为（方位向,距离向,后向反射系数）： 1200 0 1 1250 -50 1 1250 50 1 1150 -50 1

. .

1150 50 1 图13是仿真的雷达回波信号图

仿真出来的信号50100150200距离向250300350400450500100200300400500600方位向7008009001000

图13 仿真出来的SAR回波信号

图14是经过第一次相位校正之后.通过距离向压缩后的距离时域-方位时域信号图（Chirp Scaling算法的七个步骤中并不包含该信号.该信号是将步骤2之后的信号通过方位向傅里叶逆变换.再进行距离向压缩得到的.只为了验证原理）。按照理论.该图中所有点的距离徙动都应该一样。从图中大致可看出.五个点的距离徙动是差不多的。

Chirp Scaling后、经过距离向压缩，距离徙动一致50100150200距离向250300350400450500100200300400500600方位向7008009001000

图14 Chirp Scaling第2步之后、经过距离向压缩得到的图

图15为步骤5之后.信号距离压缩.距离徙动校正之后的距离多普勒域中的信号图。

. .

消除距离徙动后的信号50100150200距离向250300350400450500100200300400500600方位向7008009001000

图15 距离徙动校正之后的图

图16为步骤6之后.消除相位偏移的图。

相位校正后的信号50100150200距离向250300350400450500100200300400500600方位向7008009001000

图16 消除相位偏移的图

图17为通过Chirp Scaling算法生成的点目标图像

. .

生成的点目标50100X: 384 Y: 170Index: 175.4RGB: 0.635, 0.635, 0.635X: 2 Y: 170Index: 120.5RGB: 0.381, 0.381, 0.381150200距离向X: 513 Y: 2Index: 195.8RGB: 0.73, 0.73, 0.73250300X: 386 Y: 340Index: 155.4RGB: 0., 0., 0.X: 1 Y: 340Index: 153.7RGB: 0., 0., 0.350400450500100200300400500方位向6007008009001000

图17通过Chirp Scaling算法生成的点目标图像

6.4 几种算法比较

本文讨论了距离多普勒算法和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; %载频1GHz lambda=C/Fc; %波长 %%目标区域参数

Xmin=0; %目标区域方位向范围[Xmin,Xmax] Xmax=50;

Yc=10000; %成像区域中线

Y0=500; %目标区域距离向范围[Yc-Y0,Yc+Y0] %成像宽度为2*Y0 %%轨道参数

V=100; %SAR的运动速度100 m/s

. .

H=5000; %高度 5000 m R0=sqrt(Yc^2+H^2); %最短距离 %%天线参数

D=4; %方位向天线长度

Lsar=lambda*R0/D; %SAR合成孔径长度.《合成孔径雷达成像——算法与实现》P.100 Tsar=Lsar/V; %SAR照射时间 %%慢时间域参数

Ka=-2*V^2/lambda/R0; %多普勒频域调频率P.93 Ba=abs(Ka*Tsar); %多普勒频率调制带宽

PRF=Ba; %脉冲重复频率.PRF其实为多普勒频率的采样率.又为复频率.所以等于Ba.P.93

PRT=1/PRF; %脉冲重复时间

ds=PRT; %慢时域的时间步长 Nslow=ceil((Xmax-Xmin+Lsar)/V/ds); %慢时域的采样数.ceil为取整函数.结合P.76的图理解

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; %脉冲持续时间5us

Br=30e6; %chirp频率调制带宽为30MHz 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,1

. .

Xmin+25*DX,Yc,1 Xmin,Yc+50*DY,1

Xmin+50*DX,Yc+50*DY,1]; disp('Parameters:') %参数显示

disp('Sampling Rate in fast-time domain');disp(Fsr/Br) disp('Sampling Number in fast-time domain');disp(Nfast) disp('Sampling Rate in slow-time domain');disp(PRF/Ba) disp('Sampling Number in slow-time domain');disp(Nslow) disp('Range Resolution');disp(DY)

disp('Cross-range Resolution');disp(DX) disp('SAR integration length');disp(Lsar) disp('Position of targets');disp(Ptarget)

%%================================================================ %%生成回波信号

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.96

Srnm=Srnm+sigma*exp(j*phase).*(0%%================================================================ %%距离-多普勒算法开始 %%距离向压缩 tic;

tr=tm-2*Rmin/C;

Refr=exp(j*pi*Kr*tr.^2).*(0Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr)))); Gr=abs(Sr);

%%开始进行距离弯曲补偿正侧视没有距离走动项 主要是因为斜距的变化引起回波包络的徙动

%%补偿方法：最近邻域插值法.具体为：先变换到距离多普勒域.分别对单个像素点计算出距离徙动量.得到距离徙动量与距离分辨率的比值.

%%该比值可能为小数.按照四舍五入的方法近似为整数.而后在该像素点上减去徙动量

. .

%%方位向做fft处理 再在频域做距离弯曲补偿

Sa_RD = ftx(Sr); % 方位向FFT 变为距离多普域进行距离弯曲校正 %距离徙动运算,由于是正侧视 .fdc=0,只需要进行距离弯曲补偿。 Kp=1; %计算或者预设预滤波比 h = waitbar(0,'最近邻域插值'); %%首先计算距离迁移量矩阵

for n=1:N %总共有N个方位采样

for m=1: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);%距离徒动量的小数部分 if m+round(RMC)>M %判断是否超出边界 Sa_RD(n,m)=Sa_RD(n,M/2); else

if 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)); end end end

waitbar(n/N) end

close(h)

%========================

Sr_rmc=iftx(Sa_RD); %%距离徙动校正后还原到时域 Ga = abs(Sr_rmc); %%方位向压缩 ta=sn-Xmin/V;

Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)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('距离向压缩.未校正距离徙动的图像'), subplot(212);

imagesc(row,col,255-Ga); %距离向压缩.校正距离徙动后的图像 axis([Yc-Y0,Yc+Y0,Xmin-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 %%作者：徐一凡

%%功能：合成孔径雷达距离多普勒算法点目标成像

%%================================================================ clear;clc;close all;

%%================================================================ %%常数定义

C=3e8; %光速 %%雷达参数

Fc=1e9; %载频1GHz lambda=C/Fc; %波长 %%目标区域参数

Xmin=0; %目标区域方位向范围[Xmin,Xmax] Xmax=50;

Yc=10000; %成像区域中线

Y0=500; %目标区域距离向范围[Yc-Y0,Yc+Y0] %成像宽度为2*Y0 %%轨道参数

V=100; %SAR的运动速度100 m/s H=5000; %高度 5000 m R0=sqrt(Yc^2+H^2); %最短距离 %%天线参数

D=4; %方位向天线长度

Lsar=lambda*R0/D; %SAR合成孔径长度.《合成孔径雷达成像——算法与实现》P.100 Tsar=Lsar/V; %SAR照射时间 %%慢时间域参数

Ka=-2*V^2/lambda/R0; %多普勒频域调频率P.93 Ba=abs(Ka*Tsar); %多普勒频率调制带宽

PRF=Ba; %脉冲重复频率.PRF其实为多普勒频率的采样率.

. .

又为复频率.所以等于Ba.P.93

PRT=1/PRF; %脉冲重复时间

ds=PRT; %慢时域的时间步长 Nslow=ceil((Xmax-Xmin+Lsar)/V/ds); %慢时域的采样数.ceil为取整函数.结合P.76的图理解

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; %脉冲持续时间5us

Br=30e6; %chirp频率调制带宽为30MHz 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,1 Xmin+25*DX,Yc,1 Xmin,Yc+50*DY,1

Xmin+50*DX,Yc+50*DY,1]; disp('Parameters:') %参数显示

disp('Sampling Rate in fast-time domain');disp(Fsr/Br) disp('Sampling Number in fast-time domain');disp(Nfast) disp('Sampling Rate in slow-time domain');disp(PRF/Ba) disp('Sampling Number in slow-time domain');disp(Nslow) disp('Range Resolution');disp(DY)

disp('Cross-range Resolution');disp(DX)

. .

disp('SAR integration length');disp(Lsar) disp('Position of targets');disp(Ptarget)

%%================================================================ %%生成回波信号

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.96

Srnm=Srnm+sigma*exp(j*phase).*(0%%================================================================ %%距离-多普勒算法开始 %%距离向压缩 tic;

tr=tm-2*Rmin/C;

Refr=exp(j*pi*Kr*tr.^2).*(0Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr)))); Gr=abs(Sr);

%%开始进行距离弯曲补偿正侧视没有距离走动项 主要是因为斜距的变化引起回波包络的徙动

%%补偿方法：最近邻域插值法.具体为：先变换到距离多普勒域.分别对单个像素点计算出距离徙动量.得到距离徙动量与距离分辨率的比值.

%%该比值可能为小数.按照四舍五入的方法近似为整数.而后在该像素点上减去徙动量 %%方位向做fft处理 再在频域做距离弯曲补偿

Sa_RD = ftx(Sr); % 方位向FFT 变为距离多普域进行距离弯曲校正 %距离徙动运算,由于是正侧视 .fdc=0,只需要进行距离弯曲补偿。 Kp=1; %计算或者预设预滤波比 %%第二种方法进行截断sinc插值进行距离徒动校正 h = waitbar(0,'Sinc插值'); P=4;%4点sinc插值

RMCmaxtix = zeros(N,M); for n=1:N for m=P:M

. .

delta_R = (1/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); %距离徒动量的小数部分 for i = -P/2:P/2-1

if m+RMC+i>M %判断是否超出边界

RMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,M)*sinc(pi*(-i+RMC)); else

RMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,m+round(RMC)+i)*sinc(pi*(-i+delta_RMC));

end end end

waitbar(n/N) end

close(h)

%========================

Sr_rmc=iftx(RMCmaxtix); %%距离徙动校正后还原到时域 Ga = abs(Sr_rmc); %%方位向压缩 ta=sn-Xmin/V;

Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)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('距离向压缩.未校正距离徙动的图像'), subplot(212);

imagesc(row,col,255-Ga); %距离向压缩.校正距离徙动后的图像 axis([Yc-Y0,Yc+Y0,Xmin-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('方位向压缩后的图像'),

%% Chirp Scaling算法 %徐一凡

clear all;clc;

%%距离向参数range:x domain Tr=200;%时宽200m Br=1;%带宽1

Kr=Br/Tr;%调频斜率 Fc=4;%载频4

Nfast=512;%为了快速运算

Xc=1200;X0=150;%定义距离向范围

x=Xc+linspace(-X0,X0,Nfast);%x域序列:Xc-X0~Xc+X0 dx=2*X0/Nfast;%定义步长

kx=linspace(-1/dx/2,1/dx/2,Nfast);%kx域序列

%%方位向参数cross-range:y domain Ta=300;%时宽300m,合成孔径长度 Ba=1;%带宽1(1/m)

Ka=Fc/Xc;%调频斜率 Ka=Ba/Ta=Fc/Xc Nslow=1024;%为了快速运算 Y0=200;

y=linspace(-Y0,Y0,Nslow);%y域序列:-Y0~Y0 dy=2*Y0/Nslow;

ky=linspace(-1/dy/2,1/dy/2,Nslow);%ky域序列

%%目标几何关系target geometry %x坐标,y坐标,复后向散射系数 Ptar=[Xc,0,1+0j Xc+50,-50,1+0j Xc+50,50,1+0j Xc-50,-50,1+0j Xc-50,50,1+0j];

disp('Position of targets');disp(Ptar) %%生成SAR正交解调后的回波数据 Srnm=zeros(Nfast,Nslow); N=size(Ptar,1);%目标个数

h = waitbar(0,'SAR回波生成'); for i=1:1:N

. .

xn=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)waitbar(i/N) end

close(h) tic;

%%数据准备

phi0=-x'*sqrt(Fc^2-ky.^2);

phi1=-Fc*x'*(1./sqrt(Fc^2-ky.^2));

phi2=1/2*x'*(ky.^2./(Fc^2-ky.^2).^1.5); Cs=ones(Nfast,1)*(Fc./sqrt(Fc^2-ky.^2)-1); Ks=1./(1/Kr-2*phi2);

%%CSA:7步 开始

s_xky=fftshift(fft(fftshift(Srnm).')).';%方位向FFT(步骤1)

scs_xky=s_xky.*exp(j*pi*Cs.*Ks.*(x'*ones(1,Nslow)-Xc*(1+Cs)).^2);%Chirp Scaling（步骤2）

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)p0_y=exp(-j*pi*Ka*y.^2).*(abs(y)s_kxy=fftshift(fft(fftshift(s1)));%距离向FFT sxc_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).')).';%距离压缩后.距离-多普勒域 %%结果显示 figure(1)

colormap(gray);

imagesc(255-abs(Srnm));

xlabel('方位向'),ylabel('距离向'), title('仿真出来的信号');

figure(2)

colormap(gray);

imagesc(255-abs(sxc_xy));

xlabel('方位向'),ylabel('距离向'),

title('Chirp Scaling后、经过距离向压缩.距离徙动一致');

figure(3)

colormap(gray);

imagesc(255-abs(srmc_xky));

xlabel('方位向'),ylabel('距离向'), title('消除距离徙动后的信号');

figure(4)

colormap(gray);

imagesc(255-abs(f_xky));

xlabel('方位向'),ylabel('距离向'), title('相位校正后的信号');

figure(5)

colormap(gray);

imagesc(255-abs(f_xy));

xlabel('方位向'),ylabel('距离向'), title('生成的点目标');

. .