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

1、SAR图像点目标仿真报告徐一凡1 SAR 原理简介 合成孔径雷达(Synthetic Aperture Radar .简称 SAR)是一种高分辨率成像雷达技术。vaa . 式中Baa 表示雷达的方位分辨率Ba 表示雷达方位向多谱勒带宽 . va 表示方位向它利用脉冲压缩技术获得高的距离向分辨率 . 利用合成孔径原理获得高的方位向分辨率 . 从 而获得大面积高分辨率雷达图像。SAR回波信号经距离向脉冲压缩后 . 雷达的距离分辨率由雷达发射信号带宽决定：r 2CBr . 式中r 表示雷达的距离分辨率 . Br 表示雷达发射信号带宽 . C 表示光速。同样.SAR 回波信号经方位向合成孔径后 . 雷

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

3、.设其坐标为 (xT , yT , zT ) ；由几何关系 .目标与 SAR平台的斜距为：uuurR PT(x xT )2 (y yT)2 (z zT )2(1)由图可知： y 0, z H ,zT 0 ；令 x v s. 其中 v 为平台速度 .s 为慢时间变量( slow.rtime ).假设 xT vs.其中 s表示 SAR平台的 x 坐标为 xT 的时刻；再令 r H2 yT2 表示目标与 SAR的垂直斜距 . 重写 (1) 式为：uuurPTR(s;r)r 2 v2 (s s0)2(2)R(s;r)就表示任意时刻 s时.目标与雷达的斜距。 一般情况下 .v s s0 叶技术展开 .可

4、将(2) 式可近似写为：2R(s;r) r2 v2 (s s0)2 r v (s s0)22rr . 于是通过傅里(3)可见.斜距是 s和r的函数.不同的目标 . r也不一样.但当目标距 SAR较远时.在观测带内 .可近似认为 r 不变 . 即 r R0 。图 2：空间几何关系 (a) 正视图 (b) 侧视图图 2(a) 中. Lsar表示合成孔径长度 . 它和合成孔径时间 Tsar的关系是 Lsar vTsar。(b) 中. 为雷达天线半功率点波束角 . 为波束轴线与 Z 轴的夹角 . 即波束视角 . Rmin 为近 距点距离 . Rmax为远距点距离 .W 为测绘带宽度 . 它们的关系为：

5、Rmin H tg(2)Rmax H tg(2)(4)W Rmax Rmin3 SAR的回波信号模型SAR在运动中以一定的周期 ( 1/ PRF )发射和接收信号 .具体过程如图 3 所示。发射机以载机发射第 n 个脉冲时 . 电磁波再次回到载机时的延时n2* RC(s;r ) .带入式(6)中得：l 的时间发射啁啾脉冲 . 然后切换天线开关接收回波信号当雷达不处于发射状态时 . 它接收 3 反射回波。发射和接收回波的时间序列如图 4 所示。 在机载情况下 . 每个回波可以在脉冲发射间隔内直接接收到。 但是在星载情况下 . 由于距离过 大.某个脉冲的回波要经过 610个脉冲间隔才能接收到。这里

6、仿真为了方便 . 默认为机载情 况。图4 脉冲雷达的发射与接收周期假设Tr为 chirp 信号持续时间 .下标 r 表示距离向； PRF为重复频率 .PRT为重复周期 等于1/ PRF 。接收序列中 . n 2* R( s; r )表示发射第 i个脉冲时 .目标回波相对于发射序 nC列的延时。雷达的发射序列数学表达式为式 (5):5)p(t)rect( t )ej Krt2ej2 fct Tr式中 . rect(g) 表示矩形信号s(t) p(t n* PRT) nKr 为距离向的 chirp 信号调频率 . f c为载频。雷达回波信号由发射信号波形 .天线方向图 .斜距 .目标 RCS.环境

7、等因素共同决定 .若不考虑环境因素 . 则单点目标雷达回波信号可写成式 (6) 所示：sr (t)wp(t n PRT n )n(6)其中 . 表示点目标的雷达散射截面 . w表示点目标天线方向图双向幅度加权 . n表示t n PRT 2R(s;r)/Csr (t)w rect( )nTrexpj Kr(t n PRT 2R(s;r)/C)2(7)4exp- j 4 R(s;r) exp j2 fc(t n PRT n)式(7)就是单点目标回波信号模型 .其中. exp j Kr(t n PRT 2R(s;r)/C)2是4chirp 分量.它决定距离向分辨率； exp- j R(s;r) 为多

