合成孔径雷达(SAR)RD算法点目标成像与分析Matlab仿真
发布时间
阅读量：
阅读量
这是一个用于高分辨率合成孔径雷达（SAR）仿真系统的Matlab脚本文件。该文件的主要功能是通过模拟目标回波信号的生成、处理和成像过程，验证系统性能以及优化算法参数。以下是文件的主要功能模块和技术流程：
全局变量初始化
定义雷达工作频率（Fr）、中心频率（Fc）、 PRI码长度（Nrg）、码元周期（Tr）、采样率（Fs）、基线距离（Lg）等基本参数。
设置雷达的工作模式为“宽基线模式”并计算相关几何参数（斜距偏移角φrc, 坐标原点位置等）。目标配置
定义三个目标点A、B、C的位置坐标。
计算每个目标点到雷达平台的距离R0及其斜距偏移角φrc。
生成接收信号S_echo，并添加噪声干扰以模拟实际测量条件。主仿真流程
伪噪声信号生成
- 使用伪噪声码生成多径传播路径，并对各路径进行幅度衰减和相位补偿。
- 计算多径传播路径的角度差Δθr = (2h/c) cos(θrc) (cos(θrc) - cos(θb_t))。
距离压缩- 对接收到的信号进行频域去噪和平移校正。
- 应用二次距离补偿方法改进相位误差，并使用二维傅里叶变换进行距离压缩。
- 对压缩后的数据进行逆傅里叶变换以获得时间域的距离剖面。
方位向处理- 对压缩后的数据进行二维傅里叶变换以分离距离和方位信息。
- 进行角度校正和平移校正以提高方位分辨率。
- 利用二次相位补偿方法进一步消除余弦项引起的相位误差影响。
成像合成- 切割回波图并对其进行插值处理以提高采样密度。
- 绘制幅度包络切片图及其轮廓图以观察目标回波特性。
- 绘制分段信噪比（SNR）曲线评估回波质量。
结果验证- 输出关键性能指标如相对定位误差和信噪比SNR等结果。
- 比较不同阶段的数据对比结果以验证算法的有效性。
关键技术与算法
伪噪声码与多径传播
使用伪噪声码实现多径传播路径的选择与抵消，并通过几何模型计算各路径的角度差Δθ_r.
距离压缩技术
使用二维傅里叶变换分离距离和方位信息；通过二次距离补偿方法改进相位误差；应用非均匀插值提高信噪比SNR.
方位向处理
使用二维傅里叶变换分离空间信息；通过角度校正和平移校正提高方向分辨能力；利用二次相位补偿方法消除余弦项影响.
图像合成与可视化
切割成像窗口并插
文章目录
-
- 一、概述
- 二、仿真思路
-
- 1.概述
- 2.高分3号简介与基本参数
-
第三章 回波生成
- 第一部分 卫星运行速度计算
-
第二部分 几何学分析
-
第三部分 信号参数及时间轴构建
- （1）信号参数设置
- （2）时间轴构建
-
4.点目标回波生成
-
- (1)点目标坐标设置
-
- (2)回波生成
- 第一部分 卫星运行速度计算
-
第四章 低斜角度处理方案
4.1 距离收缩
4.2 基于方位的傅里叶变换
4.3 距离偏移校正
4.4 方位收缩
4.5 升采样操作
(1) 全局流程概述
(2) 升采样操作通过在频域中补充零点完成
**(3) 割裂视图生成'
第五部分 大斜视角处理
* 六、完整代码
* * 1.低斜视角处理
* 2.大斜视角处理一、概述
本文旨在利用Ian G. Cumming所著《合成孔径雷达成像算法及其在第六章中介绍的距离多普勒算法进行Matlab仿真研究，并结合高分三号 Synthetic Aperture Radar（SAR）卫星平台提供的信号参数展开分析。文章将从基本回波生成原理入手进行阐述，并通过代码实现与图形展示相结合的方式深入探讨相关技术要点。
需要注意的是：文中涉及的知识体系并非全面罗列所有相关内容，在完成前四章理论学习的基础上读者应具备基础了解。
此外：作为SAR成像领域的入门学习者作者在实际操作过程中难免会遇到诸多问题与不足 因此所得成果仅供参考
二、仿真思路
1.概述
即指在假设环境中对特定点的目标反射波进行数值模拟。
基本的仿真包括两个阶段：
回波生成 * 距离向包络Wr
* 方位向包络Wa
* 相位Phase
基于数据处理的小斜视角经过优化后能够适用于大斜视角问题（被视为一个升级版本）。然而，在这种情况下, 用于解决大斜角度问题的方法无法反向应用于解决小坡度角度问题
勘误：图中的256与320应分别为Naz与Nrg，下文再解释。
针对不同场景中的多维仿真建模 ，我们在模拟回波过程中对各个维度进行精确计算。具体而言，在这一环节中我们需分别针对不同类型的信号源计算其对应的参数值R_0、R_η以及波束中心穿越时刻T_η_c；然而，在后续的数据分析阶段则面临一定的限制——我们无法直接获取所有空间位置上的完整信息。这种限制也反映了该算法的主要缺陷——即仅能基于单个参考点的信息来进行推断工作。值得注意的是，在实际应用中选择合适的参考位置至关重要——近邻参考节点能够提供更为理想的成像质量。
2.高分3号简介与基本参数
GF-3号是中国高分辨率 Earth Observing System（HES）专项工程的重要组成部分之一。
该遥感平台具备1米分辨率雷达遥感能力。
它也是中国首个达到1米分辨率水平的C频段多极化合成孔径雷达（CSAR）系统。
该遥感成像系统由中国航天科技集团有限公司独立研制。
已知的基本参数如下：
%% GF3系统参数
% 几何参数
H = 755e3; % 轨道高度
% 天线参数
La = 15; % 天线尺寸
% 姿态参数（针对坐标原点的点A）
phi = 20*pi/180; % 俯仰角-20° - +20°
incidence= 20.5*pi/180; % 入射角
% 信号参数
c =3e8;%光速
f0 = 5.4e9;%雷达工作频率
Tr = 40e-6; % 距离向周期
Br = 3*5.995849160000000e+06; % 发射信号带宽这些参数与书本上的符号基本对应，接下来将基于上面的参数进行仿真
三、回波生成
1.卫星运行速度计算
通常采用直线几何模型在RD算法中。即采用等效雷达速度V_r（即V_g和V_s），并且满足关系式： V_r = V_g = V_s （适用于直线几何）。具体来说就是基于高中物理知识。利用卫星轨道高度H和地球半径R结合万有引力常量G进行计算。从而计算得出。
% 卫星轨道速度Vr计算
EarthMass = 6e24; %地球质量(kg)
EarthRadius = 6.37e6; %地球半径6371km
Gravitational = 6.67e-11; %万有引力常量
% 姿态参数
phi = 20 * pi / 180; % 俯仰角+20°
incidence = 20.5 * pi / 180; % 入射角
%计算等效雷达速度(卫星做圆周运动的线速度)
Vr = sqrt(Gravitational*EarthMass/(EarthRadius + H)); %第一宇宙速度2.几何
点目标仿真中涉及到的几何运算并不复杂，但比较重要，例如Rη瞬时斜距可以说是整个仿真当中最重要的参数，回波的计算会影响到成像的质量。
首先从下面的这一副星载的几何关系开始，雷达沿着方位向(Y轴)进行运动，在它的运动过程中，通过波束斜视角刚好能够照射到点目标的斜距即为景中心斜距；雷达与点目标之间的瞬时距离就是瞬时斜距，而最近斜距R0是固定的，通过俯仰角和高度就可以计算出来。
几何关系中的三个重要的角度，还有一个合成角(方位向包络)：
- 波束倾角Rc
- 入射倾角
- 仰视角度phi
- 合成倾角θ_bw（未在图中标注）
%景中心斜距R_eta_c和最近斜距R0，斜视角theta_rc由轨道高度、俯仰角、入射角计算得出
R_eta_c = H / cos(incidence); %景中心斜距
R0 = H / cos(phi);
theta_r_c = acos(R0/R_eta_c); %斜视角，单位为弧度;斜视角为4.6° % 合成孔径参数
rho_r = c / (2 * Fr); % 距离向分辨率
rho_a = Vr/Fa; % 距离向分辨率|La / 2
theta_bw = 0.886 * lambda / Lr; % 方位向3dB波束宽度
theta_syn = Vs / Vg * theta_bw; % 合成角宽度(斜面上的合成角)
Ls = R_eta_c * theta_syn; % 合成孔径长度3.信号参数与时间轴生成
(1)信号参数
在SAR仿真中采用的是二维（距离方向、方位方向）的线性调频信号；这些维度各自配置了一组信号参数；其中基础参数包括工作频率f₀、光速c等，并被共享使用。
%% 信号参数设置
% 电磁波参数
c = 3e+8; % 光速
Vs = Vr; % 卫星平台速度
Vg = Vr; % 波束扫描速度
La = 15; %方位向天线长度->椭圆的长轴
Lr=1.5;%距离向天线尺寸——>椭圆的短轴
f0 = 5.4e+9; % 雷达工作频率
lambda = c / f0; %电磁波波长
% 距离向信号参数
Tr = 40e-6; % 发射脉冲时宽
Br = 2.8 * 6e6; % 距离向信号带宽
Kr = Br / Tr; % 距离向调频率
alpha_os_r = 1.2; % 距离过采样率
Nrg = 2500; % 距离线采样点数
Fr = alpha_os_r * Br; % 距离向采样率
% 方位向信号参数
alpha_os_a = 1.23; % 方位过采样率
Naz = 1600; % 距离线数
delta_f_dop = 2 * 0.886 * Vr * (cos(theta_r_c)) / La; % 多普勒带宽
Fa = alpha_os_a * delta_f_dop; % 方位向采样率
Ta = 0.886 * lambda * R_eta_c / (La * Vg * cos(theta_r_c)); %目标照射时间
% 合成孔径参数
rho_r = c / (2 * Fr); % 距离向分辨率
rho_a = Vr/Fa; % 距离向分辨率|La / 2
theta_bw = 0.886 * lambda / Lr; % 方位向3dB波束宽度
theta_syn = Vs / Vg * theta_bw; % 合成角宽度(斜面上的合成角)
Ls = R_eta_c * theta_syn; % 合成孔径长度
fprintf("距离向分辨率:%.2f,方位向分辨率:%.2f\n\n",rho_r,rho_a)在本段中需注意的是Nrg与Naz这两个参数。它们决定了信号的距离向点数与方位向点数。例如当Naz=1600而Nrg=2500时将会生成一个164行×258列的二维阵列数据。值得注意的是 Naz与Nrg均与后续的时间轴相关联 从而引出了Taz与Trg这两个概念 这两个概念要特别注意不要与其他方向上的时间变量混淆 在进行包络计算时应分别使用Tr Ta作为参考基准 而在涉及时间轴相关的运算处理中则应采用Trg Taz作为运算基准
(2)时间轴生成
时间轴由时域轴和频域轴构成，在基本参数设置方面表现相似；然而需要注意的是每个轴的核心位置必须特别关注。
在距离向和方位向的时间轴上分别对应着时间参数2R_eta_c/c 和 time_eta_c 。通过分析几何关系图可知，在这两个参数中分别代表了雷达在距离方向上的电磁波运动中点以及方位方向上的雷达运动中点。
在距离向的频域轴上中点位置为0，在方位向的频域轴上中点位置即为多普勒中心频率(f_eta_c) 。在此基础上要求回波信号的频谱经过搬移后归零。
对于各个轴的数据我们都进行了扩展处理，并确保后续处理更加便捷。
%% 时间轴参数
Trg = Nrg / Fr;Taz = Naz / Fa;
%距离向/方位向采样时间间隔
Gap_t_tau = 1 / Fr;Gap_t_eta = 1 / Fa;
%距离向/方位向采样频率间隔
Gap_f_tau = Fr / Nrg;Gap_f_eta = Fa / Naz;
% 时间轴变量
time_tau_r = 2 * R_eta_c / c + (-Trg / 2:Gap_t_tau:Trg / 2 - Gap_t_tau); % 距离时间变量
time_eta_a = time_eta_c + (-Taz / 2:Gap_t_eta:Taz / 2 - Gap_t_eta); % 方位时间变量
% 随着距离向时间变化的最近斜距；c/2是因为距离向上一个时间包含两次电磁波来回
R0_tau_r = (time_tau_r * c / 2) * cos(theta_r_c);
Ext_R0_tau_r = repmat(R0_tau_r, Naz, 1); %扩展R0，用于计算变量Ka
% 频率变量
f_tau = (-Fr / 2:Gap_f_tau:Fr / 2 - Gap_f_tau); % 距离频率变量
f_tau=f_tau-(round(f_tau/Fr))*Fr;%混叠方程
f_eta = f_eta_c + (-Fa / 2:Gap_f_eta:Fa / 2 - Gap_f_eta); % 方位频率变量
f_eta=f_eta-(round((f_eta-f_eta_c)/Fa))*Fa;
% 时间轴
[Ext_time_tau_r, Ext_time_eta_a] = meshgrid(time_tau_r, time_eta_a); % 设置距离时域-方位时域二维网络坐标
% 频率轴
[Ext_f_tau, Ext_f_eta] = meshgrid(f_tau, f_eta); % 设置频率时域-方位频域二维网络坐标其中有一个变量需要特别说明,即Ext_R0_tau_r,它是通过扩展R0_tau_r得到的.即计算随着距离门(距离向点数)而变化的最近斜距,可理解为将距离转化为时间τ,每个τ值都对应一个R₀(同样也存在一个η_c值),这将用于后续计算方位向调频率K_a(K_a会根据R₀的变化而变化)
数学公式:
{R_0}' = \tau \times \frac{c}{2} \times \cos ({\theta _{r,c}})
4.点目标回波生成
(1)点目标坐标设置
基于前面获得的信号信息以及基础几何配置，在生成回波之前需要预先设定用于点目标仿真的各点坐标位置，并标明其以标准单位米为计量单位。
- 原点A的起始位置设定在坐标系的零点
- 近距参考站B的位置确定在x=500、y=500处
- 远程基准站C的定位参数中，其方位基准方向与近距站B一致；为了计算其至参考站的距离值，则需参考相关图表。（需要注意的是，在实际应用中可以选择不同的坐标设定）
%% 点目标(三个)坐标设置
% 设置目标点相对于景中心之间的距离
xA = 0;yA = 0; %A=(0,0)
xB = xA + 500;yB = xA + 500; %（500,500）
xC = H*tan(phi+theta_bw/2)-R0*sin(phi);%计算的C点距离向坐标
yC = xA + 500;
Position_x_r = [xA, xB, xC];
Position_y_a = [yA, yB, yC]; %点目标的坐标矩阵表示(2)回波生成
基于先前设定的基本信号参数与坐标系配置，在后续步骤中将进行回波生成过程的具体实施。具体而言，在回波生成过程中，在针对每个点目标计算距离向包络、方位向包络以及相位信息的基础上展开乘积运算以得到最终的结果值；由于各点目标具有不同的坐标位置这一特点，在特定时刻其瞬时斜距Rη、波束中心穿越时刻等关键参数也随之发生变化；经过对所有采样点的目标回波进行累加求和处理后即可模拟雷达对场景上各点目标进行采样操作的效果。
%% 生成回波S_echo
Target_num = 3; %目标数量
S_echo = zeros(Naz, Nrg);
for i = 1:Target_num
R0_Target = sqrt((R0 * sin(phi) + Position_x_r(i)) ^2+H^2); %对每个目标计算瞬时斜距
time_eta_c_Target = (Position_y_a(i)-R0_Target * tan(theta_r_c)) / Vr; %波束中心穿越时刻
% 计算目标点的瞬时斜距
R_eta = sqrt((R0_Target)^2+(Vr^2)*((Ext_time_eta_a - Position_y_a(i) / Vr).^2));
%R_eta=sqrt(H^2+(Position_x_r(i)-(-R0*sin(phi)))^2+(Position_y_a(i)-Ext_time_eta_a*Vr).^2);
% 距离向包络
Wr = (abs(Ext_time_tau_r-2*R_eta/c) <= Tr / 2);
% 方位向包络
% Wa = sinc((La*atan(Vg*(Ext_time_eta_a - time_eta_c_Target)./R0_Target)/lambda).^2);
Wa = ((La*atan(Vg*(Ext_time_eta_a - time_eta_c_Target)./(R0 * sin(phi) + Position_x_r(i)))/lambda).^2)<=Ta/2;
% Wa = abs(Ext_time_eta_a-time_eta_c_Target) <= Ta / 2;
% 相位
Phase = exp(-1j*4*pi*f0*R_eta/c) .* exp(+1j*pi*Kr*(Ext_time_tau_r - 2 * R_eta / c).^2);
% 接收信号叠加
S_echo_Target = Wr .* Wa .* Phase;
S_echo = S_echo + S_echo_Target;
%%%%%%%%%%%%%%格式化输出%%%%%%%%%%%%%%%%%%
R_eta_c_Target=R0_Target/cos(theta_r_c);
fprintf("当前目标:%d,坐标:(%.2f,%.2f),\n最近斜距R0=%.2f,景中心斜距R_eta_c=%.2f,波束中心穿越时刻=%.4f\n", i, Position_x_r(i), Position_y_a(i), R0_Target, R_eta_c_Target,time_eta_c_Target)
Time_tau_Point = round(((2*R0_Target)/(c*cos(theta_r_c))-time_tau_r(1))/Gap_t_tau);%目标的距离向标准坐标(距离门),Rηc的坐标
Time_tau_Point_RCMC = ceil(((2*R0_Target)/(c)-time_tau_r(1))/Gap_t_tau);%RCMC后点数坐标,R0的坐标
Time_eta_Point = Naz / 2 + (Position_y_a(i) / Vr) /Gap_t_eta;
fprintf("徙动校正前:脉冲点数坐标应为%d列(距)\n",Time_tau_Point);
fprintf("徙动校正后:点数坐标应为%d行(方),%d列\n\n",round(Time_eta_Point),Time_tau_Point_RCMC)
end
S_echo = S_echo .* exp(-2j*pi*f_eta_c*Ext_time_eta_a); %多普勒中心校正在回波生成的计算当中，有几个问题需要注意：
在计算过程中，在计算最近斜距R₀时，
采用以下公式：
{R_0}_{\eta} = \sqrt {{{({R_0} \cdot \sin (\phi) + x)}^2} + {H^2}}
其中，
R_0 \cdot \sin (\phi) 这一项表示沿方向的距离值，
即雷达飞行轨迹与地面法平面之间的距离，
该法平面是指在无仰角的情况下，
部分代码直接取值为 R_0 。
该法平面即为我们成像后得到的俯视图中最左边的位置。
基于几何关系进行推导分析，
给出了两个计算方法：
1）针对场景建模阶段进行几何推导，
求得当前点的目标瞬时斜距。
R(\eta ) = \sqrt {{R_0}{{_{t\arg et}}^2} + V{r^2}{{(\eta - \frac{y}{{Vr}})}^2}} 2) 建议读者认真考虑采用三维坐标系，并完成欧氏空间距离计算。
还有一个非常关键的概念是波束中心穿越时刻（beam center crossing time），即具有斜视角的一组雷达波束恰好准确地照射到了目标点的那个特定时刻。需要注意的是，在这种情况下计算得到的时间值是一个负数，在雷达飞行轨迹上表现为比飞至最近斜距的时间要早（即当雷达飞行至最近斜距时，则该时刻定义为零）；在时间轴构建过程中，默认是以波束中心穿越原点位置的时间作为基准值进行计算；但为了精确求解各维度上的包络值以及相应的相位信息，则必须基于具体的目标回波情况重新计算精确的时间参数。
为了完成这一过程，请先分别求取距离方向上的包络、方位方向上的包络以及相位相关的公式，并将它们依次相乘；接着对每个独立的目标点进行处理后累加所有的回波信号样本；最后通过以下运算式完成多普勒中心频移校正：% 多普勒中心频移校正：将频谱搬移到零频处。
S_echo = S_echo .* exp(-2j * pi * f_eta_c * Ext_time_eta_a);
回波的可视化以及设置的过采样率等因素均与之相关联，在此附有示例回波图以供参考（其变化主要受信号参数的影响较大）：
四、低斜视角处理
基于本文的仿真采用了高分三号（星载）雷达参数，并且设置有4.6度的倾斜角参数设置。因此，在理论上必须采用大倾斜角处理方法进行数据校正工作。本节将系统地分析低倾斜角处理方案的效果，并揭示其局限性，并在此基础上引出改进的大倾斜角处理方法作为后续研究的基础
1.距离压缩
匹配滤波也被称为压缩的原理，则是因为它的功能机制与降维打击这一概念具有相似性。具体而言，在实际应用中，则需要按照以下三个关键步骤来进行操作：第一步……第二步……第三步……
将信号从时域转换为距离频域 fft ，
通过与匹配滤波器相乘来处理信号，
将处理后的信号进行傅里叶反变换以恢复二维时域中的信息 ifft 。
Hf = (abs(Ext_f_tau) <= Br / 2) .* exp(+1j*pi*Ext_f_tau.^2/Kr);%滤波器
% 匹配滤波
S1_ftau_eta = fftshift(fft(fftshift(S_echo, 2), Nrg, 2), 2);
S1_ftau_eta = S1_ftau_eta .* Hf;
S1_tau_eta = fftshift(ifft(fftshift(S1_ftau_eta, 2), Nrg, 2), 2);所得到的效果，很像是将上面的每一个回波包络，看成一个面包，将它沿着两端向内进行压缩，最终得到一条线的效果。（后面的方位压缩，其实就是在方位向上进行所谓的"降维打击"）
2.方位向傅里叶变换
注：改写说明
%% 方位向傅里叶变换
S2_tau_feta = fftshift(fft(fftshift(S1_tau_eta, 1), Naz, 1), 1);在对信号进行前向和后向快速傅里叶变换之前及之后均执行快速傅里叶位移（FFTSHIFT）操作，在理论上等同于分别对其进行平移（对称）处理。由于该方法的理论基础较为基础，在现有文献中已有较多介绍不再赘述细节过程
需要注意的是，在完成距离压缩之后所得出的信号将在两个维度上均为时域序列；在此基础上执行方位向傅里叶变换；经过随后的距离迁移校正与方位压缩步骤；随后再次执行方位向傅里叶反变换；从而实现了对点目标的成像效果。
3.距离徙动校正
关于距离徙动校正的原理，可以参考下面的这篇博客
普遍采用SINCH内插算法来进行数据重构。该方法具有较高的计算效率和良好的重构性能。然而，在实际应用中存在一定的局限性：一方面其计算复杂度较高，在数据量较大的情况下可能导致系统响应时间变长；另一方面，在处理噪声干扰方面表现欠佳。因此，在这种情况下我们转而采用另一种常用的方法是采用相位预处理技术来进行数据重构。该方法的有效性略逊于SINCH内插算法（即sin插值），但在实现难度上具有较大的优势：操作过程相对更为简便，并且能够快速完成相关运算。
%% 距离徙动校正RCMC:采用相位补偿法
%虽然Ka是随着R0变化的，但是在相位补偿时需要假设R0是不变的，这也是相位补偿本身的局限性
delta_R = (((lambda * Ext_f_eta).^2) .* R0) ./ (8 * (Vr^2)); %距离徙动表达式
G_rcmc = exp((+4j * pi .* Ext_f_tau .* delta_R)./c); %补偿相位
S3_ftau_feta = fftshift(fft(fftshift(S2_tau_feta, 2), Nrg, 2), 2); %在方位向傅里叶变换的基础上进行距离向傅里叶变换
S3_ftau_feta = S3_ftau_feta .* G_rcmc; %与补偿相位相乘
S3_tau_feta_RCMC = fftshift(ifft(fftshift(S3_ftau_feta, 2), Nrg, 2), 2); %距离向傅里叶逆变换需要注意的关键点在于：上面所述之 R_₀ 指的是原始位置A₀处所对应的最近横向间距 R_₀ 这一特定数值；由于此时接收方向上的信号处于频域特性范围内，在进行方位压缩运算时其相关参数仅取决于固定不变的基础特性——即 R_ ₀ 本身；因此，在常规情况下，默认情况下以原始位置A₀处作为最佳成像效果的位置基准；然而，在实际应用中若采用相位补偿法进行处理，则必须重新定义新的参考基准——即更换当前处理对象的位置参数到其他目标位置（例如C）上进行相应处理
此外，在本研究中可以通过计算点目标的脉冲所在位置来校正移动带来的影响以验证校正效果。在本研究中，在距离方向上约有1570个脉冲计数值作为参考依据请根据以下公式：
R0 = {\rm{808390}}{\rm{.36}} < = > [\tau (1570)] \times (\frac{c}{2}) = {\rm{ 8}}{\rm{.0843e + 05}}
可以计算出，斜距上的误差约为35.66m
4.方位压缩
方位压缩与距离压缩具有相似性。值得注意的是，在此之前已经完成了方位向傅里叶变换。在此阶段，则可以直接应用匹配滤波器，并最终完成整个过程。
%% 方位压缩
% 根据变化的R0计算出相应的Ka矩阵(距离向变化，方位向不变)
Ka = 2 * Vr^2 * cos(theta_r_c)^2 ./ (lambda * Ext_R0_tau_r);
% 方位向匹配滤波器
Haz = exp(-1j*pi*Ext_f_eta.^2./Ka);
Offset = exp(-1j*2*pi*Ext_f_eta.*time_eta_c);%偏移滤波器，将原点搬移到Naz/2的位置，校准坐标
% 匹配滤波
S4_tau_feta = S3_tau_feta_RCMC .* Haz.*Offset;
S4_tau_eta = fftshift(ifft(fftshift(S4_tau_feta, 1), Naz, 1), 1);需要注意的是：
在方位滤波器处理之后的偏移量是将成像的目标点重新定位至Naz/2的坐标位置，并对应于方位方向的时间轴中心点。经过方位压缩处理后，在方位方向上将信号压缩为一个点。最终输出即为该目标点的空间定位信息。
可以看到基本符合我们前面的点目标坐标设置
5.升采样
(1)总体步骤
什么是上采样？实际上它是借助傅里叶变换的一个特性，在频域向信号末尾添加零点以实现时域中的更详细描绘。其核心概念等同于插值过程；需要注意的是关于补零的操作——教科书中指出：‘必须谨慎对待这一操作’。因为如果将零添加到特定位置会影响脉冲形态，则可能导致时域效果受损。总结一下——升采样的具体步骤如下：首先执行初始采样；其次在频谱中填充足够的零点；最后通过逆变换恢复高分辨率的时间序列。
- 切片：从Naz×Nrg维数中提取所需信号段落，并以32×32作为示例。
- 升采样（频域补零）：通过FFT和FFTshift将信号转换到二维频域空间中，在不影响频谱的前提下插入零数组。
- 角度校正：由于斜视角的影响，在点目标成像图中竖直方向呈现直线而水平方向倾斜的情况出现。我们需要调用imrotate函数对其进行旋转校正。
- 剖面：在获得包络后利用find和max函数确定峰值点作为中心点，在其周围分别沿着横纵坐标轴提取线条形成剖面；随后对横纵坐标分别计算绝对值并进行dB转换以获得最终结果。
%% 目标的升采样切片
%采用二维频域补零的方式进行升采样操作;升采样为了提高目标的切片细节丰富程度，便于观察
CutResolution = 32; %切片尺寸
Profile_Position = [901, 1570]; %切片的中心点位置
fprintf("\n点目标(C点)坐标对应的距离门到雷达距离为:%.2f\n",(time_tau_r(Profile_Position(2)))*(c/2))
%切片的升采样倍数：10*CutResolution；也就是在二维频域补多少个零
S5_tau_eta_Cut = S4_tau_eta(Profile_Position(1)-CutResolution/2:Profile_Position(1)+CutResolution/2, Profile_Position(2)-CutResolution/2:Profile_Position(2)+CutResolution/2);%切片
[S_zero_fill]=fft_zero_fill(S5_tau_eta_Cut,10*CutResolution);%补零函数
%由于斜视角的存在，在对升采样切片进行剖面分析前，将方位向和距离向都通过角度旋转校正到便于剖面的方向上
S5_tau_eta_Cut_UP_Azi=imrotate(S_zero_fill,rad2deg(theta_r_c),'bilinear', 'crop');%方位向切片无需角度校正
S5_tau_eta_Cut_UP_Ran = imrotate(S_zero_fill, rad2deg(theta_r_c), 'bilinear', 'crop'); %角度校正
%求解升采样切片包络的abs最大值坐标，用于下面的剖面
[UP_Profile_Position_Ran,UP_Profile_Position_Azi] = find(max(max(abs(S5_tau_eta_Cut_UP_Ran)))==abs(S5_tau_eta_Cut_UP_Ran));
%幅度db化+搬移峰值至0dB
%基于上面的升采样插值结果->获得剖面
Abs_S5_Azi = abs(S5_tau_eta_Cut_UP_Azi(:, UP_Profile_Position_Azi)); %方位向剖面
Abs_S5_Azi = Abs_S5_Azi / max(Abs_S5_Azi); %移动峰值点
Abs_S5_Ran = abs(S5_tau_eta_Cut_UP_Ran(UP_Profile_Position_Ran, :)); %距离向剖面
Abs_S5_Ran = Abs_S5_Ran / max(Abs_S5_Ran);
可以看到经过升采样，点目标的细节变得更加丰富，便于我们后续分析
(2)升采样(频域补零)
如图所示，在上述代码中调用了名为fft_zero_fill的函数，并且这是一个自编写的补零函数（完整的代码位于本节末尾）。具体思路如下：我们首先分析二维频谱（需要注意的是，在此过程中应避免调用 fftshift ）：
我们关注的频谱范围主要集中在黄色区域中，并且这些区域包含了重要信号的信息部分。为了避免全零数组对这些重要频段造成干扰，在执行信号处理时需要特别注意这一点。我打算按照图中所示的方向执行补零操作，并且通过这种方法能够显著提升数据处理后的结果质量。
观察到能量被很多0元素"压缩"到了四角的位置。这使得在进行二维时域反变换时变得更为简便。
function [S_FFT]=fft_zero_fill(x,nums)%进行补零；补零前一定要观察二维频谱，避免将补零的数组插到有能量的(黄色)区域，破坏原本的频谱！
S_FFT=fft(x,[],1);
S_FFT=fft(S_FFT,[],2);%二维频谱；这里不使用fftshift，避免频谱被搬移到中心
figure('name',"二维频域补零前");
imagesc(abs(S_FFT));%将二维频谱可视化，确定补零的范围(中点)
Y_Insert=zeros(nums,33);%在Y方向插入(方位向)
S_FFT=[S_FFT(1:21,:);Y_Insert;S_FFT(22:33,:)];%插入补零的数组;肉眼观察二维频谱的补零数组插入位置！！!
X_Insert=zeros(size(S_FFT,1),nums);%在X方向插入
S_FFT=[S_FFT(:,1:21),X_Insert,S_FFT(:,22:33)];%插入补零的数组;肉眼观察二维频谱的补零数组插入位置！！!
figure('name',"二维频域补零后");
imagesc(abs(S_FFT));%将二维频谱可视化，确定补零的范围(中点)
S_FFT=(ifft(ifft(S_FFT,[],1),[],2));%反变换，回到二维时域
%
% figure;
% imagesc(abs(ifft(ifft(S_FFT,[],1),[],2)))
end(3)剖面
经过了前面的升采样，我们也通过类似于这个代码（实际上很简单）：
Abs_S5_Azi = abs(S5_tau_eta_Cut_UP_Azi(:, UP_Profile_Position_Azi)); %方位向剖面完成了剖面的求解、幅值化、dB化的工作，得到了下面的剖面图：
在图中所示的红线位置对应着-13dB这一数值。由此可见，在剖面图中我们观察到峰值旁瓣值并未达至-13dB之下这一要求（这是必须满足的关键指标），这表明该算法在应用于星载SAR成像任务时存在一定的局限性，并且无法达到高分三号卫星所需的技术标准。因此建议采用大斜视角处理方案来提升图像质量。
五、大斜视角处理
大斜视角处理本质上只是对这几个关键环节进行了优化。其核心架构并未发生根本性变化。作为学习者，在初步了解相关知识时。
1.距离压缩->二次距离压缩(改进)
所谓二次距离压缩其实非常简单,其核心在于通过引入Km和Ksrc两个参数来解决调频偏移问题,基本上只需复制教材中的相关参数计算公式即可
%% 距离压缩(方式3)
D_feta_Vr=sqrt(1-((lambda*Ext_f_eta).^2)/(4*(Vr^2)));%徙动因子
K_src=(2*(Vr^2)*(f0^3)*(D_feta_Vr.^3))./(c*R0*(Ext_f_eta.^2));
Km=Kr./(1-Kr./K_src);%改进的调频率
Hf = (abs(Ext_f_tau) <= Br / 2) .* exp(+1j*pi*(Ext_f_tau.^2)./Km);%匹配滤波器
% 匹配滤波
S1_ftau_eta = fftshift(fft(fftshift(S_echo, 2), Nrg, 2), 2);
S1_ftau_eta = S1_ftau_eta .* Hf;%距离频域进行匹配滤波
S1_tau_eta = fftshift(ifft(fftshift(S1_ftau_eta, 2), Nrg, 2), 2);
2.方位向傅里叶变换
方位向傅里叶变换没有任何不同，这里不赘述了
3.距离徙动校正->引入新的徙动量(改进)
提出了一个新的距离迁移量，并对我们的相乘过程进行了优化。基本流程及其效果保持不变。
%% 距离徙动校正RCMC:采用相位补偿法
%虽然Ka是随着R0变化的，但是在相位补偿时需要假设R0是不变的
% delta_R = (((lambda * Ext_f_eta).^2) .* R0) ./ (8 * (Vr^2)); %距离徙动表达式
delta_R = R0*((1-D_feta_Vr)./D_feta_Vr); %距离徙动表达式
G_rcmc = exp((+4j * pi .* Ext_f_tau .* delta_R)./c); %补偿相位
S3_ftau_feta = fftshift(fft(fftshift(S2_tau_feta, 2), Nrg, 2), 2); %在方位向傅里叶变换的基础上进行距离向傅里叶变换
S3_ftau_feta = S3_ftau_feta .* G_rcmc; %与补偿相位相乘
S3_tau_feta_RCMC = fftshift(ifft(fftshift(S3_ftau_feta, 2), Nrg, 2), 2); %距离向傅里叶逆变换
%距离徙动校正结束
上图中红线代表点目标的多普勒中心频率，可以忽略
3.方位压缩->引入新的滤波器(改进)
在方位方向上采用了书上所述的新型匹配滤波器，并实现了更为显著的方位压缩效果。
%% 方位压缩
% 根据变化的R0计算出相应的Ka矩阵(距离向变化，方位向不变)
Ka = 2 * Vr^2 * cos(theta_r_c)^2 ./ (lambda * Ext_R0_tau_r);
% 方位向匹配滤波器
Haz = exp(-1j*pi*Ext_f_eta.^2./Ka);
Haz_BT=exp(+4j*pi*(Ext_R0_tau_r.*D_feta_Vr*f0)/c);%改进的方位滤波器
Offset = exp(-1j*2*pi*Ext_f_eta.*time_eta_c);%偏移滤波器，将原点搬移到Naz/2的位置，校准坐标
ABS_offset=exp(-2j*pi*Ext_f_eta*(29/Fa));%上面的Offset校准不够精确，观察方位向上原点与Naz/2的差值，二次校准
% 匹配滤波
S4_tau_feta = S3_tau_feta_RCMC .* Haz_BT.*Offset.*ABS_offset;
S4_tau_eta = fftshift(ifft(fftshift(S4_tau_feta, 1), Naz, 1), 1);ABS_offset代表校准过程中的一个固定偏差量。根据个人需求对基准点进行调整至Naz/2的位置其具体数值可暂不考虑。
得到了点目标成像结果，能够看出来变得更加规则一点了
4.升采样结果
使用与前面的小斜视角相同的升采样方法，得到的结果如下：
从图上可以看出：
该区域的轮廓图与等高线图呈现出高度美观度显著提升的状态，并且整体呈现出高度有序性特征。
该剖面轮廓整体衰减水平均在-13dB以下，并且已经满足成像要求。
六、完整代码
1.低斜视角处理
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 高分3号卫星参数RDA点目标仿真 %
% %
% 日期：2023年9月22日 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear;close all;
%代码格式化命令：MBeautify.formatCurrentEditorPage()
%% 卫星轨道参数
H = 755e3; %卫星轨道高度
% 卫星轨道速度Vr计算
EarthMass = 6e24; %地球质量(kg)
EarthRadius = 6.37e6; %地球半径6371km
Gravitational = 6.67e-11; %万有引力常量
% 姿态参数
phi = 20 * pi / 180; % 俯仰角+20°
incidence = 20.5 * pi / 180; % 入射角
%计算等效雷达速度(卫星做圆周运动的线速度)
Vr = sqrt(Gravitational*EarthMass/(EarthRadius + H)); %第一宇宙速度
%景中心斜距R_eta_c和最近斜距R0，斜视角theta_rc由轨道高度、俯仰角、入射角计算得出
R_eta_c = H / cos(incidence); %景中心斜距
R0 = H / cos(phi);
theta_r_c = acos(R0/R_eta_c); %斜视角，单位为弧度;斜视角为4.6°
%% 信号参数设置
% 电磁波参数
c = 3e+8; % 光速
Vs = Vr; % 卫星平台速度
Vg = Vr; % 波束扫描速度
La = 15; %方位向天线长度->椭圆的长轴
Lr=1.5;%距离向天线尺寸——>椭圆的短轴
f0 = 5.4e+9; % 雷达工作频率
lambda = c / f0; %电磁波波长
% 距离向信号参数
Tr = 40e-6; % 发射脉冲时宽
Br = 2.8 * 6e6; % 距离向信号带宽
Kr = Br / Tr; % 距离向调频率
alpha_os_r = 1.2; % 距离过采样率
Nrg = 2500; % 距离线采样点数
Fr = alpha_os_r * Br; % 距离向采样率
% 方位向信号参数
alpha_os_a = 1.7; % 方位过采样率(高过采样率避免鬼影目标)
Naz = 1600; % 距离线数
delta_f_dop = 2 * 0.886 * Vr * (cos(theta_r_c)) / La; % 多普勒带宽
Fa = alpha_os_a * delta_f_dop; % 方位向采样率
Ta = 0.886 * lambda * R_eta_c / (La * Vg * cos(theta_r_c)); %目标照射时间
% 景中心点(原点)的参数
time_eta_c = -R_eta_c * sin(theta_r_c) / Vr; % 波束中心穿越时刻
f_eta_c = 2 * Vr * sin(theta_r_c) / lambda; % 多普勒中心频率
% 合成孔径参数
rho_r = c / (2 * Fr); % 距离向分辨率
rho_a = Vr/Fa; % 距离向分辨率|La / 2
theta_bw = 0.886 * lambda / Lr; % 方位向3dB波束宽度
theta_syn = Vs / Vg * theta_bw; % 合成角宽度(斜面上的合成角)
Ls = R_eta_c * theta_syn; % 合成孔径长度
%% 时间轴参数
Trg = Nrg / Fr;Taz = Naz / Fa;
%距离向/方位向采样时间间隔
Gap_t_tau = 1 / Fr;Gap_t_eta = 1 / Fa;
%距离向/方位向采样频率间隔
Gap_f_tau = Fr / Nrg;Gap_f_eta = Fa / Naz;
% 时间轴变量
time_tau_r = 2 * R_eta_c / c + (-Trg / 2:Gap_t_tau:Trg / 2 - Gap_t_tau); % 距离时间变量
time_eta_a = time_eta_c + (-Taz / 2:Gap_t_eta:Taz / 2 - Gap_t_eta); % 方位时间变量
% 随着距离向时间变化的最近斜距；c/2是因为距离向上一个时间包含两次电磁波来回
R0_tau_r = (time_tau_r * c / 2) * cos(theta_r_c);
Ext_R0_tau_r = repmat(R0_tau_r, Naz, 1); %扩展R0，用于计算变量Ka
% 频率变量
f_tau = (-Fr / 2:Gap_f_tau:Fr / 2 - Gap_f_tau); % 距离频率变量
f_tau=f_tau-(round(f_tau/Fr))/Fr;%混叠方程
f_eta = f_eta_c + (-Fa / 2:Gap_f_eta:Fa / 2 - Gap_f_eta); % 方位频率变量
f_eta=f_eta-(round(f_eta/Fa))/Fa;
% 时间轴
[Ext_time_tau_r, Ext_time_eta_a] = meshgrid(time_tau_r, time_eta_a); % 设置距离时域-方位时域二维网络坐标
% 频率轴
[Ext_f_tau, Ext_f_eta] = meshgrid(f_tau, f_eta); % 设置频率时域-方位频域二维网络坐标
%% 点目标(三个)坐标设置
% 设置目标点相对于景中心之间的距离
xA = 0;yA = 0; %A=(0,0)
xB = xA + 500;yB = xA + 500; %（500,500）
xC = H*tan(phi+theta_bw/2)-R0*sin(phi);%计算的C点距离向坐标
yC = xA + 500;
Position_x_r = [xA, xB, xC];
Position_y_a = [yA, yB, yC]; %点目标的坐标矩阵表示
%% 生成回波S_echo
Target_num = 3; %目标数量
S_echo = zeros(Naz, Nrg);
for i = 1:Target_num
R0_Target = sqrt((R0 * sin(phi) + Position_x_r(i)) ^2+H^2); %对每个目标计算瞬时斜距
time_eta_c_Target = (Position_y_a(i)-R0_Target * tan(theta_r_c)) / Vr; %波束中心穿越时刻
% 计算目标点的瞬时斜距
R_eta = sqrt((R0_Target)^2+(Vr^2)*((Ext_time_eta_a - Position_y_a(i) / Vr).^2));
%R_eta=sqrt(H^2+(Position_x_r(i)-(-R0*sin(phi)))^2+(Position_y_a(i)-Ext_time_eta_a*Vr).^2);
% 距离向包络
Wr = (abs(Ext_time_tau_r-2*R_eta/c) <= Tr / 2);
% 方位向包络
%Wa = sinc((La*atan(Vg*(Ext_time_eta_a - time_eta_c_Target)./R0_Target)/lambda).^2);
Wa = (La*atan(Vg*(Ext_time_eta_a - time_eta_c_Target)./(R0 * sin(phi) + Position_x_r(i))/lambda).^2)<=Ta/2;
%Wa = abs(Ext_time_eta_a-time_eta_c_Target) <= Ta / 2;
% 相位
Phase = exp(-1j*4*pi*f0*R_eta/c) .* exp(+1j*pi*Kr*(Ext_time_tau_r - 2 * R_eta / c).^2);
% 接收信号叠加
S_echo_Target = Wr .* Wa .* Phase;
S_echo = S_echo + S_echo_Target;
%%%%%%%%%%%%%%格式化输出%%%%%%%%%%%%%%%%%%
R_eta_c_Target=R0_Target/cos(theta_r_c);
fprintf("当前目标:%d,坐标:(%.2f,%.2f),\n最近斜距R0=%.2f,景中心斜距R_eta_c=%.2f,波束中心穿越时刻=%.4f\n", i, Position_x_r(i), Position_y_a(i), R0_Target, R_eta_c_Target,time_eta_c_Target)
Time_tau_Point = round(((2*R0_Target)/(c*cos(theta_r_c))-time_tau_r(1))/Gap_t_tau);%目标的距离向标准坐标(距离门)
Time_eta_Point = Naz / 2 + (Position_y_a(i) / Vr) /Gap_t_eta;
fprintf("徙动校正前,点数坐标应为%d行(方),%d列(距)\n\n",round(Time_eta_Point),Time_tau_Point);
end
S_echo = S_echo .* exp(-2j*pi*f_eta_c*Ext_time_eta_a); %校正时间轴
%% 距离压缩
Hf = (abs(Ext_f_tau) <= Br / 2) .* exp(+1j*pi*Ext_f_tau.^2/Kr);%滤波器
% 匹配滤波
S1_ftau_eta = fftshift(fft(fftshift(S_echo, 2), Nrg, 2), 2);
S1_ftau_eta = S1_ftau_eta .* Hf;
S1_tau_eta = fftshift(ifft(fftshift(S1_ftau_eta, 2), Nrg, 2), 2);
%% 方位向傅里叶变换
S2_tau_feta = fftshift(fft(fftshift(S1_tau_eta, 1), Naz, 1), 1);
%% 距离徙动校正RCMC:采用相位补偿法
%虽然Ka是随着R0变化的，但是在相位补偿时需要假设R0是不变的
delta_R = (((lambda * Ext_f_eta).^2) .* R0) ./ (8 * (Vr^2)); %距离徙动表达式
G_rcmc = exp((+4j * pi .* Ext_f_tau .* delta_R)./c); %补偿相位
S3_ftau_feta = fftshift(fft(fftshift(S2_tau_feta, 2), Nrg, 2), 2); %在方位向傅里叶变换的基础上进行距离向傅里叶变换
S3_ftau_feta = S3_ftau_feta .* G_rcmc; %与补偿相位相乘
S3_tau_feta_RCMC = fftshift(ifft(fftshift(S3_ftau_feta, 2), Nrg, 2), 2); %距离向傅里叶逆变换
%距离徙动校正结束
%% 方位压缩
% 根据变化的R0计算出相应的Ka矩阵(距离向变化，方位向不变)
Ka = 2 * Vr^2 * cos(theta_r_c)^2 ./ (lambda * Ext_R0_tau_r);
% 方位向匹配滤波器
Haz = exp(-1j*pi*Ext_f_eta.^2./Ka);
Offset = exp(-1j*2*pi*Ext_f_eta.*time_eta_c);%偏移滤波器，将原点搬移到Naz/2的位置，校准坐标
% 匹配滤波
S4_tau_feta = S3_tau_feta_RCMC .* Haz.*Offset;
S4_tau_eta = fftshift(ifft(fftshift(S4_tau_feta, 1), Naz, 1), 1);
%% 目标的升采样切片
%采用二维频域补零的方式进行升采样操作;升采样为了提高目标的切片细节丰富程度，便于观察
CutResolution = 32; %切片尺寸
Profile_Position = [901, 1570]; %切片的中心点位置
fprintf("\n点目标(C点)坐标对应的距离门到雷达距离为:%.2f\n",(time_tau_r(Profile_Position(2)))*(c/2))
%切片的升采样倍数：10*CutResolution；也就是在二维频域补多少个零
S5_tau_eta_Cut = S4_tau_eta(Profile_Position(1)-CutResolution/2:Profile_Position(1)+CutResolution/2, Profile_Position(2)-CutResolution/2:Profile_Position(2)+CutResolution/2);%切片
[S_zero_fill]=fft_zero_fill(S5_tau_eta_Cut,10*CutResolution);%补零函数
%由于斜视角的存在，在对升采样切片进行剖面分析前，将方位向和距离向都通过角度旋转校正到便于剖面的方向上
S5_tau_eta_Cut_UP_Azi=imrotate(S_zero_fill,rad2deg(theta_r_c),'bilinear', 'crop');%方位向切片无需角度校正
S5_tau_eta_Cut_UP_Ran = imrotate(S_zero_fill, rad2deg(theta_r_c), 'bilinear', 'crop'); %角度校正
%求解升采样切片包络的abs最大值坐标，用于下面的剖面
[UP_Profile_Position_Ran,UP_Profile_Position_Azi] = find(max(max(abs(S5_tau_eta_Cut_UP_Ran)))==abs(S5_tau_eta_Cut_UP_Ran));
%幅度db化+搬移峰值至0dB
%基于上面的升采样插值结果->获得剖面
Abs_S5_Azi = abs(S5_tau_eta_Cut_UP_Azi(:, UP_Profile_Position_Azi)); %方位向剖面
Abs_S5_Azi = Abs_S5_Azi / max(Abs_S5_Azi); %移动峰值点
Abs_S5_Ran = abs(S5_tau_eta_Cut_UP_Ran(UP_Profile_Position_Ran, :)); %距离向剖面
Abs_S5_Ran = Abs_S5_Ran / max(Abs_S5_Ran);
%% 回波图可视化
figure('name', "回波可视化")
subplot(2, 2, 1);
imagesc(real(S_echo));
title('原始仿真信号实部');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 2);
imagesc(imag(S_echo));
title('原始仿真信号虚部');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 3);
imagesc(abs(S_echo));
title('原始仿真信号幅度');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 4);
imagesc(angle(S_echo));
title('原始仿真信号相位');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
%% 距离压缩可视化
%距离频谱可视化
figure('name', "回波频谱可视化")
subplot(2, 2, 1);
imagesc(real(S1_ftau_eta));
title('回波距离向实部');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 2);
imagesc(imag(S1_ftau_eta));
title('回波距离向虚部');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 3);
imagesc(abs(S1_ftau_eta));
title('回波距离向频谱幅度');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 4);
imagesc(angle(S1_ftau_eta));
title('回波距离向相位');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
%距离压缩结果
figure('name', "距离压缩时域结果")
subplot(1, 2, 1);
imagesc(real(S1_tau_eta));
title('实部');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(1, 2, 2);
imagesc(abs(S1_tau_eta));
title('幅度');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
%% 方位向傅里叶变换可视化
figure('name', "方位向傅里叶变换结果")
subplot(1, 2, 1);
imagesc(real(S2_tau_feta));
title('实部');
xlabel('距离向时间\tau');
ylabel('方位向频谱点数f_\eta');
subplot(1, 2, 2);
imagesc(abs(S2_tau_feta));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向频谱点数f_\eta');
%% 距离徙动校正可视化
figure('name', "距离徙动校正结果")
subplot(1, 2, 1);
imagesc(real(S3_tau_feta_RCMC));
title('实部');
xlabel('距离向时间\tau');
ylabel('方位向频率点数 f_\eta');
subplot(1, 2, 2);
imagesc(abs(S3_tau_feta_RCMC));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向频率点数 f_\eta');
%% 回波成像
figure('name', "点目标成像结果")
subplot(1, 2, 1);
imagesc(abs(S4_tau_eta));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向时间 \eta');
subplot(1, 2, 2);
imagesc(abs(S5_tau_eta_Cut)); %切片
title('点目标C的32×32切片(角度校正前)');
xlabel('距离向时间\tau');
ylabel('方位向时间\eta');
%% 升采样成像
figure('name', "升采样成像结果")
subplot(1, 2, 1);
imagesc(abs(S5_tau_eta_Cut_UP_Ran));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向时间 \eta');
subplot(1, 2, 2);
%contour((1:Nrg),(1:Naz),abs(S5_tau_eta_Cut_UP_Ran),20);
contour(abs(S5_tau_eta_Cut_UP_Ran), 30);
title('升采样轮廓图');
xlabel('距离向时间\tau');
ylabel('方位向时间\eta');
%% 切片包络幅度(dB)剖面
figure('name', "目标C的距离和方位幅度包络切片(插值后)")
subplot(1, 2, 1);
plot(mag2db(Abs_S5_Ran));
title('目标A距离向包络剖面(有角度校正)');
xlabel('距离向时间\tau');
ylabel('包络剖面幅度(dB)');
yticks(-50:2:50); % 设置刻度间隔为0.1
%画一条-13dB的参考线
x = linspace(0, 500, 100); % 生成 x 坐标的数据
y = repmat(-13,100); % 生成对应的 y 坐标的数据，全都是-13
hold on;
line(x, y,'Color','red'); % 绘制直线
subplot(1, 2, 2);
plot(mag2db(Abs_S5_Azi));
title('目标A方位向包络剖面(无角度校正)');
xlabel('方位向时间\eta');
ylabel('包络剖面幅度(dB)');
yticks(-50:2:50); % 设置刻度间隔为0.1
%画一条-13dB的参考线
x = linspace(0, 500, 100); % 生成 x 坐标的数据
y = repmat(-13,100); % 生成对应的 y 坐标的数据，全都是-13
hold on;
line(x, y,'Color','red'); % 绘制直线
%% 切片相位剖面
figure('name', "目标C的距离和方位幅度包络切片(插值后)")
subplot(1, 2, 1);
plot(angle(S5_tau_eta_Cut_UP_Ran(UP_Profile_Position_Ran, :)));
title('目标A距离向相位剖面(有角度校正)');
xlabel('距离向时间\tau');
ylabel('包络剖面幅度(dB)');
subplot(1, 2, 2);
plot(angle(S5_tau_eta_Cut_UP_Azi(:, UP_Profile_Position_Azi)));
title('目标A方位向相位剖面(无角度校正)');
xlabel('方位向时间\eta');
ylabel('包络剖面幅度(dB)');
%% 手动补零
function [S_FFT]=fft_zero_fill(x,nums)%进行补零；补零前一定要观察二维频谱，避免将补零的数组插到有能量的(黄色)区域，破坏原本的频谱！
S_FFT=fft(x,[],1);
S_FFT=fft(S_FFT,[],2);%二维频谱；这里不使用fftshift，避免频谱被搬移到中心
figure('name',"二维频域补零前");
imagesc(abs(S_FFT));%将二维频谱可视化，确定补零的范围(中点)
Y_Insert=zeros(nums,33);%在Y方向插入(方位向)
S_FFT=[S_FFT(1:21,:);Y_Insert;S_FFT(22:33,:)];%插入补零的数组;肉眼观察二维频谱的补零数组插入位置！！!
X_Insert=zeros(size(S_FFT,1),nums);%在X方向插入
S_FFT=[S_FFT(:,1:21),X_Insert,S_FFT(:,22:33)];%插入补零的数组;肉眼观察二维频谱的补零数组插入位置！！!
figure('name',"二维频域补零后");
imagesc(abs(S_FFT));%将二维频谱可视化，确定补零的范围(中点)
S_FFT=(ifft(ifft(S_FFT,[],1),[],2));%反变换，回到二维时域
%
% figure;
% imagesc(abs(ifft(ifft(S_FFT,[],1),[],2)))
end2.大斜视角处理
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 高分3号卫星参数RDA点目标仿真(大斜视角方法) %
% 作者：CYAN %
% 日期：2023年10月11日 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear;close all;
%代码格式化命令：MBeautify.formatCurrentEditorPage()
%基于Gaofen3_Rda.m：由于星载的斜视角超过3.5°，需要当作大斜视角进行处理，故此仿真脚本将进行大斜视角数据处理
% 使用参数Km用于匹配滤波解决调频率失配问题->二次距离压缩
% 大斜视角距离徙动校正
% 方位压缩的改进
%% 卫星轨道参数
H = 755e3; %卫星轨道高度
% 卫星轨道速度Vr计算
EarthMass = 6e24; %地球质量(kg)
EarthRadius = 6.37e6; %地球半径6371km
Gravitational = 6.67e-11; %万有引力常量
% 姿态参数
phi = 20 * pi / 180; % 俯仰角+20°
incidence = 20.5 * pi / 180; % 入射角
%计算等效雷达速度(卫星做圆周运动的线速度)
Vr = sqrt(Gravitational*EarthMass/(EarthRadius + H)); %第一宇宙速度
%景中心斜距R_eta_c和最近斜距R0，斜视角theta_rc由轨道高度、俯仰角、入射角计算得出
R_eta_c = H / cos(incidence); %景中心斜距
R0 = H / cos(phi);
theta_r_c = acos(R0/R_eta_c); %斜视角，单位为弧度;斜视角为4.6°
%% 信号参数设置
% 电磁波参数
c = 3e+8; % 光速
Vs = Vr; % 卫星平台速度
Vg = Vr; % 波束扫描速度
La = 15; %方位向天线长度->椭圆的长轴
Lr=1.5;%距离向天线尺寸——>椭圆的短轴
f0 = 5.4e+9; % 雷达工作频率
lambda = c / f0; %电磁波波长
% 距离向信号参数
Tr = 40e-6; % 发射脉冲时宽
Br = 2.8 * 6e6; % 距离向信号带宽
Kr = Br / Tr; % 距离向调频率
alpha_os_r = 1.2; % 距离过采样率
Nrg = 2500; % 距离线采样点数
Fr = alpha_os_r * Br; % 距离向采样率
% 方位向信号参数
alpha_os_a = 1.23; % 方位过采样率
Naz = 1600; % 距离线数
delta_f_dop = 2 * 0.886 * Vr * (cos(theta_r_c)) / La; % 多普勒带宽
Fa = alpha_os_a * delta_f_dop; % 方位向采样率
Ta = 0.886 * lambda * R_eta_c / (La * Vg * cos(theta_r_c)); %目标照射时间
% 景中心点(原点)的参数
time_eta_c = -R_eta_c * sin(theta_r_c) / Vr; % 波束中心穿越时刻
f_eta_c = 2 * Vr * sin(theta_r_c) / lambda; % 多普勒中心频率
% 合成孔径参数
rho_r = c / (2 * Fr); % 距离向分辨率
rho_a = Vr/Fa; % 距离向分辨率；书上另一定义为La / 2
theta_bw = 0.886 * lambda / Lr; % 方位向3dB波束宽度
theta_syn = Vs / Vg * theta_bw; % 合成角宽度(斜面上的合成角)
Ls = R_eta_c * theta_syn; % 合成孔径长度
fprintf("距离向分辨率:%.2f,方位向分辨率:%.2f\n\n",rho_r,rho_a)
%% 时间轴参数
Trg = Nrg / Fr;Taz = Naz / Fa;%采样的每个时间片为1/Fr或1/Fa;乘以点数计算出总的时长
%距离向/方位向采样时间间隔
Gap_t_tau = 1 / Fr;Gap_t_eta = 1 / Fa;
%距离向/方位向采样频率间隔
Gap_f_tau = Fr / Nrg;Gap_f_eta = Fa / Naz;
% 时间轴变量
time_tau_r = 2 * R0 / c + (-Trg / 2:Gap_t_tau:Trg / 2 - Gap_t_tau); % 距离时间变量
time_eta_a = time_eta_c + (-Taz / 2:Gap_t_eta:Taz / 2 - Gap_t_eta); % 方位时间变量
% 随着距离向时间变化的最近斜距；c/2是因为距离向上一个时间包含两次电磁波来回
R0_tau_r = (time_tau_r * c / 2) * cos(theta_r_c);
Ext_R0_tau_r = repmat(R0_tau_r, Naz, 1); %扩展R0，用于计算变量Ka
% 频率变量
f_tau = (-Fr / 2:Gap_f_tau:Fr / 2 - Gap_f_tau); % 距离频率变量
f_tau=f_tau-(round(f_tau/Fr))*Fr;%混叠方程
f_eta = f_eta_c + (-Fa / 2:Gap_f_eta:Fa / 2 - Gap_f_eta); % 方位频率变量
f_eta=f_eta-(round((f_eta-f_eta_c)/Fa))*Fa;
% 时间轴
[Ext_time_tau_r, Ext_time_eta_a] = meshgrid(time_tau_r, time_eta_a); % 设置距离时域-方位时域二维网络坐标
% 频率轴
[Ext_f_tau, Ext_f_eta] = meshgrid(f_tau, f_eta); % 设置频率时域-方位频域二维网络坐标
%% 点目标(三个)坐标设置
% 设置目标点相对于景中心之间的距离
xA = 0;yA = 0; %A=(0,0)
xB = xA + 500;yB = xA + 500; %（500,500）
xC = H*tan(phi+theta_bw/2)-R0*sin(phi);%计算的C点距离向坐标
yC = xA + 500;
Position_x_r = [xA, xB, xC];
Position_y_a = [yA, yB, yC]; %点目标的坐标矩阵表示
%% 生成回波S_echo
Target_num = 3; %目标数量
S_echo = zeros(Naz, Nrg);
for i = 1:Target_num
R0_Target = sqrt((R0 * sin(phi) + Position_x_r(i)) ^2+H^2); %对每个目标计算瞬时斜距
time_eta_c_Target = (Position_y_a(i)-R0_Target * tan(theta_r_c)) / Vr; %波束中心穿越时刻
% 计算目标点的瞬时斜距
R_eta = sqrt((R0_Target)^2+(Vr^2)*((Ext_time_eta_a - Position_y_a(i) / Vr).^2));
%R_eta=sqrt(H^2+(Position_x_r(i)-(-R0*sin(phi)))^2+(Position_y_a(i)-Ext_time_eta_a*Vr).^2);
% 距离向包络
Wr = (abs(Ext_time_tau_r-2*R_eta/c) <= Tr / 2);
% 方位向包络
% Wa = sinc((La*atan(Vg*(Ext_time_eta_a - time_eta_c_Target)./R0_Target)/lambda).^2);
Wa = ((La*atan(Vg*(Ext_time_eta_a - time_eta_c_Target)./(R0 * sin(phi) + Position_x_r(i)))/lambda).^2)<=Ta/2;
% Wa = abs(Ext_time_eta_a-time_eta_c_Target) <= Ta / 2;
% 相位
Phase = exp(-1j*4*pi*f0*R_eta/c) .* exp(+1j*pi*Kr*(Ext_time_tau_r - 2 * R_eta / c).^2);
% 接收信号叠加
S_echo_Target = Wr .* Wa .* Phase;
S_echo = S_echo + S_echo_Target;
%%%%%%%%%%%%%%格式化输出%%%%%%%%%%%%%%%%%%
R_eta_c_Target=R0_Target/cos(theta_r_c);
fprintf("当前目标:%d,坐标:(%.2f,%.2f),\n最近斜距R0=%.2f,景中心斜距R_eta_c=%.2f,波束中心穿越时刻=%.4f\n", i, Position_x_r(i), Position_y_a(i), R0_Target, R_eta_c_Target,time_eta_c_Target)
Time_tau_Point = round(((2*R0_Target)/(c*cos(theta_r_c))-time_tau_r(1))/Gap_t_tau);%目标的距离向标准坐标(距离门),Rηc的坐标
Time_tau_Point_RCMC = ceil(((2*R0_Target)/(c)-time_tau_r(1))/Gap_t_tau);%RCMC后点数坐标,R0的坐标
Time_eta_Point = Naz / 2 + (Position_y_a(i) / Vr) /Gap_t_eta;
fprintf("(仅参考)徙动校正前:脉冲点数坐标应为%d列(距)\n",Time_tau_Point);
fprintf("徙动校正后:点数坐标应为%d行(方),%d列\n\n",round(Time_eta_Point),Time_tau_Point_RCMC)
end
S_echo = S_echo .* exp(-2j*pi*f_eta_c*Ext_time_eta_a); %多普勒中心校正
%% 距离压缩(方式3)
D_feta_Vr=sqrt(1-((lambda*Ext_f_eta).^2)/(4*(Vr^2)));%徙动因子
K_src=(2*(Vr^2)*(f0^3)*(D_feta_Vr.^3))./(c*R0*(Ext_f_eta.^2));
Km=Kr./(1-Kr./K_src);%改进的调频率
Hf = (abs(Ext_f_tau) <= Br / 2) .* exp(+1j*pi*(Ext_f_tau.^2)./Km);%匹配滤波器
% 匹配滤波
S1_ftau_eta = fftshift(fft(fftshift(S_echo, 2), Nrg, 2), 2);
S1_ftau_eta = S1_ftau_eta .* Hf;%距离频域进行匹配滤波
S1_tau_eta = fftshift(ifft(fftshift(S1_ftau_eta, 2), Nrg, 2), 2);
%% 二次距离压缩相位误差计算（判断是否小于pi/2）
D_feta_Vr_O=sqrt(1-((lambda*f_eta_c).^2)/(4*(Vr^2)));
K_src_O=(2*(Vr^2)*(f0^3)*(D_feta_Vr_O.^3))./(c*R0*(Ext_f_eta.^2));
Km_O=Kr./(1-Kr./K_src_O);
delta_phi_srcf=pi*abs(Km_O-Km)*(Tr/2)^2<pi/2;%相位误差应小于pi/2;小于则数组全为1
%% 方位向傅里叶变换
S2_tau_feta = fftshift(fft(fftshift(S1_tau_eta, 1), Naz, 1), 1);
%% 距离徙动校正RCMC:采用相位补偿法
%虽然Ka是随着R0变化的，但是在相位补偿时需要假设R0是不变的
% delta_R = (((lambda * Ext_f_eta).^2) .* R0) ./ (8 * (Vr^2)); %距离徙动表达式
delta_R = R0*((1-D_feta_Vr)./D_feta_Vr); %距离徙动表达式
G_rcmc = exp((+4j * pi .* Ext_f_tau .* delta_R)./c); %补偿相位
S3_ftau_feta = fftshift(fft(fftshift(S2_tau_feta, 2), Nrg, 2), 2); %在方位向傅里叶变换的基础上进行距离向傅里叶变换
S3_ftau_feta = S3_ftau_feta .* G_rcmc; %与补偿相位相乘
S3_tau_feta_RCMC = fftshift(ifft(fftshift(S3_ftau_feta, 2), Nrg, 2), 2); %距离向傅里叶逆变换
%距离徙动校正结束
%% 方位压缩
% 根据变化的R0计算出相应的Ka矩阵(距离向变化，方位向不变)
Ka = 2 * Vr^2 * cos(theta_r_c)^2 ./ (lambda * Ext_R0_tau_r);
% 方位向匹配滤波器
Haz = exp(-1j*pi*Ext_f_eta.^2./Ka);
Haz_BT=exp(+4j*pi*(Ext_R0_tau_r.*D_feta_Vr*f0)/c);%改进的方位滤波器
Offset = exp(-1j*2*pi*Ext_f_eta.*time_eta_c);%偏移滤波器，将原点搬移到Naz/2的位置，校准坐标
ABS_offset=exp(-2j*pi*Ext_f_eta*(29/Fa));%上面的Offset校准不够精确，观察方位向上原点与Naz/2的差值，二次校准
% 匹配滤波
S4_tau_feta = S3_tau_feta_RCMC .* Haz_BT.*Offset.*ABS_offset;
S4_tau_eta = fftshift(ifft(fftshift(S4_tau_feta, 1), Naz, 1), 1);
%% 目标的升采样切片
%采用二维频域补零的方式进行升采样操作;升采样为了提高目标的切片细节丰富程度，便于观察
CutResolution = 32; %切片尺寸
Profile_Position = [872, 1914]; %切片的中心点位置
% fprintf("\n点目标(C点)坐标对应的距离门到雷达距离为:%.2f\n",(time_tau_r(Profile_Position(2)))*(c/2))
%切片的升采样倍数：10*CutResolution；也就是在二维频域补多少个零
S5_tau_eta_Cut = S4_tau_eta(Profile_Position(1)-CutResolution/2:Profile_Position(1)+CutResolution/2, Profile_Position(2)-CutResolution/2:Profile_Position(2)+CutResolution/2);%切片
[S_zero_fill]=fft_zero_fill(S5_tau_eta_Cut,10*CutResolution);%补零函数
%由于斜视角的存在，在对升采样切片进行剖面分析前，将方位向和距离向都通过角度旋转校正到便于剖面的方向上
S5_tau_eta_Cut_UP_Azi=imrotate(S_zero_fill,rad2deg(theta_r_c),'bilinear', 'crop');%方位向切片无需角度校正
S5_tau_eta_Cut_UP_Ran = imrotate(S_zero_fill, rad2deg(theta_r_c), 'bilinear', 'crop'); %角度校正
%求解升采样切片包络的abs最大值坐标，用于下面的剖面
[UP_Profile_Position_Ran,UP_Profile_Position_Azi] = find(max(max(abs(S5_tau_eta_Cut_UP_Ran)))==abs(S5_tau_eta_Cut_UP_Ran));
%幅度db化+搬移峰值至0dB
%基于上面的升采样插值结果->获得剖面
Abs_S5_Azi = abs(S5_tau_eta_Cut_UP_Azi(:, UP_Profile_Position_Azi)); %方位向剖面
Abs_S5_Azi = Abs_S5_Azi / max(Abs_S5_Azi); %移动峰值点
Abs_S5_Ran = abs(S5_tau_eta_Cut_UP_Ran(UP_Profile_Position_Ran, :)); %距离向剖面
Abs_S5_Ran = Abs_S5_Ran / max(Abs_S5_Ran);
%% 回波图可视化
figure('name', "回波可视化")
subplot(2, 2, 1);
imagesc(real(S_echo));
title('原始仿真信号实部');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 2);
imagesc(imag(S_echo));
title('原始仿真信号虚部');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 3);
imagesc(abs(S_echo));
title('原始仿真信号幅度');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 4);
imagesc(angle(S_echo));
title('原始仿真信号相位');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
%% 距离压缩可视化
%距离频谱可视化
figure('name', "回波频谱可视化")
subplot(2, 2, 1);
imagesc(real(S1_ftau_eta));
title('回波距离向实部');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 2);
imagesc(imag(S1_ftau_eta));
title('回波距离向虚部');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 3);
imagesc(abs(S1_ftau_eta));
title('回波距离向频谱幅度');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 4);
imagesc(angle(S1_ftau_eta));
title('回波距离向相位');
xlabel('距离向频谱点数f_\tau');
ylabel('方位向时间 \eta');
%距离压缩结果
figure('name', "距离压缩时域结果")
subplot(1, 2, 1);
imagesc(real(S1_tau_eta));
title('实部');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
subplot(1, 2, 2);
imagesc(abs(S1_tau_eta));
title('幅度');
xlabel('距离向时间 \tau');
ylabel('方位向时间 \eta');
%% 方位向傅里叶变换可视化
hold on;
figure('name', "方位向傅里叶变换结果")
subplot(1, 2, 1);
imagesc(real(S2_tau_feta));
title('实部');
xlabel('距离向时间\tau');
ylabel('方位向频谱点数f_\eta');
subplot(1, 2, 2);
imagesc(abs(S2_tau_feta));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向频谱点数f_\eta');
pix_fnc1 =find(f_eta==f_eta_c);
pin_fnc2 = 1053;
pix_fnc3 = 1053;
line([1, Nrg],[pix_fnc1, pix_fnc1], 'Color', 'red', 'LineWidth', 2);
line([1, Nrg],[pin_fnc2, pin_fnc2], 'Color', 'red', 'LineWidth', 2);
line([1, Nrg],[pix_fnc3, pix_fnc3], 'Color', 'red', 'LineWidth', 2);
hold off;
%% 距离徙动校正可视化
hold on;
figure('name', "距离徙动校正结果")
subplot(1, 2, 1);
imagesc(real(S3_tau_feta_RCMC));
title('实部');
xlabel('距离向时间\tau');
ylabel('方位向频率点数 f_\eta');
subplot(1, 2, 2);
imagesc(abs(S3_tau_feta_RCMC));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向频率点数 f_\eta');
pix_fnc1 =find(f_eta==f_eta_c);
pin_fnc2 = 1053;
pix_fnc3 = 1053;
line([1, Nrg],[pix_fnc1, pix_fnc1], 'Color', 'red', 'LineWidth', 2);
line([1, Nrg],[pin_fnc2, pin_fnc2], 'Color', 'red', 'LineWidth', 2);
line([1, Nrg],[pix_fnc3, pix_fnc3], 'Color', 'red', 'LineWidth', 2);
hold off;
%% 回波成像
figure('name', "点目标成像结果")
subplot(1, 2, 1);
imagesc(abs(S4_tau_eta));
title('幅度');
xlabel('距离向时间\tau');
ylabel('方位向时间 \eta');
subplot(1, 2, 2);
imagesc(abs(S5_tau_eta_Cut)); %切片
title('点目标C的32×32切片(角度校正前)');
xlabel('距离向时间\tau');
ylabel('方位向时间\eta');
%% 升采样成像
figure('name', "升采样成像结果")
subplot(2, 2, 1);
imagesc(abs(S5_tau_eta_Cut_UP_Ran));
title('幅度(距离向)');
xlabel('距离向时间\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 2);
contour(abs(S5_tau_eta_Cut_UP_Ran), 30);
title('升采样轮廓图(距离向)');
xlabel('距离向时间\tau');
ylabel('方位向时间\eta');
subplot(2, 2, 3);
imagesc(abs(S5_tau_eta_Cut_UP_Azi));
title('幅度(方位向)');
xlabel('距离向时间\tau');
ylabel('方位向时间 \eta');
subplot(2, 2, 4);
contour(abs(S5_tau_eta_Cut_UP_Azi), 30);
title('升采样轮廓图(方位向)');
xlabel('距离向时间\tau');
ylabel('方位向时间\eta');
%% 切片包络幅度(dB)剖面
figure('name', "目标C的距离和方位幅度包络切片(插值后)")
subplot(1, 2, 1);
plot(mag2db(Abs_S5_Ran));
title('目标A距离向包络剖面(有角度校正)');
xlabel('距离向时间\tau');
ylabel('包络剖面幅度(dB)');
yticks(-50:2:50); % 设置刻度间隔为0.1
%画一条-13dB的参考线
x = linspace(0, 500, 100); % 生成 x 坐标的数据
y = repmat(-13,100); % 生成对应的 y 坐标的数据，全都是-13
hold on;
line(x, y,'Color','red'); % 绘制直线
subplot(1, 2, 2);
plot(mag2db(Abs_S5_Azi));
title('目标A方位向包络剖面(无角度校正)');
xlabel('方位向时间\eta');
ylabel('包络剖面幅度(dB)');
yticks(-50:2:50); % 设置刻度间隔为0.1
%画一条-13dB的参考线
x = linspace(0, 500, 100); % 生成 x 坐标的数据
y = repmat(-13,100); % 生成对应的 y 坐标的数据，全都是-13
hold on;
line(x, y,'Color','red'); % 绘制直线
%% 切片相位剖面
figure('name', "目标C的距离和方位幅度包络切片(插值后)")
subplot(1, 2, 1);
plot(angle(S5_tau_eta_Cut_UP_Ran(UP_Profile_Position_Ran, :)));
title('目标A距离向相位剖面(有角度校正)');
xlabel('距离向时间\tau');
ylabel('包络剖面幅度(dB)');
subplot(1, 2, 2);
plot(angle(S5_tau_eta_Cut_UP_Azi(:, UP_Profile_Position_Azi)));
title('目标A方位向相位剖面(无角度校正)');
xlabel('方位向时间\eta');
ylabel('包络剖面幅度(dB)');
%% 距离与方位验算
%距离向
RelaDist_AandC_Ran=(1917-1251)*(1/Fr)*(c/2)*cos(theta_r_c)-(sqrt((R0*sin(phi)+14114.18)^2+H^2)-R0);
Ground_RelaDist_AandC_Ran=RelaDist_AandC_Ran*sin(phi+theta_bw/2);
fprintf("\n(距离向)A点与C点相对的斜距误差:%.2fm;映射到地面上为:%.2fm\n",RelaDist_AandC_Ran,Ground_RelaDist_AandC_Ran)
RelaDist_AandB_Ran=(1274-1251)*(1/Fr)*(c/2)*cos(theta_r_c)-(sqrt((R0*sin(phi)+500)^2+H^2)-R0);
fprintf("\n(距离向)A点与B点相对的斜距误差:%.2fm\n",RelaDist_AandB_Ran)
%方位向
RelaDist_AandC_Azi=(872-800)*(Vr/Fa)-500;
fprintf("\n(方位向)A点与C点(B点)相对的方位向误差:%.2fm\n",RelaDist_AandC_Azi)
%% 手动补零
function [S_FFT]=fft_zero_fill(x,nums)%进行补零；补零前一定要观察二维频谱，避免将补零的数组插到有能量的(黄色)区域，破坏原本的频谱！
S_FFT=fft(x,[],1);
S_FFT=fft(S_FFT,[],2);%二维频谱；这里不使用fftshift，避免频谱被搬移到中心
figure('name',"二维频域补零前");
imagesc(abs(S_FFT));%将二维频谱可视化，确定补零的范围(中点)
Y_Insert=zeros(nums,33);%在Y方向插入(方位向);肉眼观察二维频谱的补零数组插入位置！！!
S_FFT=[S_FFT(1:21,:);Y_Insert;S_FFT(22:33,:)];%插入补零的数组
X_Insert=zeros(size(S_FFT,1),nums);%在X方向插入;肉眼观察二维频谱的补零数组插入位置！！!
S_FFT=[S_FFT(:,1:22),X_Insert,S_FFT(:,23:33)];%插入补零的数组
figure('name',"二维频域补零后");
imagesc(abs(S_FFT));%将二维频谱可视化，确定补零的范围(中点)
S_FFT=(ifft(ifft(S_FFT,[],1),[],2));%反变换，回到二维时域
%
% figure;
% imagesc(abs(ifft(ifft(S_FFT,[],1),[],2)))
end
%%位置验证
function [temp]=PositionVal()
% 创建一个示例图像
image = [1 2 3; 4 5 6; 7 8 9];
% 绘制图像
imagesc(image);
colormap gray;
colorbar;
% 在图像上画一条二维直线
hold on;
line([1, 3], [1, 3], 'Color', 'red', 'LineWidth', 2);
hold off;
end全部评论 (0)
还没有任何评论哟~