SAR面目标的回波仿真

%% 0 参数
fc   = 10e9;            % 10 GHz
c    = 3e8;
lambda = c/fc;
B    = 150e6;           % 150 MHz
fs   = 180e6;           % 采样率
Tp   = 1e-6;            % 脉冲宽度
PRF  = 500;             % 500 Hz
v    = 150;             % 机速 150 m/s
h    = 8000;            % 高度 8 km
R0   = 800e3;           % 最近斜距
thetaDep = 30*pi/180;   % 下视角 30°
​
%% 1 生成目标网格（二选一）
% [node,elem,rcs] = target_flatRoof(20,15,12);   % 长20 宽15 高12 m
[node,elem,rcs] = target_triCone(10,8);          % 底径10 m 高8 m
​
%% 2 计算二维频域回波
[Echo,f_range,f_az] = sarEcho3D(node,elem,rcs, ...
                                fc,B,PRF,v,h,R0,thetaDep);
​
%% 3 RD 成像
[Image,dRange,dAz] = rdImaging(Echo,fc,B,Tp,fs,PRF,v,h,R0);
​
%% 4 可视化
plotResults(Echo,Image,node,elem,dRange,dAz);

function [node,elem,rcs] = target_flatRoof(L,W,H)
% 将平顶楼拆为 6 个矩形面 -> 三角面片
% 输出：node (N×3)  elem (M×3)  rcs (M×1) 复散射系数
[x,y] = meshgrid(0:W/20:W, 0:L/20:L);
% 顶面
zTop = H*ones(size(x));
[nodeTop,elemTop] = rect2tri(x,y,zTop);
% 前后立面
zFront = linspace(0,H,size(x,1))';
[nodeFront,elemFront] = rect2tri(x(1,:),y(1,:),zFront);
% ... 左右底面同理，此处省略
node = [nodeTop; nodeFront; ...];
elem = [elemTop; elemFront+size(nodeTop,1); ...];
% 物理光学散射系数（简化版）
A = triArea(node,elem);              % 面片面积
rcs = 1j*4*pi/lambda * A .* cos(30*pi/180); % PO 近似
end

function [Echo,f_r,f_a] = sarEcho3D(node,elem,rcs,fc,B,PRF,v,h,R0,thetaDep)
% 面片→二维频域相位历史（PO + 频域叠加）
Nr   = 512;              % 距离采样
Na   = 1024;             % 方位采样
f_r  = (-Nr/2:Nr/2-1)*(B/Nr);          % 距离频率
f_a  = (-Na/2:Na/2-1)*(PRF/Na);        % 方位频率
Echo = zeros(Nr,Na);
[kx,ky,kz] = sarWaveVector(fc,f_r,f_a,v,h,R0,thetaDep);
for k = 1:size(elem,1)
    % 面片重心
    cent = mean(node(elem(k,:),:),1);
    % 波矢量
    kVec = [kx(:), ky(:), kz(:)];
    % 相位 = 2*k·r
    phase = 2 * (kVec * cent');
    % 叠加
    Echo = Echo + rcs(k)*exp(1j*phase);
end
end

function [Img,dr,daz] = rdImaging(Echo,fc,B,Tp,fs,PRF,v,h,R0)
% 1) 距离向匹配滤波
Nr = size(Echo,1);
fr = (-Nr/2:Nr/2-1)*(fs/Nr);
Hr = exp(1j*pi*fr.^2*Tp/B);           %  chirp 匹配
Echo_r = ifft(fft(Echo,[],1) .* Hr.', [],1);
​
% 2) 方位向 FFT（RD 算法）
Img = fftshift(fft(Echo_r,[],2),2);
% 像素尺寸
dr   = c/2/B;                         % 距离
daz  = lambda*R0/2/size(Echo,2)/v;    % 方位
end