任务介绍
在本任务中，您将使用可用的测量值和运动模型来递归估计车辆沿轨迹的位置。车辆有了非常简单的LIDAR传感器，可以返回与环境中各个地标相对应的范围（range）和方位测量值（bearing）。假定地标的所有位置是已知的并假设他们已知的数据关联，即哪个度量属于哪个界标。
数据链接：https://pan.baidu.com/s/10jQlD0fOyJ81NWP2FOoxAA
提取码：dzkl
运动和测量模型
运动模型
车辆运动模型将线速度和角速度里程计读数作为输入，并输出车辆的状态（即2D状态）：
x k = x k − 1 + T [ cos θ k − 1 0 sin θ k − 1 0 0 1 ] ( [ v k ω k ] + w k ) , w k = N ( 0 , Q ) \begin{aligned}\mathbf{x}_{k} &= \mathbf{x}_{k-1} + T\begin{bmatrix}\cos\theta_{k-1} &0 \\\sin\theta_{k-1} &0 \\0 &1 \end{bmatrix}\left(\begin{bmatrix}v_k \\\omega_k\end{bmatrix} + \mathbf{w}_k\right)\, , \, \, \, \, \, \mathbf{w}_k = \mathcal{N}\left(\mathbf{0}, \mathbf{Q}\right)\end{aligned} xk=xk−1+T⎣⎡cosθk−1sinθk−10001⎦⎤([vkωk]+wk),wk=N(0,Q)
- x k = [ x y θ ] T \mathbf{x}_k = \left[ x \, y \, \theta \right]^T xk=[xyθ]T 是当前车辆的状态
- v k v_k vk 和 ω k \omega_k ωk 是线速度和角速度的里程表读数，并作为模型的输入值
过程噪声 w k \mathbf{w}_k wk 服从正太分布（均值为0）和常数协方差 Q \mathbf{Q} Q.
测量模型
测量模型将车辆的当前姿态与激光雷达范围和方位测量值相关联 y k l = [ r ϕ ] T \mathbf{y}^l_k = \left[r \, \phi \right]^T ykl=[rϕ]T.
y k l = [ ( x l − x k − d cos θ k ) 2 + ( y l − y k − d sin θ k ) 2 a r c t a n ( x l − x k − d cos θ k , y l − y k − d sin θ k ) − θ k ] + n k l , n k l = N ( 0 , R ) \begin{aligned} \mathbf{y}^l_k = \begin{bmatrix} \sqrt{(x_l - x_k - d\cos\theta_{k})^2 + (y_l - y_k - d\sin\theta_{k})^2} \\ arctan\left(x_l - x_k - d\cos\theta_{k},y_l - y_k - d\sin\theta_{k}\right) - \theta_k \end{bmatrix} + \mathbf{n}^l_k \, , \, \, \, \, \, \mathbf{n}^l_k = \mathcal{N}\left(\mathbf{0}, \mathbf{R}\right) \end{aligned} ykl=[(xl−xk−dcosθk)2+(yl−yk−dsinθk)2arctan(xl−xk−dcosθk,yl−yk−dsinθk)−θk]+nkl,nkl=N(0,R)
- x l x_l xl 和 y l y_l yl是 标记值 l l l 的真实值
- x k x_k xk 、 y k y_k yk 和 θ k \theta_{k} θk 表示当前车辆的状态
- d d d 机器人中心和前轮激光测距仪（laser rangefinde)的距离 (LIDAR)
标记值的测量噪声 n k l \mathbf{n}^l_k nkl 服从正太分布（均值为0）和常数协方差 R \mathbf{R} R.
准备
由于上述模型是非线性的，因此我们建议使用扩展卡尔曼滤波器（EKF）作为状态估计器。具体来说，需要实现两个环节：
- 预测环节，使用测距法测量和运动模型在给定的时间步长产生状态和协方差估计值
- 修正环节，使用范围和方位角激光雷达提供的测量值以纠正状态和他们的协方差估计
加载数据
import pickle
import numpy as np
import matplotlib.pyplot as plt
import math
with open('data/data.pickle', 'rb') as f:
data = pickle.load(f)
t = data['t'] # timestamps [s]
x_init = data['x_init'] # initial x position [m]
y_init = data['y_init'] # initial y position [m]
th_init = data['th_init'] # initial theta position [rad]
# input signal
v = data['v'] # translational velocity input [m/s]
om = data['om'] # rotational velocity input [rad/s]
# bearing and range measurements, LIDAR constants
b = data['b'] # bearing to each landmarks center in the frame attached to the laser [rad]
r = data['r'] # range measurements [m]
l = data['l'] # x,y positions of landmarks [m]
d = data['d'] # distance between robot center and laser rangefinder [m]
注意这里的参数d，即机器人中心到LIDAR的距离放到了d变量中。
Ground Truth
如果实现了该任务，相应的轨迹和角度为：
注意方向值以弧度变成了 [ − π , π ] \left[-\pi,\pi\right] [−π,π]的范围.
初始化参数
现在已经加载了数据，我们可以开始为求解程序进行设置了。设计滤波器的最重要方面之一是确定输入和测量噪声协方差矩阵，以及初始状态和协方差值，设置值如下：
v_var = 0.01 # translation velocity variance
om_var = 0.1 # rotational velocity variance
r_var = 0.1 # range measurements variance
b_var = 0.1 # bearing measurement variance
Q_km = np.diag([v_var, om_var]) # input noise covariance
cov_y = np.diag([r_var, b_var]) # measurement noise covariance
x_est = np.zeros([len(v), 3]) # estimated states, x, y, and theta
P_est = np.zeros([len(v), 3, 3]) # state covariance matrices
x_est[0] = np.array([x_init, y_init, th_init]) # initial state
P_est[0] = np.diag([1, 1, 0.1]) # initial state covariance
注意: 有必要调整测量噪声的方差 r_var, b_var 使得滤波器表现更好
为了使方向的估计值和方位测量值一致, 有必要把所有测量的 θ \theta θ 值 换到 ( − π , π ] (-\pi , \pi] (−π,π] 范围.
# Wraps angle to (-pi,pi] range
def wraptopi(x):
if x > np.pi:
x = x - (np.floor(x / (2 * np.pi)) + 1) * 2 * np.pi
elif x < -np.pi:
x = x + (np.floor(x / (-2 * np.pi)) + 1) * 2 * np.pi
return x
这个函数可以将一个单独的值换算到范围内，还考虑过这种写法
def wraptopi1(x):
x[x>np.pi]=x[x>np.pi] - (np.floor(x[x>np.pi] / (2 * np.pi)) + 1) * 2 * np.pi
x[x<-np.pi]= x[x<-np.pi] + (np.floor(x[x<-np.pi] / (-2 * np.pi)) + 1) * 2 * np.pi
return x
这样定义函数，就可以把数组内的所有值换到该范围内。以下的函数都是基于第一个定义方法写的。
修正环节
- 在 x ˇ k \mathbf{\check{x}}_{k} xˇk计算测量模型的雅可比矩阵
y k l = h ( x k , n k l ) H k = ∂ h ∂ x k ∣ x ˇ k , 0 , M k = ∂ h ∂ n k ∣ x ˇ k , 0 . \begin{aligned} \mathbf{y}^l_k = &\mathbf{h}(\mathbf{x}_{k}, \mathbf{n}^l_k) \\\\ \mathbf{H}_{k} = \frac{\partial \mathbf{h}}{\partial \mathbf{x}_{k}}\bigg|_{\mathbf{\check{x}}_{k},0}& \, , \, \, \, \, \mathbf{M}_{k} = \frac{\partial \mathbf{h}}{\partial \mathbf{n}_{k}}\bigg|_{\mathbf{\check{x}}_{k},0} \, . \end{aligned} ykl=Hk=∂xk∂h∣∣∣∣xˇk,0h(xk,nkl),Mk=∂nk∂h∣∣∣∣xˇk,0.
- 计算卡尔曼增益
K k = P ˇ k H k T ( H k P ˇ k H k T + M k R k M k T ) − 1 \begin{aligned} \mathbf{K}_k &= \mathbf{\check{P}}_k \mathbf{H}_k^T \left(\mathbf{H}_k \mathbf{\check{P}}_k \mathbf{H}_k^T + \mathbf{M}_k \mathbf{R}_k \mathbf{M}_k^T \right)^{-1} \end{aligned} Kk=PˇkHkT(HkPˇkHkT+MkRkMkT)−1
- 修正预测状态
y ˇ k l = h ( x ˇ k , 0 ) x ^ k = x ˇ k + K k ( y k l − y ˇ k l ) \begin{aligned} \mathbf{\check{y}}^l_k &= \mathbf{h}\left(\mathbf{\check{x}}_k, \mathbf{0}\right) \\ \mathbf{\hat{x}}_k &= \mathbf{\check{x}}_k + \mathbf{K}_k \left(\mathbf{y}^l_k - \mathbf{\check{y}}^l_k\right) \end{aligned} yˇklx^k=h(xˇk,0)=xˇk+Kk(ykl−yˇkl) - 修正协方差
P
^
k
=
(
I
−
K
k
H
k
)
P
ˇ
k
\begin{aligned} \mathbf{\hat{P}}_k &= \left(\mathbf{I} - \mathbf{K}_k \mathbf{H}_k \right)\mathbf{\check{P}}_k \end{aligned}
P^k=(I−KkHk)Pˇk
把车辆想象成一根棍子，
(
x
k
,
y
k
)
(x_{k},y_{k})
(xk,yk)是车辆的中心位置，车辆最前面是雷达的位置，
(
x
l
,
y
l
)
(x_{l},y_{l})
(xl,yl)是车辆要取的位置，即landmark标记点，那么很容易根据正切值算出来这个变化的角度。
因此，这里详细推导一下
H
k
H_{k}
Hk和
M
k
M_{k}
Mk，
H
k
H_{k}
Hk是h对
x
k
x_{k}
xk求偏导，而h就是
y
k
l
{y}^l_k
ykl,
y
k
l
{y}^l_k
ykl是(2,1)的矩阵，分别对里面的
x
k
x_{k}
xk、
y
k
y_{k}
yk、
θ
k
\theta_{k}
θk求偏导那么可以得到(2,3)的矩阵，即
H
k
H_{k}
Hk的维数是(2,3)
y
k
l
=
[
(
x
l
−
x
k
−
d
cos
θ
k
)
2
+
(
y
l
−
y
k
−
d
sin
θ
k
)
2
a
r
c
t
a
n
(
x
l
−
x
k
−
d
cos
θ
k
,
y
l
−
y
k
−
d
sin
θ
k
)
−
θ
k
]
+
n
k
l
,
n
k
l
=
N
(
0
,
R
)
\begin{aligned} \mathbf{y}^l_k = \begin{bmatrix} \sqrt{(x_l - x_k - d\cos\theta_{k})^2 + (y_l - y_k - d\sin\theta_{k})^2} \\ arctan\left(x_l - x_k - d\cos\theta_{k},y_l - y_k - d\sin\theta_{k}\right) - \theta_k \end{bmatrix} + \mathbf{n}^l_k \, , \, \, \, \, \, \mathbf{n}^l_k = \mathcal{N}\left(\mathbf{0}, \mathbf{R}\right) \end{aligned}
ykl=[(xl−xk−dcosθk)2+(yl−yk−dsinθk)2arctan(xl−xk−dcosθk,yl−yk−dsinθk)−θk]+nkl,nkl=N(0,R)
为了表示方便，采用三个中间变量A、B、C来简化过程。
A
=
x
l
−
x
k
−
d
cos
θ
k
A=x_l - x_k - d\cos\theta_{k}
A=xl−xk−dcosθk
B
=
y
l
−
y
k
−
d
sin
θ
k
B=y_l - y_k - d\sin\theta_{k}
B=yl−yk−dsinθk
C
=
(
x
l
−
x
k
−
d
cos
θ
k
)
2
+
(
y
l
−
y
k
−
d
sin
θ
k
)
2
=
A
2
+
B
2
C=\sqrt{(x_l - x_k - d\cos\theta_{k})^2 + (y_l - y_k - d\sin\theta_{k})^2}=\sqrt{A^2 + B^2}
C=(xl−xk−dcosθk)2+(yl−yk−dsinθk)2=A2+B2
因此求导后的
H
k
H_{k}
Hk为
H k = [ − A / C − B / C ( d / C ) ∗ ( A sin θ − B cos θ ) B / C 2 − A / C 2 ( − d / C 2 ) ∗ cos θ ( A + B ) − 1 ] \begin{aligned} \mathbf{H}^k = \begin{bmatrix} -A/C & -B/C & (d/C) *(A\sin\theta-B\cos\theta)\\ B/C^2 & -A/C^2 & (-d/C^2)*\cos\theta(A+B)-1 \end{bmatrix} \end{aligned} Hk=[−A/CB/C2−B/C−A/C2(d/C)∗(Asinθ−Bcosθ)(−d/C2)∗cosθ(A+B)−1]
矩阵的维度
思考一下矩阵的维度：
| 变量 | 维度 |
|---|---|
| x_k | (3,1) |
| y_k | (2,1) |
| K_k | (3,2) |
| F_k | (3,3) |
| L_k | (3,2) |
| Q_k | (2,2) |
| P_k | (3,3) |
| H_k | (2,3) |
| M_k | (2,2) |
| X_est | (501,3) |
| P_est | (501,3,3) |
def measurement_update(lk, rk, bk, P_check, x_check):
theta = x_check[:,2]
bk = wraptopi(bk)
theta= wraptopi(theta)
A = lk[0] - x_check[:,0] - d * math.cos(theta)
B = lk[1] - x_check[:,1] - d * math.sin(theta)
C = (A ** 2 + B ** 2) ** (0.5)
# 1. Compute measurement Jacobian
M = np.identity(2)
H = np.array([[-A/C,-B/C,(d/C)*(A*math.sin(theta)- B*math.cos(theta))],
[B/(C**2),-A/(C**2),((-d/(C**2))*math.cos(theta)*(A+B))-1]])
H = H.reshape(2,3) ## the dimension of H is (2,3,1),so you need to reshape it.
# 2. Compute Kalman Gain
K = P_check @ H.T @ np.linalg.inv(H @ P_check @ H.T + M @cov_y @ M.T )
# 3. Correct predicted state (remember to wrap the angles to [-pi,pi])
h = math.atan2(B,A) - theta
h = wraptopi(h)
ss = K @ np.array([rk-C,bk-h])
x_check = x_check + ss.reshape(1,3)
theta = x_check[:,2]
theta= wraptopi(theta)
# 4. Correct covariance
P_check = (np.identity(3)-K @ H) @ P_check
return x_check, P_check
预测环节
现在，实现主滤波器循环，使用提供的运动模型定义EKF的预测步骤：
x
ˇ
k
=
f
(
x
^
k
−
1
,
u
k
−
1
,
0
)
P
ˇ
k
=
F
k
−
1
P
^
k
−
1
F
k
−
1
T
+
L
k
−
1
Q
k
−
1
L
k
−
1
T
.
\begin{aligned} \mathbf{\check{x}}_k &= \mathbf{f}\left(\mathbf{\hat{x}}_{k-1}, \mathbf{u}_{k-1}, \mathbf{0} \right) \\ \mathbf{\check{P}}_k &= \mathbf{F}_{k-1}\mathbf{\hat{P}}_{k-1}\mathbf{F}_{k-1}^T + \mathbf{L}_{k-1}\mathbf{Q}_{k-1}\mathbf{L}_{k-1}^T \, . \end{aligned}
xˇkPˇk=f(x^k−1,uk−1,0)=Fk−1P^k−1Fk−1T+Lk−1Qk−1Lk−1T.
其中
F
k
−
1
=
∂
f
∂
x
k
−
1
∣
x
^
k
−
1
,
u
k
,
0
,
L
k
−
1
=
∂
f
∂
w
k
∣
x
^
k
−
1
,
u
k
,
0
.
\begin{aligned} \mathbf{F}_{k-1} = \frac{\partial \mathbf{f}}{\partial \mathbf{x}_{k-1}}\bigg|_{\mathbf{\hat{x}}_{k-1},\mathbf{u}_{k},0} \, , \, \, \, \, \mathbf{L}_{k-1} = \frac{\partial \mathbf{f}}{\partial \mathbf{w}_{k}}\bigg|_{\mathbf{\hat{x}}_{k-1},\mathbf{u}_{k},0} \, . \end{aligned}
Fk−1=∂xk−1∂f∣∣∣∣x^k−1,uk,0,Lk−1=∂wk∂f∣∣∣∣x^k−1,uk,0.
那么同样求一下
F
k
−
1
{F}_{k-1}
Fk−1,由于f’就是
x
k
x_{k}
xk,类似地，f是(3,1)的矩阵，对里面的
x
k
x_{k}
xk、
y
k
y_{k}
yk、
θ
k
\theta_{k}
θk求偏导那么可以得到(3,3)的矩阵，即
F
k
−
1
F_{k-1}
Fk−1的维数是(3,3)
x
k
=
x
k
−
1
+
T
[
cos
θ
k
−
1
0
sin
θ
k
−
1
0
0
1
]
(
[
v
k
ω
k
]
+
w
k
)
,
w
k
=
N
(
0
,
Q
)
\begin{aligned}\mathbf{x}_{k} &= \mathbf{x}_{k-1} + T\begin{bmatrix}\cos\theta_{k-1} &0 \\\sin\theta_{k-1} &0 \\0 &1 \end{bmatrix}\left(\begin{bmatrix}v_k \\\omega_k\end{bmatrix} + \mathbf{w}_k\right)\, , \, \, \, \, \, \mathbf{w}_k = \mathcal{N}\left(\mathbf{0}, \mathbf{Q}\right)\end{aligned}
xk=xk−1+T⎣⎡cosθk−1sinθk−10001⎦⎤([vkωk]+wk),wk=N(0,Q)
由于
x
k
−
1
{x}_{k-1}
xk−1是
(
x
k
−
1
,
y
k
−
1
,
θ
k
−
1
)
T
({x}_{k-1},{y}_{k-1},{\theta}_{k-1})^T
(xk−1,yk−1,θk−1)T,那么带入后
x
k
=
[
x
k
−
1
+
T
cos
θ
k
−
1
∗
(
v
k
+
w
k
)
y
k
−
1
+
T
sin
θ
k
−
1
∗
(
v
k
+
w
k
)
T
∗
(
v
k
+
w
k
)
]
\begin{aligned}\mathbf{x}_{k} &= \begin{bmatrix}{x}_{k-1}+T\cos\theta_{k-1} *(v_k+ \mathbf{w}_k) \\ {y}_{k-1}+T\sin\theta_{k-1} *(v_k+ \mathbf{w}_k) \\T*(v_k+ \mathbf{w}_k)\end{bmatrix}\end{aligned}
xk=⎣⎡xk−1+Tcosθk−1∗(vk+wk)yk−1+Tsinθk−1∗(vk+wk)T∗(vk+wk)⎦⎤
那么每一列分别对
x
k
−
1
,
y
k
−
1
,
θ
k
−
1
{x}_{k-1},{y}_{k-1},{\theta}_{k-1}
xk−1,yk−1,θk−1求导可以得到
F
k
−
1
F_{k-1}
Fk−1
F
k
−
1
=
[
1
0
−
T
∗
v
k
∗
sin
θ
k
−
1
0
1
T
∗
v
k
∗
cos
θ
k
−
1
0
0
1
]
\begin{aligned}\mathbf{F}_{k-1} &= \begin{bmatrix}1&0&-T*v_{k}*\sin\theta_{k-1}\\ 0&1&T*v_{k}*\cos\theta_{k-1} \\0&0&1\end{bmatrix}\end{aligned}
Fk−1=⎣⎡100010−T∗vk∗sinθk−1T∗vk∗cosθk−11⎦⎤
类似地，对
w
k
\mathbf{w}_k
wk求导显而易见,就是前面的系数
L
k
−
1
=
T
[
cos
θ
k
−
1
0
sin
θ
k
−
1
0
0
1
]
\begin{aligned}\mathbf{L}_{k-1} &=T\begin{bmatrix}\cos\theta_{k-1} &0 \\\sin\theta_{k-1} &0 \\0 &1 \end{bmatrix}\end{aligned}
Lk−1=T⎣⎡cosθk−1sinθk−10001⎦⎤
#### 5. Main Filter Loop #######################################################################
for k in range(1, len(t)): # start at 1 because we've set the initial prediciton
delta_t = t[k] - t[k - 1] # time step (difference between timestamps)
x_check = x_est[k - 1,:]
P_check = P_est[k - 1,:,:]
theta = x_check[2]
theta = wraptopi(theta)
# 1. Update state with odometry readings (remember to wrap the angles to [-pi,pi])
cc = delta_t*np.array([[math.cos(theta),0],[math.sin(theta),0],[0,1]])@ np.array([[v[k]],[om[k]]])
x_check = x_check.reshape(1,3)
x_check = x_check + cc.T
theta = x_check[0,2]
theta = wraptopi(theta)
# 2. Motion model jacobian with respect to last state
F_km = np.array([[1,0,-delta_t * v[k] *math.sin(theta)],[0,1,delta_t * v[k]* math.cos(theta)],[0,0,1]])
# 3. Motion model jacobian with respect to noise
L_km = delta_t * np.array([[math.cos(theta),0],[math.sin(theta),0],[0,1]])
# 4. Propagate uncertainty
P_check = F_km @ P_check @ F_km.T + L_km @ Q_km @L_km.T
# 5. Update state estimate using available landmark measurements
for i in range(len(r[k])):
x_check, P_check = measurement_update(l[i], r[k, i], b[k, i], P_check, x_check)
# Set final state predictions for timestep
x_est[k, 0] = x_check[0,0]
x_est[k, 1] = x_check[0,1]
x_est[k, 2] = x_check[0,2]
P_est[k, :, :] = P_check
那么显示下结果
e_fig = plt.figure()
ax = e_fig.add_subplot(111)
ax.plot(x_est[:, 0], x_est[:, 1])
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_title('Estimated trajectory')
plt.show()
e_fig = plt.figure()
ax = e_fig.add_subplot(111)
ax.plot(t[:], x_est[:, 2])
ax.set_xlabel('Time [s]')
ax.set_ylabel('theta [rad]')
ax.set_title('Estimated trajectory')
plt.show()
不好意思走错了，应该是这样的
如果想实时显示运动轨迹，可以运行一下代码
ax = [] # 定义一个 x 轴的空列表用来接收动态的数据
ay = [] # 定义一个 y 轴的空列表用来接收动态的数据
plt.ion() # 开启一个画图的窗口
for i in range(501): # 遍历0-99的值
ax.append(x_est[i, 0]) # 添加 x 到 x 轴的数据中
ay.append(x_est[i, 1]) # 添加 y 到 y 轴的数据中
plt.clf() # 清除之前画的图
plt.plot(ax,ay) # 画出当前 ax 列表和 ay 列表中的值的图形
plt.pause(0.1) # 暂停一秒
plt.ioff() # 关闭画图的窗口
同样，随时间变化的角度也可以显示
ax = [] # 定义一个 x 轴的空列表用来接收动态的数据
ay = [] # 定义一个 y 轴的空列表用来接收动态的数据
plt.ion() # 开启一个画图的窗口
for i in range(501): # 遍历0-99的值
ax.append(i) # 添加 时间t 到 x 轴的数据中
ay.append(x_est[i, 2]) # 添加 theta 到 y 轴的数据中
plt.clf() # 清除之前画的图
plt.plot(ax,ay) # 画出当前 ax 列表和 ay 列表中的值的图形
plt.pause(0.1) # 暂停一秒
plt.ioff() # 关闭画图的窗口
结果路径差不太多，但是噪声太多了，那么按照一下思路，更改一下方差
v_var不做更改（v_var=0.01）- 15>
om_var>5 - 1>
r_var>0.001 - 1>
b_var>0.001
如取 om_var=5，r_var=0.001 ，b_var=0.001
总结
每次更新完，角度应该变换在[-pi,pi]之间，而且一些矩阵的维数（主要是一维）应该用reshape命令，防止出错，尤其是一下情况：
如定义一个矩阵
x
1
x_{1}
x1,
import numpy as np
x1 = np.array([50,0.69,1.63])
x1
结果为：array([50. , 0.69, 1.63])
那么看一下该向量（vector）的维数，结果为：(3,)
x1.shape
如果此时定义另一个变量 x 2 x_{2} x2
x2=np.array([[10,30,40]]).T
x2
那么结果为,其维数为(3, 1)
array([[10],
[30],
[40]])
做矩阵乘法，结果为：array([585.9])
np.dot(x1,x2)
(3,)的向量可以和(3,3)矩阵相乘,（矩阵维数为(1,3)肯定可以和(3,3)相乘,略）但有时候你想定义一个(1,3)的矩阵，区分这两种写法：
a=np.array([[1,2,3]])
b=np.array([1,2,3])
根据上面的演示，shape分别为：(1, 3), (3,)
如果忘记了，只需要c=b.reshape(1,3),便可以把矩阵变成(1, 3)
而且，如果你让a+b，可以得到array([[2, 4, 6]])，便以a的维数为基础扩充，但如果你想以b的维度为基础，可以让b进行b=b.reshape(3,)操作。
在上述的H矩阵,(2,3,1)这样的,你需要的是(2,3),因此需要reshape一下）,出现该问题的原因是有一些是数，有一些是包含一个元素的数组，相加的时候变成了含有一个元素的数组，这样会导致维数增加，解决办法一是在定义的时候多检查想要取出的是元素还是数组，例如
x_check[:,0]是array([50.])而lk[0]是-20，可以通过x_check[:,0][0]来取出50，或者就不考虑数组和元素的区别，在计算后强制改变reshape数组的维数。
总之，为了防止出错，多用reshape处理一维的向量、元素、数组运算，或者运行后将一些参数强制reshape一下和你推导的矩阵维数对应。
扫一扫
1267
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
