wk算法的仿真

这是我的第一篇博客，写的初衷是因为当时学做的时候找了很多相关资料，但没有直接能参考的代码，有很大阻力，所以把我完成的这一点东西放这里和大家交流讨论，也希望我能帮到需要的人做入门。

我入门的时候找到一片RDA算法仿真的文章，大佬的代码很工整，井井有条，给我入门雷达点目标仿真提供了不少帮助，本文中仿真信号生成部分也是几乎完全照抄了大佬的代码逻辑，这部分不懂的可以看看原文代码，链接如下合成孔径雷达基于RDA算法的点目标仿真_rda算法点目标_风雨中的小小鸟的博客-CSDN博客

wk算法有时候又被称f-k偏移算法，主要是在频率域做信号处理完成的雷达成像，原理可以看《Digital Processing of Synthetic Aperture Radar Data Algorithms and Implementation》，翻译本是由电子工业出版社出版的《合成孔径雷达成像算法与实现》，工作流程就是

 参考函数相乘后信号特点：参考距离处完全聚焦，其他点离参考距离越远，聚焦效果越差

Stolt插值后信号特点：无论离参考距离处远近的点都能完全聚焦

上传的代码我做过一些修改和批注，应该能解决大部分的问题。当初最困扰我的是87行的那个线性补偿，是找了很多资料做了方位向和距离向补偿后，调试时把方位向补偿去掉突然解决的，真要解释清楚可能有点问题，我认为的方向是我做傅里叶变换时有时shift有时没shift回来造成的相位不对齐，所以需要做这样一个补偿。但是我目前没有解释这个问题的压力了，这个部分留给你们去探讨研究啦（开摆）！

另一个大家可能会有疑问的部分是Stolt插值里面为什么计算这么多频率，在这里解释一下：

f_tao2all：所有可能映射到的新频率fτ‘

f_tao2：实际存在的新频率fτ‘

fr_new：由新频率fτ’映射回去得对应的距离向频率，用于插值函数中完成Stolt插值（这里利用的是matlab提供的interp1（）函数，这里不是说interp1（）自带了Stolt插值的功能，这个只是个matlab自带的插值函数，千万不要以为Stolt插值由这个函数直接完成了，要好好关注之前的那个频率映射，那才是Stolt插值的核心，interp1（）只是用来实现的）

由上一步描述，f_tao2实际舍去了一些f_tao2all的点，为什么需要这么做呢？

举个例子，假如方位向频率和距离向频率都是1~100，公式不好打且计算量大，就随便以

fτ’^2=fτ^2+fη^2为例：fτ‘的全部范围应该在[,100]，但是反映射回去时，当fτ’为100,fη为1时，fτ‘会存在的原域并不存在的解，所以需要舍去一部分。

（上面可能讲的有点乱，多理理频率映射的关系和各个频率的概念就能懂了，感谢大家担待）

有人提醒了我：

程序里的iftx和ifty是自定义的函数，比如iftx（）就是在运行目录新建一个.m文件并命名为iftx，代码为下面的两行，这样可以省去fftshift和ifftshift占用的地方，这一点也是我刚学习的时候从某篇文章得到的思路，具体链接找不到了，感谢这位作者。

function s=iftx(fs)
s=ifftshift(ifft(fftshift(fs)));

完整的仿真代码如下，跑出的效果我认为还行吧

clear ;close all;
%% 仿真雷达参数
c = 3e8;%光速
Rc = 2e4; %景中心斜距
Vr = 100; %雷达等效速度
Tr = 5e-6; %发射脉冲时宽
Kr = 6e12; %(Hz/s)距离向调频率
f0 = 1e9; %雷达工作频率
lamda = c/f0; %雷达工作段波长，计算调频率需要
H = 0; %雷达轨道距地面高度
squint_c = 9/180*pi;%波束斜视角
fc = 2*Vr*sin(squint_c)/lamda;%多普勒中心频率

%% 仿真目标参数
Lmax=50;Lmin =0;%波束照射方位向范围
Wid_azi = Lmax - Lmin;%波束方位向宽度
LD_c = 10000;LD = 500;%距离向中心点LD_c；范围在LD_c±LD
Ymin = LD_c - LD;Ymax = LD_c + LD;%用于距离向分割
R0 = sqrt(LD_c^2 + H^2);%景中心距离，用作参考函数相乘部分的参考点
La = 4;%方位向天线长度
Lsar=lamda*R0/La; %合成孔径长度
Tsar=Lsar/Vr;%合成孔径内照射时间（某个目标点受照射的时间）

Rmin = sqrt(H^2 + Ymin^2);%探测范围内最短距离点
Rmax = sqrt(H^2 + Ymax^2 + (Lsar/2)^2);%探测范围内最远距离点

K=3;%仿真点数
x1 = 0;y1 = 10000; x2 = 0; y2 = 9800; x3 = 40; y3 = 10300;
target=[x1,y1,1;
        x2,y2,1;
        x3,y3,1];
    
%回波信号仿真
Ka = -2*(Vr^2)/lamda/R0;
N = (Wid_azi+Lsar)/Vr*2*abs(Ka*Tsar);       %方位向采样数
M = 2*(Rmax - Rmin)/c*2*Kr*Tr + Tr*2*Kr*Tr; %距离向采样数

N = 2^nextpow2(N);
M = 2^nextpow2(M);%将N和M化为2的幂次方便于fft

