SAR成像波数域WK成像算法

系列：

（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博客

WK算法

区别于RD,CS:在二维频域通过一种特别的方式来校正距离方位耦合与距离时间和方位频率的依赖关系。使用精确的双曲线模型，而非二阶近似，具有对宽孔径或者大斜视角数据的处理能力。

(1)处理流程

• 通过二维FFT将数据变到二维频域
• 参考函数相乘，这是WKA的第一个关键聚焦步骤。参考函数根据选定的距离(通常为测绘中心)来计算，它补偿了该距离处包括距离向频率调制、距离徙动、距离方位耦合和方位向频率调制在内的各种相位。相乘后参考点处得到了完全聚焦，而非参考点得到了部分聚焦。
• Stolt插值，这是WKA的第二个关键步骤。它在距离频域用插值方式完成其它目标的聚焦。
• 通过二维IFFT将信号变回到时域

  (3)实现细节：

①参考函数相乘

二维频域中未经压缩的基带信号为：

   

对于距离为的目标而言，相位为：

             

注意由于多普勒中心频率依赖于雷达波长，因此会随着距离频率变化，这意味着信号二维频谱在方位向总是扭曲的。

因此RFM滤波器的相位为：

          

经过处理后，二维频域中的残余相位近似为：

②Stolt插值：完成残余RCMC、残余SRC以及残余方位压缩

  变量代换：

                                     

③二维傅里叶逆变换得到压缩后的图像

(4)仿真

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]
%以合成孔径中心为原点，距离向为x轴，方位向为y轴
Xc = sqrt(R_etac^2-H^2);
Yc = 0;
Xo = 1000;%应该看不了这么大区域，设置大主要为了看非参考处的散焦
Yo =300;
Rmin=sqrt(H^2+(Xc-Xo)^2);%观测场景距飞机的最近距离
Rmax=sqrt(H^2+(Xc+Xo)^2);%观测场景距飞机的最远距离
Ra=Ls+2*Yo;%正侧视时雷达在方位向行走距离
R_ref=R_etac;%参考距离
%% 目标位置
target = [Xc,Yc;
          Xc-800,Yc+100;
          Xc-800,Yc-200;
          Xc+500,Yc-200;
          Xc+800,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=sqrt((tao*c/2).^2-H^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_AFRF=fftshift(fft2(signal_receive));
%% 参考函数相乘（一致压缩）
f_tao=linspace(-Fr/2,Fr/2,Nr);%距离频率轴
f_tao_mtx = ones(Na,1)*f_tao;%距离频率轴矩阵
f_eta_ref=2*Vr*sin(sita*pi/180)/lamda;%参考中心频率，取多普勒中心频率f_etac 4.34
f_eta=linspace(-Fa/2,Fa/2,Na);%方位频率轴
f_eta_mtx=f_eta_ref+f_eta.'*ones(1,Nr);%方位频率轴矩阵
Signal_compress=Signal_AFRF.*exp(im*4*pi*R_ref/c...
    *sqrt((f0+f_tao_mtx).^2-c^2*f_eta_mtx.^2/(4*Vr^2))+im*pi*f_tao_mtx.^2/Kr);
figure;
mesh(abs((ifft2(Signal_compress))));
view(0,90);
xlabel('距离向');
ylabel('方位向');
zlabel('幅度'); title('未进行stolt插值的成像结果（非参考距离处散焦明显）');
%% 距离频域进行stolt插值
Signal_stolt=zeros(Na,Nr);
P=6;%截断sinc插值的核函数的点数
f_tao_pie=sqrt((f0+f_tao_mtx).^2-(c*f_eta_mtx).^2/(4*Vr^2))-f0;
delta_n=((f_tao_pie-f_tao_mtx)/Fr*Nr/2);%除以2中心就不散焦了，但我感觉应该除以1
for m=1:Na
     for n=1:Nr
        for i=-P/2:1:P/2-1
            if n+i<=0
               Signal_stolt(m,n)= Signal_stolt(m,n)+Signal_compress(m,n)*sinc(delta_n(m,n)-i);
            else
                if n+i<=Nr
                   Signal_stolt(m,n)= Signal_stolt(m,n)+Signal_compress(m,n+i)*sinc(delta_n(m,n)-i);
                else
                   Signal_stolt(m,n)= Signal_stolt(m,n)+Signal_compress(m,n)*sinc(delta_n(m,n)-i);
                end
             end
        end
     end
end
%% 参考距离平移
Signal_translation=Signal_stolt.*exp(-im*4*pi*R_ref/c*f_tao_mtx);
%% 点目标成像结果
Signal_ATRT=ifft2((Signal_translation));
figure;
mesh(r,y,abs(Signal_ATRT));
view(0,90);xlim([27400 29100]);
xlabel('距离向');
ylabel('方位向');
zlabel('幅度'); title('点目标成像结果');%中心散焦了，证明程序有问题（大概率是插值），有懂哥解决了，麻烦不吝赐教

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