點目標(biāo)仿真_SAR

1、合成孔徑雷達(SAR)的點目標(biāo)仿真一 SAR原理簡介合成孔徑雷達(Synthetic Aperture Radar ，簡稱SAR)是一種高分辨率成像雷達技術(shù)。它利用脈沖壓縮技術(shù)獲得高的距離向分辨率，利用合成孔徑原理獲得高的方位向分辨率，從而獲得大面積高分辨率雷達圖像。SAR回波信號經(jīng)距離向脈沖壓縮后，雷達的距離分辨率由雷達發(fā)射信號帶寬決定：，式中表示雷達的距離分辨率，表示雷達發(fā)射信號帶寬，表示光速。同樣，SAR回波信號經(jīng)方位向合成孔徑后，雷達的方位分辨率由雷達方位向的多譜勒帶寬決定：，式中表示雷達的方位分辨率，表示雷達方位向多譜勒帶寬，表示方位向SAR平臺速度。二 SAR的成像模式和空間幾何關(guān)

2、系根據(jù)SAR波束照射的方式，SAR的典型成像模式有Stripmap(條帶式)，Spotlight(聚束式)和Scan(掃描模式)，如圖2.1。條帶式成像是最早研究的成像模式，也是低分辨率成像最簡單最有效的方式；聚束式成像是在一次飛行中，通過不同的視角對同一區(qū)域成像，因而能獲得較高的分辨率；掃描模式成像較少使用，它的信號處理最復(fù)雜。 圖2.1：SAR典型的成像模式這里分析SAR點目標(biāo)回波時，只討論正側(cè)式Stripmap SAR，正側(cè)式表示SAR波束中心和SAR平臺運動方向垂直，如圖2.2，選取直角坐標(biāo)系XYZ為參考坐標(biāo)系，XOY平面為地平面；SAR平臺距地平面高h(yuǎn)，沿X軸正向以速度V勻速飛行；P

3、點為SAR平臺的位置矢量，設(shè)其坐標(biāo)為(x,y,z)； T點為目標(biāo)的位置矢量，設(shè)其坐標(biāo)為；由幾何關(guān)系，目標(biāo)與SAR平臺的斜距為： (2.1)由圖可知：；令，其中為平臺速度，s為慢時間變量（slow time），假設(shè)，其中表示SAR平臺的x坐標(biāo)為的時刻；再令，表示目標(biāo)與SAR的垂直斜距，重寫2.1式為： (2.2) 就表示任意時刻時，目標(biāo)與雷達的斜距。一般情況下，于是2.2式可近似寫為： (2.3)可見，斜距是的函數(shù)，不同的目標(biāo)，也不一樣，但當(dāng)目標(biāo)距SAR較遠(yuǎn)時，在觀測帶內(nèi)，可近似認(rèn)為不變，即。 圖2.2：空間幾何關(guān)系 (a)正視圖 (b)側(cè)視圖圖2.2(a)中，表示合成孔徑長度，它和合成孔徑時間

4、的關(guān)系是。(b)中，為雷達天線半功率點波束角，為波束軸線與Z軸的夾角，即波束視角，為近距點距離，為遠(yuǎn)距點距離，W為測繪帶寬度，它們的關(guān)系為： (2.4)三 SAR的回波信號模型SAR在運動過程中，以一定的PRT(Pulse Repitition Time,脈沖重復(fù)周期)發(fā)射和接收脈沖，天線波束照射到地面上近似為一矩形區(qū)域，如圖2.2（a），區(qū)域內(nèi)各散射元（點）對入射波后向散射，這樣，發(fā)射脈沖經(jīng)目標(biāo)和天線方向圖的調(diào)制，攜帶目標(biāo)和環(huán)境信息形成SAR回波。從時域來看，發(fā)射和接收的信號都是一時間序列。 圖3.1：SAR發(fā)射和接收信號圖3.1表示SAR發(fā)射和接收信號的時域序列。發(fā)射序列中，為chirp信

5、號持續(xù)時間，下標(biāo)表示距離向(Range)；PRT為脈沖重復(fù)周期；接收序列中，表示發(fā)射第個脈沖時，目標(biāo)回波相對于發(fā)射序列的延時；陰影部分表示雷達接收機采樣波門，采樣波門的寬度要保證能罩住測繪帶內(nèi)所有目標(biāo)的回波。雷達發(fā)射序列的數(shù)學(xué)表達式為： (3.1)式中，表示矩形信號，為距離向chirp信號的調(diào)頻斜率，為載頻。雷達回波信號由發(fā)射信號波形，天線方向圖，斜距，目標(biāo)RCS，環(huán)境等因素共同決定，若不考慮環(huán)境因素，則單點目標(biāo)雷達回波信號可寫成： (3.2)式中，為點目標(biāo)的雷達散射截面,表示點目標(biāo)天線方向圖雙向幅度加權(quán)，表示載機發(fā)射第n個脈沖時，電磁波在雷達與目標(biāo)之間傳播的雙程時間，代入3.2式 (3.3)

