SAR学习笔记-代码部分

$\left\{\begin{matrix} H_{0}:x\left ( t \right )=n\left ( t \right )\; \; \; 0\leq t\leq T \\ H_{1}:x\left ( t \right )=s\left ( t \right )+n\left ( t \right )\; \; \; 0\leq t\leq T \end{matrix}\right.$

为了进一步简化检测模型，先暂时不考虑接收到的回波信号相对雷达发射信号的时延以及能量衰减，即认为接收的回波信号即为雷达的发射信号 $s\left ( t \right )$ 。根据接收机的热噪声特性， $n\left ( t \right )$ 为高斯白噪声。

从信号处理的角度看，此时目标检测就是：根据雷达接收机接收到的信号，判断接收信号中是否含有目标信号。

当不存在干扰 $n\left ( t \right )$ 的情况下，雷达接收的信号 $x\left ( t \right )$ 只有两种情况下：有目标情况的 $s\left ( t \right )$ ；没有目标情况的0，此时在目标有无确定的情况下，雷达接收的信号是确定的，可以根据雷达接收的信号很容易判断目标的有无。

当干扰 $n\left ( t \right )$ 存在的情况下，无论目标有无，雷达接收的信号都是不确定的，是随机信号，直接根据雷达接收的信号判断雷达探测区域内的目标有无是比较困难的。针对这个问题，可以从概率论的角度来分析。

对接收的信号 $x\left ( t \right )$ 进行采样得到 $x_{k}=x\left ( t_{k} \right )$ ,考虑接收机热噪声 $n_{k}=n\left ( t_{k} \right )$ 是高斯白噪声，则随机变量 $x_{k}$ 是个高斯变量，且目标的有无不影响随机变量 $x_{k}$ 的方差（方差为噪声方差），只影响其均值（无目标，均值为0；有目标，均值为发射信号在该时刻的采样值）。此时 $x_{k}$ 在目标有无下的概率密度函数分别为：

$f\left ( x_{k}|H_{0} \right )=\frac{1}{\sqrt{2\pi \sigma }}e^{-\frac{x^{2}_{k}}{2\sigma ^{2}}}$

$f\left ( x_{k}|H_{1} \right )=\frac{1}{\sqrt{2\pi \sigma }}e^{-\frac{\left ( x_{k} -s_{k}\right )^{2}}{2\sigma ^{2}}}$

由于不同采样点处的随机变量 $x_{k}$ 相互独立，考虑 $0\leq t\leq T$ 时间内的采样点数为N，则N个采样样本的联合概率密度为，即似然函数：

$f\left ( x_{1},x_{2},......, x_{N}|H_{0} \right )=\left ( \frac{1}{\sqrt{2\pi \sigma }} \right )^{N}e^{-\frac{1}{N_{0}}\sum_{k=1}^{N}x^{2}_{k}\Delta t}\; \; \; \; \; \; \; \; \; \; (1)$

$f\left ( x_{1},x_{2},......, x_{N}|H_{1} \right )=\left ( \frac{1}{\sqrt{2\pi \sigma }} \right )^{N}e^{-\frac{1}{N_{0}}\sum_{k=1}^{N}\left ( x_{k} -s_{k}\right )^{2}\Delta t}\; \; \; \; \; \; \; \; \; \; \left ( 2 \right )$

其中 $N_{0}$ 噪声功率谱密度（单边谱）， $\Delta t$ 为采样时间间隔， $\sigma ^{2}$ 为噪声功率。当雷达接收到某个信号，并得到对应的采样值，将采样值带入上述公式，可以得到目标存在以及目标不存在两种情况下，雷达接收到此信号的概率 $f\left ( x_{1},x_{2},......, x_{N}|H_{0} \right )$ 以及 $f\left ( x_{1},x_{2},......, x_{N}|H_{1} \right )$ 。由此，可以根据哪种情况下，雷达接收到此信号的概率最大，判断目标的有无。即：

$f\left ( x_{1},x_{2},......, x_{N}|H_{1} \right )\overset{\overset{H_{0}}{<}}{\underset{H_{1}}{>}}f\left ( x_{1},x_{2},......, x_{N}|H_{0} \right )$

这种判断目标有无的方法（准则）称为最大似然准则，将上式变形得到：

$l\left ( x \right )=\frac{f\left ( x_{1},x_{2},......, x_{N}|H_{1} \right )}{f\left ( x_{1},x_{2},......, x_{N}|H_{0} \right )}\overset{\overset{H_{0}}{<}}{\underset{H_{1}}{>}}l_{0}=1\; \; \; \; \; \; \; \; \; \; \left ( 3 \right )$

此时判决表达式左侧式子 $l\left ( x \right )$ 称为似然比，是雷达接收信号 $x=\left [ x_{1},x_{2},......, x_{N} \right ]$ 的函数。右侧数值 $l_{0}$ 是判决统计量 $l\left ( x \right )$ 在最大似然准则下的检测门限。

事实上，不同的检测准则，其检测统计量一般相同，为似然比 $l\left ( x \right )$ ，唯一不同的是检测门限 $l_{0}$ 。下面是常见检测准则下的检测门限值。

最大似然检测下

$l_{0}=1$

最大后验概率准则下

$l_{0}=\frac{P\left ( H_{0} \right )}{P\left ( H_{1} \right )}$

其中 $P\left ( H_{1}\right ),P\left ( H_{0}\right )$ 分别为目标存在与否的先验概率；

最小错误准则下

$l_{0}=\frac{P\left ( H_{0} \right )}{P\left ( H_{1} \right )}$

最小风险Bayes准则下

$l_{0}=\frac{\left ( C_{10}-C_{00} \right )P\left ( H_{0} \right )}{\left ( C_{01}-C_{11} \right )P\left ( H_{1} \right )}$

其中 $C_{00},C_{01},C_{10},C_{11}$ 分别对应的代价因子，如 $C_{10}$ 为 $H_{0}$ 情况下判为 $H_{1}$

奈曼-皮尔逊准则下

$l_{0}=\lambda \left ( P_{f} \right )$

极大极小准则下

$l_{0}=\frac{\left ( C_{10}-C_{00} \right )q_{0}}{\left ( C_{01}-C_{11} \right )\left ( 1- q_{0}\right )$

其中， $q_{0}$ 满足：

$C_{00}+\left ( C_{10}-C_{00}\right )\alpha \left ( q_{0} \right )=C_{11}+\left ( C_{01}-C_{11}\right )\beta \left ( q_{0} \right )$

而：

$\alpha \left ( q \right )=P\left ( D_{1}|H_{0} \right );\beta \left ( q \right )=P\left ( D_{0}|H_{1} \right )$

将式（1），（2）代入（3）化简可得新的判决表达式：

$I\left ( x \right )=\sum_{k=1}^{N}x_{k}s_{k}\Delta t\overset{\overset{H_{0}}{<}}{\underset{H_{1}}{>}}\beta=\frac{N_{0}}{2}lnl_{0}+\frac{E}{2}\; \; \; \; \; \; \; \; \; \; \left ( 4 \right )$

其中

$E=\sum_{k=1}^{N}x^{2}_{k}\Delta t$

考虑到 $x_{k}$ 为高斯变量，则判决统计量 $I=I\left ( x \right )$ 为高斯变量，并且可以获得两种情况下 $I$ 的概率密度函数：

$f\left ( I|H_{0} \right )=\frac{1}{\sqrt{2\pi \sigma_{I} }}e^{-\frac{I^{2}}{2\sigma_{I} ^{2}}}$

$f\left ( I|H_{1} \right )=\frac{1}{\sqrt{2\pi \sigma_{I} }}e^{-\frac{\left ( I-E \right )^{2}}{2\sigma_{I} ^{2}}}$

其中

$\sigma _{I}^{2}=\frac{N_{0}E}{2}$

由此获得虚警概率：

$P_{f}=P\left ( D_{1}|H_{0} \right )=\int_{\beta }^{\infty }f\left ( I|H_{0} \right )dI=\frac{1}{2}\left [ 1-erf\left ( \left ( \frac{\beta }{E} \right )\sqrt{\frac{E}{N_{0}}} \right ) \right ]$

检测概率：

$P_{d}=P\left ( D_{1}|H_{1} \right )=\int_{\beta }^{\infty }f\left ( I|H_{1} \right )dI=\frac{1}{2}\left [ 1-erf\left ( \left ( \frac{\beta }{E} -1\right )\sqrt{\frac{E}{N_{0}}} \right ) \right ]$

2. 代码实现

考虑雷达发射信号是chirp信号，信号模型为：

$s\left ( t \right )=cos\left ( \pi Kt^{2} \right )\; \; \; \; \; -\frac{T}{2}\leq t\leq \frac{T}{2}$

信号参数见仿真代码。考虑奈曼-皮尔逊检测准则（恒虚警检测），根据理论虚警概率计算表达式由虚警概率计算得到检测门限，通过对比每次蒙特卡洛实验下的检测统计量与检测门限，仿真目标检测场景，统计虚警概率、漏警概率，并与对应的理论值对比，代码如下：

clear all;
close all;
clc;
%% 发射信号参数（以线性调频信号为例）
B = 100e6;              % 发射信号带宽
T = 1e-6;               % 发射信号时宽
K = B/T;                % 发射信号调频率
c = 3e8;                % 光速

%% 回波信号（不考虑时延以及信号能量衰减）
fs = 500e6;                % 采样率
dt = 1/fs;
t = -T/2:dt:T/2;      
st = cos(pi*K*t.^2);   

figure;
plot(st)
title('信号波形');

%% 噪声功率
E = st*st'*dt;
SNR = 10;                  % 信噪比,单位dB    
N0 = E/(10^(SNR/10));
sigma2 = N0/2*fs;
Monte = 100000;

%% 门限(恒虚警检测）

% 设置的虚警概率
pf = 1e-3;
% 检测门限
beta = erfinv(1-2*pf)*sqrt(E*N0);
% 理论计算虚警概率
pf1 = 1/2*(1-erf(beta/sqrt(E*N0)));
% 理论计算漏警概率
pm = 1/2*(1+erf((beta-E)/sqrt(E*N0)));
% 理论计算检测概率
pd = 1-pm;
%% 虚警概率

k01 = 0;
% 虚警概率
pf_monte = zeros(1,Monte);
for i = 1:Monte
    w = sqrt(sigma2)*randn(1,length(t));
    xt = w;
    I = xt*st'*dt;
    if I>beta
        k01=k01+1;
    end   
    pf_monte(i) = k01/i;
end


%% 检测概率
k02 = 0;

% 漏警概率
pm_monte = zeros(1,Monte);
% 检测概率
pd_monte = zeros(1,Monte);

for ii = 1:Monte
    w = sqrt(sigma2)*randn(1,length(t));
    xt = st+w;
    I = xt*st'*dt;
    if I>beta
        k02=k02+1;
    end    
    pd_monte(ii) = k02/ii;
    pm_monte(ii) = 1-k02/ii;

end

figure;
plot(pf_monte,'k');
hold on;
plot(pf*ones(1,Monte),'r');
hold off;
legend('蒙特卡罗统计值','理论值');
xlabel('蒙特卡罗实验次数');
ylabel('虚警频率');
title(['虚警概率为',num2str(pf)]);

figure;
plot(pm_monte,'k');
hold on;
plot(pm*ones(1,Monte),'r');
hold off;
legend('蒙特卡罗统计值','理论值');
xlabel('蒙特卡罗实验次数');
ylabel('漏警频率');
title(['漏警概率为',num2str(pm)]);

figure;
plot(pd_monte,'k');
hold on;
plot(pd*ones(1,Monte),'r');
hold off;
legend('蒙特卡罗统计值','理论值');
xlabel('蒙特卡罗实验次数');
ylabel('检测频率');
title(['检测概率为',num2str(pd)]);

仿真结果如下：

从仿真结果可以看出，由统计得到的虚警概率、漏警概率以及检测概率随着实验次数的增加，逐渐趋于理论值。

由此，通过检测概率计算表达式可以分析恒虚警的检测概率随信噪比的变化规律，代码如下：

clear all;
close all;
clc;
%% 发射信号参数（以线性调频信号为例）
B = 100e6;              % 发射信号带宽
T = 1e-6;               % 发射信号时宽
K = B/T;                % 发射信号调频率
c = 3e8;                % 光速

%% 回波信号（不考虑时延以及信号能量衰减）
fs = 500e6;                % 采样率
dt = 1/fs;
t = -T/2:dt:T/2;      
st = cos(pi*K*t.^2);   

figure;
plot(st)
title('信号波形');

%% 噪声功率
E = st*st'*dt;
SNR = -20:20;                  % 信噪比,单位dB    
N0 = E./(10.^(SNR/10));

%% 门限(恒虚警检测）

% 设置的虚警概率
pf= 1e-3;
% 检测门限
beta = erfinv(1-2*pf).*sqrt(E*N0);
% 理论计算漏警概率
pm = 1/2*(1+erf((beta-E)./sqrt(E*N0)));
% 理论计算检测概率
pd = 1-pm;

%%

figure;
plot(SNR,pm,'k');
hold on
plot(SNR,pd,'r');
hold off
legend('漏警概率','检测概率')
xlabel('信噪比/dB');
ylabel('漏警概率/检测概率');
title(['恒虚警检测，虚警概率为',num2str(pf)]);

仿真结果如下

由仿真结果可以看出，随着信噪比的增加，目标的检测概率逐渐增加，在特定的信噪比范围内，信噪比的增加对检测概率的增加影响很大，当信噪比增加到一定值时，检测概率很高，但继续增加信噪比，检测概率增加量很小。

二、一维距离像

1. 信号模型

当检测到雷达接收的信号存在目标回波时，下一步任务就是对目标距离的测量。雷达的测距原理是通过雷达接收的回波信号相对于雷达发射信号的时延来确定的。

假设雷达发射信号：

$s\left ( t \right )\; \; \; \; \; \; \; 0\leq t\leq T$

只考虑接收的热噪声，则距离为 $R$ 的目标回波信号模型为：

$x\left ( t \right )=as\left ( t -\frac{2R}{c}\right )+n\left ( t \right )$

为保障雷达探测威力，发射信号时宽较大，直接测量时延困难。可以通过设计大带宽信号，利用脉冲压缩原理，通过匹配滤波的方式形成时域窄脉冲，实现时延的精准估计（在噪声环境中，测距分辨率也会影响测距精度），进而得到目标与雷达的距离。所设计的匹配滤波器与雷达发射信号有关，滤波器表达式为：

$h\left ( t \right )=s^{*}\left ( -t \right )$

利用匹配滤波器 $h\left ( t \right )$ 对雷达接收信号 $x\left ( t \right )$ 匹配滤波得到：

$y\left ( t \right )=x\left ( t \right )\bigotimes h\left ( t \right )$

时域卷积运算为：将匹配滤波器 $h\left ( t \right )$ 关于 $t=0$ 的反转得到 $h_{1}\left ( t \right )=h\left ( -t \right )$ ，然后将 $h_{1}\left ( t \right )$ 滑动并与接收信号 $x\left ( t \right )$ 点乘求和运算（也可利用时域卷积对应频域相乘进行卷积的快速计算），得到不同时延 $t_{0}$ 的发射信号 $s\left ( t-t_{0} \right )$ 与接收信号 $x\left ( t \right )$ 的相关值。

当时延 $t_{0}$ 满足：

$t_{0}=\frac{2R}{c}$

不考虑信号能量的衰减，即 $a=1$ 。则 $t_{0}$ 时延下的相关值为匹配滤波输出模值最大所对应时刻，且最大值即为目标检测中的检测统计量的值

2. 代码实现

考虑雷达发射信号是chirp信号，信号模型为：

$s\left ( t \right )=exp\left ( j\pi K\left (t-\frac{T}{2} \right )^{2} \right )\; \; \; \; \; 0\leq t\leq T$

仿真参数见代码部分。代码如下：

clear all;
close all;
clc;
%% 仿真参数

% 信号参数
B = 100e6;              % 发射信号带宽
T = 1e-6;               % 发射信号时宽
K = B/T;                % 发射信号调频率
c = 3e8;                % 光速

% 目标参数
R1 = 750;               % 目标1距离
R2 = 1500;              % 目标2距离

% 雷达探测参数
Rmin = c*T/2;           % 由信号脉宽决定
Rmax = 2000;            % 由雷达探测威力确定

%% 发射信号
fs = 500e6;                        % 采样率
dt = 1/fs;
tt = 0:1/fs:T;         
st = exp(1j*pi*K*(tt-T/2).^2);     % 发射信号
Ns = length(st);

figure;
plot(real(st))
title('发射信号波形')
%% 噪声功率
E = st*st'*dt;
SNR = 20;                        % 信噪比,单位dB    
N0 = E/(10^(SNR/10));
siama2 = N0/2*fs;                % 噪声功率

%% 接收信号
Tmin = 2*Rmin/c;
Tmax = 2*Rmax/c;
t = Tmin:1/fs:Tmax;    
xt = zeros(1,length(t));
t1 = round((2*R1/c-Tmin)*fs)+1;
t2 = round((2*R2/c-Tmin)*fs)+1;
xt(t1:t1+Ns-1) = xt(t1:t1+Ns-1)+st;
xt(t2:t2+Ns-1) = xt(t2:t2+Ns-1)+st;
w = sqrt(siama2)*randn(1,length(t));
xt = xt+w;

figure;
plot(t,real(xt));
title('接收信号波形')

%% 匹配滤波

ht = conj(fliplr(st));
r_mf = conv(xt,ht);
r_mf = r_mf(Ns:Ns+length(t)-1);

figure;
plot(c*t/2,abs(r_mf));
xlabel('距离/m')
ylabel('匹配输出信号强度');
title('一维距离像');

仿真结果：

三、二维距离像

1. 方位搜索

一维距离像能够显示探测区域目标的距离分布，但由于目标的方位未知，仅根据估计的距离并不能实现对目标位置的估计，这一节引入方位向的估计并联合距离向实现目标的二维定位。代码分两个场景：自编的二维平面定位场景，方位向定位通过相控阵波束扫描得到，为实现方位向的高精度定位，一般采用空间普估计技术，典型的算法如MUSIC和ESPRIT；利用MATLAB中的phased工具箱编写的三维定位场景，为了搜索的简便，方位向的搜索是在俯仰角已知的情况下进行的（之后看情况会对phased工具箱做一期介绍）。

波束形成

波束形成器是滤波器的空间形式，它用来接收沿空间某一方向传播的期望信号，同时抑制来自其它方向的干扰和随机噪声。

常规波束形成

如图是波束形成器的一般结果，信号输出为：

$y = W^{H}\left ( \theta \right )X$

其中

$y=\begin{bmatrix} y\left ( 0 \right ) & y\left ( 1 \right )& \cdots & y\left ( N-1 \right ) \end{bmatrix}$

$W\left ( \theta \right )=\begin{bmatrix} 1 & e^{-j2\pi d/\lambda sin\vartheta }& \cdots & e^{-j2\pi (M-1)d/\lambda sin\vartheta }\end{bmatrix}^{T}$

$X=\begin{bmatrix} x_{1}& x_{2}& \cdots & x_{M} \end{bmatrix}^{T}$

$x_{i}=\begin{bmatrix} x_{i}\left ( 0 \right )& x_{i}\left ( 1 \right )& \cdots & x_{i}\left ( N-1\right ) \end{bmatrix}^{T},i=1,2,\cdots ,M$

对于不同的搜索方向矢量 $a^{H}\left ( \theta \right )=W^{H}\left ( \theta \right )$ ,则不同角度下的输出功率为：

$P_{CBF}\left ( \theta \right )=a^{H}\left ( \theta \right )\left ( XX^{H} \right )a^{H}\left ( \theta \right )$

根据输出功率最大准则，此时最佳空间滤波器为：

${W}_{opt}=\underset{a\left ( \theta \right )}{max} \, \, P_{CBF}\left ( \theta \right )$

Capon波束形成器

常规波束形成虽然空间滤波输出信号功率最大，但信噪比不一定是最大的，由此引入线性约束最小方差准则（LCMV），优化方程如下：

$\left\{\begin{matrix} \underset{W}{min}P=\, \, W^{H}RW\\ W^{H}a\left ( \theta \right )=1 \end{matrix}\right.$

其中 $R$ 为不含期望信号的阵列协方差矩阵。由此得到此准则下的最佳空间滤波器：

${W}_{opt}=\frac{R^{-1}a\left ( \theta \right )}{a^{H}\left ( \theta \right )R^{-1}a\left ( \theta \right )}$

对于的阵列方向响应，或者波束形成图为：

$g\left ( \theta \right )={W}^{H}_{opt}a\left ( \theta \right )$

MUSIC多重信号分类法

基本思路：利用信号子空间和噪声子空间的正交性来估计信号的方位，即利用不同信号源之间，信号源与噪声之间的弱相干性。

窄带远场信号模型：

$X = AS+n$

其中

$X=\begin{bmatrix} x_{1}& x_{2}& \cdots & x_{M} \end{bmatrix}^{T}$

$x_{i}=\begin{bmatrix} x_{i}\left ( 0 \right )& x_{i}\left ( 1 \right )& \cdots & x_{i}\left ( N-1\right ) \end{bmatrix}^{T},i=1,2,\cdots ,N$

$A=\begin{bmatrix} a_{1} & a_{2}& \cdots & a_{K}\end{bmatrix}$

$a_{k}=\begin{bmatrix} 1 & e^{-j2\pi d/\lambda sin\vartheta_{k} }& \cdots & e^{-j2\pi (M-1)d/\lambda sin\vartheta_{k} }\end{bmatrix}^{T},k=1,2,\cdots ,K$

$S=\begin{bmatrix} s_{1}& s_{2}& \cdots & x_{K} \end{bmatrix}^{T}$

$s_{k}=\begin{bmatrix} s_{k}\left ( 0 \right )& s_{k}\left ( 1 \right )& \cdots & s_{k}\left ( N-1\right ) \end{bmatrix}^{T},k=1,2,\cdots ,K$

M为阵元数目，K为信号源数目，N每个阵元时域采样点数。数据 $X$ 的相关矩阵为：

$R=E\left [ XX^{H} \right ]=AR_{s}A+\sigma ^{2}I$

对 $R$ 特征分解，其中K个最大特征值对应的特征矢量构成的是信号子空间，其他特征矢量构成的子空间为噪声子空间 $V$ 。由此搜索矢量 $a\left ( \theta \right )$ 对于噪声子空间每个特征矢量正交，由此可以得到如下谱峰：

$S\left ( \theta \right )=\frac{1}{a^{H}\left ( \theta \right )VV^{H}a\left ( \theta \right )}$

由此得到信号源的方位估计。

2. 代码实现

常规波束形成与MUSIC对比

利用相控阵天线方位扫描特性对目标方位进行搜索，并将波束形成方法与MUSIC方法结果进行对比，代码如下：

clc; close all;clear; 

%%
% 相控阵天线扫描
%% 仿真参数

c = 3e8;
fc = 300e6;                                  % 载频
lambda = c/fc;                               % 载波波长
d = 0.5*lambda;                              % 阵元空间间距
M = 18;                                      % 阵元数目（天线通道数）
N = 100;                                     % 采样点数
thetas1 = 30;                                % 信号1入射角
Px1 = 1;
thetas2 = 50;                                % 信号2入射角度 
Px2 = 1;
m = 0:M-1;
fs = 1e9;                                    % 采样频率
t = (0:1:N-1)/fs; 
% 信号生成
vs1 = exp(-1j*2*pi*m'*d/lambda*sind(thetas1)); % 信号1方向矢量 
vs2 = exp(-1j*2*pi*m'*d/lambda*sind(thetas2)); % 信号2方向矢量 
x1 = sqrt(Px1)*vs1*(randn(1,N).*exp(1j*2*pi*fc*t));                    
x2 = sqrt(Px2)*vs2*(randn(1,N).*exp(1j*2*pi*fc*t));                     
x = x1+x2;                                     % M个阵元接收信号

snr = 10;
x=awgn(x,snr,'measured');
%% 波束形成
Rx = x*x'/N;                                 % 协方差矩阵
sita = 80*(-1:0.001:1);                      % 扫描方向范围
v = exp(-1j*2*pi*m'*d/lambda*sind(sita));    % 扫描导向矢量
Py = abs(diag(v'*Rx*v));                     % 各扫描方向接收的信号功率

figure;
plot(sita,10*log10(Py/max(Py)),'k'); 
xlabel('角度/degree');
ylabel('功率/dB'); 
title('波束形成方位向扫描');
grid on 
axis([-90 90 -50 0]); 
hold off

%% MUSIC

[EV,D]=eig(Rx);%%%%  按照 特征值 分解  
theta = -90:0.5:90;
SP = zeros(1,length(theta));
Num = 2; 
for i = 1:length(theta)
        a=exp(-1j*2*pi*m'*d/lambda*sind(theta(i)));
        En=EV(:,1:M-Num); 
        an = a'*En;
        SP(i)=1/(an*an');
end
   
SP=abs(SP);
SPmax=max(SP);
SP=10*log10(SP/SPmax);

figure;
plot(theta,SP);
xlabel('角度 (degree)')
ylabel('幅度 (dB)')
title('MUSIC方位估计');
axis([-90 90 -60 0])
grid on

仿真结果：

Capon波束形成器

上述波束形成的加权系数简单，并未考虑噪声的影响，下面是线性约束最小方差准则波束形成代码：

clc; close all;clear; 

%%
% 相控阵天线扫描
%% 仿真参数

c = 3e8;
fc = 300e6;                                  % 载频
lambda = c/fc;                               % 载波波长
d = 0.5*lambda;                              % 阵元空间间距
M = 18;                                      % 阵元数目（天线通道数）
N = 100;                                     % 采样点数
thetas = 30;                                % 信号1入射角
m = 0:M-1;
fs = 1e9;                                    % 采样频率
t = (0:1:N-1)/fs; 
% 信号生成
vs = exp(-1j*2*pi*m'*d/lambda*sind(thetas)); % 信号1方向矢量 
x = vs*(randn(1,N).*exp(1j*2*pi*fc*t));                    
snr = 10;
x=awgn(x,snr,'measured');
%% 波束形成
Rn = cov(x.');                                  % 噪声协方差矩阵
wop1=(Rn\vs)/(vs'*(Rn\vs));                    
sita=80*(-1:0.001:1);                           % 扫描方向范围
v=exp(-1j*2*pi*m'*d/lambda*sind(sita));         % 扫描方向矢量 
B=abs(wop1'*v); 
plot(sita,20*log10(B/max(B)),'k'); 
title('波束图');xlabel('角度/degree');ylabel('波束图/dB'); 
grid on 
axis([-90 90 -50 0]); 
hold off

仿真结果：

利用MATLAB中phased工具实现二维距离成像

考虑雷达发射信号是chirp信号，基带信号模型为：

$s\left ( t \right )=exp\left ( j\pi K\left (t-\frac{T}{2} \right )^{2} \right )\; \; \; \; \; 0\leq t\leq T$

仿真参数见代码部分，代码如下：

clear;
close all;
clc;
%%
%
%--------------------------------回波信号生成-------------------------
%% 参数设置
fc = 300e6;
c = 3e8;
lambda = c/fc;
% 目标参数
TgtModel = phased.RadarTarget;
tgaz = 50;                                      %% 目标方位
tgel = 0;                                       %% 目标俯仰
Rtg = 30e3;                                     %% 目标距离
vtg = 100;                                      %% 雷达速度
tgtpos = Rtg*[cosd(tgel)*cosd(tgaz);cosd(tgel)*sind(tgaz);sind(tgel)];     
tgtvel = vtg*[cosd(tgel)*cosd(tgaz);cosd(tgel)*sind(tgaz);sind(tgel)];      

% 阵列天线参数
antenna = phased.ULA;
antenna.NumElements = 16;
antenna.ElementSpacing = 0.3*lambda;
cosineElement = phased.CosineAntennaElement;
antenna.Element = cosineElement;
% 波形参数
waveform = phased.LinearFMWaveform;
waveform.PRF = 1000;
waveform.PulseWidth = 1e-4;
prf = waveform.PRF;
nSamples = waveform.SampleRate/prf;
% 发射机参数
TX = phased.Transmitter;
TX.Gain = 20;
% 平台参数
PlatformModel = phased.Platform;
PlatformModel.InitialPosition = tgtpos;
PlatformModel.Velocity = tgtvel;

% 信道参数
ChannelModel = phased.FreeSpace;
ChannelModel.TwoWayPropagation = true;

% 发射与接收参数

txArray = phased.Radiator;
txArray.Sensor = antenna;
txArray.OperatingFrequency = fc;
rxArray = phased.Collector;
rxArray.Sensor = antenna;
rxArray.OperatingFrequency = fc;
rxPreamp = phased.ReceiverPreamp;
rxPreamp.Gain = 10;
rxPreamp.NoiseFigure = 5;
% 变量设置
radarPos = [0;0;0];
radarVel = [0;0;0];

nPulses = 32;
tgtAng = zeros(2,nPulses);
tgtAngcopy = zeros(2,nPulses);
data = complex(zeros(nSamples,antenna.NumElements,nPulses));


%% 雷达脉冲生成
for ii=1:nPulses
    wf = step(waveform);                                               % 生成波形
    [tgtPos, tgtVel] = step(PlatformModel,1/prf);                      % 更新目标位置
    [tgtRng, tgtAng] = rangeangle(tgtPos, radarPos);                   % 计算对应的距离与角度
    tgtAngcopy(:,ii) = tgtAng;
    s0 = step(TX, wf);                                                 % 信号发射增益
    s1 = step(txArray,s0, tgtAng);                                     % 从阵列发射信号
    s2 = step(ChannelModel, s1, radarPos, tgtPos, radarVel, tgtVel);   % 雷达与目标之间信号往返传播
    s3 = step(TgtModel, s2);                                           % 信号从目标反射
    s4 = step(rxArray,s3,tgtAng);                                      % 阵列接收信号
    s5 = step(rxPreamp,s4);                                            % 加入接收机噪声
    data(:,:,ii) = s5(:,:);                                            % 信号逐脉冲生成
end
t = (0:nPulses*nSamples-1)/waveform.SampleRate;
y = abs(data(:,1,:));

figure;
plot(t,y(:));
xlabel('时间/s');
ylabel('幅度');
title('回波信号(阵元1)');

%--------------------------------波束形成-------------------------
%% 波束形成
% 波束形成参数
beamformer = phased.PhaseShiftBeamformer;
beamformer.SensorArray = antenna;
beamformer.DirectionSource = 'Input port';
beamformer.WeightsOutputPort = true;
beamformer.WeightsNormalization = 'Preserve power';


[bf0,w0] = beamformer(data(:,:,1),[0;0]);
fieldval0=pattern(antenna,fc,-180:180,0,'Weights',w0,...
    'PropagationSpeed',physconst('LightSpeed'),'Normalize',false,...
    'Type','powerdb','CoordinateSystem','rectangular');

[bfaz, waz] = step(beamformer,data(:,:,1),[tgaz;tgel]);
fieldvalaz  = pattern(antenna,fc,-180:180,tgel,'Weights',waz,...
    'PropagationSpeed',physconst('LightSpeed'),'Normalize',false,...
    'Type','powerdb','CoordinateSystem','rectangular');
figure;
pattern(antenna,fc,'Weights',w0);
title('朝（0，0）方向的波束合成方向图')
figure;
pattern(antenna,fc,'Weights',waz);
title(['朝（',num2str(tgaz),'，',num2str(tgel),'）方向的波束合成方向图']);
figure;
subplot(2,2,1);
plot(-180:180,fieldval0);
xlabel('方位角（度）');
ylabel('功率（dB）');
title(' (0,0)合成方向方位截面');

subplot(2,2,2);
plot(abs(bf0));
xlabel('时间 /ms');
ylabel('幅度');
title(' (0,0)方向波束形成回波');

subplot(2,2,3);
plot(-180:180,fieldvalaz);
xlabel('方位角（度）');
ylabel('功率（dB）');
title(['（',num2str(tgaz),'，',num2str(tgel),'）合成方向方位截面 ']);

subplot(2,2,4);plot(abs(bfaz));
xlabel('时间 /ms');
ylabel('幅度');
title(['（',num2str(tgaz),'，',num2str(tgel),'）方向波束形成回波 ']);

%% 所有脉冲波束形成

beamformed=complex(zeros(nSamples,nPulses));
for ii=1:nPulses
    beamformed(:,ii)=beamformer(data(:,:,ii),tgtAngcopy(:,ii));
end
t = (0:nPulses*nSamples-1)/waveform.SampleRate;
y = abs(reshape(beamformed,nPulses*nSamples,1));
dc=abs(reshape(data(:,1,:),nPulses*nSamples,1));

figure;
subplot(2,1,1);
plot(t,dc);
xlabel('时间/s');
ylabel('幅度');
title('目标回波');

subplot(2,1,2);
plot(t,y);
xlabel('时间/s');
ylabel('幅度');
title('超目标方向波束形成回波');

%--------------------------------参数估计-------------------------
%% 匹配滤波
b = getMatchedFilter(waveform);
matchedfilter = phased.MatchedFilter(...
    'Coefficients',b,...
    'SpectrumWindow','Hamming');

matchFiltered = step(matchedfilter,beamformed);
[m,ind] = max(abs(matchFiltered(:,nPulses/2)));  % ind为中间脉冲脉冲压缩峰值对应的距离单元
t = (0:nSamples-1)/waveform.SampleRate;

figure;
subplot(1,2,1);
plot(real(wf));
xlabel('时间/s');
ylabel('幅度');
title('雷达波形');
subplot(1,2,2);
plot(real(b));
xlabel('时间/s');
ylabel('幅度');
title('匹配滤波器');

figure;
subplot(1,2,1);
plot(t,abs(bfaz));
xlabel('时间 /ms');
ylabel('幅度')
title(['（',num2str(tgaz),'，',num2str(tgel),'）方向波束形成回波 ']);
subplot(1,2,2);
plot(t,abs(matchFiltered(:,nPulses/2)));
xlabel('时间/s');
ylabel('幅度');
title('回波的脉冲压缩输出');

%% 利用 RangeDopplerResponse进行多普勒处理
rangeDoppler=phased.RangeDopplerResponse('PropagationSpeed',physconst('LightSpeed'),...
    'SampleRate',waveform.SampleRate,...
    'DopplerFFTLengthSource','Property',...
    'DopplerFFTLength',512, ...
    'DopplerWindow','Kaiser',...
    'DopplerSidelobeAttenuation',30,...
    'DopplerOutput','Speed', ...
    'OperatingFrequency',fc);
[Rv,range,vecter] = rangeDoppler(beamformed,b);

figure;
mesh(vecter,range,abs(Rv));
xlabel('速度 (m/s)');
ylabel('距离 (m)');
title('距离-速度图像');

[~,Irv] = max(abs(Rv(:)));
Rtg_est = range(mod(Irv-1,length(range))+1);
vtg_est = vecter(ceil(Irv/length(range)));

fprintf('距离估计Rtg=%.6fm\n',Rtg_est);
fprintf('速度估计vtg=%.6fm/s\n',vtg_est);

%%  利用 AngleDopplerResponse 进行多普勒处理

matched = complex(zeros(nSamples,nPulses,antenna.NumElements));
for ii=1:antenna.NumElements
    for jj=1:nPulses
        matched(:,jj,ii) = filter(b,1,data(:,ii,jj));
    end
end

dd = squeeze(matched(ind,:,:)).';
angleDopplerResp = phased.AngleDopplerResponse('SensorArray',antenna,...
    'OperatingFrequency',3e8, ...
    'PropagationSpeed',physconst('LightSpeed'),...
    'PRF',prf, 'ElevationAngle',tgel,...
    'NumAngleSamples',512,'NumDopplerSamples',512);

[AnDo,angle,dop] = angleDopplerResp(dd);
figure;
mesh(angle,dop,abs(AnDo));
xlabel('角度 (度)');
ylabel('多普勒频率 (Hz)');
title('方位-多普勒图像');

[~,Iad] = max(abs(AnDo(:)));
fd_est = dop(mod(Iad-1,length(dop))+1);
tgaz_est = angle(ceil(Iad/length(dop)));
fprintf('多普勒估计fd=%.6fHz\n',fd_est);
fprintf('方位估计tgaz=%.6f度\n',tgaz_est);

%% 定位
figure;
polar(tgaz*pi/180,Rtg,'o');
hold on;
polar(tgaz_est*pi/180,Rtg_est,'x');
hold off;

仿真结果：

估计结果：

注：估计的速度为负值的原因：目标初始速度为正（目标远离雷达的速度为正），因此最后得到的多普勒频率为负值。由于phased工具向箱中速度与频率的转化关于速度的方向与场景设置的速度方向不一致，因此得到的速度估计与真实值差个负号。

四、SAR成像

SAR的公式推导较为复杂，这节内容不做详细分析，之后看情况会出一期SAR公式详细推导的内容，这节重点是SAR代码的编写，代码部分注释比较详细，对照SAR相关书籍不难看懂。下面简单介绍RD（距离多普勒）算法的主要思路以及步骤。

1. RD算法

一维距离像介绍了通过匹配滤波器实现距离维的脉冲成像，这部分原理与SAR成像中的距离维成像相对应；波束形成原理本质上可以理解为利用不同方位信号的弱相干实现方位对位的，这部分与SAR成像的方位维成像对应，与波束形成不同的是，SAR是通过平台的移动得到不同方位的弱相干信号的。

上图是仿真的场景图，仿真的斜视角为零度，所用方法为RD算法。之所以利用RD算法由于相同距离不同方位上的目标多普勒历程一致，只是时延不一样，所以可以在距离多普勒域统一处理相同距离处目标的成像，提升计算效率。基本上，包含RD算法在内的SRA成像算法的主要步骤包括：距离维脉冲压缩；距离徙动校正；方位向脉冲压缩形成SAR图像。根据SAR成像精度的需要可能需要相位补偿的步骤。

本节代码数据是根据参数仿真生成的，SAR成像算法主要考虑了距离维脉冲压缩和方位维脉冲压缩。没有考虑距离徙动校正的原因：仿真目标为点目标，徙动现象不是很明显。

2. 代码实现

RD算法代码如下：

%% RD算法
clear;close all;clc;

%% 参数设置
% 载频信号参数
c = 3e8;
fc = 1e9;                               % 信号载频
lambda = c/fc;                          % 载波波长

% 探测范围(地面范围）
Az0 = 10e3;
AL = 1000;
Azmin = Az0-AL/2;                       % 方位向范围
Azmax = Az0+AL/2;                  
Rg0 = 10e3;                             % 中心地距
RgL = 1000;                             % 测绘带宽

% 平台参数
vr = 100;                               % SAR搭载平台速度
H = 5000;                               % 平台高度
R0 = sqrt(Rg0^2+H^2);                   % 中心斜距（斜视角零度）

%天线参数
D = 4;                                  % 方位向天线长度
La = lambda*R0/D;                       % 合成孔径长度
Ta = La/vr;                             % 合成时长

% 方位维/慢时间维参数
Ka = -2*vr^2/lambda/R0;                 % 慢时间维调频率
Ba = abs(Ka*Ta);                        % 慢时间维带宽
PRF = 1.2*Ba;                           % 脉冲重复频率
Mslow = ceil((Azmax-Azmin+La)/vr*PRF);  % 慢时间维点数/脉冲数
Mslow = 2^nextpow2(Mslow);              % 用于慢时间维FFT的点数
ta = linspace((Azmin-La/2)/vr,(Azmax+La/2)/vr,Mslow);  
PRF = 1/((Azmax-Azmin+La)/vr/Mslow);    % 与慢时间维FFT点数相符的脉冲重复频率

% 距离维/快时间维参数
Tw = 5e-6;                              % 脉冲持续时间
Br = 30e6;                              % 发射信号带宽
Kr = Br/Tw;                             % 调频率
Fr = 2*Br;                              % 快时间维采样频率
Rmin = sqrt((Rg0-RgL/2)^2+H^2);
Rmax = sqrt((Rg0+RgL/2)^2+H^2+(La/2)^2);
Nfast = ceil(2*(Rmax-Rmin)/c*Fr+Tw*Fr); % 快时间维点数
Nfast = 2^nextpow2(Nfast);              % 用于快时间维FFT的点数
tr = linspace(2*Rmin/c,2*Rmax/c+Tw,Nfast);  
Fr = 1/((2*Rmax/c+Tw-2*Rmin/c)/Nfast);  % 与快时间维FFT点数相符的采样率

% 分辨率参数
Dr = c/2/Br;                            % 距离分辨率
Da = D/2;                               % 方位分辨率

% 点目标参数
Ntarget = 4;                            % 目标数量
Ptarget = [Az0-10,Rg0-20,1              % 目标位置\散射信息
         Az0+10,Rg0-10,1
         Az0-10,Rg0+10,1
         Az0+10,Rg0+20,1];  
fprintf('仿真参数：\n');     
fprintf('快时间/距离维过采样率：%.4f\n',Fr/Br);     
fprintf('快时间/距离维采样点数：%d\n',Nfast);     
fprintf('慢时间/方位维过采样率：%.4f\n',PRF/Ba);     
fprintf('慢时间/方位维采样点数：%d\n',Mslow);     
fprintf('距离分辨率：%.1fm\n',Dr);     
fprintf('距离横向分辨率：%.1fm\n',Da);     
fprintf('合成孔径长度：%.1fm\n',La);     
disp('目标方位/地距/斜距：');
disp([Ptarget(:,1),Ptarget(:,2),sqrt(Ptarget(:,2).^2+H^2)])

%% 探测范围离散化
R = Rmin:1/Fr*c/2:Rmax;
Nr = length(R);
Az = Azmin:vr/PRF:Azmax;
Na = length(Az);

%% 回波信号生成
snr = 0;                                % 信噪比
Srnm = zeros(Mslow,Nfast);
for k = 1:1:Ntarget
    sigmak = Ptarget(k,3);
    Azk = ta*vr-Ptarget(k,1);
    Rk = sqrt(Azk.^2+Ptarget(k,2)^2+H^2);
    tauk = 2*Rk/c;
    tk = ones(Mslow,1)*tr-tauk'*ones(1,Nfast);
    phasek = pi*Kr*tk.^2-(4*pi/lambda)*(Rk'*ones(1,Nfast));
    Srnm = Srnm+sigmak*exp(1i*phasek).*(0<tk&tk<Tw).*((abs(Azk)<La/2)'*ones(1,Nfast));
end                                
Srnm = awgn(Srnm,snr,'measured');

%% 距离压缩
thr = tr-2*Rmin/c;
hrc = exp(1i*pi*Kr*thr.^2).*(0<thr&thr<Tw);  % 距离维匹配滤波器
SARr =ifft((fft(Srnm,Nfast,2).*(ones(Mslow,1)*conj(fft(hrc,Nfast,2)))),Nfast,2);

figure;
imagesc(255-abs(SARr));                       
xlabel('快时间');
ylabel('慢时间');
title('快时间维脉压');
colormap(gray)

figure;
Sta = round(Ta/2*PRF);
imagesc(R,Az,255-abs(SARr(Sta:Sta+Na-1,1:Nr)));                       
xlabel('\rightarrow\it斜距/m');
ylabel('\it方位/m\leftarrow');
title('距离压缩');
colormap(gray);

figure;
col = round(linspace(1,Na,7));row =round(Nr/2)-200:round(Nr/2)+200;
waterfall(R(row),Az(col),abs(SARr(col+Sta-1,row)));
xlabel('斜距/m');
ylabel('方位/m');
title('不同方位处的距离脉压结果');
colormap(gray);

%% 方位压缩
tha = ta-Azmin/vr;
hac = exp(1i*pi*Ka*tha.^2).*(abs(tha)<Ta/2);  % 方位维匹配滤波器
SARra = ifft(fft(SARr).*(conj(fft(hac)).'*ones(1,Nfast)));
SAR = abs(SARra);

figure;
imagesc(255-SAR);                       
xlabel('快时间');
ylabel('慢时间');
title('快时间-慢时间维脉压');
colormap(gray)

figure;
imagesc(R,Az,255-SAR(1:Na,1:Nr));      
xlabel('\rightarrow\it斜距/m');
ylabel('\it方位/m\leftarrow');
title('二维压缩后SAR图像');
colormap(gray);

figure;
mesh(R,Az,SAR(1:Na,1:Nr));
xlabel('斜距');
ylabel('方位');
title('SAR图像（数据获取面）'),

%% 3dB目标图像/地距转换器
Max = max(max(SAR));
figure;
col1 = round(Na/2)-20:round(Na/2)+20;
row1 = round(Nr/2)-20:round(Nr/2)+20;
contourf(R(row1),Az(col1),SAR(col1,row1),[0.707*Max,Max],'b');grid on
xlabel('\rightarrow\it斜距/m');
ylabel('\it方位/m\leftarrow');
title('数据获取面目标分辨率');
colormap(gray);

%% 地距转换
[X,Y] = meshgrid(R,Az);
Rg = linspace(Rg0-RgL/2,Rg0+RgL/2,Nr);
R1 = sqrt(Rg.^2+H^2);
[Xq,Yq] = meshgrid(R1,Az);
SAR = interp2(X,Y,SAR(1:Na,1:Nr),Xq,Yq,'linear');
Max1 = max(max(SAR));

figure;
mesh(Rg,Az,SAR(1:Na,1:Nr));
xlabel('地距');
ylabel('方位');
title('SAR图像（投影到地面）'),

figure;
contourf(Rg(row1),Az(col1),SAR(col1,row1),[0.707*Max1,Max1],'b');grid on
xlabel('\rightarrow\it地距/m');
ylabel('\it方位/m\leftarrow');
title('投影到地面目标分辨率');
colormap(gray);

仿真结果

总结

这部分内容比较简单的介绍了雷达领域中个人认为比较基础的几个仿真问题，希望对初次接触雷达领域的人有所帮助。关于SAR的部分，这部分内容介绍比较简单，代码部分并涉及距离徙动校正，残余相位校正，后续准备充当后，再出一期关于SAR比较详细的介绍。如果内容上有所问题，希望留言指出，如果有什么建议或意见，也可以评论区留言。