当前位置: 首页 > news >正文

利用 IMU 估计人体关节轴向和位置 —— 论文推导

Title: 利用 IMU 估计人体关节轴向和位置 —— “Joint axis and position estimation from inertial measurement data by exploiting kinematic constraints” —— 论文推导


文章目录

  • I. 论文回顾
  • II. 铰接关节的约束
    • 1. 铰接关节约束的原理
    • 2. 铰接关节约束的梯度
    • 3. 铰接关节约束梯度的人工推导
  • III. 球窝关节的约束
    • 1. 球窝关节约束的原理
    • 2. 球窝关节约束的人工解读
    • 3. 球窝关节约束梯度的符号计算
  • IV. 算法实现
  • V. 小结


本篇博客中, 如在引用段落中的是 AI 生成的内容. 此处主要用的是豆包. 因提示 “服务器繁忙,请稍后再试”, 没有用 Deepseek.

在正常段落中的是由博主真人输出的.


I. 论文回顾

论文 Joint axis and position estimation from inertial measurement data by exploiting kinematic constraints 由 Thomas Seel 等人撰写。文章提出利用关节运动学约束从惯性测量数据中估计关节轴和关节位置的新方法,经仿真和实验验证有效,可应用于机器人、车辆等领域, 尤其是外骨骼控制、人体运动跟踪等场景。

研究背景:惯性测量单元(IMU)在人体运动分析等领域应用广泛,但使用中存在两个问题:传感器偏差导致位置和角度估计漂移;需知道传感器坐标系相对关节轴或安装段的方向信息。现有解决方法存在局限性,因此需要能准确估计关节位置和轴的方法。

利用运动学约束

  • 铰链关节约束:两个通过铰链关节连接的刚性段,其陀螺仪测量的角速度在关节平面投影长度相等,利用该特性可在传感器方向未知时识别铰链关节轴。通过计算相关梯度,使用高斯 - 牛顿算法可搜索满足条件的关节轴向量。此外,还需判断估计的关节轴向量是否对齐。
  • 球窝关节约束:对于球窝关节,结合加速度计读数,可利用关节中心加速度相等的约束来估计关节中心到传感器坐标系原点的偏移向量,同样通过计算梯度并使用高斯 - 牛顿算法进行估计。

建模与仿真:开发三连杆运动学仿真模型,模拟传感器测量。对铰链关节轴和球窝关节偏移向量估计算法进行实现,采用特定参数化和更新过程。仿真结果表明,算法在不同运动、采样率和噪声幅度下表现良好,估计值能快速收敛到真实值附近。

实验结果:将算法应用于基于 IMU 的步态分析实验数据,实验中使用无线运动追踪器,让受试者进行腿部和脚部环绕运动及行走,记录传感器数据。结果显示,算法能在不到二十次迭代中正确识别膝关节轴和踝关节位置,虽存在一定误差,但最小二乘法能有效处理这些不准确性。

结论与展望:提出的最小二乘法能有效估计关节轴和位置,无需传感器安装知识和积分运算,对采样率要求低。该算法可应用于机器人操纵器、链接车辆等领域,未来将关注角度相关关节轴和位置的识别等扩展研究。


II. 铰接关节的约束

1. 铰接关节约束的原理

∥ g 1 ( t ) × j 1 ∥ 2 − ∥ g 2 ( t ) × j 2 ∥ 2 = 0 (1) \left\| g_{1}(t) × j_{1}\right\| _{2}-\left\| g_{2}(t) × j_{2}\right\| _{2}=0 \tag{1} g1(t)×j12g2(t)×j22=0(1)

两个通过铰链关节连接的刚性段,其陀螺仪测量的角速度在关节平面投影长度相等,利用该特性可在传感器方向未知时识别铰链关节轴,其原理涉及几何关系和运动学知识:

  1. 铰链关节的运动特性:铰链关节只允许两个刚性段在一个特定平面内相对转动,这个平面就是关节平面,而关节轴垂直于该平面。
    两个刚性段虽可在空间自由旋转和移动,但相对运动受限于这个铰链关节的约束。
  2. 角速度的几何关系:分别安装在两个刚性段上的陀螺仪,测量得到各自的角速度 g 1 ( t ) g_{1}(t) g1(t) g 2 ( t ) g_{2}(t) g2(t)
    从几何角度看, g 1 ( t ) g_{1}(t) g1(t) g 2 ( t ) g_{2}(t) g2(t) 的差异仅在于关节角的变化速度以及随时间变化的旋转矩阵。
    因为它们围绕同一个铰链关节轴运动,所以在垂直于关节轴的平面(即关节平面)上的投影具有特殊性质。
  3. 投影长度相等的证明:设单位关节轴向量在第一段陀螺仪坐标系中为 j 1 j_{1} j1,在第二段陀螺仪坐标系中为 j 2 j_{2} j2
    根据向量叉乘的几何意义,向量 g i ( t ) g_{i}(t) gi(t) j i j_{i} ji 叉乘的结果,其模长 ∥ g i ( t ) × j i ∥ 2 \left\| g_{i}(t) × j_{i}\right\| _{2} gi(t)×ji2 表示 g i ( t ) g_{i}(t) gi(t) 在垂直于 j i j_{i} ji 方向(即关节平面)上的投影长度。
    由于两个刚性段围绕同一关节轴运动,在任意时刻 t t t g 1 ( t ) g_{1}(t) g1(t) g 2 ( t ) g_{2}(t) g2(t) 在关节平面的投影长度必然相等,即 ∥ g 1 ( t ) × j 1 ∥ 2 − ∥ g 2 ( t ) × j 2 ∥ 2 = 0 \left\| g_{1}(t) × j_{1}\right\| _{2}-\left\| g_{2}(t) × j_{2}\right\| _{2}=0 g1(t)×j12g2(t)×j22=0
  4. 识别关节轴的方法:当传感器方向未知时,在大量的陀螺仪测量数据中,搜索满足 ∥ g 1 ( t ) × j 1 ∥ 2 − ∥ g 2 ( t ) × j 2 ∥ 2 = 0 \left\| g_{1}(t) × j_{1}\right\| _{2}-\left\| g_{2}(t) × j_{2}\right\| _{2}=0 g1(t)×j12g2(t)×j22=0 的向量 j 1 ^ \hat{j_{1}} j1^ j 2 ^ \hat{j_{2}} j2^ ,就能找到可能的关节轴方向。
    通过计算相关梯度(如 d ( ∥ g i ( t ) × j i ∥ 2 ) d j i = ( g i ( t ) × j i ) × g i ( t ) ∥ g i ( t ) × j i ∥ 2 \frac{d\left(\left\| g_{i}(t) × j_{i}\right\| _{2}\right)}{d j_{i}}=\frac{\left(g_{i}(t) × j_{i}\right) × g_{i}(t)}{\left\| g_{i}(t) × j_{i}\right\| _{2}} djid(gi(t)×ji2)=gi(t)×ji2(gi(t)×ji)×gi(t) i = 1 , 2 i = 1, 2 i=1,2),并利用高斯-牛顿算法等优化方法,可在最小二乘意义下找到最符合条件的 j 1 ^ \hat{j_{1}} j1^ j 2 ^ \hat{j_{2}} j2^ ,从而实现对铰链关节轴的识别。

