Cadence Velocity Diagrams 在毫米波雷达系统中的应用
毫米波雷达凭借其在全天候条件下的探测能力、高精度的距离和速度分辨率以及对生物体非侵入性特性,正在智能监测和人体运动分析领域获得广泛应用。在这一背景下,Cadence Velocity Diagrams(步频-速度图)作为一种分析工具,为毫米波雷达在人体运动特征提取方面提供了强大的理论支持和实际应用价值。
毫米波雷达基本原理与步频检测
毫米波雷达通常工作在24GHz、60GHz或77GHz频段,通过发射电磁波并接收目标反射回波来探测环境。在人体运动分析中,毫米波雷达能够捕捉到人体各部位的微小运动,特别是步行和跑步过程中的周期性模式。
雷达信号模型与多普勒效应
当雷达信号照射到运动目标时,反射信号会因多普勒效应而产生频移。对于运动人体,其不同部位以不同速度运动,产生一系列多普勒频移。设发射信号为:
s
t
(
t
)
=
A
t
cos
(
2
π
f
c
t
+
ϕ
t
)
s_t(t) = A_t \cos(2\pi f_c t + \phi_t)
st(t)=Atcos(2πfct+ϕt)
其中,
A
t
A_t
At是发射信号幅度,
f
c
f_c
fc是载波频率,
ϕ
t
\phi_t
ϕt是初始相位。
当信号被运动速度为
v
(
t
)
v(t)
v(t)的目标反射后,接收信号可表示为:
s
r
(
t
)
=
A
r
cos
(
2
π
f
c
(
t
−
τ
(
t
)
)
+
ϕ
r
)
s_r(t) = A_r \cos(2\pi f_c (t-\tau(t)) + \phi_r)
sr(t)=Arcos(2πfc(t−τ(t))+ϕr)
其中,
τ
(
t
)
=
2
R
(
t
)
c
\tau(t) = \frac{2R(t)}{c}
τ(t)=c2R(t)是信号往返延迟,
R
(
t
)
R(t)
R(t)是目标到雷达的瞬时距离,
c
c
c是光速。考虑到目标运动,
R
(
t
)
=
R
0
+
∫
0
t
v
(
t
′
)
d
t
′
R(t) = R_0 + \int_0^t v(t')dt'
R(t)=R0+∫0tv(t′)dt′,其中
R
0
R_0
R0是初始距离。
混频后的信号包含多普勒频移:
f
d
(
t
)
=
2
v
(
t
)
λ
=
2
v
(
t
)
f
c
c
f_d(t) = \frac{2v(t)}{\lambda} = \frac{2v(t)f_c}{c}
fd(t)=λ2v(t)=c2v(t)fc
其中,
λ
\lambda
λ是雷达波长。
微多普勒分析与特征提取
人体运动产生的微多普勒特征是识别步频的关键。通过短时傅里叶变换(STFT),雷达数据可以转换为时频域表示:
STFT
(
t
,
f
)
=
∫
−
∞
∞
s
(
t
′
)
w
(
t
′
−
t
)
e
−
j
2
π
f
t
′
d
t
′
\text{STFT}(t,f) = \int_{-\infty}^{\infty} s(t') w(t'-t) e^{-j2\pi ft'} dt'
STFT(t,f)=∫−∞∞s(t′)w(t′−t)e−j2πft′dt′
其中,
w
(
t
)
w(t)
w(t)是窗函数,常用的有汉宁窗或汉明窗。结果形成微多普勒频谱图,显示不同身体部位的多普勒特征随时间的变化。
对于步行分析,腿部运动产生的微多普勒特征最为明显。一个完整的步态周期包含摆动相和支撑相,这在微多普勒频谱上表现为具有特定周期性的模式。步频可以通过以下方法提取:
- 频谱峰值追踪: f peak ( t ) = arg max f ∣ STFT ( t , f ) ∣ f_{\text{peak}}(t) = \arg\max_f |\text{STFT}(t,f)| fpeak(t)=argfmax∣STFT(t,f)∣
- 周期性分析:应用自相关函数或周期图分析: R x x ( τ ) = ∫ − ∞ ∞ f peak ( t ) f peak ( t + τ ) d t R_{xx}(\tau) = \int_{-\infty}^{\infty} f_{\text{peak}}(t) f_{\text{peak}}(t+\tau) dt Rxx(τ)=∫−∞∞fpeak(t)fpeak(t+τ)dt
- 步频估计: C = 1 T stride C = \frac{1}{T_{\text{stride}}} C=Tstride1 其中, T stride T_{\text{stride}} Tstride是通过自相关峰值间距确定的步态周期。
Cadence-Velocity 关系的数学模型
基于多普勒频移的速度估计
毫米波雷达可以直接通过多普勒频移估计目标速度:
v
(
t
)
=
f
d
(
t
)
⋅
c
2
f
c
v(t) = \frac{f_d(t) \cdot c}{2f_c}
v(t)=2fcfd(t)⋅c
考虑到人体不同部位的多普勒贡献,整体速度可以表示为加权和:
v
body
(
t
)
=
∑
i
=
1
N
σ
i
v
i
(
t
)
∑
i
=
1
N
σ
i
v_{\text{body}}(t) = \frac{\sum_{i=1}^{N} \sigma_i v_i(t)}{\sum_{i=1}^{N} \sigma_i}
vbody(t)=∑i=1Nσi∑i=1Nσivi(t)
其中,
σ
i
\sigma_i
σi是第
i
i
i个散射点的雷达截面积(RCS),
v
i
(
t
)
v_i(t)
vi(t)是对应的速度,
N
N
N是散射点总数。
非线性步频-速度关系模型
在毫米波雷达应用中,步频与速度的关系通常需要考虑非线性因素。扩展二次模型可表示为:
C
(
v
)
=
a
0
+
a
1
v
+
a
2
v
2
+
a
3
v
3
+
ϵ
(
v
)
C(v) = a_0 + a_1v + a_2v^2 + a_3v^3 + \epsilon(v)
C(v)=a0+a1v+a2v2+a3v3+ϵ(v)
其中,
ϵ
(
v
)
\epsilon(v)
ϵ(v)是与速度相关的误差项,可以通过自回归移动平均(ARMA)模型描述:
ϵ
(
v
)
=
∑
i
=
1
p
ϕ
i
ϵ
(
v
−
Δ
v
i
)
+
∑
j
=
1
q
θ
j
η
(
v
−
Δ
v
j
)
+
η
(
v
)
\epsilon(v) = \sum_{i=1}^{p} \phi_i \epsilon(v-\Delta v_i) + \sum_{j=1}^{q} \theta_j \eta(v-\Delta v_j) + \eta(v)
ϵ(v)=i=1∑pϕiϵ(v−Δvi)+j=1∑qθjη(v−Δvj)+η(v)
ϕ
i
\phi_i
ϕi和
θ
j
\theta_j
θj分别是自回归和移动平均参数,
η
(
v
)
\eta(v)
η(v)是白噪声过程。
状态空间表示与卡尔曼滤波
为了实时跟踪和预测步频与速度的关系,可以使用状态空间模型和卡尔曼滤波:
定义状态向量
x
(
t
)
=
[
v
(
t
)
,
C
(
t
)
,
v
˙
(
t
)
,
C
˙
(
t
)
]
T
\mathbf{x}(t) = [v(t), C(t), \dot{v}(t), \dot{C}(t)]^T
x(t)=[v(t),C(t),v˙(t),C˙(t)]T,状态转移方程为:
x
(
t
+
Δ
t
)
=
F
x
(
t
)
+
w
(
t
)
\mathbf{x}(t+\Delta t) = \mathbf{F} \mathbf{x}(t) + \mathbf{w}(t)
x(t+Δt)=Fx(t)+w(t)
其中,
F
\mathbf{F}
F是状态转移矩阵:
F
=
[
1
0
Δ
t
0
0
1
0
Δ
t
0
0
1
0
0
0
0
1
]
\mathbf{F} = \begin{bmatrix} 1 & 0 & \Delta t & 0 \ 0 & 1 & 0 & \Delta t \ 0 & 0 & 1 & 0 \ 0 & 0 & 0 & 1 \end{bmatrix}
F=[10Δt0 010Δt 0010 0001]
w
(
t
)
\mathbf{w}(t)
w(t)是过程噪声,假设服从协方差矩阵为
Q
\mathbf{Q}
Q的零均值高斯分布。
观测方程为:
z
(
t
)
=
H
x
(
t
)
+
v
(
t
)
\mathbf{z}(t) = \mathbf{H} \mathbf{x}(t) + \mathbf{v}(t)
z(t)=Hx(t)+v(t)
其中,
H
\mathbf{H}
H是观测矩阵,
v
(
t
)
\mathbf{v}(t)
v(t)是观测噪声,服从协方差矩阵为
R
\mathbf{R}
R的零均值高斯分布。
卡尔曼滤波的预测步骤:
x
^
−
(
t
)
=
F
x
^
(
t
−
Δ
t
)
\hat{\mathbf{x}}^-(t) = \mathbf{F} \hat{\mathbf{x}}(t-\Delta t)
x^−(t)=Fx^(t−Δt)
P
−
(
t
)
=
F
P
(
t
−
Δ
t
)
F
T
+
Q
\mathbf{P}^-(t) = \mathbf{F} \mathbf{P}(t-\Delta t) \mathbf{F}^T + \mathbf{Q}
P−(t)=FP(t−Δt)FT+Q
更新步骤:
K
(
t
)
=
P
−
(
t
)
H
T
[
H
P
−
(
t
)
H
T
+
R
]
−
1
\mathbf{K}(t) = \mathbf{P}^-(t) \mathbf{H}^T [\mathbf{H} \mathbf{P}^-(t) \mathbf{H}^T + \mathbf{R}]^{-1}
K(t)=P−(t)HT[HP−(t)HT+R]−1
x
^
(
t
)
=
x
^
−
(
t
)
+
K
(
t
)
[
z
(
t
)
−
H
x
^
−
(
t
)
]
\hat{\mathbf{x}}(t) = \hat{\mathbf{x}}^-(t) + \mathbf{K}(t) [\mathbf{z}(t) - \mathbf{H} \hat{\mathbf{x}}^-(t)]
x^(t)=x^−(t)+K(t)[z(t)−Hx^−(t)]
P
(
t
)
=
[
I
−
K
(
t
)
H
]
P
−
(
t
)
\mathbf{P}(t) = [I - \mathbf{K}(t) \mathbf{H}] \mathbf{P}^-(t)
P(t)=[I−K(t)H]P−(t)
这提供了对步频-速度关系的平滑估计和预测。
毫米波雷达中的微多普勒特征与步频分析
微多普勒特征分解
人体运动在微多普勒频谱上形成复杂模式,可以通过奇异值分解(SVD)进行特征提取:
S
=
U
Σ
V
∗
\mathbf{S} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^*
S=UΣV∗
其中,
S
\mathbf{S}
S是微多普勒频谱矩阵,
U
\mathbf{U}
U和
V
\mathbf{V}
V包含左右奇异向量,
Σ
\mathbf{\Sigma}
Σ是奇异值对角矩阵。主要特征可以通过选择前
k
k
k个最大奇异值对应的成分来表示。
步频特征的小波变换分析
小波变换提供了更精细的时频分析能力,特别适合分析步频的非平稳特性:
CWT
(
a
,
b
)
=
1
a
∫
−
∞
∞
s
(
t
)
ψ
∗
(
t
−
b
a
)
d
t
\text{CWT}(a,b) = \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} s(t) \psi^*\left(\frac{t-b}{a}\right) dt
CWT(a,b)=a1∫−∞∞s(t)ψ∗(at−b)dt
其中,
ψ
(
t
)
\psi(t)
ψ(t)是小波函数,
a
a
a是尺度参数,
b
b
b是平移参数。对于步频检测,Morlet小波是常用选择:
ψ
(
t
)
=
1
π
f
b
e
2
π
i
f
c
t
e
−
t
2
/
f
b
\psi(t) = \frac{1}{\sqrt{\pi f_b}} e^{2\pi if_c t} e^{-t^2/f_b}
ψ(t)=πfb1e2πifcte−t2/fb
其中,
f
c
f_c
fc是中心频率,
f
b
f_b
fb控制时频分辨率。
步频可以通过分析小波系数的局部周期性来提取:
C
(
t
)
=
1
a
^
(
t
)
C(t) = \frac{1}{\hat{a}(t)}
C(t)=a^(t)1
其中,
a
^
(
t
)
\hat{a}(t)
a^(t)是在时间
t
t
t处小波能量谱中占主导的尺度参数。
深度学习增强的特征提取
最近的研究将深度学习应用于微多普勒数据处理。一种方法是使用卷积神经网络(CNN)自动提取步频特征。网络架构通常包括多个卷积层、池化层和全连接层,训练目标是最小化预测步频与真实值之间的损失:
L
=
1
N
∑
i
=
1
N
(
C
pred
,
i
−
C
true
,
i
)
2
+
λ
∣
W
∣
2
\mathcal{L} = \frac{1}{N} \sum_{i=1}^{N} (C_{\text{pred},i} - C_{\text{true},i})^2 + \lambda |\mathbf{W}|^2
L=N1i=1∑N(Cpred,i−Ctrue,i)2+λ∣W∣2
其中,
λ
\lambda
λ是正则化参数,
W
\mathbf{W}
W表示网络权重。
毫米波雷达中步频-速度关系的实际应用
行人意图识别
基于步频-速度特征,可以分析行人的运动意图。例如,加速、减速或转向等意图可以通过步频-速度曲线的变化来预测。这对自动驾驶系统和智能交通管理至关重要。定义意图概率模型:
P
(
I
∣
C
,
v
,
θ
)
=
P
(
C
,
v
,
θ
∣
I
)
P
(
I
)
P
(
C
,
v
,
θ
)
P(I|C,v,\theta) = \frac{P(C,v,\theta|I)P(I)}{P(C,v,\theta)}
P(I∣C,v,θ)=P(C,v,θ)P(C,v,θ∣I)P(I)
其中,
I
I
I表示意图类别,
θ
\theta
θ是行人朝向角,
P
(
C
,
v
,
θ
∣
I
)
P(C,v,\theta|I)
P(C,v,θ∣I)可以通过训练高斯混合模型(GMM)来估计:
P
(
C
,
v
,
θ
∣
I
)
=
∑
k
=
1
K
π
k
N
(
[
C
,
v
,
θ
]
∣
μ
k
,
Σ
k
)
P(C,v,\theta|I) = \sum_{k=1}^{K} \pi_k \mathcal{N}([C,v,\theta]|\mu_k,\Sigma_k)
P(C,v,θ∣I)=k=1∑KπkN([C,v,θ]∣μk,Σk)
π
k
\pi_k
πk、
μ
k
\mu_k
μk和
Σ
k
\Sigma_k
Σk分别是第
k
k
k个混合成分的权重、均值和协方差。
跌倒检测与健康监测
毫米波雷达结合步频-速度分析可以有效检测异常步态和跌倒事件。正常步行的步频-速度关系相对稳定,而跌倒前通常表现为步频-速度关系的突变。跌倒风险评估函数可以定义为:
R
fall
(
t
)
=
α
∣
d
C
(
t
)
d
t
∣
+
β
∣
d
v
(
t
)
d
t
∣
+
γ
∣
d
2
v
(
t
)
d
t
2
∣
R_{\text{fall}}(t) = \alpha \left|\frac{dC(t)}{dt}\right| + \beta \left|\frac{dv(t)}{dt}\right| + \gamma \left|\frac{d^2v(t)}{dt^2}\right|
Rfall(t)=α
dtdC(t)
+β
dtdv(t)
+γ
dt2d2v(t)
其中,
α
\alpha
α、
β
\beta
β和
γ
\gamma
γ是权重系数,可通过机器学习方法优化。当
R
fall
(
t
)
R_{\text{fall}}(t)
Rfall(t)超过预设阈值时,系统触发警报。
生物识别与身份验证
个体的步频-速度特征具有生物识别潜力。每个人的步态模式具有独特性,可以通过步频-速度关系的参数化表示进行量化:
f
gait
=
[
a
0
,
a
1
,
a
2
,
.
.
.
,
σ
C
2
,
ρ
C
v
,
.
.
.
]
T
\mathbf{f}{\text{gait}} = [a_0, a_1, a_2, ..., \sigma_C^2, \rho{Cv}, ...]^T
fgait=[a0,a1,a2,...,σC2,ρCv,...]T
其中包含步频-速度模型参数和统计特征(如步频方差
σ
C
2
\sigma_C^2
σC2、步频与速度的相关系数
ρ
C
v
\rho_{Cv}
ρCv等)。身份验证可以通过计算当前观测特征与参考模板之间的马氏距离实现:
d
M
(
f
,
f
ref
)
=
(
f
−
f
ref
)
T
Σ
−
1
(
f
−
f
ref
)
d_M(\mathbf{f}, \mathbf{f}{\text{ref}}) = \sqrt{(\mathbf{f} - \mathbf{f}{\text{ref}})^T \mathbf{\Sigma}^{-1} (\mathbf{f} - \mathbf{f}_{\text{ref}})}
dM(f,fref)=(f−fref)TΣ−1(f−fref)
其中,
Σ
\mathbf{\Sigma}
Σ是特征空间的协方差矩阵。
多普勒频谱分析
时频分布与Cohen类分析
Cohen类时频分布提供了比STFT更灵活的时频分析工具:
P
C
(
t
,
f
)
=
∭
−
∞
∞
e
j
2
π
[
θ
(
t
′
−
t
)
−
ν
τ
−
f
(
t
′
′
−
t
)
]
ϕ
(
θ
,
τ
)
s
∗
(
t
′
−
τ
/
2
)
s
(
t
′
+
τ
/
2
)
d
θ
d
τ
d
t
′
P_C(t,f) = \iiint_{-\infty}^{\infty} e^{j2\pi[\theta(t'-t) - \nu\tau - f(t''-t)]} \phi(\theta,\tau) s^*(t'-\tau/2) s(t'+\tau/2) d\theta d\tau dt'
PC(t,f)=∭−∞∞ej2π[θ(t′−t)−ντ−f(t′′−t)]ϕ(θ,τ)s∗(t′−τ/2)s(t′+τ/2)dθdτdt′
其中,
ϕ
(
θ
,
τ
)
\phi(\theta,\tau)
ϕ(θ,τ)是核函数,不同选择产生不同的分布。例如,Wigner-Ville分布(WVD)的核函数为
ϕ
(
θ
,
τ
)
=
1
\phi(\theta,\tau) = 1
ϕ(θ,τ)=1:
WVD
(
t
,
f
)
=
∫
−
∞
∞
s
(
t
+
τ
2
)
s
∗
(
t
−
τ
2
)
e
−
j
2
π
f
τ
d
τ
\text{WVD}(t,f) = \int_{-\infty}^{\infty} s\left(t+\frac{\tau}{2}\right) s^*\left(t-\frac{\tau}{2}\right) e^{-j2\pi f\tau} d\tau
WVD(t,f)=∫−∞∞s(t+2τ)s∗(t−2τ)e−j2πfτdτ
平滑伪Wigner-Ville分布(SPWVD)引入时间和频率平滑:
SPWVD
(
t
,
f
)
=
∫
−
∞
∞
∫
−
∞
∞
g
(
s
)
h
(
τ
)
s
(
t
−
s
+
τ
2
)
s
∗
(
t
−
s
−
τ
2
)
e
−
j
2
π
f
τ
d
s
d
τ
\text{SPWVD}(t,f) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(s)h(\tau) s\left(t-s+\frac{\tau}{2}\right) s^*\left(t-s-\frac{\tau}{2}\right) e^{-j2\pi f\tau} ds d\tau
SPWVD(t,f)=∫−∞∞∫−∞∞g(s)h(τ)s(t−s+2τ)s∗(t−s−2τ)e−j2πfτdsdτ
其中,
g
(
s
)
g(s)
g(s)和
h
(
τ
)
h(\tau)
h(τ)分别是时间和频率平滑窗。
稀疏表示与压缩感知
微多普勒信号通常具有稀疏性,可以利用压缩感知技术进行高效处理:
y
=
Φ
x
+
n
\mathbf{y} = \mathbf{\Phi} \mathbf{x} + \mathbf{n}
y=Φx+n
其中,
y
\mathbf{y}
y是测量向量,
Φ
\mathbf{\Phi}
Φ是测量矩阵,
x
\mathbf{x}
x是稀疏信号,
n
\mathbf{n}
n是噪声。重构问题可以表述为:
min
x
∣
x
∣
1
subject to
∣
y
−
Φ
x
∣
2
≤
ϵ
\min_{\mathbf{x}} |\mathbf{x}|_1 \quad \text{subject to} \quad |\mathbf{y} - \mathbf{\Phi} \mathbf{x}|_2 \leq \epsilon
xmin∣x∣1subject to∣y−Φx∣2≤ϵ
解决方案包括正交匹配追踪(OMP)和基追踪(BP)等算法。
联合时频分析与经验模态分解
经验模态分解(EMD)将信号分解为多个固有模态函数(IMF):
s
(
t
)
=
∑
i
=
1
N
c
i
(
t
)
+
r
N
(
t
)
s(t) = \sum_{i=1}^{N} c_i(t) + r_N(t)
s(t)=i=1∑Nci(t)+rN(t)
其中,
c
i
(
t
)
c_i(t)
ci(t)是IMF,
r
N
(
t
)
r_N(t)
rN(t)是残差。对每个IMF应用Hilbert变换得到分析信号:
z
i
(
t
)
=
c
i
(
t
)
+
j
H
c
i
(
t
)
=
a
i
(
t
)
e
j
ϕ
i
(
t
)
z_i(t) = c_i(t) + j \mathcal{H}{c_i(t)} = a_i(t) e^{j\phi_i(t)}
zi(t)=ci(t)+jHci(t)=ai(t)ejϕi(t)
瞬时频率定义为:
f
i
(
t
)
=
1
2
π
d
ϕ
i
(
t
)
d
t
f_i(t) = \frac{1}{2\pi} \frac{d\phi_i(t)}{dt}
fi(t)=2π1dtdϕi(t)
Hilbert-Huang变换(HHT)提供了Hilbert频谱:
H
(
ω
,
t
)
=
∑
i
=
1
N
a
i
(
t
)
δ
(
ω
−
ω
i
(
t
)
)
H(\omega, t) = \sum_{i=1}^{N} a_i(t) \delta(\omega - \omega_i(t))
H(ω,t)=i=1∑Nai(t)δ(ω−ωi(t))
这种方法特别适合分析非平稳的步频变化。