8、普勒分量 . 它决定方位向分辨率。对于任意一个脉冲 . 回波信号可表示为式 (8) 所示：sr(t,s) A0wr( 2R(s;r)/ C)wa(s sc) exp j4 f0R(s;r)/Cr 0 r a 2 c 0(8)exp j Kr( 2R(s;r) /C)2我们知道 .由于 R(s;r)随慢时间 s 的变化而变化 .所以计算机记录到的回波数据存储形式如图 5 所示：图 5 目标照射时间内 . 单个点目标回波能量在信号处理器的二维存储器中的轨迹4 距离徙动及校正根据图 2可知 .在倾斜角为零或很小的时候 .目标与雷达的瞬时距离为 R(s;r). 根据几何关系可知 . R(s;r) r2

9、. 根据泰勒级数展开可得：v2 (s s0)222 2 2 v 2R(s;r) r v (s s0)r (s s0)(9)2r由式(9) 可知. 不同慢时间对应着不同的 R(s;r). 并且是一个双曲线形式或者近似为一 个二次形式。如图 5所示. 同一目标的回波存储在计算机里不在同一直线上.存在距离徙动。从而定义距离徙动量：2v2R(s,r) (s s0)2(10)2r为了进行方位向的压缩 .方位向的回波数据必须在同一条直线上. 也就是说必须校正距离徙动 R(s, r ) 。由式( 10)可知 . 不同的最近距离 r 对应着不同的 R(s,r). 因此在时域 处理距离徙动会非常麻烦。因此 .

10、对方位向进行傅里叶变换 . 对距离向不进行变换 . 得到新的 域。由于方位向的频率即为多普勒频率 . 所以这个新的域也称为距离多普勒域。将斜距 R写成多普勒 fa 的函数 . 即 R( fa,r)。众所周知 . 对最近距离为 r 的点目标 P.2V回波多普勒 fa 是倾斜角 的函数 .即 fasin .斜距 R(fa,r) r / cos . 于是R(fa,r)r /cos r / 1 sin 2x1281(V)2rfa2(11)1 所以距离多普勒域中的我距离徙动为R(fa,r)=1( )2rfa2. 可发现它不随慢时间变换 . 同一最短距离 r 对应着相同大小的距离徙动。因此在距离多普勒域对

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

12、徙动可以直接通过平移消除 . 对于小数部分则通过四舍五入的方 法变为 0 或者 1. 这样就可以得到较为精确的插值。Sinc 插值原理如下：在基带信号下 . 卷积核是 sinc 函数h(x) sinc(x) sin( x)(12)插值信号为g(x)g d ( i )sin c(x i)(13)i即为所有输入样本的加权平均。可通过频域来理解 .如图 6所示.采样信号 gd(i)的频谱 Gd( f )等于以采样率重复的信号频谱。为了重建信号 g(x) .只需要一个周期频谱(如基带周期).因此需要理想矩形低通滤波器在频域中提取基带频谱(如图6)所示。已知该理想滤波器在时域中是sinc 函数。17)5

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

14、/ cwa(s sc)2 (14)exp- j4 f0R(s)/cexp j Kr(t-2R(s)/c)2src (t,s)距离向压缩后的信号为：IFFTtS0( ft,s)H(ft)(15)A0 rt 2R(s)/cwa(s sc)exp j4 f0R(s) /cH( ft)rectf f 2|Kf|Texp j fK exp j2 ft016)2. 通过方位 FFT将数据变换至距离多普勒域 . 多普勒中心频率估计以及大部分后续操作 都在该域进行。方位向傅里叶变换后信号为：S1(t, fs) FFTs src (t,s)Kfsa2A0prt 2Rrd(fs)Wa(fs fsc )ccexp-

