一种基于运动轨迹的乒乓球旋转估计方法及其应用

1.本发明公开了一种基于运动轨迹的乒乓球旋转估计方法，属于机器人技术领域。


背景技术：

2.乒乓球作为我国的国球，对于我国的体育事业具有重大的影响。随着计算机、光学、自动控制等新型技术的兴起，除了人人对打这一传统的竞技模式以外，人机对打(vr游戏、人
ꢀ‑
机器人对打)逐渐步入了我们的生活中，大大增强了乒乓球运动的趣味性。因此，如何尽可能地将真实的乒乓球运动结合到虚拟世界中，如何高效地感知乒乓球的运动性质，提升人类
‑ꢀ
机器人乒乓球对战的效果是急需解决的问题。
3.旋转是乒乓球运动的精髓，而高速旋转是乒乓球运动相比于其他球类运动所没有的，因此，准确地估计旋转是实现人机对打的关键。乒乓球的飞行轨迹和飞行线速度、旋转角速度、空气阻力系数、马格努斯力系数有关，且呈高度非线性。球的飞行方式也确定了球上商标的运动特性。如何利用有效信息(如球轨迹、球商标)准确估算球的旋转参数是至关重要的一步。
4.现有方法虽然能在某种程度上解决旋转估计的问题，然而仍有一定的局限性，具体表现如下：
5.1、基于商标识别的旋转估计很大程度上依赖于长焦镜头、高帧率的相机和快速的视觉实时处理，对硬件的要求较高，且需要保证球能处在相机的合适视野范围内，大大增加了工作量。
6.2、现有的基于运动轨迹的旋转估计中，较为常见的做法是将空气阻力系数(kd)、马格努斯力系数(km)看作预先设好的常值。这其中有两点局限性：a.kd,km与空气湿度、温度、球的线速度和旋转角速度等参数有关，直接将其定为常值无法囊括不同球速、不同空气温度的情况；b.事实上kd,km与球的线速度和旋转角速度的比值呈非线性关系，所以球在飞行过程中随着球线速度的衰减，kd,km也势必有变化，不符合“常数”的假设。
7.现有公开号为cn110941795a的专利提出了一种基于随机采样一致算法的旋转角速度估计方法。该方法首先通过函数拟合得到球的飞行轨迹，并采用递归最小二乘法获取球的速度和加速度，再计算各时刻马格努斯力的单位向量，最终基于采样一致算法得到乒乓球的旋转角速度。但该方法认为马格努斯力系数提前已知，未考虑其随球线速度、角速度的变化，无法保证旋转估计的精度。
8.现有公开(公告)号为cn110458281a(cn110458281b)的专利提供了一种基于深度强化学习的乒乓球旋转速度预测方法。该方法首先将来球位置进行归一化，再将该归一化序列输入到深度lstm网络中，最终将得到的lstm的状态向量输入到旋转估计网络中得到入射旋转速度。然而，该方法同样未考虑变化的马格努斯力系数对旋转估计的影响，无法保证旋转估计的精度。且这类基于网络的方法需要大量的数据样本，前期工作复杂耗时。


技术实现要素：

