Rehabilitation robot motion state digital twinborn monitoring and dynamic compensation system

Abstract

The invention belongs to the technical field of digital twinning, and discloses a digital twinning monitoring and dynamic compensation system for the motion state of a rehabilitation robot. The system obtains a limb surface deformation image and joint torque data of a patient through a multi-view depth camera, constructs a dynamic deformation field, analyzes and screens a high-stability deformation area based on deformation direction consistency and amplitude fluctuation, extracts deformation features and joint torque data to construct dynamic mechanical boundary conditions, and obtains a dynamic mechanical model. And updating the digital twin model in real time by adopting a low-order finite element grid. According to the system, deviation vectors are calculated based on a digital twin model, a dynamic compensation instruction is generated in combination with impedance control parameters, and precise adjustment of the motion trail of the rehabilitation robot is achieved. According to the invention, the robot can make an instant response to the state change of the patient, and the safety and effectiveness of rehabilitation training are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the field of digital twin technology, and more specifically, to a digital twin monitoring and dynamic compensation system for the motion state of a rehabilitation robot. Background Art

[0002] As an important means of assisting patients in recovering their motor function, rehabilitation robots have been widely used in recent years in fields such as neurological rehabilitation and orthopedic rehabilitation. With the deep integration of robotics and rehabilitation medicine, rehabilitation robots have evolved from early passive assistive devices to interactive systems with a degree of intelligence. Currently, rehabilitation robots in clinical use primarily include upper limb exoskeletons, lower limb gait training robots, and hand function rehabilitation robots. These robots help patients recover their motor function by providing precise motion guidance and appropriate mechanical assistance.

[0003] During actual rehabilitation training, soft tissues such as muscles and skin in a patient's limbs undergo highly nonlinear and individualized deformations with movement, exhibiting complex mechanical properties such as strain rate dependence, anisotropy, and viscoelasticity. These deformations not only affect the distribution of contact forces applied by the patient on the robot but also alter the kinematic and dynamic properties of the limb. However, existing technologies still primarily rely on rigid-body assumptions or simplified elastic models, failing to accurately capture these soft tissue deformations. When a patient performs elbow flexion and extension exercises, the contraction and relaxation of the forearm muscles cause significant changes in local tissue stiffness. Unable to monitor these changes in real time, traditional systems often apply inappropriate assistive forces, leading to patient discomfort and even secondary injury. Furthermore, existing digital twin modeling approaches either rely on overly simplified mass-spring models, resulting in insufficient accuracy, or employ high-order nonlinear finite element models, resulting in excessive computational complexity. Solving a typical detailed full-limb model takes minutes or even hours, far from meeting the millisecond-level response requirements of real-time rehabilitation control. In stroke rehabilitation, in particular, muscle tone and spasticity can change instantly, and traditional models lack the ability to provide real-time dynamic boundary condition updates, preventing the robot from adjusting its assistive strategy in a timely manner. This technical defect not only reduces the rehabilitation effect, but may also aggravate the patient's psychological resistance.

[0004] In view of this, the present invention proposes a digital twin monitoring and dynamic compensation system for the motion state of a rehabilitation robot to solve the above problems. Summary of the Invention

[0005] In order to overcome the above-mentioned defects of the prior art and achieve the above-mentioned objectives, the present invention provides the following technical solution: a rehabilitation robot motion state digital twin monitoring and dynamic compensation system, comprising: An acquisition module is used to obtain the surface deformation image and joint torque data of the patient's limb at each moment during the movement of the rehabilitation robot; a deformation field construction module, configured to construct a dynamic deformation field of the patient's limb soft tissue based on the deformation gradient distribution and local curvature change between pixels in the surface deformation image; a stability analysis module for obtaining the deformation stability of each local area in the dynamic deformation field at each moment based on the consistency of deformation direction and deformation amplitude fluctuation of the same local area in the dynamic deformation field at adjacent moments within the analysis period of each moment; and screening out high-stability deformation areas based on the deformation stability; A boundary condition module is used to construct the dynamic mechanical boundary conditions of the digital twin model based on the deformation characteristics of each high-stability deformation area in the dynamic deformation field at each moment and the corresponding joint torque data; A model updating and compensation module is used to update the motion state of the digital twin model in real time based on the dynamic mechanical boundary conditions and a preset low-order finite element grid, and to generate dynamic compensation instructions for the motion trajectory of the rehabilitation robot according to the motion state; The modules are connected via wired and / or wireless means to achieve data transmission between modules.

[0006] The technical effects and advantages of the rehabilitation robot motion state digital twin monitoring and dynamic compensation system of the present invention are as follows: The present invention achieves precise adaptation to individual differences in patients, transforming rehabilitation training from a "one-size-fits-all" model to truly personalized precision intervention, which not only improves the training effect, but also significantly shortens the rehabilitation cycle. During functional training, abnormal muscle tone and spasticity can be sensed in real time and trigger corresponding adjustments, effectively avoiding the risk of secondary injury that may be caused by traditional systems, and further improving patient training compliance. The system's real-time response capability enables the rehabilitation robot to exhibit "flexibility" similar to that of a human therapist, and can intelligently adjust the auxiliary strength and training difficulty when the patient is tired or uncomfortable, creating a safe and comfortable rehabilitation environment. Especially for patients with complex pathological conditions, such as those with dystonia or involuntary movements, the system can adapt to their highly nonlinear and unpredictable movement characteristics, reducing the risk of mechanical injuries common in traditional rehabilitation. In addition, the precise digital twin model generated by the system provides doctors with quantitative evaluation and prognosis prediction tools, which makes up for the lack of subjectivity in traditional evaluation methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0007] Figure 1 Schematic diagram of the digital twin monitoring and dynamic compensation system for the rehabilitation robot's motion state of the present invention. DETAILED DESCRIPTION

