CS算法（三）—— SAR卫星回波数据成像处理

SAR成像专栏目录

本期将使用CS算法对哨兵1号和RADARSAT-1的回波数据进行SAR成像处理，展示成像结果，并附上MATLAB代码，CS算法的核心代码见本专栏的《CS算法（二）—— 斜视SAR点目标仿真》

1. 哨兵1号SAR回波数据成像处理

1.1 成像结果

成像区域为美国旧金山，成像结果如下。由于数据较大（62510*34856）切分为了两次进行成像：

对应位置的光学图如下，SAR图像未做几何校正，所以不是正南正北的，与光学图有个夹角。

SAR图像局部放大如下，这张图里拍到了金门大桥。

对应光学图为：

1.2 成像代码

close all; 
clc;clear;

%% read file
load C:\xxx\echo2.mat
%% bu 0
din=double(echo(30001:end, :));
N=ceil(para.Tc*para.Fs);
[Nslow0,Nfast0]=size(din);
Nfast=Nfast0+2*N;
bu0=zeros(Nslow0,N);
din=[bu0,din,bu0];
bu00_num=2000;
Nslow=Nslow0+bu00_num;
bu0=zeros(bu00_num/2,Nfast);
echo=[bu0;din;bu0];
figure; imagesc(abs(echo));
clear bu0 din din2;
fprintf('回波数据导入完毕\n\n');

%% 参数
Vs=sqrt(para.Vx.^2+para.Vy.^2+para.Vz.^2);
Vs=Vs(10);
RR=sqrt(para.Px.^2+para.Py.^2+para.Pz.^2);
RR=RR(10);
% Re=6378.13649e3; %Semi-major axis of the earth
% Re=6356.752e3; %Semi-minor axis of the earth
Re=6371.393e3;  %Mean radius of earth

Vg=Vs*Re/RR;
para.V =sqrt(Vs*Vg);

R0=(para.delay_pulse_num*para.PRT+para.sample_delay)*para.c/2;
r_gate=para.c/para.Fs/2;
para.Rmin=R0-N*r_gate;
para.bw=para.bw_a;
para.fdc=0;
para.PRF=1/para.PRT;
height=RR-Re;
alpha=acosd((height/(R0+para.Nr/2*r_gate)));
para.Kr =      -para.Kr;
para.Vst =     para.V ;    
para.Rcc =    para.Rmin+Nfast/2*r_gate ;  
para.theta =   0/180*pi ; 

%%
disp('开始成像......');
tic
Sa = CSA_linear(echo,para);
toc

figure;imagesc(abs(Sa));
fprintf('CS成像完毕\n');

SSa = abs(flip(Sa,2));


%% 存储tiff图像
ra_length=round(max(ra_L1)); rr_length=round((max(rr_L1)-min(rr_L1)));
image1=imresize(SSa, [round(Nslow*ra_length/rr_length), Nslow], 'nearest');

gen_tiff_img(image1, 16, 'C:\xxx\image.tiff');
clear image1;

使用普通电脑和MATALB成像每景图大约用时2分钟，使用RD算法需要的时间预计得超1小时，所以使用CS算法可以更快地对更大的区域进行精确的成像，这是CS算法的初衷。

由于哨兵1号卫星的多普勒导引做的比较好，SAR回波数据的多普勒中心就在0附近，不需要做距离走动校正，不足以验证我们CS算法的好坏，接下来我们再以RADARSAT-1的回波数据为例进行CSA成像处理：

2. RADARSAT-1 SAR回波数据成像处理

2.1 成像结果

我们这次把之前加拿大温哥华地区的回波数据全部成像出来。

对应光学图为：

局部放大为：

2.2 成像代码

clear ;	
close all; 
clc;

%% parameter
c=299792458;                    	%电磁波传播速度
fc=5.3e+09  ;                   	%雷达频率,GHz
lambda=c/fc;

Fs=32.317e6; Ts=1/Fs;                 	%采样率,MHz
Br=30.111e6;
Kr=-0.72135e12;                	%距离调频率,MHz/us
Tc=41.74e-6 ;                  	%脉宽,us
PRF=1.25698e3 ;  PRT=1/PRF;       	%脉冲重复频率, Hz
Nr=round(Fs*Tc/2)*2;
tstart=6.5956e-3;              	%数据窗开始时间,ms
first_rg_cell=1850;
r_gate=c*Ts/2;
R0=tstart*c/2;%+first_rg_cell*r_gate;
V=7062;                        	%有效雷达速率,m/s
fd=2*V/lambda;
Ka=1733;                        	%方位调频率,Hz/s
fdc=-6900;                       	%多普勒中心频率,Hz
bw=2.9/180*pi;                      %波束宽度，弧度   
D=0.886*lambda/bw;
theta=40.15;

%% load data
% load CDdata1.mat
load C:\xxx\echo.mat
din=double(echo);

%% bu 0
[Nslow0,Nfast0]=size(din);
Nfast=Nfast0+2*Nr;
bu0=zeros(Nslow0,Nr);
din2=[bu0,din,bu0];
bu00_num=1000;
Nslow=Nslow0+bu00_num;
bu0=zeros(bu00_num/2,Nfast);
echo=[bu0;din2;bu0];
figure; imagesc(abs(echo));




%% 
para.c=c;
para.lambda=lambda;
para.Fs  =Fs  ; 
para.Kr=Kr;
para.V   =V   ;
para.PRF =PRF ;
para.Rmin=R0-Nr*r_gate;
para.bw=bw;
para.theta=asin(fdc/(2*V/lambda));

SS=CSA_linear(echo,para);


fprintf('方位向脉压结束\n\n');

SS=flip(SS(:,Nr/2:end-Nr/2));


%% 存储tiff图像

ra_length=round(max(ra_L1)); rr_length=round((max(rr_L1)-min(rr_L1)));
image1=imresize(img, [round(Nslow*ra_length/rr_length), Nslow], 'nearest');

gen_tiff_img(image1, 16, 'C:\Desktop\MATH\image\image0712-6.tiff');
clear image0 image1;

最后给大家附上生成tiff图的代码：

function gen_tiff_img(img_in, bits, path)


mmax=max(abs(img_in),[],'all');
scale=mmax/2^bits;
image=img_in/scale;
if bits==8
    data_tiff = uint8(image);
elseif bits==16
    data_tiff = uint16(image);
% elseif bits==24
%     data_tiff = uint24(image);
elseif bits==32
    data_tiff = uint32(image);
end
% data_tiff=single(image);

t = Tiff(path,'w');
setTag(t,'Photometric',Tiff.Photometric.MinIsBlack);
setTag(t,'Compression',Tiff.Compression.None);
% % setTag(t,'YCbCrSubSampling',[2 2]);
setTag(t,'BitsPerSample',bits);
setTag(t,'SamplesPerPixel',1);
setTag(t,'SampleFormat',Tiff.SampleFormat.UInt); %IEEEFP
setTag(t,'ImageLength',size(data_tiff,1));
setTag(t,'ImageWidth',size(data_tiff,2));
% % setTag(t,'TileLength',32);
% % setTag(t,'TileWidth',32);
setTag(t,'PlanarConfiguration',Tiff.PlanarConfiguration.Chunky);
% % setTag(t,'JPEGColorMode',Tiff.JPEGColorMode.RGB);
% % setTag(t,'JPEGQuality',75);

write(t,data_tiff);
close(t);