www.pudn.com > real-pictures-imaging.zip > SAR_nichuli.m, change:2012-03-07,size:8582b


%SAR回拨模拟及CS算法成像,参考书:合成孔径雷达成像原理(皮亦鸣 
clc  
clear all;  
close all;  
clc  
clear all;  
close all;  
%-----------------------------目标位置-----------------------------------%  
a=imread('test5.jpg');             % 载入图片 
b=imresize(a,1); 
I=rgb2gray(b); 
hang=size(I,1); 
lie=size(I,2); 
 
k=1; 
for m=1:6:hang 
    for n=1:6:lie 
        if I(m,n)~=255 
            X0(k)=m; 
            R0(k)=n; 
            k=k+1; 
        end 
    end 
end 
mid_X0=(X0(1)+X0(k-1))/2; 
%-----------------------------------------------------------------------%  
Npt_num = length(X0);                 %目标点数 =2 
%---------------------------雷达平台参数-----------------------------------%  
 
% H=input('输入雷达飞行高度值:');               % 初始雷达高度(z轴:m),要求高度值在3000~5000 
% V=input('输入雷达飞行速度值:');               % 设置雷达飞行平台的速度(x轴:m/s),要求速度值100~300 
% ThetaCenterdeg=input('输入雷达下视角度数:');  % 正侧视SAR系统的下视角,设置雷达照射区域中心,要求角度值在30~60 
% ThetaCenter=ThetaCenterdeg/180*pi;            % 转换为弧度值 
% DX=input('输入雷达方位向分辨率:');            % 设置用户要求的方位向分辨率,要求分辨率达到0.5~1 
% DY=input('输入雷达距离向分辨率:');            % 设置用户要求的距离向分辨率,要求分辨率达到0.5~1 
% fc=input('输入雷达工作频率即载频');           % 根据用户要求设置雷达工作频率,选X波段8GHz~12GHz 
% Tr=input('输入雷达信号脉冲宽度');             % 脉冲宽度,1.5us 
 
H=10000;                                 %雷达高度 m  
Theta=30*pi/180;                        %波束入射角  
Theta_rc=0*pi/180;                      %波束斜视角  
fc=4e9;                                 %载频  
resR=1;                                 %距离向分辨率 
Va=300;                                 %雷达速度  
resA=1;                                 %方位向分辨率  
Tp=3e-6;                                %脉冲宽度  
 
Ra=H/cos(Theta);                        %景中心斜距   
R=H*tan(Theta);                         %景中心斜距地面投影  
c=3e8;                                  %光速  
lambda=c/fc;                            %载波波长 
 
%---------------------------距离向参数-----------------------------------% 
 
FsFactor=1.2;                           %过采样因子  
Br=c/2/resR;                            %发射信号带宽  
Kr=Br/Tp;                               %距离向线性调频率  
Fs=Br*FsFactor;                         %距离向抽样频率  
gal=c/2/Fs;                             %距离门  
nr=floor(Fs*Tp)                         %脉冲抽样点数,floor函数即向负无穷方向取整 
Nr=nr+2000;                             %距离向点数  
PTheta=3*0.886*lambda/(2*resR);         %合成孔径角  
  
Rmax=H/cos(Theta+PTheta/2);  
Rmin=H/cos(Theta-PTheta/2);  
nnr=fix((Rmax-Rmin)/gal)                %fix向0靠拢取整 
  
RRmax=H*tan(Theta+PTheta/2);  
RRmin=H*tan(Theta-PTheta/2);  
%gale=gal*sin(Theta); 
gale=(RRmax-RRmin)/nnr;  
%-----------------------------------------------------------------------%  
  
%---------------------------方位向参数-----------------------------------%  
 
La=2*resA;                              %天线长度  
Ls=0.886*Ra*lambda/La;                  %合成孔径长度 Ra=H/cos(Theta);载波波长 
Ts=Ls/Va;                               %合成孔径时间  
Theta_syn=0.886*lambda/La;              %合成孔径角  
  
  
Ta=Ts;                                  %测量时间  
Xs=Ta*Va;                               %测量长度  
fdc=2*Va*sin(Theta_rc)/lambda;          %多普勒中心频率,波束斜视角 
Ka=2*Va^2*(cos(Theta_rc))^3/lambda/Ra;  %多普勒调频率  
Bd=abs(Ka)*Ts                           %多普勒带宽  
PRF=Bd*2;                               %脉冲重复频率  
Na=fix(PRF*Ta);                         %方位向采样点数Na  
Na=Na+mod(Na,2); 
%-----------------------------------------------------------------------%  
  
Xt=linspace(-Ta/2,Ta/2,Na);  
ra=zeros(Npt_num,Na);  
for i=1:Npt_num  
        ra(i,:)=sqrt((R+R0(i)).^2+((Va)*Xt-(X0(i)-mid_X0)).^2+H^2);  
end  
rmax=max(max(ra));						% max. slant range   最大斜距  
rmin=min(min(ra));						% min. slant range   最小斜距  
rmc=fix((rmax-Rmin)/gal);			    % range migration	最大距离徙动  
rg=fix((ra-Rmin)/gal+1);				% range gate beause of range migration 距离徙动  
rgmax=max(max(rg));  
rgmin=min(min(rg));  
  
  
tr=linspace(-Tp/2,Tp/2,nr);								% fast time along range  距离快时间  
sig=zeros(Na,Nr);  
for i=1:Na			  
	for k=1:Npt_num  
        RM=sqrt(H^2+(Ls/2)^2+(R+R0(k))^2);             %照射时间内目标与雷达的最大斜距  
        if ra(k,i)<RM  
            phase=(tr).^2;  
            sig(i,rg(k,i):rg(k,i)+nr-1)=sig(i,rg(k,i):rg(k,i)+nr-1)+exp(-j*4*pi*ra(k,i)/lambda)*exp(-j*pi*Kr*phase);  
%   		sig(i,rg(k,i):rg(k,i)+nr)=sig(i,rg(k,i):rg(k,i)+nr)+exp(-j00*4*pi/lamda*ra(k,i))*exp(-j00*pi*Kr*(tr).^2);  
        end  
    end    
end  
disp('end of echo return modelling')                %回波建模  
figure(1)  
imshow(real(sig))       %显示一张二值图像,imshow是matlab中显示图像的函数 
  
%-----------------------------方位向FFT------------------------------------%  
for i=1:Nr  
   sig(:,i)=fftshift(fft(sig(:,i)));		  
end  
disp('end of azimuth fft')                                   
figure(2)  
% mesh(abs(sig))  
imagesc(abs(sig)),colormap(gray); %imagesc(C)将输入变量C显示为图像,colormap设定和获取当前的颜色 
%plot(real(sig(:,436)));  
  
%-----------------------------Chirp Scaling--------------------------------%  
rref=Ra;                                            %参考距离  
fa=linspace(-PRF/2,PRF/2,Na);                       %方位向频率  
Kref=Kr./(1+Kr*rref*2*lambda*(lambda*fa/2/Va).^2/c^2./(1-(lambda*fa/2/Va).^2).^(3/2));      %等效的距离调频率  
Cs=1./sqrt(1-(lambda*fa/2/Va).^2)-1;                   %弯曲因子  
Rref=rref*(1+Cs);                                   %参考距离处的距离曲线  
t=zeros(1,Nr);  
for i=1:Nr  
    t(i)=2*sqrt((RRmin+(i-1)*gale)^2+H^2)/c;  
     %t(i)=2*sqrt((Rmin+(i-1)*gal*2)^2)/c;  
end  
CSf=zeros(Na,Nr);  
for i=1:Nr  
    CSf(:,i)=exp(-j*pi*Kref.*Cs.*(t(i)-2*Rref/c).^2);  
end  
sig=sig.*CSf;  
disp('end of CS')    
  
%--------------------------------距离向FFT-----------------------------------%  
for i=1:Na  
   sig(i,:)=fftshift(fft(sig(i,:)));		  
end  
disp('end of range fft')                 
  
  
%-------------------------------------距离压缩和距离徙动校正------------------%  
RMCR=zeros(Na,Nr);  
fr=linspace(-Fs/2,Fs/2,Nr);  
for i=1:Nr  
    RMCR(:,i)=exp(-j*pi*fr(i)^2./(Kref.*(1+Cs))).*exp(j*4*pi*fr(i)*rref*Cs/c);  
end  
sig=sig.*RMCR;  
disp('end of RMCR')   
  
%---------------------------距离向IFFT---------------------------------------%  
for i=1:Na  
   sig(i,:)=(ifft(sig(i,:)));		  
end  
disp('end of range ifft')     
figure(3)  
imagesc(abs(sig)),colormap(gray);  
  
%------------------------------------方位压缩---------------------------------%  
tt=ones(Na,1)*t;  
DT=4*pi/c^2*Kref'.*(1+Cs)'.*Cs'*(t*c/2-rref);  
ACP=exp(j*4*pi/lambda*sqrt(1-(lambda*fa'/2/Va).^2)*t*c/2).*exp(-j*DT).*exp(-j*2*pi*c/lambda*tt);  
sig=sig.*ACP;  
disp('end of amizuth compress')     
  
%-----------------------------------方位向IFFT-------------------------------%  
for i=1:Nr  
   sig(:,i)=(ifft(sig(:,i)));		  
end  
disp('end of azimuth fft')         
figure(4)  
imagesc(abs(sig)),colormap(gray);  
 
 
  
%imshow和image: 图像的显示是最为重要的,用imshow和image都可以显示图像,但是有一定的区别。用的不对,就会象我最初一样,老是出错,或者得到一张空白图或者是彩色图显示成颗粒状、反相黑白图等等。image是用来显示附标图像,即显示的图像上有x,y坐标轴的显示,可以看到图像的像素大小。imshow只是显示图像。它们都可以用subplot来定位图像显示的位置,用colormap来定义图像显示用的颜色查找表,比如用colormap(pink),可以把黑白图像显示成带粉红色的图像,很有趣的。在这里最值得注意的是要显示的图像像素矩阵的数据类型。显示真彩色图像像素三维矩阵X,如果是uint8 类型,要求矩阵的数据范围为0-255,如果是double型,则其数据范围为0-1,要不就会出错或者出现空白页。类型转换很简单,如果你原来的数值是 uint8,在运算中转换为double后,实际要显示的数值没有改变的话,只要用uint8(X)就可转换为uint8型,如果不想转换频繁,也可在显示时用X/255来转换为符合0-1double类型范围要求的数值显示。如果显示索引图像(二维矩阵),因为不同数据类型对应颜色查找表colormap的基点不同,会有所区别,如果不对的话,会出现很多意外的显示效果的。如果索引图像像素数值是double型,则它的取值范围为1-length(colormap),数值起点为1,则矩阵中数值为1的对应colormap中第一行数据,如果索引图像像素数值是uint8,则取值范围为0-255,数值起点为0,则矩阵中数值为0的对应colormap中第一行数据,所以索引图像这两个数据类型之间的转换,要考虑到+1或-1。直接用uint8或double转换则会查找移位,产生失真情况。uint16数据类型与uint8类似,取值范围为0-65536。