【转载】从零开始学习「张氏相机标定法」

【WHO：张氏标定法发明人】

先来简单介绍一下我们的主角：张正友博士。他是世界著名的计算机视觉和多媒体技术的专家，ACM Fellow，IEEE Fellow。现任微软研究院视觉技术组高级研究员。他在立体视觉、三维重建、运动分析、图像配准、摄像机标定等方面都有开创性的贡献。

「张氏标定法」是张正友博士在1999年发表在国际顶级会议ICCV上的论文《Flexible Camera Calibration By Viewing a Plane From Unknown Orientations》中，提出的一种利用平面棋盘格进行相机标定的实用方法。该方法介于摄影标定法和自标定法之间，既克服了摄影标定法需要的高精度三维标定物的缺点，又解决了自标定法鲁棒性差的难题。标定过程仅需使用一个打印出来的棋盘格，并从不同方向拍摄几组图片即可，任何人都可以自己制作标定图案，不仅实用灵活方便，而且精度很高，鲁棒性好。因此很快被全世界广泛采用，极大的促进了三维计算机视觉从实验室走向真实世界的进程。

在介绍「张氏标定法」之前，我们先来搞清楚一些基本问题。

【WHY：为什么要进行相机标定？】

相机标定的目的是：建立相机成像几何模型并矫正透镜畸变。这句话有点拗口，下面分别对其中两个关键部分进行解释。

建立相机成像几何模型：计算机视觉的首要任务就是要通过拍摄到的图像信息获取到物体在真实三维世界里相对应的信息，于是，建立物体从三维世界映射到相机成像平面这一过程中的几何模型就显得尤为重要，而这一过程最关键的部分就是要得到相机的内参和外参（后续文有具体解释）。

矫正透镜畸变：我们最开始接触到的成像方面的知识应该是有关小孔成像的，但是由于这种成像方式只有小孔部分能透过光线就会导致物体的成像亮度很低，于是聪明的人类发明了透镜。虽然亮度问题解决了，但是新的问题又来了：由于透镜的制造工艺，会使成像产生多种形式的畸变，于是为了去除畸变（使成像后的图像与真实世界的景象保持一致），人们计算并利用畸变系数来矫正这种像差。虽然理论上可以设计出不产生畸变的透镜，但其制造工艺相对于球面透镜会复杂很多，所以相对于复杂且高成本的制造工艺，人们更喜欢用数学来解决问题。

【HOW：相机标定的原理】

前面已经说过，相机标定的目的之一是为了建立物体从三维世界到成像平面上各坐标点的对应关系，所以首先我们需要定义这样几个坐标系来为整个过程做好铺垫：

世界坐标系(world coordinate system)：用户定义的三维世界的坐标系，为了描述目标物在真实世界里的位置而被引入。单位为m。

相机坐标系(camera coordinate system)：在相机上建立的坐标系，为了从相机的角度描述物体位置而定义，作为沟通世界坐标系和图像/像素坐标系的中间一环。单位为m。

图像坐标系(image coordinate system)：为了描述成像过程中物体从相机坐标系到图像坐标系的投影透射关系而引入，方便进一步得到像素坐标系下的坐标。单位为m。

像素坐标系(pixel coordinate system)：为了描述物体成像后的像点在数字图像上（相片）的坐标而引入，是我们真正从相机内读取到的信息所在的坐标系。单位为个（像素数目）。

一下子定义出来四个坐标系可能有点晕，下图可以更清晰地表达这四个坐标系之间的关系：

世界坐标系：Xw、Yw、Zw。相机坐标系： Xc、Yc、Zc。图像坐标系：x、y。像素坐标系：u、v。

其中，相机坐标系的轴与光轴重合，且垂直于图像坐标系平面并通过图像坐标系的原点，相机坐标系与图像坐标系之间的距离为焦距f（也即图像坐标系原点与焦点重合）。像素坐标系平面u-v和图像坐标系平面x-y重合，但像素坐标系原点位于图中左上角（之所以这么定义，目的是从存储信息的首地址开始读写）。

在这里我们先引入「棋盘」的概念：

棋盘是一块由黑白方块间隔组成的标定板，我们用它来作为相机标定的标定物（从真实世界映射到数字图像内的对象）。之所以我们用棋盘作为标定物是因为平面棋盘模式更容易处理（相对于复杂的三维物体），但与此同时，二维物体相对于三维物体会缺少一部分信息，于是我们会多次改变棋盘的方位来捕捉图像，以求获得更丰富的坐标信息。如下图所示，是相机在不同方向下拍摄的同一个棋盘图像。

下面将依次对刚体进行一系列变换，使之从世界坐标系进行仿射变换、投影透射，最终得到像素坐标系下的离散图像点，过程中会逐步引入各参数矩阵。

