Chirp Scaling SAR成像算法(CS算法)

（1）SAR距离多普勒成像（RD）算法_Xc Lbb的博客-CSDN博客

（2）Chirp Scaling SAR成像算法(CS算法)_chirp scaling算法_Xc Lbb的博客-CSDN博客

（3）SAR成像波数域WK成像算法_Xc Lbb的博客-CSDN博客

（4）SAR后向投影(BP)成像算法_Xc Lbb的博客-CSDN博客

Chirp Scaling 算法

(1)原理：通过对Chirp信号进行频率调制，实现对该信号的尺度变换(变标)或平移。基于这种原理，可以通过相位相乘的方法代替距离时域插值来校正随距离变化的RCM，此外，由于需要在频域进行数据处理，CS算法还能解决SRC对方位频率依赖的问题。

由于通过频率调制实现变标或者平移不能太大，否则会引起信号的中心频率以及带宽发生改变。因此一般通过两步实现RCM的校正，第一步校正不同距离门上的RCM，使所有信号具有一致的RCM；第二步则可以在二维频域中通过相位相乘完成RCMC。

(2) Chirp Scaling实现步骤

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

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

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

4．通过与参考函数进行相位相乘，同时完成距离压缩、SRC和一致RCMC。这是第二步相位相乘。

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

6．通过与随距离变化的匹配滤波器进行相位相乘，实现方位压缩。此外，由于步骤2中的Chirp Scaling操作，相位相乘中还需要附加一项相位校正。这是第三步（最后一步）相位相乘。

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

(3)具体处理步骤

  ①雷达原始数据：

解调后单个点目标的基带信号可以表示为：

                                              

其中——距离向快时间轴(光速)

       ——方位向慢时间轴(波束照射区的推移)

    ——多普勒中心频率

   ——距离向包络(矩形窗函数)

  ——方位向包络(sinc平方型函数)

   ——载波频率

 ——距离向调频率

瞬时斜距： (WK算法一般用到二阶近似)

②方位向傅里叶变换

忽略较小的距离包络调制，距离多普勒域中的信号频谱可以表示为：

                                     

徙动参数:

                                  

新的距离向调频率:

                                 

③补余RCMC中的Chirp Scaling操作：

       在满足发射脉冲是线性调频的，雷达等效速度Vr以及改变后的调频率Km不随距离改变，则变标方程是线性调频的，基于以上假设变标方程为：

                

其中新的距离时间:

                                          

将距离多普勒域中的信号与变标方程相乘：

                                 

④利用POSP计算傅里叶积分，得到二维频域中的表达形式：

                                                              

                               

                               

                                

五个指数项的解释如下：

第一个指数项包含方位调制，可近似为随距离变化的关于方位频率的二次函数。方位向处理时补偿掉。

第二个指数项表示变标后的距离调制，其为的二次函数。由于因子，它微弱的依赖于距离和方位。当将其变换回距离时域后，会发现其中的 表征的是变标因子1+α，并且包含了距离方位耦合。

第三个指数项表征的是点目标位置的线性相位，距离压缩后目标峰值将出现在此处。

第四个指数项是一致RCM，可近似为的二次函数。

第五个指数项是一个附加相位，它是距离和方位的函数，在方位处理时进行补偿。

⑤参考函数相乘用于距离压缩、SRC和一致RCMC(参考函数为第二个和第四个指数项的复共轭)

⑥距离向傅里叶逆变换                         

                                

⑦方位压缩及附加相位校正

⑧方位向傅里叶逆变换得到压缩数据

(4)仿真

数据来源于表4.1，典型的机载数据

close all; clear; clc;
%% 仿真参数
%参数来源表4.1，P98
%2023.6.27  lbb
R_etac=30e3;%景中心斜距
H=10e3;%飞行高度
Tr=10e-6;%脉冲宽度
B=100e6;%信号带宽
Kr=B/Tr;%距离脉冲调频率
Fr=1.2*B;%距离采样率
Vr=250;%雷达有效速度
f0=9.4e9;%载波频率
c=3e8;%光速
lamda=c/f0;%波长
Ka=2*Vr^2/lamda/R_etac;%方位向调频率
La=1;%天线真实孔径
Ls=0.886*R_etac*lamda/La;%合成孔径长度
Ta=Ls/Vr;%目标照射时间
Bw_doppler=0.886*2*Vr/La;%多普勒带宽
Fa=600;%方位向采样率
im=sqrt(-1);%虚数单位
sita=0;%斜视角度，正侧视
%% 成像区域[Xc-X0,Xc+X0; Yc-Y0,Yc+Y0]
Xc = sqrt(R_etac^2-H^2);
Yc = 0;
Xo = 500;
Yo =300;
Rmin=sqrt(H^2+(Xc-Xo)^2);%观测场景距飞机的最近距离
Rmax=sqrt(H^2+(Xc+Xo)^2);%观测场景距飞机的最远距离
Ra=Ls+2*Yo;%正侧视时雷达在方位向行走距离
%% 目标位置
target = [Xc,Yc;
          Xc-300,Yc+100;
          Xc-300,Yc-200
          Xc-100,Yc-100];  