在推导方程 (1) 时,向量 g i ( t ) g_{i}(t) gi(t) j i j_{i} ji 叉乘的模长公式是基于向量叉乘的几何意义和性质得到的。

  1. 向量叉乘的几何意义: 向量 a ⃗ \vec{a} a b ⃗ \vec{b} b 叉乘 a ⃗ × b ⃗ \vec{a}\times\vec{b} a ×b ,其结果是一个向量,该向量的模长 ∣ a ⃗ × b ⃗ ∣ \vert\vec{a}\times\vec{b}\vert a ×b 等于 ∣ a ⃗ ∣ ∣ b ⃗ ∣ sin ⁡ θ \vert\vec{a}\vert\vert\vec{b}\vert\sin\theta a ∣∣b sinθ,其中 θ \theta θ a ⃗ \vec{a} a b ⃗ \vec{b} b 之间的夹角。
    在当前问题中,对于向量 g i ( t ) g_{i}(t) gi(t) 与单位关节轴向量 j i j_{i} ji ∣ j i ∣ = 1 \vert j_{i}\vert = 1 ji=1 ),它们叉乘的模长 ∥ g i ( t ) × j i ∥ 2 \left\| g_{i}(t) × j_{i}\right\| _{2} gi(t)×ji2,根据上述几何意义,就等于 ∣ g i ( t ) ∣ ∣ j i ∣ sin ⁡ θ = ∣ g i ( t ) ∣ sin ⁡ θ \vert g_{i}(t)\vert\vert j_{i}\vert\sin\theta = \vert g_{i}(t)\vert\sin\theta gi(t)∣∣jisinθ=gi(t)sinθ,这里 θ \theta θ g i ( t ) g_{i}(t) gi(t) j i j_{i} ji 的夹角。
  2. 投影长度的表示: 由于 j i j_{i} ji 垂直于关节平面, g i ( t ) g_{i}(t) gi(t) j i j_{i} ji 叉乘的模长 ∥ g i ( t ) × j i ∥ 2 \left\| g_{i}(t) × j_{i}\right\| _{2} gi(t)×ji2,实际上就是 g i ( t ) g_{i}(t) gi(t) 在垂直于 j i j_{i} ji 方向(即关节平面)上的投影长度。
    从几何直观上理解,在以 j i j_{i} ji 为法向量的关节平面上, g i ( t ) g_{i}(t) gi(t) 在该平面的投影长度可以通过上述叉乘模长公式来表示,这是由向量运算的几何性质所决定的。

图形表示:


2. 铰接关节约束的梯度