[0008] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0009] See also Figure 1 The present invention provides a rehabilitation robot motion state digital twin monitoring and dynamic compensation system, including: An acquisition module is used to obtain the surface deformation image and joint torque data of the patient's limb at each moment during the movement of the rehabilitation robot; The deformation field construction module is used to construct the dynamic deformation field of the patient's limb soft tissue based on the deformation gradient distribution and local curvature changes between pixels in the surface deformation image; The stability analysis module is used to obtain the deformation stability of each local area in the dynamic deformation field at each moment based on the consistency of deformation direction and deformation amplitude fluctuation of the same local area in the dynamic deformation field at adjacent moments within the analysis period of each moment; and to screen out high-stability deformation areas based on deformation stability; The boundary condition module is used to construct the dynamic mechanical boundary conditions of the digital twin model based on the deformation characteristics of each high-stability deformation area in the dynamic deformation field at each moment and its corresponding joint torque data; The model update and compensation module is used to update the motion state of the digital twin model in real time based on dynamic mechanical boundary conditions and preset low-order finite element meshes, and generate dynamic compensation instructions for the motion trajectory of the rehabilitation robot based on the motion state.

[0010] The modules are connected via wired and / or wireless means to achieve data transmission between modules.

[0011] The present invention captures the surface deformation images and joint torque data of the patient's limbs in real time, constructs a dynamic deformation field to analyze the deformation characteristics of the limb soft tissue, screens high-stability deformation areas as reliable observation points, constructs dynamic mechanical boundary conditions to drive the update of the digital twin model, and generates dynamic compensation instructions for the rehabilitation robot based on the motion state of the digital twin model, thereby realizing precise adjustment of the robot's motion trajectory during the rehabilitation process, effectively adapting to individual differences among patients and dynamic changes during the rehabilitation process, and improving the safety and effectiveness of rehabilitation training.

[0012] In an embodiment of the present invention, obtaining the deformation stability of each local area in the dynamic deformation field at each moment includes: Track the deformation vectors of the pixels in the surface deformation image at all times, and obtain the deformation vector of each pixel in the surface deformation image at each moment; Count the number of times the direction angle of the deformation vectors of the same local area in the dynamic deformation field at adjacent moments within the analysis period at each moment is less than the preset direction threshold, and record it as the deformation direction consistency of each local area in the dynamic deformation field at each moment; Calculate the standard deviation of the amplitude of the deformation vector of the same local area in the dynamic deformation field at all moments within the analysis period at each moment as the deformation amplitude fluctuation of each local area in the dynamic deformation field at each moment; The deformation direction consistency and the inverse of the deformation amplitude fluctuation are weightedly fused to obtain the deformation stability of each local area in the dynamic deformation field at each moment.

[0013] In the embodiment of the present invention, the "same local area" refers to a rectangular or circular area of ​​a preset size centered on a pixel point in the surface deformation image, and the specific size is determined according to the actual application requirements. For example, in this embodiment, the local area can be defined as a rectangular area with a side length of 10×10 pixels centered on a certain pixel point, or a circular area with a radius of 5 pixels. The same local area maintains a consistent spatial position in the dynamic deformation field at consecutive moments, that is, it corresponds through the spatial coordinates of the pixel points. The division of the local area covers the entire surface deformation image, ensuring that each pixel point belongs to at least one local area.

[0014] In this embodiment, the surface deformation images of the continuous time series are first analyzed by pixel tracking technology, the temporal motion trajectory of each pixel is calculated, and the deformation vector field at each moment is generated; the analysis period is usually set to 0.5-2 seconds, and the directional change characteristics of the deformation vectors of the same local area are counted within this time window; the directional angle between the deformation vectors of the same local area at adjacent moments is calculated, and a preset directional threshold (usually 15°-30°) is set; the number of times the directional angle is less than the preset threshold is counted, and the higher the value, the more stable the deformation direction of the area, which is recorded as Deformation direction consistency; at the same time, the standard deviation of the deformation vector amplitude of the local area during the analysis period is calculated. The smaller the value, the smaller the deformation amplitude fluctuation, which is recorded as deformation amplitude fluctuation. The deformation direction consistency and the inverse of the deformation amplitude fluctuation are combined through weighted fusion. The typical weighting coefficients are α=0.6 and β=0.4. The calculation formula is: deformation stability = α×deformation direction consistency + β×(1 / deformation amplitude fluctuation). The higher the deformation stability value, the more stable the deformation characteristics of the local area in continuous time, and the more suitable it is as a reliable observation point for subsequent analysis.

[0015] To further clarify the "deformation vector of the same local area" and its relationship to the deformation vector of the pixel point, in this embodiment, the deformation vector of the same local area refers to the set of deformation vectors of all pixels in the local area. Specifically, for a local area, after obtaining the deformation vector of each pixel point in the area, the direction angle and amplitude standard deviation of these deformation vectors are calculated. The direction angle is calculated as follows: for adjacent time points t and t+1, the deformation vectors of the i-th pixel point in the local area are v_i(t) and v_i(t+1) respectively; Calculate the angle between two vectors , where · represents the vector dot product and |·| represents the vector modulus. s is the adjustment coefficient, a very small real number that ensures the denominator is non-zero. The number of times the included angle of all pixels in the local region is less than the preset directional threshold is counted as the deformation direction consistency. The amplitude standard deviation is calculated as follows: for all moments in the analysis period, the deformation vector amplitude of the i-th pixel in the local region is |v_i(t)|. The standard deviation of the deformation vector amplitudes of all pixels in the local region at all moments is calculated as the deformation amplitude fluctuation.

[0016] This implementation achieves precise identification of stable regions within a dynamic deformation field by comprehensively considering both the directional and amplitude stability of deformation. Compared to traditional methods that focus solely on deformation amplitude or analyze a single moment in time, this method captures the dynamic characteristics of deformation through time series analysis, improving the accuracy and robustness of deformation stability assessment and providing a reliable foundation for subsequent screening of highly stable deformation regions.

[0017] In an embodiment of the present invention, a method for screening a high-stability deformation region includes: The deformation stability of each local area in the dynamic deformation field at each moment is sorted, and the local area with deformation stability greater than a preset stability threshold is selected as the candidate area; The mean boundary curvature of each candidate region in the surface deformation image is calculated, and candidate regions with a mean boundary curvature greater than a preset curvature threshold are eliminated to obtain highly stable deformation regions.