15、 j 4 f0R0exp jc3. 在距离多普勒域进行随距离时间及方位频率变化的RCMC该. 域中同一距离上的一组目标轨迹相互重合。 RCMC将距离徙动曲线拉直到与方位频率轴平行的方向。这里可以采用 最近邻点插值法或者 sinc 插值法 . 具体插值方法见前面。假设 RCMC插值是精确的 . 信号变 为：S2(t, fs) A0 pr (t0)Wa( fs fsc)c4 f Rf 2( 18)exp- j 4 f0R0exp j fs cKa4. 通过每一距离门上的频域匹配滤波实现方位压缩。为进行方位压缩. 将 RCMC后的S2(t, fs) 乘以频域匹配滤波器 Haz(fs)。fs2Haz(

16、 fs) exp j s (19)KaS3(t, fs) S2(t, fs)Haz(fs)4 f0R0 (20)A0pr(t 2R0/ c)Wa(fs fsc)exp j 0 0c5. 最后通过方位 IFFT 将数据变换回时域 . 得到压缩后的复图像。复原后的图像为：sac(t,s) IFFTsS3(t, fs)A0 pr (t - 2R0 / c) pa(s)(21)4 f Rexp- j 0 0exp j2 fsc s图 8 距离多普勒算法流程图5.2 Chirp Scaling算法距离多普勒算法具有诸多优点 . 但是距离多普勒算法有两点不足：首先 .当用较长的核 函数提高距离徙动校正(

17、RCM)C精度时 . 运算量较大；其次 . 二次距离压缩( SRC)对方位频 率的依赖性问题较难解决 . 从而限制了其对某些大斜视角和长孔径SAR的处理精度。Chirp Scaling 算法避免了 RCMC中的插值操作 .通过对 Chirp 信号进行频率调制 . 实现 了对该信号的尺度变换或平移。图 8 显示了 Chirp Scaling 算法处理流程。这里主要讨论小倾斜角及短孔径下的基本 CSA处理框图。主要步骤包括四次 FFT 和三次相位相乘。1. 通过方位向 FFT 将数据变换到距离多普勒域。2. 通过相位相乘实现 Chirp Scaling 操作 . 使所有目标的距离徙动轨迹一致化。这

18、是第 一步相位相乘。用以改变线调频率尺度的 Chirp Scaling 二次相位函数为：2R( fa;Rs) 2H1(t, fa;Rs) exp j ( fa;RB)a( fa)(ta s )2(22)c3. 通过距离向 FFT 将数据变到二维频域。4. 通过与参考函数进行相位相乘 . 同时完成距离压缩、 SRC和一致 RCMC。这是第二步相 位相乘。用于距离压缩 . 距离徙动校正的相位函数写为：H2( fr, fa;Rs) exp j(fa;RB)1 a( fa)exp j 4 Rsa( fa)frfr223)5. 通过距离向 IFFT 将数据变回到距离多普勒域。6. 通过与随距离变化的匹配

19、滤波器进行相位相乘 .实现方位压缩。此外 .由于步骤 2 中 的 Chirp Scaling 操作 . 相位相乘中还需要附加一项相位校正。这是第三步相位相乘。补偿 由 Chirp Scaling 引起的剩余相位函数是：2H2(tr, fa;RB) exp j 2V RB faM2 fa 2 exp j (fa;RB)(24)7. 最后通过方位向 IFFT 将数据变回到二维时域 . 即 SAR图像域。图 8Chirp Scaling 算法流程图简而言之 .R-D 算法是将徙动曲线逐一校正 .CS算法是以某一徙动曲线为参考 . 在 Doppler 域内消除不同距离门的徙动曲线的差异 . 令这些曲线

20、成为一组相互“平行”的曲 线. 然后在二维频率域内统一的去掉距离徙动。通俗一点就是.RD 算法是将弯曲的信号一根根掰直 .而 CS算法是先把所有信号都掰得一样弯 .然后再统一掰直。6 仿真结果6.1 使用最近邻点插值的距离多普勒算法仿真结果本文首先对 5 个点目标的回波信号进行了仿真 .5 个点目标构成了矩形的 4 个顶点和中 心. 其坐标分别如下 . 格式为（方位向 , 距离向 , 后向反射系数） ：09750110097501501000010102501100102501图 9 的上图是距离向压缩后的图像 .从图中可以看到 5 条回波信号（其中有几条部分重 合. 但仍能看出来）目标回波信

21、号存在明显的距离徙动. 需要进行校正。图 9 的下图是通过最近邻点插值法校正后的图像 . 可以看出图像基本被校正为直线。-40000200向位方4000.-40000200向位方4000.距离 向 压 缩 ，未 校正 距离 徙 动 的 图 像距离 向 压 缩 ，校 正距 离徙 动 后 的 图 像1.054x 101.054x 10图 9 距离向压缩后最近邻点插值的结果图 10 为进行方位向压缩后形成的图像 . 可以明显看出 5 个点目标 . 并且 5 个点目标构成 了矩形的四个顶点及其中心。方位向 压 缩后的 图像向 位 方距 离 向4x 10图 10 通过最近邻点插值生成的点目标图像6.2

22、使用最近邻点插值的距离多普勒算法仿真结果图 11 上图为通过距离压缩后的图像 . 图 11 的下图为通过 sinc 插值法校正后的图像。距离向压缩，未校正距离徙动的图像-400 -200向位 方02004000.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 1.05距 离 向x 10 4距 离 向x 10 4向位图 11 距离向压缩后 sinc 插值的结果-400-300-200-100向0位方100方位向 压 缩后的 图像2003004000.95 0.96 0.97 0.980.991距离向1.011.02 1.031.04x 101.054图

