CN112244833A - 一种基于协作机械臂的人体上肢多维末端刚度测量方法 - Google Patents

一种基于协作机械臂的人体上肢多维末端刚度测量方法 Download PDF

Info

Publication number
CN112244833A
CN112244833A CN202011032974.XA CN202011032974A CN112244833A CN 112244833 A CN112244833 A CN 112244833A CN 202011032974 A CN202011032974 A CN 202011032974A CN 112244833 A CN112244833 A CN 112244833A
Authority
CN
China
Prior art keywords
disturbance
upper limb
mechanical arm
dimensional
function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011032974.XA
Other languages
English (en)
Other versions
CN112244833B (zh
Inventor
张琴
刘济汉
熊蔡华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN202011032974.XA priority Critical patent/CN112244833B/zh
Publication of CN112244833A publication Critical patent/CN112244833A/zh
Application granted granted Critical
Publication of CN112244833B publication Critical patent/CN112244833B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/22—Ergometry; Measuring muscular strength or the force of a muscular blow
    • A61B5/221—Ergometry, e.g. by using bicycle type apparatus
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/22—Ergometry; Measuring muscular strength or the force of a muscular blow
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7239—Details of waveform analysis using differentiation including higher order derivatives

Abstract

本发明公开了一种基于协作机械臂的人体上肢多维末端刚度测量方法，属于肢体刚度测量领域。首先，产生特定频带范围内的多维随机微扰动，使用多自由度协作机械臂对人体上肢末端施加产生的扰动，并采用多输入多输出系统的参数辨识方法解算多维的末端刚度参数。本发明提出了利用多自由度协作机械臂产生扰动来进行人体上肢末端刚度测量，无需专门开发扰动设备就能产生满足幅值与频带范围要求的随机扰动，省去了开发专用型测量装置的复杂过程，降低了人体上肢末端刚度测量的难度，结合多维随机扰动的产生与多输入多输出系统的参数辨识方法，在机械臂工具端能实现足够多维度的运动的情况下，可以根据测量需求实现一维、二维或三维的人体上肢末端刚度测量。

Description

一种基于协作机械臂的人体上肢多维末端刚度测量方法
技术领域
本发明属于肢体刚度测量领域，更具体的，涉及一种基于协作机械臂的人体上肢多维末端刚度测量方法。
背景技术
人体的肢体刚度描述了肢体抵抗外界干扰的强烈程度，肢体的刚度调节特性对维持肢体与未知环境的稳定交互起着至关重要的作用。由于人体上肢末端执行了人体大部分的复杂交互任务，上肢末端刚度受到广泛的关注。获取上肢末端刚度的方法是对上肢末端施加微小的扰动并采集扰动引发的末端回复力，再通过末端动力学模型解算得到刚度参数。尽管扰动的幅值与频带范围设定会影响刚度测量结果，但为了保证测量的精准度，所施加的扰动需要满足三个基本条件：1)幅度在10mm以内以保证受试者上肢姿态近似不变，2)具有随机性以避免受试者产生自发反应，3)在期望研究的频带范围内有均衡的扰动功率。
传统的方案大多是根据研究需求而专门开发连杆式机电设备来用于末端刚度测量，所需的开发时间长，所产生的扰动特征强烈依赖于测量设备的设计特点而难以同时满足扰动所需条件，且大多只能应用于二维末端刚度测量。
因此，需要开发一种新型的方法或者设备，能提供符合要求的扰动，实现多维的末端刚度测量，并且测量精度较高。
发明内容
针对现有技术的缺陷，本发明的目的在于，提供一种基于协作机械臂的人体上肢多维末端刚度测量方法，利用多自由度协作机械臂对上肢末端施加规划的扰动，最后经过多输入多输出系统参数辨识过程实现了人体上肢多维末端刚度测量。
为实现上述目的，本发明提供了一种基于协作机械臂的人体上肢多维末端刚度测量方法，其包括依次执行的如下步骤：
(1)生成扰动控制信息，以用于发送给机械臂工具端，运动控制信息包括位置序列ds与对应的时间序列ts，
(2)受试者以上肢握住机械臂工具端的握把，按照ts中的时间信息向机械臂发送对应的ds中的目标位置信息，使机械臂执行扰动，收集上肢末端回复力数据f和上肢末端的位移数据x，
(3)根据上肢末端回复力数据f和上肢末端的位移数据x，进行刚度解算，获得人体上肢多维末端刚度值。
进一步的，步骤(1)中，为使产生的扰动能满足频带范围的需求，构建双边的理想幅值谱函数，为满足N维测量需求，分别构建N个服从区间(-π,π)上均匀分布的随机相位谱函数，将理想幅值谱函数与随机相位谱函数相结合，得到期望频谱函数，对期望频谱函数进行反快速傅里叶变换得到对应的时域函数，时域函数为离散的位置序列，对时域函数执行缩放，以满足扰动幅值条件，将缩放后时域函数的关键点作为扰动位置目标，从而获得位置序列ds与对应的时间序列ts。
进一步的，步骤(1)的具体步骤为：
设所使用机械臂的控制频率为fc，产生扰动的总持续时间为ta，总控制点数num＝tafc，num为偶数，
(1.1)构建理想幅值谱函数
设期望的频带范围的上限是fh，下限是fl，构建双边的理想幅值谱函数A(n)如下：
Figure BDA0002704328600000021
其中，n为控制点序号，n＝1,2,…,num，ta为扰动的总持续时间，num为总控制点数，
(1.2)构建随机相位谱函数
设需测量维度为N的末端刚度，分别构建N个服从区间(-π,π)上均匀分布的单边随机相位谱函数
Figure BDA0002704328600000031
如下：
Figure BDA0002704328600000032
其中，i代表维度序号，i＝1,2,…,N，m为控制点序号，m＝1,2,…,num/2，每个维度上的
Figure BDA0002704328600000033
充分随机且相互独立的，保证产生的扰动具有很强的随机性，则双边的随机相位谱函数
Figure BDA0002704328600000034
如下：
Figure BDA0002704328600000035
其中i代表维度序号，n为控制点序号，num为总控制点数，
(1.3)变换为时域位置函数
将理想幅值谱函数与随机相位谱函数相结合，得到N个期望频谱函数Di(n)：
Figure BDA0002704328600000036
其中，n为控制点序号，
Figure BDA0002704328600000037
表示虚数单位，对期望频谱函数进行反快速傅里叶变换得到对应的时域函数di(k)如下：
di(k)＝ifft(Di(n))
di(k)是离散的位置序列，k为位置序号，k＝1,2,…,num，
(1.4)扰动幅值缩放
设定实际施加扰动的最大幅值为dmax，依据位置序列di(k)与扰动中心点的欧氏距离
Figure BDA0002704328600000038
的最大值，对di(k)进行缩放以满足幅值条件，缩放后的位置序列
Figure BDA0002704328600000039
如下：
Figure BDA00027043286000000310
其中k为位置序号，
(1.5)生成扰动控制信息
将缩放后位置序列的关键点作为扰动位置目标，从而获得位置序列ds与对应的时间序列ts。
进一步的，确定距离量ρ(k)中所有局部最大值对应的序号和所有局部最小值对应的序号，进而确定缩放后位置序列的关键点，组成序号集合Id如下：
Id＝{k|ρ(k)≥ρ(k-1)andρ(k)≥ρ(k+1)}∪{k|ρ(k)≤ρ(k-1)andρ(k)≤ρ(k+1)}
得到位置序列
Figure BDA0002704328600000041
与对应的时间序列
Figure BDA0002704328600000042
其中p代表关键点的位置序号，ds与ts共同构成了最终发送给机械臂工具端的扰动控制信息。
进一步的，步骤(2)中，受试者保持固定的上肢姿态，实验过程中机械臂工具端固定的力传感器收集上肢因受到扰动而引发的回复力数据f，外部的光学三维捕捉系统收集上肢末端的位移数据x。
进一步的，上肢末端动力学为以末端位移x(t)＝[x1(t) ... xN(t)]T为输入，以末端回复力f(t)＝[f1(t) … fN(t)]T为输出的多输入多输出的系统，其中t为时间，经过多输入多输出系统参数辨识过程实现人体上肢多维末端刚度解算。
进一步的，刚度解算具体过程为：
构建多维空间内的上肢末端系统方程为：
Figure BDA0002704328600000043
其中，
Figure BDA0002704328600000044
代表末端速度，
Figure BDA0002704328600000045
代表末端加速度，I、B、K分别为末端的惯性、阻尼和刚度矩阵，t为时间，
(3.1)传递函数求解
设f为频率，X(f)＝[X1(f) ... XN(f)]T与F(f)＝[F1(f) ... FN(f)]T分别是末端位移x(t)与末端回复力f(t)的频域表达形式，令传递函数矩阵H如下：
Figure BDA0002704328600000046
其中，
Figure BDA0002704328600000051
是联系输入xi与输出fj的单输入单输出子系统传递函数，i、j代表维度序号，i＝1,2,…,N，j＝1,2,…,N，设
Figure BDA0002704328600000052
为xi与fj的互功率谱，
Figure BDA0002704328600000053
为xi的自功率谱或者
Figure BDA0002704328600000054
为xi与xo的互功率谱，o代表维度序号，o＝1,2,…,N：
Figure BDA0002704328600000055
令Gxf为包含所有
Figure BDA0002704328600000056
的输入输出互功率谱矩阵、Gxx为包含所有
Figure BDA0002704328600000057
的输入自/互功率谱矩阵：
Figure BDA0002704328600000058
计算
Figure BDA0002704328600000059
得到N*N个子系统各自的传递函数
Figure BDA00027043286000000517
(3.2)刚度参数辨识
每个子系统的传递函数可被表达为：
Figure BDA00027043286000000510
其中，s代表复频率，f代表频率，利用最小二乘方法求出使
Figure BDA00027043286000000511
最小的惯性参数
Figure BDA00027043286000000512
阻尼参数
Figure BDA00027043286000000513
与刚度参数
Figure BDA00027043286000000514
组合所有子系统的
Figure BDA00027043286000000515
最终得到多维末端刚度矩阵K如下：
Figure BDA00027043286000000516
通过本发明所构思的以上技术方案，与现有技术相比，能够取得下列有益效果：
1、提出了利用多自由度协作机械臂产生扰动来进行人体上肢末端刚度测量的方法，省去了开发专用型测量装置的复杂过程，降低了人体上肢末端刚度测量的难度。
2、产生的扰动能很好地满足上肢末端刚度测量所需的扰动幅值、频带范围以及随机性要求，易于根据需求对扰动的幅值、频带范围进行调整，且最终产生的运动控制信息能较容易地应用在任何支持工具端位置控制的协作机械臂平台上。
3、刚度解算方法从功率谱角度出发，能针对带有随机性输入的系统求解其传递函数，并在指定频带范围内辨识传递函数中的未知参数，有较高的刚度解算精度。
4、结合了多维随机扰动的产生与多输入多输出系统的参数辨识方法，在机械臂工具端能实现足够多维度的运动的情况下，可以根据测量需求实现一维、二维或三维的人体上肢末端刚度测量。
附图说明
图1是本发明实施例中基于协作机械臂的人体上肢末端刚度测量方法的原理图；
图2是本发明实施例中扰动实验的示意图；
图3是本发明实施例中产生的三维扰动与扰动功率谱，其包括(a)三个维度上扰动的轨迹，(b)三个维度上扰动的功率谱，(c)三维扰动的轨迹距中心点的欧氏距离曲线；
图4是本发明实施例中子系统传递函数拟合效果，其中每个子图都对应一个单输入单输出的末端动力学子系统，各子图的行数与列数分别对应输出与输入的维度数；
图5是本发明实施例中实测回复力与估计回复力局部对比图，其种三个子图分别表达了三个维度上的实测回复力与估计回复力对比；
图6是本发明实施例中测量得到的空间刚度椭球，其中三幅人像图分别为(a)正视图、(b)左视图、(c)俯视图，子图(d)表现了空间刚度椭球中的长轴、次长轴和短轴的长度与空间方向。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明提出了特定频带范围的多维随机微扰动产生方法，并利用多自由度协作机械臂对上肢末端施加规划的扰动，最后经过多输入多输出系统参数辨识过程实现了人体上肢多维末端刚度测量。所发明的方法在机械臂工具端有足够多运动维度的情况下，能对一维、二维或三维的上肢末端刚度进行测量。其中技术方案包括扰动产生、扰动实验、刚度解算三部分，具体细节如下：
(1)扰动产生
假设所使用机械臂的控制频率为fc，产生扰动的总持续时间为ta，则总控制点数num＝tafc，为了计算方便一般保证num为偶数。
(1.1)构建理想幅值谱函数
为了使产生的扰动能满足频带范围的需求，构建一个理想的幅值谱函数。假设期望的频带范围的上限是fh，下限是fl，则构建双边的理想幅值谱函数A(n)如下：
Figure BDA0002704328600000071
其中，n为控制点序号，n＝1,2,…,num，ta为扰动的总持续时间，num为总控制点数。
(1.2)构建随机相位谱函数
假设需要对维度为N(N∈{1,2,3})的末端刚度进行测量，则使用随机数生成器，分别构建N个服从区间(-π,π)上均匀分布的单边随机相位谱函数
Figure BDA0002704328600000072
如下：
Figure BDA0002704328600000073
其中i代表维度序号，i＝1,2,…,N，m为控制点序号，m＝1,2,…,num/2。每个维度上的
Figure BDA0002704328600000074
充分随机且是相互独立的，保证产生的扰动具有很强的随机性。则双边的随机相位谱函数
Figure BDA0002704328600000081
如下：
Figure BDA0002704328600000082
其中，i代表维度序号，n为控制点序号，num为总控制点数。
(1.3)变换为时域位置函数
将理想幅值谱函数与随机相位谱函数相结合，得到N个期望频谱函数Di(n)如下：
Figure BDA0002704328600000083
其中n为控制点序号，
Figure BDA0002704328600000084
表示虚数单位。对期望频谱函数进行反快速傅里叶变换得到对应的时域函数di(k)如下：
di(k)＝ifft(Di(n))
di(k)即是离散的位置序列，k为位置序号，根据快速傅里叶变换不改变数据个数的性质，有k＝1,2,…,num。
(1.4)扰动幅值缩放
设定实际施加扰动的最大幅值为dmax，依据di(k)与扰动中心点(即原点)的欧氏距离
Figure BDA0002704328600000085
的最大值，对di(k)进行缩放以满足幅值条件，缩放后的位置序列
Figure BDA0002704328600000086
如下：
Figure BDA0002704328600000087
其中k为位置序号，
(1.5)生成运动控制信息
每个
Figure BDA0002704328600000088
都是一个离散位置序列，带有很强的随机性，但若直接将
Figure BDA0002704328600000089
设定为机械臂末端的时序位置目标，则会在每个控制点都会给机械臂发送运动指令，这将导致机械臂运动速度与加速度的频繁变化，进而可能引起机械臂保护性停止或运行故障。为了避免上述问题的发生，采用
Figure BDA00027043286000000810
中的某些关键路径点来替代完整的
Figure BDA00027043286000000811
作为运动位置目标，具体方法为找到距离量ρ(k)中所有局部最大值对应的序号和所有局部最小值对应的序号，组成序号集合Id如下：
Id＝{k|ρ(k)≥ρ(k-1)andρ(k)≥ρ(k+1)}∪{k|ρ(k)≤ρ(k-1)andρ(k)≤ρ(k+1)}
进而得到位置序列
Figure BDA0002704328600000091
与对应的时间序列
Figure BDA0002704328600000092
其中p代表关键点的位置序号，ds与ts共同构成了最终发送给机械臂工具端的运动控制信息。
(2)扰动实验
扰动实验中，先使受试者握住机械臂工具端所固定的握把，随后按照时间序列ts中的时间信息向机械臂发送对应的位置序列ds中的目标位置信息，使机械臂执行扰动。整个过程中受试者保持固定的上肢姿态。实验过程中机械臂工具端固定的力传感器收集上肢末端因收到扰动而引发的回复力数据f，外部的光学三维捕捉系统收集上肢末端的位移数据X。
(3)刚度解算
多维空间内的上肢末端动力学是一个以末端位移x(t)＝[x1(t) ... xN(t)]T为输入，以末端回复力f(t)＝[f1(t) … fN(t)]T为输出的多输入多输出的系统，构建其系统方程为：
Figure BDA0002704328600000093
其中t为时间，
Figure BDA0002704328600000094
代表末端速度，
Figure BDA0002704328600000095
代表末端加速度，I、B、K分别为末端的惯性、阻尼和刚度矩阵。由于其系统输入是从频域角度产生的，因此也从频域角度来求解其中的未知参数，具体细节如下：
(3.1)传递函数求解
设f为频率，X(f)＝[X1(f) ... XN(f)]T与F(f)＝[F1(f) ... FN(f)]T分别是末端位移x(t)与末端回复力f(t)的频域表达形式，令传递函数矩阵H如下：
Figure BDA0002704328600000096
其中
Figure BDA0002704328600000097
是联系输入xi与输出fj的单输入单输出子系统传递函数，i、j代表维度序号，i＝1,2,…,N，j＝1,2,…,N，。则有：
F＝HX
由于该系统输入信号有强随机性，因此其没有确定的频谱，故系统传递函数需要根据功率谱来求解。设
Figure BDA0002704328600000101
为xi与fj的互功率谱，
Figure BDA0002704328600000102
为xi的自功率谱(i＝o时)或者xi与xo的互功率谱(i≠o时)，其中o代表维度序号，o＝1,2,…,N，则有：
Figure BDA0002704328600000103
令Gxf为包含所有
Figure BDA0002704328600000104
的输入输出互功率谱矩阵、Gxx为包含所有
Figure BDA0002704328600000105
的输入自/互功率矩阵：
Figure BDA0002704328600000106
则有：
Gxf＝GxxHT
因此根据
Figure BDA0002704328600000107
即可以得到N*N个子系统各自的传递函数
Figure BDA0002704328600000108
(3.2)刚度参数辨识
每个子系统的传递函数可被表达为：
Figure BDA0002704328600000109
其中s代表复频率，f代表频率。利用最小二乘方法求出使
Figure BDA00027043286000001010
最小的惯性参数
Figure BDA00027043286000001011
阻尼参数
Figure BDA00027043286000001012
与刚度参数
Figure BDA00027043286000001013
组合所有子系统的
Figure BDA00027043286000001014
即可最终得到多维末端刚度矩阵K如下：
Figure BDA00027043286000001015
下面以更为具体的实施例进一步阐述本发明方法，本发明测量实例是使用一种6自由度的协作型机械臂来对人体上肢末端施加三维随机扰动，根据扰动过程中收集的末端位移与回复力数据，计算得到三维的上肢末端刚度。整体测量流程如图1，图1为人体上肢末端刚度测量流程图，由图可知，其主要的流程分为三大步骤，扰动产生、扰动实验和刚度解算。扰动产生包括理想幅值谱函数与随机相位谱函数的构建、反快速傅里叶变换、幅值缩放和生成运动控制信息等过程。扰动实验包括使用协作机械臂施加扰动和收集末端位置、回复力信息等过程。刚度解算包括计算子系统传递函数和刚度参数辨识等过程。
本发明的实施例具体包含以下步骤：
(1)扰动的产生
所使用的协作型机械臂的控制频率为fc＝125Hz，计划产生扰动的总持续时间为ta＝20s，则总控制点数num＝tafc。
(1.1)构建理想幅值谱函数
设定期望研究的频带范围的上限是fh＝5Hz，下限是fl＝0Hz，则根据构建双边的理想幅值谱函数A(n)如下：
Figure BDA0002704328600000111
其中n为控制点序号，n＝1,2,…,num，ta为扰动的总持续时间，num为总控制点数。
(1.2)构建随机相位谱函数
测量维度N＝3，分别构建3个服从区间(-π,π)上均匀分布的单边随机相位谱函数：
Figure BDA0002704328600000112
其中i代表维度序号，i＝1,2,3，m为控制点序号，m＝1,2,…,num/2。则构建双边的随机相位谱函数
Figure BDA0002704328600000113
如下：
Figure BDA0002704328600000114
其中，i代表维度序号，n为控制点序号，num为总控制点数。
(1.3)变换为时域位置函数
将理想幅值谱函数与随机相位谱函数相结合，得到3个期望频谱函数Di(n)如下：
Figure BDA0002704328600000121
其中n为控制点序号，
Figure BDA0002704328600000122
表示虚数单位。对期望频谱函数进行反快速傅里叶变换得到对应的时域函数di(k)如下：
di(k)＝ifft(Di(n))
其中k为位置序号，k＝1,2,…,num，di(k)即是初始的离散位置序列。
(1.4)扰动幅值缩放
计算di(k)与扰动中心点(即原点)的欧氏距离
Figure BDA0002704328600000123
设定实际施加扰动的最大幅值为dmax＝10mm，则对di(k)进行缩放以满足幅值条件，缩放后的位置序列
Figure BDA0002704328600000124
如下：
Figure BDA0002704328600000125
其中k为位置序号，
(1.5)生成运动控制信息
若直接将
Figure BDA0002704328600000126
设定为机械臂末端的时序位置目标可能引起机械臂保护性停止或运行故障，因此挑选
Figure BDA0002704328600000127
中的关键路径点来代替完整
Figure BDA0002704328600000128
具体方法为找到距离量ρ(k)的所有局部最大值对应的序号和所有局部最小值对应的序号，组成序号集合Id如下：
Id＝{k|ρ(k)≥ρ(k-1)andρ(k)≥ρ(k+1)}∪{k|ρ(k)≤ρ(k-1)andρ(k)≤ρ(k+1)}
进而得到最终发送给机械臂的运动控制信息，包括末端位置序列
Figure BDA0002704328600000129
与对应的发送时间序列
Figure BDA00027043286000001210
其中p代表关键点的位置序号。
(2)扰动实验
图2为扰动实验示意图，由图可知，受试者坐在一硬质椅子上，6自由度的协作型机械臂的握把位于受试者前方，受试者握住机械臂工具端所固定的握把。实验中，按照ts中的时间信息向机械臂发送对应的ds中的目标位置信息，使机械臂执行扰动。在此过程中，受试者保持固定的上肢姿态。实验过程中机械臂工具端固定的力传感器收集上肢因收到扰动而引发的回复力数据f，外部的光学三维捕捉系统收集上肢末端的位移数据X。实际产生的三维扰动及其功率谱如图3，图3中(a)给出了X、Y、Z三个维度上扰动的轨迹,(b)给出了X、Y、Z三个维度上扰动的功率谱，(c)给出了三维扰动的轨迹距中心点的欧氏距离曲线，由图可知，扰动具有很强的随机性，其功率谱在5Hz以下基本保持平坦，三维扰动幅值在10mm以内，因此满足设计需求。
(3)刚度解算
三维空间内的上肢末端动力学是一个以末端位移x(t)＝[x1(t) … x3(t)]T为输入，以末端回复力f(t)＝[f1(t) … f3(t)]T为输出的多输入多输出的系统，构建其系统方程为：
Figure BDA0002704328600000131
其中t为时间，
Figure BDA0002704328600000132
代表末端速度，
Figure BDA0002704328600000133
代表末端加速度，I、B、K分别为末端的惯性、阻尼和刚度参数，解算I、B、K的具体细节如下：
(3.1)传递函数求解
计算xi与fj的互功率谱
Figure BDA0002704328600000134
计算xi的自功率谱(i＝o时)或者xi与xo的互功率谱(i≠o时)
Figure BDA0002704328600000135
其中i、j、o都代表维度序号，o＝1,2,3，j＝1,2,3，令Gxf为包含所有
Figure BDA0002704328600000136
的输入输出互功率谱矩阵、Gxx为包含所有
Figure BDA0002704328600000137
的输入自/互功率谱矩阵：
Figure BDA0002704328600000138
则根据下面的解算方程：
Figure BDA0002704328600000141
即可以得到3*3＝9个子系统各自的传递函数
Figure BDA0002704328600000142
(3.2)刚度参数辨识
每个子系统的传递函数可被表达为：
Figure BDA0002704328600000143
其中，s代表复频率，f代表频率。利用最小二乘方法在0～10Hz范围内辨识出使
Figure BDA0002704328600000144
最小的惯性参数
Figure BDA0002704328600000145
阻尼参数
Figure BDA0002704328600000146
与刚度参数
Figure BDA0002704328600000147
为了评估辨识的效果，图4展示了对各子系统传递函数
Figure BDA0002704328600000148
的拟合效果，平均决定系数R2＝0.962，图4中每个子图都对应一个单输入单输出的末端动力学子系统，9个子图的横坐标都为频率且刻度范围相同，纵坐标都为幅值且刻度范围相同，因此采用简化的表达方法省略了部分子图的坐标刻度信息，由图可知，各子系统的拟合传递函数与实测传递函数基本重合，展现了良好的拟合效果。图5展示对实际测得的局部末端回复力f的估计效果，平均决定系数R2＝0.825，图5中三个子图分别对应X、Y、Z三个维度，横坐标都为时间且刻度范围相同，因此采用简化的表达方法省略了部分子图的横坐标刻度信息，由图可知，三个维度上估计回复力与实测回复力都非常接近，展现了较高的刚度解算精度。组合所有子系统的
Figure BDA0002704328600000149
即可得到三维末端刚度矩阵K如下：
Figure BDA00027043286000001410
三维末端刚度矩阵K可用一空间椭球来表达其各向异性，椭球面上任意一点距球心的距离即代表该点与球心连线方向上的刚度大小，故椭球长轴的方向即是末端刚度最大的方向，椭球长短轴的方向即是末端刚度最小的方向。椭球长轴、次长轴、短轴的长度与空间方向信息由刚度矩阵K的奇异值分解得到，最终得到的刚度椭球如图6，图6中分别从正视、左视、俯视角度展现了三维末端刚度椭球的平面投影，并从三维空间视角展现出其长轴、次长轴、短轴的长度与空间方向，充分表现出实测三维末端刚度的空间各向异性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，其包括依次执行的如下步骤：
(1)生成扰动控制信息，以用于发送给机械臂工具端，运动控制信息包括位置序列ds与对应的时间序列ts，
(2)受试者以上肢握住机械臂工具端的握把，按照ts中的时间信息向机械臂发送对应的ds中的目标位置信息，使机械臂执行扰动，收集上肢末端回复力数据f和上肢末端的位移数据x，
(3)根据上肢末端回复力数据f和上肢末端的位移数据x，进行刚度解算，获得人体上肢多维的末端刚度值。
2.如权利要求1所述的一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，步骤(1)中，为使产生的扰动能满足频带范围的需求，构建双边的理想幅值谱函数，为满足N维测量需求，分别构建N个服从区间(-π,π)上均匀分布的随机相位谱函数，将理想幅值谱函数与随机相位谱函数相结合，得到期望频谱函数，对期望频谱函数进行反快速傅里叶变换得到对应的时域函数，时域函数为离散的位置序列，对时域函数执行缩放，以满足扰动幅值条件，将缩放后时域函数的关键点作为扰动位置目标，从而获得位置序列ds与对应的时间序列ts。
3.如权利要求2所述的一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，步骤(1)的具体步骤为：
设所使用机械臂的控制频率为fc，产生扰动的总持续时间为ta，总控制点数num＝tafc，num为偶数，
(1.1)构建理想幅值谱函数
设期望的频带范围的上限是fh，下限是fl，构建双边的理想幅值谱函数A(n)如下：
Figure FDA0002704328590000021
其中，n为控制点序号，n＝1,2,…,num，ta为扰动的总持续时间，num为总控制点数，
(1.2)构建随机相位谱函数
设需测量维度为N的末端刚度，分别构建N个服从区间(-π,π)上均匀分布的单边随机相位谱函数
Figure FDA0002704328590000027
如下：
Figure FDA0002704328590000028
其中，i代表维度序号，i＝1,2,…,N，m为控制点序号，m＝1,2,…,num/2，每个维度上的
Figure FDA0002704328590000029
充分随机且相互独立的，保证产生的扰动具有很强的随机性，则双边的随机相位谱函数
Figure FDA00027043285900000210
如下：
Figure FDA0002704328590000022
其中i代表维度序号，n为控制点序号，num为总控制点数，
(1.3)变换为时域位置函数
将理想幅值谱函数与随机相位谱函数相结合，得到N个期望频谱函数Di(n)：
Figure FDA0002704328590000023
其中，n为控制点序号，
Figure FDA0002704328590000024
表示虚数单位，对期望频谱函数进行反快速傅里叶变换得到对应的时域函数di(k)如下：
di(k)＝ifft(Di(n))
di(k)是离散的位置序列，k为位置序号，k＝1,2,…,num，
(1.4)扰动幅值缩放
设定实际施加扰动的最大幅值为dmax，依据位置序列di(k)与扰动中心点的欧氏距离
Figure FDA0002704328590000025
的最大值，对di(k)进行缩放以满足幅值条件，缩放后的位置序列
Figure FDA0002704328590000026
如下：
Figure FDA0002704328590000031
其中k为位置序号，
(1.5)生成扰动控制信息
将缩放后位置序列的关键点作为扰动位置目标，从而获得位置序列ds与对应的时间序列ts。
4.如权利要求3所述的一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，
确定距离量ρ(k)中所有局部最大值对应的序号和所有局部最小值对应的序号，进而确定缩放后位置序列的关键点，组成序号集合Id如下：
Id＝{k|ρ(k)≥ρ(k-1)andρ(k)≥ρ(k+1)}∪{k|ρ(k)≤ρ(k-1)andρ(k)≤ρ(k+1)}
得到位置序列
Figure FDA0002704328590000032
与对应的时间序列
Figure FDA0002704328590000033
其中p代表关键点的位置序号，ds与ts共同构成了最终发送给机械臂工具端的扰动控制信息。
5.如权利要求4所述的一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，
步骤(2)中，受试者保持固定的上肢姿态，实验过程中机械臂工具端固定的力传感器收集上肢因受到扰动而引发的回复力数据f，外部的光学三维捕捉系统收集上肢末端的位移数据x。
6.如权利要求5所述的一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，上肢末端动力学为以末端位移x(t)＝[x1(t) … xN(t)]T为输入，以末端回复力f(t)＝[f1(t) … fN(t)]T为输出的多输入多输出的系统，其中t为时间，经过多输入多输出系统参数辨识过程实现人体上肢多维末端刚度解算。
7.如权利要求6所述的一种基于协作机械臂的人体上肢多维末端刚度测量方法，其特征在于，刚度解算具体过程为：
构建多维空间内的上肢末端系统方程为：
Figure FDA0002704328590000041
其中，
Figure FDA0002704328590000042
代表末端速度，
Figure FDA0002704328590000043
代表末端加速度，I、B、K分别为末端的惯性、阻尼和刚度矩阵，t为时间，
(3.1)传递函数求解
设f为频率，X(f)＝[X1(f) … XN(f)]T与F(f)＝[F1(f) … FN(f)]T分别是末端位移x(t)与末端回复力f(t)的频域表达形式，令传递函数矩阵H如下：
Figure FDA0002704328590000044
其中，
Figure FDA0002704328590000045
是联系输入xi与输出fj的单输入单输出子系统传递函数，i、j代表维度序号，i＝1,2,…,N，j＝1,2,…,N，设
Figure FDA0002704328590000046
为xi与fj的互功率谱，
Figure FDA0002704328590000047
为xi的自功率谱或者
Figure FDA0002704328590000048
为xi与xo的互功率谱，o代表维度序号，o＝1,2,…,N：
Figure FDA0002704328590000049
令Gxf为包含所有
Figure FDA00027043285900000410
的输入输出互功率谱矩阵、Gxx为包含所有
Figure FDA00027043285900000411
的输入自/互功率谱矩阵：
Figure FDA00027043285900000412
计算
Figure FDA00027043285900000413
得到N*N个子系统各自的传递函数
Figure FDA00027043285900000414
(3.2)刚度参数辨识
每个子系统的传递函数可被表达为：
Figure FDA00027043285900000415
其中，s代表复频率，f代表频率，利用最小二乘方法求出使
Figure FDA00027043285900000416
最小的惯性参数
Figure FDA00027043285900000417
阻尼参数
Figure FDA00027043285900000418
与刚度参数
Figure FDA00027043285900000419
组合所有子系统的
Figure FDA00027043285900000420
最终得到多维末端刚度矩阵K如下：
Figure FDA0002704328590000051
CN202011032974.XA 2020-09-27 2020-09-27 一种基于协作机械臂的人体上肢多维末端刚度测量方法 Active CN112244833B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011032974.XA CN112244833B (zh) 2020-09-27 2020-09-27 一种基于协作机械臂的人体上肢多维末端刚度测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011032974.XA CN112244833B (zh) 2020-09-27 2020-09-27 一种基于协作机械臂的人体上肢多维末端刚度测量方法

