sar合成孔徑雷達圖像點目標仿真報告附matlab代碼

1、SAR圖像點目標仿真報告徐一凡1SAR原理簡介合成孔徑雷達(SyntheticApertureRadar，簡稱SAR)是一種高分辨率成像雷達技術。它利用脈沖壓縮技術獲得高的距離向分辨率，利用合成孔徑原理獲得高的方位向分辨率，從而獲得大面積高分辨率雷達圖像。SAR回波信號經距離向脈沖壓縮后，雷達的距離分辨率由雷達發(fā)射信號帶寬決定：CP二，式中P表示雷達的距離分辨率，B表示雷達發(fā)射信號帶寬，C表示光速。同r2Brrr樣，SAR回波信號經方位向合成孔徑后，雷達的方位分辨率由雷達方位向的多譜勒帶寬決定：P二?，式中P表示雷達的方位分辨率，B表示雷達方位向多譜勒帶寬，v表示aBaaaa方位向SAR平臺速

2、度。在小斜視角的情況下，方位分辨率近似表示為p二D，其中D為a2方位向合成孔徑的長度。2SAR的幾何關系雷達位置和波束在地面覆蓋區(qū)域的簡單幾何模型如圖1所示。此次仿真考慮的是正側視的條帶式仿真，也就是說傾斜角為零，SAR波束中心和SAR平臺運動方向垂直的情況。圖1雷達數(shù)據(jù)獲取的幾何關系建立坐標系XYZ如圖2所示，其中XOY平面為地平面；SAR平臺距地平面高H,以速度V沿X軸正向勻速飛行；P點為SAR平臺的位置矢量，設其坐標為(x,y,z);T點為目標的位置矢量，設其坐標為(兀丁,打,zj；由幾何關系，目標與SAR平臺的斜距為：R二PT=J(x-x)2+(y-y)2+(z-z)2ITtt(1)由

3、圖可知：y=°，z=H,Zf=0;令x=vs-，其中v為平臺速度,s為慢時間變量(slowtime),假設需=vs，其中s表示SAR平臺的x坐標為需的時刻；再令廠=込2+晉，r表示目標與SAR的垂直斜距，重寫(1)式為：PT=R(s;r)=r2+v2-(s-s)2(2)R(s;r)就表示任意時刻s時，目標與雷達的斜距。一般情況下，v|s-s0|<<r，于是通過傅里葉技術展開，可將(2)式可近似寫為：v2R(s;r)=r2+v2-(s一s)2qr+(s一s)2(3)丫02r0可見，斜距是sWr的函數(shù)，不同的目標，r也不一樣，但當目標距SAR較遠時，在觀測帶內，可近似認為r不

4、變，即r=R0。圖2:空間幾何關系(a)正視圖(b)側視圖圖2(a)中，Lsa廠表示合成孔徑長度，它和合成孔徑時間Tsa廠的關系是Lsar=vTsar.(b)中，&為雷達天線半功率點波束角，°為波束軸線與Z軸的夾角，即波束視角，Rmin為近距點距離，Rmax為遠距點距離，w為測繪帶寬度，它們的關系為：Rmin二H-tg(0-Rmax二H-tg(0+辱；)(4)W二Rmax-Rmin3SAR的回波信號模型SAR在運動中以一定的周期(1/PR巧發(fā)射和接收信號，具體過程如圖3所示。發(fā)射機以工/的時間發(fā)射啁啾脈沖，然后切換天線開關接收回波信號。圖3雷達發(fā)射脈沖串的時序當雷達不處于發(fā)射

5、狀態(tài)時,它接收3反射回波。發(fā)射和接收回波的時間序列如圖4所示。在機載情況下，每個回波可以在脈沖發(fā)射間隔內直接接收到。但是在星載情況下，由于距離過大，某個脈沖的回波要經過610個脈沖間隔才能接收到。這里仿真為了方便，默認為機載情況。圖4脈沖雷達的發(fā)射與接收周期假設T為chirp信號持續(xù)時間，下標r表示距離向；PRF為重復頻率，PRT為重復周期，r等于1/PRF。接收序列中，工二2*號s;r)表示發(fā)射第i個脈沖時，目標回波相對于發(fā)射nC序列的延時。雷達的發(fā)射序列數(shù)學表達式為式(5):s(t)=遠p(t-n*PRT)n=s(5)P(t)=rect()ej冗Krt2ej2kftr式中，WCtH表示矩形

6、信號，Kr為距離向的chirp信號調頻率，f為載頻。雷達回波信號由發(fā)射信號波形，天線方向圖，斜距，目標RCS，環(huán)境等因素共同決定，若不考慮環(huán)境因素，則單點目標雷達回波信號可寫成式(6)所示：s(t)=£Qwp(t-n-PRT-T)(6)rnn其中，Q表示點目標的雷達散射截面，w表示點目標天線方向圖雙向幅度加權，T表n示載機發(fā)射第n個脈沖時，電磁波再次回到載機時的延時T二2*R('；廠),帶入式中得:nC8/nPRT2R(s;r)/Cs(t)二w-rect()-rTn=8rexpj兀K(tn-PRT2R(s;r)/C)2-r4兀exp-jR(s;r)-expj2兀f(tn-PR

7、TT)人cn式就是單點目標回波信號模型，其中，expjK(tn-PRT2R(s;r)/C)2是r4兀chirp分量，它決定距離向分辨率；exp-jR(s;r)為多普勒分量，它決定方位向分辨率。入對于任意一個脈沖，回波信號可表示為式(8)所示：s(t,s)=Aw(t2R(s;r)/C)w(ss)xexpj4兀fR(s;r)/C-r0rac'.0(8)expj兀K(t2R(s;r)/C)2r我們知道，由于R(s;r)隨慢時間s的變化而變化，所以計算機記錄到的回波數(shù)據(jù)存儲形式如圖5所示：圖5目標照射時間內，單個點目標回波能量在信號處理器的二維存儲器中的軌跡4距離徙動及校正根據(jù)圖2可知，在傾斜

8、角為零或很小的時候，目標與雷達的瞬時距離為R(s;r)，根據(jù)幾何關系可知，R(s;丫)=丫2+v2-(ss0)2，根據(jù)泰勒級數(shù)展開可得:R(s;r)=Jr2+v2-(ss)2沁r+(ss)2(9)v02r0q由式(9)可知，不同慢時間對應著不同的R(s;廠)，并且是一個雙曲線形式或者近似為一個二次形式。如圖5所示,同一目標的回波存儲在計算機里不在同一直線上，存在距離徙動。從而定義距離徙動量：v2AR(s,r)=(ss)2(10)2r0為了進行方位向的壓縮，方位向的回波數(shù)據(jù)必須在同一條直線上，也就是說必須校正距離徙動AR(s,r)。由式(10)可知，不同的最近距離r對應著不同的AR(s,r)，因

9、此在時域處理距離徙動會非常麻煩。因此，對方位向進行傅里葉變換，對距離向不進行變換，得到新的域。由于方位向的頻率即為多普勒頻率，所以這個新的域也稱為距離多普勒域。將斜距R寫成多普勒fa的函數(shù)，即R(f,r)。眾所周知，對最近距離為r的點目標P,a2V回波多普勒f是傾斜角0的函數(shù)，即/=vsin0，斜距R(f,r)=r/cos9，于是aa九aR(f,r)=r/cos9=r/、；1一sin292(11)a1(九)£r+8(V)2fa1a所以距離多普勒域中的我距離徙動為AR(f,r)=-(-)2rf2，可發(fā)現(xiàn)它不隨慢時間變換，a8Va同一最短距離r對應著相同大小的距離徙動。因此在距離多普勒域

10、對一個距離徙動校正就是對一組具有相同最短距離的點目標的距離徙動校正，這樣可以節(jié)省運算量。為了對距離徙動進行校正，需要得到距離徙動單元，即距離徙動體現(xiàn)在存儲單元中的移動數(shù)值，距離徙動單元可以表示為AR(f,r)/p，這個值通常為一個分數(shù)，由于存儲單元ar都是離散的，所以不同通過在存儲單元簡單的移動得到準確的值。為了得到準確的徙動校正值，通常需要進行插值運算。本仿真采用了兩種插值方法最近鄰點插值和sinc插值，下面分別進行介紹。最近鄰點插值法的優(yōu)點是簡單而快速，缺點是不夠精確。AR(f,r)/p=N+n，其中N為整數(shù)部分，arn為小數(shù)部分，整數(shù)部分徙動可以直接通過平移消除，對于小數(shù)部分則通過四舍五

11、入的方法變?yōu)?或者1,這樣就可以得到較為精確的插值。Sinc插值原理如下：在基帶信號下，卷積核是sinc函數(shù)sin(兀x)h(x)=sinc(x)=(12)兀x，、_-'插值信號為g(x)=工g(i)sinc(x-i)(13)di即為所有輸入樣本的加權平均。可通過頻域來理解，如圖6所示，采樣信號g(i)的頻譜G(f)等于以采樣率重復的信dd號頻譜。為了重建信號g(x)，只需要一個周期頻譜(如基帶周期)，因此需要理想矩形低通濾波器在頻域中提取基帶頻譜(如圖6)所示。已知該理想濾波器在時域中是sinc函數(shù)。由于頻域相乘相當于時域卷積，故插值可以通過與sinc核的卷積來實現(xiàn)。圖6理想低通濾波

12、器怎樣對采樣信號進行插值5點目標成像matlab仿真5.1距離多普勒算法距離多普勒算法(RDA)是在1976年至1978年為民用星載SAR提出的，它兼顧了成熟、簡單、高效和精確等因素，至今仍是使用最廣泛的成像算法。它通過距離和方位上的頻域操作，到達了高效的模塊化處理要求，同時又具有了一維操作的簡便性。圖7示意了RDA的處理流程。這里主要討論小傾斜角及短孔徑下的基本RDA處理框圖。1.當數(shù)據(jù)處在方位時域時，可通過快速卷積進行距離壓縮。也就是說，距離FFT后隨即進行距離向匹配濾波，再利用距離IFFT完成距離壓縮。回波信號為：14)s(t,s)二Awt-2R(s)/cw(ss)00racxexp-j

