A gait analysis method and system that fuses depth camera and inertial sensors

Abstract

The application discloses a gait analysis method and system fusing a depth camera and an inertial sensor, and the method comprises the following steps: collecting human gait data; pre-processing the original human gait data; performing gait-related calculation by using the processed data; performing sensor data fusion by using a fusion algorithm; and performing gait analysis. The application fuses the noisy accelerometer and the gyro with drift in the inertial measurement unit by using a complementary filter, so that more accurate joint angle results are obtained, the IMU is installed on a limb in an arbitrary direction, does not depend on a uniform magnetic field, avoids noise and errors caused by magnetic field interference, is suitable for single-limb rotation, and combines the joint angle results calculated by using the depth camera data, fuses the depth camera and the IMU data, and when a motion with a larger measurement error is detected, the depth camera and the IMU are used to mutually correct the joint angle, so that the calculation precision is improved without reducing the operation efficiency.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of inertial sensors and vision sensors, and particularly relates to a gait analysis method and system fusing a depth camera and an inertial sensor. BACKGROUND

[0002] Traditional gait analysis requires a special laboratory and expensive equipment, thus limiting its application range. With the progress of science and technology, more and more reliable and economical gait analysis sensors appear, which can be used outside the laboratory environment and in daily activities. Vision sensors and inertial measurement units (IMU) are two low-cost sensors, which are increasingly applied to gait analysis systems for home rehabilitation exercises in recent years.

[0003] Depth cameras are widely used in many fields due to their low cost, multiple uses and strong functions. However, when used as a skeletal tracking system, the depth camera often loses skeletal data due to occlusion and other problems, and therefore needs to be corrected and technically breakthrough. IMU is a lightweight, energy-saving and inexpensive inertial measurement unit, which is widely used in the field of gait analysis. The “Wearable inertial device human gait information acquisition and calculation method” disclosed in Chinese patent application CN 108245164 A and the “Human gait acquisition and analysis system and method based on inertial sensor information fusion” disclosed in Chinese patent application CN 108245164 A both use two IMUs and an additional accelerometer to obtain knee joint angle information, and the installation direction of the sensor is limited, and the measurement using the IMU often drifts over time, affecting the calculation accuracy. Therefore, through multi-sensor fusion technology, the advantages of various sensors are combined, which is an inevitable development trend of gait analysis research. SUMMARY

[0004] The purpose of the present application is to provide a gait analysis method and system fusing a depth camera and an inertial sensor.

[0005] The technical solution for achieving the purpose of the present application is as follows: In a first aspect, the present application provides a gait analysis method fusing a depth camera and an inertial sensor, comprising the following steps:

[0006] Step S1: collecting human gait data;

[0007] Step S2: preprocessing the original human gait data;

[0008] Step S3: using the data obtained in step S2 to perform gait-related calculations;

[0009] Step S4: fusing the sensor data in step S3;

[0010] Step S5: gait analysis using the results from steps S3, S4.

[0011] In a second aspect, the application provides a gait analysis system fusing a depth camera and an inertial sensor, comprising:

[0012] A human gait data acquisition module: for acquiring gait-related human activity data of the sensor;

[0013] A data preprocessing module: for preprocessing the original human gait data;

[0014] A gait-related calculation module: for calculating human joint angles, bone lengths, and gait parameters;

[0015] A sensor fusion module: for fusing sensor data to obtain fused joint angles.

[0016] Compared with the prior art, the application has the following beneficial effects:

[0017] (1) The application designs a depth camera data acquisition system, pre-processes the collected data based on a One Euro filter, and uses the processed joint tracking data to calculate joint angles, bone lengths, and gait parameters, thereby improving the accuracy and real-time performance of the calculation;

[0018] (2) The application designs a joint position estimation method based only on an accelerometer and a gyroscope, and uses a complementary filter to fuse the noisy accelerometer and the drifting gyroscope for joint angle calculation, which does not depend on a uniform magnetic field, avoids noise and errors caused by magnetic field interference, and is suitable for single-limb rotation; the IMU installation direction is not limited; the subject does not need to accurately perform a predefined posture and action;

[0019] (3) The application designs an extended Kalman filter scheme based on quaternions to correct the data results with large errors or missing data from the depth camera and the IMU, and calculates the joint angles by fusing the depth camera and the IMU, in order to solve the problem of missing bone data due to occlusion and the like when the depth camera is used as a bone tracking system, and the problem of drift over time when using an IMU alone to measure joint angles. BRIEF DESCRIPTION OF DRAWINGS

[0020] Figure 1 A flowchart is designed for the depth camera acquisition system.

[0021] Figure 2 A flowchart is designed for the gait cycle calculation based on the depth camera.

[0022] Figure 3 A schematic diagram of the calculation method for measuring joint angles based on the IMU.

[0023] Figure 4 Flow chart for fusing depth camera and IMU to calculate joint angle.

[0024] Figure 5 Flow chart for extended Kalman filter.

[0025] Figure 6 Schematic diagram of the invented gait analysis system fusing depth camera and inertial sensor. DETAILED DESCRIPTION

