SAR成像学习(三)距离方向成像matlab代码解析 1

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/qiupingzhao/article/details/51254850


本文将结合matlab代码讲解SAR距离向成像问题。
本文只研究距离向,且是正侧视情况。
文中以同一方位向坐标上四个目标点的成像为例,这四个目标的关系如下:



目标的相关信息:

% 关于目标
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xc=2.e3;                 % Range distance to center of target area
X0=50;                   % target area in range is within [Xc-X0,Xc+X0]
ntarget=4;                        % number of targets
%%%%%%%%%%%%% Targets' parameters  %%%%%%%%%%%%%%%%%%
%
% xn: range;               fn: reflectivity 发射系数
%
xn(1)=0;                   fn(1)=1;
xn(2)=.7*X0;               fn(2)=.8;
xn(3)=xn(2)+2*dx;          fn(3)=1.;
xn(4)=-.5*X0;              fn(4)=.8;
%注意这里的xn是相对于中间的Xc的位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
发射信号是线性调频信号:
% 关于发射信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
colormap(gray(256))
cj=sqrt(-1);
pi2=2*pi;
%
c=3e8;                   % Propagation speed
B0=100e6;                % Baseband bandwidth is plus/minus B0 注意这里的带宽是2B0
w0=pi2*B0;
fc=1e9;                  % Carrier frequency
wc=pi2*fc;

Tp=.1e-6;               % Chirp pulse duration

alpha=w0/Tp;             % Chirp rate
wcm=wc-alpha*Tp;         % Modified chirp carrier 即是式子中的beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

1 问题的描述

成像的过程只要分为两步: 发射信号到接收信号;后处理,即接收信号到影像 如图:



第一步是个正问题,主要由硬件完成,第二步是个逆问题,主要由软件完成。
第一个问题的输入f 0 (x) 是地面目标的理想函数,其与发射信号p(t) 的联合起来之后得到输c出:回波信号s(t) 
第二个问题的输入是s(t) ,通过与p  (t) 匹配滤波,得到输出f(x) 。这个输出就是我们想要的。如何得到f(x) 就是我们的核心问题。
理想情况下f(x)=f 0 (x) ,这意味着我们得到了地面目标函数,也即是我们的影像完全真实地反映了地面的情况。但是,这是不可能实现的。不过,通过合理地解这个逆问题,使得f(x) 接近f 0 (x) 是成像的关键问题,也是我们不断努力的目标。

2 问题的分析和解决



其实问题的解决办法在上一个图中已经暗示了,他就是匹配滤波技术,这个技术我们已经在“急救箱系列”分析过了。但是面对实际问题,还是有许多问题需要说明,尤其是像我这样的新手。

2.1 再谈匹配滤波



匹配滤波的实施公式:
S M (t)=F 1 (ω) [S(ω)P  (ω)] 

其中S(ω),P(ω) 分别为回波信号s(t) 和参考(发射)信号p(t) 的傅里叶变换,p(t) s(t) S(ω) 的形式分别如下:
p(t)=a(t)exp(jβt+jαt 2 ) 

s(t)= n σ n p(t2x n /c) 

S(ω)=P(ω) n σ n exp(jω2x n c ) 

接着推导匹配滤波的实施公式得到:
s M (t) =F 1 (ω) [ n σ n |P(ω)| 2 exp(jω2x n c )]= n σ n Psf t (t2x n c )  

其中Psf t =F 1 (ω) [|P(ω)| 2 ] 表示时间域距离成像点扩散函数。
这样就得到了我们的问题的答案:
f(x)=s M (t)= n σ n Psf t [2c (xx n )] 

这里的转换时用到了时间和距离的比例关系:
x=c2 t 

利用这个关系同样可以得到距离域的点扩散函数。

2.2 采样


