%分析SAR成像时的PSF
%编写者:林月冠
%时间:2010-7-12 21:26:26
%%================================================================
clear;clc;close all;
%%================================================================
%%雷达参数
lambda=1.245e-2; %载波波长
C=3e8; %光速
Fc=C/lambda; %载波频率
%%参数--目标区域
Xmin=-100; %方位向目标区域[Xmin,Xmax]
Xmax=100;
Yc=850e3; %图像中心
Y0=100; %地距目标区域[Yc-Y0,Yc+Y0]
%图像宽度 2*Y0
%%参数--平台
Vr=7090; %平台速度100 m/s
H=0; %平台高度 5000 m
R0=sqrt(Yc^2+H^2); %图像中心斜距
%%参数--天线
D=10; %方位向天线长度
Lsar=0.886*lambda*R0/D; %SAR 合成孔径宽度
Tsar=2*Lsar/Vr; %SAR 合成孔径时间
%参数--慢时间域
Ka=2*Vr^2/lambda/R0; %多普勒调频率
Ba=abs(Ka*Tsar); %多普勒带宽
PRF=1.2*Ba/3; %脉冲重复频率
PRT=1/PRF; %脉冲重复时间
dat=PRT; %方位向采样时间间隔
Na=ceil((Xmax-Xmin+Lsar)/Vr/dat); %方位向采样数
sn=linspace((Xmin-Lsar/2)/Vr,(Xmax+Lsar/2)/Vr,Na);%方位向离散时间采样序列
df=-PRF/2:PRF/Na:(PRF/2-PRF/Na); %方位向频谱采样处
dat=PRT;
%%%%参数--快时间域
Br=30e6; %距离向带宽
Fsr=1.1*Br; %距离采样率
Tau=1e-6; %脉冲持续时间
Kr=-Br/Tau; %距离向调频率
drt=1/Fsr; %距离向采样时间间隔
Rmin=sqrt((Yc-Y0)^2+H^2); %距离采样最近斜距
Rmax=sqrt((Yc+Y0)^2+H^2+(Lsar/2)^2); %距离采样的最远斜距
Nr=ceil(2*(Rmax-Rmin)/C/drt+Tau/drt); %距离向的采样数
tm=linspace(2*Rmin/C,2*Rmax/C+Tau,Nr); %距离向离散时间采样序列
fu=-Fsr/2:(Fsr/Nr):(Fsr/2 - Fsr/Nr); %距离向频率采样处
%%参数--分辨率
DY=C/2/Br; %距离向分辨率
DX=D/2; %方位向分辨率
Ntarget=25; %目标数
%%%%现实参数信息和目标属性
disp('参数显示 单位m');disp(' ');
disp('距离向分辨率:');disp(DY)
disp('方位向分辨率:');disp(DX)
disp('SAR合成孔径长度:');disp(Lsar)
%%================================================================
%%基信号
Ia=31;Ir=41;
S0=zeros(Na,Nr);
nn=ceil(numel(S0)*0.1);
Basis=zeros(nn,Ia*Ir);
for k_Ia=1:1:Ia
for k_Ir=1:1:Ir
Dslow=sn*Vr-(k_Ia-16)*PRT*Vr; %各个方位时刻平台相对于目标的方位距离序列
R=sqrt(Dslow.^2+(Yc-(k_Ir-21)*drt*C/2)^2+H^2); %各个方位时刻平台相对于目标的斜距
tau=2*R/C; %各个方位时刻平台到目标的来回时间序列
Dfast=ones(Na,1)*tm-tau'*ones(1,Nr)-Tau/2; %调制时间
phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,Nr)); %见书本P156
S0=exp(j*phase).*(abs(Dfast)<Tau/2).*((abs(Dslow)<Lsar/2)'*ones(1,Nr));
inde=randperm(numel(S0));
ind=sort(inde(1:nn));
Basis(:,(k_Ia-1)*Ir+k_Ir)=S0(ind);
end
end
%%===============================================================
figure,imagesc(abs(Basis));
%%% 计算PSF
PSF=zeros(1,Ia*Ir);
for k=1:Ia*Ir
PSF(k)=abs(Basis(:,636)'*Basis(:,k))/norm(Basis(:,636))/norm(Basis(:,k));
end
hh=figure;
semilogy(PSF,'k');axis tight;grid on;hold on
dd=PSF;dd((dd==max(dd)))=0;dd=max(dd)*ones(size(PSF));
semilogy(dd,'r','LineWidth',2);
title('PSF');ylim([10^-2.5 1])
set(hh,'units','centimeters');
set(hh,'position',[5 5 11 7]);
%%===========================CS重建===============================
Ns=Ia*Ir; %场景网格的点数
Nt=10; %场景中的目标点数
x=zeros(Ns,1);
inde=randperm(Ns);
ind=sort(inde(1:Nt));
x(ind)=1;
y=Basis*x; %回波数据,Basis为观测矩阵
% 重建观测场景
tic
Theta=[real(Basis) -imag(Basis);imag(Basis) real(Basis)];
Y=[real(y);imag(y)];
sigma2 = std(Y)^2/1e2;
eta = 1e-20;
lambda_init = [];
adaptive = 0;
optimal = 1;
[weights,used]=FastLaplace(Theta,Y,sigma2,eta,lambda_init, adaptive,optimal);
out=zeros(1,2*Ns);out(used)=weights;
xr=out(1:Ns)+1i*out(Ns+1:end);
scene=reshape(xr,31,41); %按31*41排列原来的1271*1的向量
toc
hh=figure;
imagesc(abs(scene));axis tight;axis off
colormap(copper)
set(hh,'units','centimeters');
set(hh,'position',[5 5 11 7]);
- 1
- 2
前往页