角色动画研究 —— 使用雅可比矩阵(Jacobian)来结算IK

这篇博客介绍了如何通过雅可比行列式来进行IK的结算。

---

雅可比矩阵

在IK结算中，如果骨骼的架构比较复杂，那么往往几何分析的结算方法就不太现实了。此时通常使用数值方法让其自己构建出来 —— 在每个时间步长的时候进行对应的计算来得到此时每一个骨节的角度应该进行怎样的改变。在这些数值方法中比较好用的就是使用雅可比矩阵(Jacobian)来进行结算。

从数学意义上来解释雅可比矩阵，我们可以想象有6个函数，每个函数对应着有6个变量。那么针对每个输入变量xi，就会能够得到对应的yi.

y1=f1(x1,x2,x3,x4,x5,x6)

y2=f2(x1,x2,x3,x4,x5,x6)

y3=f3(x1,x2,x3,x4,x5,x6)

y4=f4(x1,x2,x3,x4,x5,x6)

y5=f5(x1,x2,x3,x4,x5,x6)

y6=f6(x1,x2,x3,x4,x5,x6)

因此yi的导数可以被写成：

dyi=∂fi∂x1dx1+∂fi∂x2dx2+∂fi∂x3dx3+∂fi∂x4dx4+∂fi∂x5dx5+∂fi∂x6dx6

因此结合上面的方程，我们可以将上面的方程写为向量的形式：

dY=∂F∂XdX

函数F对于X的偏导数矩阵，就被称为雅可比矩阵(Jacobian)。换句话说，雅可比矩阵就是X的速度到Y的速度的映射：

Y˙=J(X)X˙

在某个时间步长，雅可比矩阵其实也就是针对于xi的函数。在下一个时间步长的时候，X已经改变了，因此雅可比矩阵也进行了改变。

雅可比矩阵与动力学

雅可比矩阵的解释

当雅可比矩阵被用在机械臂&动力学&铰链结构时，此时输入的向量xi就变成了此时各个关节的对应值，而输出向量yi就变成了终端关节的位置与朝向：

Y=[px py pz αx αy αz]T

针对于雅可比矩阵与机械臂的分析，在百度文库上有一篇讲述的很详细的PPT，传送门。

在这种情况下，雅可比矩阵实际上真正要表述的是各个关节的角度变化率θ˙与终端关节的位置朝向Y˙的关系，这里划重点敲黑板。

V=Y˙=J(θ)θ˙

向量V就代表终端骨节的线速度和角速度，而且对应的线速度和角速度都与其当前的位置和朝向和目标位置相关。这些速度都是三维空间的，因此每个向量都有x,y,z的元素。

V=[vx vy vz ωx ωy ωz]

而θ˙指代的是各个关节的速度，在上面的方程中，是未知的。

θ˙=[θ1˙ θ2˙ .... θn˙ ]T

J就是雅可比矩阵，指代的就是在当前的pose下，各个关节的变化导致的终端骨节的变化的矩阵。

J=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∂px∂θ1∂py∂θ1∂pz∂θ1...∂αz∂θ1∂px∂θ2∂py∂θ2∂pz∂θ2...∂αz∂θ2...............∂px∂θn∂py∂θn∂pz∂θn...∂αz∂θn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

针对于n关节的骨骼架构，雅可比矩阵应该是6*n的形式。

因此针对于不同的骨骼架构，雅可比矩阵的形式也不一样。例如针对铰链关节(revolute joint)，终端的朝向变化率基本上就是关节基于回转轴的角速度ω；针对移动关节，终端关节的朝向实际上与关节的回转无关；而针对于旋转关节，终端的线速度实际上要通过旋转轴向量与关节到旋转关节的向量的叉乘来实现。

雅可比矩阵示例

考虑一个如下图的三臂铰链关节，所有的关节都约束在xy平面上，所有的约束轴向量都是(0,0,1)。

在这个示例中，终端关节E被尝试移动到目标位置G，我们不考虑终端的目标的朝向应该是哪里。

