三坐标在拟合一条直线或者一个平面的时候,最小二乘法是常见的拟合方法。在我们国家标准或者ISO标准中,对被测要素或理想要素(甚至基准)的拟合,也有用最小二乘法拟合的规范,比如下图:

可是,在用最小二乘法拟合直线或者平面的时候,最小二乘法常见有两种方法:- 普通最小二乘法( OLS,Ordinary Least Squares)
- 总体最小二乘法( TLS,Total Least Squares)
那么三坐标在拟合直线或者平面的时候,采用的是哪种最小二乘法呢?今天我们就来讨论一下这个问题,并详细介绍一下总体最小二乘法(TLS)的计算方法。今天的文章,数学味道比较重,有不少线性代数的内容,实在没有耐心看的小伙伴,在这里我直接给出结论: 三坐标采用的是总体最小二乘法(TLS) , 并亲测有效。我已经做成EXCEL表格(仅针对直线),有兴趣的小伙伴可以直接跳到最后,联系助理老师Judy老师(本文最后有她微信联系方式),她会发您TLS和OLS的EXCEL表格。
现在,我们已经采集了一堆离散的点(P1-P5),见下图,我们要将这些点,采用最小二乘法拟合成一条直线(y=kx+b)。图3 一堆离散的点
对最小二乘法熟悉的小伙伴都应该知道,最小二乘法的核心思想是“残差的平方和最小”,可是问题是,哪个残差的平方和最小呢,或者说哪个方向的残差平方和最小呢?这是问题的关键。OLS(普通最小二乘法),在它的眼里,上面的5个离散的点Pi(xi, yi)中,自变量xi是值得信任的,也就是的xi坐标值是精确的或者可控的,因变量yi的坐标是存在测量误差(或者观察误差)的,不值得信任的。而在TLS(总体最小二乘法)看来,上面5个点Pi(xi, yi)中,xi和yi都是有测量或者观察误差的,都不值得信任。OLS既然认为y坐标不值得信任,那么它就要权衡,尽量找到一条直线 y=kx+b,使得y方向的的残差di(或者说估计偏差)的平方和最小,见下图:见上图中红色的d1-d5, 残差的平方和D, 即为:D=d12 + d22 + d32 + d42 +d52所以OLS的目的是, 寻求合适的k和b, 使得上面的D值最小化 。见下面动画:上面的动画中,只要人们能够寻求一个合适的k和b, 使得D值最小,这个时候,这条OLS最小二乘直线,就被拟合出来了。关于OLS的公式和推导方法,德辉学堂往期的公众号文章已经详细讨论过这个话题,有兴趣的小伙伴,可以点击文章最后的链接,这里就不再赘述。因为在TLS看来,观察点P的x值和y值都不可信,这个时候,它的残差就不再是单纯的y方向的偏差,同时也考虑到x方向的偏差。所以对于TLS来说,我们把观察点P到估计直线的 垂直距离d 定义为残差。见下图:D = d12 + d22 + d32 + d42 +d52TLS的目的,就是在寻找合适的k和b, 使得D值最小。见下面动画:如动画2中显示,只要我们能求出合适的k和b, 使得残差平方和D值最小,这样对应的k和b,就能确定这条直线了。稍微总结一下,OLS和TLS的目的,都是为了使得残差d的平方和最小, 只不过两者残差的方向不一样。以直线为例,OLS的残差是y轴方向,而TLS是垂直于直线方向。同样的原理,可以扩展到平面,超平面或者更高维度,这里就不再啰嗦。OLS的计算公式在我之前的公众号文章中讨论过,这里不再讨论。那么TLS的公式是什么呢?我们接下来就来推导一下TLS的公式(有恐数学症的小伙伴可以直接跳过)。假设有n个离散的点Pi(xi, yi),也就是我们用三坐标,在零件表面上测量出来实际的坐标点, 现在我们准备采用TLS(总体最小二乘法)把它们拟合成一条直线。为了方便推导,我不采用斜距式来表达直线,用一般式来表达直线:上面的等式(1)中,因为等式两边可以任意除以一个不为0的常数,所以,我们还可以做一个限制,假设:为啥这样假设?因为这样做,在计算点到直线距离的时候,可以把分母去掉,有利于下一步的推导(同时该限制也是避免残差的平方和趋向无穷小)。根据点到直线的距离公式,这样点Pi(xi, yi)到公式(1)中直线的距离di(即残差)为:再对(3)式求平方得(平方后,绝对值符号可以去掉):大家不要被公式(5)吓到,那个侧翻的M,它就是一个求和公式而已。注意上面的公式中, xi, yi都是已知的测量坐标值,a,b,c是未知的。我把公式(5)再转化一下:现在TLS的任务,就是在函数Q(a,b,c)中,求合适的a,b,c,使得Q(a,b,c)或者D值最小(还记得那个“平方和最小”的宗旨吗?)。接下来用到大学的数学知识,先想办法消去C。我们先对C求偏导,并令其为0(因为极值会发生在这个位置)。注意,上面的公式不管怎么折腾,大家记住,下标是i的量都是已知的(它们是观察数据),比如ui, vi都是已知值,x和y的平均值也是已知的,他们是由实际的坐标值倒腾过来的,a和b才是我们真正关心的(它们是一个未知量)。我们的问题就转化为,求的合适的a,b,使下面公式中D最小:接下来引入矩阵,因为我们最终要用到矩阵的二次型来解决这个极值问题。上式中 wT 中的T表示矩阵的转置。因为 wT *zi 是一个常数,不是一个矩阵,可以任意转置,所以有:wTzi=(wTzi)T =ziT *w (13)对于求和公式(14)来说,w或者wT 都是常数,而且和观察点i无关,可以单独提到求和公式外边,于是将公式(14)再整理可得:在公式15中,括号里的那一坨是矩阵,先整理括号里的一坨:上面这一坨矩阵看起来吓人,其实很简单,他里边都是常数,无非就是一个2x2的一个对称矩阵了,没啥吓人的。再说一遍,等式(17)中的矩阵,我们叫系数矩阵,是由常数构成,是由实测的坐标值倒腾出来常数矩阵,里边没有可怕的变量。我们再来观察一下公式(18), 其中的a,b是未知数,当然有个约束前提(a2+b2=1),中间的矩阵是一个系数矩阵,而且看这个架势,这不是大学线性代数最后一章,那个特别莫名其妙的二次型矩阵么?到线性代数最后一章的时候,矩阵的特征值,特征向量,矩阵的二次型,合同矩阵之类的概念,把人搞得云里雾里。事实上这些概念,在实际工作中非常有用,尤其是在算极值的时候非常有用。公式(18)中,我们从右往左看,中间的矩阵和向量(a, b)相乘,意味着把向量(a,b)经过线性 变换成一个新向量 (这个新向量长度和方向都有可能会发生改变),然后这个新向量再和原向量(a, b)进行内积。然后,我们的问题是,在模长为1(因为a,b的平方和开根号为1)的所有向量中,哪个倒霉的向量,经过公式(18)的转换,结果最小?稍微直观的解释一下,我把公式(18)整理一下, 用(x,y) 代替(a,b),用z代替D, 注意,这里的z就是我们前面讲的平方和这个值,现在我们要看它的最小值在哪里:公式(19)实际上是一个椭圆的抛物面(有兴趣的小伙伴,展开就可以发现),见下图红色抛物面,我们现在在 x2+y2=1的约束条件下,寻求z的最小值。几何上解释,我们用一个圆筒: x2 +y2 =1的一个圆柱(见下图蓝色的圆柱),和这个椭圆抛物面去交,交出来的最低点L(见下图黄点),这个时候最低点L对应的坐标(x,y)就是我们寻求的,满足条件(x2 +y2 =1)下,残差的平方和最小的点。图6中,我们可以发现最低点发生的位置是L点,那么L点发生在哪个位置呢?上图中红色的是一个标准的椭圆抛物面,凭直觉,我们可以发现,这个L点肯定发生在椭圆的大径方向上。我们垂直于z轴做一系列平面,和红色的抛物面进行相交,得到一系列椭圆,我们再来观察这个最低点L发生的位置:见图7,如果我们能找到L点处,得到对应的坐标值(a,b), 然后再计算出c, 我们想要的TLS直线就被拟合出来了。这里的关键点就在于公式(19)中的系数矩阵!我再把公式写在下面:上面A由Sxx, Syy和Sxy构成的系数矩阵,它有一系列的特点。因为篇幅原因,不能给出推导过程,这里只能说结论:![]()
![]()
![]()
![]()
![]()
![]()
- 公式(19)中取得最小值时(约束条件:x2+y2=1),(x,y)一定在系数矩阵的最小特征值对应的特征向量方向上,也就是说向量(a, b), 一定是最小特征值对应的特征向量集中的一个,它一定在椭圆的大径方向上。见图(8)现在的关键点在于,要计算出系数矩阵A的特征值和特征向量。再回到大学时代吧,先令特征值为λ,则矩阵A的特征多项式为:好了,再来求特征向量 p ,我们假设特征向量p为(a1, b1)T , 因为根据矩阵A的特征值和特征向量λ的关系,有以下公式:设 b1=1(因为特征向量是一个集合,不是唯一的,其中元素可以随便设,当然不是0,只要其它元素等比就行),则:然后将 (a1,b1) 归一化(就是单位化,就是使得a2 +b2 =1):好了,直线方程ax+by+c=0中的a,b,c就这样被痛苦的给弄出来了,终于,欧了。整个计算过程比较繁琐,容易迷失,我们再快速梳理一下整个主要计算过程:
整过计算过程,我已经弄成EXCEL表格,只要输入原始数据,自动会生成直线。如果对推导过程不感兴趣的小伙伴,可以联系助理老师索要Excel表格。我用EXCEL表格,将一堆离散的点拟合成直线,然后再计算每个点到该直线的距离,三坐标做同样的操作(因为三坐标里找不到直线方程的表达式),我比较了海克斯康, 蔡司,OGP,国产POWERDMIS等软件的计算结果,结果一模一样。所以我最终的结论是,关于最小二乘法的拟合,三坐标采用的是TLS,也就是总体最小二乘法。或许有小伙伴会问,我怎么知道应该采用哪一种最小二乘法呢?这里涉及到最小二乘法OLS和TLS两种方法的应用场景。简单总结一下就是:- 若自变量是人为设定 / 精确记录 的可控变量(如预算、温度、配比),或者自变量x的误差远远小于因变量y的误差,优先用 OLS。
- 若自变量与因变量 均含显著测量误差 (如坐标测量、坐标转换),则切换至 TLS。
举个例子,比如我们知道一家公司,在某个平台每个月投放的广告金额x和销售额y的数据,我想知道广告金额x和销售额y之间线性关系,因为自变量,投入的广告金额x是精确可控的,而销售额y受各种因素的影响,误差会比较大,这种拟合计算就应该采用OLS。而三坐标在采点的时候,因为x和y方向(或者z方向)的坐标值都会有误差,所以在拟合的时候,应该采用TLS,如同本案例所示。好了,本期文章的内容就到这里了,希望对您有所启发。本期文章讲了3节内容,第一节内容介绍了两种最小二乘法OLS和TLS的拟合原理;第二节内容详细介绍了TLS的具体计算和公式推导;第三节内容将TLS的计算结果和各大三坐标品牌软件的计算进行了比对,并介绍了OLS和TLS的各自的适用情况。