[0026] The present invention proposes a gait analysis method and system fusing depth camera and inertial sensor, using a complementary filter to fuse the noisy accelerometer and the drifting gyroscope in the inertial measurement unit, obtaining more accurate joint angle results, the IMU being installed on the limb in any direction, not relying on uniform magnetic field, avoiding noise and error caused by magnetic field interference, suitable for single-limb rotation, and combining the joint angle results calculated through depth camera data, using an extended Kalman filter algorithm based on quaternion to fuse the data of depth camera and IMU, when detecting motion with large measurement error, using depth camera and IMU to correct joint angle, improving the calculation accuracy without reducing the operation efficiency.

[0027] The gait analysis method fusing depth camera and inertial sensor of the present invention comprises the following steps:

[0028] Step S1: collecting human gait data;

[0029] Step S2: preprocessing the original human gait data;

[0030] Step S3: using the data obtained in step S2 to perform gait-related calculation;

[0031] Step S4: using the sensor data in step S3 for fusion;

[0032] Step S5: using the results obtained in steps S3 and S4 to perform gait analysis.

[0033] Further, step S1 comprises the following steps:

[0034] Step S11: designing a depth camera data acquisition system, obtaining human depth image coordinates through depth camera bone tracking technology, determining human joint positions, and thus collecting human gait-related data such as step length, stride, step width, step speed and bone length;

[0035] Step S12: using the IMU sensor containing the accelerometer, gyroscope, collecting the acceleration and angular velocity signals of the sensor on the human body; the IMU is installed on two body segments connected by the measured human joint, does not rely on the magnetometer, does not need to manually measure any distance, does not need the subject to accurately perform the predefined posture or action, and is suitable for single-limb rotation.

[0036] Further, step S2 comprises the following steps:

[0037] Step S21: using the One Euro filtering algorithm to filter the depth camera joint coordinate data obtained in step S11, so as to reduce jitter while reducing lag, so that the data is more stable while meeting rapidity;

[0038] Step S22: using the Butterworth low-pass filter to denoise the accelerometer and gyroscope signals collected by the IMU, to obtain more stable data.

[0039] Further, step S3 comprises the following steps:

[0040] Step S31: using the depth camera joint tracking data obtained in step S2 to calculate the joint angle, bone length and gait parameters;

[0041] Step S32: using the data of the IMU sensor to identify the joint axis and the position of the joint axis;

[0042] Step S33: identifying the coordinates of the joint axis obtained in step S32;

[0043] Step S34: symbol matching the joint axis coordinates in step S33;

[0044] Step S35: identifying the coordinates of the joint position obtained in step S32;

[0045] Step S36: using the IMU sensor to calculate the joint angle based on the complementary filter.

[0046] Further, step S4 comprises the following steps:

[0047] Step S41: unifying the coordinate systems of the depth camera and the IMU, converting the depth camera coordinates to obtain the depth camera data information in the global coordinate system, and when the data of the IMU and the depth camera are mainly referenced in the global coordinate system, the data fusion of the two sensors can be performed;

[0048] Step S42: designing an extended Kalman filtering scheme based on quaternion, calibrating the data results with large errors or loss from the depth camera and the IMU, and calculating the joint angle through the fusion of the depth camera and the IMU.

[0049] Further, the step S5 comprises the following steps:

[0050] The step S51: the results obtained by the steps S3 and S4 can be used for gait analysis, gait symmetry test, disease prevention, rehabilitation evaluation and other gait-related researches.

[0051] In order to make the technical solutions, objectives and advantages of the present application clearer, the present application will be described in detail below with reference to the drawings, wherein the drawings form a part of the present application and are used together with the embodiments of the present application to explain the principles of the present application, but not to limit the scope of the present application.

[0052] The present application designs a gait analysis method fusing a depth camera and an inertial sensor, fuses sensor information in multiple aspects, and provides a new idea for further improving the accuracy of gait-related calculations. The specific embodiments will be described in detail as follows:

[0053] The embodiment provides a gait analysis method fusing a depth camera and an inertial sensor, and the method comprises the following steps.

[0054] The step S1: collecting human gait data, specifically comprising the following steps.

[0055] The step S11: combining Figure 1 , a depth camera data acquisition system is designed, first, defining a Kinect device, reading depth, color and skeleton data frame variables, then performing coordinate conversion, processing depth, color and skeleton data frame variables, saving feature data, and finally drawing human skeleton. The depth data, color data and skeleton data obtained from the depth camera are read and processed to collect the position data of each skeletal joint;

[0056] The step S12: a plurality of subjects are selected to perform the gait analysis experiment, and the length of the experimental path is set to 6 meters and the width is set to 2 meters, and external interference is excluded. The depth camera sensor is placed flat at one end of the path, with a height of 0.5 meters from the ground, so as to correctly detect the frame of the human skeleton and its movement as much as possible, thereby minimizing the error of the recorded data. The subjects stand in the capture range in front of the depth camera in a natural posture, and before collecting data, the operator needs to calibrate the distance and angle parameters of the depth camera.

[0057] The step S13: an IMU sensor is arranged on the lateral side of the middle segment of the thigh, the lateral side of the middle segment of the lower leg and the lateral side of the foot of the subject using an elastic band, without being limited to a specific position or direction. The IMU acquisition system is started, and the subject stands at attention for 5 seconds, and waits for the calibration and correction of the IMU data.