这个问题是必然会碰见的问题,因为我们要进行的是“数字处理”。这里的关键问题有:


  • 采样间隔,这个由信号带宽2B 0  决定:
    δt<=1/2B 0 =π/ω 0  

  • 采样起止时刻:
    T s =2(X c X 0 )c  

    T f =2(X c +X 0 )c + T p  

    其中,T p  表示信号持续的时间,最后加上这一项是为了避免出现回波信号中含有部分波的情况。

  • 采样的时间段:
    t[T s ,T f ] 

    T f T s =4X 0 c + T p  

  • 采样点数:
    n=2ceil((.5(TfTs))/dt) 

    取这么多采样点是考虑到离散傅里叶变换的问题。

  • 空间域和频率域采样序列:
    x i =ct i 2 =(X c X 0 )+(i1)cδt2  

    ω i =ω c +(in2 1)δω 

    δω=2πnδt  


该部分的代码是:

% 采样问题
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=pi/(2*alpha*Tp);      % Time domain sampling (gurad band plus minus
                         % 50 per) or use dt=1/(2*B0) for a general
                         % radar signal
%
Ts=(2*(Xc-X0))/c;        % Start time of sampling
Tf=(2*(Xc+X0))/c+Tp;     % End time of sampling

% Measurement parameters
n=2*ceil((.5*(Tf-Ts))/dt);   % Number of time samples
t=Ts+(0:n-1)*dt;             % Time array for data acquisition
dw=pi2/(n*dt);               % Frequency domain sampling
w=wc+dw*(-n/2:n/2-1);        % Frequency array (centered at carrier)
x=Xc+.5*c*dt*(-n/2:n/2-1);   % range bins (array); reference signal is
                             % for target at x=Xc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2.3 结果

四个目标的回波信号对应的代码是:
% 获得回波信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
s=zeros(1,n);              % Initialize echoed signal array
%暂时可以忽略
na=8;                      % Number of harmonics in random phase               
ar=rand(1,na);             % Amplitude of harmonics
ter=2*pi*rand(1,na);       % Phase of harmonics

for i=1:ntarget;
   td=t-2*(Xc+xn(i))/c;
   %注意这里的时间是对采样时间序列的一个延迟,另外有一个平移
   %这个要跟后面的参考信号也就是匹配滤波器对应上,而且他会影响恢复出来的f(x)的解译

   pha=wcm*td+alpha*(td.^2);         % Chirp (LFM) phase
   for j=1:na;                       % Loop for CPM harmonics 相位调制!可以暂时忽略
      pha=pha+ar(j)*cos((w(n/2+1+j)-wc)*td+ter(j));
   end;
   s=s+fn(i)*exp(cj*pha).*(td >= 0 & td <= Tp); %回波信号 注意是各个目标的回波信号的叠加 
   %特别注意这里的Tp,是脉冲持续的时间
end;
匹配滤波用到的参考信号对应的代码是:
% Reference echoed signal
%
td0=t-2*(Xc+0)/c; 
%注意这里的时间参考点
%这个要跟前面的回波信号的时间对应上,而且他会影响恢复出来的f(x)的解译
pha0=wcm*td0+alpha*(td0.^2);         % Chirp (LFM) phase
for j=1:na;                          % Loop for CPM harmonics
   pha0=pha0+ar(j)*cos((w(n/2+1+j)-wc)*td0+ter(j));
end;
s0=exp(cj*pha0).*(td0 >= 0 & td0 <= Tp);
解调及显示:
% Baseband conversion 解调 
% 这里和该系列前面讲的正交解调还是不一样的,当时讨论过
% 为什么会出来了复数图像,现在看来,这个问题有点幼稚。
% 线性调频信号本身就可以表示成复数,回波信号也是。
% 任何一个数只要我们愿意,都可以在形式上表示成复数(无非实数的虚数部分为零嘛)。
% 不应该觉得复数是什么奇怪的事情,而应该觉得它实实在在,就在那里,这是见得少了罢了。
%
sb=s.*exp(-cj*wc*t);
sb0=s0.*exp(-cj*wc*t);