[0018] In this embodiment, first, all local regions in the dynamic deformation field at each moment are sorted from high to low according to the deformation stability value; a preset stability threshold is set (typical value is 0.75-0.85, and the full score is 1), and local regions with deformation stability higher than the threshold are selected as preliminary candidate regions; boundary feature analysis is performed on each candidate region to extract its boundary contour in the surface deformation image; the local curvature distribution of the boundary contour is calculated, including calculating the local curvature of each point on the boundary and obtaining the mean curvature of the entire boundary; a preset curvature threshold is set (usually 0.08-0.15mm -1), and eliminate candidate areas with a mean boundary curvature greater than the threshold, because high curvature boundaries usually indicate complex or irregular regional shapes, which may lead to unstable calculations of mechanical properties; after double screening, the retained areas have both high deformation stability and good boundary shape characteristics, and are identified as high-stability deformation areas; these high-stability deformation areas will serve as key observation points for constructing the dynamic mechanical boundary conditions of the digital twin model.

[0019] In this embodiment, the deformation characteristics of the high stability deformation area include the deformation gradient tensor and the deformation velocity. Specifically, the deformation gradient tensor F1 is obtained by calculating the deformation displacement vectors of all pixels in the area. For the set of pixels in the high stability deformation area, the deformation displacement vector u(x) represents the displacement of the pixel in space, and the deformation gradient tensor F1 is defined as ,in is the Jacobian matrix of the displacement gradient, u(x,t) is the deformation displacement vector of pixel point x at time t, and for each pixel point x, the spatial derivative of the deformation displacement vector u(x,t) at time t is calculated to obtain I is the identity matrix. The deformation velocity v is calculated by the time difference of the deformation displacement vector at consecutive moments, that is, v = Δu / Δt, where Δu is the displacement difference between adjacent moments and Δt is the time interval.

[0020] This implementation utilizes a dual screening mechanism based on deformation stability and boundary curvature to ensure that the selected regions possess both stable deformation characteristics and geometric properties suitable for mechanical analysis. Compared to traditional methods that rely solely on a single metric, this method significantly improves the reliability and applicability of screening results, avoids the numerical instability that can result from irregular boundary regions, and provides high-quality observation areas for the subsequent construction of dynamic mechanical boundary conditions.

[0021] In an embodiment of the present invention, a method for constructing dynamic mechanical boundary conditions of a digital twin model includes: Extract the deformation features of each high-stability deformation area in the dynamic deformation field at each moment. The deformation features include deformation gradient tensor and deformation velocity. Calculate the equivalent stiffness coefficient of each high-stability deformation area based on the deformation gradient tensor and deformation velocity; According to the joint torque data, the torque distribution vector corresponding to each high-stability deformation area is obtained; The equivalent stiffness coefficient and the moment distribution vector are matrix mapped to obtain the dynamic mechanical boundary conditions corresponding to the high stability deformation area.

[0022] In this embodiment, the key deformation features of each high-stability deformation area are first extracted from the dynamic deformation field; the deformation features mainly include the deformation gradient tensor (describing the direction and amplitude of local deformation) and the deformation velocity (describing the time rate of deformation); the deformation gradient tensor F1 is obtained by calculating the spatial derivative of the displacement field; the deformation velocity v is obtained by calculating the time difference of the displacement field at consecutive moments; based on the theory of continuum mechanics, a hyperelastic model such as Neo-Hookean or Mooney-Rivlin is used to combine the deformation gradient tensor and deformation velocity to establish a stress-strain relationship; each high-stability deformation is calculated by least squares fitting or neural network regression. The equivalent stiffness coefficient K of the region is determined, which reflects the mechanical properties of the local tissue. At the same time, based on the joint torque data M measured by the sensor, combined with the robot kinematic model and the human biomechanical model, the decomposition vector of the torque in each spatial direction is calculated. A mapping relationship between the high-stability deformation region and the joint torque is established to obtain the torque distribution vector T_m corresponding to each region. Finally, the equivalent stiffness coefficient K and the torque distribution vector T_m are mapped through matrix operations to generate dynamic mechanical boundary conditions in the form of B=K·T_m. These boundary conditions include various forms such as displacement constraints, force constraints, or mixed constraints, which fully describe the relationship between tissue deformation and external forces.

[0023] In this embodiment, based on the Neo-Hookean hyperelastic model, the relationship between the stress tensor σ and the deformation gradient tensor F1 is σ=λtr(E)I+2μE, where E is the Green-Lagrange strain tensor, , λ and μ are Lamé constants, tr represents the matrix trace, and T is the transpose of the matrix or vector. The damping effect is introduced by the deformation velocity v, and the damping stress σ_d = η·v, where η is the damping coefficient. The stress tensor σ and the damping stress σ_d are superimposed to obtain the total stress σ_total = σ + σ_d. The equivalent stiffness coefficient K is obtained by fitting the linear relationship between the total stress and the deformation gradient tensor using the least squares method.

[0024] It needs to be explained that in actual applications, the deformation of the patient's limb soft tissue during rehabilitation training is usually small (i.e., F1 is close to the unit matrix I, and E is small).

[0025] Therefore, the complex nonlinear stress-strain relationship is approximated as a linear relationship: σtotal≈K·some deformation; the “some deformation” here needs to be associated with the deformation gradient tensor F1 and the deformation velocity v.

[0026] Specifically, an equivalent deformation D1 is defined such that σtotal ≈ K·D1 (equivalent to a certain deformation mentioned above); D1 here comprehensively considers the deformation gradient tensor F1 and the deformation velocity v; Specifically, based on F1, the equivalent strain ε_eq= , represents the magnitude of deformation; extract the equivalent velocity from v , which represents the amplitude of the velocity; the equivalent deformation D1=α1·ε_eq+β1·v_eq is obtained comprehensively, where α1 and β1 are weighting coefficients used to balance the contributions of elastic effect and damping effect. The typical values ​​are α1=0.7 and β1=0.3.

[0027] By using the least squares method to solve K, the error between σtotal and K·D1 is minimized; thus, the equivalent stiffness coefficient K is obtained.

[0028] The obtained K is a scalar that represents the equivalent stiffness characteristics of the high-stability deformation region. It integrates the elastic stiffness (determined by F and the Neo-Hookean model) and the damping stiffness (determined by v and the damping coefficient η), reflecting the comprehensive mechanical behavior of this region in dynamic deformation.