T_azimuth = linspace((Lmin-Lsar/2)/Vr,(Lmax+Lsar/2)/Vr,N);%方位向时间
T_tao = linspace(2*Rmin/c - 0.5*Tr,2*Rmax/c + 0.5*Tr,M);  %接收器开启时间段（时间起点可以自己定，总长为2*（Rmax-Rmin）/c + Tr即可）
sig_re = zeros(N,M);
for k = 1:1:K
    A_refle = target(k,3);             %反射系数（对应目标点电介质常数）
    L_bet = T_azimuth*Vr - target(k,1);%接收端与目标的方位向距离
    r0 = sqrt(H^2 + target(k,2)^2);    %瞬时水平平面距离
    R_inta = sqrt(r0^2 + L_bet.^2);     %瞬时斜距
    T_behind = 2*R_inta/c;              %被减去的2R/c延时
    t = ones(N,1)*T_tao - T_behind'*ones(1,M);
    phase = pi*(Kr*(t.^2) - 4*f0/c*R_inta'*ones(1,M));
    sig_re = sig_re + A_refle*exp(1i*phase).*(0<t&t<Tr).*((abs(L_bet)<Lsar/2)'*ones(1,M));
    %%仿真信号的生成依据书<合成孔径雷达成像算法与实现>page105的（5.2）式完成
    %%英文原著为《Digital Processing of Synthetic Aperture Radar Data Algorithms and Implementation》
    %%上述页码是由电子工业出版社出版的翻译本，名字《合成孔径雷达成像算法与实现》
end

%% 点目标仿真图像
figure(1) 
subplot(211) 
imagesc(abs(sig_re));title('SAR 3point') 
subplot(212) 
imagesc(angle(sig_re)) 

row=T_tao*c/2;  %% 行相当于距离（不是y）
col=T_azimuth*Vr;    %% 列

%% 二维傅里叶变换
Sig_re = fftshift(fft(ifftshift(sig_re,1),N,1),1);
Sig_re = fftshift(fft(ifftshift(Sig_re,2),M,2),2);
%% 参考函数相乘
%频域分点
%注：此处命名有不规范，因为是一开始照着RDA那篇文章自己写的，概念有模糊，但是懒得
%%%改了，一个是方位向时间和频率，不应该用PRT命名；一个是距离向时间和频率
PRT_amu = (Lmax - Lmin + Lsar)/Vr/N;
PRF_amu = 1/PRT_amu;
f_inta = linspace(-PRF_amu/2,PRF_amu/2,N) ;%+ fc;

PRT_dis = (2*Rmax/c - 2*Rmin/c + Tr)/M;
PRF_dis = 1/PRT_dis;
f_tao = linspace(-PRF_dis/2,PRF_dis/2,M);

%参考函数构造
Rref = 9500;%取参考距离点为9500处
phase_ref = 4*pi*Rref/c*sqrt(ones(N,1)*(f0+f_tao).^2 - c^2*(f_inta'*ones(1,M)).^2/4/(Vr^2)) + pi*ones(N,1)*(f_tao).^2/Kr;
%参考函数的构造推导详情看书<合成孔径雷达成像算法与实现>中第八章wk算法，结果在page105的（8.3）式
phase_ref = phase_ref - 2*pi*f_tao*2*R0/c;%这里出来的有点偶然，应该是一个线性相位延迟，R0是景中心频率
H_ref = exp(1i*phase_ref);

%参考函数相乘
Sig_RFM = Sig_re.*H_ref;
WinForR = ones(N,1)*kaiser(M,2.5)';%凯瑟窗
Sig_RFM = Sig_RFM.*WinForR;

%成像看效果，方便debug
sig_ref = iftx(ifty(Sig_RFM));

figure(2)
colormap(gray)
imagesc(row,col,255-abs(sig_ref));
title('一致压缩'),xlabel('距离向'),ylabel('方位向');

figure(3)
contour(row,col,abs(sig_ref));
title('一致压缩等高线图'),xlabel('距离向'),ylabel('方位向');

%% stolt插值
Sstolt_phase_NaN = zeros(N,M);Sstolt_abs_NaN = zeros(N,M);Sig_stolt = zeros(N,M); 
%依据page206的（8.5）式计算出f_tao'用以替代f_tao从而完成距离徙动的修正；
f_tao2all = sqrt((f0+f_tao).^2 - (c*f_inta'/2/Vr).^2) - f0;%%f_tao映射到f_tao'

f2_min = max(min(transpose(f_tao2all)));
f2_max = min(max(transpose(f_tao2all)));
f_tao2 = linspace(f2_min,f2_max,M) ;%这三行选f_tao在范围内的f_tao'

fr_new = sqrt((f0+f_tao2).^2 + (c*f_inta'/2/Vr).^2) - f0;%%f_tao'映射到f_tao

for i1 = 1:1:N
       Sstolt_phase_NaN(i1,:) = interp1(f_tao(:),unwrap(angle(Sig_RFM(i1,:))),fr_new(i1,:),'spline',0); 
       Sstolt_abs_NaN(i1,:) = interp1(f_tao(:),abs(Sig_RFM(i1,:)),fr_new(i1,:),'spline',0); 
       %注：这里插值也可以直接对复数信号Sig_RFM做插值，这里是因为debug的时候不出结果，所以对幅角和幅度分开插值
end
Sig_stoltNaN = Sstolt_abs_NaN.*exp(1i*Sstolt_phase_NaN);
%下述步骤用于清除一部分不在理论幅度范围内的信号，用补0处理
num_find = find(~isnan(Sig_stoltNaN));
Sig_stolt(num_find) = Sig_stoltNaN(num_find);
sig_stolt = iftx(ifty(Sig_stolt));

figure(4)
contour(row,col,abs(sig_stolt));
title('stolt插值后等高线图'),xlabel('距离向'),ylabel('方位向');

figure(5)
colormap(gray)
imagesc(row,col,255-abs(sig_stolt));
title('stolt插值后幅度图'),xlabel('距离向'),ylabel('方位向');
deta = abs(sig_stolt - sig_ref);

谢谢大家看到这里，有问题或者补充欢迎讨论呀