[0058] Step S14: After the IMU completes data calibration and correction, the depth camera acquisition system is started. When the indicator light of the depth camera is on, the subject walks in the direction of the depth camera at a natural pace in a relaxed state and performs the following actions: walks straight until the coordinates of the human body joint cannot be identified, and stops; each subject tests several times;

[0059] Step S2, pre-process the original human gait data, specifically including the following steps:

[0060] Step S21: The One Euro filtering algorithm is used to filter the depth camera joint coordinate data obtained in step S11, to reduce jitter while reducing lag, so that the data is more stable while meeting the speed requirement.

[0061] Step S22: The Butterworth low-pass filter is used to denoise the accelerometer and gyroscope signals collected by the IMU to obtain more stable data.

[0062] Step S3, use the data obtained in step S2 to perform gait-related calculations, the specific steps are as follows:

[0063] Step S31: Use the depth camera joint tracking data obtained in step S2 to calculate joint angles, bone lengths, and gait parameters.

[0064] After collecting the spatial coordinate sequence of the human skeleton points and filtering the data, because the depth camera camera cannot be completely placed with the X axis parallel to the ground, and the zero view of the Y axis cannot be completely regarded as the landing point, the landing point of the two feet cannot be directly determined according to the returned coordinates to determine the cycle.

[0065] During walking, the human center point shows periodic oscillation in the up-down and left-right directions. The lowest point of each oscillation is selected as a special point, and the frame when the special point appears is the special frame (from the knee Y axis position image, it is the wave trough position). Gait cycle (T) is the time from one side of the heel landing to the heel landing again, that is, the time required for three consecutive special frames to appear.

[0066] T=t k+2 -t k #(1)

[0067] In formula (1), t k is the time when the kth special frame appears, and t k+2 is the time when the (k+2)th special frame appears.

[0068] Combining Figure 2 , the gait cycle can be calculated by finding the position of the special frame.

[0069] Bone length measurement is to calculate the distance between two joints connected by a bone. The three-dimensional coordinate information of the bone points is obtained by a depth camera, a joint vector is generated, and the bone length is calculated. In the 3D coordinate system of the depth camera, the 3D coordinates of the adjacent joints of the bone are retrieved by the depth camera, the quaternion is used, and the following formula (2) is used to calculate the vector of each segment:

[0070]

[0071] wherein E(X e , Y e , Z e ), F(X f , Y f , Z f ) represent the positions of two adjacent joints, respectively.

[0072] Then, the bone length is calculated:

[0073]

[0074] The step length is the distance from the position of one side of the heel landing to the position of the other side of the heel landing during walking, which is expressed by the distance formula of the spatial point coordinates:

[0075]

[0076] In formula (4), s pi represents the step length of the i-th step, X k , Y k , Z k represent the spatial coordinates at the special point number k, and X k+1 , Y k+1 , Z k+1 represent the spatial coordinates at the special point number k+1.

[0077] The stride length is the distance from the position of one side of the heel landing to the position of the other side of the heel landing during walking, which is expressed by the distance formula of the spatial point coordinates:

[0078] s di = s pi + s p(i+1) #(5)

[0079] The step width is the distance from the position of one side of the heel landing to the position of the other side of the heel landing during walking, which is expressed by the distance formula of the spatial point coordinates:

[0080] The step width is the distance from the position of one side of the heel landing to the position of the other side of the heel landing during walking, which is expressed by the distance formula of the spatial point coordinates:

[0081] AX + BY + CZ = D (7)

[0082] The distance of the left foot landing point to the plane, i.e. the step width, is represented as follows:

[0083]

[0084] where (X0, Y0, Z0) is the coordinate of the right foot landing point.

[0085] The step velocity is the instantaneous velocity of the bone point captured in each data collection. Because the collection frequency is 30 frames / s, the time interval of adjacent frames is 0.033 s, and the average velocity of the adjacent three points is used to replace the instantaneous velocity of the center point, and the instantaneous velocity of the center point is:

[0086]

[0087] where f is the depth camera collection frequency, and δ is the distance between point B and point B. c t-1 t-1 t

[0088] Acceleration: the acceleration is the change rate of adjacent frames, i.e.

[0089]

[0090] According to the collected m+1 point coordinates, the time parameter sequence can be determined, that is, the displacement sequence x1, x2, x3,... x m , the velocity sequence v1, v2, v3,... v m-1 , and the acceleration sequence a1, a2, a3,... a m-3 .

[0091] The depth camera uses its depth sensor to retrieve information in 3D coordinates. It uses its own position as a reference point to retrieve the joint points on the skeleton. The depth camera uses the relative difference between the joint points and the reference point to calculate the three-axis coordinates of the human body and the motion tracking information and the corresponding 3D joint angles. Since the depth camera can only provide a set of 3D coordinates for the human body limbs, it is necessary to calculate the joint angles according to the 3D coordinates of the back-end system.

[0092] The three-dimensional coordinate information of the bone points is obtained by the depth camera, the joint vector is generated, and the joint control angle is calculated by calculating the included angle of the joint vector. In the 3D coordinate system of the depth camera, the depth camera is used to retrieve the 3D coordinates of the adjacent joints of the measured joint, and the quaternion is used to calculate the vector of each segment using formula (11):

[0093] ​​​​

[0094] Then, the joint angle is calculated using the following equation:

[0095]

