合成孔径雷达成像——CS算法实现

CS（Chirp Scaling）算法通过相位相乘的形式避免了插值操作，相比于RD算法极大的提高了SAR的成像效率。

仿真代码如下

close all; clear; clc;
c=3e8;%光速
%% 仿真参数
theta = 45;                            %下视角
thetac = 0;                            %斜视角
pd = 3;                                %地距分辨力为3mxy均为3m
pr0 = 3;                               %方位向分辨力
pr = pd*cosd(theta);                   %距离分辨力c/2B=>B>c/2pr
h = 40e3;                              %飞行高度h=10e3;
gr = h*tand(theta);                    %地距
Rc = h/cosd(theta)/cosd(thetac);       %中心穿越时刻距离（距离场景中心距离，对于波束内部的其他点有所差异）
% widthr = 1000;                         %测绘带宽为1000m,y方向
% widtha = 1000;                         %方位向带宽x方向1000m
vs = 7500;
vr = vs/1.06;                          %卫星比等效速度大6% 参考第90页vr=250;%雷达有效速度
vg = vr/1.06;
vr=vs;


% Rc=30e3;%景中心斜距

Tp=10e-6;%脉冲宽度
B=150e6;%信号带宽
Kr=B/Tp;%距离脉冲调频率
fs=1.2*B;%距离采样率
fc=10e9;%载波频率
lambda=c/fc;%波长
% Ka=2*vr^2/lambda/Rc;%方位向调频率


la=2*pr0;%天线真实孔径
fdwidth=0.886*2*vr/la;%多普勒带宽
prf = 5e3;                                   %脉冲重复频率
im=sqrt(-1);%虚数单位
thetac=0;%斜视角度，正侧视
fd0=2*vr*sin(thetac*pi/180)/lambda;%参考中心频率，取多普勒中心频率f_etac 4.34

ls=0.886*Rc*lambda/la;%合成孔径长度
Ta=ls/vr;%目标照射时间
%% 成像区域[Xc-X0,Xc+X0; Yc-Y0,Yc+Y0]
Xc = sqrt(Rc^2-h^2);
Yc = 0;
Xo = 500;
Yo = 500;
Rmin=sqrt(h^2+(Xc-Xo)^2);%观测场景距飞机的最近距离
Rmax=sqrt(h^2+(Xc+Xo)^2);%观测场景距飞机的最远距离
Ra=ls+2*Yo;%正侧视时雷达在方位向行走距离
centerpos=[gr,0,0];
%% 目标位置

target=[centerpos;gr+50,centerpos(2)+50,0;gr-30,centerpos(2)-20,0;gr+50,centerpos(2),0;gr+150,centerpos(2),0]; %点目标位置(第一个为方位向，第二个为 距离向)


widthr = 1000;widtha=1000;
xylim = [gr-widthr/2,gr+widthr/2;-widtha/2,widtha/2;0,0];%成像区域范围

figure;plot(target(:,1),target(:,2),'r.','Markersize',10);
xlim([xylim(1,1),xylim(1,2)]);ylabel('方位向距离/m');xlabel('距离向距离/m'); ylim([xylim(2,1),xylim(2,2)]);title('目标位置')
%% 生成回波
tm=-Ra/vr/2:1/prf:Ra/vr/2-1/prf;%慢时间轴
tk=2*Rmin/c-Tp/2:1/fs:2*Rmax/c+Tp/2-1/fs;%快时间轴
r=sqrt((tk*c/2).^2-h^2);%距离向横坐标
x=r;
nm=length(tm);%方位向采样点数
nk=length(tk);%距离向采样点数
echoall=zeros(nm,nk);%回波
y=vr*tm;%飞机的位置(注意：慢时间轴不同时表达式不同)
Rtm=zeros(size(target,1),nm);
A0=1;%幅度