6、3.3式就是單點目標(biāo)回波信號模型。其中，為chirp分量，它決定距離向分辨率，為doppler分量，它決定方位向分辨率。距離向變量遠(yuǎn)大于方位向變量t(典型相差量級)，于是一般可以假設(shè)SAR滿足“停走?！蹦Ｊ?，即SAR在發(fā)射和接收一個脈沖信號中間，載機未發(fā)生運動。為了理論分析方便，稱為慢時間變量(slow time)，稱t為快時間變量(fast time)于是，一維回波信號可以寫成二維形式，正交解調(diào)去除載波后，單點目標(biāo)的回波可寫成： (3.3) 圖3.2：單點目標(biāo)回波二維分布示意圖 在方位向(慢時間域)是離散的，其中V是SAR的速度，是0時刻目標(biāo)在參考坐標(biāo)系中的x坐標(biāo)。為了作數(shù)字信號處理，在距離

7、向(快時間域)也要采樣，假設(shè)采樣周期為Tr,則，如圖3.2，方位向發(fā)射N個脈沖，距離向采樣得到M個樣值點，則SAR回波為一矩陣，K個理想點目標(biāo)的回波經(jīng)采樣后的表達式為： (3.4)上式用Matlab語言可表示為：%*%Generate the raw signal dataK=Ntarget; %number of targetsN=Nslow; %number of vector in slow-time domainM=Nfast; %number of vector in fast-time domainT=Ptarget; %position of targetsSrnm=zeros(

8、N,M);for k=1:1:K sigma=T(k,3); Dslow=sn*V-T(k,1); R=sqrt(Dslow.2+T(k,2)2+H2); tau=2*R/C; Dfast=ones(N,1)*tm-tau*ones(1,M); phase=pi*Kr*Dfast.2-(4*pi/lambda)*(R*ones(1,M);Srnm=Srnm+sigma*exp(j*phase).*(0Dfast&DfastTr).*(abs(Dslow)Lsar/2)*ones(1,M);end%*四 SAR的信號系統(tǒng)模型從信號與系統(tǒng)的角度看，SAR回波可看作目標(biāo)的散射特性通過一個二維線性系統(tǒng)

9、的輸出。點目標(biāo)的信號與系統(tǒng)模型如圖4.1：圖4.1：點目標(biāo)信號與系統(tǒng)模型模型的數(shù)學(xué)表達式為： (4.1)式中，表示點目標(biāo)的散射特性，表示等效系統(tǒng)，設(shè)為發(fā)射的chirp信號，則： (4.2)4.2式表明只在維是線性時不變(LTI)的，在維是時變的，相同的, 不同的，響應(yīng)不一樣。但通常情況下可近似認(rèn)為不變，即，這時，系統(tǒng)等效為一個二維LTI系統(tǒng)。五 點目標(biāo)SAR的成像處理算法仿真SAR的回波數(shù)據(jù)不具有直觀性，不經(jīng)處理人無法理解它，如圖5.1。從原理上講，SAR成像處理的過程是從回波數(shù)據(jù)中提取目標(biāo)區(qū)域散射系數(shù)的二維分布，本質(zhì)上是一個二維相關(guān)處理過程,因此最直接的處理方法是對回波進行二維匹配濾波，但其

10、運算量很大，再加上SAR的數(shù)據(jù)率本來就高，這使得實時處難于實現(xiàn)。通常，可以把二維過程分解成距離向和方位向兩個一維過程，Range-Dopper Algorithm（簡稱RD算法）就是采用這種思想的典型算法，這里也只討論RD算法。 圖5.1：SAR回波數(shù)據(jù) (a)未經(jīng)處理 (b)處理后RD算法通過距離遷移(Range Migration)矯正，消除距離和方位之間的耦合。在滿足聚焦深度的前提小，將成像處理分解成兩個一維的LTI系統(tǒng)進行相關(guān)處理，并采用頻域快速相關(guān)算法提高了速度。RD算法已非常成熟，并成為衡量其它算法優(yōu)劣的標(biāo)準(zhǔn)。RD算法典型的數(shù)字處理流程如圖5.2。 圖5.2: Range-Dopp

11、er 算法流程 預(yù)處理這是對SAR回波處理的第一步，一般在SAR平臺(衛(wèi)星，飛機)上實時處理，包括解調(diào)和數(shù)字化。雷達信號的載頻較高(GHz)，不宜直接采樣數(shù)字化處理，常常通過正交解調(diào)方式解調(diào)出基帶信號，再對基帶信號(MHz)采樣數(shù)字化，然后存儲或傳到地面做進一步處理。采樣后的數(shù)據(jù)常采用矩陣形式存儲，假設(shè)方位向發(fā)射(采樣)N個脈沖，距離向采樣得到M每個采樣值(圖3.2)，則待處理數(shù)據(jù)是一個的矩陣，如圖5.3。實際處理時，要在方位向上加窗截斷，因此，在方位向上的開始和結(jié)束的一段數(shù)據(jù)(圖中影陰區(qū)所示)是不充分的，對應(yīng)的長度均為，表示SAR的合成孔徑長度。仿真時，這個數(shù)據(jù)陣是程序根據(jù)3.4式產(chǎn)生的。圖

12、5.3：待處理數(shù)據(jù) 距離壓縮距離向信號是典型的Chirp信號，相關(guān)算法是在頻域利用FFT進行的。Matlab語句為：Refr=exp(j*pi*Kr*tr.2).*(0tr&trTr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr); 和分別是Chirp信號的調(diào)頻斜率和脈沖持續(xù)時間，Refr表示參考信號，fty,ifty是對矩陣的行(對應(yīng)距離向)進行FFT和IFFT運算的子程序。例如，fty的代碼為： %FFT in row of matrixfunction fs=fty(s);fs=fftshift(fft(fftshift(s.).; 距離遷移矯

13、正距離遷移是SAR信號處理中必然出現(xiàn)的現(xiàn)象，它的大小隨系統(tǒng)參數(shù)不同而變化，并不總需要補償。點目標(biāo)仿真時，可以先不考慮。 方位壓縮方位向的處理是SAR成像處理算法最核心的部分。正側(cè)式點目標(biāo)(圖2.2)情況下，回波經(jīng)距離壓縮后在方位向也是一Chirp信號，因此其壓縮處理同距離壓縮處理類似，只是壓縮因子不同。仿真中，調(diào)頻斜率已知，因此不需要進行Doppler參數(shù)估計。 SAR參數(shù)SAR平臺：水平速度V=100m/s高度H=5000m天線等效孔徑D=4mSAR平臺與測繪帶的垂直斜距R0=11180m(計算結(jié)果)發(fā)射信號： 載波頻率=1GHz Chirp信號持續(xù)時間=5usChirp信號調(diào)頻帶寬=30M

14、Hz Chirp信號調(diào)頻斜率=(計算結(jié)果) 脈沖重復(fù)頻率PRF=57.6Hz(計算結(jié)果)Doppler調(diào)頻帶寬=50Hz(計算結(jié)果)Doppler調(diào)頻斜率=5.96(計算結(jié)果)分辨率： 距離向分辨率DY=5m 方位向分辨率DX=2m目標(biāo)位置: 距離向Y=Yc-Y0,Yc+Y0=9500,10500m 方位向X=Xmin,Xmax=0,50m 目標(biāo)個數(shù)Ntarget=3 目標(biāo)位置矩陣: 格式 x坐標(biāo),y坐標(biāo),目標(biāo)散射系數(shù)Ptarget=Xmin,Yc,1 Xmin,Yc+10*DY,1 Xmin+20*DX,Yc+50*DY,1stripmapSAR.m程序(見附錄)實現(xiàn)了仿真功能，圖5.4到圖

15、5.7為仿真結(jié)果。運行程序，在Command Window 中列出了仿真的參數(shù)：Parameters:Sampling Rate in fast-time domain 3.0996Sampling Number in fast-time domain 1024Sampling Rate in slow-time domain 1.1525Sampling Number in slow-time domain 512Range Resolution 5Cross-range Resolution 2SAR integration length 838.5255Position of targe

16、ts 0 10000 1 0 10050 1 40 10250 1當(dāng)然，這些參數(shù)可以改變以得到不同的結(jié)果，但值得注意的是，采樣點數(shù)不宜過大，否則數(shù)據(jù)量過大將導(dǎo)致程序運行時間過長，甚至計算機因內(nèi)存耗盡而死機。本例采用的是5121024個點。 圖5.4：SAR的點目標(biāo)仿真結(jié)果 圖5.5:兩點目標(biāo)的回波仿真3D圖 圖5.6：兩點目標(biāo)距離向壓縮后的3D圖圖5.7：兩點目標(biāo)距離向和方位向壓縮后的3D圖 圖5.8：兩點目標(biāo)壓縮后的3dB等高線圖附錄：SAR的點目標(biāo)仿真Matlab程序主程序：stripmapSAR.m%=clear;clc;close all;%=%Parameter-constantC=

17、3e8; %propagation speed%Parameter-radar characteristicsFc=1e9; %carrier frequency 1GHzlambda=C/Fc; %wavelength %Parameter-target areaXmin=0; %target area in azimuth is withinXmin,XmaxXmax=50; Yc=10000; %center of imaged areaY0=500; %target area in range is withinYc-Y0,Yc+Y0 %imaged width 2*Y0%Parame

18、ter-orbital informationV=100; %SAR velosity 100 m/sH=5000; %height 5000 mR0=sqrt(Yc2+H2);%Parameter-antennaD=4; %antenna length in azimuth directionLsar=lambda*R0/D; %SAR integration lengthTsar=Lsar/V; %SAR integration time%Parameter-slow-time domainKa=-2*V2/lambda/R0; %doppler frequency modulation

19、rateBa=abs(Ka*Tsar); %doppler frequency modulation bandwidthPRF=Ba; %pulse repitition frequencyPRT=1/PRF; %pulse repitition timeds=PRT; %sample spacing in slow-time domainNslow=ceil(Xmax-Xmin+Lsar)/V/ds); %sample number in slow-time domainNslow=2nextpow2(Nslow); %for fftsn=linspace(Xmin-Lsar/2)/V,(X

20、max+Lsar/2)/V,Nslow);%discrete time array in slow-time domainPRT=(Xmax-Xmin+Lsar)/V/Nslow; %refreshPRF=1/PRT;ds=PRT;%Parameter-fast-time domainTr=5e-6; %pulse duration 10usBr=30e6; %chirp frequency modulation bandwidth 30MHzKr=Br/Tr; %chirp slopeFsr=3*Br; %sampling frequency in fast-time domaindt=1/

21、Fsr; %sample spacing in fast-time domainRmin=sqrt(Yc-Y0)2+H2);Rmax=sqrt(Yc+Y0)2+H2+(Lsar/2)2);Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%sample number in fast-time domainNfast=2nextpow2(Nfast); %for ffttm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast); %discrete time array in fast-time domaindt=(2*Rmax/C+Tr-2*Rmin

22、/C)/Nfast; %refreshFsr=1/dt;%Parameter-resolutionDY=C/2/Br; %range resolutionDX=D/2; %cross-range resolution%Parameter-point targetsNtarget=2; %number of targets%format x, y, reflectivityPtarget=Xmin,Yc,1 %position of targets Xmin,Yc+10*DY,1 Xmin+20*DX,Yc+50*DY,1; disp(Parameters:)disp(Sampling Rate

23、 in fast-time domain);disp(Fsr/Br)disp(Sampling Number in fast-time domain);disp(Nfast)disp(Sampling Rate in slow-time domain);disp(PRF/Ba)disp(Sampling Number in slow-time domain);disp(Nslow)disp(Range Resolution);disp(DY)disp(Cross-range Resolution);disp(DX) disp(SAR integration length);disp(Lsar)