%% 生成回波
eta=-Ra/Vr/2:1/Fa:Ra/Vr/2-1/Fa;%慢时间轴
tao=2*Rmin/c-Tr/2:1/Fr:2*Rmax/c+Tr/2-1/Fr;%快时间轴
r=((tao*c/2).^2-H^2).^(1/2);%距离向横坐标
Na=length(eta);%方位向采样点数
Nr=length(tao);%距离向采样点数
signal_receive=zeros(Na,Nr);%回波
y=Vr*eta;%飞机的位置(注意：慢时间轴不同时表达式不同)
R_eta=zeros(size(target,1),Na);
A0=1;%幅度
for i=1:size(target,1)
    R_eta(i,:)=sqrt(target(i,1)^2+(target(i,2)-y).^2+H^2);%瞬时斜距
     for j=1:Na
         signal_receive(j,:)=A0*rectpuls(tao-2*R_eta(i,j)/c,Tr).*(abs(target(i,2)-y(j))<Ls/2).*...
         exp(-im*4*pi*f0*R_eta(i,j)/c).*exp(im*pi*Kr*(tao-2*R_eta(i,j)/c).^2)+signal_receive(j,:);
    end
end
%% 方位向傅里叶变换
Signal_azimuth_FFT=zeros(Na,Nr);
for i=1:Nr
   Signal_azimuth_FFT(:,i)= fftshift(fft((signal_receive(:,i)),Na));%转换到多普勒域
end
figure;
imagesc(abs(Signal_azimuth_FFT));
%% Chirp Scaling 操作
tao_mtx = ones(Na,1)*tao;    %距离时间轴矩阵
eta_mtx = eta.'*ones(1,Nr);  %方位时间轴矩阵
f_tao_mtx = ones(Na,1)*linspace(-Fr/2,Fr/2,Nr); %距离频率轴矩阵
f_eta_ref=2*Vr*sin(sita*pi/180)/lamda;%参考中心频率，取多普勒中心频率f_etac 4.34
f_eta=f_eta_ref+linspace(-Fa/2,Fa/2,Na).'*ones(1,Nr);%方位频率轴矩阵
D_feta_Vr=sqrt(1-(c*f_eta).^2/(4*Vr^2*f0^2));%徙动参数7.17
Km=Kr./(1-(Kr*c*R_etac*f_eta.^2)./(2*Vr^2*f0^3*D_feta_Vr.^3));%改变后的距离向调频率7.18
D_fetaref_Vrref=sqrt(1-(c*f_eta_ref).^2/(4*Vr^2*f0^2));%7.18
tao_pie=tao_mtx-2*R_etac./(c*D_feta_Vr);%新的距离时间7.27
Signal_scaling=exp(im*pi*Km.*(D_fetaref_Vrref./D_feta_Vr-1).*tao_pie.^2);%变标方程7.30
Signal_scaled=Signal_azimuth_FFT.*Signal_scaling;%变标后
%% 距离向傅里叶变换
Signal_AfRf=zeros(Na,Nr);
for i=1:Na
    Signal_AfRf(i,:)=fftshift(fft((Signal_scaled(i,:))));%二维频域
end
%% 距离压缩，SRC，一致RCMC
f_tao=ones(Na,1)*linspace(-Fr/2,Fr/2,Nr);
Signal_RCMC=zeros(Na,Nr);
Signal_RCMC=Signal_AfRf.*exp(im*pi.*D_feta_Vr./(D_fetaref_Vrref.*Km).*f_tao_mtx.^2).*...
            exp(im*4*pi/c*(1./D_feta_Vr-1./D_fetaref_Vrref)*R_etac.*f_tao_mtx);%7.32
%% 距离向傅里叶逆变换
Signal_AFRT=zeros(Na,Nr);
for i=1:Na
     Signal_AFRT(i,:)=ifft(fftshift(Signal_RCMC(i,:)));%多普勒域
end
%% 方位向匹配滤波以及附加相位校正
Signal_azimuth_compensations=zeros(Na,Nr);
Signal_azimuth_compensations= Signal_AFRT.*exp(im*4*pi*R_etac*f0*D_feta_Vr/c).*...
                  exp(im*4*pi*Km/c^2.*(1-D_feta_Vr./D_fetaref_Vrref).*(R_etac./...
                  D_feta_Vr-R_etac./D_feta_Vr).^2);%方位压缩以及附加相位补偿7.32
figure;
[R,Y] = meshgrid(r,y);
mesh(R,Y,abs(( Signal_azimuth_compensations)));view(0,90);xlim([27800 28400]);
%% SAR成像结果 
signal_result=zeros(Na,Nr);
for i = 1:Nr
    signal_result(:,i) = ifft(Signal_azimuth_compensations(:,i)); %变回SAR图像域   
end
figure;
[R,Y] = meshgrid(r,y);
mesh(R,Y,abs(( signal_result)));view(0,90);xlim([27800 28300]);ylim([-300,300]);

仿真结果：