辨识实例:最小二乘法拟合圆的标准方程
以下以拟合圆的方程来演示最小二乘二阶辨识的实例,已知圆的通式为
(
x
−
a
)
2
+
(
y
−
b
)
2
=
r
2
\begin{equation} \begin{split} (x-a)^2 + (y-b)^2 = r^2 \end{split} \end{equation}
(x−a)2+(y−b)2=r2
由于最小二乘法不能直接辨识二阶以上的方程,所以我们需要对该方程进行处理
(
x
−
a
)
2
+
(
y
−
b
)
2
=
r
2
x
2
−
2
a
x
+
a
2
+
y
2
−
2
b
y
+
b
2
=
r
2
x
2
+
y
2
−
2
a
x
−
2
b
y
+
a
2
+
b
2
−
r
2
=
0
\begin{equation} \begin{split} (x-a)^2 + (y-b)^2 &= r^2 \\ x^2 - 2ax + a^2 + y^2 - 2by + b^2 &= r^2 \\ x^2 + y^2 - 2ax - 2by + a^2+ b^2 - r^2 &= 0 \end{split} \end{equation}
(x−a)2+(y−b)2x2−2ax+a2+y2−2by+b2x2+y2−2ax−2by+a2+b2−r2=r2=r2=0
上式中
a
2
a^2
a2、
b
2
b^2
b2以及
r
2
r^2
r2在应用时都是常数项,所以可以将其当作一个常数项
c
c
c来看待。待辨识方程即为
f
(
a
,
b
,
c
)
=
x
2
+
y
2
−
2
a
x
−
2
b
y
+
c
\begin{equation} \begin{split} f(a,b,c) = x^2 + y^2 - 2ax - 2by + c \end{split} \end{equation}
f(a,b,c)=x2+y2−2ax−2by+c
这里将上式化简为最小二乘的矩阵求解形式如下
f
(
θ
)
=
Z
m
=
[
−
2
x
−
2
y
E
]
[
a
b
c
]
=
H
m
θ
^
\begin{equation} \begin{split} f(\theta) = Z_m &= \begin{bmatrix} -2x & -2y & E \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} \\ &=H_m \hat \theta \end{split} \end{equation}
f(θ)=Zm=[−2x−2yE]
abc
=Hmθ^
式中
E
E
E为单位矩阵,其他字母均为矩阵,需要保证维度一致。
使用最小二乘的求解公式对其进行求解
θ
^
=
[
a
^
b
^
c
^
]
=
H
m
−
1
Z
m
\begin{equation} \begin{split} \hat \theta = \begin{bmatrix} \hat a \\ \hat b \\ \hat c \end{bmatrix}= H_m^{-1}Z_m \end{split} \end{equation}
θ^=
a^b^c^
=Hm−1Zm
求解得到的
a
^
\hat a
a^与
b
^
\hat b
b^即是辨识得到的圆的圆心坐标
(
x
,
y
)
(x,y)
(x,y),而半径
r
r
r则由于我们对其进行了合并处理,所以需要进一步处理才能得到,如下
a
2
+
b
2
+
r
2
=
c
r
2
=
c
−
a
2
−
b
2
r
=
c
−
a
2
−
b
2
\begin{equation} \begin{split} a^2 + b^2 + r^2 &= c \\ r^2 &= c - a^2 - b^2 \\ r &= \sqrt{c - a^2 - b^2} \end{split} \end{equation}
a2+b2+r2r2r=c=c−a2−b2=c−a2−b2
这里以圆心为
(
0
,
0
)
(0,0)
(0,0),半径
r
=
3
r = 3
r=3m作为圆的真实参数。以下是该实例对应的matlab代码
% 生成没有噪音的圆的轨迹点
t = linspace(0, 2*pi, 100);
r = 3;
x = r*cos(t);
y = r*sin(t);
figure;
plot(x,y); hold on;
axis equal
% 采集带噪音的真实坐标数据
z_x = x + 0.2*randn(1,length(t));
z_y = y + 0.2*randn(1,length(t));
plot(z_x,z_y,'r.');
% 构造最小二乘形式
zm = -(z_x.^2 + z_y.^2)';
H = [-2*z_x', -2*z_y', ones(1,length(z_x))'];
% 最小二乘求解
hat = H\zm;
x_hat = hat(1);
y_hat = hat(2);
r_hat = sqrt(-hat(3)-hat(1)^2-hat(2)^2);
fprintf("圆心位置为:(%f, %f)\n",x_hat,y_hat);
fprintf("圆心半径为:%f m\n",r_hat);
%画图验证
x_fit = x_hat + r_hat*cos(t);
y_fit = y_hat + r_hat*sin(t);
plot(x_fit, y_fit,'r-');
legend('真实轨迹','采集数据','拟合轨迹','Location','best')
从结果我们可以看出
- 最小二乘并不是无偏的,它是将所有的误差平摊到每个采集的数据点上,它只能起到类似滤波器的效果,对于那些对精度要求较高的参数(小数点后两位以上),使用最小二乘较难满足精度要求
- 最小二乘是最简单的辨识方法,它是基于假设噪音为高斯(正态)分布时的白噪音时才能获得较好的效果,如果不是白噪音,则可能需要使用到卡尔曼滤波、神经网络辨识、粒子群辨识、粒子滤波等高级辨识方法