文章目录
前言
暑假开始入坑通信感知一体化这个方向,最近在学习一些雷达相关的知识,发现Capon算法和APES算法的资料主要集中在文献和教材,并没有太多的博客对其进行讲解。所以打算把自己最近看到的东西记录一下。
参考教材:MIMO RADAR SIGNAL PROCESSING
在这里打个广告:https://mr-bulijiojio.github.io/
虽然我基本上看不懂,但我感受到了自己和你电本科生雷达天花板的差距。
一、MIMO雷达场景建模
如上图所示,MIMO雷达有
M
t
M_t
Mt个发射天线和
M
r
M_r
Mr个接收天线,同时假设环境中的目标个数为
K
K
K。设
M
r
M_r
Mr根接收天线接收到的信号用
y
(
t
)
=
[
y
1
(
t
)
,
y
2
(
t
)
,
…
…
,
y
M
r
(
t
)
]
T
y(t)=[y_1(t),y_2(t),……,y_{Mr}(t)]^{T}
y(t)=[y1(t),y2(t),……,yMr(t)]T来表示。对于有效可实现的算法,模型可实现的前提条件是信号离散化。因此,假设采样点数为
N
N
N,则接收信号的离散形式为
y
(
n
)
=
[
y
1
(
n
)
,
y
2
(
n
)
,
…
…
,
y
M
r
(
n
)
]
T
y(n)=[y_1(n),y_2(n),……,y_{Mr}(n)]^{T}
y(n)=[y1(n),y2(n),……,yMr(n)]T。
同理,假设发射信号(有
M
t
M_t
Mt根发射天线)的离散形式为
x
(
n
)
=
[
x
1
(
n
)
,
x
2
(
n
)
,
…
…
,
x
M
t
(
n
)
]
T
x(n)=[x_1(n),x_2(n),……,x_{Mt}(n)]^{T}
x(n)=[x1(n),x2(n),……,xMt(n)]T。
这里先介绍多天线阵列信号处理中一个常见的概念:导向矢量,文献中往往翻译成steering vector。如下图所示,考虑均匀分布的多天线阵列,同时假设电磁波到达每个天线的角度相同(远场模型),那么相邻天线之间波传播的路程差即为 d ⋅ s i n θ d \cdot sin\theta d⋅sinθ,根据波的基本性质,可以推出相邻天线之间的相位差 φ k = 2 π d ⋅ s i n θ k λ \varphi_k=\frac{2 \pi d \cdot sin\theta_k}{\lambda} φk=λ2πd⋅sinθk。故当有M根天线的时候,steering vector的形式如下: [ e j ⋅ φ 1 , e j ⋅ φ 2 , … … , e j ⋅ φ M ] T [e^{j \cdot\varphi_1},e^{j \cdot\varphi_2},……,e^{j \cdot\varphi_M}]^{T} [ej⋅φ1,ej⋅φ2,……,ej⋅φM]T(有些参考教材也定义 e − j ⋅ φ 1 , e − j ⋅ φ 2 e^{-j \cdot\varphi_1},e^{-j \cdot\varphi_2} e−j⋅φ1,e−j⋅φ2这种形式,取上述向量的共轭即可)
有了steering vector的定义之后,现在来进行MIMO雷达模型的建模,先作出如下定义(此处统一声明:上标*代表共轭转置,上标c代表共轭,上标T代表转置):
a
(
θ
)
=
[
e
j
2
π
f
0
τ
1
(
θ
)
e
j
2
π
f
0
τ
2
(
θ
)
⋯
e
j
2
π
f
0
τ
M
t
(
θ
)
]
T
\mathbf{a}(\theta)=\left[e^{j 2 \pi f_{0} \tau_{1}(\theta)} \quad e^{j 2 \pi f_{0} \tau_{2}(\theta)} \cdots e^{j 2 \pi f_{0} \tau_{M_{t}(\theta)}}\right]^{T}
a(θ)=[ej2πf0τ1(θ)ej2πf0τ2(θ)⋯ej2πf0τMt(θ)]T
其中
τ
m
\tau_m
τm代表传播延迟,而在目标位置处的信号(多个支路叠加)可以写成:
∑
m
=
1
M
t
e
−
j
2
π
f
0
⋅
τ
m
(
θ
)
⋅
x
m
(
n
)
=
a
∗
(
θ
)
⋅
x
(
n
)
\sum_{m=1}^{M t} e^{-j 2 \pi f 0 \cdot \tau m(\theta)} \cdot x_{m}(n)=a^{*}(\theta) \cdot x(n)
m=1∑Mte−j2πf0⋅τm(θ)⋅xm(n)=a∗(θ)⋅x(n)
我们仔细观察就会不难发现,这个
a
(
θ
)
\mathbf{a}(\theta)
a(θ)的构建思想和steering vector的思想是一样的,也体现了相邻天线的相位差对合并信号的影响。同理,在接收端定义如下的向量:
b
(
θ
)
=
[
e
j
2
π
f
0
t
1
(
θ
)
e
j
2
π
f
0
t
2
(
θ
)
⋯
e
j
2
π
f
0
t
M
r
(
θ
)
]
T
\mathbf{b}(\theta)=\left[e^{j 2 \pi f_{0} t_{1}(\theta)} \quad e^{j 2 \pi f_{0} t_{2}(\theta)} \cdots e^{j 2 \pi f_{0} t_{M_{r}(\theta)}}\right]^{T}
b(θ)=[ej2πf0t1(θ)ej2πf0t2(θ)⋯ej2πf0tMr(θ)]T
我们就可以将接收端的信号写成如下形式:
y
(
n
)
=
∑
k
=
1
K
β
k
⋅
b
c
(
θ
k
)
⋅
a
∗
(
θ
k
)
⋅
x
(
n
)
+
ε
(
n
)
y(n)=\sum_{k=1}^{K} \beta_{k} \cdot b^{c}\left(\theta_{k}\right) \cdot a^{*}\left(\theta_{k}\right) \cdot x(n)+\varepsilon(n)
y(n)=k=1∑Kβk⋅bc(θk)⋅a∗(θk)⋅x(n)+ε(n)
其中
ε
(
n
)
\varepsilon(n)
ε(n)代表噪声矢量,
β
k
\beta_{k}
βk代表对于第
k
k
k个目标的复幅度(体现了幅度和相位信息),其他变量在之前已经声明过,不再重复声明。
二、待解决的问题
事实上,有很多问题可以去求解,例如如何通过优化的方法得出能够正常检测到的目标 K K K的最大值,但我们此处关注的主要是如何去做到达角的估计(AOA estimation)。对于到达角估计算法,最常见的即为MUSIC算法,但是在此处主要讨论Capon算法和APES算法。
对于Capon算法和APES算法,我们主要关注以下两个问题的实现:
1.如何估计出角度
θ
\theta
θ
2.当已知角度
θ
\theta
θ后,如何去估计幅度
β
\beta
β的值
三、Capon算法
1.标准Capon算法
我们首先关注第一个问题,如何去估计角度
θ
\theta
θ。Capon算法的关键步骤就是通过beamforming的方式优化波束设计。此处的beamforming在矢量运算的角度就是在接收信号
Y
Y
Y的左边乘上beamforming vector
w
∗
w^{*}
w∗,关键在于如何设置一个合适的beamforming vector
w
w
w。在Capon算法中,优化问题的建立如下:
min
w
P
=
w
∗
⋅
R
^
y
y
⋅
w
,
s.t
:
ω
∗
⋅
b
c
(
θ
)
=
1
\min _{w} P=w^{*} \cdot \hat{\mathbf{R}}_{y y} \cdot w, \text { s.t }: \omega^{*} \cdot b^c(\theta)=1
wminP=w∗⋅R^yy⋅w, s.t :ω∗⋅bc(θ)=1
如何理解优化问题的意义呢?乍一看确实有些费解,因为传统的beamforming的意义就是让波束方向汇聚,减少损耗,提高功率,但是此处却要求信号的功率达到最低。
事实上,我们可以从类比滤波器的角度去了解,右边的约束条件确保了从一个特定方向来的信号的增益是固定的,那么如果总的功率达到最小,那么就可以证明这个beamforming方式的角度选择性最好。
我们不加证明的给出该优化问题的解:
w
^
capon
=
R
^
y
y
−
1
⋅
b
c
(
θ
)
b
T
(
θ
)
⋅
R
^
y
y
−
1
⋅
b
c
(
θ
)
\widehat{w}_{\text {capon }}=\frac{\hat{R}_{y y} ^{-1} \cdot b^{c}(\theta)}{b^{T}(\theta) \cdot \hat{R}_{y y}^{-1} \cdot b^{c}(\theta)}
w
capon =bT(θ)⋅R^yy−1⋅bc(θ)R^yy−1⋅bc(θ)
而将
w
^
capon
\widehat{w}_{\text {capon }}
w
capon 的值代入之前的等式
P
=
w
∗
⋅
R
^
y
y
⋅
w
P=w^{*} \cdot \hat{\mathbf{R}}_{y y} \cdot w
P=w∗⋅R^yy⋅w,即可得到功率
P
P
P的表达式:
P
=
1
b
T
(
θ
)
⋅
R
^
y
y
−
1
⋅
b
c
(
θ
)
P=\frac{1}{b^{T}(\theta) \cdot \hat{\mathbf{R}}_{y y}^{-1} \cdot b^c(\theta) }
P=bT(θ)⋅R^yy−1⋅bc(θ)1
由上式可知,不同的到达角对应不同的功率,且功率的值和beamforming vector无关,进而可通过和MUSIC类似的方法,寻找功率表达式的极值计算到达角。
我们现在关注第二个问题,已知角度
θ
\theta
θ后,如何估计幅度
β
\beta
β的值。在多目标场景下,不同的
k
k
k对应不同的β,为了简化计算,我们此处考虑单目标场景下
β
\beta
β的估计。则接收信号表达式可以写成:
Y
=
b
c
(
θ
)
β
(
θ
)
a
∗
(
θ
)
X
+
Z
\mathbf{Y}=\mathbf{b}^{c}(\theta) \beta(\theta) \mathbf{a}^{*}(\theta) \mathbf{X}+\mathbf{Z}
Y=bc(θ)β(θ)a∗(θ)X+Z
在经过了beamforming之后,接收信号的表达式可以表示为:
w
∗
⋅
Y
=
w
∗
⋅
b
c
(
θ
)
⋅
β
(
θ
)
⋅
a
∗
(
θ
)
⋅
X
+
ω
∗
⋅
Z
w^{*} \cdot Y=w^{*} \cdot b^{c}(\theta) \cdot \beta(\theta) \cdot a^{*}(\theta) \cdot X+\omega^{*} \cdot Z
w∗⋅Y=w∗⋅bc(θ)⋅β(θ)⋅a∗(θ)⋅X+ω∗⋅Z
将
w
^
capon
\widehat{w}_{\text {capon }}
w
capon 的值代入beamforming后的接收信号表达式,得到如下结果:
b
T
(
θ
)
R
^
y
y
−
1
Y
b
T
(
θ
)
R
^
y
y
−
1
b
c
(
θ
)
=
β
(
θ
)
a
∗
(
θ
)
X
+
b
T
(
θ
)
R
^
y
y
−
1
Z
b
T
(
θ
)
R
^
y
y
−
1
b
c
(
θ
)
\frac{\mathbf{b}^{T}(\theta) \hat{\mathbf{R}}_{y y}^{-1} \mathbf{Y}}{\mathbf{b}^{T}(\theta) \hat{\mathbf{R}}_{y y}^{-1} \mathbf{b}^{c}(\theta)}=\beta(\theta) \mathbf{a}^{*}(\theta) \mathbf{X}+\frac{\mathbf{b}^{T}(\theta) \hat{\mathbf{R}}_{y y}^{-1} \mathbf{Z}}{\mathbf{b}^{T}(\theta) \hat{\mathbf{R}}_{y y}^{-1} \mathbf{b}^{c}(\theta)}
bT(θ)R^yy−1bc(θ)bT(θ)R^yy−1Y=β(θ)a∗(θ)X+bT(θ)R^yy−1bc(θ)bT(θ)R^yy−1Z
事实上,由于
β
\beta
β是一个标量,最简单的估计方法就是忽略噪声的影响,并且用LS的思想(least square)直接将
β
\beta
β解出。基本的矢量运算的过程就不详细推导了,给出最终得到的结果(分为不考虑beamforming vector的估计值和考虑beamforming vector的估计值,个人理解为beamforming vector的引入进一步优化了对
β
\beta
β的估计)。
不考虑beamforming vector时,得到的估计结果如下:
β
^
L
S
(
θ
)
=
b
T
(
θ
)
Y
X
∗
a
(
θ
)
N
∥
b
(
θ
)
∥
2
[
a
∗
(
θ
)
R
^
x
x
a
(
θ
)
]
\hat{\beta}_{\mathrm{LS}}(\theta)=\frac{\mathbf{b}^{T}(\theta) \mathbf{Y} \mathbf{X}^{*} \mathbf{a}(\theta)}{N\|\mathbf{b}(\theta)\|^{2}\left[\mathbf{a}^{*}(\theta) \hat{\mathbf{R}}_{x x} \mathbf{a}(\theta)\right]}
β^LS(θ)=N∥b(θ)∥2[a∗(θ)R^xxa(θ)]bT(θ)YX∗a(θ)
其中:
R
^
x
x
=
1
N
X
X
∗
,
R
^
y
y
=
1
N
Y
Y
∗
\hat{\mathbf{R}}_{x x}=\frac{1}{N} \mathbf{X} \mathbf{X}^{*},\hat{\mathbf{R}}_{y y}=\frac{1}{N} \mathbf{Y} \mathbf{Y}^{*}
R^xx=N1XX∗,R^yy=N1YY∗
考虑beamforming vector时,得到的估计结果如下:
β
^
Capon
(
θ
)
=
b
T
(
θ
)
R
^
y
y
−
1
Y
X
∗
a
(
θ
)
N
[
b
T
(
θ
)
R
^
y
y
−
1
b
c
(
θ
)
]
[
a
∗
(
θ
)
R
^
x
x
a
(
θ
)
]
\hat{\beta}_{\text {Capon }}(\theta)=\frac{\mathbf{b}^{T}(\theta) \hat{\mathbf{R}}_{y y}^{-1} \mathbf{Y} \mathbf{X}^{*} \mathbf{a}(\theta)}{N\left[\mathbf{b}^{T}(\theta) \hat{\mathbf{R}}_{y y}^{-1} \mathbf{b}^{c}(\theta)\right]\left[\mathbf{a}^{*}(\theta) \hat{\mathbf{R}}_{x x} \mathbf{a}(\theta)\right]}
β^Capon (θ)=N[bT(θ)R^yy−1bc(θ)][a∗(θ)R^xxa(θ)]bT(θ)R^yy−1YX∗a(θ)
2.Capon算法的不足
Capon算法的不足之处主要可以归结为两点:
1.Capon的数值估计的准确性十分依赖于steering vector估计的准确性,如果steering vector的估计值和真实值有一定的差异,Capon算法的性能会受到较大的影响。但实际应用中,往往只能将steering vector确定在一个范围之间,因此,其改进算法即为Robust Capon beamforming。限于篇幅,此处不讲解去具体原理,大家有兴趣可以去Google scholar上搜一下引用量较高的一篇经典文章。
2.即使steering vector的估计和真实值完全一致,但Capon算法在估计
β
\beta
β的时候也忽略了噪声项的影响(其实是一种zero forcing的思想)。在SNR条件较好的情况下,这样不会带来太高的误差,但当SNR条件较差时,噪声项直接忽略势必会对估计结果的准确性带来影响。因此,在低信噪比条件下,Capon方法的性能会比MUSIC方法差。
四、APES算法(amplitude and phase estimation)
由于干扰项的存在,我们很自然的想到,如果要提高估计精度,我们可以直接将优化问题构造成使得beamforming之后的噪声项功率最小。但是,这样是无法直接求解的,因为我们并不知道噪声信号的方差,但是,接收端的 Y Y Y信号是我们已知的,我们可以通过构造另一个含义相同的优化表达式来表达这个思想:
min
w
,
β
∥
w
∗
⋅
Y
−
β
(
θ
)
⋅
a
∗
(
θ
)
⋅
X
∥
2
,
s.t
:
w
∗
⋅
b
c
(
θ
)
=
1
\min _{w, \beta}\left\|w^{*} \cdot Y-\beta(\theta) \cdot a^{*}(\theta) \cdot X\right\|^{2}, \text { s.t }: w^{*} \cdot b^{c}(\theta)=1
w,βmin∥w∗⋅Y−β(θ)⋅a∗(θ)⋅X∥2, s.t :w∗⋅bc(θ)=1
这是一个关于
w
w
w和
β
\beta
β的双变量优化问题,再对比经过beamforming后接收信号的表达式:
w
∗
⋅
Y
=
w
∗
⋅
b
c
(
θ
)
⋅
β
(
θ
)
⋅
a
∗
(
θ
)
⋅
X
+
ω
∗
⋅
Z
w^{*} \cdot Y=w^{*} \cdot b^{c}(\theta) \cdot \beta(\theta) \cdot a^{*}(\theta) \cdot X+\omega^{*} \cdot Z
w∗⋅Y=w∗⋅bc(θ)⋅β(θ)⋅a∗(θ)⋅X+ω∗⋅Z
很容易发现其本质就是使得干扰项的影响降到最低。
对于该优化问题的求解,可以先求出 w w w一定时 β \beta β的最优解,将双变量优化问题转化为单变量优化问题,再解出 w w w的最优解并代入 β \beta β,得到 β ^ APES ( θ ) \hat{\beta}_{\text {APES }}(\theta) β^APES (θ)。
β
^
APES
(
θ
)
\hat{\beta}_{\text {APES }}(\theta)
β^APES (θ)用
w
w
w表示的结果如下:
β
^
A
P
E
S
(
θ
)
=
w
∗
Y
X
∗
a
(
θ
)
N
a
∗
(
θ
)
R
^
x
x
a
(
θ
)
\hat{\beta}_{\mathrm{APES}}(\theta)=\frac{\mathbf{w}^{*} \mathbf{Y} \mathbf{X}^{*} \mathbf{a}(\theta)}{N \mathbf{a}^{*}(\theta) \hat{\mathbf{R}}_{x x} \mathbf{a}(\theta)}
β^APES(θ)=Na∗(θ)R^xxa(θ)w∗YX∗a(θ)
则可转化为如下的单变量优化问题:
min
w
w
∗
Q
^
w
\min _{\mathbf{w}} \mathbf{w}^{*} \hat{\mathbf{Q}} \mathbf{w} \quad
minww∗Q^w subject to
w
∗
b
c
(
θ
)
=
1
\quad \mathbf{w}^{*} \mathbf{b}^{c}(\theta)=1
w∗bc(θ)=1
Q
^
=
R
^
y
y
−
Y
X
∗
a
(
θ
)
a
∗
(
θ
)
X
Y
∗
N
2
a
∗
(
θ
)
R
^
x
x
a
(
θ
)
\hat{\mathbf{Q}}=\hat{\mathbf{R}}_{y y}-\frac{\mathbf{Y} \mathbf{X}^{*} \mathbf{a}(\theta) \mathbf{a}^{*}(\theta) \mathbf{X} \mathbf{Y}^{*}}{N^{2} \mathbf{a}^{*}(\theta) \hat{\mathbf{R}}_{x x} \mathbf{a}(\theta)}
Q^=R^yy−N2a∗(θ)R^xxa(θ)YX∗a(θ)a∗(θ)XY∗
有趣的是,这个时候优化问题的结构就和之前的Capon算法一模一样了,因此也对应一样的解法。但是需要说明的是,虽然这两个优化问题可以变换为一样的结构,这两种算法本身的思想和效果仍存在较大的差异。
可以解出
w
w
w的最优解如下:
w
^
A
P
E
S
=
Q
^
−
1
b
c
(
θ
)
b
T
(
θ
)
Q
^
−
1
b
c
(
θ
)
\hat{\mathbf{w}}_{\mathrm{APES}}=\frac{\hat{\mathbf{Q}}^{-1} \mathbf{b}^{c}(\theta)}{\mathbf{b}^{T}(\theta) \hat{\mathbf{Q}}^{-1} \mathbf{b}^{c}(\theta)}
w^APES=bT(θ)Q^−1bc(θ)Q^−1bc(θ)
代入到
β
\beta
β的表达式,得到
β
^
A
P
E
S
(
θ
)
\hat{\beta}_{\mathrm{APES}}(\theta)
β^APES(θ)的值为:
β
^
A
P
E
S
(
θ
)
=
b
T
(
θ
)
Q
^
−
1
Y
X
∗
a
(
θ
)
N
[
b
T
(
θ
)
Q
^
−
1
b
c
(
θ
)
]
[
a
∗
(
θ
)
R
^
x
x
a
(
θ
)
]
\hat{\beta}_{\mathrm{APES}}(\theta)=\frac{\mathbf{b}^{T}(\theta) \hat{\mathbf{Q}}^{-1} \mathbf{Y} \mathbf{X}^{*} \mathbf{a}(\theta)}{N\left[\mathbf{b}^{T}(\theta) \hat{\mathbf{Q}}^{-1} \mathbf{b}^{c}(\theta)\right]\left[\mathbf{a}^{*}(\theta) \hat{\mathbf{R}}_{x x} \mathbf{a}(\theta)\right]}
β^APES(θ)=N[bT(θ)Q^−1bc(θ)][a∗(θ)R^xxa(θ)]bT(θ)Q^−1YX∗a(θ)
总结
以上就是MIMO雷达场景建模的过程,Capon算法和APES算法的基本思想和解法了。当然这些只是较为基础的知识,在更具体更复杂的科研场景中,算法的形式也会或多或少的发生改变。鉴于博主还是学生,且能力有限、时间仓促,如有谬误,欢迎各位大佬指正交流。