[0029] In order to further clarify the physical meaning of B=K·T_m and the specific content of the dynamic mechanical boundary conditions, in this embodiment, B represents the dynamic mechanical boundary condition vector, and its physical meaning is the mechanical response of the high-stability deformation area under the action of external forces, including displacement constraints, force constraints or mixed constraints. K is the equivalent stiffness coefficient, T_m is the torque distribution vector, and B=K·T_m means mapping the torque distribution to the mechanical constraints at the boundary through the stiffness characteristics. Specifically, the displacement constraint is expressed as the displacement value of the boundary node, the force constraint is expressed as the external force value of the boundary node, and the mixed constraint is expressed as a combination of displacement and external force. The decomposition of the torque distribution vector T_m in each spatial direction is realized through the robot kinematic model, that is, T_m=[T_mx,T_my,T_mz] T , where T_mx, T_my, and T_mz are the components of the torque in the x, y, and z directions respectively. The dynamic boundary condition B in different spatial directions is expressed as B=[B_x,B_y,B_z] T , where each component corresponds to the mechanical constraint in the corresponding direction.

[0030] This example establishes precise dynamic mechanical boundary conditions by extracting the mechanical characteristics of highly stable deformation regions and combining them with measured joint torque data. Compared to traditional methods that rely on simplified boundary assumptions or empirical parameters, this method, based on actual observational data, constructs boundary conditions that more closely resemble the true mechanical behavior of human tissue, accurately reflecting individual differences and dynamic changes, and providing reliable driving conditions for high-precision updates of digital twin models.

[0031] In an embodiment of the present invention, based on dynamic mechanical boundary conditions and a preset low-order finite element grid, the motion state of the digital twin model is updated in real time, including: Constructing a preset low-order finite element mesh, which includes adaptive mesh elements based on high-stability deformation regions; Mapping dynamic mechanical boundary conditions to boundary nodes of low-order finite element meshes to generate boundary constraint equations; By iteratively solving the boundary constraint equations, the stress and displacement distributions of the low-order finite element mesh are obtained, and the motion state of the digital twin model is updated.

[0032] In this embodiment, a preset low-order finite element grid is first constructed, using low-order unit types such as tetrahedron units (for volume deformation analysis) or shell units (for surface analysis); the grid density is higher near the high-stability deformation area and relatively sparse away from these areas, forming an adaptive grid structure; the use of low-order units (such as linear tetrahedrons, linear hexahedrons, etc.) reduces the computational complexity while maintaining sufficient simulation accuracy, which is particularly suitable for real-time computing requirements; the dynamic mechanical boundary conditions are mapped to the corresponding boundary nodes of the finite element grid, including converting displacement constraints, force constraints or mixed constraints into boundary conditions on the nodes; for displacement boundary conditions, the displacement value of the node is directly specified; for force boundary conditions, the force is converted into a node force; for mixed boundary conditions, a penalty method or Lagrange multiplication is used. The boundary constraint equation is established in the form of K1·U=F2, where K1 is the stiffness matrix, U is the node displacement vector, and F2 is the node force vector. Taking into account the nonlinear characteristics of human soft tissue, iterative methods such as Newton-Raphson are used to solve the boundary constraint equation. In each iteration, the tangent stiffness matrix is ​​updated and the residual is calculated until convergence. The complete displacement field and stress field on the finite element mesh are obtained by solving the problem. The displacement field describes the geometric characteristics of tissue deformation, and the stress field reflects the mechanical state inside the tissue. The motion state of the digital twin model is updated according to the calculation results, including complete information such as limb posture, joint angle, and soft tissue deformation. The entire calculation process is optimized for real-time performance, and a complete update can usually be completed within 50-100ms.

[0033] In this embodiment, the node displacement vector U represents the displacement of each node in the finite element mesh, and its initial value is directly specified by the displacement constraint in the dynamic mechanics boundary condition B, or is gradually updated through iterative solution. The node force vector F2 represents the external force acting on each node, and its value is directly specified by the force constraint in the dynamic mechanics boundary condition B, or is obtained by interpolation calculation of the torque distribution vector T_m on the mesh nodes. Specifically, the torque distribution vector T_m is decomposed into force components in the spatial direction through the robot kinematic model, and is distributed to the boundary nodes through interpolation method to form the node force vector F2. The stiffness matrix K1 is expanded from the equivalent stiffness coefficient K. Specifically, the stiffness characteristics of each mesh unit are assembled through the finite element method to obtain the global stiffness matrix K1.

[0034] In this embodiment, the Newton-Raphson iterative method is used to solve the boundary constraint equation K1·U=F2. The specific steps are: (1) Initialize the node displacement vector U_0, which is usually set to a zero vector or a displacement value based on the previous moment; (2) Calculate the residual vector R=F2-K1·U_0, where the residual vector R represents the error between the current solution and the actual constraint; (3) Calculate the tangent stiffness matrix K_t, which is a linearized approximation of the stiffness matrix K1 at the current displacement U_0, usually calculated by numerical methods; (4) Solve the linear equation K_t·ΔU=R to obtain the displacement increment ΔU; (5) Update the displacement vector U_1=U_0+ΔU; (6) Repeat steps (2)-(5) until the modulus length |R| of the residual vector R is less than the preset convergence threshold (typical value is 10 -6 ), the iteration is considered to have converged.

[0035] In this embodiment, the displacement distribution is directly used to update the posture of the limbs and the deformation of the soft tissue in the digital twin model, that is, the coordinates of each node in the model are adjusted through the displacement vector U to reflect the geometric changes of the limbs; the stress distribution is used to update the angle of the joint, that is, the torque change at the joint is calculated through the stress distribution, and the adjustment amount of the joint angle is inferred in combination with the robot kinematic model, thereby updating the joint state of the digital twin model.

[0036] It should be noted that the equivalent stiffness coefficient K is a scalar, which represents the local stiffness characteristics of the high-stability deformation area; the stiffness matrix K1 is a matrix, which represents the global stiffness characteristics of the finite element grid and is assembled from the local stiffness coefficients of each grid unit.

