本章內容

重點

•SAR成像幾何關係；SAR回波模型；SAR成像模型；

•RD成像演算法；SAR成像質量；

•距離徙動；距離校正；

•實時成像；回波模擬；

要求

•掌握SAR回波模擬的過程；

•理解距離徙動的概念；

•掌握RD成像演算法的過程。

一、成像模型

對合成孔徑原理的理解

l方位向虛擬陣列天線合成大孔徑

l方位向多普勒頻率變化形成LFM訊號

》》》雷達與目標發生相對運動，引起兩者相對位置的改變，這是實現合成孔徑的物理基礎；

  SAR成像主要圍繞以上相對位置（方位向+距離向）的變化展開研究工作。

• 回波訊號中的兩個時間變數

l距離向時間τ，一般為us

l方位向時間t，一般為ms量級，慢時間。

• 根據時間的定義，回波接收與平臺運動過程可假設為START-STOP過程

• 即雷達在每一個脈衝發射與接收過程中方位向保持靜止，這便於理解回波形成過程，也便於後續成像演算法的設計與理解。

% 名稱： 基於RD演算法的雷達成像模擬
% 作者： _Sure_  
% 時間： 2016.05.25 

%%產生Stripmap SAR的回波
close all;  
clear all;  
clc;  
  %{  
        close all 是關閉所有視窗(程式執行產生的，不包括命令窗，editor窗和幫助窗）  
        clear all是清除所有工作空間中的變數  
        clc 是清除命令列  
        clear 是清除那一堆存在記憶體裡的資料的，  
        close 是關閉打開了的檔案，一般檔案開頭寫  
 %}  
 

%---------------------------------start-----------------------------------%  
%-------------------------------基本引數設定-------------------------------%  
thetaT=0;                              %T平臺波束斜視角
thetaT=thetaT*pi/180;                  %rad弧度
thetaR=0;                              %R平臺波束斜視角
thetaR=thetaR*pi/180;                  %rad弧度
c=3e8;                                 %光速
fc=1.5e9;                    %載波頻率1.5G
lambda=c/fc;                 %載波的波長
 

%%-----------測繪帶區域----------------
X0=？;          %方位向[-X0,X0]，範圍有自己根據實際情況確定
Rtc=3000;       %載波發射距離
Rrc=3000;       %載波接收距離
Rc=(Rtc+Rrc)/2; %載波距離向
R0=150;         %距離向[Rc-R0,Rc+R0]

%%-----------距離向(Range),r/t domain---
Tr=1.5e-6;  %LFM訊號脈寬 1.5us (200m)
%{
  概念普及：
  LFM訊號：(Linear-freuency-modulation)
  線性調頻訊號
%}
Br=150e6;     %LFM訊號頻寬 150MHz
Kr=Br/Tr;     %調頻斜率
Nr=512;       %快時間取樣點數

%+++++距離域序列+++++
r=Rc+linspace(-R0,R0,Nr);    %距離域序列
%{
  概念普及：
  linspace(-A,A,M)：生成一些從—A到+A的M份等間距的數值
  舉例：
     linspace(-5,5,6)
     ans =
        -5    -3    -1     1     3     5
%}
t=2*r/c;                       %距離域t值對應
dt=R0*4/c/Nr;                  %快時間取樣週期
%+++++頻率域序列+++++
f=linspace(-1/2/dt,1/2/dt,Nr); %f域序列

%%----------方位向(Azimuth,Cross-Range),x/u domain----
v=100;           %SAR 平臺移動速度
%{
  概念普及：
  合成孔徑雷達——SAR 
               ——Synthetic Aperture Radar 的縮寫
  合成孔徑觀點：  用運動的小天線合成一個大天線
  SAR在方位向是真正的微波全息;
  通過PRF進行方位向全息圖的取樣;
  取樣間隔必須固定不變，因此要求PRF與地速成正比.
%}
Lsar=300;             %合成孔徑長度
Na=1024;              % 慢時間取樣點數
%+++++u域序列+++++
x=linspace(-X0,X0,Na);         %u域序列
u=x/v;                         %u域序列t值對應
du=2*X0/v/Na;                  %慢時間取樣間隔
%+++++fu域序列+++++
fu=linspace(-1/2/du,1/2/du,Na);%fu域序列
ftdc=v*sin(thetaT);           %SAR-T平臺波束速度
ftdr=-(v*cos(thetaT))^2/lambda/Rtc;
frdc=v*sin(thetaR);           %SAR-R平臺波束速度 
frdr=-(v*cos(thetaR))^2/lambda/Rrc;
fdc=ftdc+frdc;                %Doppler調頻中心頻率
fdr=ftdr+frdr;                %Doppler調頻斜率
%-------目標位置-----------
Ntar=3;%目標個數
Ptar=[  Rrc ,    0 ,  1 %引數對應：距離向座標,方位向座標,sigma             
      Rrc+50,   -50,  1
      Rrc+50,    50,  1];

%%-----------LFM產生回波---------
s_ut=zeros(Nr,Na);  %生成m×n的double類零矩陣，設定數值記憶體
%{
  概念普及：
  zeros:
  zeros(m,n)產生m×n的double類零矩陣，zeros(n)產生n×n的全0方陣。
  舉例
    zeros(2,4)
    ans =
         0     0     0     0
         0     0     0     0
    zeros(4)
    ans =
         0     0     0     0
         0     0     0     0
         0     0     0     0
         0     0     0     0
%}
U=ones(Nr,1)*u;            %擴充為矩陣
T=t'*ones(1,Na);           % 將快時間拓展為 Na 列
%{
  概念普及：
  ones:
      ones函式——生成全1陣
  zeros(m,n)產生m×n的double類零矩陣，zeros(n)產生n×n的全0方陣。
  舉例
    ones(1,6)
    ans =
         1     1     1     1     1     1
%}
for i=1:1:Ntar
    rn=Ptar(i,1);          % 目標距離向座標
    xn=Ptar(i,2);          % 目標方位向座標
    sigma=Ptar(i,3);       % 目標RCS
    rtn=rn+Rtc-Rrc;        
    RT=sqrt(rtn^2+(rtn*tan(thetaT)+xn-v*U).^2);     % 發射 目標斜距
    RR=sqrt(rn^2+(rn*tan(thetaT)+xn-v*U).^2);       % 接收 目標斜距
    R=RT+RR;
    DT=T-R/c;
    phase=pi*Kr*DT.^2-2*pi/lambda*R;
    s_ut=s_ut+sigma*exp(j*phase).*(abs(DT)<Tr/2).*(abs(v*U-xn)<Lsar/2);
end;
%-------------------------------end---------------------------------------% 

%%-------------距離壓縮-------------
%參考訊號
p0_t=exp(j*pi*Kr*(t-2*Rc/c).^2).*(abs(t-2*Rc/c)<Tr/2);   %距離向LFM訊號
p0_f=fftshift(fft(fftshift(p0_t)));     % 距離向LFM參考訊號的快速傅立葉變換
%{
  概念普及：Matlab fftshift 詳解

一. 實訊號情況
因為實訊號以fs為取樣速率的訊號在 fs/2 處混疊，所以實訊號fft的結果中前半部分對應
[0, fs/2],後半部分對應[ -fs/2, 0]
1）實訊號fft的結果前半部分對應[0, fs/2]是正頻率的結果,後半部分對應[ -fs/2, 0]是
   負頻率的結果。大於fs/2的部分的頻譜實際上是實訊號的負頻率加fs的結果。故要得到
   正確的結果，只需將視在頻率減去fs即可得到頻譜對應的真實負頻率
2）如果要讓實訊號fft的結果與[-fs/2, fs/2]對應，則要fft後fftshift一下即可，
   fftshift的操作是將fft結果以fs/2為中心左右互換
3）如果實訊號fft的繪圖頻率f從[-fs/2, fs/2]，並且沒有fftshift，則fft正頻譜對應
   f在[0, fs/2]的結果將混疊到(f - fs/2)的位置;fft負頻譜對應f在[-fs/2, 0]的結果
   混疊到 f + fs - fs/2 的位置，注意這裡f為負值，也就是說此種情況下fft負頻譜對
    應的視在頻率減去fs/2即可得到頻譜對應的真實負頻率
 
二. 覆信號情況
1）覆信號沒有負頻率，以fs為取樣速率的訊號，fft的頻譜結果是從[0, fs]的。
2）在 f > fs/2 時，對覆信號的fft結果進行fftshift會產生頻率混疊(將下面的示例2中
   的頻率從f=15改為f=85可以驗證f=85的譜線在fftshift後跑到 f = -15 = 85 - fs =
   85 - 100的位置了)，所以覆信號也一般要求 f <= fs/2
3）在對雷達的慢時間維(覆信號)進行fft後，由於要用doppler = ((0:LFFT-1)/LFFT 
   - 0.5)*PRF; 計算多普勒頻率，所以對該慢時間訊號fft後要fftshift下，以便和正確
   的頻率單元相對應。注意多普勒頻率fd < = PRF/2 時才測的準！
%}

%--------------------距離向壓縮---------------
s_uf=fftshift(fft(fftshift(s_ut)));                 %距離向FFT
%{
  概念普及： fftshift
    作用：將零頻點移到頻譜的中間
    用法：
    Y=fftshift(X)
    Y=fftshift(X，dim)
    描述：fftshift移動零頻點到頻譜中間，重新排列fft,fft2和fftn的輸出結果。將
         零頻點放到頻譜的中間對於觀察傅立葉變換是有用的。
%}
src_uf=s_uf.*(conj(p0_f).'*ones(1,Na));             %距離壓縮
src_ut=fftshift(ifft(fftshift(src_uf)));      % IFFT後得到距離壓縮後的訊號 

%--------------------方位向壓縮--------------
src_fut=fftshift(fft(fftshift(src_ut).')).';        %距離多普勒域
%%二次距離壓縮,距離遷移校正原理模擬
src_fuf=fftshift(fft(fftshift(src_uf).')).';        %距離壓縮後的二維頻譜
F=f'*ones(1,Na);%擴充為矩陣
FU=ones(Nr,1)*fu;
p0_2f=exp(j*pi/fc^2/fdr*(FU.*F).^2+j*pi*fdc^2/fc/fdr*F-j*pi/fc/fdr*FU.^2.*F);
s2rc_fuf=src_fuf.*p0_2f;
s2rc_fut=fftshift(ifft(fftshift(s2rc_fuf)));%距離多普勒域
%%

%方位壓縮
p0_2fu=exp(j*pi/fdr*(FU-fdc).^2);%方位向壓縮因子
s2rcac_fut=s2rc_fut.*p0_2fu;%方位壓縮
s2rcac_fuf=fftshift(fft(fftshift(s2rcac_fut)));%距離方位壓縮後的二維頻譜
s2rcac_ut=fftshift(ifft(fftshift(s2rcac_fut).')).';%方位向IFFT
%%

%%
%%畫圖顯示結果
subplot(221)
G=20*log10(abs(s_ut)+1e-6);
gm=max(max(G));
gn=gm-40;%顯示動態範圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(x,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(a)原始訊號')

subplot(222)
G=20*log10(abs(src_fut)+1e-6);
gm=max(max(G));
gn=gm-40;%顯示動態範圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(fu,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(b)距離多普勒域頻譜')

subplot(223)
G=20*log10(abs(s2rc_fut)+1e-6);
gm=max(max(G));
gn=gm-40;%顯示動態範圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(fu,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(c)RMC後的RD域頻譜')

感謝各位讀者的支援，作者會努力提高部落格水準，歡迎讀者對錯誤或有待改進的地方提出建議和意見 .