SAR chirp scaling(CSA)算法仿真_sar成像csa算法、

合成孔径雷达，雷达成像技术（保铮院士）SAR chirp scaling(CSA)算法仿真
参考书为 雷达成像技术（保铮等）
代码有注释，非常容易理解
最终成像为绝对坐标

%%          CS 算法
%           正侧视
%%================================================================
clear;clc;close all;
%% 基本参数
c=3e8;
Fc=1e9;                                  %载波
lambda=c/Fc;
Xmin=-50;                                  %[Xmin,Xmax]方位向
Xmax=50;
Yc=10000;                                %场景中心线
Y0=500;                                  %[Yc-Y0,Yc+Y0]距离向
%2*Y0
V=100;                                   %SAR velosity 100 m/s 正侧视
H=5000;                                  %height 5000 m
R0=sqrt(Yc^2+H^2);                       %目标离飞机的最近距离
D=4;                                     %方位向天线长度
Lsar=lambda*R0/D;                        %合成孔径长度
Tsar=Lsar/V;                             %合成孔径时间
%% 方位向参数
Ka=-2*V^2/lambda/R0;                      %方位向多普勒斜率
Ba=abs(Ka*Tsar);                          %方位向多普勒调频范围即多普勒带宽
PRF=Ba;                                   %方位向发射序列频率就等于方位向多普勒频带宽度
PRT=1/PRF;                                %计算方位向发射序列间隔
ds=PRT;                                   %脉冲重复时间
Nslow=ceil((Xmax-Xmin+Lsar)/V/ds);        %慢时间采样点数，但是考虑到边缘，须多加一个合成孔径的采样时间（方位向变换时间/脉冲重复时间）
Nslow=2^nextpow2(Nslow);                  %for fft
tm=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);     %产生方位向信号离散序列，以时间为量纲
PRT=(Xmax-Xmin+Lsar)/V/Nslow;            %方位向单位间隔
PRF=1/PRT;                               %采样频率
ds=PRT;                                  %采样间隔（时间域）
%% 距离向参数
Tr=5e-6;
B=30e6;
miu=B/Tr;                                %调频斜率，正扫频
Fsr=3*B;                                 %采样频率
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);                 %for fft

t=linspace(2*Rmin/c,2*Rmax/c+Tr,Nfast);  %产生距离向信号离散序列，以时间为量纲
dt=(2*Rmax/c+Tr-2*Rmin/c)/Nfast;         %距离向采样间隔
Fsr=1/dt;                                %距离向采样频率

% 分辨率
DY=c/2/B;                                %距离向分辨率
DX=D/2;                                  %方位向分辨率


%% 这几行为了figure 的坐标 对算法没有影响
fa=linspace(-PRF/2,+PRF/2,Nslow);        %方位向
fr=linspace(-Fsr/2,+Fsr/2,Nfast);        %距离向
r=t/2*c;                                 %距离RB  t=linspace(2*Rmin/c,2*Rmax/c+Tr,Nfast);
real_r = sqrt(r.^2-H^2);                 %r是斜距，转换成水平距离  相当于 数据录取平面 到水平面的转化
%% 目标位置
Num_target=3;                            %三个点目标
%format [x, y, reflectivity]
Ptarget=[20,    Yc-Y0/2,    0.1
          0,    Yc,         0.2
          -20,  Yc+Y0/2,    0.3];
%% 回波 
Srnm=zeros(Nslow,Nfast);