[0037] This embodiment achieves real-time updates of the digital twin model through a low-order finite element mesh and an efficient solution algorithm. Compared to the computational complexity and slow update limitations of traditional high-order finite element methods, this method significantly improves computational efficiency while maintaining accuracy, meeting the stringent real-time requirements of rehabilitation robot control. The adaptive mesh strategy and boundary condition mapping mechanism further enhance the model's ability to accurately describe key areas, ensuring that the digital twin model accurately reflects the actual motion state of the patient's limbs.

[0038] In an embodiment of the present invention, generating a dynamic compensation instruction for the motion trajectory of the rehabilitation robot according to the motion state includes: According to the motion state of the digital twin model, the deviation vector between the desired and actual poses of the rehabilitation robot end effector is calculated; Generate dynamic compensation instructions for the rehabilitation robot's motion trajectory based on the deviation vector and preset impedance control parameters; The driving joint angles of the rehabilitation robot are adjusted through dynamic compensation instructions to achieve real-time compensation of the motion trajectory.

[0039] In this embodiment, the real motion state of the patient's limbs, including joint position, soft tissue deformation, etc., is first obtained through the digital twin model; based on the human biomechanical model and the requirements of the rehabilitation task, the ideal expected posture of the end effector of the rehabilitation robot (such as a manipulator or a strap) is calculated at this time; at the same time, the actual posture of the end effector is obtained through the forward kinematics and sensor feedback of the robot; the deviation vector Δx between the expected posture and the actual posture is calculated, including the position deviation and the posture deviation; the preset impedance control parameters are set, including the impedance stiffness matrix K2 (usually a diagonal matrix, corresponding to the stiffness coefficients in different directions), the damping matrix D and the inertia matrix M; based on the impedance control model F=K2×Δx+D×Δx1+M×Δx2, the compensation force F to be applied is calculated, where Δx1 and Δx2 are the first-order and second-order derivatives of the deviation, respectively; the compensation force F is reversely mapped to the joint space through the Jacobian matrix J of the robot, and the required joint torque τ=J T F; Calculates joint angle adjustment Δθ based on joint torque and the robot's dynamics model; Generates formatted dynamic compensation instructions containing control parameters such as the target angle, angular velocity, and acceleration of each driven joint; Executes the dynamic compensation instructions through the robot's underlying controller, adjusting the robot's motion trajectory in real time; The compensation process is typically completed within 10-30ms, ensuring real-time tracking and adaptation to the patient's movements.

[0040] In this embodiment, the adjustment amount Δθ of the joint angle is calculated by the robot dynamics model, specifically: τ=M(θ)·Δθ2+C(θ,θ1)·Δθ1+G(θ), where τ is the joint torque vector, representing the torque acting on each joint of the robot. Its value is the joint torque calculated by the impedance control model. Δθ2 is the joint angle acceleration vector, representing the acceleration of each joint angle. Its value is the unknown quantity to be solved and is obtained by solving the equation; M(θ) is the inertia matrix, representing the inertial characteristics of the robot joint, calculated by the robot dynamics model, and depends on the current joint angle θ; Δθ1 is the joint angle velocity vector, representing the velocity of each joint angle; C(θ,θ1) is the Coriolis force and centrifugal force matrix, representing the force generated by velocity coupling in the robot motion, which depends on the current joint angle θ and the joint angular velocity θ1; G(θ) is the gravity term, representing the torque of the robot joint affected by gravity.

[0041] By solving this equation (such as the Newton-Euler method or the Lagrange method), the joint angle adjustment Δθ is obtained; In the actual solution, the joint torque τ is a known quantity calculated by the impedance control model; M(θ), C(θ,θ1), and G(θ) are known matrices or vectors calculated by the robot dynamics model, which depend on the current joint angle θ and the joint angle velocity vector Δθ1.

[0042] Δθ1 is the joint velocity at the current moment, which is usually provided by the system state (or sensor measurement). Δθ1 at the initial moment is known, and Δθ1 at subsequent moments is updated through numerical integration.

[0043] Therefore, the only real unknown in the equation is Δθ2 (joint angle acceleration vector), which can be solved directly by linear algebra methods: ; The relationship from Δθ2 to Δθ can be established through numerical integration: Update of joint velocity: Δθ1(t+Δt)=Δθ1(t)+Δθ2(t)·Δt; Update of joint angles: Δθ(t+Δt)=Δθ(t)+Δθ1(t)·Δt; Where Δt is the time step, Δθ(t) is the joint angle adjustment at time t, Δθ1(t) is the joint angle velocity vector at time t, Δθ2(t) is the joint angle acceleration vector at time t, Δθ(t+Δt) is the joint angle adjustment at time t+Δt, and Δθ1(t+Δt) is the joint angle velocity vector at time t+Δt; The target angle of each driven joint is θ_target = θ + Δθ, where θ is the current joint angle; the target angular velocity θ1_target and target acceleration θ2_target are calculated by the first-order and second-order time derivatives of Δθ, respectively.

[0044] It should be noted that the impedance stiffness matrix K2 is a matrix in the impedance control parameters, which represents the stiffness characteristics of the robot end effector and is used to control the magnitude of the compensation force.

[0045] This embodiment achieves dynamic compensation of the rehabilitation robot's motion trajectory through an impedance control strategy based on a digital twin model. Compared to the rigidity of traditional fixed trajectory control and the limitations of a predefined trajectory library, this method dynamically adjusts the robot's behavior based on the patient's real-time status, providing personalized rehabilitation assistance. The impedance control mechanism enables the robot to exhibit appropriate "compliance," enabling it to guide patients through standard movements while adapting to their motor skills and comfort needs, significantly improving the safety and effectiveness of rehabilitation training.

[0046] In an embodiment of the present invention, a method for constructing a dynamic deformation field includes: Perform feature point matching on the pixels in the surface deformation image to obtain the deformation displacement vector of each pixel; According to the deformation displacement vector, the deformation gradient tensor between any two pixels in the surface deformation image is calculated; According to the deformation gradient tensor and the local curvature change of the pixel points in the surface deformation image, the dynamic deformation field at the corresponding moment is constructed.

