合成孔径雷达(SAR)的点目标仿真(附件带代码程序).doc

合成孔径雷达(SAR)的点目标仿真(附件带代码程序)
-
合成孔径雷达(SAR)的点目标仿真 { %]imf|g.  
一．    SAR原理简介 kmt1vV.9  
合成孔径雷达(Synthetic Aperture Radar ，简称SAR)是一种高分辨率成像雷达技术。它利用脉冲压缩技术获得高的距离向分辨率，利用合成孔径原理获得高的方位向分辨率，从而获得大面积高分辨率雷达图像。 NIbK3`1  
SAR回波信号经距离向脉冲压缩后，雷达的距离分辨率由雷达发射信号带宽决定： ，式中 表示雷达的距离分辨率， 表示雷达发射信号带宽， 表示光速。同样，SAR回波信号经方位向合成孔径后，雷达的方位分辨率由雷达方位向的多谱勒带宽决定： ，式中 表示雷达的方位分辨率， 表示雷达方位向多谱勒带宽， 表示方位向SAR平台速度。 (1
GU  
二．    SAR的成像模式和空间几何关系 y5+%8#3  
根据SAR波束照射的方式，SAR的典型成像模式有Stripmap(条带式)，Spotlight(聚束式)和Scan(扫描模式)，如图2.1。条带式成像是最早研究的成像模式，也是低分辨率成像最简单最有效的方式；聚束式成像是在一次飞行中，通过不同的视角对同一区域成像，因而能获得较高的分辨率；扫描模式成像较少使用，它的信号处理最复杂。 btq 4diW  
*x;4::'Jn  
                  图2.1：SAR典型的成像模式 u Z-ZZE C  
这里分析SAR点目标回波时，只讨论正侧式Stripmap SAR，正侧式表示SAR波束中心和SAR平台运动方向垂直，如图2.2，选取直角坐标系XYZ为参考坐标系，XOY平面为地平面；SAR平台距地平面高h，沿X轴正向以速度V匀速飞行；P点为SAR平台的位置矢量，设其坐标为(x,y,z)； T点为目标的位置矢量，设其坐标为 ；由几何关系，目标与SAR平台的斜距为： m}=E$zPbO  
                      (2.1) b"N!#&O
]  
由图可知： ；令 ，其中 为平台速度，s为慢时间变量（slow time），假设 ，其中 表示SAR平台的x坐标为 的时刻；再令 ， 表示目标与SAR的垂直斜距，重写2.1式为： }(ORh2Ri  
                                       (2.2)     ?&~q^t?u  
就表示任意时刻 时，目标与雷达的斜距。一般情况下， ，于是2.2式可近似写为： @&WHX#  
                 (2.3) }shxEsq  
可见，斜距是 的函数，不同的目标， 也不一样，但当目标距SAR较远时，在观测带内，可近似认为 不变，即 。           1H/I-

 
/+SLq`'u)  
               图2.2：空间几何关系 (a)正视图 (b)侧视图 +W*~=*h|  
图2.2(a)中， 表示合成孔径长度，它和合成孔径时间 的关系是 。(b)中， 为雷达天线半功率点波束角， 为波束轴线与Z轴的夹角，即波束视角， 为近距点距离， 为远距点距离，W为测绘带宽度，它们的关系为： 'Uf?-t*LT@  
                                                     (2.4) #H>{>0q  
SXOAa1表示SAR发射和接收信号的时域序列。发射序列中， 为chirp信号持续时间，下标 表示距离向(Range)；PRT为脉冲重复周期；接收序列中， 表示发射第 个脉冲时，目标回波相对于发射序列的延时；阴影部分表示雷达接收机采样波门，采样波门的宽度要保证能罩住测绘带内所有目标的回波。 l:!L+t*}6  
雷达发射序列的数学表达式为： E6A"Xo  
                           (3.1) oY18a*_>M1  