以下是对梯度公式 d ( ∥ g i ( t ) × j i ∥ 2 ) d j i = ( g i ( t ) × j i ) × g i ( t ) ∥ g i ( t ) × j i ∥ 2 \frac{d\left(\left\| g_{i}(t) × j_{i}\right\| _{2}\right)}{d j_{i}}=\frac{\left(g_{i}(t) × j_{i}\right) × g_{i}(t)}{\left\| g_{i}(t) × j_{i}\right\| _{2}} djid(gi(t)×ji2)=gi(t)×ji2(gi(t)×ji)×gi(t) 的推导过程:

  1. 明确符号与基本定义
    g i ( t ) g_{i}(t) gi(t) j i j_{i} ji 为三维向量, × \times ×表示向量的叉乘运算, ∥ ⋅ ∥ 2 \left\| \cdot \right\| _{2} 2表示向量的2 - 范数(即欧几里得范数)。对于向量 a ⃗ = ( a 1 , a 2 , a 3 ) \vec{a}=(a_1,a_2,a_3) a =(a1,a2,a3),其 2 - 范数 ∥ a ⃗ ∥ 2 = a 1 2 + a 2 2 + a 3 2 \left\|\vec{a}\right\| _{2}=\sqrt{a_1^2 + a_2^2 + a_3^2} a 2=a12+a22+a32

  2. 对向量叉乘的模长函数求导
    h ( j i ) = ∥ g i ( t ) × j i ∥ 2 h(j_{i})=\left\| g_{i}(t) × j_{i}\right\| _{2} h(ji)=gi(t)×ji2,我们要计算 d h ( j i ) d j i \frac{dh(j_{i})}{d j_{i}} djidh(ji)
    根据向量叉乘的定义,若 g i ( t ) = ( g i 1 ( t ) , g i 2 ( t ) , g i 3 ( t ) ) g_{i}(t)=(g_{i1}(t),g_{i2}(t),g_{i3}(t)) gi(t)=(gi1(t),gi2(t),gi3(t)) j i = ( j i 1 , j i 2 , j i 3 ) j_{i}=(j_{i1},j_{i2},j_{i3}) ji=(ji1,ji2,ji3),则 g i ( t ) × j i = ( g i 2 ( t ) j i 3 − g i 3 ( t ) j i 2 , g i 3 ( t ) j i 1 − g i 1 ( t ) j i 3 , g i 1 ( t ) j i 2 − g i 2 ( t ) j i 1 ) g_{i}(t) × j_{i}=(g_{i2}(t)j_{i3}-g_{i3}(t)j_{i2},g_{i3}(t)j_{i1}-g_{i1}(t)j_{i3},g_{i1}(t)j_{i2}-g_{i2}(t)j_{i1}) gi(t)×ji=(gi2(t)ji3gi3(t)ji2,gi3(t)ji1gi1(t)ji3,gi1(t)ji2gi2(t)ji1)
    u = g i ( t ) × j i u = g_{i}(t) × j_{i} u=gi(t)×ji,那么 h ( j i ) = ∥ u ∥ 2 = u 1 2 + u 2 2 + u 3 2 h(j_{i})=\left\|u\right\| _{2}=\sqrt{u_1^2 + u_2^2 + u_3^2} h(ji)=u2=u12+u22+u32
    根据复合函数求导法则,我们先对 h ( j i ) h(j_{i}) h(ji) 关于 u u u 求导,再乘以 u u u 关于 j i j_{i} ji 的导数。

  • h ( j i ) = u 1 2 + u 2 2 + u 3 2 h(j_{i})=\sqrt{u_1^2 + u_2^2 + u_3^2} h(ji)=u12+u22+u32 关于 u u u 求导:
    根据求导公式 ( x ) ′ = 1 2 x (\sqrt{x})^\prime=\frac{1}{2\sqrt{x}} (x )=2x 1,对 h ( j i ) h(j_{i}) h(ji) 关于 u k u_k uk k = 1 , 2 , 3 k = 1,2,3 k=1,2,3)求偏导数可得:
    ∂ h ( j i ) ∂ u k = u k ∥ u ∥ 2 \frac{\partial h(j_{i})}{\partial u_k}=\frac{u_k}{\left\|u\right\| _{2}} ukh(ji)=u2uk,所以 ∂ h ( j i ) ∂ u = u ∥ u ∥ 2 \frac{\partial h(j_{i})}{\partial u}=\frac{u}{\left\|u\right\| _{2}} uh(ji)=u2u(这里 u u u 是向量形式)。
  • u = g i ( t ) × j i u = g_{i}(t) × j_{i} u=gi(t)×ji 关于 j i j_{i} ji 求导:
    我们知道 u k u_k uk 是关于 j i 1 , j i 2 , j i 3 j_{i1},j_{i2},j_{i3} ji1,ji2,ji3 的线性函数。以 u 1 = g i 2 ( t ) j i 3 − g i 3 ( t ) j i 2 u_1 = g_{i2}(t)j_{i3}-g_{i3}(t)j_{i2} u1=gi2(t)ji3gi3(t)ji2 为例,对 u 1 u_1 u1 关于 j i 1 j_{i1} ji1 求偏导数为 0 0 0,对 u 1 u_1 u1 关于 j i 2 j_{i2} ji2 求偏导数为 − g i 3 ( t ) -g_{i3}(t) gi3(t),对 u 1 u_1 u1关于 j i 3 j_{i3} ji3 求偏导数为 g i 2 ( t ) g_{i2}(t) gi2(t)
    同理可得 u 2 u_2 u2 u 3 u_3 u3 关于 j i 1 , j i 2 , j i 3 j_{i1},j_{i2},j_{i3} ji1,ji2,ji3 的偏导数。
    经过整理可以发现, ∂ u ∂ j i = g i ( t ) × \frac{\partial u}{\partial j_{i}}=g_{i}(t)\times jiu=gi(t)×(这里“ × \times ×”表示一种与叉乘相关的矩阵 - 向量运算结构,从结果上看等同于 g i ( t ) g_{i}(t) gi(t) 再次与 u u u 叉乘)。
  1. 应用链式法则
    根据链式法则 d h ( j i ) d j i = ∂ h ( j i ) ∂ u ⋅ ∂ u ∂ j i \frac{dh(j_{i})}{d j_{i}}=\frac{\partial h(j_{i})}{\partial u}\cdot\frac{\partial u}{\partial j_{i}} djidh(ji)=uh(ji)jiu,将前面求得的结果代入:
    d ( ∥ g i ( t ) × j i ∥ 2 ) d j i = u ∥ u ∥ 2 ⋅ ( g i ( t ) × ) = ( g i ( t ) × j i ) × g i ( t ) ∥ g i ( t ) × j i ∥ 2 \frac{d\left(\left\| g_{i}(t) × j_{i}\right\| _{2}\right)}{d j_{i}}=\frac{u}{\left\|u\right\| _{2}}\cdot(g_{i}(t)\times)=\frac{\left(g_{i}(t) × j_{i}\right) × g_{i}(t)}{\left\| g_{i}(t) × j_{i}\right\| _{2}} djid(gi(t)×ji2)=u2u(gi(t)×)=gi(t)×ji2(gi(t)×ji)×gi(t)

综上,得到了给定的梯度计算公式。


3. 铰接关节约束梯度的人工推导

曾经人工推导过如下. (AI 时代, 人工涂鸦看起来可能更温暖一点, 我们何去何从?)

[notes]Joint-1 [notes]Joint-2-page2

III. 球窝关节的约束