[0047] In this embodiment, the surface deformation image at continuous moments is first processed by a computer vision algorithm to identify and track feature points in the image; a feature descriptor method such as SIFT, SURF or ORB is used to extract stable feature points in the surface deformation image; a feature matching algorithm such as KLT (Kanade-Lucas-Tomasi) optical flow method or a descriptor-based matching method is applied to establish a pixel-level correspondence between consecutive frames; the displacement difference between each matching point pair is calculated to generate a dense deformation displacement vector field; sparse areas in the deformation displacement vector field are filled using methods such as bilinear interpolation or radial basis function (RBF) interpolation to ensure the continuity of the deformation field; and the obtained The deformation displacement vector is used to calculate the deformation gradient tensor F1 between any two adjacent pixels in the surface deformation image; at the same time, the local curvature change of each pixel in the surface deformation image is calculated, including features such as principal curvature, mean curvature and Gaussian curvature; the local curvature change reflects the local change characteristics of the surface shape and provides additional geometric constraints for the deformation field; the deformation gradient tensor is fused with the local curvature change information to construct a complete dynamic deformation field; the dynamic deformation field is represented by a tensor field, and each spatial point contains complete deformation information, including multiple dimensions such as displacement, strain, and velocity; the update frequency of the deformation field is synchronized with the image acquisition frequency, usually 25-30Hz, to ensure that transient changes in limb movement are captured.

[0048] This embodiment constructs a complete dynamic deformation field describing the deformation of the patient's limb soft tissue through high-precision feature point matching and deformation gradient calculation, combined with local curvature analysis. Compared with the limitations of traditional methods that only focus on surface contours or joint positions, this method captures the detailed characteristics of soft tissue deformation, providing a rich information foundation for subsequent deformation stability analysis and high-stability area screening. The dynamic deformation field acts as a bridge connecting image data and mechanical models, enabling the system to infer the internal mechanical state of the tissue from visual observations, providing key support for the construction of digital twin models.

[0049] In an embodiment of the present invention, a method for dividing a low-order finite element grid into adaptive grid units includes: Generate initial mesh elements based on the boundary shape of the high-stability deformation area; Calculate the mean deformation gradient of each initial grid cell and record it as the deformation complexity of the grid cell; The initial grid cells whose deformation complexity is greater than the preset complexity threshold are subdivided twice to obtain adaptive grid cells.

[0050] In this embodiment, first, according to the boundary shape of the identified high-stability deformation area, an initial grid unit is generated using algorithms such as Delaunay triangulation or tetrahedron partitioning; the initial grid maintains a relatively uniform unit size and shape quality, usually using tetrahedron or hexahedron units; for each initial grid unit, the deformation gradient data of the area covered by the unit is extracted from the dynamic deformation field; the mean of the deformation gradient tensor in the unit is calculated, and the mean reflects the complexity of the deformation of the area; the calculation of the deformation gradient mean takes into account all components of the tensor, usually using the Frobenius norm or eigenvalue analysis method; the calculated deformation gradient mean is defined as the deformation complexity of the grid unit, and the higher the deformation complexity, the more complex the deformation of the area; a preset complexity threshold (typically ) is set. 0.15-0.25), which is dynamically adjusted according to application requirements and computing resources; an adaptive subdivision strategy is applied to initial grid cells whose deformation complexity is greater than a preset threshold; subdivision methods include uniform subdivision (uniformly dividing a cell into multiple subcells) and non-uniform subdivision (directional subdivision according to the gradient direction); secondary subdivision is generally used, that is, the original cell is divided into 8 (tetrahedron) or 27 (hexahedron) subcells; for particularly complex areas, multi-level subdivision can be performed, but usually no more than three levels, to balance accuracy and computational efficiency; the subdivided grid cells are dense in high deformation complexity areas and sparse in low deformation complexity areas, forming an adaptive grid structure; the number of adaptive grid cells finally generated is usually controlled between 5,000 and 20,000 to ensure real-time computing performance.

[0051] In this embodiment, the boundary shape of the high-stability deformation region is determined by the following steps: (1) extracting the boundary pixel points of the high-stability deformation region in the surface deformation image, where the boundary pixel points are the dividing points between the pixels within the region and the pixels outside the region; (2) performing curve fitting on the boundary pixel points to generate a continuous boundary contour curve, usually using a B-spline curve or a polynomial curve; (3) the boundary shape is the geometric description of the boundary contour curve, including features such as curvature, length, and direction.

[0052] In this embodiment, the deformation gradient mean of the initial grid cells is calculated by the following steps: (1) Determine the set of pixels covered by the grid cells. Specifically, the initial grid cell division is based on the geometric shape of the surface deformation image at a certain time t and is generated through a Delaunay triangulation or tetrahedronization algorithm. The set of pixels covered by each grid cell refers to all the pixels within the projection area of ​​the cell in the surface deformation image.

[0053] (2) For each pixel in the set, obtain its deformation displacement vector u(x); it should be noted that the deformation displacement vector u(x) is not directly calculated when the grid unit is divided, but is pre-calculated during the construction of the dynamic deformation field. Specifically, the dynamic deformation field is a deformation displacement vector field u1(x) calculated by using feature point matching and pixel tracking algorithms on the surface deformation images at consecutive moments. In the dynamic deformation field, the deformation displacement vector u(x, t) of each pixel x at time t represents the displacement difference of the pixel from the reference state (usually the initial time t_0) to the current time t, that is, u(x, t) = x(t) - x(t_0), where x(t) and x(t_0) are the spatial coordinates of the pixel at time t and t_0, respectively. Through feature matching algorithms (such as KLT optical flow method or deep learning-based pixel matching algorithms), pixel-level correspondences are established between consecutive frames, and the displacement difference between each matching point pair is calculated, thereby generating a dense deformation displacement vector field u1(x). Therefore, when calculating the mean deformation gradient of the initial grid cell, the deformation displacement vector u(x, t) at the time t is directly extracted from the dynamic deformation field.

[0054] (3) Calculate the deformation gradient tensor ,in is the Jacobian matrix of the displacement gradient; specifically, for each pixel point x in the set of pixel points covered by the grid unit, the spatial derivative of the deformation displacement vector u(x,t) at time t is calculated to obtain .

[0055] (4) Take the average value of the deformation gradient tensor F1 of all pixels in the grid unit to obtain the mean deformation gradient, which is usually expressed using the Frobenius norm.