13、4kfR(s)/cexpj兀K(t-2R(s)/c)20r距離向壓縮后的信號為：s(t,s)=IFFTS(f,s)H(f)t0tt二Apt2R(s)/cw(ss)expj4kfR(s)/c0rac0rc15)ff2h(打=rectfKTTexp-加丘exp-j加如16)2通過方位FFT將數(shù)據(jù)變換至距離多普勒域，多普勒中心頻率估計以及大部分后續(xù)操作都在該域進行。方位向傅里葉變換后信號為：S(t,f)=FFTs(t,s)1ssrc2R(f)=Apt“W(ff)0rca4kfRf2Xexp-jiexpj兀rcKa3.在距離多普勒域進行隨距離時間及方位頻率變化的RCMC，該域中同一距離上的一組目標軌跡

14、相互重合。RCMC將距離徙動曲線拉直到與方位頻率軸平行的方向。這里可以采用最近鄰點插值法或者sine插值法，具體插值方法見前面。假設RCMC插值是精確的，信號變?yōu)椋?17)2RS(t,f)=Ap(t)W(ff)20rcassc4kfR】rxexp-j00expjkc(18)4通過每一距離門上的頻域匹配濾波實現(xiàn)方位壓縮。匚Ka為進行方位壓縮，將RCMC后的S(t,f)乘以頻域匹配濾波器H(f)。azsf2H(f)二expjfrazsKa21)(19)S(t,f)=S(t,f)H(f)3s2sazsssc=Ap(t2R/c)W(ff)expjf0Ro)0r0asscc5最后通過方位IFFT將數(shù)據(jù)變