Publications (2)

Publication Number Publication Date
CN112244833A true CN112244833A (zh) 2021-01-22
CN112244833B CN112244833B (zh) 2021-09-07

Family

ID=74233289

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011032974.XA Active CN112244833B (zh) 2020-09-27 2020-09-27 一种基于协作机械臂的人体上肢多维末端刚度测量方法

Country Status (1)

Country Link
CN (1) CN112244833B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113084812A (zh) * 2021-04-09 2021-07-09 吉林大学 一种机器人末端刚度性能评价方法
CN113524235A (zh) * 2021-07-15 2021-10-22 华中科技大学 一种面向动态交互的对象肢体时变刚度辨识方法及装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130317627A1 (en) * 2012-05-25 2013-11-28 Kcf Technologies, Inc. Apparatuses and methods for harvesting energy from prosthetic limbs
CN108387351A (zh) * 2018-02-12 2018-08-10 华中科技大学 一种肢体末端刚度测量装置及其测量方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130317627A1 (en) * 2012-05-25 2013-11-28 Kcf Technologies, Inc. Apparatuses and methods for harvesting energy from prosthetic limbs
CN108387351A (zh) * 2018-02-12 2018-08-10 华中科技大学 一种肢体末端刚度测量装置及其测量方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
PANAGIOTIS K等: "Human Arm Impedance: Characterization and Modeling in 3D Space", 《IEEE》 *
韦成朋等: "人体上肢在推拉运动中的末端机械阻抗研究", 《工业控制计算机》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113084812A (zh) * 2021-04-09 2021-07-09 吉林大学 一种机器人末端刚度性能评价方法
CN113084812B (zh) * 2021-04-09 2022-06-21 吉林大学 一种机器人末端刚度性能评价方法
CN113524235A (zh) * 2021-07-15 2021-10-22 华中科技大学 一种面向动态交互的对象肢体时变刚度辨识方法及装置
CN113524235B (zh) * 2021-07-15 2022-03-11 华中科技大学 一种面向动态交互的对象肢体时变刚度辨识方法及装置