[0096] θ = cos -1 (cos θ) # (13)

[0097] Step S32: Identification of joint axis and joint axis position using data of the IMU sensor.

[0098] The present application does not assume any local sensor axis to coincide with the joint axis or limb segment. Therefore, the direction and position of the joint flexion axis are arbitrary. However, these coordinates can be identified from the measurement data of arbitrary motion by exploiting kinematic constraints. The first step of this identification is to collect identification data while the joint is moved in an arbitrary way around its degrees of freedom (i.e. without assuming any type of specific motion like walking or motion in a certain direction). Approximately every 0.1 s a data set S(i) is acquired, which is in the form of:

[0099]

[0100] where a1(t), a2(t) e R 3 is the acceleration over a certain sampling period At, g1(t), g2(t) e R 3 is the angular velocity, is the angular rate. The time between two data sets must be a multiple of the sampling period, i.e. Thus, N > > 1 data sets are collected, which will be used in the subsequent section for the identification of the local joint axis and position coordinates.

[0101] Step S33: Coordinate identification of the joint axis resulting from step S32.

[0102] The data sets S(i), i e [1, N], are used to identify the unit length direction vector j The joint flexion axis in the local coordinates of the two sensors. j1and j2are constants, depending only on the mounting direction of the sensors with respect to the joint. The angular rates g1(t), g2(t) differ only in the joint angular velocity vector and the (time-varying) rotation matrix measured on the joint. Therefore, they project into the joint plane (i.e. the geometric plane whose normal vector is the joint axis) with the same length for each time instant t, i.e.

[0103]

[0104] where ||•||2denotes the Euclidean norm. This constraint holds true regardless of the position and orientation of the sensors mounted on the segment. In particular, each data set S(i), i e [1, N] must satisfy (15). We can therefore identify j1and j2by minimizing the left-hand side of (15) for all data sets. j1and j2are expressed in spherical coordinates as:

[0105] j1= (cos(φ1)cos(θ1), cos(φ1)sin(θ1), sin(φ1)) T (16)

[0106] j2= (cos(φ2)cos(θ2), cos(φ2)sin(θ2), sin(φ2)) T (17)

[0107] and define the sum of squared errors as:

[0108]

[0109] where φ1, φ2denote the roll angles of the joint axes and θ1, θ2denote the yaw angles of the joint axes. Since the sign of (18) with respect to j1and j2is invariant, this cost function has four minima, corresponding to the four possible combinations of signs, (j1, j2), (-j1, j2), (j1, -j2), and (-j1, -j2). These true joint axis coordinates are determined by minimizing Ψ(φ1, φ2, θ1, θ2) in its arguments. This optimization can be achieved using the Gauss-Newton algorithm or any other standard optimization method.

[0110] Step S34: Sign matching of joint axis coordinates from step S33.

[0111] The gyro-based joint angular velocity approximation g1(t) x j1- g2(t) x j2will be used in the later joint angle calculation. It is therefore important to ensure that the signs of j1and j2match, i.e. that they point in the same direction. This is the case when the Z-axes of both sensors point approximately in the lateral direction (i.e. the coordinate axes point into the lateral half-space, which can not restrict the mounting orientation in any way). Conversely, if the local Y-axis of the first sensor points approximately in the medial direction and the local Z-axis of the second sensor points approximately in the lateral direction, the X-coordinate of j1and the Z-coordinate of j2should have opposite signs. If the mounting orientation of the sensors cannot be observed, the correct pairing of signs can also be determined from the inertial data itself. First, select a period from the identification data during which the angular velocity about the joint axis is negligible, i.e. g1(t) x j1~ 0, g2(t) x j2~ 0 (19)

[0112] Step S35: Coordinate recognition is performed on the joint position obtained in step S32.

[0113] The acceleration of each sensor is considered as the sum of the acceleration of the joint center and the acceleration due to the rotation of the sensor around the joint center. Obviously, in the two local coordinate systems, the acceleration of the joint center must be the same, up to a certain time-varying rotation matrix corresponding to the mutual rotation of the two local coordinate systems. Mathematically, this constraint is expressed by the following equation:

[0114]

[0115] where a is the measured acceleration of the sensor, and a is the acceleration of the joint center. is the radial and tangential acceleration due to rotation around the joint center. By subtracting a from a, the acceleration of the joint center is obtained. i (t) is the measured acceleration of the sensor. i

[0116] Each pair of coordinates, o1, o2, describes a point on the joint axis that satisfies the constraints of any given motion that the joint can perform. Therefore, using the same arbitrary motion data S(i), i∈[1, N], and defining another sum of squared errors:

[0117]

[0118] The minimization of S2 can be performed by the Gauss-Newton algorithm. As mentioned above, any other optimization method can also be used. As a result of this optimization, the joint axis is represented by o , which refers to an arbitrary point along the joint axis, we shift it to the position as close as possible to the sensor by applying the following command:

[0119]