[0056] This embodiment achieves efficient allocation of computing resources through an adaptive mesh subdivision strategy based on deformation complexity. Compared to the redundant calculations of uniform meshes or the subjectivity of empirical subdivision, this method automatically adjusts the mesh density based on the deformation characteristics, providing high-precision expression in key areas while maintaining computational simplicity in less important areas. This adaptive mesh strategy significantly improves the balance between accuracy and efficiency of numerical simulations, ensuring maximum simulation results with limited computing resources and meeting the requirements of real-time updates of digital twin models.

[0057] In an embodiment of the present invention, a deformation vector tracking method includes: A pixel matching algorithm based on deep learning is used to extract features from pixels in the surface deformation image and obtain the feature descriptor of each pixel. According to the feature descriptor, the matching probability of the same pixel in the surface deformation image at adjacent moments is calculated; According to the matching probability, a deformation vector is generated for each pixel.

[0058] In this embodiment, a deep learning method is first used to extract features and match pixels of the surface deformation image; a pre-trained convolutional neural network (CNN) such as ResNet, VGG or a specially designed U-Net is used as a feature extractor; a feature descriptor is extracted for each pixel in the image, and the descriptor dimension is usually 128-256, which includes the texture, gradient and semantic features of the area around the point; the feature descriptor has rotation invariance, scale invariance and partial illumination invariance, ensuring that stable features can be maintained during the deformation process; for two frames of surface deformation images at consecutive moments, the similarity between the feature descriptors of each pixel is calculated, and commonly used similarity metrics include cosine similarity or L2 distance ; Establish a feature descriptor matching matrix between the two frames of images, where each element in the matrix represents the matching probability of the corresponding pixel pair; use the matching probability as the weight, combined with the local spatial consistency constraint and the temporal continuity constraint, to optimize the pixel correspondence relationship; optimization methods include graph cut algorithm, Markov random field (MRF) or variational inference; for pixel pairs with a matching probability higher than a threshold (usually 0.7-0.8), calculate the position difference between them to form a pixel-level deformation vector; for low matching probability or many-to-one matching, use the information of the surrounding high-confidence matching points for interpolation or inference; finally generate a complete deformation vector field that covers all pixels in the image, and each vector contains displacement direction and amplitude information.

[0059] In this embodiment, the matching probability P is obtained by normalizing the feature descriptor similarity S. Specifically, for the feature descriptors f_i and f_j of pixel points i and j, the similarity S(i, j) = cos(f_i, f_j), where cos represents cosine similarity; ... , the denominator represents the sum of all possible matching pixel points j, and the probability distribution is obtained after normalization.

[0060] This embodiment achieves high-precision pixel matching and deformation vector tracking through deep learning technology. Compared to the limitations of traditional feature point methods or optical flow methods, deep learning-based methods can extract richer semantic features and contextual information, improving tracking accuracy under conditions of complex textures, large deformations, and varying illumination. High-quality deformation vectors are the foundation of deformation stability analysis, directly affecting the reliability of subsequent analysis results and providing a stable data input source for the system.

[0061] In an embodiment of the present invention, a method for acquiring a surface deformation image includes: The multi-view depth camera installed on the rehabilitation robot collects 3D point cloud data of the patient's limbs during movement. Surface reconstruction is performed on the three-dimensional point cloud data to generate a surface deformation image of the patient's limb.

[0062] In this embodiment, multiple depth cameras are first installed at appropriate positions of the rehabilitation robot to form a multi-view observation system; 3-5 depth cameras are usually configured to cover various key angles of the patient's limbs during rehabilitation training; the depth cameras using infrared structured light (such as Kinect), time of flight (ToF) or binocular stereo vision technology have a resolution of no less than 640×480 and a depth accuracy of millimeter level; the camera acquisition frequency is set to 30-60Hz to ensure that transient changes in rapid motion are captured; multiple depth cameras are strictly spatially calibrated and time synchronized to establish a unified world coordinate system; the original depth image and RGB image of the patient's limb during rehabilitation exercise are collected; 3D point cloud data are extracted from the depth image, where each point contains (x, y, z) spatial coordinates and corresponding color information; align and fuse point cloud data from different perspectives using algorithms such as iterative closest point (ICP) or feature matching; remove noise points, outliers, and overlapping redundant points in the point cloud to improve data quality; based on the streamlined point cloud data, use Poisson surface reconstruction, moving least squares (MLS) or implicit function methods to reconstruct the surface; the reconstructed surface model is usually represented by a triangular mesh or NURBS surface; the surface model is projected onto a 2D plane to generate a surface deformation image containing depth information; the surface deformation image contains rich information such as the spatial coordinates, normal vector, curvature, etc. of each pixel point; to facilitate subsequent analysis, the surface deformation image is organized into a time series data structure, retaining complete timestamps and camera parameter information.

[0063] This embodiment utilizes a multi-view depth camera and high-precision surface reconstruction technology to fully capture the surface deformation of a patient's limbs. Compared to the blind spots and information loss of single-view observations, multi-view 3D reconstruction provides complete spatial deformation information, overcoming the limitations of traditional methods in situations such as limb self-occlusion and complex postures. High-quality surface deformation images provide a reliable data foundation for subsequent dynamic deformation field construction and deformation analysis, and serve as a critical input source for the entire system.

[0064] In the embodiment of the present invention, the detailed implementation steps of tracking the deformation vectors of the pixels in the surface deformation image at all times include: Local feature enhancement preprocessing is applied to each frame of the surface deformation image to improve texture contrast and edge clarity; consecutive multiple frames of images are organized into a time series buffer to ensure temporal continuity analysis; multi-scale feature descriptors are extracted for each pixel, including local oriented gradient histogram (HOG), local binary pattern (LBP) and depth features; a recursive optical flow algorithm is applied to track the displacement of each pixel between consecutive frames, and forward-backward consistency check is used to verify tracking reliability; temporal smoothing filtering is applied to the tracking results to eliminate high-frequency noise and instantaneous jumps, generating a continuous and smooth deformation vector field; for pixels that fail to be tracked, spatial coherence constraints are used for interpolation estimation to ensure the integrity of the deformation vector field.

