本文链接：https://blog.csdn.net/weixin_39871649/article/details/111275540

本文介绍了MUSIC算法在生命体微动模型和雷达探测技术中的应用，通过仿真验证展示其在处理微弱回波信号、识别生命体心跳和呼吸频率方面的优势，对比了与傅立叶变换的区别，并提供了MATLAB实现代码链接。

摘要由CSDN通过智能技术生成

近些年来，地震、海啸等自然灾害频发，恐怖活动也时有发生。这些自然和人为的灾难往往会造成大量人员被困于复杂而封闭的环境中，因此，如何采取快速有效的探测手段迅速确定被困人员的位置和状态成为了生命探测领域中的关键问题。根据使用的传感器类型，生命探测技术可以分为音频探测技术、红外探测技术、光学探测技术、雷达探测技术^[1]。根据探测设备是否发射信号，生命探测技术又可分为有源和无源探测技术。雷达探测技术是一种最经典的有源探测技术，其向目标区域发射电磁波，根据目标体的反射信号特征来判断区域中是否存在生命体，如心跳(1~1.5 Hz)和呼吸(0.27~0.33 Hz)的频率特征。音频探测技术、红外探测技术、光学探测技术等则是无源探测技术。对于单频雷达波，主要是采取一些时频分析的方法分析回波信号的频率成分以确定微弱生命体的有无，常用的时频分析方法包括傅立叶谱分析、多重自相关、Gabor变换、SPWVD算法、小波变换、MUSIC(Multiple Signal Classification，简称MUSIC)算法^[2]、MVDR(Minimum Variance Distortionless Response，简称MVDR)算法等。在下文中，本文将主要介绍生命体微动模型、MUSIC回波处理算法和仿真验证，并在文末给出MUSIC算法的matlab实现代码。

生命体微动模型

根据雷达的工作原理，接收机与微弱生命体的距离可以表示为：上式中，是雷达接收机与微弱生命体的平均距离，是由于生命体呼吸和心跳所引起的距离变化。生命体的呼吸和心跳存在着显著的周期性，因此，可以表示为：其中，是心脏跳动信号的峰值振幅、是心脏跳动的频率(1~1.5 Hz)、是心脏跳动信号的初相，是呼吸信号的峰值振幅、是呼吸频率(0.27~0.33 Hz)、是对应的初相。由式(2)，在接收机与生命体无相对运动且之间不存在障碍物时，根据波动方程的行波解，回波信号的表达式为：是雷达发射的载波的频率，若目标区域存在障碍物，如墙，则在式(3)中加上墙体的回波项即可。需要注意的是，若接收机与生命体存在相互运动时，会发生多普勒频移，此时，回波的表达式会发生改变。

MUSIC回波处理算法

一个信号中必然包含有效成分和噪声，而有效成分的子空间与噪声子空间是正交的，因此，可以采用MUSIC算法进行谱估计，达到滤除噪声的目的。MUSIC算法本质上是一种信号频率估计的多重分类算法，利用微弱生命体信号子空间与噪声子空间的正交性构造空间谱函数，通过谱峰搜索，确定心跳和呼吸信号的频率。对于本文中回波信号处理，其主要目标是判断呼吸和心跳所对应的频率的能量是否为真实生命体产生，或者是其他信号产生的谐波。因此，此处不考虑呼吸和心跳信号的初相，假设都为0，回波信号可以表示为：其中，，是方差为σ的白噪声，与信号中的有效成分相互独立。本文中信号的有效成分数量为2(K=2)，即生命体的呼吸和心跳。

定义微弱回波信号向量为：

得到回波信号向量后，就可以计算其自相关矩阵，由于矩阵是一埃尔米特矩阵(共轭转置是其本身)，因此，矩阵可以表达为：其中，是正交的归一化特征向量，是的共轭转置，是对应的特征值。矩阵的个特征值当中，仅有(信号中有效成分的个数)个较大的特征值与有效成分有关，剩余的个特征值与噪声分量有关。噪声子空间的特征向量，可以构成矩阵，如下：其中，矩阵的尺寸为。真实的信号中，噪声分量不能完全满足高斯白噪声这一基本假设，因此，噪声子空间与有效成分子空间并不完全正交。于是，可以构造如下的扫描函数：() 其中，，是信号的频率向量，此处的是归一化圆频率，范围是，与真实的频率之间满足，是信号的采样率。综上，MUSIC算法的回波信号处理流程可以概括为： 1)计算回波信号的自相关矩阵； 2)对进行特征分解，得到个噪声子空间的归一化特征向量； 3)改变圆频率，P将在有效成分的特征频率处产生峰值，根据特征频率值可以确定信号有效成分的频率，根据频率值确定目标区域内是否存在生命体。