plot(t,real(sb))%只是显示了振幅
xlabel('Time, sec')
ylabel('Real Part')
title('Baseband Echoed Signal')
axis('square')
axis([Ts Tf 1.1*min(real(sb)) 1.1*max(real(sb))])
print P1.1a.ps
pause(1)
%
plot(t,real(sb0))
xlabel('Time, sec')
ylabel('Real Part')
title('Baseband Reference Echoed Signal')
axis('square')
axis([Ts Tf 1.1*min(real(sb0)) 1.1*max(real(sb0))])
print P1.2a.ps
pause(1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
![这里写图片描述]()




注意图中只是画出了振幅信息,而忽略了相位信息,实际上相位信息是非常有用的,特别是在有的应用中。
匹配滤波及结果:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 目标函数重建
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fourier Transform
%
fsb=fty(sb);
fsb0=fty(sb0);

% Power equalization
% 这里是对参考信号做了一个处理,暂时可以先忽略
mag=abs(fsb0);
amp_max=1/sqrt(2);     % Maximum amplitude for equalization
afsb0=abs(fsb0);
P_max=max(afsb0);
I=find(afsb0 > amp_max*P_max);
nI=length(I);
fsb0(I)=((amp_max*(P_max^2)*ones(1,nI))./afsb0(I)).* ...
        exp(cj*angle(fsb0(I)));
%
% Apply a window (e.g., power window) on fsb0 here
%
E=sum(mag.*abs(fsb0));

%
% Matched Filtering
%
fsmb=fsb.*conj(fsb0);
%
%
% Inverse Fourier Transform
%
smb=ifty(fsmb);       % Matched filtered signal (range reconstruction)
%
% Display
%
plot(x,(n/E)*abs(smb))
xlabel('Range, meters')
ylabel('Magnitude')
title('Range Reconstruction Via Matched Filtering')
axis([Xc-X0 Xc+X0 0 1.1]); axis('square')
print P1.6a.ps
pause(1)
结果:


3 如果不采用匹配滤波技术


如果不采用匹配滤波技术,可以用一种叫pulse compression的方法。但是它对信号的形式有一定的限制,比如线性调频信号。上面介绍的匹配滤波方法并不局限于线性调频信号。
代码:

% 另一种方法Time domain Compression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%解调并压缩
td0=t-2*(Xc+0)/c;
scb=conj(s).* ...
     exp(cj*wcm*td0+cj*alpha*(td0.^2));  % Baseband compressed signal
%
plot(t,real(scb))
xlabel('Time, sec')
ylabel('Real Part')
title('Time Domain Compressed Signal')
axis('square')
print TimeDomainCompressedSignal.png
pause(1)

%因为上面的信号包含了与目标的位置相关的频率成分,所以只要进行一个傅里叶变换即可得到问题的解。
fscb=fty(scb);

X=(c*(w-wc))/(4*alpha);     % Range array for time domain compression 注意参考点
plot(X+Xc,(dt/Tp)*abs(fscb)) %?

xlabel('Range, meters')
ylabel('Magnitude')
title('Range Reconstruction Via Time Domain Compression')
axis([Xc-X0 Xc+X0 0 1.1]); axis('square')
%axis([Xc-X0 Xc+X0 0 1.1*max(abs(fscb))]); axis('square')
pause(1)
结果:


说明:
观察上面的结果和匹配滤波得到的结果可以发现明显的不同,这里的原因我暂时也不是特别清楚,也可能是这段代码有问题。但是必须注意到:这肯定是由于信号的压缩造成的。这个后面再回来说。

以上所有的分析中有一个问题被我掩盖了,那就是幅宽问题,如果幅宽很大,以上的分析没有问题,但是如果幅宽相对脉冲持续时间比较小,会带来一个问题,这个问题的解决需要在匹配滤波的过程中添加一些步骤。
未完待续。

阅读更多
换一批

没有更多推荐了,返回首页