由于不考虑目标的朝向应该是哪里，那么V向量很容易可以计算出来：

V=⎡⎣⎢⎢(G−E)x(G−E)y(G−E)z⎤⎦⎥⎥

同样的，根据约束条件，可以很容易的计算出各个角速度：

J=⎡⎣⎢⎢((0,0,1)×E)x   ((0,0,1)×(E−P1))x   ((0,0,1)×(E−P2))x((0,0,1)×E)y   ((0,0,1)×(E−P1))y   ((0,0,1)×(E−P2))y((0,0,1)×E)z   ((0,0,1)×(E−P1))z   ((0,0,1)×(E−P2))z⎤⎦⎥⎥

那么根据V=Jθ˙这个方程，可以求得J的逆，然后再求出θ˙的值，也就能够算出此时各个关节的角速度了……

矩阵的求解

正如之前所说，根据方程需要求得θ˙值的话需要求得J的逆矩阵：

V=Jθ˙

因此：

J−1V=θ˙

但是就如上面所说的，通常的雅可比矩阵不一定是n*n形式的，所以此时的求解需要使用到线性代数中的广义逆（pseudoinverse）。

广义逆

矩阵的广义逆指的是在矩阵A不是n∗n形式的时候，还是有一个矩阵A+可以使得：

A+A=In       (n≤m)

或：

AA+=Im       (m≤n)

而A+的值为：

A+=(ATA)−1AT       (n≤m)

或：

A+=AT(AAT)−1       (m≤n)

具体的广义逆的解释这里就略过了，有兴趣的读者可以去自行查阅资料，传送门。

推导过程

具体的推导过程如下：

V=Jθ˙

假设此时m≤n，那么此时有：

J+V=J+Jθ˙

展开J+=JT(JJT)−1，此时有：

J+V=JT(JJT)−1Jθ˙

将(JJT)−1展开为(JT)−1J−1，会发现右边的公式直接被消掉，也就是说：

JT(JJT)−1V=θ˙

针对于上面的等式，令β=(JJT)−1V，此时针对于β可以使用LU分解，也就是说：

(JJT)β=V

最终得到：

JTβ=θ˙

骨节迭代Integration

正常的Integration方法都可以用于骨骼的迭代，例如Euler Integration或者Verlet Integration等…… 但是需要注意的是，雅可比矩阵在下一个时间步长的时候已经改变了所以需要进行一次重新的计算。这个过程一直持续到终端骨节与目标位置的差值已经到了一个可以接受的误差范围为止。

下图展示了一个臂长分别为15，10，5的机械臂。初始的pose角度分别为{π/8,π/4,π/4}，终端骨节目标点位置为{−20,5}。在21帧中，通过目标位置与终端骨节的线性插值方法来进行IK结算。下图显示了时间分别为0，5，15，20的情况：

当然，需要注意的是有时矩阵依然可能会出现奇异性的状况。有人提出一种方法是使用最小二乘法来进行优化，公式如下：

θ˙=JT(JJT+λ2I)−1V

但是这种方法可能会对迭代的收敛性有影响，因此这个方法到底是不是好也有争论。下图展示了是否使用最小二乘法来进行迭代结算的两种情况：

在这种情况下，目标位置是(−35,5)，因此造成了矩阵的奇异性。因此在这个例子中使用最小二乘法能够有更好的表现。

使用雅可比矩阵的转置而不是逆

在线性代数中求逆的工作量是很大的，因此有一种可选择的方法是使用雅可比矩阵的转置而不是逆矩阵来进行结算：

θ˙=αJTV

这种方法可以免去求逆或广义逆的计算量，但是可能存在的风险是终端骨节可能不够稳定。具体的分析这里略过，可能日后会另开博客讲解。

下图是使用转置矩阵的效果图：

额外收获

在研究雅可比矩阵的结算的过程中，我最大的额外收获应该是广义逆的使用。回想起有一段时间再研究Laplacian Mesh Editing的过程中，似乎也可以使用广义逆来进行非n*n形式矩阵的求逆呢……

<全文完>