23、12 为进行方位向压缩后形成的图像 . 可以明显看出 5 个点目标 . 并且 5 个点目标构成 了矩形的四个顶点及其中心。图 12 通过 sinc 插值生成的点目标图像6.3 Chirp Scaling 算法仿真结果同样.在 Chirp Scaling 中.对 5个点目标的回波信号进行了仿真 .5 个点目标构成了矩 形的 4 个顶点和中心 . 其坐标分别如下 . 格式为（方位向 , 距离向 , 后向反射系数） ：1200 0 11250-50112505011150-5011150501图 13 是仿真的雷达回波信号图向离距仿真出来的信号图 13 仿真出来的 SAR回波信号图 14 是经过第一

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

25、徙动后的信号50100150200250300350400450500100200 300 400 500 600 700 800 900 1000 方位向图 15 距离徙动校正之后的图 图 16 为步骤 6 之后 . 消除相位偏移的图。相位校正后的 信号向 离距350400450500100 200 300 400 500 600 700 800 900 1000 方位向图 16 消除相位偏移的图 图 17 为通过 Chirp Scaling 算法生成的点目标图像50100150Index: 175.4RGB: 0.635, 0.635, 0.635Index: 120.5RGB: 0.38

26、1, 0.381, 0.381200250300X: 386 Y: 340Index: 155.4RGB: 0.54, 0.54, 0.54X: 641 Y: 340Index: 153.7RGB: 0.54, 0.54, 0.54350400450500100 200 300 400 500 600 700 800 900 1000图 17 通过 Chirp Scaling 算法生成的点目标图像6.4 几种算法比较本文讨论了距离多普勒算法和 Chirp Scaling 算法.其中距离多普勒算法考虑了最近邻 点插值和 sinc 插值两种插值方法。距离多普勒算法兼顾了成熟、简单、高效和精确等因素

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

28、SAR.m%作者：徐一凡%功能：合成孔径雷达距离多普勒算法点目标成像%= clear;clc;close all;%=%常数定义C=3e8;%光速%雷达参数Fc=1e9;%载频1GHzlambda=C/Fc;%波长%目标区域参数Xmin=0;%目标区域方位向范围Xmin,XmaxXmax=50;Yc=10000;%成像区域中线Y0=500;%目标区域距离向范围Yc-Y0,Yc+Y0%成像宽度为 2*Y0%轨道参数V=100;%SARH=5000; %R0=sqrt(Yc2+H2); % %天线参数D=4; %的运动速度 100 m/s 高度 5000 m 最短距离方位向天线长度Lsar=lam

29、bda*R0/D; %SAR 合成孔径长度 . 合成孔径雷达成像算法与实现P.100Ka=-2*V2/lambda/R0; %多普勒频域调频率 P.93Tsar=Lsar/V; %SAR 照射时间 %慢时间域参数Ba=abs(Ka*Tsar);%PRF=Ba;%多普勒频率调制带宽脉冲重复频率 .PRF其实为多普勒频率的采样率脉冲重复时间慢时域的时间步长又为复频率 . 所以等于 Ba.P.93PRT=1/PRF;%ds=PRT;%Nslow=ceil(Xmax-Xmin+Lsar)/V/ds); %慢时域的采样数 .ceil 为取整函数 . 结合 P.76的图理解为最靠近 2的幂次函数 . 这里

30、为Nslow=2nextpow2(Nslow); %nextpow2 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;%Br=30e6;%chirpKr=Br/Tr;%chirpFsr=2*Br;%dt=1/Fsr;%Rmin=sqrt(Yc-Y0)2+H2);Rmax=sqrt(Yc+Y0)2+H2

31、+(Lsar/2)2);Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%Nfast=2nextpow2(Nfast); %脉冲持续时间 5us频率调制带宽为 30MHz 调频率快时域采样频率 . 为3倍的带宽快时域采样间隔快时域的采样数量更新为 2的幂次 . 方便进行 fft 变换DY=C/2/Br;%DX=D/2;%点目标参数设置Ntarget=5;%点目标格式 x,y,反射系数 sigmaPtarget=Xmin,Yc-50*DY,1 %dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; %Fsr=1/dt;%分辨率参数设置tm=linspace(

32、2*Rmin/C,2*Rmax/C+Tr,Nfast); % 快时域的离散时间矩阵更新间隔距离向分辨率方位向分辨率点目标的数量点目标位置 .这里设置了 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