1. 球窝关节约束的原理

  1. 场景与前提
    考虑由球窝关节连接的两个连杆,如原文中图 2 所示。球窝关节具有三个自由度,这意味着两个相连刚体可在多个方向上相对转动。由于这种多自由度的特性,两个传感器测量到的角速度之间不存在普遍的固定关系。为了利用运动学约束进行相关研究,需要引入加速度计的读数。

  2. 变量定义
    设两个传感器(分别位于两个连杆上)的加速度分别为 a 1 ( t ) a_{1}(t) a1(t) a 2 ( t ) a_{2}(t) a2(t) o 1 o_{1} o1 o 2 o_{2} o2 分别是从关节中心到第一和第二传感器坐标系原点的偏移向量。
    定义 Γ g i ( o i ) \Gamma_{g_{i}}(o_{i}) Γgi(oi) Γ g i ( o i ) = g i ( t ) × ( g i ( t ) × o i ) + g ˙ i ( t ) × o i \Gamma_{g_{i}}(o_{i}) = g_{i}(t)\times(g_{i}(t)\times o_{i})+\dot{g}_{i}(t)\times o_{i} Γgi(oi)=gi(t)×(gi(t)×oi)+g˙i(t)×oi i = 1 , 2 i = 1,2 i=1,2
    Γ g i ( o i ) \Gamma_{g_{i}}(o_{i}) Γgi(oi) 表示由于绕关节中心旋转而产生的(径向和切向)加速度。
    它综合考虑了角速度 g i ( t ) g_{i}(t) gi(t) 及其变化率 g ˙ i ( t ) \dot{g}_{i}(t) g˙i(t) 与偏移向量 o i o_{i} oi 的关系。
    通过向量叉乘运算,体现了旋转运动对加速度的贡献。

  3. 方程(3)解释
    方程 ∥ a 1 ( t ) − Γ g 1 ( o 1 ) ∥ 2 − ∥ a 2 ( t ) − Γ g 2 ( o 2 ) ∥ 2 = 0 ∀ t \left\| a_{1}(t)-\Gamma_{g_{1}}(o_{1})\right\| _{2}-\left\| a_{2}(t)-\Gamma_{g_{2}}(o_{2})\right\| _{2}=0 \ \forall t a1(t)Γg1(o1)2a2(t)Γg2(o2)2=0 t ,其含义是 ( a 1 ( t ) − Γ g 1 ( o 1 ) ) (a_{1}(t)-\Gamma_{g_{1}}(o_{1})) (a1(t)Γg1(o1)) 表示在第一个局部坐标系下关节中心的加速度, ( a 2 ( t ) − Γ g 2 ( o 2 ) ) (a_{2}(t)-\Gamma_{g_{2}}(o_{2})) (a2(t)Γg2(o2))表示在第二个局部坐标系下关节中心的加速度。
    因为球窝关节连接的两个连杆在关节中心处的加速度,在不考虑坐标系旋转差异时应相等(这里忽略了旋转矩阵的影响,只考虑加速度的大小关系),所以两者加速度的欧几里得范数(即 ∥ a 1 ( t ) − Γ g 1 ( o 1 ) ∥ 2 \left\| a_{1}(t)-\Gamma_{g_{1}}(o_{1})\right\| _{2} a1(t)Γg1(o1)2 ∥ a 2 ( t ) − Γ g 2 ( o 2 ) ∥ 2 \left\| a_{2}(t)-\Gamma_{g_{2}}(o_{2})\right\| _{2} a2(t)Γg2(o2)2 )之差为0。
    这个方程在任何时刻 t t t都成立,基于此,可以利用大量的数据集 ( a 1 ( t k ) , a 2 ( t k ) , g 1 ( t k ) , g 2 ( t k ) ) k = 1 N (a_{1}(t_{k}), a_{2}(t_{k}), g_{1}(t_{k}), g_{2}(t_{k}))_{k = 1}^{N} (a1(tk),a2(tk),g1(tk),g2(tk))k=1N 来估计 o 1 o_{1} o1 o 2 o_{2} o2

  4. 方程(4)解释
    方程 d ( ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 ) d o i = Γ g i T ( a i ( t ) − Γ g i ( o i ) ) ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 \frac{d\left(\left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2}\right)}{d o_{i}}=\frac{\Gamma_{g_{i}}^{T}(a_{i}(t)-\Gamma_{g_{i}}(o_{i}))}{\left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2}} doid(ai(t)Γgi(oi)2)=ai(t)Γgi(oi)2ΓgiT(ai(t)Γgi(oi)) i = 1 , 2 i = 1,2 i=1,2 ,这是关于 ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 \left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2} ai(t)Γgi(oi)2 o i o_{i} oi 的梯度公式。
    其中 Γ g i T ( o i ) = ( o i × g i ( t ) ) × g i ( t ) + o i × g ˙ i ( t ) \Gamma _{g_{i}}^{T}(o_{i})=(o_{i}\times g_{i}(t))\times g_{i}(t)+o_{i}\times \dot {g}_{i}(t) ΓgiT(oi)=(oi×gi(t))×gi(t)+oi×g˙i(t) ,它是 Γ g i \Gamma_{g_{i}} Γgi 矩阵乘法表示形式的转置。
    这个梯度公式在优化算法中非常重要,用于计算当 o i o_{i} oi 发生微小变化时, ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 \left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2} ai(t)Γgi(oi)2 的变化率,进而通过迭代优化方法来求解 o 1 o_{1} o1 o 2 o_{2} o2 ,以满足方程(3)的约束条件,实现从测量数据中准确估计关节中心到传感器坐标系原点的偏移向量。


2. 球窝关节约束的人工解读

绕关节中心旋转而产生的(径向和切向)加速度的人工推导:

球窝关节的约束原理

关节中心加速度的人工推导如下.

球窝关节的加速度

说明