[0120] In addition, in order to make the estimation of the joint axis more accurate, when the subject's limb 1 remains stationary and the limb 2 rotates, the IMU data set of the rotating limb is measured. Assuming that the noise of the sensor is zero, theoretically, the direction of these vectors measured by the gyroscope should be parallel to the joint axis, but in reality, the sensor is affected by noise, and the direction of each vector will deviate to a certain extent, and will not be completely parallel to the joint axis. In order to reduce or eliminate noise to some extent, these vectors can be added together. However, the measured values cannot be directly added together, because the limb rotates around the joint axis, and the direction of rotation is also positive and negative in turn, and the direction of the measured angular velocity is also positive and negative in turn, and directly adding the measured values together may cause the positive and negative angular velocity values to cancel each other out, thereby being affected by noise more. To this end, the positive and negative directions can be calculated separately, and each measured value g2(t k ​​) between the first measurement g2(t1) and the second measurement g2(t2) is calculated. If the angle is less than 90°, the measurement is kept unchanged. If the angle is greater than 90°, the opposite of the measurement is taken. The processed vectors are then added together. The resulting vector is then unitized to obtain the optimal estimate of j2. Similarly, limb 2 is kept stationary and limb 1 is rotated about the joint axis. The IMU data set on limb 1 is processed in the above manner to obtain the optimal estimate of j1.

[0121] To extend the algorithm to the case of single limb rotation, the vector from the joint center to the origin of the IMU coordinate frame on the rotating limb can be estimated by least squares when the subject's limb 1 is kept stationary and limb 2 is rotated. When a single limb (assume limb 2) is rotated about the joint center, the following equation holds:

[0122]

[0123] In equation (3.11), g is the acceleration due to gravity (9.8 m / s 2 ), Also represents the acceleration due to gravity. When only limb 2 is rotated about the joint, the optimization can be achieved using the Gauss-Newton iterative algorithm:

[0124]

[0125] Step S36: Calculate the joint angle using the IMU sensors based on a complementary filter.

[0126] A complementary filter is used to combine two measurements of a given signal, one consisting of high-frequency disturbance noise and the other consisting of low-frequency disturbance noise, to produce a single signal output measurement. Using filter gains, one can manipulate the dependence on each input and response time to drift correction, shortening the response time at the expense of greater output noise.

[0127] A two-stage complementary filter is used to obtain a combined direction estimate with smaller error components compared to using only a single sensor signal. In the following, only the Figure 3 The flexion angle is calculated from the acceleration and angular rate only. The gyroscope-based flexion angle can be calculated by integrating the difference in angular rate about the joint axis, i.e.,

[0128]

[0129] In addition, by knowing the joint axis coordinates, many constraint methods can be employed that require the sensor axis to coincide with the joint axis or segment axis, which can be extended to three-dimensional space. The measured acceleration can be transferred to the joint axis by:

[0130]

[0131] and are the same quantities measured in two different local coordinate systems, which are rotated with respect to each other around one axis. Therefore, the flexion angle can be expressed as and into the joint plane (ideally, the two angles are the same, but due to measurement inaccuracies, they are approximations). Therefore, for each local frame a pair of joint plane axes

[0132]

[0133] The accelerometer-based joint angle is calculated by

[0134]

[0135] where represents the angle between the two vectors in . Since no integration is employed to calculate it, the resulting α acc (t) is not affected by drift. The above equation is sensitive to measurement errors, if the shift acceleration, is almost collinear with the joint axis j1, j2. However, in almost all practical cases, the gravitational acceleration dominates the acceleration signals a1(t), a2(t) and Therefore, the error is significant only when the joint axis is close to vertical or during strong acceleration in the medial or lateral direction. While both of these cases are rare in walking and most other motions or daily activities. In addition, the term in equation (26) is usually small compared to the gravitational acceleration, thus the sensitivity to o1, o2 is low.

[0136] The gyro-based angle is very accurate for short periods of time, but exhibits a slow drift (the drift depends on the bias of the gyroscope). The accelerometer-based angle does not drift, but it is affected by the accelerometer noise and is not reliable at moments of large acceleration changes. Therefore, combining the two angles using standard sensor fusion tools (e.g., a complementary filter or a Kalman filter) can make the result more accurate. The result should be expressed in terms of α acc+gyr (t). For example, using a complementary filter, it is expressed as

[0137] α acc+gyr (t) = λ α acc (t) + (1 - λ) (α acc+gyr (t - Δt) + α gyr (t) - α gyr (t - Δt))#(29)

[0138] Where λ∈[0, 1].

[0139] Knee joint angle with noise but no drift angle alpha based on sensor fusion of accelerometer and gyroscope acc (t) with very accurate but drifting angle alpha gyr (t) using complementary filter to get angle alpha acc+gyr (t).

[0140] Step S4, use the sensor data in step S3 for fusion, the specific steps are as follows:

[0141] The present application aims at the problem that the depth camera often loses bone data due to occlusion and the like when used as a bone tracking system, and the problem that the joint angle measurement using IMU alone often drifts over time, and designs a quaternion-based extended Kalman filtering scheme to calibrate the data results with large errors or loss from the depth camera and IMU, and calculates the joint angle through fusion of the depth camera and IMU.

[0142] Step S41: unify the coordinate systems of the depth camera and the IMU, convert the depth camera coordinates to obtain the depth camera data information in the global coordinate system, and when the data of the IMU and the depth camera are mainly referenced to the global coordinate system, the data fusion of the two sensors can be performed. The steps of unifying the coordinate systems are as follows:

[0143]

