TomoSAR仿真

基于压缩感知的层析SAR仿真

   层析SA原理：通过多基线（多通道）实现高度向的高分辨。在得到多幅二维SAR图像后，配准、去平地相位，之后利用谱分析或压缩感知的方法重建出目标高度。

只考虑雷达与目标在同一个方位切片的情况。

参数	载频
数值	10GHz
高度	2000m
带宽	80MHz
脉冲宽度	8㎲
基线间隔	1m
通道数	5

下面展示一些 内联代码片。

close all;
H = 2000;
c=3e8;            %光速
Fc=10e9;
lamda=c / Fc;      %波长
Tr=8e-6;    %脉冲宽度
Br=80e6;
Kr=Br/Tr;   %chirp调频率
Yc = 2000;
target  = [0,2000,0;
    0,2100,10;
    0,2200,20];  %目标位置
R0 = sqrt((H-target(1,3))^2 + target(1,2)^2);
R1 = sqrt((H-target(2,3))^2 + target(2,2)^2);
R2 = sqrt((H-target(3,3))^2 + target(3,2)^2);
N = size(target,1);
Nrg = 1024;
tr =linspace(-Tr/2,Tr/2,Nrg)+2*R0/c;
r = tr*c/2;
fsr = 1/(Tr/Nrg);
a = fsr/Br;
dr = c/2/Br;
% 基线设置
del_H = 1;
Npass = 5;
theta_r = atan(2000/H);
dz = lamda * R0 / (2 *(Npass-1)* del_H * sin(theta_r))  * sin(theta_r)    % 0.4m
Hmax = (Npass-1) * dz %最大不模糊高程


%回波信号
s_echo = zeros(Npass,Nrg);
for k = 1:Npass
    for ii = 1:N
        R= sqrt(target(ii,2)^2 + (H+del_H*(k-1)-target(ii,3))^2 ) * ones(1,Nrg);
        tau = tr-2*R/c;
        s = exp(-1j*4*pi*R/lamda) .* exp(1j* pi* Kr*tau.^2) ;
        s_echo(k,:) = s_echo(k,:) + s;
    end
end
clear s
%figure;plot(r,abs(s_echo));xlabel('距离');
rr = sqrt(target(:,2).^2 + (H+del_H*(0:Npass-1)-target(:,3)).^2 );
%距离压缩
tm = tr-2*R0/c;
Refr=exp(1j*pi*Kr*tm.^2);%距离压缩参考函数
F_Refr=fft(Refr);
s_rd = zeros(Npass,Nrg);
for k = 1:Npass
    S_echo = fft(s_echo(k,:));
    s_rd(k,:) = fftshift(ifft(S_echo .* conj(F_Refr)));
end
% figure;plot(r,20*log10(abs(s_rd)));xlabel('距离');
fprintf('配准前的相关系数：');
for k=2:Npass
    cor = abs(corrcoef(abs(s_rd(1,:)),abs(s_rd(k,:))));
    fprintf('%1.2f ',cor(2,1));
end
fprintf('\n');
%% 粗配准
s_rd2 = s_rd;
%基于最大相关系数法像素级配准
s_rd3 = s_rd2;
Nx = 30;
dx = 1;
s1 = s_rd2(1,:);
Coeff = zeros(Npass,Nx);
shift = zeros(1,Npass);
for k = 2:Npass
    for n = 1:Nx
        s2 = s_rd2(k,:);
        S2 = fftshift(fft(s2));
        S2 = S2 .* exp(1j*2*pi*(-Nrg/2:Nrg/2-1)/Nrg * (n -Nx/2-1)*dx);   %搬移
        s2 = ifft(ifftshift(S2));
        Coeff(k,n) = sum(s1 .* conj(s2))/sum(abs(s1).^2) /sum(abs(s2).^2);  
    end
    [temp, indx] = max(Coeff(k,:));
    shift(k) = indx-Nx/2-1;
    s2 = s_rd2(k,:);
    S2 = fftshift(fft(s2));
    S2 = S2 .* exp(1j*2*pi*(-Nrg/2:Nrg/2-1)/Nrg * shift(k)*dx);
    s2 = ifft(ifftshift(S2));
    s_rd3(k,:) = s2;
end
% figure;plot(r,20*log10(abs(s_rd3(1:2,:))));xlabel('距离');
fprintf('粗配准后的相关系数：');
for k=2:Npass
    cor = abs(corrcoef(abs(s_rd3(1,:)),abs(s_rd3(k,:))));
    fprintf('%1.2f ',cor(2,1));
end
fprintf('\n');
%%  基于最大干涉频谱法亚像素级配准
s_rd4 = s_rd3;
Nx = 18;
dx2 = 0.1;
s1 = s_rd3(1,:);
SNR = zeros(Npass,Nx);
shift2 = zeros(1,Npass);
for k = 2:Npass
    for n = 1:Nx
        s2 = s_rd3(k,:);
        S2 = fftshift(fft(s2));
        S2 = S2 .* exp(1j*2*pi*(-Nrg/2:Nrg/2-1)/Nrg * (n -Nx/2-1)*dx2);   %搬移
        s2 = ifft(ifftshift(S2));
        tmp = angle(s1 .* conj(s2));
        tmp =fftshift(fft(exp(1j*tmp)));
        [temp , indx] = max(abs(tmp));
        tmp(indx) = 0;
        SNR(k,n) = temp/mean(abs(tmp));  
    end
    [temp, indx] = max(SNR(k,:));
    shift2(k) = indx-Nx/2-1;
    s2 = s_rd3(k,:);
    S2 = fftshift(fft(s2));
    S2 = S2 .* exp(1j*2*pi*(-Nrg/2:Nrg/2-1)/Nrg * shift2(k)*dx2);
    s2 = ifft(ifftshift(S2));
    s_rd4(k,:) = s2;
end
figure;plot(r,20*log10(abs(s_rd4)));xlabel('距离');title('精配准后');
fprintf('精配准后的相关系数：');
for k=2:Npass
    cor = abs(corrcoef(abs(s_rd4(1,:)),abs(s_rd4(k,:))));
    fprintf('%1.2f ',cor(2,1));
end
fprintf('\n');

%%  去斜
Rn0 = sqrt(ones(Npass,1) * (r.^2 - H^2)  +(H+(0:Npass-1)*del_H)' .^2 * ones(1,Nrg) ) ;
s_rd5 = s_rd4 .* exp(1j*4*pi/lamda * Rn0);

%画图
figure;
subplot(1,2,1);xlabel('距离');title('去斜后相位');
for k=1:Npass
    hold on; plot(r,phase(s_rd5(k,:)));
end
hold off;
subplot(1,2,2);xlabel('距离');title('去斜后的干涉相位');
for k=2:Npass
    hold on; plot(r,phase(s_rd5(1,:) .* conj(s_rd5(k,:))));
end
hold off;
%% CS高程重建
Nz = 128;
s_cs= zeros(Nz,Nrg);
z = linspace(-Hmax/2,Hmax/2,Nz);
z =Hmax*((-Nz/2+1:Nz/2)+20)/Nz;
 %z =((-Nz/2:Nz/2));
K = [1,2,3,4,5];
    for n=1:Nrg
        data = s_rd5(:,n);
        s_cs(:,n) = cs_omp(data,Nz);
    end
for i=1:Nz
    for j=1:Nrg
        if abs(s_cs(i,j))< T*0.4
            s_cs(i,j) = 0;
        end
    end
end
figure,
imagesc(r,z,abs(s_cs));%colormap(gray);
xlabel('斜距');ylabel('高度');axis xy;title('反演结果');