为了推导方便,我们构建的全局坐标系. 而论文中都是基于瞬时固定的 IMU 局部坐标系来推导. 事实上, 瞬时固定的局部坐标系看作为了参考坐标系, 作用类似于全局坐标系. 瞬时固定的局部坐标系并不是固定在 IMU 上的随动局部坐标系, 而是 IMU 运动瞬时在外部惯性空间中与 IMU 随动坐标系瞬时重合的外部参考坐标系.

基于瞬时固定的 IMU 局部坐标系与我们引入的全局坐标系相差一个旋转矩阵. 相应得到的速度和加速度关系也相差一个旋转矩阵.


3. 球窝关节约束梯度的符号计算

d ( ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 ) d o i = Γ g i T ( a i ( t ) − Γ g i ( o i ) ) ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 (4) \frac{d\left(\left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2}\right)}{d o_{i}}=\frac{\Gamma_{g_{i}}^{T}(a_{i}(t)-\Gamma_{g_{i}}(o_{i}))}{\left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2}} \tag{4} doid(ai(t)Γgi(oi)2)=ai(t)Γgi(oi)2ΓgiT(ai(t)Γgi(oi))(4)

简写 Δ i ≜ a i ( t ) − Γ g i ( o i ) \Delta_i \triangleq a_{i}(t)-\Gamma_{g_{i}}(o_{i}) Δiai(t)Γgi(oi). 因为 ∥ Δ i ∥ 2 = Δ i ⋅ Δ i \Vert \Delta_i \Vert_2 = \sqrt{\Delta_i\cdot \Delta_i} Δi2=ΔiΔi , 可知式 (4) 左手边为
d ( ∥ Δ i ∥ 2 ) d o i = d Δ i ⋅ Δ i d o i = 1 2 Δ i ⋅ Δ i d ( Δ i ⋅ Δ i ) d o i = 1 2 ∥ Δ i ∥ 2 d ( Δ i ⋅ Δ i ) d o i (4-a) \frac{d (\Vert \Delta_i \Vert_2 )}{d o_i} = \frac{d \sqrt{\Delta_i \cdot \Delta_i}}{d o_i} = \frac{1}{2\sqrt{\Delta_i \cdot \Delta_i}} \frac{d(\Delta_i \cdot \Delta_i)}{d o_i} = \frac{1}{2\Vert{\Delta_i }\Vert_2} \frac{d(\Delta_i \cdot \Delta_i)}{d o_i} \tag{4-a} doid(Δi2)=doidΔiΔi =2ΔiΔi 1doid(ΔiΔi)=2∥Δi21doid(ΔiΔi)(4-a)
式 (4) 右手边可以简写为
Γ g i T ( a i ( t ) − Γ g i ( o i ) ) ∥ a i ( t ) − Γ g i ( o i ) ∥ 2 = Γ g i T ( Δ i ) ∥ Δ i ∥ 2 (4-b) \frac{\Gamma_{g_{i}}^{T}(a_{i}(t)-\Gamma_{g_{i}}(o_{i}))}{\left\| a_{i}(t)-\Gamma_{g_{i}}(o_{i})\right\| _{2}} = \frac{\Gamma_{g_{i}}^{T}(\Delta_i)}{\left\| \Delta_i \right\| _{2}} \tag{4-b} ai(t)Γgi(oi)2ΓgiT(ai(t)Γgi(oi))=Δi2ΓgiT(Δi)(4-b)
所以要证明式 (4) 成立, 只要验证
1 2 d ( Δ i ⋅ Δ i ) d o i = Γ g i T ( Δ i ) (4-c) \frac{1}{2} \frac{d(\Delta_i \cdot \Delta_i)}{d o_i} = {\Gamma_{g_{i}}^{T}(\Delta_i)} \tag{4-c} 21doid(ΔiΔi)=ΓgiT(Δi)(4-c)
成立即可.

下面的推导计算比较复杂, 我们借用 Maxima 符号计算来推导. (不过结果和论文中相差一个负号)

/*Thomas Seel, etc. Joint Axis and Position Estimation from Inertial Measurement Data by Exploiting Kinematic Constraints*/
/*B. Constraints induced by sphernoidal joints*/
/* Derivation of Equation (4)*/
/*a---acceleration of IMU*/
/*o---offset vestor from the joint center to the origin of the IMU frame*/
/*g---angular velocity of the gyroscope(IMU)*/
/*dg---dg(t)/dt*/
/*Γ---the (radial and tangential) acceleration due to rotation around the joint center(spheroidal joint)*/
/*Gnorm---function to be derivated at left hand of equation (4). Derivatives with respect to vectors*/
/*DGnormX---X-component of the left hand of equation (4). Derivatives with respect to vectors*/
/*DGnormY---Y-component of the left hand of equation (4)*/
/*DGnormZ---Z-component of the left hand of equation (4)*/a:[a_x,a_y,a_z];
o:[o_x,o_y,o_z];
g:[g_x,g_y,g_z];
dg:[dg_x,dg_y,dg_z];
load("vect");
express(g~o);
Γ:(g~(g~o)+dg~o);
Γ2:express(Γ);
Γ2a:(a-Γ2);Gnorm: Γ2a.Γ2a$
DGnormX: expand(diff(Gnorm,o_x)/2);
/*numerator of the right hand of equation (4)*/:(Γ2a~g)~g+Γ2a~dg$
TΓ2:express()$/*X-component of the numerator of the right hand of equation (4)*/
expand(2[1]);DGnormY:expand(diff(Gnorm,o_y)/2);
/*Y-component of the numerator of the right hand of equation (4)*/
expand(2[2]);DGnormZ:expand(diff(Gnorm,o_z)/2);
/*Z-component of the numerator of the right hand of equation (4)*/
expand(2[3]);
eq4_1 eq4-2