9.本发明为了解决上述现有技术的不足之处，提出了一种基于运动轨迹的乒乓球旋转估计方法及其应用，以期能准确地估计球在飞行中的旋转速度，从而能提高乒乓球轨迹预测的精度。
10.为达到上述发明目的，采用如下技术方案：
11.本发明一种基于运动轨迹的乒乓球旋转估计方法的特点在于，所述方法包括：
12.步骤1、基于乒乓球飞行空气动力学，对若干组乒乓球飞行的轨迹进行求解，得到乒乓球飞行过程中的空气阻力系数、马格努斯力系数、球线速度、由马格努斯力引起的乒乓球的加速度，用于拟合马格努斯力系数与球线速度、角速度的关系；
13.步骤1.1，利用发球机发出角速度为ω的乒乓球，在乒乓球飞行过程中，通过双目视觉以采样间隔ts对球位置坐标进行采样，得到乒乓球的轨迹{wp(k):k＝1,2
…
n}；其中，k表示采样时刻数，n表示一条轨迹中采样点的个数，wp(k)表示第k个采样时刻采集到的乒乓球的位置矢量；且wp(k)＝[wx(k),wy(k),wz(k)]，wx(k),wy(k),wz(k)分别表示位置矢量wp(k)在世界坐标系的 x轴、y轴、z轴上的分量；
[0014]
步骤1.2，通过对wx(k),wy(k),wz(k)进行三次多项式拟合，得到如式(1)、式(2)、式(3)所示的乒乓球的坐标与时间t的连续关系：
[0015]wx(t)＝f
x
(wx(k),t)
ꢀꢀ
(1)
[0016]wy(t)＝f
x
(wy(k),t)
ꢀꢀ
(2)
[0017]wz(t)＝f
x
(wz(k),t)
ꢀꢀ
(3)
[0018]
式(1)、式(2)、式(3)中，f
x
(
·
,
·
)、fy(
·
,
·
)、fz(
·
,
·
)分别为世界坐标系的x轴、y轴、z轴上关于时间变量t的三次多项式函数，wx(t),wy(t),wz(t)表示拟合后t时刻乒乓球的位置矢量在世界坐标系的x轴、y轴、z轴的分量；
[0019]
通过对式(1)、式(2)、式(3)求关于t的一阶、二阶导数，并令t＝kts，从而得到如式(4)和式 (5)所示的第k个采样时刻球在世界坐标系下的球线速度wv(k)和加速度wa(k)：
[0020][0021][0022]
式(4)和式(5)中，wv(k)＝[
wvx
(k),wvy(k),wvz(k)]，wa(k)＝[
wax
(k),
way
(k),
waz
(k)]，
wvx
(k),wvy(k),wvz(k) 为wv(k)在世界坐标系的x轴，y轴，z轴的分量，
wax
(k),
way
(k),
waz
(k)为wa(k)在世界坐标系的x 轴，y轴，z轴的分量；
[0023]
步骤1.3，利用式(6)、式(7)、式(8)分别求解第k个采样时刻乒乓球的空气阻力系数kd(k)、由马格努斯力产生的加速度在世界坐标系下的加速度
wam
(k)、马格努斯力系数km(k)：
[0024]
[0025][0026][0027]
式(6)、式(7)中，wg表示在世界坐标系中的重力加速度矢量；
[0028]
步骤1.4，按照步骤1.1-1.3的过程，采集不同角速度，不同线速度，不同发射方向的n 组乒乓球轨迹并得到包含所有的km(k)、wω(k)、wv(k)的集合：y＝{{km(kj)}j}、x1＝{{wω(kj)}j}、 x2＝{{wv(kj)}j},j＝1,2
…
n,kj＝1,2
…
nj，其中，j表示轨迹索引，kj表示第j个轨迹的坐标点索引， y为包含个马格努斯力系数的集合，x1,x2为包含个角速度、线速度的集合；并将 x1，x2作为自变量，y作为自变量进行二维多项式函数拟合，得到如式(9)所示的关系式：
[0029]km
＝f(wv,wω)
ꢀꢀ
(9)。
[0030]
步骤2、基于所述马格努斯力系数与球线速度、角速度的拟合关系，利用乒乓球飞行的新轨迹构建角速度矢量方向的线性方程组并进行求解，得到乒乓球旋转角速度矢量的方向角。
[0031]
步骤3、基于所述乒乓球旋转角速度矢量的方向，利用乒乓球飞行的新轨迹构建旋转角速度模的非线性函数并求解其最小值，得到乒乓球旋转角速度的模。
[0032]
本发明所述的乒乓球旋转估计方法的特点也在于，所述步骤2包括：
[0033]
步骤2.1，给定一条待求的新轨迹td＝{wpd(kd)|kd＝1,2
…
nd},wpd(kd)表示待求的新轨迹中第 kd个采样时刻球在世界坐标系下的位置矢量，nd表示待求的新轨迹中包含的乒乓球坐标的个数；
[0034]
利用式(4)计算每一个采样时刻球在世界坐标系下的速度矢量wvd(k)＝[
wvdx
(k),
wvdy
(k), wvdz
(k)]、利用式(7)计算在世界坐标系下由马格努斯力引起的加速度矢量
wadm
(k)＝[
wadmx
(k), wadmy
(k),
wadmz
(k)]，其中，
wvdx
(k),
wvdy
(k),
wvdz
(k)为wvd(k)在世界坐标系的x轴、y轴、z轴的分量；
wadmx
(k),
wadmy
(kd),
wadmz
(k)为
wadm
(k)在世界坐标系的x轴、y轴、z轴的分量；
[0035]
步骤2.2，按照式(10)计算待求的新轨迹中第kd个时刻乒乓球的相对坐标系和世界坐标系的转换矩阵
wrdl
(k)：
[0036][0037]
式(10)中，r
d11
(k),r
d12
(k),r
d21
(k),r
d22
(k)为矩阵
wrdl
(k)中第1行第1列，第1行第2列，第 2行第1列，第2行第2列的元素；
[0038]
步骤2.3，利用式(11)和式(12)分别计算关于旋转角速度矢量方向角线性方程的左系数矩阵 ld和右系数矩阵rd：
[0039]
ld＝[a
dmz
(1) a
dmz
(2)
ꢀ…ꢀadmz
(nd)]
t
ꢀꢀ
(11)
[0040][0041]
步骤2.4，计算旋转角速度矢量的方向角
[0042]
a、通过式(13)计算θd的正切值tanθd：
[0043]
tanθd＝(l
dt
ld)-1
l
dt
rdꢀꢀ
(13)
[0044]
b、通过式(14)计算旋转方向矢量的方向角
[0045][0046]
所述步骤3包括：
[0047]
步骤3.1，根据步骤2中得到的方向角用式(15)计算第kd个采样时刻角速度矢量wωd和速度矢量wvd(kd)的夹角θ
dωv
(kd)：
[0048][0049]
式(15)中，
wedω
表示与wωd同向的单位向量，并有：
[0050][0051]
步骤3.2，利用式(17)计算观测变量y2：
[0052]
y2＝[||
wadm
(kd)|||kd＝1，2...nd]
ꢀꢀ
(17)
[0053]
步骤3.3，通过求解函数最小值求解旋转角速度矢量的模：
[0054]
a、令自变量x3＝[||wvd(kd)|||kd＝1,2
…
nd],x4＝[θ
dωv
(kd)|kd＝1,2
…
nd]；
[0055]
b、利用式(18)构建拟合函数y：
[0056]
y＝g(x3(kd)，||ωd||，x4(kd))＝f(x3(kd)，||ωd||)
·
x3(kd)||ωd||sin(x4(kd))
ꢀꢀ
(18)
[0057]
式(18)中，||ωd||为非线性回归的待定系数；
[0058]
c、利用式(19)构建非线性函数并求解最小值，从而得到乒乓球旋转角速度的模
[0059][0060]
本实施例中，一种电子设备，包括存储器以及处理器，其特点在于，所述存储器用于存储支持处理器执行所述乒乓球旋转估计方法的程序，所述处理器被配置为用于执行所述存储器中存储的程序。
[0061]
本实施例中，一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，其特点在于，所述计算机程序被处理器运行时执行所述乒乓球旋转估计方法的步骤。
[0062]
与现有技术相比，本发明的有益效果在于：
[0063]
1、本发明通过实时求解空气阻力系数、马格努斯力系数，解决了现有技术中将空气阻力系数、马格努斯力系数看作预先知道的常数，未考虑球的线速度、角速度对该常数影响的问题，从而提高了乒乓球旋转估计的精度。
[0064]
2、本发明步骤1中，通过球的运动轨迹，利用球飞行动力学的特征，实时求解空气阻力系数、马格努斯力系数，有效避免了“系数为常值”这一假设，避免了以往专利中“根据当前空气温度、湿度提前计算空气阻力系数，马格努斯力系数”这一做法，从而提高了乒乓球旋转估计的精度。
[0065]
3、本发明步骤1中，针对实时求解的马格努斯力系数和球线速度、旋转角速度，拟合出了二者的函数关系，并将其用在了旋转估计上，克服了现有技术中“未考虑马格努斯力系数和球角速度、线速度关系”的问题，从而提高了乒乓球旋转估计的精度。
附图说明
[0066]
图1是本发明实施例中描述世界坐标系、相对坐标系、球线速度、旋转角速度的示意图；
[0067]
图2是本发明实施例中描述拟合km与角速度、线速度关系这一过程的流程图；
[0068]
图3是本发明实施例中描述旋转角速度矢量方向求解过程的流程图；
[0069]
图4是本发明实施例中描述绝对坐标系与相对坐标系姿态关系的示意图；
[0070]
图5是本发明实施例中描述旋转角速度矢量模求解过程的流程图。
具体实施方式
[0071]
本实施例中，一种基于运动轨迹的乒乓球旋转估计方法，包括：
[0072]
步骤1、基于乒乓球飞行空气动力学，对若干组乒乓球飞行的轨迹进行求解，得到乒乓球飞行过程中的空气阻力系数、马格努斯力系数、球线速度、由马格努斯力引起的乒乓球的加速度，用于拟合马格努斯力系数与球线速度、角速度的关系；
[0073]
如图1所示建立世界坐标系o-xyz，y轴垂直于桌面的宽，平行于桌面且指向远离机器人的方向，z轴竖直向上，通过右手规则确定x轴方向。
[0074]
利用空气动力学和牛顿第二定律对乒乓球的飞行轨迹进行建模。乒乓球在飞行过程主要受到重力、空气阻力、马格努斯力的作用，由此得到球的飞行动力学方程如式(1)所示：
[0075]wa＝wg+kd||wv||wv+k
mw
ω
×
wvꢀꢀ
(1)
[0076]
式(1)中，wa表示球在世界坐标系下的加速度，wg表示在世界坐标系下重力加速度矢量，wv,wω表示在世界坐标系球的线速度矢量和角速度矢量，kd,km分别表示空气阻力系数和马格努斯力系数。
[0077]
步骤1.1，采集若干组乒乓球轨迹，目的是产生足够的数据样本，具体实施方式如下：利用发球机发出角速度为ω的乒乓球，在乒乓球飞行过程中，通过双目视觉以采样间隔ts对球位置坐标进行采样，得到乒乓球的轨迹{wp(k):k＝1,2
…
n}，其中，k表示采样时刻数，
n表示一条轨迹中采样点的个数，wp(k)表示第k个采样时刻采集到的乒乓球在世界坐标系的位置矢量；且wp(k)＝[wx(k),wy(k),wz(k)]，wx(k),wy(k),wz(k)分别表示位置矢量wp(k)在世界坐标系的x轴、y 轴、z轴上的分量。
[0078]
步骤1.2，为了减小双目视觉系统噪声对球速度、加速度的影响，本技术采用多项式拟合的方法对每个采样时刻的球坐标进行拟合，从而起到滤波的作用，然后求解每个采样时刻球的运动速度和运动加速度，具体步骤为，通过对wx(k),wy(k),wz(k)进行3次多项式拟合，得到如式(2)、式(3)、式(4)所示的乒乓球的坐标与时间t的连续关系：
[0079]wx(t)＝f
x
(wx(k),t)
ꢀꢀ
(2)
[0080]wy(t)＝f
x
(wy(k),t)
ꢀꢀ
(3)
[0081]wz(t)＝f
x
(wz(k),t)
ꢀꢀ
(4)
[0082]
式(2)、式(3)、式(4)中，f
x
(
·
,
·
),fy(
·
,
·
),fz(
·
,
·
)分别为世界坐标系的x轴、y轴、z轴上关于时间变量t的三次多项式函数，wx(t),wy(t),wz(t)表示拟合后t时刻乒乓球的位置矢量在世界坐标系的x轴、y轴、z轴的分量。
[0083]
通过对式(2)、式(3)、式(4)求关于t的一阶、二阶导数，并令t＝kts,从而得到如式(5)和(6) 所示的第k个采样时刻球在世界坐标系下的球线速度wv(k)和加速度wa(k)：
[0084][0085][0086]
式(5)和式(6)中，wv(k)＝[
wvx
(k),wvy(k),wvz(k)]，wa(k)＝[
wax
(k),
way
(k),
waz
(k)]，
wvx
(k),wvy(k),wvz(k) 为wv(k)在世界坐标系的x轴，y轴，z轴的分量，
wax
(k),
way
(k),
waz
(k)为wa(k)在世界坐标系的x 轴，y轴，z轴的分量。
[0087]
步骤1.3，为了减小双目视觉系统噪声对球速度、加速度的影响，本技术采用多项式拟合的方法对每个采样时刻的球坐标进行拟合，从而起到滤波的作用，然后求解每个采样时刻球的运动速度和运动加速度，具体步骤为，利用式(7)、式(8)、式(9)分别求解第k个采样时刻乒乓球的空气阻力系数kd(k)、由马格努斯力产生的加速度在世界坐标系下的加速度
wam
(k)、马格努斯力系数km(k)：
[0088][0089][0090][0091]
步骤1.4，步骤1.4，根据球飞行流体力学相关理论，km为球飞行线速度与飞行角速度的函数。为得到具体的函数表达式，本技术采取的方法为：收集若干组球飞行轨迹，利用式(7)、式(8)、式(9)计算每个采样时刻的km(k),wv(k)和||ω||，再用函数拟合的方法近似
得到km(k)与 ||wv(k)||和||ω||的函数关系式。具体实施为：如图2所示，按照步骤1.1-步骤1.3的过程，采集不同角速度，不同线速度，不同发射方向的n组乒乓球轨迹并得到得到包含所有的km(k)、wv(k)、wω(k)的集合：y＝{{km(kj)}j}与x1＝{{wω(kj)}j},x2＝{{wv(kj)}j},j＝1,2
…
n,kj＝1,2
…
nj，其中，j表示轨迹索引，kj表示第j个轨迹的坐标点索引，y为包含了个马格努斯力系数的集合，x1,x2为包含了个角速度、线速度的集合；并将x1，x2作为自变量，y作为因变量进行二维多项式函数拟合，得到如式(10)的关系：
[0092]km
＝f(wv,wω)
ꢀꢀ
(10)
[0093]
步骤2，如图3所示，给定一条来球轨迹，计算其旋转角速度矢量方向。
[0094]
如图1所示建立在球上的相对坐标系{ox
lylzl
},ox
l
为球的速度矢量在世界坐标系xwoyw平面上的投影，oz
l
竖直向上，oy
l
由右手规则决定。根据式(1)可知，球飞行过程中马格努斯力大小与飞行角速度直接相关，所以本技术利用马格努斯力的受力方向推导出球旋转矢量的方向，具体实施方式如下：
[0095]
步骤2.1，给定待求轨迹，计算每个采样时刻球的运动速度、运动加速度和由马格努斯力引起的加速度，具体实施方式如下：给定一条待求的新轨迹td＝{wpd(kd)|kd＝1,2
…
nd},wpd(kd) 表示待求的新轨迹中第kd个采样时刻球在世界坐标系下的位置矢量，nd表示待求的新轨迹中包含的乒乓球坐标的个数。
[0096]
利用式(5)计算每一个采样时刻球在世界坐标系下的速度矢量wvd(k)＝[
wvdx
(k),
wvdy
(k), wvdz
(k)]、利用式(8)计算由马格努斯力引起的在世界坐标系下的加速度矢量
wadm
(k)＝[
wadmx
(k), wadmy
(k),
wadmz
(k)]，其中，
wvdx
(k),
wvdy
(k),
wvdz
(k)为wvd(k)在世界坐标系x轴、y轴、z轴的分量。
wadmx
(k),
wadmy
(k),
wadmz
(k)为
wadm
(k)在世界坐标系下的x轴、y轴、z轴的分量；
[0097]
步骤2.2，由于球在飞行过程中相对坐标系{ox
lylzl
}的方向时时刻刻都在改变，故需求出每个采样时刻{ox
lylzl
}和{ox
wywzw
}的转换矩阵。具体实施方式为：按照式(11)计算待求的新轨迹中第kd个时刻乒乓球的相对坐标系和世界坐标系的转换矩阵
wrdl
(k)：
[0098][0099]
式(11)中，r
d11
(k),r
d12
(k),r
d21
(k),r
d22
(k)为矩阵
wrdl
(k)中第1行第1列，第1行第2列，第 2行第1列，第2行第2列的元素；式(11)的推导如下：
[0100]
如图4所示，建立新坐标系{ox
w’y
w’z
w’}，该坐标系的原点位于球心，ox
w’,oy
w’,oz
w’的方向分别与oxw,oyw,ozw的方向相同，v
xy
为球的线速度在平面{x
w’oy
w’}上的投影，w′vx
,w′
vy为 v
xy
在ox
w’,oy
w’上的投影，故有{ox
lylzl
}的三个与坐标轴同向的单位向量
weoxl
,
weoyl
,
weozl
在 {ox
w’y
w’z
w’}中的坐标为：
[0101][0102][0103][0104]
故：
[0105][0106]
步骤2.3，由于
wadm
(kd)＝k
mw
ω
×wvd(kd)恒成立，所以可以得到式(16)：
[0107]rd12
(kd)
wadmx
(kd)cosθd+r
d22
(kd)
wadmy
(kd)cosθd+
wadmz
(kd)θd＝0
ꢀꢀ
(16)
[0108]
由于在实际中乒乓球旋转主要以带有少量侧旋的上下旋为主，所以cisθd≠0恒成立，上式两边除以cosθd，得到式(17)：
[0109]rd12
(kd)
wadmx
(kd)+r
d22
(kd)
wadmy
(kd)+wad
mz
(kd)tanθd＝0
ꢀꢀ
(17)
[0110]
联立(17)，利用式(18)和式(19)分别计算关于旋转角速度矢量方向角线性方程的左系数矩阵ld和右系数矩阵rd：
[0111]
ld＝[a
dmz
(1) a
dmz
(2) ... a
dmz
(nd)]
t
ꢀꢀ
(18)
[0112][0113]
步骤2.4，利用球的运动速度和由马格努斯力引起的加速度，计算旋转角速度矢量的方向角具体实施方式如下：
[0114]
a、通过(20)计算θd的正切值tanθd：
[0115]
tanθd＝(l
dt
ld)-1
l
dt
rdꢀꢀ
(20)
[0116]
b、通过式(21)计算旋转方向矢量的方向角
[0117][0118]
步骤3，本技术利用步骤1.4中得到的km与||wv||,||wω||的函数关系，基于步骤2求得的旋转角速度矢量方向，通过每个采样时刻的马格努斯力引起的加速度求解旋转角速
度矢量的模。
[0119]
步骤3.1，根据步骤2中得到的方向角用式(22)计算第kd个采样时刻角速度矢量wωd和速度矢量wvd(kd)的夹角θ
dωv
(kd)：
[0120][0121]
式(22)中，
wedω
表示与wωd同向的单位向量，并有：
[0122][0123]
步骤3.2，利用式(24)计算观测变量y2：
[0124]
y2＝[||
wadm
(kd)|||kd＝1，2...nd]
ꢀꢀ
(24)
[0125]
步骤3.3，通过求解函数最小值求解旋转角速度矢量的模：
[0126]
a、令自变量x3＝[||wvd(kd)|||kd＝1,2
…
nd],x4＝[θ
dωv
(kd)|kd＝1,2
…
nd]；
[0127]
b、利用式(25)构建拟合函数y：
[0128]
y＝g(x3(kd)，||ωd||，x4(kd))＝f(x3(kd)，||ωd||)
·
x3(kd)||ωd||sin(x4(kd))
ꢀꢀ
(25)
[0129]
式(25)中，||ωd||为非线性回归的待定系数；
[0130]
c、利用式(26)构建非线性函数并求解最小值，从而得到乒乓球旋转角速度的模c、利用式(26)构建非线性函数并求解最小值，从而得到乒乓球旋转角速度的模
[0131]
步骤3的流程如图4所示。