仿真验证

雷达发射载波的频率为4 GHz，振幅为1，信号的采样率为20 Hz，共采集50 秒信号。在距离雷达接收机1.2 m处存在一障碍物，5 m处存在一生命体，心跳幅度为2 mm，频率为1.2 Hz，呼吸幅度为4 mm，频率为0.3 Hz。并在回波信号中加入高斯白噪声，使信号的信噪比为10。图1为加入噪声后回波信号的时域波形，单纯观察时域波形是很难分析出有效成分的特征频率的。对信号作傅立叶变换，可以得到其傅立叶域的频谱图(图2)，从中可以清晰地看到两个尖峰(0.3 Hz和1.2 Hz处)，这两个尖峰对应于生命体的呼吸和心跳，然而，仍然有其他尖峰出现，这是由于傅立叶分析很难克服信号中谐波的影响。对于实际信号，采用傅立叶变换分析，可能会误判，甚至遗漏生命体。图3是回波信号的MUSIC谱图，谱图较为干净，不存在其他尖峰，相比于经典的傅立叶分析方法，优势明显。

代码

TestMusic.m：MUSIC算法分析微弱生命体回波信号的示例程序。

%% MUSIC微弱生命体雷达回波信号分析
clear;

% 微弱生命体探测系统参数
c = 2.998e8;          % 光速
f0 = 4e9;             % 发射载波频率
Ac = 1;               % 雷达波振幅
fs = 20;              % 采样频率，Hz
n = 1000;             % 采样总点数
t = (0:n-1)/fs;       % 采样时窗为1000/20 = 50s
fft_num = 1024;       % 傅立叶变换点数

end_f = 3;            % 最大的频率边界

% 生命体参数
fr = 0.3;             % 呼吸频率
fh = 1.2;             % 心跳频率
mr = 4e-3;            % 呼吸振幅
mh = 2e-3;            % 心跳振幅

xr = mr*sin(2*pi*fr*t);  % 呼吸信号
xh = mh*sin(2*pi*fh*t);  % 心跳信号
signalnum = 2;        % 信号中有效成分的数量

x0 = 5;               % 雷达接收机到人体外部距离，单位 m
x1 = 1.2;             % 雷达接收机到障碍物(墙)的距离，单位 m

% 回波信号
signal = Ac*exp(-1i*4*pi*f0*x0/c)*exp(-1i*4*pi*f0/c*(xr+xh))+Ac*exp(-1i*4*pi*f0*x1/c);

SNR = 10;
signal=awgn(signal,SNR,'measured');  % 加入高斯白噪声

signal = signal-mean(signal);  % 去除信号中的零频成分

figure(1);
plot(t,real(signal));
title('含噪声的回波时域波形')
xlabel('时间/s')
ylabel('振幅')


% 对含噪信号作傅立叶变换，并作傅立叶谱分析
y = fft(signal,fft_num);
mag = abs(y/fft_num);
spec = mag(1:fft_num/2+1);
spec(2:end-1) = 2*spec(2:end-1);
freq = fs*(0:(fft_num/2))/fft_num;  % 傅立叶谱的频率向量

figure(2);
plot(freq,spec);
title('含噪声的回波时域波形的频谱')
axis([0 end_f 0 0.6]);
xlabel('频率/Hz')
ylabel('振幅')



% 计算加噪声信号的自相关矩阵，令M=1000;
M = 1000;
R = corrmtx(signal, M-1);
Rxx = R'*R;
[V, D] = eig(Rxx);  % 对相关矩阵作特征值分解


% 计算MUSIC算法中的扫描函数P，以确定有效成分的频率
omega = 2*pi*(0.05:0.005:3)/fs;
l = length(omega);
P = zeros(l,1);    % 扫描函数
for i=1:l
    deno = 0;
    a = exp(-1i*omega(i).*(0:M-1));
for j=1:M-signalnum
        deno = deno+abs((a*V(:,j)))^(2)/D(j,j);
end
    P(i) = 1/deno;
end

% 绘制扫描函数图
figure(3)
plot(omega/(2*pi)*fs,P);
title('含噪声的回波时域波形的MUSIC谱图')
xlabel('频率/Hz')
ylabel('振幅')