[0144] In formula (30), X w , Y w , and Z w are the coordinates of the global coordinate system in (O w , x w , y w , z w ), X k , Y k , and Z k are the coordinates of the depth camera coordinate system in (O k , x k , y k , z k ), and X * , Y * , and Z * are the translation amounts of the global coordinate system and the depth camera coordinate system. Because the two coordinate systems are parallel to each other, only the global coordinate origin (O k , x k , y k , z k ) needs to be determined in the depth camera coordinate system (O w) position. The walking gait acquisition system is calibrated to obtain the translation amount, and the depth camera coordinates are converted by formula (30) to obtain the depth camera data information in the global coordinate system. When the data of the IMU and the depth camera are mainly referenced to the global coordinate system, the data fusion of the two sensors can be performed.

[0145] Step S42: combining Figure 4 An extended Kalman filter scheme based on quaternion is designed to calibrate the data results with large errors or missing data from the depth camera and the IMU, and to calculate the joint angle through the fusion of the depth camera and the IMU.

[0146] The quaternion q can be obtained by integrating the quaternion derivative, and the state model can be written as:

[0147]

[0148] Among them

[0149]

[0150] In the formula, ω = [ω x ω y ω z ] T represents the output of the three-axis gyroscope, and [ω×] is a 3×3 skew-symmetric matrix, which is given by the following formula:

[0151]

[0152] The present application uses an extended Kalman (EKF) to process data from the IMU and the depth camera. The state vector in the designed extended Kalman filter fusion scheme is composed of a quaternion and a knee joint angular velocity. The covariance matrix of the measurement noise is a 7×7 diagonal matrix, and the measurement vector has seven elements (four quaternions and three angular velocity components). The quaternion error component is estimated by combining the Euler angle and the quaternion transformation. The calculation formula of the quaternion obtained from the Euler angle is derived as follows:

[0153]

[0154]

[0155]

[0156]

[0157] In the formula, φ is the knee joint roll angle, θ is the knee joint heading angle, for the knee pitch angle, and Δφ is the error of the knee in the roll angle. The measurement noise matrix requires the values of the current inflection point spacing and the inclination angle. For the remaining diagonal values related to the angular velocity, the IMU gyroscope accuracy is known, and the accuracy of the depth camera is calculated according to the root mean square (RMS) error of the angular velocity.

[0158] In combination Figure 5 , the input from each individual sensor is a vector y The state vector consists of four components of the normalized quaternion and three components of the angular velocity. The extended Kalman filter processes the input parameters using the evolution model in the figure, the predicted state is used to update the measurement noise covariance matrices R1 and R2, the update matrix A and the process noise covariance matrix Q, and the state quaternion is normalized after each step.

[0159] Step S5, using the results obtained in steps S3, S4 to perform gait analysis.

[0160] Step S51: The results obtained in steps S3, S4 can be used to analyze the knee angle, gait symmetry test, disease prevention, rehabilitation evaluation, and other gait-related research.

[0161] Based on the same inventive concept, the present application also provides a gait analysis system that fuses a depth camera and an inertial sensor, in combination Figure 6 , comprising the following modules:

[0162] Human gait data acquisition module: the acquisition module designs a data acquisition system, runs an experimental acquisition program to collect sensor-related human activity data;

[0163] Data preprocessing module: preprocess the original human gait data;

[0164] Gait-related calculation module: calculate the human joint angle, bone length, and gait parameters;

[0165] Sensor fusion module: fuse the depth camera and IMU sensor data to calculate the joint angle.

[0166] The specific implementation methods of the above modules are the same as the aforementioned gait analysis method that fuses a depth camera and an inertial sensor, and will not be repeated here.