33、 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(Lsa

34、r) 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+H2);%实际距离矩阵tau=2*R/C;%回波相对于发射波的延时Dfast=ones(N,1

35、)*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;tr=tm

36、-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 = ft

37、x(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

38、-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);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)en

39、dclose(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;%= %绘图colormap(gray);figure(1)subplot(211);row=tm*C/2-2008;col=sn*V-26;imagesc(row,col,255-Gr); % 距离向压缩 .

40、未校正距离徙动的图像 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(&#

41、39; 距离向压缩 . 校正距离徙动后的图像 '), 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;%

42、光速%雷达参数Fc=1e9;%载频1GHzlambda=C/Fc;%波长%目标区域参数Xmin=0;%目标区域方位向范围 Xmin,XmaxXmax=50;Yc=10000;%成像区域中线Y0=500;%目标区域距离向范围 Yc-Y0,Yc+Y0%成像宽度为 2*Y0%轨道参数V=100;%SAR的运动速度 100 m/sH=5000;%高度 5000 mR0=sqrt(Yc2+H2);%最短距离%天线参数D=4;%方位向天线长度Lsar=lambda*R0/D;%SAR 合成孔径长度 . 合成孔径雷达成像算法与实现P.100Tsar=Lsar/V;%SAR照射时间%慢时间域参数Ka=-2*V

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

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

45、(Rmax-Rmin)/C/dt+Tr/dt);% 快时域的采样数量 Nfast=2nextpow2(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 % 一个矩形以及矩形的中心距离向分辨率方位向分辨率点

46、目标的数量点目标位置 .这里设置了 5个点目标 .构成Xmin+50*DX,Yc-50*DY,1Xmin+25*DX,Yc,1Xmin,Yc+50*DY,1Xmin+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 do

47、main');disp(PRF/Ba)K=Ntarget;%N=Nslow;%M=Nfast;%T=Ptarget;%Srnm=zeros(N,M);%for k=1:1:K%sigma=T(k,3);%Dslow=sn*V-T(k,1);%R=sqrt(Dslow.2+T(k,2)2+H2);tau=2*R/C;%生成回波信号Dfast=ones(N,1)*tm-tau'*ones(1,M);%disp('Sampling Number in slow-time domain');disp(Nslow) disp('Range Resolution&#

48、39;);disp(DY) disp('Cross-range Resolution');disp(DX) disp('SAR integration length');disp(Lsar) disp('Position of targets');disp(Ptarget) %=%(t-tau).阵 .ones(N,1) 和 ones(1,M) 都是为了将其扩展为矩阵phase=pi*Kr*Dfast.2-(4*pi/lambda)*(R'*ones(1,M);%目标数目 慢时域的采样数 快时域的采样数 目标矩阵 生成零矩阵存储回波信号

49、总共 K个目标 得到目标的反射系数 方位向距离 . 投影到方位向的距离 实际距离矩阵 回波相对于发射波的延时 其实就是时间矩相位 . 公式参见 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=ift

50、y(fty(Srnm).*(ones(N,1)*conj(fty(Refr);Gr=abs(Sr);%开始进行距离弯曲补偿正侧视没有距离走动项主要是因为斜距的变化引起回波包络的徙动%补偿方法：最近邻域插值法 .具体为：先变换到距离多普勒域 . 分别对单个像素点计 算出距离徙动量 . 得到距离徙动量与距离分辨率的比值 .%该比值可能为小数 . 按照四舍五入的方法近似为整数 . 而后在该像素点上减去徙动量 %方位向做 fft 处理 再在频域做距离弯曲补偿Sa_RD = ftx(Sr); % 方位向FFT 变为距离多普域进行距离弯曲校正 %距离徙动运算 , 由于是正侧视 .fdc=0, 只需要进行距

51、离弯曲补偿。 Kp=1; % 计算或者预设预滤波比 %第二种方法进行截断 sinc 插值进行距离徒动校正 h = waitbar(0,'Sinc插值 ');P=4;%4点 sinc 插值RMCmaxtix = zeros(N,M); for n=1:Nfor m=P:Mdelta_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_

52、RMC = RMC-round(RMC); % for i = -P/2:P/2-1if m+RMC+i>M %RMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,M)*sinc(pi*(-i+RMC); 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(Sa);toc;%= %绘图colormap(gray);figure(1)subplot(211);row=tm*C/2-2008;col=sn*V-26;imagesc(row,col,255-Gr); %