for k=1:Num_target
    R=sqrt((V*tm-Ptarget(k,1)).^2+Ptarget(k,2)^2+H^2);%1*512
    tau=2*R/c;                      
    sigma=Ptarget(k,3);
    Dfast=ones(Nslow,1)*t-tau'*ones(1,Nfast);
    phase=pi*miu*Dfast.^2-(4*pi/lambda)*(R'*ones(1,Nfast));
    Srnm=Srnm+sigma*((0<Dfast&Dfast<Tr)).*exp(1j*phase);
end
%% 方位维FFT
Srnm_xfft=zeros(Nslow,Nfast);
for k = 1:Nfast
    Srnm_xfft(:,k) = fftshift(fft(Srnm(:,k))); 
end

figure(1);
set(gcf,'position',[150,250,1200,400]);
subplot(121)
mesh(real_r,fa,abs(Srnm_xfft));axis tight;
title("距离多普勒域");xlabel("距离/m");ylabel("方位/Hz");
subplot(122)
imagesc(real_r,fa,abs(Srnm_xfft));axis tight;
title("距离多普勒域");xlabel("距离/m");ylabel("方位/Hz");
%% CS 
fa=linspace(-PRF/2,PRF/2,Nslow);             %方位向
fr=linspace(-Fsr/2,Fsr/2,Nfast);            %距离向
r=t/2*c;               %距离RB  t=linspace(2*Rmin/c,2*Rmax/c+Tr,Nfast);
real_r = sqrt(r.^2-H^2);%r是斜距，转换成水平距离  相当于 数据录取平面 到水平面的转化

sin_theta=lambda*fa/2/V;
cos_theta= sqrt(1-(sin_theta).^2);
afa=1./cos_theta-1;
Rs=sqrt(H^2+Yc^2);

Ks=miu./(1-miu*Rs*2*lambda/c/c.*(sin_theta.^2)./(cos_theta.^3));
H1=exp(1j*pi*(Ks.'*ones(1,Nfast)).*(afa.'*ones(1,Nfast)).*((ones(Nslow,1)*t-2*Rs/c*(afa+1).'*ones(1,Nfast))).^2);
Srnm_cs=Srnm_xfft.*H1;
%% 距离FFT
Srnm_yfft=zeros(Nslow,Nfast);
for n = 1:Nslow
    Srnm_yfft(n,:) = fftshift(fft(Srnm_cs(n,:)));    
end

figure(2);
set(gcf,'position',[150,250,1200,400]);
subplot(121)
mesh(fr,fa,abs(Srnm_yfft));axis tight;
title("距离频域-方位频域");xlabel("距离/Hz");ylabel("方位/Hz");
subplot(122)
imagesc(fr,fa,abs(Srnm_yfft));axis tight;
title("距离频域-方位频域");xlabel("距离/Hz");ylabel("方位/Hz");
%% 一致RCMC和距离压缩
H2=exp(1j*pi*(ones(Nslow,1)*fr.^2)./((Ks.*(1+afa))'*ones(1,Nfast)) ).*exp( 1j*4*pi/c*(ones(Nslow,1)*fr)*Rs.*(afa'*ones(1,Nfast)) );%%------（7.32式）
Srnm_rmc=Srnm_yfft.*H2;                 %距离徙动校正与距离压缩
%% 距离IFFT
Srnm_yifft=zeros(Nslow,Nfast);

for n = 1:Nslow
    Srnm_yifft(n,:) = ifft(fftshift(Srnm_rmc(n,:)));    
end


figure(3);%经过距离脉压之后
set(gcf,'position',[150,250,1200,400]);
subplot(121)
mesh(real_r,fa,abs(Srnm_yifft));axis tight;
title("距离多普勒域");xlabel("距离/m");ylabel("方位/Hz");
subplot(122)
imagesc(real_r,fa,abs(Srnm_yifft));axis tight;
title("距离多普勒域");xlabel("距离/m");ylabel("方位/Hz");
%% 方位压缩与相位补偿
r_sub=ones(Nslow,1)*(r-Rs).^2;
phase_cor=(4*pi/c^2*Ks.*(1+afa).*afa)'*ones(1,Nfast).*r_sub;
H3=exp(1j*4*pi/lambda*ones(Nslow,1)*r.*( sqrt(1-((lambda*fa/2/V).^2)'*ones(1,Nfast)) )-1j*phase_cor-j*2*pi*c*ones(Nslow,1)*t/lambda   ); %%-----（7.34式）
Srnm_cor=Srnm_yifft.*H3;                %方位压缩与相位补偿
%% 方位IFFT
f_xy=zeros(Nslow,Nfast);
%方位IFFT
for n = 1:Nfast
    f_xy(:,n) = ifft(Srnm_cor(:,n));    
end

%% 成像
figure(4);
set(gcf,'position',[150,250,1200,400]);
subplot(121)
mesh(real_r,tm*V,abs(f_xy));axis tight;
title("最终成像");xlabel("距离/m");ylabel("方位/m");
subplot(122)
imagesc(real_r,tm*V,abs(f_xy));axis tight;
title("最终成像");xlabel("距离/m");ylabel("方位/m");
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155