IV. 算法实现

  1. 初始化
    设定输入数据,即获取 N N N 个数据集 ( g 1 ( t k ) , g 2 ( t k ) ) k = 1 N (g_{1}(t_{k}), g_{2}(t_{k}))_{k = 1}^{N} (g1(tk),g2(tk))k=1N,其中 N ≫ 4 N \gg 4 N4 。这些数据代表在不同时刻 t k t_{k} tk 从两个传感器获取的与关节运动相关的信息。
    初始化变量 x x x,在估计球窝关节偏移向量时, x x x o 1 o_{1} o1 o 2 o_{2} o2 的拼接。这里 o 1 o_{1} o1 o 2 o_{2} o2 分别是从关节中心到第一和第二传感器坐标系原点的偏移向量。
    (铰接关节的 x x x 由公式 (5) 第一式定义)

  2. 关节轴估计
    利用公式(5),基于当前的 x x x 值计算关节轴估计值 j 1 ^ \hat{j_{1}} j1^ j 2 ^ \hat{j_{2}} j2^
    公式(5)建立了 x x x j 1 ^ \hat{j_{1}} j1^ j 2 ^ \hat{j_{2}} j2^ 之间的数学关系,通过将 x x x 的各个分量代入公式(5),按照规定的数学运算顺序进行计算,得出 j 1 ^ \hat{j_{1}} j1^ j 2 ^ \hat{j_{2}} j2^ 的值。
    同时,为保证估计的合理性,将关节轴估计值限制为单位长度,使得估计问题转化为四维。

  3. 计算误差与梯度
    [球窝关节] 定义误差函数 e ( t ) ≜ ∥ a 1 ( t ) − Γ g 1 ( o 1 ) ∥ 2 − ∥ a 2 ( t ) − Γ g 2 ( o 2 ) ∥ 2 e(t) \triangleq \left\|a_{1}(t) - \Gamma_{g_{1}}(o_{1})\right\|_{2} - \left\|a_{2}(t) - \Gamma_{g_{2}}(o_{2})\right\|_{2} e(t)a1(t)Γg1(o1)2a2(t)Γg2(o2)2 。其中 a 1 ( t ) a_{1}(t) a1(t) a 2 ( t ) a_{2}(t) a2(t) 分别是两个传感器在时刻 t t t 的加速度测量值, Γ g 1 ( o 1 ) \Gamma_{g_{1}}(o_{1}) Γg1(o1) Γ g 2 ( o 2 ) \Gamma_{g_{2}}(o_{2}) Γg2(o2) 是与关节旋转相关的加速度项,通过此误差函数衡量当前估计的偏移向量 o 1 o_{1} o1 o 2 o_{2} o2 的准确性。
    [铰接关节] 定义误差函数 e ( k ) ≜ ∥ j 1 ^ × g 1 ( t k ) ∥ 2 − ∥ j 2 ^ × g 2 ( t k ) ∥ 2 e(k) \triangleq \left\|\hat{j_{1}} × g_{1}(t_{k})\right\|_{2}-\left\|\hat{j_{2}} × g_{2}(t_{k})\right\|_{2} e(k) j1^×g1(tk) 2 j2^×g2(tk) 2,即两个范数的差值。如果关节轴估计准确,那么对于每个时刻 t k t_{k} tk ∥ j 1 ^ × g 1 ( t k ) ∥ 2 \left\|\hat{j_{1}} × g_{1}(t_{k})\right\|_{2} j1^×g1(tk) 2 ∥ j 2 ^ × g 2 ( t k ) ∥ 2 \left\|\hat{j_{2}} × g_{2}(t_{k})\right\|_{2} j2^×g2(tk) 2 的值应该接近,此时 e ( k ) e(k) e(k) 接近零。通过分析误差向量 e e e 的值,可以评估关节轴估计的准确性,并进一步用于优化估计过程,例如在迭代算法中,根据误差向量 e e e 来调整 j 1 ^ \hat{j_{1}} j1^ j 2 ^ \hat{j_{2}} j2^ 的值,使得误差逐渐减小,从而得到更准确的关节轴估计。
    根据公式(4)计算误差 e e e 关于 x x x 的梯度 d e d x \frac{d e}{d x} dxde 。公式(4)给出了关于误差对 x x x 的导数形式,在计算过程中,需根据 x x x 的具体形式(即 o 1 o_{1} o1 o 2 o_{2} o2 的拼接)以及误差函数 e e e 的表达式,按照公式(4)的规则进行详细的数学推导和计算,得到梯度值。

  4. 更新变量 x x x
    使用公式 x = x − p i n v ( d e d x ) e x = x - pinv(\frac{d e}{d x}) e x=xpinv(dxde)e 更新变量 x x x 。这里 p i n v ( d e d x ) pinv(\frac{d e}{d x}) pinv(dxde) 表示梯度 d e d x \frac{d e}{d x} dxde 的伪逆,通过该公式,让 x x x 沿着负梯度方向(乘以误差 e e e)进行更新,步长由伪逆矩阵 p i n v ( d e d x ) pinv(\frac{d e}{d x}) pinv(dxde) 决定,使得误差 e e e 逐渐减小。

  5. 循环与终止条件
    重复步骤 2 至 4,直到满足特定的终止条件。终止条件可以设定为误差 e e e 小于某个预设的阈值,表明算法已经收敛到一个足够准确的结果;或者达到一定的迭代次数,防止算法无限循环。

  6. 噪声处理
    在计算过程中,由于 g 1 ( t ) g_{1}(t) g1(t) g 2 ( t ) g_{2}(t) g2(t) 的时间导数对计算结果至关重要,且数据可能存在噪声干扰,采用一种抗噪声的非因果低通滤波器与简单的差商近似相结合的方法。
    非因果低通滤波器用于去除数据中的高频噪声,使数据更加平滑;简单的差商近似则用于从离散的测量数据中估计 g 1 ( t ) g_{1}(t) g1(t) g 2 ( t ) g_{2}(t) g2(t) 的时间导数,从而保证在噪声环境下算法仍能准确运行。
    通过以上算法实现步骤,能够有效地估计关节轴及球窝关节的偏移向量,为后续相关问题的解决提供可靠的数据支持。