15、換回時域，得到壓縮后的復圖像。復原后的圖像為:s(t,s)=IFFTS(t,f)acs3s二Ap(t-2R/c)p(s)0r0a4kfRxexp-jooexpj2兀fscsc圖8距離多普勒算法流程圖5.2 ChirpScaling算法距離多普勒算法具有諸多優(yōu)點，但是距離多普勒算法有兩點不足：首先，當用較長的核函數(shù)提高距離徙動校正（RCMC）精度時，運算量較大；其次，二次距離壓縮（SRC）對方位頻率的依賴性問題較難解決，從而限制了其對某些大斜視角和長孔徑SAR的處理精度。ChirpScaling算法避免了RCMC中的插值操作，通過對Chirp信號進行頻率調制，實現(xiàn)了對該信號的尺度變換或平移。圖8

16、顯示了ChirpScaling算法處理流程。這里主要討論小傾斜角及短孔徑下的基本CSA處理框圖。主要步驟包括四次FFT和三次相位相乘。1通過方位向FFT將數(shù)據(jù)變換到距離多普勒域。2. 通過相位相乘實現(xiàn)ChirpScaling操作，使所有目標的距離徙動軌跡一致化。這是第一步相位相乘。用以改變線調頻率尺度的ChirpScaling二次相位函數(shù)為：H（t,f;R）=exp加丫（f;R）a（f）（t-2Rf；R）2（22）1asaBac3. 通過距離向FFT將數(shù)據(jù)變到二維頻域。4. 通過與參考函數(shù)進行相位相乘，同時完成距離壓縮、SRC和一致RCMC。這是第二步相位相乘。用于距離壓縮，距離徙動校正的相位

17、函數(shù)寫為：H(f,f;R)=expjn2ras1丫(f；r)1+af)aBaf2rXexpj4兀Rsa(f)fcr(23)5. 通過距離向IFFT將數(shù)據(jù)變回到距離多普勒域。6通過與隨距離變化的匹配濾波器進行相位相乘，實現(xiàn)方位壓縮。此外，由于步驟2中的ChirpScaling操作，相位相乘中還需要附加一項相位校正。這是第三步相位相乘。補償由ChirpScaling引起的剩余相位函數(shù)是：H（t,f;R）二expj2Rf2-f2expj0（f;R）（24）2raBVBaMaAaB7.最后通過方位向IFFT將數(shù)據(jù)變回到二維時域，即SAR圖像域。圖8ChirpScaling算法流程圖簡而言之，R-D算法

18、是將徙動曲線逐一校正，CS算法是以某一徙動曲線為參考，在Doppler域內消除不同距離門的徙動曲線的差異，令這些曲線成為一組相互，平行”的曲線，然后在二維頻率域內統(tǒng)一的去掉距離徙動。通俗一點就是，RD算法是將彎曲的信號一根根掰直，而CS算法是先把所有信號都掰得一樣彎，然后再統(tǒng)一掰直。6仿真結果6.1使用最近鄰點插值的距離多普勒算法仿真結果本文首先對5個點目標的回波信號進行了仿真，5個點目標構成了矩形的4個頂點和中心，其坐標分別如下，格式為（方位向,距離向,后向反射系數(shù)）：09750110097501501000010102501100102501圖9的上圖是距離向壓縮后的圖像，從圖中可以看到5

19、條回波信號（其中有幾條部分重合，但仍能看出來）目標回波信號存在明顯的距離徙動，需要進行校正。圖9的下圖是通過最近鄰點插值法校正后的圖像，可以看出圖像基本被校正為直線。圖9距離向壓縮后最近鄰點插值的結果圖10為進行方位向壓縮后形成的圖像，可以明顯看出5個點目標，并且5個點目標構成了矩形的四個頂點及其中心。圖10通過最近鄰點插值生成的點目標圖像6.2 使用最近鄰點插值的距離多普勒算法仿真結果圖11上圖為通過距離壓縮后的圖像，圖11的下圖為通過sine插值法校正后的圖像。圖11距離向壓縮后sine插值的結果圖12為進行方位向壓縮后形成的圖像，可以明顯看出5個點目標，并且5個點目標構成了矩形的四個頂點

20、及其中心。圖12通過sine插值生成的點目標圖像6.3 ChirpScaling算法仿真結果同樣，在ChirpSealing中，對5個點目標的回波信號進行了仿真，5個點目標構成了矩形的4個頂點和中心，其坐標分別如下，格式為（方位向,距離向，后向反射系數(shù)）：1200011250-50112505011150-5011150501圖13是仿真的雷達回波信號圖圖13仿真出來的SAR回波信號圖14是經過第一次相位校正之后，通過距離向壓縮后的距離時域-方位時域信號圖（ChirpSealing算法的七個步驟中并不包含該信號，該信號是將步驟2之后的信號通過方位向傅里葉逆變換，再進行距離向壓縮得到的，只為了驗

21、證原理）。按照理論，該圖中所有點的距離徙動都應該一樣。從圖中大致可看出，五個點的距離徙動是差不多的。圖14ChirpSealing第2步之后、經過距離向壓縮得到的圖圖15為步驟5之后，信號距離壓縮，距離徙動校正之后的距離多普勒域中的信號圖。圖15距離徙動校正之后的圖圖16為步驟6之后，消除相位偏移的圖。圖16消除相位偏移的圖圖17為通過ChirpSealing算法生成的點目標圖像圖17通過ChirpSealing算法生成的點目標圖像6.4幾種算法比較本文討論了距離多普勒算法和ChirpSealing算法，其中距離多普勒算法考慮了最近鄰點插值和sine插值兩種插值方法。距離多普勒算法兼顧了成熟、

22、簡單、高效和精確等因素，至今仍被廣泛使用，但是距離多普勒算法有兩點不足：首先，當用較長的核函數(shù)提高距離徙動校正RCMC）精度時，運算量較大；其次，二次距離壓縮（SRC）對方位頻率的依賴性問題較難解決，從而限制了其對某些大斜視角和長孔徑SAR的處理精度。最近鄰點插值的優(yōu)點是速度快，該插值的運行時間為2.267137秒，缺點是不夠精確；sine插值的優(yōu)點是精確，該方法的運行時間為29.148728秒，缺點是速度慢；ChirpSealing算法避免了插值運算，提高了速度，運行時間為0.323327秒，但是其算法較為復雜。%=%文件名：NearSAR.m%作者：徐一凡%功能：合成孔徑雷達距離多普勒算法

23、點目標成像%=clear;clc;closeall;%=%常數(shù)定義C=3e8;%雷達參數(shù)Fc=1e9;lambda=C/Fc;%目標區(qū)域參數(shù)Xmin=0;Xmax=50;Yc=10000;Y0=500;%光速%載頻1GHz%波長%目標區(qū)域方位向范圍Xmin,Xmax%成像區(qū)域中線%目標區(qū)域距離向范圍Yc-YO,Yc+YO%成像寬度為2*Y0%軌道參數(shù)V=100;H=5000;R0=sqrt(Yc"2+2);%天線參數(shù)D=4;%SAR的運動速度100m/s%高度5000m%最短距離%方位向天線長度Lsar=lambda*RO/D;%SAR合成孔徑長度，合成孔徑雷達成像算法與實現(xiàn)P.1O

24、OTsar=Lsar/V;%慢時間域參數(shù)Ka=-2*V"2/lambda/R0;Ba二abs(Ka*Tsar);PRF二Ba;又為復頻率，所以等于PRT=1/PRF;ds=PRT;%SAR照射時間%多普勒頻域調頻率P.93%多普勒頻率調制帶寬%脈沖重復頻率，PRF其實為多普勒頻率的采樣率,%脈沖重復時間%慢時域的時間步長Nslow=ceil(Xmax-Xmin+Lsar)/V/ds);%慢時域的采樣數(shù)，ceil為取整函數(shù)，結合P.76的圖理解%nextpow2為最靠近2的幕次函數(shù)，這里為Nslow=2nextpow2(Nslow);fft變換做準備sn=linspace(Xmin-L

25、sar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢時間域的時間矩陣PRT=(Xmax-Xmin+Lsar)/V/Nslow;%由于Nslow改變了，所以相應的一些參數(shù)也需要更新，周期減小了PRF=1/PRT;ds=PRT;%快時間域參數(shù)設置%脈沖持續(xù)時間5usTr=5e-6;Br=30e6;%chirp頻率調制帶寬為30MHzKr二Br/Tr;%chirp調頻率Fsr=2*Br;%快時域采樣頻率，為3倍的帶寬dt=1/Fsr;%快時域采樣間隔Rmin=sqrt(Yc-Y0廠2+2);Rmax=sqrt(Yc+Y0廠2+2+(Lsar/2廠2);Nfast=ceil(2*(Rm

26、ax-Rmin)/C/dt+Tr/dt);%快時域的采樣數(shù)量Nfast=2"nextpow2(Nfast);%更新為2的幕次，方便進行fft變換tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);%快時域的離散時間矩陣dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;%更新間隔Fsr=1/dt;%分辨率參數(shù)設置%距離向分辨率%方位向分辨率''%點目標的數(shù)量%點目標位置，這里設置了5個點目標，構DY=C/2/Br;DX=D/2;%點目標參數(shù)設置Ntarget=5;%點目標格式x,y,反射系數(shù)sigmaPtarget=Xmin,Y

27、c-50*DY,l成一個矩形以及矩形的中心Xmin+50*DX,Yc-50*DY,1Xmin+25*DX,Yc,1Xmin,Yc+50*DY,1Xmin+50*DX,Yc+50*DY,1;disp('Parameters:')%參數(shù)顯示%=%生成回波信號K=Ntarget;N=Nslow;M=Nfast;T=Ptarget;Srnm=zeros(N,M);fork=1:1:Ksigma=T(k,3);Dslow=sn*V-T(k,1);disp('SamplingRateinfast-timedomain');disp(Fsr/Br)disp('Samp

28、lingNumberinfast-timedomain');disp(Nfast)disp('SamplingRateinslow-timedomain');disp(PRF/Ba)disp('SamplingNumberinslow-timedomain');disp(Nslow)disp('RangeResolution');disp(DY)disp('Cross-rangeResolution');disp(DX)disp('SARintegrationlength');disp(Lsar)disp(

29、'Positionoftargets');disp(Ptarget)%目標數(shù)目%慢時域的采樣數(shù)%快時域的采樣數(shù)%目標矩陣%生成零矩陣存儲回波信號%總共K個目標%得到目標的反射系數(shù)%方位向距離，投影到方位向的距離%實際距離矩陣R=sqrt(Dslow.2+T(k,2)2+H2);tau=2*R/C;%回波相對于發(fā)射波的延時Dfast=ones(N,l)*tm-tau'*ones(l,M);%(t-tau),其實就是時間矩陣，ones(N,l)和ones(1,M)都是為了將其擴展為矩陣phase二pi*Kr*Dfast."2-(4*pi/lambda)*(R

30、9;*ones(1,M);%相位，公式參見P.96Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*(abs(Dslow)<Lsar/2)'*ones(1,M);%由于是多個目標反射的回波，所以此處進行疊加end%=%距離-多普勒算法開始%距離向壓縮tic;tr=tm-2*Rmin/C;Refr=exp(j*pi*Kr*tr."2).*(0tr&trTr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr);Gr二abs(Sr);%開始進行距離彎曲補償正側

31、視沒有距離走動項主要是因為斜距的變化引起回波包絡的徙動%補償方法：最近鄰域插值法，具體為：先變換到距離多普勒域，分別對單個像素點計算出距離徙動量，得到距離徙動量與距離分辨率的比值，%該比值可能為小數(shù)，按照四舍五入的方法近似為整數(shù)，而后在該像素點上減去徙動量%方位向做fft處理再在頻域做距離彎曲補償Sa_RD=ftx(Sr);%方位向FFT變?yōu)榫嚯x多普域進行距離彎曲校正%距離徙動運算，由于是正側視，fdc=0,只需要進行距離彎曲補償。Kp=1;%計算或者預設預濾波比h=waitbar(0,'最近鄰域插值');%首先計算距離遷移量矩陣forn=1:N%總共有N個方位采樣form=1

32、:M%每個方位采樣上有M個距離采樣delta_R=(1/8)*(lambda/V廠2*(R0+(m-M/2)*C/2/Fsr)*(n-N/2)*PRF/N廠2;%距離遷移量P.160；(R0+(m-M/2)*C/2/Fsr):每個距離向點m的R0更新；(n-N/2)*PRF/N：不同方位向的多普勒頻率不一樣RMC=2*delta_R*Fsr/C;%此處為delta_R/DY，距離徒動了幾個距離單元delta_RMC=RMC-round(RMC);%距離徒動量的小數(shù)部分ifm+round(RMC)M%判斷是否超出邊界Sa_RD(n,m)=Sa_RD(n,M/2);elseifdelta_RMC=

33、0.5%五入Sa_RD(n,m)=Sa_RD(n,m+round(RMC)+1);else%四舍Sa_RD(n,m)=Sa_RD(n,m+round(RMC);endendendwaitbar(n/N)endclose(h)%=Sr_rmc=iftx(Sa_RD);%距離徙動校正后還原到時域Ga=abs(Sr_rmc);%方位向壓縮ta=sn-Xmin/V;Refa二exp(j*pi*Ka*ta."2).*(abs(ta)Tsar/2);Sa=iftx(ftx(Sr_rmc).*(conj(ftx(Refa).'*ones(1,M);Gar=abs(Sa);toc;%=%繪圖

34、colormap(gray);figure(l)subplot(211);row二tm*C/2-2008;col二sn*V-26;imagesc(row,col,255-Gr);%距離向壓縮，未校正距離徙動的圖像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel('距離向'),ylabel('方位向'),title('距離向壓縮，未校正距離徙動的圖像')，subplot(212);imagesc(row,col,255-Ga);%距離向壓縮，校正距離徙動后的圖像axis(Yc-Y0,Yc+Y0,Xmi

35、n-Lsar/2,Xmax+Lsar/2);xlabel('距離向'),ylabel('方位向'),title('距離向壓縮，校正距離徙動后的圖像')，figure(2)colormap(gray);imagesc(row,col,255-Gar);%方位向壓縮后的圖像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel('距離向'),ylabel('方位向'),title('方位向壓縮后的圖像'),%=%文件名:SincSAR.m%作者：徐一凡%功能:合

36、成孔徑雷達距離多普勒算法點目標成像%=clear;clc;closeall;%=%常數(shù)定義C=3e8;%光速%雷達參數(shù)Fc=1e9;%載頻1GHzlambda=C/Fc;%波長%目標區(qū)域參數(shù)Xmin=O;Xmax=5O;%目標區(qū)域方位向范圍Xmin,XmaxYc=1OOOO;YO=5OO;%成像區(qū)域中線%目標區(qū)域距離向范圍Yc-Y0,Yc+Y0%成像寬度為2*Y0%軌道參數(shù)V=1OO;H=5OOO;R0=sqrt(Yc"2+2);%天線參數(shù)D=4;%SAR的運動速度100m/s%高度5OOOm%最短距離%方位向天線長度Lsar=lambda*RO/D;%SAR合成孔徑長度，合成孔徑雷

37、達成像算法與實現(xiàn)P.1OOTsar=Lsar/V;%慢時間域參數(shù)Ka=-2*V"2/lambda/R0;Ba=abs(Ka*Tsar);PRF=Ba;又為復頻率，所以等于PRT=1/PRF;ds=PRT;%SAR照射時間%多普勒頻域調頻率P.93%多普勒頻率調制帶寬%脈沖重復頻率，PRF其實為多普勒頻率的米樣率，%脈沖重復時間%慢時域的時間步長Nslow二ceil(Xmax-Xmin+Lsar)/V/ds);%慢時域的采樣數(shù)，ceil為取整函數(shù)，結合P.76的圖理解Nslow=2nextpow2(Nslow);%nextpow2為最靠近2的幕次函數(shù)，這里為fft變換做準備sn=lin

38、space(Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢時間域的時間矩陣PRT=(Xmax-Xmin+Lsar)/V/Nslow;%由于Nslow改變了，所以相應的一些參數(shù)也需要更新，周期減小了PRF=1/PRT;ds=PRT;%快時間域參數(shù)設置Tr=5e-6;Br=30e6;%脈沖持續(xù)時間5us%chirp頻率調制帶寬為30MHzKr=Br/Tr;Fsr=2*Br;%chirp調頻率%快時域采樣頻率，為3倍的帶寬dt=1/Fsr;%快時域采樣間隔Rmin=sqrt(Yc-YO)2+H2);Rmax=sqrt(Yc+Y0廠2+2+(Lsar/2廠2);Nfa

39、st=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%快時域的采樣數(shù)量Nfast=2"nextpow2(Nfast);%更新為2的幕次，方便進行fft變換tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);%快時域的離散時間矩陣dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;%更新間隔Fsr=1/dt;%分辨率參數(shù)設置DY=C/2/Br;DX=D/2;%點目標參數(shù)設置Ntarget=5;%點目標格式x,y,反射系數(shù)sigmaPtarget=Xmin,Yc-50*DY,1成一個矩形以及矩形的中心%距離向分辨率%方位向分辨率%點

40、目標的數(shù)量%點目標位置，這里設置了5個點目標，構Xmin+50*DX,Yc-50*DY,1Xmin+25*DX,Yc,1Xmin,Yc+50*DY,1Xmin+50*DX,Yc+50*DY,l;disp('Parameters:')%參數(shù)顯示disp('SamplingRateinfast-timedomain');disp(Fsr/Br)disp('SamplingNumberinfast-timedomain');disp(Nfast)disp('SamplingRateinslow-timedomain');disp(PRF

41、/Ba)disp('SamplingNumberinslow-timedomain');disp(Nslow)disp('RangeResolution');disp(DY)disp('Cross-rangeResolution');disp(DX)disp('SARintegrationlength');disp(Lsar)disp('Positionoftargets');disp(Ptarget)%=%生成回波信號K二Ntarget;N=Nslow;M=Nfast;T=Ptarget;Srnm=zeros(N

42、,M);fork=1:1:Ksigma=T(k,3);Dslow=sn*V-T(k,1);R=sqrt(Dslow."2+T(k,2廠2+2);tau=2*R/C;Dfast=ones(N,1)*tm-tau'*ones(1,M);%目標數(shù)目%慢時域的采樣數(shù)%快時域的采樣數(shù)%目標矩陣%生成零矩陣存儲回波信號%總共K個目標%得到目標的反射系數(shù)%方位向距離，投影到方位向的距離%實際距離矩陣%回波相對于發(fā)射波的延時%(t-tau)，其實就是時間矩陣，ones(N,1)和ones(1,M)都是為了將其擴展為矩陣phase=pi*Kr*Dfast.2-(4*pi/lambda)*(R&

43、#39;*ones(1,M);%相位，公式參見P.96Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*(abs(Dslow)<Lsar/2)'*ones(1,M);%由于是多個目標反射的回波，所以此處進行疊加end%=%距離-多普勒算法開始%距離向壓縮tic;tr=tm-2*Rmin/C;Refr二exp(j*pi*Kr*tr."2).*(0tr&trTr);Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr);Gr=abs(Sr);%開始進行距離彎曲補償

44、正側視沒有距離走動項主要是因為斜距的變化引起回波包絡的徙動%補償方法：最近鄰域插值法，具體為：先變換到距離多普勒域，分別對單個像素點計算出距離徙動量，得到距離徙動量與距離分辨率的比值，%該比值可能為小數(shù)，按照四舍五入的方法近似為整數(shù)，而后在該像素點上減去徙動量%方位向做fft處理再在頻域做距離彎曲補償Sa_RD=ftx(Sr);%方位向FFT變?yōu)榫嚯x多普域進行距離彎曲校正%距離徙動運算，由于是正側視，fdc=O,只需要進行距離彎曲補償。Kp=1;%計算或者預設預濾波比%第二種方法進行截斷sinc插值進行距離徒動校正'h=waitbar(O,'Sinc插值')；P=4;%

45、4點sinc插值RMCmaxtix=zeros(N,M);forn=1:Nform=P:Mdelta_R=(l/8)*(lambda/V廠2*(R0+(m-M/2)*C/2/Fsr)*(n-N/2)*PRF/N廠2;%首先計算距離遷移量計算方法就是把斜距變換到距離多普勒域就知道了RMC=2*delta_R*Fsr/C;%距離徒動了幾個距離單元delta_RMC=RMC-round(RMC);%距離徒動量的小數(shù)部分fori=-P/2:P/2-1ifm+RMC+i>M%判斷是否超出邊界RMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,M)*sinc(pi*(-i+R

46、MC);elseRMCmaxtix(n,m)=RMCmaxtix(n,m)+Sa_RD(n,m+round(RMC)+i)*sinc(pi*(-i+delta_RMC);endendendwaitbar(n/N)endclose(h)%=Sr_rmc=iftx(RMCmaxtix);%距離徙動校正后還原到時域Ga=abs(Sr_rmc);%方位向壓縮ta=sn-Xmin/V;Refa=exp(j*pi*Ka*ta."2).*(abs(ta)<Tsar/2);Sa=iftx(ftx(Sr_rmc).*(conj(ftx(Refa).'*ones(1,M);Gar=abs(

47、Sa);toc;%=%繪圖colormap(gray);figure(1)subplot(211);row=tm*C/2-2008;col=sn*V-26;imagesc(row,col,255-Gr);%距離向壓縮，未校正距離徙動的圖像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabelC距離向，),ylabel(，方位向)，title('距離向壓縮，未校正距離徙動的圖像')，subplot(212);imagesc(row,col,255-Ga);%距離向壓縮，校正距離徙動后的圖像axis(Yc-Y0,Yc+Y0,Xmin-Lsa

48、r/2,Xmax+Lsar/2);xlabelC距離向'),ylabel('方位向'),title('距離向壓縮，校正距離徙動后的圖像')，figure(2)colormap(gray);imagesc(row,col,255-Gar);%方位向壓縮后的圖像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabelC距離向'),ylabel('方位向'),title('方位向壓縮后的圖像'),%ChirpScaling算法%徐一凡clearall;clc;%距離向參數(shù)range

49、:xdomainTr=200;%時寬200mBr=1;%帶寬1Kr=Br/Tr;%調頻斜率Fc=4;%載頻4Nfast=512;%為了快速運算Xc=1200;X0=150;%定義距離向范圍x=Xc+linspace(-X0,X0,Nfast);%x域序列:Xc-X0Xc+X0dx=2*X0/Nfast;%定義步長kx=linspace(-1/dx/2,1/dx/2,Nfast);%kx域序列%方位向參數(shù)cross-range:ydomainTa=300;%時寬300m,合成孔徑長度Ba=1;%帶寬1(1/m)Ka=Fc/Xc;%調頻斜率Ka=Ba/Ta=Fc/XcNslow=1024;%為了快

50、速運算Y0=200;y=linspace(-Y0,Y0,Nslow);%y域序列:-Y0Y0dy=2*Y0/Nslow;ky=linspace(-1/dy/2,1/dy/2,Nslow);%ky域序列%目標幾何關系targetgeometry%x坐標,y坐標，復后向散射系數(shù)Ptar=Xc,0,1+0jXc+50,-50,1+0jXc+50,50,1+0jXc-50,-50,1+0jXc-50,50,1+0j;disp('Positionoftargets');disp(Ptar)%生成SAR正交解調后的回波數(shù)據(jù)Srnm=zeros(Nfast,Nslow);N=size(Pta

51、r,l);%目標個數(shù)h=waitbar(0,'SAR回波生成');fori=1:1:Nxn=Ptar(i,1);yn=Ptar(i,2);sigma=Ptar(i,3);%提取每個目標的信息X=x.'*ones(1,Nslow);%擴充為矩陣Y=ones(Nfast,1)*y;%擴充為矩陣DX=sqrt(xn"2+(Y-yn)."2);%中間變量phase二pi*Kr*(X-DX)."2-2*pi*Fc*DX;%回波相位Srnm二Srnm+sigma*exp(j*phase).*(abs(X-DX)Tr/2).*(abs(Y-yn)Ta/2);%回波累加waitbar(i/N)endclose(h)tic;%數(shù)據(jù)準備phi0=-x'*sqrt(Fc“2-ky.“2);phi1二-Fc*x'*(1./sqrt(Fc“2-ky.“2);phi2=1/2*x'*(ky2./(F"2-ky2).J.5);Cs=ones(Nfast,1)*(Fc./sqrt(Fc"2-