[0065] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art will be able to modify the technical solutions described in the foregoing embodiments or to substitute equivalents for some of the technical features. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention shall be included within the scope of protection of the present invention.

[0066] It should be noted that the formulas in this manual are all dimensionless and calculated numerically. The formulas are obtained by collecting a large amount of data and performing software simulation to obtain the most recent real situation. The preset parameters and thresholds in the formulas are set by technicians in this field according to actual conditions.

[0067] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to the embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the claims and their equivalents.

Claims

1. The digital twin monitoring and dynamic compensation system for the motion state of rehabilitation robots is characterized by: include: An acquisition module is used to obtain the surface deformation image and joint torque data of the patient's limb at each moment during the movement of the rehabilitation robot; a deformation field construction module, configured to construct a dynamic deformation field of the patient's limb soft tissue based on the deformation gradient distribution and local curvature change between pixels in the surface deformation image; a stability analysis module, configured to obtain the deformation stability of each local area in the dynamic deformation field at each moment based on the consistency of deformation directions and deformation amplitude fluctuations of the same local area in the dynamic deformation field at adjacent moments within an analysis period at each moment; Screening out high-stability deformation regions based on the deformation stability; A boundary condition module is used to construct the dynamic mechanical boundary conditions of the digital twin model based on the deformation characteristics of each high-stability deformation area in the dynamic deformation field at each moment and the corresponding joint torque data; A model updating and compensation module is used to update the motion state of the digital twin model in real time based on the dynamic mechanical boundary conditions and a preset low-order finite element grid, and to generate dynamic compensation instructions for the motion trajectory of the rehabilitation robot according to the motion state; The modules are connected via wired and / or wireless means to achieve data transmission between modules.

2. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 1 is characterized in that: The obtaining of the deformation stability of each local area in the dynamic deformation field at each moment includes: Tracking deformation vectors of pixels in the surface deformation image at all moments to obtain the deformation vector of each pixel in the surface deformation image at each moment; Counting the number of times that the direction angle of the deformation vectors of the same local area in the dynamic deformation field at adjacent moments within the analysis period at each moment is less than a preset direction threshold, and recording this as the deformation direction consistency of each local area in the dynamic deformation field at each moment; Calculating the standard deviation of the amplitude of the deformation vector of the same local area in the dynamic deformation field at all moments within the analysis period at each moment as the deformation amplitude fluctuation of each local area in the dynamic deformation field at each moment; The deformation direction consistency and the inverse of the deformation amplitude fluctuation are weightedly fused to obtain the deformation stability of each local area in the dynamic deformation field at each moment.

3. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 1 is characterized in that: The method for screening the high-stability deformation region comprises: sorting the deformation stability of each local area in the dynamic deformation field at each moment, and selecting a local area whose deformation stability is greater than a preset stability threshold as a candidate area; The mean value of the boundary curvature of each candidate region in the surface deformation image is calculated, and candidate regions whose mean value of the boundary curvature is greater than a preset curvature threshold are eliminated to obtain highly stable deformation regions.

4. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 1 is characterized in that: The method for constructing the dynamic mechanical boundary conditions of the digital twin model includes: Extracting deformation features of each high-stability deformation region in the dynamic deformation field at each moment, wherein the deformation features include a deformation gradient tensor and a deformation velocity; Calculating an equivalent stiffness coefficient of each high-stability deformation region according to the deformation gradient tensor and the deformation velocity; Obtaining a torque distribution vector corresponding to each high-stability deformation region according to the joint torque data; The equivalent stiffness coefficient and the moment distribution vector are subjected to matrix mapping to obtain the dynamic mechanical boundary conditions corresponding to the high stability deformation region.

5. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 1, characterized in that: The real-time updating of the motion state of the digital twin model based on the dynamic mechanical boundary conditions and the preset low-order finite element grid includes: Constructing a preset low-order finite element grid, wherein the low-order finite element grid includes adaptive grid units based on high-stability deformation region division; Mapping the dynamic mechanical boundary conditions to the boundary nodes of the low-order finite element mesh to generate boundary constraint equations; By iteratively solving the boundary constraint equations, the stress distribution and displacement distribution of the low-order finite element grid are obtained, and the motion state of the digital twin model is updated.

6. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 1, characterized in that: The step of generating a dynamic compensation instruction for the motion trajectory of the rehabilitation robot according to the motion state includes: Calculating a deviation vector between a desired pose and an actual pose of an end effector of the rehabilitation robot according to a motion state of the digital twin model; Generating a dynamic compensation instruction for the motion trajectory of the rehabilitation robot according to the deviation vector and a preset impedance control parameter; The driving joint angles of the rehabilitation robot are adjusted by the dynamic compensation instructions to achieve real-time compensation of the motion trajectory.

7. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 2, characterized in that: The method for constructing the dynamic deformation field comprises: Performing feature point matching on pixel points in the surface deformation image to obtain a deformation displacement vector for each pixel point; Calculating a deformation gradient tensor between any two pixel points in the surface deformation image according to the deformation displacement vector; A dynamic deformation field at a corresponding moment is constructed according to the deformation gradient tensor and the local curvature change of the pixel points in the surface deformation image.

8. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 5, characterized in that: The method for dividing the adaptive grid units of the low-order finite element grid includes: generating an initial grid cell according to the boundary shape of the high-stability deformation region; Calculate the mean deformation gradient of each initial grid cell and record it as the deformation complexity of the grid cell; The initial grid cells whose deformation complexity is greater than the preset complexity threshold are subdivided twice to obtain adaptive grid cells.

9. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 2, characterized in that: The deformation vector tracking method includes: Using a deep learning-based pixel matching algorithm, feature extraction is performed on the pixels in the surface deformation image to obtain a feature descriptor for each pixel; Calculating, based on the feature descriptor, the matching probability of the same pixel point in the surface deformation image at adjacent moments; A deformation vector for each pixel is generated according to the matching probability.

10. The rehabilitation robot motion state digital twin monitoring and dynamic compensation system according to claim 1, characterized in that: The method for acquiring the surface deformation image comprises: The multi-view depth camera installed on the rehabilitation robot collects 3D point cloud data of the patient's limbs during movement. Surface reconstruction is performed on the three-dimensional point cloud data to generate a surface deformation image of the patient's limb.

Images