V. 小结

以上推导了利用 IMU 估计人体关节轴向与位置的数学原理. 该算法可用于外骨骼控制、人体运动跟踪等场景.

写这边博文 AI 帮了很大忙, 提高了效率. 但我在这篇博文中的作用却降低了吧?


版权声明:本文为博主原创文章,遵循 CC 4.0 BY 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/woyaomaishu2/article/details/145503633
本文作者:wzf@robotics_notes

相关文章:

利用 IMU 估计人体关节轴向和位置 —— 论文推导

Title: 利用 IMU 估计人体关节轴向和位置 —— “Joint axis and position estimation from inertial measurement data by exploiting kinematic constraints” —— 论文推导 文章目录 I. 论文回顾II. 铰接关节的约束1. 铰接关节约束的原理2. 铰接关节约束的梯度3. 铰接关节约…...

脚本一键生成管理下游k8s集群的kubeconfig

一、场景 1.1 需要管理下游k8s集群的场景。 1.2 不希望使用默认的cluster-admin权限的config. 二、脚本 **重点参数: 2.1 配置变量。 1、有单独namespace的权限和集群只读权限。 2、自签名的CA证书位置要正确。 2.2 如果配置错误,需要重新…...

数据库系统概念第六版记录 三

外码约束(Foreign Key Constraint) 外码(Foreign Key, FK)是关系数据库中的一个约束,它用于保证表之间的引用完整性。外码的值必须: 要么存在于被引用表的主键列中,要么为空(NULL&…...

YOLOv11-ultralytics-8.3.67部分代码阅读笔记-files.py

files.py ultralytics\utils\files.py 目录 files.py 1.所需的库和模块 2.class WorkingDirectory(contextlib.ContextDecorator): 3.def spaces_in_path(path): 4.def increment_path(path, exist_okFalse, sep"", mkdirFalse): 5.def file_age(path__fi…...

微信小程序案例1——制作猫眼电影底部标签导航栏

文章目录 一、项目步骤1 新建一个无AppID的movie项目2将准备好的底部标签导航图标拷贝到movie项目下面(将图标文件夹image放到项目文件夹里)3 打开App.json配置文件,在pages数组里添加4个页面路径:电影“pages/movie/movie”、影院“pages/cinema/cinema…...

【大数据技术】搭建完全分布式高可用大数据集群(Kafka)

搭建完全分布式高可用大数据集群(Kafka) kafka_2.13-3.9.0.tgz注:请在阅读本篇文章前,将以上资源下载下来。 写在前面 本文主要介绍搭建完全分布式高可用集群 Kafka 的详细步骤。 注意: 统一约定将软件安装包存放于虚拟机的/software目录下,软件安装至/opt目录下。 安…...

【服务器知识】如何在linux系统上搭建一个nfs

文章目录 NFS网络系统搭建**1. 准备工作****2. 服务器端配置****(1) 安装 NFS 服务****(2) 创建共享目录****(3) 配置共享规则****(4) 生效配置并启动服务****(5) 防火墙配置** **3. 客户端配置****(1) 安装 NFS 客户端工具****(2) 创建本地挂载点****(3) 挂载 NFS 共享目录***…...

图片画质增强:轻松提升画质

前言: 今天给大家推荐一款非常实用的图片画质增强软件,它无需联网即可使用,完全离线操作,这款软件基于先进的深度学习技术,能够对模糊图片进行强大的高清处理,效果令人惊艳。 图片画质增强:一…...

vscode快速接入deepseek 实践操作

背景说明 在deepseek快速火爆的情况下,也想自己体验一把。看看在vscode中集成进来,方便平时的脚本开发。对于年纪大的人还是非常方便的。 操作过程 安装continue 打开vscode进入扩展市场,搜索安装 安装完成就是上面的样子,会…...

mapbox进阶,添加绘图扩展插件,绘制圆形

👨‍⚕️ 主页: gis分享者 👨‍⚕️ 感谢各位大佬 点赞👍 收藏⭐ 留言📝 加关注✅! 👨‍⚕️ 收录于专栏:mapbox 从入门到精通 文章目录 一、🍀前言1.1 ☘️mapboxgl.Map 地图对象1.2 ☘️mapboxgl.Map style属性1.3 ☘️MapboxDraw 绘图控件二、🍀添加绘图扩…...

Cursor 与多语言开发:全栈开发的利器

引言 全栈开发要求开发者跨越前端、后端、数据库甚至数据科学等多个技术领域,而不同技术栈往往需要切换工具和思维方式。Cursor 作为一款 AI 驱动的智能编程助手,凭借其对 20 编程语言 和主流框架的深度支持,正在成为全栈开发的“瑞士军刀”…...

2025 CCF BDCI|“基于TPU平台的OCR模型性能优化”一等奖作品

2024年12月,中国计算机学会在海南博鳌成功举办了第十二届CCF大数据与计算智能大赛(简称2024 CCF BDCI)。本届比赛的算能赛道吸引了1748名选手报名,经过激烈角逐,北京航空航天大学的“常务副SOTA”团队脱颖而出&#xf…...

FPGA的IP核接口引脚含义-快解

疑问 手册繁琐,怎样快速了解IP核各输入输出接口引脚的含义。 答疑 不慌不慌,手册确实比较详细但繁琐,如何快速知晓该部分信息,涛tao道长给你们说,简单得很,一般新入门的道友有所不知,往往后面…...

数据库高安全—审计追踪:传统审计统一审计

书接上文数据库高安全—角色权限:权限管理&权限检查,从权限管理和权限检查方面解读了高斯数据库的角色权限,本篇将从传统审计和统一审计两方面对高斯数据库的审计追踪技术进行解读。 4 审计追踪 4.1 传统审计 审计内容的记录方式通…...

机器学习 - 需要了解的条件概率、高斯分布、似然函数

似然函数是连接数据与参数的桥梁,通过“数据反推参数”的逆向思维,成为统计推断的核心工具。理解它的关键在于区分“参数固定时数据的概率”与“数据固定时参数的合理性”,这种视角转换是掌握现代统计学和机器学习的基础。 一、在学习似然函…...

Spring Boot Web 入门

目录 Spring Boot Web 是 Spring Boot 框架的一个重要模块,它简化了基于 Spring 的 Web 应用程序的开发过程。以下是一个 Spring Boot Web 项目的入门指南,涵盖了项目创建、代码编写、运行等关键步骤。 1. 项目创建 使用 Spring Initializr 使用 IDE …...

神经网络|(八)概率论基础知识-二项分布及python仿真

【1】引言 前序已经学习了古典概型、条件概率、全概率公式和贝叶斯公式,它们作为基础,解释了事件发生及其概率的对应关系,相关文章链接为: 神经网络|(四)概率论基础知识-古典概型-CSDN博客 神经网络|(五)概率论基础知识-条件概…...

【面试场景】MySQL分布式主键选取

文章目录 一. MySQL的自增主键二. UUID三. 雪花ID(推荐) 我的博客地址 一. MySQL的自增主键 适合单表的情况, 在分布式分库分表下可能会有一些问题 主键冲突问题 在分布式系统中,多个数据库节点独立生成自增主键,很容易出现重复的主键值。例如&#xff…...

执行git stash drop stash@{x} 时出现error: unknown switch `e‘ 的解决方式

原因: 在 PowerShell 或某些 Shell 中,{} 是特殊符号,stash{0} 会被解析成 stash 0,而 后的字符可能被误认为选项(如 -e),使 Git 收到意外的 -e 参数,导致报错 unknown switch ‘e’。 解决方…...

链表和 list

一、单链表的模拟实现 1.实现方式 链表的实现方式分为动态实现和静态实现两种。 动态实现是通过 new 申请结点,然后通过 delete 释放结点的形式构造链表。这种实现方式最能体 现链表的特性; 静态实现是利用两个数组配合来模拟链表。一个表示数据域&am…...

2025年能源电力系统与流体力学国际会议 (EPSFD 2025)

2025年能源电力系统与流体力学国际会议(EPSFD 2025)将于本年度在美丽的杭州盛大召开。作为全球能源、电力系统以及流体力学领域的顶级盛会,EPSFD 2025旨在为来自世界各地的科学家、工程师和研究人员提供一个展示最新研究成果、分享实践经验及…...

Admin.Net中的消息通信SignalR解释

定义集线器接口 IOnlineUserHub public interface IOnlineUserHub {/// 在线用户列表Task OnlineUserList(OnlineUserList context);/// 强制下线Task ForceOffline(object context);/// 发布站内消息Task PublicNotice(SysNotice context);/// 接收消息Task ReceiveMessage(…...

DAY 47

三、通道注意力 3.1 通道注意力的定义 # 新增:通道注意力模块(SE模块) class ChannelAttention(nn.Module):"""通道注意力模块(Squeeze-and-Excitation)"""def __init__(self, in_channels, reduction_rat…...

UE5 学习系列(三)创建和移动物体

这篇博客是该系列的第三篇,是在之前两篇博客的基础上展开,主要介绍如何在操作界面中创建和拖动物体,这篇博客跟随的视频链接如下: B 站视频:s03-创建和移动物体 如果你不打算开之前的博客并且对UE5 比较熟的话按照以…...

在 Nginx Stream 层“改写”MQTT ngx_stream_mqtt_filter_module

1、为什么要修改 CONNECT 报文? 多租户隔离:自动为接入设备追加租户前缀,后端按 ClientID 拆分队列。零代码鉴权:将入站用户名替换为 OAuth Access-Token,后端 Broker 统一校验。灰度发布:根据 IP/地理位写…...

Java面试专项一-准备篇

一、企业简历筛选规则 一般企业的简历筛选流程:首先由HR先筛选一部分简历后,在将简历给到对应的项目负责人后再进行下一步的操作。 HR如何筛选简历 例如:Boss直聘(招聘方平台) 直接按照条件进行筛选 例如&#xff1a…...

IP如何挑?2025年海外专线IP如何购买?

你花了时间和预算买了IP,结果IP质量不佳,项目效率低下不说,还可能带来莫名的网络问题,是不是太闹心了?尤其是在面对海外专线IP时,到底怎么才能买到适合自己的呢?所以,挑IP绝对是个技…...

基于Springboot+Vue的办公管理系统

角色: 管理员、员工 技术: 后端: SpringBoot, Vue2, MySQL, Mybatis-Plus 前端: Vue2, Element-UI, Axios, Echarts, Vue-Router 核心功能: 该办公管理系统是一个综合性的企业内部管理平台,旨在提升企业运营效率和员工管理水…...

三分算法与DeepSeek辅助证明是单峰函数

前置 单峰函数有唯一的最大值,最大值左侧的数值严格单调递增,最大值右侧的数值严格单调递减。 单谷函数有唯一的最小值,最小值左侧的数值严格单调递减,最小值右侧的数值严格单调递增。 三分的本质 三分和二分一样都是通过不断缩…...

MySQL的pymysql操作

本章是MySQL的最后一章,MySQL到此完结,下一站Hadoop!!! 这章很简单,完整代码在最后,详细讲解之前python课程里面也有,感兴趣的可以往前找一下 一、查询操作 我们需要打开pycharm …...