《合成孔径雷达成像算法与实现》仿真代码：点目标仿真（书5.6节）

一、仿真代码及结果

%% 教材5.6节 点目标仿真 by ZZQclearvars;close all;clc;%% 仿真参数Retac = 20E3; %景中心斜距 mVr = 150; %等效雷达速度 m/sTr = 25E-6; %发射脉宽时宽 sKr = 0.25E12; %距离向调频率 Hz/sf0 = 5.3E9; %雷达工作频率 HzBWdop = 80; %多普勒带宽 HzFr = 7.5E6; % 距离向采样率 HzFa = 104; %方位向采样率 HzNaz = 256; %距离线数Nrg = 256; %距离线采样点数 BW = Tr*Kr; %发射信号带宽 Hzc = 3E8; %光速 m/slambda = c/f0; %发射载波波长 mthetasqc = deg2rad(3.5); %波束斜视角etac = -Retac*sin(thetasqc)/Vr;fetac = 2*Vr*sin(thetasqc)/lambda;La = 0.886*2*Vr*cos(thetasqc)/BWdop; % 方位向天线长度，根据公式（4.36）betabw = 0.886*lambda/La; %方位向波束宽度R0 = Retac*cos(thetasqc); %最近斜距%% 回波信号生成trgVec = (round(2*Retac/c*Fr)-Nrg/2:1:round(2*Retac/c*Fr)+Nrg/2-1)/Fr; %快时间序列tazVec = (round(etac*Fa)-Naz/2:1:round(etac*Fa)+Naz/2-1)/Fa; %慢时间序列radarposVec = tazVec*Vr; %雷达方位向位置序列，以零多普勒位置为原点RtargetVec = sqrt((Retac*cos(thetasqc))^2+radarposVec.^2); %雷达和目标间距离序列delayVec = 2*RtargetVec/c; %雷达到目标的延迟序列 winazVec = sinc(0.886*(asin(-radarposVec./RtargetVec)-thetasqc)/betabw).^2;%方位向窗函数序列，winazVec(i)表示tazVec(i)慢时间下的目标回波的天线增益trgMat = repmat(trgVec,Naz,1);RtargetMat = repmat(RtargetVec',1,Nrg);delayMat = repmat(delayVec',1,Nrg);winrgMat = (abs(trgMat-delayMat)/Tr <= 0.5);winazMat = repmat(winazVec',1,Nrg);signalMat = winrgMat.*winazMat.*...exp(-1i*2*pi*f0*delayMat).*exp(1i*pi*Kr*(trgMat-delayMat).^2);%回波信号序列Naz*Nrg，其中signal(i,j)表示慢时间为tazVec(i)，快时间为trgVec(j)的回波信号采样值%% 时域信号figure('Name','时域信号')subplot(1,2,1);imagesc(abs(signalMat));xlabel('距离向（采样点）');ylabel('方位向（采样点）');title('幅度');subplot(1,2,2);imagesc(angle(signalMat));xlabel('距离向（采样点）');ylabel('方位向（采样点）');title('相位');%% 距离多普勒域频谱Taz = Naz/Fa; %孔径时间tazBase = round(etac/Taz)*Taz;idtaz = find(tazVec-tazBase>=0,1);SRDMat = fft(circshift(signalMat,-idtaz+1,1),[],1);% SRDMat = fft(signalMat,[],1);figure('Name','距离多普勒域信号')subplot(1,2,1);imagesc(abs(SRDMat));xlabel('距离向（采样点）');ylabel('方位频率（采样点）');title('幅度');subplot(1,2,2);imagesc(angle(SRDMat));xlabel('距离向（采样点）');ylabel('方位频率（采样点）');title('相位');%% 二维频域频谱Trg = Nrg/Fr; %快时间采样时间trgBase = round(2*Retac/c/Trg)*Trg;idtrg = find(trgVec-trgBase>=0,1);S2dfMat = fft(circshift(SRDMat,-idtrg+1,2),[],2);% S2dfMat = fft(SRDMat,[],2);figure('Name','二维频域信号')subplot(1,2,1);imagesc(abs(S2dfMat));xlabel('距离频率（采样点）');ylabel('方位频率（采样点）');title('幅度');subplot(1,2,2);imagesc(angle(S2dfMat));xlabel('距离频率（采样点）');ylabel('方位频率（采样点）');title('相位');%% 比较仿真二维频谱相位和理论推导二维频谱相位% ----------------------------------------------------% 参照书5.3节的公式(5.29)获得理论上精确的二维频谱相位% ----------------------------------------------------Mamb = floor(fetac/Fa); fazVec = (0:Naz-1)*Fa/Naz+Mamb*Fa;fetacId = 0; %多普勒中心频率在fazVec中的索引值for i = 1:Nazif (fazVec(i)-fetac < -Fa/2)fazVec(i) = fazVec(i)+Fa;elseif (fazVec(i)-fetac >= Fa/2)fazVec(i) = fazVec(i)-Fa;elseif (fazVec(i)>fetac && fazVec(i)-fetac<Fa/Naz)fetacId = i;endend% 根据多普勒中心频率估计绝对多普勒频率frgVec = fftshift((-Nrg/2:Nrg/2-1)*Fr/Nrg);DfetaVrVec = sqrt(1-c^2/(4*Vr^2*f0^2)*fazVec.^2);S2dfphaseMat = -pi/Kr*repmat(frgVec.^2,Naz,1)...-4*pi*R0*f0/c*sqrt(repmat((DfetaVrVec.^2).',1,Nrg)+repmat(2*frgVec/f0+frgVec.^2/(f0^2),Naz,1));S2dfphaseMat = angle(exp(1i*S2dfphaseMat));figure('Name','二维频域相位')subplot(1,2,1);imagesc(angle(S2dfMat));xlabel('距离频率（采样点）');ylabel('方位频率（采样点）');title('仿真相位');subplot(1,2,2);imagesc(S2dfphaseMat);xlabel('距离频率（采样点）');ylabel('方位频率（采样点）');title('理论相位');