从世界坐标系到相机坐标系

刚体从世界坐标系转换到相机坐标系的过程，可以通过旋转和平移来得到，我们将其变换矩阵由一个旋转矩阵和平移向量组合成的齐次坐标矩阵（为什么要引入齐次坐标可见后续文章）来表示：

其中，R为旋转矩阵，t为平移向量，因为假定在世界坐标系中物点所在平面过世界坐标系原点且与Zw轴垂直（也即棋盘平面与Xw-Yw平面重合，目的在于方便后续计算），所以zw=0，可直接转换成式1的形式。其中变换矩阵

即为前文提到的外参矩阵，之所称之为外参矩阵可以理解为只与相机外部参数有关，且外参矩阵随刚体位置的变化而变化。

下图表示了用R，t将上述世界坐标系转换到相机坐标系的过程。

从相机坐标系到理想图像坐标系（不考虑畸变）

这一过程进行了从三维坐标到二维坐标的转换，也即投影透视过程（用中心投影法将物体投射到投影面上，从而获得的一种较为接近视觉效果的单面投影图，也就是使我们人眼看到景物近大远小的一种成像方式）。我们还是拿针孔成像来说明（除了成像亮度低外，成像效果和透镜成像是一样的，但是光路更简单）。

成像过程如图二所示：针孔面（相机坐标系）在图像平面（图像坐标系）和物点平面（棋盘平面）之间，所成图像为倒立实像。

但是为了在数学上更方便描述，我们将相机坐标系和图像坐标系位置对调，变成图三所示的布置方式（没有实际的物理意义，只是方便计算）：

此时，假设相机坐标系中有一点M，则在理想图像坐标系下（无畸变）的成像点P的坐标为（可由相似三角形原则得出）：

将上式化为齐次坐标表示形式为：

从理想图像坐标系到实际图像坐标系（考虑畸变）

透镜的畸变主要分为径向畸变和切向畸变，还有薄透镜畸变等等，但都没有径向和切向畸变影响显著，所以我们在这里只考虑径向和切向畸变。

径向畸变是由于透镜形状的制造工艺导致。且越向透镜边缘移动径向畸变越严重。下图所示是径向畸变的两种类型：桶形畸变和枕形畸变。

实际情况中我们常用r=0处的泰勒级数展开的前几项来近似描述径向畸变。矫正径向畸变前后的坐标关系为：

由此可知对于径向畸变，我们有3个畸变参数需要求解。

切向畸变是由于透镜和CMOS或者CCD的安装位置误差导致。因此，如果存在切向畸变，一个矩形被投影到成像平面上时，很可能会变成一个梯形。切向畸变需要两个额外的畸变参数来描述，矫正前后的坐标关系为：

由此可知对于切向畸变，我们有2个畸变参数需要求解。

综上，我们一共需要5个畸变参数（k1、k2、k3、p1和p2 ）来描述透镜畸变。

从实际图像坐标系到像素坐标系

由于定义的像素坐标系原点与图像坐标系原点不重合，假设图像坐标系原点在像素坐标系下的坐标为（u0，v0），每个像素点在图像坐标系x轴、y轴方向的尺寸为：dx、dy，且像点在实际图像坐标系下的坐标为（xc，yc），于是可得到像点在像素坐标系下的坐标为：

化为齐次坐标表示形式可得：

公式2中(xp, yp)与公式5中(xc, yc)相同，都是图像坐标系下的坐标。

若暂不考虑透镜畸变，则将式2与式5的转换矩阵相乘即为内参矩阵M：

之所以称之为内参矩阵可以理解为矩阵内各值只与相机内部参数有关，且不随物体位置变化而变化。

最后用一幅图来总结从世界坐标系到像素坐标系（不考虑畸变）的转换关系：

从零开始学习「张氏相机标定法」（二）单应矩阵

前面文章《从零开始学习「张氏相机标定法」（一）成像几何模型》中我们已经得到了像素坐标系和世界坐标系下的坐标映射关系：

其中，u、v表示像素坐标系中的坐标，s表示尺度因子，fx、fy、u0、v0、γ（由于制造误差产生的两个坐标轴偏斜参数，通常很小）表示5个相机内参，R,t表示相机外参，Xw、Yw、Zw（假设标定棋盘位于世界坐标系中Zw=0的平面）表示世界坐标系中的坐标。

单应性概念的引出

我们在这里引入一个新的概念：单应性（Homography）变换。可以简单的理解为它用来描述物体在世界坐标系和像素坐标系之间的位置映射关系。对应的变换矩阵称为单应性矩阵。在上述式子中，单应性矩阵定义为：

其中，M是内参矩阵

从单应矩阵定义式子来看，它同时包含了相机内参和外参。在进一步介绍相机标定知识之前，我们重点来了解一下单应性，这有助于深入理解相机标定。因为在计算机视觉领域，单应性是一个非常重要的概念。

为了不让读者一上来就淹没在公式的汪洋大海中失去兴趣，我们颠倒一下顺序，先来看看单应性到底有什么用，然后再介绍单应矩阵的估计方法。

单应性在计算机视觉中的应用

单应性在计算机视觉领域是一个非常重要的概念，它在图像校正、图像拼接、相机位姿估计、视觉SLAM等领域有非常重要的作用。

图像校正

用单应矩阵进行图像矫正的例子如下图所示，最少需要四个对应点对（后面会给出原因）就可以实现。

视角变换

单应矩阵用于视角变换的例子如下图所示，可以方便地将左边普通视图转换为右图的鸟瞰图。

图像拼接

既然单应矩阵可以进行视角转换，那我们把不同角度拍摄的图像都转换到同样的视角下，就可以实现图像拼接了。如下图所示，通过单应矩阵H可以将image1和image2都变换到同一个平面。

单应矩阵用于图像拼接的例子如下所示。

增强现实（AR）

平面二维标记图案（marker）经常用来做AR展示。根据marker不同视角下的图像可以方便的得到虚拟物体的位置姿态并进行显示，如下图所示。

如何估计单应矩阵？

了解了上述单应性的部分应用后，我们就有很大的动力来学习单应矩阵的推导和计算了。首先，我们假设两张图像中的对应点对齐次坐标为(x',y',1)和(x,y,1)，单应矩阵H定义为：

则有：

矩阵展开后有3个等式，将第3个等式代入前两个等式中可得：

也就是说，一个点对对应两个等式。在此插入一个讨论：单应矩阵H有几个自由度？

或许有人会说，9个啊，H矩阵不是9个参数吗？从h11到h33总共9个。真的是这样吗？实际上并不是，因为这里使用的是齐次坐标系，也就是说可以进行任意尺度的缩放。比如我们把hij乘以任意一个非零常数k并不改变等式结果：

所以实际上单应矩阵H只有8个自由度。8自由度下H计算过程有两种方法。

第一种方法：直接设置 h33=1，那么上述等式变为：

第二种方法：将H添加约束条件，将H矩阵模变为1，如下：

以第2种方法（用第1种也类似）为例继续推导，我们将如下等式（包含||H||=1约束）:

乘以分母展开，得到：

整理，得到：

假如我们得到了两幅图片中对应的N个点对（特征点匹配对），那么可以得到如下线性方程组：

写成矩阵形式：

由于单应矩阵H包含了||H||=1约束，因此根据上图的线性方程组，8自由度的H我们至少需要4对对应的点才能计算出单应矩阵。这也回答了前面图像校正中提到的为何至少需要4个点对的根本原因。

但是，以上只是理论推导，在真实的应用场景中，我们计算的点对中都会包含噪声。比如点的位置偏差几个像素，甚至出现特征点对误匹配的现象，如果只使用4个点对来计算单应矩阵，那会出现很大的误差。因此，为了使得计算更精确，一般都会使用远大于4个点对来计算单应矩阵。另外上述方程组采用直接线性解法通常很难得到最优解，所以实际使用中一般会用其他优化方法，如奇异值分解、Levenberg-Marquarat（LM）算法（后续文章会介绍）等进行求解。

如何根据标定图得到单应矩阵？

经过前面一系列的介绍，我们应该大致明白如何根据打印的棋盘标定图和拍摄的照片来计算单应矩阵H。我们来总结一下大致过程。

1、打印一张棋盘格标定图纸，将其贴在平面物体的表面。

2、拍摄一组不同方向棋盘格的图片，可以通过移动相机来实现，也可以移动标定图片来实现。

3、对于每张拍摄的棋盘图片，检测图片中所有棋盘格的特征点（角点，也就是下图中黑白棋盘交叉点，中间品红色的圆圈内就是一个角点）。我们定义打印的棋盘图纸位于世界坐标系Zw=0的平面上，世界坐标系的原点位于棋盘图纸的固定一角（比如下图中黄色点）。像素坐标系原点位于图片左上角。

4、因为棋盘标定图纸中所有角点的空间坐标是已知的，这些角点对应在拍摄的标定图片中的角点的像素坐标也是已知的，如果我们得到这样的N>=4个匹配点对（越多计算结果越鲁棒），就可以根据LM等优化方法得到其单应矩阵H。当然计算单应矩阵一般不需要自己写函数实现，OpenCV中就有现成的函数可以调用，对应的c++函数是:

Mat findHomography(InputArray srcPoints, InputArray dstPoints, int method=0, double ransacReprojThreshold=3, OutputArray mask=noArray() )

从函数定义来看，只要输入匹配点对，指定具体计算方法即可输出结果。

至此，我们已经搞清楚单应矩阵的概念、推导和应用了，下一篇文章我们继续学习张氏相机标定法。

从零开始学习「张氏相机标定法」（三）推导求解

前面文章《从零开始学习「张氏相机标定法」（一）成像几何模型》中我们学习了相机成像几何模型，知道了如何将世界坐标系中的三维坐标和像素坐标系中的二维坐标联系起来，根据《从零开始学习「张氏相机标定法」（二）单应矩阵》，我们根据标定棋盘图纸及其对应的照片已经可以得到单应矩阵H了。如下所示：

下一步如何求相机内外参数呢？

我们知道H是内参矩阵和外参矩阵的混合体，而我们想要最终分别获得内参和外参。所以需要想个办法，先把内参求出来（先求内参是因为更容易），得到内参后，外参也就随之解出了。

我们先不考虑镜头畸变，来看看如何求解内参和外参。求解思路是利用旋转向量的约束关系，以下是具体推导，建议自己演算一遍，加深理解。

为了利用旋转向量之间的约束关系，我们先将单应性矩阵H化为3个列向量，即H=[h1 h2 h3]，则有

按元素对应关系可得：

因为旋转向量在构造中是相互正交的，即r1和r2相互正交，由此我们就可以利用“正交”的两个含义，得出每个单应矩阵提供的两个约束条件：

约束条件1：旋转向量点积为0（两垂直平面上的旋转向量互相垂直），即：

约束条件2：旋转向量长度相等（旋转不改变尺度），即：

所以一个单应性矩阵H可以提供上述两个约束条件。那么如何利用上述两个约束条件求解内参或者外参呢？我们一步一步来看，由前面可知内参矩阵M：

记：

我们看到B为对称矩阵，真正有用的元素只有6个（主对角线任意一侧的6个元素）。把B带入前面两个约束条件后可转化为：

上面两约束中的式子均可写为通式

的形式，定义3X3的单应矩阵H=[h1 h2 h3]的第i列列向量：

将如下表达式代入上述的约束单项式：

为了简化表达形式，令：

则有：

由此，两约束条件最终可以转化为如下形式：

如果我们拍摄了n张不同角度的标定图片，因为每张图片我们都可以得到一组（2个）上述的等式。其中，v12,v11,v22可以通过前面已经计算好的单应矩阵得到，因此是已知的，而b中6个元素是待求的未知数。因此，至少需要保证图片数 n>=3，我们才能解出b。

根据n张不同角度的标定图片，最终我们得到了一个矩阵集合 Vb=0 ，其中V是一个 (2n x 6) 的矩阵。如果 n>=3，就可以得到唯一解b（带有一个比例因子）。

如果 n=2，也就是只有两张标定图片，那么我们可以设置内参中的γ=0（γ表示由于制造误差产生的两个坐标轴偏斜参数，通常很小，可忽略），将前面式子（搬运到下图）中γ=0可以看到对应 B12=0，换句话说，就是增加了一个约束条件：[0, 1, 0, 0, 0, 0]b = 0。

如果n=1，只能假设u0, v0已知（位于图像中心）且 γ=0，这样只能解出fx, fy两个参数。

前面说到，B中包含一个尺度因子λ，即：

假设我们已经根据前面计算得到了矩阵B元素的值，那么根据已知的矩阵B很容易解出内参，如下：

得到内参数后，内参矩阵M也已知。单应矩阵H也已知，因此可继续求得外参数：

其中又由旋转矩阵性质有

则可得

实际情况下，数据中是存在噪音的，所以计算得到的旋转矩阵R并不一定能满足旋转矩阵的性质。所以通常根据奇异值分解来得到旋转矩阵R。

上述的推导结果是基于理想情况下的解，从理论上证明了张氏标定算法的可行性。但在实际标定过程中，一般使用最大似然估计进行优化。假设我们拍摄了n张标定图片，每张图片里有m个棋盘格角点。三维空间点X在图片上对应的二维像素为x，三维空间点经过相机内参M，外参R，t变换后得到的二维像素为x'，假设噪声是独立同分布的，我们通过最小化x, x'的位置来求解上述最大似然估计问题：

现在我们来考虑透镜畸变的影响，由于径向畸变的影响相对较明显，所以主要考虑径向畸变参数，根据经验，通常只考虑径向畸变的前两个参数k1,k2就可以（增加更多的参数会使得模型变的复杂且不稳定）。实际求解中，通常把k1,k2也作为参数加入上述函数一起进行优化，待优化函数如下所示