bar = waitbar(0,"回波生成...");
radar = [zeros(length(tm),1),y',h*ones(length(tm),1)];      %卫星坐标
% Rtm1 = sqrt(sum((target(1,:)-[0,y(280),h]).^2)); plot(rectpuls(tk-2*Rtm1/c-Tp/2,Tp));
for k=1:size(target,1)  %目标数
    echo=[];
        Rtm = sqrt((target(k,1).*ones(length(tm),1)-radar(:,1)).^2+...
            (target(k,2).*ones(length(tm),1)-radar(:,2)).^2+...
            (target(k,3).*ones(length(tm),1)-radar(:,3)).^2)*ones(1,length(tk));
        echo(:,:,k)=(abs(target(k,2)-y')*ones(1,length(tk))<ls/2).*(abs((ones(length(tm),1)*tk-2.*Rtm/c)/2)<=Tp/2).*...
                    exp(1j*2*pi*fc.*(-2*Rtm./c)+1j*pi*Kr*(tk-2.*Rtm./c).^2);  %回波模型(停走模型)
    waitbar(k/size(target,1),bar)
    echoall = echo(:,:,k) + echoall;
end
close(bar)
% for i=1:size(target,1)

tk_mat = ones(nm,1)*tk;    %距离时间轴矩阵
ta_mtx = tm.'*ones(1,nk);  %方位时间轴矩阵
ft_mat = ones(nm,1)*linspace(-fs/2,fs/2,nk); %距离频率轴矩阵
fdc=2*vr*sin(thetac*pi/180)/lambda;%参考中心频率，取多普勒中心频率f_etac 4.34
fd1=fdc+linspace(-prf/2,prf/2,nm).'*ones(1,nk);%方位频率轴矩阵
f=linspace(-fs/2,fs/2,nk);
HFk = rectpuls(f,Kr*Tp).*exp(-1j*pi*(f+0).^2./Kr); %PSP计算已经下变频不需要-fc,若以及移动Tp/2要注意移动+B/2(移动Tp/2)
fd = fdc+linspace(-prf/2,prf/2,nm);

figure,imagesc(tk,tm,abs(echoall)) ;title('原始回波幅度')          

figure,imagesc(tk,tm,angle(echoall)) ;title('原始回波相位')         
caxis([0 1]);axis xy                         

%% 方位向傅里叶变换
echoallrd=zeros(nm,nk);
% for i=1:nk
%    echoallrd(:,i)= fftshift(fft((echoall(:,i)),nm));%转换到多普勒域
% end
   echoallrd= fftshift(fft(echoall,nm,1),1);%转换到多普勒域
figure;imagesc(x,fd,abs(echoallrd));title('原始回波幅度距离多普勒域');xlim([xylim(1,1),xylim(1,2)]); 

% figure;
% imagesc(abs(echoallrd));
%% Chirp Scaling 操作

Dfd=sqrt(1-(c*fd1).^2/(4*vr^2*fc^2));%徙动参数7.17
Km=Kr./(1-(Kr*c*Rc*fd1.^2)./(2*vr^2*fc^3*Dfd.^3));%改变后的距离向调频率7.18
Dfdc=sqrt(1-(c*fdc).^2/(4*vr^2*fc^2));%7.18
% tao_pie=tk_mat-2*Rc./(c*Dfd);%新的距离时间7.27
srcm=exp(im*pi*Km.*(Dfdc./Dfd-1).*(tk_mat-2*Rc./(c*Dfd)).^2);%变标方程7.30
echorcm1rd=echoallrd.*srcm;%变标后
%% 距离向傅里叶变换
echorcm1ff=zeros(nm,nk);

echorcm1ff = fftshift(fft(echorcm1rd,[],2),2);%二维频域
%% 距离压缩，SRC，一致RCMC
f_tao=ones(nm,1)*linspace(-fs/2,fs/2,nk);
echorcm2ff=zeros(nm,nk);
echorcm2ff=echorcm1ff.*exp(im*pi.*Dfd./(Dfdc.*Km).*ft_mat.^2).*...
            exp(im*4*pi/c*(1./Dfd-1./Dfdc)*Rc.*ft_mat);%7.32
        
%% 距离向傅里叶逆变换
echorcm2=zeros(nm,nk);


echorcm2=ifft(fftshift(echorcm2ff,2),[],2);%多普勒域

%% 方位向匹配滤波以及附加相位校正
echorcmc1=zeros(nm,nk);
echorcmc1= echorcm2.*exp(im*4*pi*Rc*fc*Dfd/c).*...
                  exp(im*4*pi*Km/c^2.*(1-Dfd./Dfdc).*(Rc./...
                  Dfd-Rc./Dfdc).^2);%方位压缩以及附加相位补偿7.32

%% SAR成像结果 
res_image=zeros(nm,nk);
% for i = 1:nk
%     res_image(:,i) = ifft(echorcmc1(:,i)); %变回SAR图像域   
% end
res_image  = ifft(echorcmc1,[],1); %变回SAR图像域
figure;[R,Y] = meshgrid(r,y);
mesh(R,Y,abs(( res_image)));view(0,90);xlim([xylim(1,1),xylim(1,2)]);
ylabel('方位向距离/m');xlabel('距离向距离/m'); ylim([xylim(2,1),xylim(2,2)]);title('成像结果')

参考资料《合成孔径雷达成像算法与实现》——lan G. Cumming, et al.