空间球拟合与椭球拟合

实验科学离不开拟合,考虑到二次曲面拟合中球拟合和椭球拟合较常使用,并且球拟合有简单的最小二乘方法,因此这里稍微整理一下二者拟合相关的数学工具与处理思路。

球拟合

处理球的拟合时不会遇到处理椭球问题中的非线性问题,可以直接使用最小二乘法来得出结果。

三维空间中的任意一个球面可以表示成:

$$ (x-x_0)^2+(y-y_0)^2+(z-z_0)^2=r^2 $$

对该式展开并移相得到:

$$ x^2+y^2+z^2=2x_0x+2y_0y+2z_0z+(r^2-x_0^2-y_0^2-z_0^2) $$

不妨记:

$$ A = \begin{pmatrix} 2x_1 & 2y_1 & 2z_1 & 1 \\ 2x_2 & 2y_2 & 2z_2 & 1 \\ \vdots & \vdots & \vdots & \vdots \\ 2x_n & 2y_n & 2z_n & 1 \\ \end{pmatrix}, \vec{b} = \begin{pmatrix} x_0 \\ y_0 \\ z_0 \\ r^2 - x_{0}^2 - y_{0}^2 - z_{0}^2 \end{pmatrix}, \vec{Y} = \begin{pmatrix} x_1^2 + y_1^2 + z_1^2 \\ x_2^2 + y_2^2 + z_2^2 \\ \vdots \\ x_{n}^2 + y_{n}^2 + z_{n}^2 \end{pmatrix} $$

其中$(x_1,y_1,z_1),(x_2,y_2,z_2),...,(x_n,y_n,z_n)$分别是$n$个散点的坐标。这样问题就转化为了求下述超定方程的解:

$$ A\vec{b}=\vec{Y} $$

由于这是一个线性问题,我们可以直接写出该超定方程的解:

$$ \vec{b_0}=(A^\intercal A)^{-1}A^\intercal Y $$

椭球拟合

三维坐标变换

现在我们考虑一个空间中的任意椭球,为了方便计算误差函数我们需要对它进行坐标变换,使其中心与坐标原点对齐;假设该变换形式为:先平移一个向量$\vec{s}=-(x_0,y_0,z_0)^\intercal$,后分别绕$z,y,x$轴旋转$\theta_z,\theta_y,\theta_x$。将其写成矩阵形式为:

$$ \vec{x'}=R_xR_yR_z(\vec{x}+\vec{s})\tag{2.1} $$

此时椭球方程可以转换为标准形式:

$$ \frac{x'^2}{a^2}+\frac{y'^2}{b^2}+\frac{z'^2}{c^2}=1 $$

在这种情况下,我们可以简单定义最小距离平方和为:

$$ S_{err}^2=\sum_i\frac{1}{\frac{x_i'^2}{a^2}+\frac{y_i'^2}{b^2}+\frac{z_i'^2}{c^2}} $$

椭球曲面拟合

求解使$S_{err}^2(\theta_x,\theta_y,\theta_z,a,b,c,x_0,y_0,z_0)$最小的拟合曲面解$\vec{s_0}=(\theta_x,\theta_y,\theta_z,a,b,c,x_0,y_0,z_0)^\intercal$是一个非线性问题,无法使用常规线性代数方法解决,不过由于函数值在任意处数值可求,因此在计算机上可以使用高斯-牛顿方法、梯度下降等方法进行求解得到区域极小值;

例如我们可以用一个试解用$\vec{s}$代入误差函数,然后在试解$\vec{s}$附近对误差函数求方向导数$\vec{p}$,然后用新的试解$\vec{p}+\vec{s}$作为下一次尝试,直到两次试解的误差函数相差小于某个特定的值退出计算。

这种方法会遇到两种困难:

1、参数$\vec{s}$和对应坐标的变换表示形式非常复杂,并且目前的数学工具下暂时找不到更加简洁、有效的表达方式;

2、该方法可能会遇到数值困难(在椭球向球体退化的时候),点偏差会造成角度的巨大变化(虽然我们可以把这种情况单独判断,并把对应最小二乘法的Jacobi矩阵秩从9改成6)

我们考虑另一种二次曲面的一般表示形式:

$$ Ax^2+By^2+Cz^2+Dxy+Exz+Fyz+Gx+Hy+Kz+L=0\tag{2.4} $$

该方程有10个元,比之前多了一个元,因此我们需要加一个约束上去,或者消去一个元;

我们将式(2.1)做适当改写:

$$ (R_1R_2R_3)^{-1}\begin{pmatrix} x'\\y'\\z' \end{pmatrix}= \begin{pmatrix} x\\y\\z \end{pmatrix}+ \begin{pmatrix} -x_0\\-y_0\\-z_0 \end{pmatrix} $$

将$(x',y',z')^\intercal$写成极坐标形式:

$$ \begin{pmatrix} x\\y\\z \end{pmatrix}= (R_1R_2R_3)^{-1}\begin{pmatrix} a\cos u\cos v\\b\sin u\cos v\\c\sin v \end{pmatrix}+ \begin{pmatrix} x_0\\y_0\\z_0 \end{pmatrix} $$

则标准椭球方程可以写成:

$$ \begin{aligned} \frac{1}{a^2}&[(x-x_0)\cos\theta_z\cos\theta_y+(y-y_0)(\sin\theta_z\cos\theta_x+\cos\theta_z\sin\theta_y\sin\theta_x)+\\ &(z-z_0)(\sin\theta_z\sin\theta_x-\cos\theta_z\sin\theta_y\cos\theta_z)]^2+\\ \frac{1}{b^2}&[-(x-x_0)\sin\theta_z\cos\theta_y+(y-y_0)(\cos\theta_x\cos\theta_z-\sin\theta_x\sin\theta_y\sin\theta_z)+\\ &(z-z_0)(\cos\theta_z\sin\theta_x+\sin\theta_y\sin\theta_z\cos\theta_x)]^2+\\ \frac{1}{c^2}&[(x-x_0)\sin\theta_y-(y-y_0)\cos\theta_y\sin\theta_x+(z-z_0)\cos\theta_y\cos\theta_x]^2=1 \end{aligned} $$

对比(2.4)式,我们可以得到:

$$ \begin{aligned} A=&\frac{\cos^2\theta_y\cos^2\theta_z}{a^2}+\frac{\cos^2\theta_y\sin^2\theta_z}{b^2}+\frac{\sin^2\theta_y}{c^2}\\ B=&\frac{(\cos\theta_x\sin\theta_z+\sin\theta_x\sin\theta_y\cos\theta_z)^2}{a^2}+ \frac{(\cos\theta_x\cos\theta_z-\sin\theta_x\sin\theta_y\sin\theta_z)^2}{b^2}+ \frac{(\sin\theta_x\cos\theta_y)^2}{c^2}\\ C=&\frac{(\sin\theta_x\sin\theta_y\cos\theta_z-\sin\theta_x\sin\theta_z)^2}{a^2}+ \frac{(\sin\theta_x\cos\theta_z+\cos\theta_x\sin\theta_y\sin\theta_z)^2}{b^2}+ \frac{(\sin\theta_x\cos\theta_y)^2}{c^2}\\ \end{aligned} $$

显然,$A,B,C$是严格正的;相应的我们也可以写出$D,E,...,L$的表达式,不过他们没有什么特殊的性质因此不在此列出。

由于$A,B,C$的性质,我们做如下变形:

$$ \begin{aligned} \frac{1}{3}&[(A+B+C)(x^2+y^2+z^2)+(A-C)(x^2+y^2-2z^2)+(A-B)(x^2-2y^2+z^2)]+\\ &Dxy+Exz+Fyz+Gx+Hy+Kz+L=0 \end{aligned} $$

在上式两边同乘$\frac{3}{A+B+C}$可以得到:

$$ x^2+y^2+z^2-U(x^2+y^2-2z^2)-V(x^2-2y^2+z^2)-4Mxy-2Nxz-2Pyz-Qx-Ry-Sz-T=0 $$

该方程求解虽然也是一个非线性问题,但相比直接求解,不需要每一步都进行那么多次的矩阵运算,极大减小了计算量,提高了计算速度。

遗传算法

既然我们没有线性的方法来求解空间椭球面拟合的方法,那么最后问题基本上都只能转换成数值求解,既然绕不开这个,不如想一些骚一点的办法来解决问题,大佬S同学作为一名技术高超的炼丹师告诉我可以尝试下遗传算法,我一看,感觉这个很靠谱……

具体情况是这样的,从上面的讨论我们可以看出,一般椭球面有九个元,如果我们写出关于这九个元的误差函数,这个函数在空间中的外形是很难确定的,尤其是这九个元当中涉及到三个空间旋转变换,我们所知道的信息几乎就只有连续,以及有下确界,以及在全局最小的解附近函数应该是非常平坦的;与此同时误差函数可能会有很多的区域极小抑或是极度平坦的小区域,使用普适的函数优化方法可能会在这些地方效率低下甚至找到错误的解(即使其对应的误差可能很小,但仍令人不是很满意)。

遗传算法给了我们这么一种看起来在任何情况下都可以比较靠谱的解决问题的方法,形象上的理解是这样的:我们在求解空间里先随机撒一些点作为一个个的个体,这些点的坐标对应着参数,我们可以用某种规则计算出每一个点的“适应度”,比如在椭球拟合问题中,我们可以直接用误差函数来表示“适应度”,在遗传算法中,个体之间采用某种方式交配并产生后代(我们可以理解成使用父母的坐标参数,通过某种算法来计算出后代的坐标参数,使得后代保留父母的部分(坐标)特点),并通过“适应度”来决定让一些很差的个体(比如后代的误差变大了很多)go die;这是一个概率事件;同时我们引入突变函数,使得后代的坐标参数在一个适当的范围内有一个随机的变化,保证群体的演化不至于落入一个区域极小之中。整个计算过程非常类似于大自然中的物种演化过程,因此被称作遗传算法。

遗传算法的优点在于,通过“适应度”和交配,我们可以让物种向我们需要的方向演化,而突变函数则保证,如果有足够丰富和数量的个体和足够长时间的演化,我们几乎最后总能找到全局最优的解。这似乎能打消我们之前的顾虑——我们不用担心椭球的误差函数的实际形状,只需要让遗传算法跑起来就行——但同时我们也知道,任何事情的成功都需要付出代价,遗传算法的代价就是空间和时间;在实际操作中,我们必须小心地调整突变函数的概率和突变程度,即使这样,有时候演化也不尽如人意,在某些极端情况下几乎得不到正确的结果(比如椭球面退化成一个平面),当然这时候之前的方法也没法work了。。

暂时放了个代码在Github上,考虑有空把遗传算法和一般的函数优化问题结合一下,看看能不能提高计算速度和鲁棒性。抽了个空随手把遗传算法部分写了一下,有空再完善了……

Last modification:June 15th, 2019 at 12:37 am
If you think my article is useful to you, please feel free to appreciate

Leave a Comment