上述非线性优化问题通常用Levenberg-Marquardt（LM）算法进行迭代求解。一般将k1,k2初值设为0。

至此，张氏相机标定理论部分已经介绍完毕。下一篇会详细介绍一下LM算法。

从零开始学习「张氏相机标定法」（四）优化算法前传

原创 sixgod 计算机视觉life 2018-03-15

坚持原创，点击上方蓝字关注我！

前一篇说到LM优化算法。在介绍LM优化算法前，得先说说它是怎么一步步发展来的。因此，优化算法部分的介绍顺序是：梯度下降、牛顿法、高斯牛顿法、LM法。

一般待优化的目标函数都是由误差的平方和组成，也就是我们常见的最小二乘问题：

这里的 f 一般是非线性函数。我们在微积分课程里学过，如果想要求目标函数的极值（可能是极大值、极小值或者鞍点处的值），可以通过令目标函数的导数（如果可导）为零，求解析解得到。但是实际应用中函数 f 一般是非常复杂的非线性函数，没有解析解。所以一般最常用的优化方法就是迭代求解。

如何迭代呢？

迭代求解

以迭代求解极小值为例（如果求极大值，在前面加个负号转化为求极小值问题）。通用的方法就是先给定一个初值 θ0，对应目标函数值为 f(θ0)，这个初值可以是经验估计的或者随机指定（随机的话收敛会慢一些），然后我们改变（增大或减少） θ0 的值，得到一个新的值 θ1，如果 f(θ1)<f(θ0)，那么说明我们迭代的方向是朝着目标函数值减小的方向，离我们期待的极小值更近了一步，继续朝着这个方向改变 θ1的值，否则朝着相反的方向改变 θ1的值。如此循环迭代多次，直到达到终止条件结束迭代过程。这个终止条件可以是：相邻几次迭代的目标函数值差别在某个阈值范围内，也可以是达到了最大迭代次数等。我们认为迭代终止时的目标函数值就是极小值。

因此迭代的一般过程如下：

1、对于目标函数 f(x)，给定自变量某个初始值x0。这个初值可以是经验估计的或者随机指定。

2、根据采用的具体方法（梯度下降、牛顿、高斯牛顿、LM等）确定一个增量△xk

3、计算目标函数添加了增量后的值如下

4、如果达到迭代终止条件（达到最大迭代次数或函数值/自变量变化非常小），则迭代结束，可以认为此时对应的目标函数值就是最小值。

5、如果没有达到迭代终止条件，按如下方式更新自变量，并返回第2步。

因此，不同的优化算法的不同主要体现在增量的更新方式上。如果采用不合适的更新方式，其实很容易陷入局部最小值。这里给出一个关于局部最小值和全局最小值的直观的理解。如下图所示：

通常我们希望迭代得到的是全局最小值（只有一个）而不是局部最小值（可能有多个）。下图表示的是不同初始值对迭代结果的影响。如果我们初值设的是红色的θ0，最后找到的极小值位于 θ∞，在这个例子中它是全局最小值。但是如果我们初值设的是蓝色的 θ'0，最后我们找到的极小值位于 θ'∞，显然只是局部最小值，不是全局最小值。

实际上按照上述迭代方法，很多时候我们找到的所谓「全局最小值」其实只是局部最小值，这并不是我们期望的。有没有保证局部最小值一定是全局最小值的情况呢？

答案是有，当然这就需要对目标函数有一定的限制。如果我们的目标函数是凸函数，那么可以保证最终找到的极小值就是最小值。很多人可能忘记了什么是凸函数，下面简单介绍一下。

凸函数的概念

凸函数（convex function）在优化算法中是一个非常重要的概念，其定义如下：

如果上面定义没有看懂也没关系，我们可以发现，凸函数其实就是一个「碗」状曲线，碗口向上。下图列举了几种凸函数和非凸函数的例子。如果一个函数是凸函数，那么它只有一个极小值点，也就是最小值点，从下图中也可以看出来。下图从左到右从上到下分别展示了：只有一个极小值的凸函数、只有一个极小值的非凸函数、有多个极值的非凸函数、加了噪音的非凸函数。

如何判断一个函数是否是凸函数？

对于一元二次可微函数，如果它的二阶导数是正的，那么该函数就是严格凸函数。对于多元二次可微函数，当它的Hessian矩阵（不懂没关系，后面会讲）是正定的时候，就是严格凸函数。好了，凸函数就介绍这么些，对凸函数做个总结吧。

1、如果一个函数是凸函数，那么它有且只有一个极值点，且这个极值点是最值点。