技术特征：
1.一种基于运动轨迹的乒乓球旋转估计方法，其特征在于，所述方法包括：步骤1、基于乒乓球飞行空气动力学，对若干组乒乓球飞行的轨迹进行求解，得到乒乓球飞行过程中的空气阻力系数、马格努斯力系数、球线速度、由马格努斯力引起的乒乓球的加速度，用于拟合马格努斯力系数与球线速度、角速度的关系；步骤1.1，利用发球机发出角速度为ω的乒乓球，在乒乓球飞行过程中，通过双目视觉以采样间隔t
s
对球位置坐标进行采样，得到乒乓球的轨迹{
w
p(k):k＝1,2
…
n}；其中，k表示采样时刻数，n表示一条轨迹中采样点的个数，
w
p(k)表示第k个采样时刻采集到的乒乓球的位置矢量；且
w
p(k)＝[
w
x(k),
w
y(k),
w
z(k)]，
w
x(k),
w
y(k),
w
z(k)分别表示位置矢量
w
p(k)在世界坐标系的x轴、y轴、z轴上的分量；步骤1.2，通过对
w
x(k),
w
y(k),
w
z(k)进行三次多项式拟合，得到如式(1)、式(2)、式(3)所示的乒乓球的坐标与时间t的连续关系：
w
x(t)＝f
x
(
w
x(k),t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
w
y(t)＝f
x
(
w
y(k),t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
w
z(t)＝f
x
(
w
z(k),t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)式(1)、式(2)、式(3)中，f
x
(
·
,
·
)、f
y
(
·
,
·
)、f
z
(
·
,
·
)分别为世界坐标系的x轴、y轴、z轴上关于时间变量t的三次多项式函数，
w
x(t),
w
y(t),
w
z(t)表示拟合后t时刻乒乓球的位置矢量在世界坐标系的x轴、y轴、z轴的分量；通过对式(1)、式(2)、式(3)求关于t的一阶、二阶导数，并令t＝kt
s
，从而得到如式(4)和式(5)所示的第k个采样时刻球在世界坐标系下的球线速度
w
v(k)和加速度
w
a(k)：a(k)：式(4)和式(5)中，
w
v(k)＝[
w
v
x
(k),
w
v
y
(k),
w
v
z
(k)]，
w
a(k)＝[
w
a
x
(k),
w
a
y
(k),
w
a
z
(k)]，
w
v
x
(k),
w
v
y
(k),
w
v
z
(k)为
w
v(k)在世界坐标系的x轴，y轴，z轴的分量，
w
a
x
(k),
w
a
y
(k),
w
a
z
(k)为
w
a(k)在世界坐标系的x轴，y轴，z轴的分量；步骤1.3，利用式(6)、式(7)、式(8)分别求解第k个采样时刻乒乓球的空气阻力系数k
d
(k)、由马格努斯力产生的加速度在世界坐标系下的加速度
w
a
m
(k)、马格努斯力系数k
m
(k)：(k)：(k)：式(6)、式(7)中，
w
g表示在世界坐标系中的重力加速度矢量；步骤1.4，按照步骤1.1-1.3的过程，采集不同角速度，不同线速度，不同发射方向的n组乒乓球轨迹并得到包含所有的k
m
(k)、
w
ω(k)、
w
v(k)的集合：y＝{{k
m
(k
j
)}
j
}、x1＝{{
w
ω
(k
j
)}
j
}、x2＝{{
w
v(k
j
)}
j
},j＝1,2
…
n,k
j
＝1,2
…
n
j
，其中，j表示轨迹索引，k
j
表示第j个轨迹的坐标点索引，y为包含个马格努斯力系数的集合，x1,x2为包含个角速度、线速度的集合；并将x1，x2作为自变量，y作为自变量进行二维多项式函数拟合，得到如式(9)所示的关系式：k
m
＝f(
w
v,
w
ω)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)步骤2、基于所述马格努斯力系数与球线速度、角速度的拟合关系，利用乒乓球飞行的新轨迹构建角速度矢量方向的线性方程组并进行求解，得到乒乓球旋转角速度矢量的方向角；步骤3、基于所述乒乓球旋转角速度矢量的方向，利用乒乓球飞行的新轨迹构建旋转角速度模的非线性函数并求解其最小值，得到乒乓球旋转角速度的模。2.根据权利要求1中所述的乒乓球旋转估计方法，其特征在于，所述步骤2包括：步骤2.1，给定一条待求的新轨迹t
d
＝{
w
p
d
(k
d
)|k
d
＝1,2
…
n
d
},
w
p
d
(k
d
)表示待求的新轨迹中第k
d
个采样时刻球在世界坐标系下的位置矢量，n
d
表示待求的新轨迹中包含的乒乓球坐标的个数；利用式(4)计算每一个采样时刻球在世界坐标系下的速度矢量
w
v
d
(k)＝[
w
v
dx
(k),
w
v
dy
(k), w
v
dz
(k)]、利用式(7)计算在世界坐标系下由马格努斯力引起的加速度矢量
w
a
dm
(k)＝[
w
a
dmx
(k), w
a
dmy
(k),
w
a
dmz
(k)]，其中，
w
v
dx
(k),
w
v
dy
(k),
w
v
dz
(k)为
w
v
d
(k)在世界坐标系的x轴、y轴、z轴的分量；
w
a
dmx
(k),
w
a
dmy
(k
d
),
w
a
dmz
(k)为
w
a
dm
(k)在世界坐标系的x轴、y轴、z轴的分量；步骤2.2，按照式(10)计算待求的新轨迹中第k
d
个时刻乒乓球的相对坐标系和世界坐标系的转换矩阵
w
r
dl
(k)：式(10)中，r
d11
(k),r
d12
(k),r
d21
(k),r
d22
(k)为矩阵
w
r
dl
(k)中第1行第1列，第1行第2列，第2行第1列，第2行第2列的元素；步骤2.3，利用式(11)和式(12)分别计算关于旋转角速度矢量方向角线性方程的左系数矩阵l
d
和右系数矩阵r
d
：l
d
＝[a
dmz
(1) a
dmz
(2) ... a
dmz
(n
d
)]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)步骤2.4，计算旋转角速度矢量的方向角
a、通过式(13)计算θ
d
的正切值tanθ
d
：tanθ
d
＝(l
dt
l
d
)-1
l
dt
r
d
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)b、通过式(14)计算旋转方向矢量的方向角b、通过式(14)计算旋转方向矢量的方向角3.根据权利要求2中所述的乒乓球旋转估计方法，其特征在于，所述步骤3包括：步骤3.1，根据步骤2中得到的方向角用式(15)计算第k
d
个采样时刻角速度矢量
w
ω
d
和速度矢量
w
v
d
(k
d
)的夹角θ
dωv
(k
d
)：式(15)中，
w
e
dω
表示与
w
ω
d
同向的单位向量，并有：步骤3.2，利用式(17)计算观测变量y2：y2＝[||
w
a
dm
(k
d
)|||k
d
＝1，2...n
d
]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(17)步骤3.3，通过求解函数最小值求解旋转角速度矢量的模：a、令自变量x3＝[||
w
v
d
(k
d
)|||k
d
＝1,2
…
n
d
],x4＝[θ
dωv
(k
d
)|k
d
＝1,2
…
n
d
]；b、利用式(18)构建拟合函数y：y＝g(x3(k
d
)，||ω
d
||，x4(k
d
))＝f(x3(k
d
)，||ω
d
||)
·
x3(k
d
)||ω
d
||sin(x4(k
d
))
ꢀꢀ
(18)式(18)中，||ω
d
||为非线性回归的待定系数；c、利用式(19)构建非线性函数并求解最小值，从而得到乒乓球旋转角速度的模c、利用式(19)构建非线性函数并求解最小值，从而得到乒乓球旋转角速度的模4.一种电子设备，包括存储器以及处理器，其特征在于，所述存储器用于存储支持处理器执行权利要求1或2或3所述乒乓球旋转估计方法的程序，所述处理器被配置为用于执行所述存储器中存储的程序。5.一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，其特征在于，所述计算机程序被处理器运行时执行权利要求1或2或3所述乒乓球旋转估计方法的步骤。

技术总结
本发明公开了一种基于运动轨迹的乒乓球旋转估计方法及其应用，该方法包括：1、拟合马格努斯力系数与球线速度、角速度的关系；2、基于所述马格努斯力系数与球线速度、角速度的拟合关系，利用乒乓球飞行的新轨迹构建角速度矢量方向的线性方程组并进行求解，得到乒乓球旋转角速度矢量的方向角；3、基于所述乒乓球旋转角速度矢量的方向，利用乒乓球飞行的新轨迹构建旋转角速度模的非线性函数并求解其最小值，得到乒乓球旋转角速度的模。本发明能准确地估计球在飞行中的旋转速度，从而能提高乒乓球轨迹预测的精度。迹预测的精度。迹预测的精度。


技术研发人员：宋博 王昱欣 孙智涌 张强 成二康
受保护的技术使用者：中国科学院合肥物质科学研究院
技术研发日：2022.09.26
技术公布日：2023/1/6