式中， 表示矩形信号， 为距离向chirp信号的调频斜率， 为载频。 Yv;s3>r  
雷达回波信号由发射信号波形，天线方向图，斜距，目标RCS，环境等因素共同决定，若不考虑环境因素，则单点目标雷达回波信号可写成： ?0VR2Yb${b  
                        (3.2) vjGQ!xF  
式中， 为点目标的雷达散射截面, 表示点目标天线方向图双向幅度加权， 表示载机发射第n个脉冲时，电磁波在雷达与目标之间传播的双程时间， ，代入3.2式 I#M3cI!X?  
            (3.3) W p=3heCa6  
3.3式就是单点目标回波信号模型。其中， 为chirp分量，它决定距离向分辨率， 为doppler分量，它决定方位向分辨率。 O  
距离向变量 远大于方位向变量t(典型相差 量级)，于是一般可以假设SAR满足“停－走－停”模式，即SAR在发射和接收一个脉冲信号中间，载机未发生运动。为了理论分析方便，称 为慢时间变量(slow time)，称t为快时间变量(fast time)于是，一维回波信号可以写成二维形式，正交解调去除载波后，单点目标的回波可写成： m|@H`=`d  
          (3.3) nVb@sI{{k  
SA x9cjj +  
_!\d?]Ya  
          图3.2：单点目标回波二维分布示意图 X W)TI  
   在方位向(慢时间域)是离散的， ，其中V是SAR的速度， 是0时刻目标在参考坐标系中的x坐标。为了作数字信号处理，在距离向(快时间域)也要采样，假设采样周期为Tr,则 ，如图3.2，方位向发射N个脉冲，距离向采样得到M个样值点，则SAR回波为一 矩阵，K个理想点目标的回波经采样后的表达式为： %@/^ UE:  
vT'Bs;QR  
          (3.4) 0+1wi4wy/  
VM"cpC_8  
 gV_/
t+jI  
%%Generate the raw signal data ^NW[)Dq1<  
K=Ntarget;                             %number of targets @=Kq99=\U  
N=Nslow;                              %number of vector in slow-time domain ^4y]7 p  
M=Nfast;                               %number of vector in fast-time domain ;t`  ?|  
T=Ptarget;                              %position of targets (J hX:1  
Srnm=zeros(N,M); =BsV`p7rU  
for k=1:1:K ^ k^y|\UtZ  
    sigma=T(k,3); adON&<  
    Dslow=sn*V-T(k,1); Z O}Og&%  
    R=sqrt(Dslow.^2+T(k,2)^2+H^2); 0bNvmZ$  
    tau=2*R/C; 'nJF:+30ZH  
    Dfast=ones(N,1)*tm-tau'*ones(1,M); HDyf]2N*N  
    phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M)); mgxIxusR  
