Singer模型与CT模型状态转移矩阵的求解
前言
回想起来,第一次接触Singer模型与CT模型时的状态转移矩阵时,对求解过程一知半解。现在,将推导过程记录下来,温故知新。若有纰漏,欢迎讨论。
状态方程
当我们对一个线性时不变系统的输入输出与状态变量感兴趣时,状态空间方程是一种常用的手段。状态空间方程由状态方程与输出方程构成,可以表述为
{
x
˙
=
A
x
+
B
u
状态方程
y
=
C
x
+
D
u
输出方程
\left\{\begin{matrix} \dot{ \mathbf{x}} = \mathbf{Ax}+\mathbf{Bu}& 状态方程\\ \mathbf{y} = \mathbf{Cx} + \mathbf{Du} & 输出方程\\ \end{matrix}\right.
{x˙=Ax+Buy=Cx+Du状态方程输出方程
其中,状态方程是一个一阶微分方程。
x
\mathbf{x}
x与
x
˙
\dot{\mathbf{x}}
x˙表示状态向量及其一阶导数。
u
\mathbf{u}
u为系统的外部输入。
A
\mathbf{A}
A表示系统矩阵。
B
\mathbf{B}
B表示输入矩阵。
y
\mathbf{y}
y为系统输出。
C
\mathbf{C}
C为输出矩阵。
D
\mathbf{D}
D为直接输出矩阵。
我们只对状态方程感兴趣。现在为了简化问题,我们并不考虑外部输入对系统的影响,因此,状态方程简化为
x
˙
=
A
x
\dot{\mathbf{x}} = \mathbf{A}\mathbf{x}
x˙=Ax
使用拉普拉斯变换对上述一阶微分方程进行求解。这里补充几个常用的拉普拉斯变换对
L ( x ( t ) ) = X ( s ) L ( x ˙ ( t ) ) = s X ( s ) − x ( 0 ) L ( e a t ) = 1 s − a L ( 1 ) = 1 s L ( t ) = 1 s 2 L ( c o s ω t ) = s s 2 + ω 2 L ( s i n ω t ) = ω s 2 + ω 2 L ( e − α t s i n ( ω t ) ) = ω ( s + α ) 2 + ω 2 \begin{aligned} \mathcal{L}(x(t)) &= X(s) \\ \mathcal{L}(\dot{x}(t)) &= sX(s) - x(0) \\ \mathcal{L}(e^{at}) &= \frac{1}{s-a} \\ \mathcal{L}(1) &= \frac{1}{s} \\ \mathcal{L}(t) &= \frac{1}{s^2} \\ \mathcal{L}(cos\omega t) &= \frac{s}{s^2 + \omega ^2} \\ \mathcal{L}(sin\omega t) &= \frac{\omega}{s^2 + \omega ^2} \\ \mathcal{L}(e^{-\alpha t}sin(\omega t)) &= \frac{\omega}{(s + \alpha)^2 + \omega ^2} \end{aligned} L(x(t))L(x˙(t))L(eat)L(1)L(t)L(cosωt)L(sinωt)L(e−αtsin(ωt))=X(s)=sX(s)−x(0)=s−a1=s1=s21=s2+ω2s=s2+ω2ω=(s+α)2+ω2ω
现在对简化的状态方程的等式两端进行拉普拉斯变换并进行变换,可得状态方程的解
L
(
x
˙
(
t
)
)
=
L
(
A
x
(
t
)
)
s
X
(
s
)
−
x
(
0
)
=
A
X
(
s
)
(
s
−
A
)
X
(
s
)
=
x
(
0
)
X
(
s
)
=
x
(
0
)
s
−
A
L
−
1
(
X
(
s
)
)
=
L
−
1
(
x
(
0
)
s
−
A
)
x
(
t
)
=
e
A
t
x
(
0
)
\begin{aligned} \mathcal{L}\left (\dot{\mathbf{x}}(t) \right )&= \mathcal{L}\left (\mathbf{A}\mathbf{x}(t) \right )\\ s\mathbf{X}(s) - \mathbf{x}(0) &= \mathbf{A}\mathbf{X}(s) \\ (s - \mathbf{A})\mathbf{X}(s) &=\mathbf{x}(0) \\ \mathbf{X}(s) &=\frac{\mathbf{x}(0)}{s - \mathbf{A}} \\ \mathcal{L}^{-1}\left (\mathbf{X}\left (s\right ) \right ) &=\mathcal{L}^{-1}\left (\frac{\mathbf{x}(0)}{s - \mathbf{A}}\right ) \\ \mathbf{x}(t) &= e^{\mathbf{A}t}\mathbf{x}(0) \end{aligned}
L(x˙(t))sX(s)−x(0)(s−A)X(s)X(s)L−1(X(s))x(t)=L(Ax(t))=AX(s)=x(0)=s−Ax(0)=L−1(s−Ax(0))=eAtx(0)
现在我们得到了线性时不变系统的状态方程的解
x
(
t
)
=
e
A
t
x
(
0
)
\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}(0)
x(t)=eAtx(0)
该解表述了
t
t
t时刻的状态向量
x
(
t
)
\mathbf{x}(t)
x(t)与初始状态向量
x
(
0
)
\mathbf{x}(0)
x(0)的关系。
e
A
t
e^{\mathbf{A}t}
eAt为矩阵指数函数。该解是在连续时间系统下的解。通常,我们更关心在离散时间系统下的解。假设离散时间系统的时间维度采样率为
Δ
T
\Delta T
ΔT,则
T
T
T时刻与
T
+
Δ
T
T+\Delta T
T+ΔT时刻的状态向量分别为
x ( T ) = e A T x ( 0 ) x ( T + Δ T ) = e A ( T + Δ T ) x ( 0 ) \begin{aligned} \mathbf{x}(T) &= e^{\mathbf{A}T}\mathbf{x}(0) \\ \mathbf{x}(T+\Delta T) &= e^{\mathbf{A}(T+\Delta T)}\mathbf{x}(0) \end{aligned} x(T)x(T+ΔT)=eATx(0)=eA(T+ΔT)x(0)
对上述两式求比值可得
x
(
T
)
{x}(T)
x(T)与
x
(
T
+
Δ
T
)
{x}(T+\Delta T)
x(T+ΔT)的关系
x
(
T
+
Δ
T
)
=
e
A
Δ
T
x
(
T
)
\mathbf{x}(T + \Delta T) = e^{\mathbf{A}{\Delta T}}\mathbf{x}(T)
x(T+ΔT)=eAΔTx(T)
这就是我们熟知的状态转移方程,其中
e
A
Δ
T
e^{\mathbf{A}{\Delta T}}
eAΔT为状态转移矩阵。
矩阵指数函数
现在我们已经推导出了离散时间系统下的状态转移方程,但是矩阵指数函数的存在妨碍我们在实际应用中的使用。矩阵指数函数的解法有许多种,这里简单介绍其中的两种。
泰勒展开
实数的泰勒展开我们已经十分熟悉,
e
a
t
=
1
+
(
a
t
)
+
1
2
!
(
a
t
)
2
+
1
3
!
(
a
t
)
3
+
1
4
!
(
a
t
)
4
+
.
.
.
.
.
.
e^{at} = 1 + (at) + \frac{1}{2!}(at)^2+ \frac{1}{3!}(at)^3+ \frac{1}{4!}(at)^4 + ......
eat=1+(at)+2!1(at)2+3!1(at)3+4!1(at)4+......
矩阵指数函数的泰勒展开的形式类似
e
A
t
=
1
+
(
A
t
)
+
1
2
!
(
A
t
)
2
+
1
3
!
(
A
t
)
3
+
1
4
!
(
A
t
)
4
+
.
.
.
.
.
.
e^{\mathbf{A}t} = 1 + (\mathbf{A}t) + \frac{1}{2!}(\mathbf{A}t)^2+ \frac{1}{3!}(\mathbf{A}t)^3+ \frac{1}{4!}(\mathbf{A}t)^4 + ......
eAt=1+(At)+2!1(At)2+3!1(At)3+4!1(At)4+......
注意,如果用有限项求和的方式代替矩阵指数函数的解,则只能得到近似解。
拉普拉斯变换
使用拉普拉斯变换的形式,同样可以得到矩阵指数函数的解。
e A t = L − 1 ( 1 s I − A ) e^{\mathbf{A}t} = \mathcal{L}^{-1}(\frac{1}{s\mathbf{I}-\mathbf{A}}) eAt=L−1(sI−A1)
Singer模型
常规的运动模型包括CV模型与CA模型。但是,当目标发生机动运动时(简单理解为目标加速度的大小与方向发生了变化),CV模型与CA模型便不再适用。Singer引入了机动频率
α
\alpha
α来修正传统的模型。当机动频率
α
\alpha
α越高,意味着目标加速度变化越快。反之亦然。当机动频率为零时,目标不再处于机动状态。此时,可以根据当前目标加速度的大小,使用CV模型或CA模型对其进行描述。Singer认为目标的加速度的状态方程可以描述为
a
˙
(
t
)
=
−
α
a
(
t
)
+
v
~
(
t
)
\dot{a}(t) = -\alpha a(t) + \tilde{v}(t)
a˙(t)=−αa(t)+v~(t)
其中,
a
(
t
)
a(t)
a(t)与
a
˙
(
t
)
\dot{a}(t)
a˙(t)为t时刻的加速度与加加速度。
v
~
\tilde{v}
v~为一零均值的高斯白噪声。可以发现,状态方程的描述与我们对机动频率的阐述在数值变化趋势上是一致的。当机动频率固定时,目标加加速度的大小取决于当前加速度的大小。由于该方程是一个一阶微分方程,我们同样使用拉普拉斯变换的方法对方程求解。这里暂且忽略噪声的影响,由此得到解的形式为
a
(
t
)
=
e
−
α
t
a
(
0
)
a(t) = e^{-\alpha t}a(0)
a(t)=e−αta(0)
其中,
a
(
0
)
a(0)
a(0)为初始加速度大小。根据上式,加速度的自相关函数为
R
(
τ
)
=
E
[
a
(
t
)
a
(
t
+
τ
)
]
=
σ
e
−
α
∣
τ
∣
R(\tau) = E[a(t)a(t+\tau)] = \sigma e^{-\alpha|\tau|}
R(τ)=E[a(t)a(t+τ)]=σe−α∣τ∣
其中,
σ
\sigma
σ为加速度方差。由此可见,Singer模型下的加速度自相关函数是一个指数函数。
Singer模型的状态方程可以表述为
X
˙
=
A
X
+
V
~
\dot{\mathbf{X}} = \mathbf{A}\mathbf{X}+\tilde{\mathbf{V}}
X˙=AX+V~
其中,状态向量
X
=
[
x
,
x
˙
,
x
¨
]
′
\mathbf{X} = [x,\dot{x},\ddot{x}]^{'}
X=[x,x˙,x¨]′,
x
x
x表示位置信息。系统矩阵为
A
=
[
0
1
0
0
0
1
0
0
−
α
]
\mathbf{A} = \begin{bmatrix} 0 & 1& 0\\ 0 & 0& 1\\ 0 & 0& -\alpha\\ \end{bmatrix}
A=
00010001−α
其中,
α
\alpha
α为机动频率。过程噪声矩阵为
V
~
=
[
0
0
v
~
]
\tilde{\mathbf{V}} = \begin{bmatrix} 0\\ 0\\ \tilde{v}\\ \end{bmatrix}
V~=
00v~
因此,状态方程可以重写为
[
x
˙
x
¨
a
˙
]
=
[
0
1
0
0
0
1
0
0
−
α
]
[
x
x
˙
x
¨
]
+
[
0
0
v
~
]
\begin{bmatrix} \dot{x}\\ \ddot{x}\\ \dot{a}\\ \end{bmatrix} = \begin{bmatrix} 0 & 1& 0\\ 0 & 0& 1\\ 0 & 0& -\alpha\\ \end{bmatrix} \begin{bmatrix} x\\ \dot{x}\\ \ddot{x}\\ \end{bmatrix} + \begin{bmatrix} 0\\ 0\\ \tilde{v}\\ \end{bmatrix}
x˙x¨a˙
=
00010001−α
xx˙x¨
+
00v~
其中,
x
¨
=
a
\ddot{x} = a
x¨=a。使用拉普拉斯变换求解状态方程的解,可得
e
A
t
=
L
−
1
(
1
s
I
−
A
)
=
L
−
1
(
[
s
−
1
0
0
s
−
1
0
0
s
+
α
]
−
1
)
=
L
−
1
(
[
1
s
1
s
2
1
s
2
(
s
+
α
)
0
1
s
1
s
(
s
+
α
)
0
0
1
s
+
α
]
)
=
[
1
t
α
t
−
1
+
e
−
α
t
α
2
0
1
1
−
e
−
α
t
α
0
0
e
−
α
t
]
\begin{aligned} e^{\mathbf{A}t} &= \mathcal{L}^{-1}(\frac{1}{s\mathbf{I}-\mathbf{A}}) \\ &=\mathcal{L}^{-1}(\begin{bmatrix} s & -1& 0\\ 0 & s& -1\\ 0 & 0& s+\alpha\\ \end{bmatrix}^{-1}) \\ &=\mathcal{L}^{-1}( \begin{bmatrix} \frac{1}{s} & \frac{1}{s^2}& \frac{1}{s^2(s+\alpha)}\\ 0 & \frac{1}{s}& \frac{1}{s(s+\alpha)}\\ 0 & 0& \frac{1}{s+\alpha}\\ \end{bmatrix}) \\ &=\begin{bmatrix} 1 & t& \frac{\alpha t - 1 + e^{-\alpha t}}{\alpha ^ 2}\\ 0 & 1& \frac{1-e^{-\alpha t}}{\alpha}\\ 0 & 0& e^{-\alpha t}\\ \end{bmatrix} \end{aligned}
eAt=L−1(sI−A1)=L−1(
s00−1s00−1s+α
−1)=L−1(
s100s21s10s2(s+α)1s(s+α)1s+α1
)=
100t10α2αt−1+e−αtα1−e−αte−αt
在离散时间系统下,假设采样间隔为
T
T
T,则状态转移矩阵为
F
=
[
1
T
α
T
−
1
+
e
−
α
T
α
2
0
1
1
−
e
−
α
T
α
0
0
e
−
α
T
]
\mathbf{F} = \begin{bmatrix} 1 & T& \frac{\alpha T - 1 + e^{-\alpha T}}{\alpha ^ 2}\\ 0 & 1& \frac{1-e^{-\alpha T}}{\alpha}\\ 0 & 0& e^{-\alpha T}\\ \end{bmatrix}
F=
100T10α2αT−1+e−αTα1−e−αTe−αT
CT模型
二维目标运动学方程为
{
x
˙
(
t
)
=
v
(
t
)
c
o
s
ϕ
(
t
)
y
˙
(
t
)
=
v
(
t
)
s
i
n
ϕ
(
t
)
v
˙
(
t
)
=
a
t
(
t
)
ϕ
˙
(
t
)
=
a
n
(
t
)
v
(
t
)
\left\{\begin{matrix} \dot{x}(t) = v(t)cos\phi(t)\\ \dot{y}(t) = v(t)sin\phi(t)\\ \dot{v}(t) = a_{t}(t)\\ \dot{\phi}(t) = \frac{a_{n}(t)}{v(t)}\\ \end{matrix}\right.
⎩
⎨
⎧x˙(t)=v(t)cosϕ(t)y˙(t)=v(t)sinϕ(t)v˙(t)=at(t)ϕ˙(t)=v(t)an(t)
其中,
x
,
y
x,y
x,y为笛卡尔坐标系下的坐标,
v
v
v为空间速度,
a
t
(
t
)
a_{t}(t)
at(t)为切向加速度(沿着速度方向),
a
n
(
t
)
a_{n}(t)
an(t)为法相加速度,
ϕ
\phi
ϕ为航向角。若令
a
n
(
t
)
≠
0
,
a
t
(
t
)
=
0
a_{n}(t)\neq 0,a_{t}(t) = 0
an(t)=0,at(t)=0,则此时目标进行匀速曲线运动。若此时
a
n
(
t
)
a_{n}(t)
an(t)为一常数,则此时目标进行匀速转弯运动。令状态向量
X
=
[
x
,
x
˙
,
y
,
y
˙
]
′
\mathbf{X} = [x,\dot{x},y,\dot{y}]^{'}
X=[x,x˙,y,y˙]′,且角速度
ω
=
ϕ
˙
\omega = \dot{\phi}
ω=ϕ˙。
x
x
x维度的加速度可以表示为
x
¨
(
t
)
=
a
t
(
t
)
c
o
s
ϕ
(
t
)
−
ω
v
(
t
)
s
i
n
ϕ
(
t
)
=
−
ω
y
˙
(
t
)
\begin{aligned} \ddot{x}(t) &= a_{t}(t)cos\phi(t) - \omega v(t)sin\phi(t) \\ &=-\omega\dot{y}(t) \end{aligned}
x¨(t)=at(t)cosϕ(t)−ωv(t)sinϕ(t)=−ωy˙(t)
同理
y
¨
(
t
)
=
ω
x
˙
(
t
)
\ddot{y}(t) = \omega \dot{x}(t)
y¨(t)=ωx˙(t)
因此系统矩阵可以写为
A
=
[
0
1
0
0
0
0
0
−
ω
0
0
0
1
0
ω
0
0
]
\mathbf{A} = \begin{bmatrix} 0 & 1& 0& 0\\ 0 & 0& 0& -\omega\\ 0 & 0& 0& 1\\ 0 & \omega& 0& 0\\ \end{bmatrix}
A=
0000100ω00000−ω10
使用拉普拉斯变换求解状态方程的解,可得
e
A
t
=
L
−
1
(
1
s
I
−
A
)
=
L
−
1
(
[
s
−
1
0
0
0
s
0
ω
0
0
s
−
1
0
−
ω
0
s
]
−
1
)
=
L
−
1
(
[
1
s
1
s
2
+
ω
2
0
−
ω
s
3
+
s
ω
2
0
s
s
2
+
ω
2
0
−
ω
s
2
+
ω
2
0
ω
s
3
+
s
ω
2
1
s
1
s
2
+
ω
2
0
ω
s
2
+
ω
2
0
s
s
2
+
ω
2
]
)
=
[
1
s
i
n
ω
t
ω
0
−
1
−
c
o
s
ω
t
ω
0
c
o
s
ω
t
0
−
s
i
n
ω
t
0
1
−
c
o
s
ω
t
ω
1
s
i
n
ω
t
ω
0
s
i
n
ω
t
0
c
o
s
ω
t
]
\begin{aligned} e^{\mathbf{A}t} &= \mathcal{L}^{-1}(\frac{1}{s\mathbf{I}-\mathbf{A}}) \\ &=\mathcal{L}^{-1}(\begin{bmatrix} s & -1& 0& 0\\ 0 & s& 0& \omega\\ 0 & 0& s& -1\\ 0 & -\omega& 0& s\\ \end{bmatrix}^{-1}) \\ &= \mathcal{L}^{-1}(\begin{bmatrix} \frac{1}{s} & \frac{1}{s^2+\omega^2}& 0& \frac{-\omega}{s^3+s\omega^2}\\ 0 & \frac{s}{s^2+\omega^2}& 0& \frac{-\omega}{s^2+\omega^2}\\ 0 & \frac{\omega}{s^3+s\omega^2}& \frac{1}{s}& \frac{1}{s^2+\omega^2}\\ 0 & \frac{\omega}{s^2+\omega^2}& 0& \frac{s}{s^2+\omega^2}\\ \end{bmatrix}) \\ &= \begin{bmatrix} 1 & \frac{sin\omega t}{\omega}& 0& -\frac{1-cos\omega t}{\omega}\\ 0 & cos\omega t& 0& -sin\omega t\\ 0 & \frac{1-cos\omega t}{\omega}& 1&\frac{sin\omega t}{\omega}\\ 0 & sin\omega t& 0& cos\omega t\\ \end{bmatrix} \end{aligned}
eAt=L−1(sI−A1)=L−1(
s000−1s0−ω00s00ω−1s
−1)=L−1(
s1000s2+ω21s2+ω2ss3+sω2ωs2+ω2ω00s10s3+sω2−ωs2+ω2−ωs2+ω21s2+ω2s
)=
1000ωsinωtcosωtω1−cosωtsinωt0010−ω1−cosωt−sinωtωsinωtcosωt
在离散时间系统下,假设采样间隔为
T
T
T,则状态转移矩阵为
F
=
[
1
s
i
n
ω
T
ω
0
−
1
−
c
o
s
ω
T
ω
0
c
o
s
ω
T
0
−
s
i
n
ω
T
0
1
−
c
o
s
ω
T
ω
1
s
i
n
ω
T
ω
0
s
i
n
ω
T
0
c
o
s
ω
T
]
\mathbf{F} = \begin{bmatrix} 1 & \frac{sin\omega T}{\omega}& 0& -\frac{1-cos\omega T}{\omega}\\ 0 & cos\omega T& 0& -sin\omega T\\ 0 & \frac{1-cos\omega T}{\omega}& 1&\frac{sin\omega T}{\omega}\\ 0 & sin\omega T& 0& cos\omega T\\ \end{bmatrix}
F=
1000ωsinωTcosωTω1−cosωTsinωT0010−ω1−cosωT−sinωTωsinωTcosωT