2、几个非负凸函数的和仍然是凸函数。

3、如果一个函数不是凸函数，那么它可能只有一个极值点，也可能有多个极值点（局部极值点和全局极值点）

凸函数讲完了，我们继续说说迭代求解最小值法。从最简单易懂的梯度下降法说起。

梯度下降（Gradient Descent）法

先来假设一个场景：一个人去探险被困在一座山上的某个位置，他必须赶在太阳下山之前回到山底的营地，否则会有危险。但是他对这座山不熟悉，而且视野有限无法看清哪里是山底，时间不多了，因此他必须采取比较激进的最快下山路线。由于可视距离有限，比如说10米。那么他需要以当前位置（坐标x0）为中心，以10米为半径寻找周围最陡峭的路（假设没有障碍），盯着这条10米长的路终点（坐标x1）走，当走到坐标为x1位置时，需要重新观察以当前位置为中心，以可视距离10米为半径区域重新寻找最陡峭的路，盯着这条路的终点（坐标x2）走，如此往复，那么最终他能够走到山底吗？

上述过程就是梯度下降法的实例。梯度下降法是求解无约束优化问题最简单和最古老的方法之一，虽然现在已经不具有实用性，但是许多有效算法都是以它为基础进行改进和修正而得到的。

梯度下降法，其本质是一阶优化算法，有时候也称为最速下降法（Steepest Descent）。与其对应的是梯度上升法(Gradient ascent)，用来求函数的极大值，两种方法原理一样，只是计算的过程中正负号不同而已。本文后续都是以梯度下降法为例进行介绍。

在微积分里，对函数求导数（如果是多元函数则求偏导）得到的就是梯度。具体来说，对于函数 f(x,y)，在点 (x0,y0)，沿着梯度向量的方向就是(∂f/∂x0, ∂f/∂y0)，它的方向是 f(x,y) 增加最快的方向。或者说，沿着梯度向量的方向，更加容易找到函数的极大值。反过来说，沿着梯度向量相反的方向，也就是 -(∂f/∂x0, ∂f/∂y0) 的方向，梯度减少最快，也就是更加容易找到函数的极小值。

以一元函数为例，假设目标函数为F(x)，梯度的符号为∇，那么在梯度下降中，自变量x每次按照如下方式迭代更新：

其中 λ 称为步长或者学习率。按照上述方式迭代更新自变量，直到满足迭代终止条件，此时认为自变量更新到函数最小值点。下图是梯度下降法的示意图。

梯度下降法看起来很好，但是在具体使用过程中，会存在不少问题。

问题1：初始值选取

如下图所示，是二维函数优化的一个例子，不同颜色表示不同高度。我们初始值在蓝色的加号位置，最小值位于绿色的加号位置。图中展示了从不同的初始值开始迭代所需要的迭代次数，可以看到不同初始值会产生非常大的差异。(a) 图只经过5次迭代就收敛了， (b) 图经过了22次迭代才收敛。

那么 (b)图里究竟为什么会多出那么多次的迭代呢？将局部细节放大后如下图所示，这部分「锯齿」状的震荡式下降是因为它经历了一个高度逐渐下降的「山谷」，而每次下降方向都与上一次垂直。换句话说，这种下降方式太「短视」，只关注了一个极小范围内的变化，而忽略了长远的趋势（牛顿法就是根据这点进行的改进，具体见后介绍），并不是一种「平滑」的下降。

既然初始值影响这么大，那么一般采用随机多次选取初始值进行多次迭代，然后将多次迭代结果进行比较，找出其中最小的那个作为最小值。这样会导致计算量急剧变大，最终确定的值也不能保证最小，并且因为随机选取，会导致每次最小值都不同，结果不稳定。

问题2：步长选取

步长，就是每一步走多远，这个参数如果设置的太大，那么很容易就在最小值附近徘徊；相反，如果设置的太小，则会导致收敛速度过慢。所以步长的选取也是一个非常关键的因素。

总结一下梯度下降法：

1、实际使用中使用梯度下降法找到的通常是局部极小值，而不是全局最小值。

2、具体找到的是哪个极小值和初始点位置有很大关系。

3、如果函数是凸函数，那么梯度下降法可以保证收敛到全局最优解（最小值）

4、梯度下降法收敛速度通常比较慢

5、梯度下降法中搜索步长 λ 很关键，步长太长会导致找不到极值点甚至震荡发散，步长太小收敛非常慢。

因此，前面的例子中，迷路的探险者可能花费了大量时间找到的却是一个假的山脚（山腰的某个低洼处，局部最小值）。

从零开始学习「张氏相机标定法」（五）优化算法正传

原创 sixgod 计算机视觉life 2018-03-26