斜视角为22.8°时，时域仿真信号：

斜视角为22.8°时，仿真信号距离多普勒域：

斜视角为22.8°时，仿真信号二维频域：

斜视角为22.8°时，仿真信号二维频域和理论二维频域比较：

​ 通过仿真信号二维频域和理论二维频域的比较可以基本证实仿真的正确性。

二、仿真参数设置及实现思路

​ 仿真参数按书中表5.2设置，其中成像场景按书中图4.1进行，为了简化分析，画出斜距平面（目标和传感器路径形成的平面）的几何关系如下图

图1 斜距平面成像几何关系

三、注意事项

以下内容为个人观点，如有错误敬请谅解，也麻烦大家批评指正！

1.时间轴设置

​ 为了使仿真信号与书中理论推导保持一致，慢时间轴零点应设置为零多普勒时刻（雷达在上图中 P 2 P_2 P2​位置的时刻），因为距离方程为

R ( η ) = R 0 2 + V r 2 η 2 \begin{equation} \tag {3.1} \begin{aligned} R(\eta)=\sqrt{R_0^2+V_r^2\eta^2} \end{aligned} \end{equation} R(η)=R02​+Vr2​η2 ​​​(3.1)​

式(3.1)中 R 0 R_0 R0​为上图中 P 2 T 1 P_2T_1 P2​T1​的长度， V r η V_r\eta Vr​η等于在 η \eta η时刻雷达位置 P η P_\eta Pη​到 P 2 P_2 P2​的距离。

​ 仿真代码中的27、28行设置了快时间轴和慢时间轴。

2.二维频域的准确获取

​ 由于FFT和连续傅里叶变换的区别，为了获得能正确二维频谱，在每次FFT前都进行了循环移位操作（代码中56-59行和74-77行）。如果不进行这种操作相当于在时域上进行了平移，会在频域产生附加的线性相位，对于成像来说，虽然不会影响成像结果的质量，但会使图像产生整体的平移，影响最终目标绝对位置的判断。

3.多普勒绝对频率的估计

​ 代码中第95-106行利用多普勒中心频率对多普勒绝对频率进行了估计。这也是由于FFT输出的频率轴相当于从零频到 F s ( 采样率 ) Fs(采样率) Fs(采样率)频率的周期循环，而实际频谱应该是以多普勒频率为中心 ± F s / 2 \pm Fs/2 ±Fs/2 范围内的频率，所以需要上述操作。