[0167] The above-described embodiments are merely specific embodiments of the present application, which are used to illustrate the technical solutions of the present application, but not to limit the present application. The protection scope of the present application is not limited to this. Although the present application is described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily think of changes to the technical solutions recorded in the foregoing embodiments, or make equivalent replacements to some technical features therein, within the technical scope disclosed by the present application. The modifications, changes or replacements do not make the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application, and should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A gait analysis method fusing depth camera and inertial sensors, characterized in that, The method comprises the following steps: Step S1: collecting human gait data, specifically: Step S11: designing a depth camera data collection system, obtaining human depth image coordinates through a depth camera skeletal tracking method, determining human joint positions, and thus collecting human gait-related data, the human gait-related data including step length, stride, step width, step speed, and bone length; Step S12: using an IMU sensor containing an accelerometer and a gyroscope to collect acceleration and angular velocity signals of the sensor on the human body; the IMU sensor is installed on two body segments connecting the measured human joints; Step S2: preprocessing the original human gait data, which has: Step S21: filtering the depth camera joint coordinate data obtained in step S11 using the One Euro filtering algorithm; Step S22: denoising the accelerometer and gyroscope signals collected by the IMU through a Butterworth low-pass filter; Step S3: using the data obtained in step S2 to perform gait-related calculations, specifically: Step S31: using the depth camera joint tracking data obtained in step S2 to calculate joint angles, bone length, and gait parameters; After collecting the spatial coordinate sequence of the human skeletal points and filtering the data, because the depth camera camera cannot be completely placed with the X axis parallel to the ground, and the zero view of the U axis cannot be completely regarded as the landing point, the landing point of the two feet cannot be directly determined according to the returned coordinates to determine the cycle; During walking, the human center point presents periodic oscillation in the up-down and left-right directions; the lowest point of each oscillation is selected as a special point, and the frame when the special point appears is the special frame; the gait cycle T is the time taken from the landing of one side of the heel to the landing of the same side of the heel again during walking, that is, the time required for the appearance of three consecutive special frames; T = t k+2 -t k (1) In formula (1), t k is the time of occurrence of the kth special frame, t k+2 is the time of occurrence of the k+2th special frame; The gait cycle can be calculated by finding the position of the special frame; Bone length measurement is to calculate the distance between two joints connected by a bone; three-dimensional coordinate information of skeletal points is obtained through a depth camera, a joint vector is generated, and bone length is calculated; in the 3D coordinate system of the depth camera, the 3D coordinates of adjacent joints of the bone are retrieved using the depth camera, and the following formula (2) is used to calculate the vector of each segment using quaternions: where E(X e ,Y e ,Z e ), F(X f ,Y f ,Z f ) represent the positions of two adjacent joints, respectively. Then, the bone length is calculated: Step length is the distance from the position where one side of the heel lands to the position where the other side of the heel lands next during walking, which is expressed using the distance formula of spatial point coordinates: In formula (4), s pi represents the step length of the i-th step, X k , Y k , Z k represents the spatial coordinates when the special point sequence number is k, X k+1 , Y k+1 , Z k+1 represents the spatial coordinates when the special point sequence number is k+1; Stride is the distance from the position where one side of the heel lands to the position where the same side of the heel lands next during walking, which is expressed using the distance formula of spatial point coordinates: s di = s pi + s p(i+1) (5) Step width is the distance from the landing point of the other side of the foot to the plane established by the two adjacent landing points of the right foot and a hovering point; the above three points through which the right foot moves are selected to establish an empty space, and the two landing points and one hovering point are respectively R1(X1, Y1, Z1), R2(X2, Y2, Z2), and R3(X3, Y3, Z3); the space plane equation is: AX+BY+CZ=D (7) The distance from the left foot landing point to the plane, i.e. the step width, is expressed as follows: where (X0, Y0, Z0) is the coordinate of the right foot landing point; Step speed is the instantaneous speed of the captured bone point in each data collection; the average speed of the adjacent three points is used to replace the instantaneous speed of the central point, and the instantaneous speed of the central point is: where f c is the frequency of the depth camera acquisition, δ t-1 is the distance between point B t-1 and point B t . Acceleration: the acceleration is the change rate of adjacent frames, i.e.: From the m+1 acquired point coordinates, the time parameter sequence, i.e. the displacement sequence X1, X2, X3,... X m , the velocity sequence V1, V2, V3,... V m-1 , and the acceleration sequence a1, a2, a3,... a m-3 , can be determined. The depth camera uses its depth sensor to retrieve information in 3D coordinates; it uses its own position as a reference point to retrieve the joint points on the skeleton; the depth camera uses the relative difference between the joint points and the reference point to calculate the three-axis coordinates of the human body and the motion tracking information and the corresponding 3D joint angles; since the depth camera can only provide a set of 3D coordinates for the human body limbs, it is necessary to calculate the joint angles according to the 3D coordinates of the back-end system; The three-dimensional coordinate information of the bone points is obtained by the depth camera, the joint vector is generated, and then the joint control angle is obtained by calculating the included angle of the joint vector; in the 3D coordinate system of the depth camera, the 3D coordinates of the adjacent joints of the measured joint are retrieved using the depth camera, and the vectors of each segment are calculated using quaternions and formula (11): Then, the joint angle is calculated using the following formula: cos -1 (cosθ) (13) Step S32: using the data of the IMU sensor to identify the joint axis and the position of the joint axis; Without assuming that any local sensor axis coincides with the joint axis or the limb segment, the direction and position of the joint flexion and extension axis are arbitrary; while the joint moves in an arbitrary manner around its degrees of freedom, the identification data is collected, and a data set S(i) is obtained every 0.1s, which is in the form of: where a1(t), a2(t) e R 3 is the acceleration at some sampling period Δt, g1(t), g2(t) e R 3 is the angular velocity, is the angular rate; the time between the two data sets must be a multiple of the sampling period, i.e. Thus, N»1 data sets are collected, which will be used in the subsequent sections to identify the local joint axes and position coordinates; Step S33: identifying the joint axis coordinates based on the results obtained in step S32; Data set S(i), i e [1, N], for identifying unit length direction vectors Joint flexion-extension axis in local coordinates of the two sensors; j1 and j2 are constants, depending only on the mounting direction of the sensors with respect to the joint; angular rates g1(t), g2(t) differ only in that the joint angular velocity vector and the rotation matrix are measured on the joint; therefore, they project into the joint plane with the same length for each time instant, i.e.: where ||·||2 represents the Euclidean norm; this constraint is valid regardless of the position and direction of the sensor mounted on the segment; each data set S(i), i∈[1, N] must satisfy (15); therefore, j1 and j2 are identified by minimizing the left side of formula (15) of all data sets; j1 and j2 are represented in spherical coordinates as: j1= (cos(φ1)cos(θ1), cos(φ1)sin(θ1), sin(φ1)) T (16) j2= (cos(φ2)cos(θ2), cos(φ2)sin(θ2), sin(φ2)) T (17) and the sum of squared errors is defined as: where φ1, φ2 represent the roll angles of the joint axis, and θ1, θ2 represent the heading angles of the joint axis; since the sign of formula (18) relative to j1 and j2 is unchanged, this cost function has four minimum values, corresponding to the four possible combinations of signs, (j1, j2), (-j1, j2), (j1, -j2), and (-j1, -j2); Step S34: performing sign matching on the joint axis coordinates in step S33; Step S35: identifying the coordinates of the joint position obtained in step S32; Step S36: using the IMU sensor to calculate the joint angle based on the complementary filter; Step S4: using the sensor data in step S3 for fusion, specifically: Step S41: unifying the coordinate systems of the depth camera and the IMU, converting the depth camera coordinates to obtain depth camera data information in the global coordinate system; when the data of the IMU and the depth camera are mainly referenced in the global coordinate system, the data of the two sensors are fused; the steps of unifying the coordinate systems are as follows: In formula (30), X w , Y w , and Z w are coordinates of the global coordinate system in (O w , x w , y w , z w ), X k , Y k , and Z k are coordinates of the depth camera coordinate system in (O k , x k , y k , z k ), and X * , Y * , and Z * are translation amounts of the global coordinate system and the depth camera coordinate system coinciding with each other; because the two coordinate systems are parallel to each other, only the position of the global coordinate origin (O k ) needs to be determined in the depth camera coordinate system (O k , x k , y k , z w ); the walking gait acquisition system is calibrated to obtain the translation amount, the depth camera coordinates are converted by using formula (30) to obtain the depth camera data information in the global coordinates; and when the data of the IMU and the depth camera are mainly in the global coordinate system as the reference system, the data fusion of the two sensors can be performed. Step S42: using the quaternion-based extended Kalman filter method, calibrating the data results from the depth camera and the IMU, calculating the joint angle through the fusion of the depth camera and the IMU; The quaternion q is obtained by integrating the quaternion derivative, and the state model is written as: Wherein where ω = [ωxωyωz]T x ω y ω z ] T denotes the output of the tri-axial gyroscope, [ωx] is a 3 x 3 skew-symmetric matrix given by: The data from the IMU and the depth camera are processed using the extended Kalman filter; the state vector in the designed extended Kalman filter fusion scheme is composed of the quaternion and the knee joint angular velocity; the covariance matrix of the measurement noise is a 7*7 diagonal matrix, and the measurement vector has seven elements: four components of the normalized quaternion and three components of the angular velocity; the quaternion error component is estimated by combining the Euler angle and the quaternion transformation; the calculation formula of the quaternion is obtained from the Euler angle, and part of the formula is derived as follows: where φ is the knee roll angle, θ is the knee yaw angle, is the knee pitch angle, and Δφ is the knee roll angle error; the measurement noise matrix requires the current knee pitch and roll values; for the remaining diagonal values related to angular velocity, the IMU gyroscope accuracy is known, and the depth camera accuracy is calculated from the root mean square error of angular velocity; The input from each individual sensor is a vector y, and the state vector is composed of four components of the normalized quaternion and three components of the angular velocity; the extended Kalman filter uses the evolution model in the figure to process the input parameters, and the predicted state is used to update the measurement noise covariance matrix R1 and R2, the update matrix A and the process noise covariance matrix Q, and the state quaternion is normalized after each step; Step S5: using the results obtained in steps S3 and S4 for gait analysis.