24、 disp(Position of targets);disp(Ptarget)%=%Generate the raw signal dataK=Ntarget; %number of targetsN=Nslow; %number of vector in slow-time domainM=Nfast; %number of vector in fast-time domainT=Ptarget; %position of targetsSrnm=zeros(N,M);for k=1:1:K sigma=T(k,3); Dslow=sn*V-T(k,1); R=sqrt(Dslow.2+T

25、(k,2)2+H2); tau=2*R/C; Dfast=ones(N,1)*tm-tau*ones(1,M); phase=pi*Kr*Dfast.2-(4*pi/lambda)*(R*ones(1,M); Srnm=Srnm+sigma*exp(j*phase).*(0Dfast&DfastTr).*(abs(Dslow)Lsar/2)*ones(1,M);end%=%Range compressiontr=tm-2*Rmin/C;Refr=exp(j*pi*Kr*tr.2).*(0tr&trTr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr);

26、Gr=abs(Sr);%Azimuth compressionta=sn-Xmin/V;Refa=exp(j*pi*Ka*ta.2).*(abs(ta)Tsar/2);Sa=iftx(ftx(Sr).*(conj(ftx(Refa).*ones(1,M);Ga=abs(Sa);%=%graw the intensity image of signalcolormap(gray);figure(1)subplot(211);row=tm*C/2-2008;col=sn*V-26;imagesc(row,col,255-Gr); %intensity image of Sraxis(Yc-Y0,Y

27、c+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel(rightarrowitRange in meters),ylabel(itAzimuth in metersleftarrow),title(Stripmap SAR after range compression),subplot(212);imagesc(row,col,255-Ga); %intensity image of Saaxis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel(rightarrowitRange in meters),ylabel(itAzimuth in metersleftarrow),title(Stripmap SAR after range and azimuth compression),%=%draw 3D picturefigure(2)waterfall(real(Srnm(200:205),:);axis tight