坚持原创，点击上方蓝字关注我！

前一篇《从零开始学习「张氏相机标定法」（四）优化算法前传》花了不少篇幅介绍了最简单的求解无约束优化问题的梯度下降法，这篇我们从算法的演进方向分别介绍牛顿法、高斯牛顿法、Levenberg-Marquardt（LM）法。

牛顿法

介绍牛顿法之前，我们先来看看它到底有什么优势？能够解决什么问题？这样才有兴趣去了解它。

前面文章中介绍了梯度下降法在山谷中容易出现锯齿状的迭代，如下图所示。从不同的初始值开始迭代所需要的迭代次数会产生非常大的差异。

这种锯齿状的迭代，可以通过牛顿法来避免。对于同样的目标函数，使用牛顿法在不同的初始值位置（下图中蓝色十字），只需要3、4次迭代即可找到最小值（下图中绿色十字）。

知道了牛顿法相较于梯度下降法的优势了，那么牛顿法是怎么来的？原理是什么？

首先回顾一下，泰勒展开的通式，如下：

待优化的目标函数和前面一样：

我们将目标函数在x处进行二阶泰勒展开。对照上述通式，目标函数泰勒展开如下：

其中 J(x) 是 ||f(x)||2 关于x的Jacobian矩阵（一阶导数），H(x)则是Hessian矩阵（二阶导数）。

什么是Jacobian矩阵和Hessian矩阵呢？

Jacobian矩阵

假设函数 f = {f1(x1,...,xn), .., fm((x1,...,xn))} 有m维空间组成，对应n个自变量。那么函数 f 的一阶偏导数（如果存在）可以组成一个m行n列的矩阵，称为Jacobian（译为雅克比）矩阵

Hessian矩阵

假设有一个实数函数 f(x1,...,xn)，如果函数 f 的所有二阶偏导数存在，那么称如下矩阵为Hessian矩阵，它是一个n x n 的方阵。

而梯度下降法可以认为是泰勒展开式中只取到了一阶项，二阶项后面的部分丢弃。而在牛顿法中，泰勒展开中取到二阶导数，我们的目标函数变为：

对其关于△x求导，并令其为0，可以得到步长

类似于梯度下降法，我们用上述步长来迭代优化，最终收敛到极值点。这个过程就称为牛顿法。

牛顿法相比于梯度下降法的优点是什么呢? 为什么会有这样的优点？我们来直观的理解一下。

如下图所示。红点和蓝点的梯度（一阶导数）是一样的，但是红点出的二阶导数比蓝点大，也就是说在红点处梯度变化比蓝点处更快，那么最小值可能就在附近，因此步长就应该变小；而蓝点处梯度变化比较缓慢，也就是说这里相对平坦，多走一点也没事，那么可以大踏步往前，因此步长可以变大一些。所以我们看到牛顿法迭代优化时既利用了梯度，又利用了梯度变化的速度（二阶导数）的信息。

从几何上说，牛顿法就是用一个二次曲面去拟合你当前所处位置的局部曲面，而梯度下降法是用一个平面去拟合当前的局部曲面，通常情况下，二次曲面的拟合会比平面更好，所以牛顿法选择的下降路径会更符合真实的最优下降路径。

从本质上去看，牛顿法是二阶收敛，梯度下降是一阶收敛，所以牛顿法收敛更快。比如你想找一条最短的路径走到一个盆地的最底部，梯度下降法每次只从你当前所处位置选一个坡度最大的方向走一步，牛顿法在选择方向时，不仅会考虑坡度是否够大，还会考虑你走了一步之后，坡度是否会变得更大。所以，可以说牛顿法比梯度下降法看得更远一点，能更快地走到最底部。也就是说梯度下降法更贪心，只考虑了眼前，反而容易走出锯齿形状走了弯路；而牛顿法目光更加长远，所以少走弯路。

有人会问了，牛顿法既然那么好，为啥还有后续的高斯牛顿法、LM法？牛顿法肯定有什么缺点吧？

是的，牛顿法确实有比较明显的缺点。前面推导过程，可以看到：牛顿法中包含了Hessian 矩阵的计算。而在高维度下计算Hessian 矩阵需要消耗很大的计算量，很多时候甚至无法计算。

下面介绍高斯牛顿法，看看它是怎么改进牛顿法这个缺点的。

高斯牛顿法

高斯牛顿（Gauss-Newton）法是对牛顿法的一种改进，它用雅克比矩阵的乘积近似代替牛顿法中的二阶Hessian 矩阵，从而省略了求二阶Hessian 矩阵的计算。下面来看看高斯牛顿法是怎么一步步推导来的。