2. A gait analysis system fusing a depth camera with inertial sensors, characterized in that, The system for implementing the method of claim 1 comprises: A human gait data acquisition module: for acquiring gait-related human activity data of the sensor; the human gait data acquisition module obtains human depth image coordinates by a depth camera bone tracking method, determines the position of human joints, and thus acquires gait-related human data, which includes step length, stride, step width, step speed and bone length; an IMU sensor containing an accelerometer and a gyroscope is used to acquire acceleration and angular velocity signals of the sensor on the human body; the IMU sensor is installed on two body segments connected to the measured human joints; A data preprocessing module: for preprocessing the original human gait data; the data preprocessing module filters the depth camera joint coordinate data using the One Euro filter algorithm; the accelerometer and gyroscope signals collected by the IMU are denoised by a Butterworth low-pass filter; A gait-related calculation module: for calculating human joint angles, bone lengths and gait parameters; the specific implementation of the gait-related calculation module is as follows: The depth camera joint tracking data is used to calculate the joint angle, bone length and gait parameter; The data of the IMU sensor are used to identify the joint axis and the position of the joint axis, and the result is identified as the joint axis coordinates; The joint axis coordinates are symbol matched; The joint position is identified as coordinates; Based on the complementary filter, the IMU sensor is used to calculate the joint angle; The sensor fusion module is used for fusing sensor data to obtain fused joint angles, and the specific implementation manner of the sensor fusion module is as follows: unifying the coordinate systems of the depth camera and the IMU, converting the depth camera coordinates to obtain depth camera data information in the global coordinate system, fusing data of the two sensors when the data of the IMU and the depth camera are mainly referenced in the global coordinate system, and adopting an extended Kalman filtering method based on a quaternion to calibrate data results from the depth camera and the IMU, and calculating the joint angle through fusion of the depth camera and the IMU.

Images