Also Published As

Publication number Publication date
CN112244833B (zh) 2021-09-07

Similar Documents

Publication Publication Date Title
CN112244833B (zh) 一种基于协作机械臂的人体上肢多维末端刚度测量方法
CN110561421B (zh) 机械臂间接拖动示教方法及装置
Hussein A review on vision-based control of flexible manipulators
US9092737B2 (en) Systems, methods, and apparatus for 3-D surface mapping, compliance mapping, and spatial registration with an array of cantilevered tactile hair or whisker sensors
CN108594660B (zh) 一种时不变结构的工作模态参数识别方法和系统
Tan et al. Realtime simulation of thin-shell deformable materials using CNN-based mesh embedding
Hirota et al. Interaction with virtual object using deformable hand
Ruppel et al. Simulation of the SynTouch BioTac sensor
Chen et al. AI based gravity compensation algorithm and simulation of load end of robotic arm wrist force
CN112926152A (zh) 一种数字孪生驱动的薄壁件装夹力精准控制与优化方法
Zhai et al. An intelligent control system for robot massaging with uncertain skin characteristics
Du et al. Natural human–machine interface with gesture tracking and cartesian platform for contactless electromagnetic force feedback
Jin et al. Motion prediction of human wearing powered exoskeleton
Besari et al. Feature-based egocentric grasp pose classification for expanding human-object interactions
Yang et al. Model-free 3-d shape control of deformable objects using novel features based on modal analysis
CN110900608B (zh) 基于最优测量构型选择的机器人运动学标定方法
Messay et al. Gpgpu acceleration of a novel calibration method for industrial robots
CN113524235B (zh) 一种面向动态交互的对象肢体时变刚度辨识方法及装置
EP3536467A1 (en) Action transfer device, action transfer method, and non-temporary computer readable medium having action transfer program stored thereon
CN111002292B (zh) 基于相似性度量的机械臂仿人运动示教方法
Vertechy et al. Real-time direct position analysis of parallel spherical wrists by using extra sensors
Meng et al. Spring-IMU Fusion-Based Proprioception for Feedback Control of Soft Manipulators
JP5083992B1 (ja) 把持姿勢生成装置、把持姿勢生成方法及び把持姿勢生成プログラム
Yang et al. Analysis of kinematic parameter identification method based on genetic algorithm
Dash et al. A Inverse Kinematic Solution Of A 6-DOF Industrial Robot Using ANN

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant