ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

RDAC成像-Vancouver airport

2020-08-21 19:00:35  阅读:254  来源: 互联网

标签:采样 RDAC Rrd Nrg 距离 airport Vancouver fn 时域


1.程序

  1 %% RDA成像仿真
  2 clear;close all;clc;
  3 %% 得到可以进行后续信号处理的原始数据data(s_echo)
  4 % 载入参数
  5 load CD_run_params;
  6 load data_1;                % 分区1的数据
  7 s_echo = data_1;            % 原始数据记为s_echo,用于后续成像。
  8 %% 图1-原始数据
  9 figure;
 10 subplot(2,2,1);
 11 imagesc(real(s_echo));
 12 title('(a)实部');
 13 xlabel('距离时域(采样点)');
 14 ylabel('方位时域(采样点)');
 15 text(300,-60,'图1,原始数据');       % 给图1进行文字说明
 16 
 17 subplot(2,2,2);
 18 imagesc(imag(s_echo));
 19 title('(b)虚部');
 20 xlabel('距离时域(采样点)');
 21 ylabel('方位时域(采样点)');
 22 
 23 subplot(2,2,3);
 24 imagesc(abs(s_echo));
 25 title('(c)幅度');
 26 xlabel('距离时域(采样点)');
 27 ylabel('方位时域(采样点)');
 28 
 29 subplot(2,2,4);
 30 imagesc(angle(s_echo));
 31 title('(d)相位');
 32 xlabel('距离时域(采样点)');
 33 ylabel('方位时域(采样点)');
 34 %% 图2-原始数据频谱
 35 figure;
 36 subplot(2,2,1);
 37 imagesc(abs(fft(s_echo,[],1)));
 38 % imagesc(abs(fft(ifftshift(s_echo,1),[],1)));
 39 title('RD 频谱幅度');
 40 text(300,-60,'图2,原始数据频谱');
 41 subplot(2,2,2);
 42 imagesc(angle(fft(s_echo,[],1)));
 43 % imagesc(angle(fft(ifftshift(s_echo,1),[],1)));
 44 title('RD 频谱相位');
 45 subplot(2,2,3);
 46 imagesc(abs(fft2(s_echo)));
 47 % imagesc(abs(fft2(ifftshift(s_echo,1))));
 48 title('二维频谱幅度');
 49 subplot(2,2,4);
 50 imagesc(angle(fft2(s_echo)));
 51 % imagesc(angle(fft2(ifftshift(s_echo,1))));
 52 title('二维频谱相位');
 53 %% 定义一些参数
 54 Kr = -Kr;                       % 将调频率Kr改成负值
 55 BW_range = 30.111e+06;          % 脉冲宽度
 56 Vr = 7062;                      % 有效雷达速率
 57 Ka = 1733;                      % 方位调频率
 58 fnc = -6900;                    % 多普勒中心频率
 59 Fa = PRF;                       % 方位向采样率
 60 lamda = c/f0;                   % 波长
 61 T_start = 6.5959e-03;           % 数据窗开始时间
 62 
 63 Nr = round(Tr*Fr);              % 线性调频信号采样点数
 64 Nrg = Nrg_cells;                % 距离线采样点数
 65 Naz = Nrg_lines_blk;         % 每一个数据块的距离线数
 66 
 67 NFFT_r = Nrg;                   % 距离向FFT长度
 68 NFFT_a = Naz;                   % 方位向FFT长度
 69 %% 距离(方位)向时间,频率相关定义
 70 tr = T_start + ( -Nrg/2 : (Nrg/2-1) )/Fr;                % 距离时间轴
 71 fr = ( -NFFT_r/2 : NFFT_r/2-1 )*( Fr/NFFT_r );          % 距离频率轴
 72 ta = ( -Naz/2: Naz/2-1 )/Fa;                            % 方位时间轴
 73 fa = fnc + fftshift( -NFFT_a/2 : NFFT_a/2-1 )*( Fa/NFFT_a );    % 方位频率轴
 74 %% 距离压缩-距离向傅里叶变换
 75 S_range = fft(s_echo,NFFT_r,2);     % 进行距离向傅里叶变换,零频在两端。
 76 %% 图3-距离频域,方位时域,频谱(未距离压缩)
 77 figure;
 78 subplot(1,2,1);
 79 imagesc(real(S_range));
 80 title('(a)实部');
 81 xlabel('距离频域(采样点)');
 82 ylabel('方位时域(采样点)');
 83 text(280,-60,'图3,距离频域');
 84 text(340,-10,'未压缩');
 85 
 86 subplot(1,2,2);
 87 imagesc(abs(S_range));
 88 title('(b)幅度');
 89 xlabel('距离频域(采样点)');
 90 ylabel('方位时域(采样点)');
 91 %% 距离压缩-采用方式2-生成距离向匹配滤波器
 92 % 时域复制脉冲,末端补零,fft,再取复共轭。
 93 t_ref = ( -Nr/2 : (Nr/2-1) )/Fr;    % 用来生成距离MF的距离时间轴
 94 t_ref_mtx = ones(Naz,1)*t_ref;      % 矩阵形式
 95 w_ref = kaiser(Nr,2.5);             % 距离向,构建Kaiser窗,此为列向量。
 96 w_ref = ones(Naz,1)*(w_ref.');      % 构成矩阵形式,每一行都相同的加窗。
 97 
 98 s_ref = exp((1j*pi*Kr).*((t_ref_mtx).^2)); % 复制(发射)脉冲,未加窗 S3.26
 99 % s_ref = w_ref.*exp((1j*pi*Kr).*((t_ref_mtx).^2)); % 复制(发射)脉冲,加了窗。
100 
101 s_ref = [s_ref,zeros(Naz,Nrg-Nr)];      % 对复制脉冲,后端补零。
102 
103 S_ref = fft(s_ref,NFFT_r,2);            % 复制脉冲的距离傅里叶变换,零频在两端。
104 H_range = conj(S_ref);                  % 距离向匹配滤波器,零频在两端。
105 S_range_c = S_range.*H_range;           % 匹配滤波器,零频在两端。
106 s_rc = ifft(S_range_c,[],2);            % 完成距离压缩,回到二维时域。
107 % s_rc的长度为:Naz*Nrg。未去除弃置区。
108 %% 图4-距离频域,方位时域,频谱(已距离压缩)
109 figure;
110 subplot(1,2,1);
111 imagesc(real(S_range_c));
112 title('(a)实部');
113 xlabel('距离频域(采样点)');
114 ylabel('方位时域(采样点)');
115 text(280,-60,'图4,距离频域');       % 给图3进行文字说明
116 text(340,-10,'已压缩');
117 
118 subplot(1,2,2);
119 imagesc(abs(S_range_c));
120 title('(b)幅度');
121 xlabel('距离频域(采样点)');
122 ylabel('方位时域(采样点)');
123 %% 图5-二维时域(完成距离压缩)
124 figure;
125 subplot(1,2,1);
126 imagesc(real(s_rc));  % 这及其以下,都直接使用去除弃置区后的结果
127 title('(a)实部');
128 xlabel('距离时域(采样点)');
129 ylabel('方位时域(采样点)');
130 text(150,-60,'图5,二维时域');       % 给图4进行文字说明
131 text(172,-10,'完成压缩');
132 
133 subplot(1,2,2);
134 imagesc(abs(s_rc));
135 title('(b)幅度');
136 xlabel('距离时域(采样点)');
137 ylabel('方位时域(采样点)');
138 %% 变换到二维频域,进行SRC
139 s_rc = s_rc.*exp(-1j*2*pi*fnc.*(ta.'*ones(1,Nrg)));    % 数据搬移 零斜视fnc=0
140 S_2df = fft(s_rc,NFFT_a,1);             % 方位向傅里叶变换,到距离多普勒域
141 S_2df = fft(S_2df,Nrg,2);               % 距离向傅里叶变换,到二维频域
142 % !!!注意:距离向零频在两端。
143 % ====================================================================
144 % 设置方位频率轴——这是关键点
145 fa = fnc + fftshift(-NFFT_a/2:NFFT_a/2-1)/NFFT_a*Fa;     % 方位频率轴如此设置。
146 % =====================================================================
147 D_fn_Vr = sqrt(1-lamda^2.*(fa.').^2./(4*Vr^2));         % 大斜视角下的徙动因子
148 K_src = 2*Vr^2*f0^3.*D_fn_Vr.^3./(c*R0*(fa.').^2);      % 列向量
149 K_src_1 = 1./K_src;             % 列向量。为了后面能使用矩阵乘法,这里先求倒数
150 H_src = exp(-1j*pi.*K_src_1*(fr.^2)); % 二次距离压缩滤波器。距离向,零频在中间。
151 % 这是矩阵,大小Naz*Nrg
152 H_src = fftshift(H_src,2);      % (左右半边互换)距离向,零频在两端。 !!!这很关键!!!
153 
154 S_2df_src = S_2df.*H_src;       % 这一步点乘时,要注意两者的距离向频率轴应该对应上,不然会出错!!
155 % 这就是为什么上面的 H_src 要 fftshift 的原因!!
156 
157 S_rd = ifft(S_2df_src,[],2);       % 完成二次距离压缩(SRC),回到距离多普勒域。
158 
159 disp('在二维频域,采取方式2完成SRC,并变换到距离多普勒域');
160 %% 作图
161 figure;
162 imagesc(abs(S_rd));
163 title('完成距离压缩,完成SRC后,距离多普勒域(未RCMC)');
164 %% 变换到距离多普勒域,RCMC
165 % S_rd = fft(s_rc,NFFT_a,1);
166 % fa = fnc + fftshift(-NFFT_a/2:NFFT_a/2-1)/NFFT_a*Fa;    % 方位频率轴 [1 NFFT_a] 0频增大至最大正频,突变至最大负频增大至0
167 tr_RCMC = T_start + ( -Nrg/2 : (Nrg/2-1) )/Fr;   % 在新的距离线长度下的时间轴 [1 Nrg]
168 
169 R0_RCMC = (c/2).*tr_RCMC;       % 随距离线变化的R0,用来计算RCM和Ka [1 Nrg]
170 delta_Rrd_fn = lamda^2.*((fa.').^2)*(R0_RCMC)/(8*Vr^2);   %需要校正的RCM,各行缓慢增大,各列由0先增大再减小至0,式6.11 [NFFT_a 1]*[1 Nrg]=[NFFT_a Nrg]
171 % delta_Rrd_fn = ((1-D_fn_Vr)./D_fn_Vr)*R0_RCMC;      % 大斜视角下的RCM
172 
173 num_range = c/(2*Fr);   % 一个距离采样单元对应的长度
174 delta_Rrd_fn_num = delta_Rrd_fn./num_range; % 需要校正的RCM:各个方位频率下的RCM对应的距离采样单元数 [NFFT_a Nrg]
175 
176 R = 8;  % sinc插值核长度
177 S_rd_rcmc = zeros(NFFT_a,Nrg); % 用来存放RCMC后的值
178 for p = 1 : NFFT_a
179     for q = 1 : Nrg   % 此时距离向的长度是 (Nrg-Nr+1)=Nrg
180         delta_Rrd_fn_p = delta_Rrd_fn_num(p,q);
181         Rrd_fn_p = q + delta_Rrd_fn_p; %各方位频率、时域点的采样对应的距离单元数
182         Rrd_fn_p_zheng = ceil(Rrd_fn_p);        % ceil,向上取整
183         ii = ( Rrd_fn_p-(Rrd_fn_p_zheng-R/2):-1:Rrd_fn_p-(Rrd_fn_p_zheng+R/2-1)  );
184         rcmc_sinc = sinc(ii);   %式2.59
185         rcmc_sinc = rcmc_sinc/sum(rcmc_sinc);   % 插值核的归一化
186         % ii 是sinc插值过程的变量;
187         % g(x)=sum(h(ii)*g_d(x-ii)) = sum(h(ii)*g_d(ll));
188         
189         % 由于S_rd只有整数点取值,且范围有限。因此插值中要考虑它的取值溢出边界问题。
190         % 这里我采取循环移位的思想,用来解决取值溢出问题。
191         if (Rrd_fn_p_zheng-R/2) >Nrg     % 全右溢
192             ll = (Rrd_fn_p_zheng-R/2-Nrg:1:Rrd_fn_p_zheng+R/2-1-Nrg);
193         else
194             if (Rrd_fn_p_zheng+R/2-1) > Nrg    % 部分右溢
195                 ll_1 = (Rrd_fn_p_zheng-R/2:1:Nrg);
196                 ll_2 = (1:1:Rrd_fn_p_zheng+R/2-1-Nrg);
197                 ll = [ll_1,ll_2];
198             else
199                 if (Rrd_fn_p_zheng+R/2-1) < 1    % 全左溢(不可能发生,但还是要考虑)
200                     ll = (Rrd_fn_p_zheng-R/2+Nrg:1:Rrd_fn_p_zheng+R/2-1+Nrg);
201                 else
202                     if (Rrd_fn_p_zheng-R/2) < 1       % 部分左溢
203                         ll_1 = (Rrd_fn_p_zheng-R/2+Nrg:1:Nrg);
204                         ll_2 = (1:1:Rrd_fn_p_zheng+R/2-1);
205                         ll = [ll_1,ll_2];
206                     else
207                         ll = (Rrd_fn_p_zheng-R/2:1:Rrd_fn_p_zheng+R/2-1);
208                     end
209                 end
210             end
211         end
212         rcmc_S_rd = S_rd(p,ll);
213         S_rd_rcmc(p,q) = sum( rcmc_sinc.*rcmc_S_rd );
214     end
215 end
216 % S_rd_rcmc 就是RCMC后的距离多普勒域频谱。
217 %% 图6——距离多普勒域(未RCMC)
218 figure;
219 subplot(1,2,1);
220 imagesc(real(S_rd));
221 title('(a)实部');
222 xlabel('距离时域(采样点)');
223 ylabel('方位频域(采样点)');
224 text(150,-60,'图6,距离多普勒域');
225 text(172,-10,'未RCMC');
226 subplot(1,2,2);
227 imagesc(abs(S_rd));
228 title('(b)幅度');
229 xlabel('距离时域(采样点)');
230 ylabel('方位频域(采样点)');
231 %% 图7——距离多普勒域,RCMC后的结果
232 figure;
233 subplot(1,2,1);
234 imagesc(real(S_rd_rcmc));
235 title('(a)实部');
236 xlabel('距离时域(采样点)');
237 ylabel('方位频域(采样点)');
238 text(150,-60,'图7,距离多普勒域');
239 text(172,-10,'已RCMC');
240 
241 subplot(1,2,2);
242 imagesc(abs(S_rd_rcmc));
243 title('(b)幅度');
244 xlabel('距离时域(采样点)');
245 ylabel('方位频域(采样点)');
246 %% 方位压缩
247 % D_fn_Vr = sqrt(1-lamda^2.*(fa.').^2./(4*Vr^2));
248 Haz = exp(1j*4*pi.*(D_fn_Vr*R0_RCMC).*f0./c);   % 改进的方位向MF
249 S_rd_c = S_rd_rcmc.*Haz;            % 乘以匹配滤波器
250 s_ac = ifft(S_rd_c,[],1);           % 完成方位压缩,变到图像域。结束。
251 %% 图8——成像结果
252 sout = (s_ac);
253 sout = abs(sout)/max(max(abs(sout)));
254 G = 20*log10(sout+eps);                         % db显示
255 clim = [-55,0];                                 % 动态显示范围
256 figure;
257 imagesc(((0:Nrg-1)+first_rg_cell)/Fr*c/2+R0,((0:Naz-1)+first_rg_line)/Fa*Vr,G,clim);
258 axis xy;
259 title('RADARSAT-1数据,使用RD算法,成像结果')
260 xlabel('Range(m)')
261 ylabel('Azimuth(m)')
262 % colormap(gray);

2..文件

 

3.结果

 

 

 

 

 

 

标签:采样,RDAC,Rrd,Nrg,距离,airport,Vancouver,fn,时域
来源: https://www.cnblogs.com/heheda-jl/p/13542874.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有