牛顿法中我们是将目标函数 f(x+△x)^2 进行泰勒展开，而在高斯牛顿法中，我们对 f(x+△x) 进行一阶泰勒展开，如下：

前面我们说过，我们要确定一个步长，使得目标函数

达到最小，将一阶泰勒展开代入，问题转化为如下最小二乘问题：

展开后可得：

对△x求导并令导数为0，可得：

则有：

其中我们把J^TJ作为牛顿法中Hessian矩阵的近似，从而用求一阶雅克比矩阵代替了直接求二阶Hessian矩阵的过程。上述就是高斯牛顿法中步长的推导过程。

这个结果完美吗？No!

我们知道牛顿法中，要求Hessian矩阵是可逆的并且正定。但在高斯牛顿法中，用来近似Hessian矩阵的J^TJ可能是奇异矩阵或者病态的，此时会导致稳定性很差，算法不收敛。

另外一个不可忽视的问题是，由于前面我们采用二阶泰勒展开来进行的推导，而泰勒展开只是在一个较小的范围内的近似，因此如果高斯牛顿法计算得到的步长较大的话，上述的近似将不再准确，也会导致算法不收敛。

肿么办？LM算法可以修正上述问题。

LM法

Levenberg-Marquardt（LM）法在一定程度上修正了高斯牛顿法的缺点，因此它比高斯牛顿法更加鲁棒，不过这是以牺牲一定的收敛速度为代价的--它的收敛速度比高斯牛顿法慢。

下面来看看LM算法到底怎么修正高斯牛顿法的缺点的。

LM算法中定义的步长为：

其中I是单位矩阵，u是一个非负数。如果 u 取值较大时，uI 占主要地位，此时的LM算法更接近一阶梯度下降法，说明此时距离最终解还比较远，用一阶近似更合适。反之，如果 u 取值较小时，H 占主要地位，说明此时距离最终解距离较近，用二阶近似模型比较合适，可以避免梯度下降的“震荡”，容易快速收敛到极值点。因此参数 u 不仅影响到迭代的方向还影响到迭代步长的大小。

设x初值为x0，根据经验可以设置u的初值u0为：

其中，τ可以自己指定，aii表示矩阵A对角线元素。

LM采用的搜索方法是信赖域（Trust Region）方法，和梯度下降、牛顿法、高斯牛顿法采用的线性搜索（Line Search）方法不同。

为什么要用信赖域呢？这是因为高斯牛顿法中采用近似二阶泰勒函数只在展开点附近有较好的近似效果，如果步长太大近似就不准确，因此我们应该给步长加个信赖区域，在信赖区域里，我们认为近似是有效的，出了这个区域，近似会出问题。

那么如何确定信赖区域的范围呢？一个比较好的方法是根据我们的近似模型跟实际函数之间的差异来确定。使用如下因子来判断泰勒近似是否足够好：

其中，分子是实际函数迭代下降的值，分母是近似模型下降的值。如果 ρ 接近于1，认为近似比较准确，可以扩大信赖范围；如果 ρ 远小于1，说明实际减小的值和近似减少的值差别很大，也就是说近似比较差，需要缩小信赖范围。

下面给出LM算法优化的一个实例流程。

1、给定初始点 x0 ，计算初始信赖域半径 u0

2、对于第 k 次迭代，根据前面式子求出步长 △xk，并计算得到 ρ

3、若 ρ ≤0.25 ，说明步子迈得太大了，应缩小信赖域半径，令

4、若 ρ ≥0.75 ，说明这一步近似比较准确，可以尝试扩大信赖域半径，令

5、若 0.25<ρ<0.75 ，说明这一步迈出去之后，处于“可信赖”和“不可信赖”之间，可以维持当前的信赖域半径，令

6、若 ρ≤0 ，说明函数值是向着上升而非下降的趋势变化了（与最优化的目标相反），这说明这一步迈得错得“离谱”了，这时不应该走到下一点，而应“原地踏步”，即

并且和上面 ρ≤0.25的情况一样缩小信赖域。反之，在 ρ>0 的情况下，都可以走到下一点，即

上面优化过程中的阈值参数只是作为示例使用的经验值，也可以自己指定。

LM算法可以一定程度避免系数矩阵的非奇异和病态问题，可以提供更鲁棒、更准确的步长。因此LM算法在相机标定、视觉SLAM等领域中应用非常广泛。

LM算法的C/C++实现见：http://users.ics.forth.gr/~lourakis/levmar/

https://mp.weixin.qq.com/mp/appmsgalbum?__biz=MzIxOTczOTM4NA==&action=getalbum&album_id=1362922280875311105&scene=173&from_msgid=2247484931&from_itemidx=1&count=3&nolastread=1#wechat_redirect