Srnm=Srnm+sigma*exp(j*phase).*(0 8$Q`wRt(%  
61\u{@
o$  
四．    SAR的信号系统模型 l,wlxh$}(  
b`0tfXzS5  
从信号与系统的角度看，SAR回波可看作目标的散射特性通过一个二维线性系统的输出。点目标的信号与系统模型如图4.1： ! VwU=5  
.^YxhUH,G  
图4.1：点目标信号与系统模型 &|SWy 2 N  
模型的数学表达式为： h!"2Ux3!x  
                   (4.1) IQQv+af5  
式中， 表示点目标的散射特性， 表示等效系统，设 为发射的chirp信号，则： 6fwNlC/9  
                             (4.2) C|{Sj`,XG  
4.2式表明 只在 维是线性时不变(LTI)的，在 维是时变的，相同的 , 不同的 ，响应 不一样。但通常情况下可近似认为 不变，即 ，这时，系统等效为一个二维LTI系统。 37bM e@W  
D{9a'0J  
OqGp|`  
五．    点目标SAR的成像处理算法仿真 |T#cq!  
SAR的回波数据不具有直观性，不经处理人无法理解它，如图5.1。从原理上讲，SAR成像处理的过程是从回波数据中提取目标区域散射系数的二维分布，本质上是一个二维相关处理过程,因此最直接的处理方法是对回波进行二维匹配滤波，但其运算量很大，再加上SAR的数据率本来就高，这使得实时处难于实现。通常，可以把二维过程分解成距离向和方位向两个一维过程，Range-Dopper Algorithm（简称RD算法）就是采用这种思想的典型算法，这里也只讨论RD算法。 9#E *o~1  
UBo0c?,4  
                 图5.1：SAR回波数据 (a)未经处理 (b)处理后 ITPE2x  
RD算法通过距离迁移(Range Migration)矫正，消除距离和方位之间的耦合。在满足聚焦深度的前提小，将成像处理分解成两个一维的LTI系统进行相关处理，并采用频域快速相关算法提高了速度。RD算法已非常成熟，并成为衡量其它算法优劣的标准。RD算法典型的数字处理流程如图5.2。 X)g X9DA  
(}#&HE
<  
                  图5.2: Range-Dopper 算法流程 =^"Sx??V  
●    预处理 lM-\:Q!
  
这是对SAR回波处理的第一步，一般在SAR平台(卫星，飞机)上实时处理，包括解调和数字化。雷达信号的载频较高(～GHz)，不宜直接采样数字化处理，常常通过正交解调方式解调出基带信号，再对基带信号(～MHz)采样数字化，然后存储或传到地面做进一步处理。采样后的数据常采用矩阵形式存储，假设方位向发射(采样)N个脉冲，距离向采样得到M每个采样值(图3.2)，则待处理数据是一个 的矩阵，如图5.3。实际处理时，要在方位向上加窗截断，因此，在方位向上的开始和结束的一段数据(图中影阴区所示)是不充分的，对应的长度均为 ， 表示SAR的合成孔径长度。仿真时，这个数据阵是程序根据3.4式产生的。 \Z*:l(  
L03I:IJ  
图5.3：待处理数据 F6{g{ B  
●    距离压缩 He_O+[sc  
距离向信号是典型的Chirp信号，相关算法是在频域利用FFT进行的。Matlab语句为： #v c+;`X  
Refr=exp(j*pi*Kr*tr.^2).*(0  
Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr)))); @Y ~gd K  
    和 分别是Chirp信号的调频斜率和脉冲持续时间，Refr表示参考信号，fty,ifty是对矩阵的行(对应距离向)进行FFT和IFFT运算的子程序。例如，fty的代码为： t%@iF U;}  
    %%FFT in row of matrix S iNgV\('U  
function fs=fty(s); 4'faE="1)S  
fs=fftshift(fft(fftshift(s.'))).'; e`Yj}i*bx]  
●    距离迁移矫正 $`xpn#l z  
距离迁移是SAR信号处理中必然出现的现象，它的大小随系统参数不同而变化，并不总需要补偿。点目标仿真时，可以先不考虑。 m
bA
zn  
●    方位压缩 jS,Pu%fR  
方位向的处理是SAR成像处理算法最核心的部分。正侧式点目标(图2.2)情况下，回波经距离压缩后在方位向也是一Chirp信号，因此其压缩处理同距离压缩处理类似，只是压缩因子不同。仿真中，调频斜率 已知，因此不需要进行Doppler参数估计。 clDn=k<  
●    SAR参数 
sta/i?n  
SAR平台： y0O(n/  
水平速度V=100m/s >MJ %6A>  
高度H=5000m \Vv)(/q{  
天线等效孔径D=4m }[%F  
SAR平台与测绘带的垂直斜距R0=11180m(计算结果) [\F:NLjiUy  
发射信号： (muJ-~CJk  
   载波频率 =1GHz :WH{wm|  
   Chirp信号持续时间 =5us !:d\A  
Chirp信号调频带宽 =30MHz
EAY+#>L*  
   Chirp信号调频斜率 = (计算结果) +VIEDV+   
   脉冲重复频率PRF=57.6Hz(计算结果) i Q3wi  
Doppler调频带宽 =50Hz(计算结果) qZQm*q(jM  
Doppler调频斜率 =－5.96(计算结果) h*4wi.-  
分辨率： wQ-BY"cK\  
   距离向分辨率DY=5m ABE EJQ  
   方位向分辨率DX=2m r
I)&.5^  
目标位置: % [,^2s  
   距离向Y=[Yc-Y0,Yc+Y0]=[9500,10500]m _#stripmapSAR.m程序(见附录)实现了仿真功能，图5.4到图5.7为仿真结果。 Io;26F""  
运行程序，在Command Window 中列出了仿真的参数： g.EKdvY"%H  
Parameters: E<_+Tc  
Sampling Rate in fast-time domain )3_g&&  
    3.0996 fI(H : N  
u`X}AKC  
Sampling Number in fast-time domain Eh;~y*k\  
    1024 |7T!rnr  
4~r=[|(aY  
Sampling Rate in slow-time domain \+w -{"u$  
    1.1525 SY 2B\TV  
dTP$7nfe
 
Sampling Number in slow-time domain QSy#k~  
   512 (5N&b h`E  
e~we YGK  
Range Resolution oZ~M`
yOz.  
     5 /$x6//0If  
NLF{W|X  
Cross-range Resolution U=%S6uL\bx  
     2 ^IM;D)X&:  
rReZ$U  
SAR integration length :/Nz' n  
  838.5255 82V xk  
diqG8KaK  
Position of targets =J?<M?ugf  
           0       10000           1 8'#%7+ "=!  
           0       10050           1 c#`Z[  
          40       10250           1 IEeh)aj[  
当然，这些参数可以改变以得到不同的结果，但值得注意的是，采样点数不宜过大，否则数据量过大将导致程序运行时间过长，甚至计算机因内存耗尽而死机。本例采用的是512 1024个点。 |2c'0Ibu  
   F/RV{} 17E  
                         图5.4：SAR的点目标仿真结果 (g~&$&pa  
iAWPE`u4  
                         图5.5:两点目标的回波仿真3D图 NY1olnI  
0Km{fZYq7;  
                        图5.6：两点目标距离向压缩后的3D图 Z10#6v  
ga~rllm;i  
(h|l$OL/  
图5.7：两点目标距离向和方位向压缩后的3D图 pt|u?T_+  
H-0A&oG  
                图5.8：两点目标压缩后的3dB等高线图 nvsuF)%9hZ  
附录：SAR的点目标仿真Matlab程序 WQ=C5^u  
主程序：stripmapSAR.m zn>*^h0B  
%%======================================================== Z]":xl\7  
clear;clc;close all; +!dIEt).U  
%%==== *w59BO&M4  
%%Parameter--constant TpuN[Y  
C=3e8;                           %propagation speed 2|kx:^D p  
%%Parameter--radar characteristics .& B_\*  
Fc=1e9;                          %carrier frequency 1GHz Pfe&wA't  
lambda=C/Fc;                     %wavelength ;:hyW,J  
%%Parameter--target area `\(F ax  
Xmin=0;                          %target area in azimuth is within[Xmin,Xmax] ;Q:^|Fw!F  
Xmax=50;                   +k
d88Fx  
Yc=10000;                        %center of imaged area G5Y 8]N  
Y0=500;                          %target area in range is within[Yc-Y0,Yc+Y0] ~}z{RE($v  
                                 %imaged width 2*Y0 CY? ]o4IV  
%%Parameter--orbital information ~2beVQ(U  
V=100;                           %SAR velosity 100 m/s @c,}\"(  
H=5000;                          %height 5000 m =4/LixsV|  
R0=sqrt(Yc^2+H^2); @ZK#Y){  
%%Parameter--antenna {  KE[8n  
D=4;                            %antenna length in azimuth direction lNs 'jaD  
Lsar=lambda*R0/D;               %SAR integration length p{"p["T+L`  
%%Parameter--slow-time domain yK+76\} I  
Ka=-2*V^2/lambda/R0;            %doppler frequency modulation rate rteViq+|.  
Ba=abs(Ka*Tsar);                %doppler frequency modulation bandwidth )NGBA."t  
PRF=Ba;                        %pulse repitition frequency J*F-tRuEw  
PRT=1/PRF;                     %pulse repitition time [Qa0uM#SU  
ds=PRT;                         %sample spacing in slow-time domain j wlmWO6  
Nslow=ceil((Xmax-Xmin+Lsar)/V/ds); %sample number in slow-time domain AMA :hQ  
Nslow=2^nextpow2(Nslow);          %for fft c5>&~^~>Tx  
sn=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%discrete time array in slow-time domain %?J\P@  
PRT=(Xmax-Xmin+Lsar)/V/Nslow;    %refresh - q[T0^e S  
PRF=1/PRT; #dE#w#=r  
ds=PRT; lSC3m=4g  
%%Parameter--fast-time domain 64z9Yr@  
Tr=5e-6;                         %pulse duration 10us SE)_
5|k*  
Br=30e6;                        %chirp frequency modulation bandwidth 30MHz H@bra~k-  
Kr=Br/Tr;                        %chirp slope @y[Zr 6\z