A method for detecting abnormal motion trajectory by fusing multi-source positioning data

Abstract

The application relates to the technical field of sports health monitoring, and discloses a sports trajectory anomaly detection method fusing multi-source positioning data, which comprises the following steps: synchronously collecting ultra-wideband ranging values, double-foot inertia data and foot bottom pressure touch ground states; constructing a zero-speed virtual constraint for a foot according to foot bottom pressure to identify a touch ground moment; calculating a theoretical displacement of a waist label from a supporting foot based on a lower limb motion chain subjected to joint limiting; integrating the ranging, inertia pre-integration, zero-speed and motion chain constraint into a factor graph optimization solution space to solve a space state; projecting a trajectory to a local coordinate system of a ski blanket; and determining motion abnormality through a phase diagram and a dynamic residual error. The application can distinguish between positioning jumps caused by ultra-wideband non-line-of-sight errors and real action out-of-control of athletes by integrating limb motion chain constraints based on muscle-skeletal inverse dynamics into factor graph fusion positioning and constructing consistency verification by using gait touch zero-speed states.

Description

Technical Field

[0001] This invention relates to the field of sports health monitoring technology, specifically to a method for detecting abnormal motion trajectories by integrating multi-source positioning data. Background Technology

[0002] In professional sports training using indoor ski simulators, wearable sensor networks are typically used to acquire athletes' kinematic data. Commonly used sensing methods include ultra-wideband positioning systems (UWBS), inertial measurement units (IMUs), and plantar pressure sensors. UWBS calculates distance by measuring the bidirectional flight time of electromagnetic pulses between a tag and a fixed base station, and then calculates spatial coordinates to achieve high temporal resolution positioning and tracking. IMUs integrate a three-axis accelerometer and a three-axis gyroscope, which can independently measure the linear acceleration and angular velocity of the vehicle, and after integration, can continuously calculate the vehicle's attitude and displacement changes. Plantar pressure sensors are distributed in an array within the insole, collecting pressure values ​​in different areas of the foot in real time to determine whether the foot is in contact with the support surface and the distribution of force.

[0003] In the confined environment of an indoor ski simulator, athletes stand on a reciprocatingly rolling ski mat and perform specialized maneuvers such as turns and carving. Their spatial position relative to fixed indoor reference points changes repeatedly within a limited range. However, when athletes perform carving postures such as large body tilts and deep squats, their torso and lower limbs create a physical obstruction structure against the positioning tag worn on their waist. This obstruction causes the ultra-wideband signal to be blocked in its propagation path or to reach the receiving end after reflection and diffraction, resulting in non-line-of-sight ranging errors. This causes outliers in the calculated real-time trajectory, which are commonly referred to in the industry as "outliers" in the positioning trajectory.

[0004] Existing motion trajectory anomaly detection methods, when applied to general scenarios such as pedestrian navigation and vehicle positioning, largely rely on analyzing abrupt changes in the geometric shape, sudden speed shifts, or deviations from preset geofence boundaries in the positioning trajectory to determine whether an anomaly has occurred. These methods assume that the absolute spatial coordinates collected by sensors faithfully reflect the actual motion state of the vehicle. However, in indoor ski simulator training scenarios, the simulator's own motion couples the absolute position reference frame with the athlete's movements, and frequent body occlusion continuously introduces non-line-of-sight errors into the positioning data. Existing methods struggle to distinguish whether coordinate deviations in the positioning trajectory stem from measurement distortion caused by signal occlusion or from actual loss of balance and movement errors by the athlete. This introduces numerous false alarms due to sensor interference into the anomaly detection results, limiting the accurate capture of disturbances in specific sports techniques. Therefore, the following solution is proposed to address these issues. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides a method for detecting motion trajectory anomalies by fusing multi-source positioning data, comprising the following steps: The system collects raw ranging values ​​from the ultra-wideband positioning system, acceleration and angular velocity data from the inertial measurement unit, and plantar pressure sensor data, and synchronizes them in time. Using plantar pressure sensor data to identify the zero-velocity state upon ground contact and extracting zero-velocity virtual observation constraints; Based on the lower limb bone length calibration value and the angular velocity data of the inertial measurement unit, a limb kinematic chain model constrained by the joint angle range is constructed, and the spatial displacement of the waist relative to the supporting foot is derived. The original ranging value, the pre-integrated quantity of the inertial measurement unit, the zero-velocity virtual observation constraint, and the spatial displacement are incorporated into the factor graph for nonlinear optimization to solve the waist positioning trajectory; The waist positioning trajectory is projected onto the local coordinate system formed by the propulsion direction and normal vector of the moving vehicle to obtain the tangential displacement and normal displacement. Based on the phase diagram closure degree of the tangential displacement and normal displacement, as well as the residual statistical characteristics between the spatial displacement and the displacement increment obtained by factor diagram optimization, the motion mode is determined to be abnormal.

[0006] Preferably, identifying the zero-velocity state upon ground contact includes: processing plantar pressure sensing data using a continuous sampling window; when the proportion of ground contact marker frames within the window to the total number of frames is greater than 0.9, and the pressure value of the first frame of the window increases by more than 30% compared to the pressure value of the previous frame, the corresponding foot is determined to be in a zero-velocity state, and a velocity observation residual is established. In the formula, For a moment The zero-velocity virtual observation residual vector; For a moment feet The instantaneous velocity vector to be estimated; It is a three-dimensional zero vector; For physical time variables.

[0007] Preferably, the lower limb bone length calibration value is obtained in the following way: when the subject is in a zero-speed state of supporting the body with one leg upright, the pitch angle of the supporting leg is obtained by the inertial measurement unit, combined with the pre-input height parameters, the initial estimated values ​​of thigh length and calf length are calculated according to the human leg-to-body ratio coefficient of 0.45 to 0.53, and then the thigh length calibration value and calf length calibration value are obtained by solving the geometric relationship of limb projection in the supporting state. The joint angle range constraints are as follows: the range of motion of the knee joint flexion and extension angle is limited to 0° to 130°, and the range of motion of the hip joint flexion and extension angle is limited to -15° to 120°. The range of motion is a hard boundary preset based on the physiological limit angle of the human lower limb skeletal structure in the sagittal plane without dislocation or fracture. When deriving the spatial displacement of the waist relative to the supporting foot, penalties are imposed on the state variables that exceed the range of motion.

[0008] Preferably, the penalty is implemented by constructing a cost function, the specific form of which is: in the factor graph optimization objective function, the knee joint angle is... Introducing penalty items Regarding the hip joint angle Introducing penalty items ,in This is a function of the square of the one-sided excess; it takes the value of zero when the angle is within the range of motion and the square of the excess when it exceeds the boundary; penalty coefficient. and The method for determining the penalty coefficient is as follows: based on the angular velocity noise variance of the inertial measurement unit and the passive resistance torque characteristics of human joints at extreme angles, the penalty coefficient is set so that the increment of the penalty term generated when the angle exceeds the constraint boundary by 1° is on the same order of magnitude as the typical variance value of the zero-velocity observation residual, thereby ensuring the reasonable weight of physiological constraints in optimization.

[0009] Preferably, the factor plot includes: UWB ranging factor, whose residual form is as follows: In the formula, For the reason about the Ranging residual scalar of a UWB base station; For waist label at all times The position vector, For the first The location vectors of each base station. The UWB ranging value in step S1; And the kinematic chain consistency factor, whose residual is the difference between the spatial displacement and the position increment of the waist state node at the current moment relative to the waist state node at the last zero-velocity update moment.

[0010] Preferably, a kernel function based on kinematic conflict is introduced into the UWB ranging factor for weight reduction processing. The mathematical expression of the kernel function is: In the formula, The residual vector of the UWB ranging factor; Here is the covariance matrix of the UWB ranging factor; The Mahalanobis distance; Threshold; threshold The method for determining the value is as follows: the Mahalanobis distance is regarded as a chi-square distribution variable, and its degree of freedom is equal to the upper 95th percentile of the number of effective base stations of the UWB ranging factor at the current time as the threshold. In this way, abnormal ranging values ​​that exceed the sensor noise level are statistically judged as kinematic conflicts and their weights are reduced in the form of an inverse proportional function.

[0011] Preferably, solving the waist positioning trajectory includes: performing nonlinear least squares optimization on the factor map, and iteratively estimating the position, velocity, attitude of the waist tag, and the zero bias parameters of the inertial measurement unit.

[0012] Preferably, the local coordinate system uses the propulsion direction of the moving vehicle and its normal vector as coordinate axes; the waist positioning trajectory is projected onto the local coordinate system to obtain the time series of tangential displacement and normal displacement.

[0013] Preferably, determining motion pattern anomalies based on the coupling characteristics of projected displacement and spatial displacement residuals includes: Calculate the phase diagram closure, the phase diagram is determined by tangential displacement. and normal displacement The calculation method for phase diagram closure is as follows: extract the trajectory line segment of a complete action cycle in the phase diagram, calculate the Euclidean distance between the first and last endpoints of the line segment, and use this distance as the closure index. A preset envelope is established by collecting data from several normal motion patterns. - Phase diagram trajectory, calculated by dividing the phase plane into angular bins and statistically averaging the radial distances. and standard deviation ,by It constitutes the normal motor envelope; The non-stationarity test is performed on the residual sequence of the motion chain consistency factor. The statistical test method for non-stationary divergence is as follows: the augmented Dickey-Fowler test is performed on the residual sequence. If the test statistic is greater than the critical value, the sequence is determined to be non-stationary. At the same time, the sliding window variance of the residual sequence is calculated. When the sliding variance increases continuously and exceeds a predetermined multiple of the initial stage variance, the residual is determined to be non-stationary divergence. An abnormal event is defined as when the phase diagram closure exceeds the preset closure threshold, the phase diagram trajectory exceeds the normal motion envelope, and the residual is simultaneously determined to be non-stationary divergence.

[0014] The present invention has the following beneficial effects: This invention introduces musculoskeletal inverse dynamics constraints into a multi-source localization fusion framework, enabling the detection of localization outliers caused by body occlusion during indoor ski simulator training. It utilizes gait ground contact events to trigger zero-velocity correction and kinetic chain consistency verification, distinguishing between sensor interference and actual motion anomalies from the underlying signal level, thus avoiding misjudgments of occlusion errors. The fusion of two-dimensional signature analysis in a projected coordinate system allows anomaly detection to focus on stable changes in motion control rather than isolated trajectory point offsets, better aligning with the needs of specialized sports assessment. This method does not rely on absolute trajectory morphology and is sensitive to relative motion in constrained environments, broadening the adaptability of localization anomaly detection in professional sports training scenarios. Attached Figure Description

[0015] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 This is a flowchart illustrating a motion trajectory anomaly detection method that integrates multi-source positioning data according to the present invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Please see Figure 1 As shown, this invention is a method for detecting motion trajectory anomalies by fusing multi-source positioning data. The detection method includes the following steps: The system collects raw ranging values ​​from the ultra-wideband positioning system, acceleration and angular velocity data from the inertial measurement unit, and plantar pressure sensor data, and synchronizes them in time. Using plantar pressure sensor data to identify the zero-velocity state upon ground contact and extracting zero-velocity virtual observation constraints; Based on the lower limb bone length calibration value and the angular velocity data of the inertial measurement unit, a limb kinematic chain model constrained by the joint angle range is constructed, and the spatial displacement of the waist relative to the supporting foot is derived. The original ranging value, the pre-integrated quantity of the inertial measurement unit, the zero-velocity virtual observation constraint, and the spatial displacement are incorporated into the factor graph for nonlinear optimization to solve the waist positioning trajectory; The waist positioning trajectory is projected onto the local coordinate system formed by the propulsion direction and normal vector of the moving vehicle to obtain the tangential displacement and normal displacement. Based on the phase diagram closure degree of the tangential displacement and normal displacement, as well as the residual statistical characteristics between the spatial displacement and the displacement increment obtained by factor diagram optimization, the motion mode is determined to be abnormal.

[0019] Identifying the zero-velocity state upon ground contact involves processing plantar pressure sensor data using a continuous sampling window. When the proportion of ground contact marker frames within the window to the total number of frames is greater than 0.9, and the pressure value in the first frame of the window increases by more than 30% compared to the pressure value in the previous frame, the corresponding foot is determined to be in a zero-velocity state, and a velocity observation residual is established. In the formula, For a moment The zero-velocity virtual observation residual vector; For a moment feet The instantaneous velocity vector to be estimated; It is a three-dimensional zero vector; For physical time variables.

[0020] The lower limb bone length calibration values ​​are obtained as follows: with the subject maintaining a zero-speed state of single-leg upright support, the pitch angle of the supporting leg is obtained by the inertial measurement unit. Combined with the pre-input height parameters, the initial estimated values ​​of thigh length and calf length are calculated according to the human leg-to-body ratio coefficient of 0.45 to 0.53. Then, the thigh length calibration value and calf length calibration value are obtained by solving the geometric relationship of limb projection in the support state. The joint angle range constraints are as follows: the range of motion of the knee joint flexion and extension angle is limited to 0° to 130°, and the range of motion of the hip joint flexion and extension angle is limited to -15° to 120°. The range of motion is a hard boundary preset based on the physiological limit angle of the human lower limb skeletal structure in the sagittal plane without dislocation or fracture. When deriving the spatial displacement of the waist relative to the supporting foot, penalties are imposed on the state variables that exceed the range of motion.

[0021] The penalty is implemented by constructing a cost function, which takes the following form: in the factor graph optimization objective function, the knee joint angle is... Introducing penalty items Regarding the hip joint angle Introducing penalty items ,in This is a function of the square of the one-sided excess; it takes the value of zero when the angle is within the range of motion and the square of the excess when it exceeds the boundary; penalty coefficient. and The method for determining the penalty coefficient is as follows: based on the angular velocity noise variance of the inertial measurement unit and the passive resistance torque characteristics of human joints at extreme angles, the penalty coefficient is set so that the increment of the penalty term generated when the angle exceeds the constraint boundary by 1° is on the same order of magnitude as the typical variance value of the zero-velocity observation residual, thereby ensuring the reasonable weight of physiological constraints in optimization.

[0022] The factor diagram includes: UWB ranging factor, whose residual form is... In the formula, For the reason about the Ranging residual scalar of a UWB base station; For waist label at all times The position vector, For the first The location vectors of each base station. The UWB ranging value in step S1; And the kinematic chain consistency factor, whose residual is the difference between the spatial displacement and the position increment of the waist state node at the current moment relative to the waist state node at the last zero-velocity update moment.

[0023] The UWB ranging factor incorporates a kernel function based on kinematic conflict for weight reduction processing. The mathematical expression of the kernel function is as follows: In the formula, The residual vector of the UWB ranging factor; Here is the covariance matrix of the UWB ranging factor; The Mahalanobis distance; Threshold; threshold The method for determining the value is as follows: the Mahalanobis distance is regarded as a chi-square distribution variable, and its degree of freedom is equal to the upper 95th percentile of the number of effective base stations of the UWB ranging factor at the current time as the threshold. In this way, abnormal ranging values ​​that exceed the sensor noise level are statistically judged as kinematic conflicts and their weights are reduced in the form of an inverse proportional function.

[0024] Solving the waist localization trajectory involves: performing nonlinear least squares optimization on the factor map, and iteratively estimating the position, velocity, attitude of the waist tag, and the zero bias parameters of the inertial measurement unit.

[0025] The local coordinate system uses the propulsion direction of the moving vehicle and its normal vector as coordinate axes; the waist positioning trajectory is projected onto this local coordinate system to obtain the time series of tangential and normal displacements.

[0026] Based on the coupling characteristics of projected displacement and spatial displacement residuals, motion pattern anomalies are identified, including: Calculate the phase diagram closure, the phase diagram is determined by tangential displacement. and normal displacement The calculation method for phase diagram closure is as follows: extract the trajectory line segment of a complete action cycle in the phase diagram, calculate the Euclidean distance between the first and last endpoints of the line segment, and use this distance as the closure index. A preset envelope is established by collecting data from several normal motion patterns. - Phase diagram trajectory, calculated by dividing the phase plane into angular bins and statistically averaging the radial distances. and standard deviation ,by It constitutes the normal motor envelope; The non-stationarity test is performed on the residual sequence of the motion chain consistency factor. The statistical test method for non-stationary divergence is as follows: the augmented Dickey-Fowler test is performed on the residual sequence. If the test statistic is greater than the critical value, the sequence is determined to be non-stationary. At the same time, the sliding window variance of the residual sequence is calculated. When the sliding variance increases continuously and exceeds a predetermined multiple of the initial stage variance, the residual is determined to be non-stationary divergence. An abnormal event is defined as when the phase diagram closure exceeds the preset closure threshold, the phase diagram trajectory exceeds the normal motion envelope, and the residual is simultaneously determined to be non-stationary divergence.

[0027] The specific application of this embodiment is as follows: Step S1: Time synchronization and acquisition of multi-source data Before training began, a UWB tag was placed on the center of the back of the subject's waist, and an IMU integrating a three-axis accelerometer and a three-axis gyroscope and a foot membrane pressure sensor array were placed under the insoles of the left and right feet, respectively.

[0028] S101: Synchronous Trigger Acquisition; After the ski simulator starts, the motor start pulse signal issued by the simulator controller simultaneously triggers the UWB positioning base station, UWB tag, IMU, and pressure sensor to start recording data; Define the physical time axis. ; S102: Acquire multi-source observation data; in each sampling period Internally, the following data is collected synchronously: The bidirectional flight time from the tag to N fixed base stations in the UWB positioning system is collected and denoted as . ,in , Based on this, the label can be calculated in real time up to the 1st. The original ranging values ​​of each base station ; IMU measurements, including triaxial acceleration, were acquired at the left and right feet. and triaxial angular velocity ; The pressure values ​​of the thin-film pressure sensor arrays on the soles of the left and right feet were collected, and the foot dynamic phase identifier of the left foot was extracted respectively. and right foot phase identifier The value can be 1 (ground support state) or 0 (swinging and airborne state).

[0029] Step S2: Zero-velocity correction and local velocity extraction based on gait ground contact events This step aims to provide a high-confidence zero-velocity virtual observation for subsequent factor plots by utilizing precise moments of foot contact and stillness.

[0030] S201: Utilizing plantar pressure data and It identifies the precise time interval between the transition of any foot from the swing phase to the support phase; it defines a length of... A continuous sampling window, when the window is within ( for or The proportion of ) exceeds the threshold Furthermore, if the pressure gradient jump condition is met in the first frame of the window, the foot is determined to be in a stable ground-touching zero-velocity state. The following is the process for determining the jump condition: Define window length and ground contact ratio threshold Specific values: The system is set to a sampling frequency of 100%. Hz. According to general biomechanical statistics, the duration of the support phase in normal human movement is typically greater than 100 milliseconds. To avoid misjudgments triggered by transient noise in the plantar pressure sensor, a continuous sampling window length is set. This corresponds to a 100-millisecond time window.

[0031] Set ground contact ratio threshold This value is based on the physical conductivity characteristics of the plantar pressure sensor array: when the pressure sensor array is stepped on, due to the arch curvature of the foot or slight deformation of the insole, the sensing points at the edge of the pressure may vibrate. A 90% overwhelming proportion confirms that the sole of the foot is fully bearing the body weight and is in a state of full-foot contact with the ground. The formula is expressed as: In the formula, For the feet exist Gait phase indicators at any given moment (1 indicates ground contact, 0 indicates airborne). For foot identifiers; This represents the start time of the current detection window.

[0032] Define the calculation method for the pressure gradient abrupt change condition and the corresponding threshold.

[0033] Define the total plantar pressure value as This value is determined by the foot. The pressure values ​​at all sensing points in the plantar diaphragm pressure sensor array are summed to obtain the result. The current time is then calculated. absolute value of time-series pressure gradient for: In the formula, The time interval between two adjacent sampling points; Sampling time Total plantar pressure value; Sampling time The total pressure value of the plantar surface.

[0034] The pressure gradient abrupt change condition is defined as follows: within three consecutive sampling points, Continuously exceeding the preset jump threshold Jump threshold System calibration is performed based on the sensor's no-load fluctuation characteristics: With the sensor in a no-load state, data is collected at 100 sampling points per second, and the mean absolute value of the pressure gradient is calculated. and standard deviation Set the jump threshold This threshold of 6 standard deviations statistically ensures that the false trigger rate of ground contact events approaches zero.

[0035] The logical condition for a jump is: In the formula, For traversing variables; This determination formula directly addresses the physical force abrupt change of the pressure sensor, defining the precise starting point of the irreversible physical contact event from when the foot is completely airborne to when it impacts the ground.

[0036] S202: After confirming the foot Extract accelerometer data from its IMU during the period when it is in a zero-speed state. and gyroscope data Construct a zero-velocity virtual observation constraint for this time period; theoretically, the velocity of the foot in the global frame should be a zero vector at this time; establish a velocity observation residual model: In the formula, For a moment The zero-velocity virtual observation residual vector; For a moment feet The instantaneous velocity vector to be estimated; It is a three-dimensional zero vector; For physical time variables.

[0037] Step S3: Construct virtual limb kinematic chain constraints based on musculoskeletal inverse dynamics This step utilizes the angular velocity information from the IMU and, based on a simplified lower limb skeletal model, derives the spatial dynamic constraints of the waist label relative to the supporting foot.

[0038] S301: Establish a simplified kinetic chain model of the human lower limbs; during system initialization, calibrate and obtain the thigh length of the subject. and calf length During the stable ground contact support phase, the center of the support foot is used as the reference point. Utilizing the angular velocity of its IMU (according to choose or (corresponding data), and the posture of the supporting leg is recursively deduced through quaternion integration; At the time of first use of this system, the subject's thigh length and calf length The calibration process is as follows: The subject stood naturally, and their height was measured externally. According to ergonomic statistics, the ratios of thigh length and calf length to height range from 0.23-0.27 and 0.22-0.26 respectively, resulting in a combined leg-to-body ratio coefficient between 0.45 and 0.53. Therefore, the initial estimate of thigh length can be calculated using the following formula. and initial estimate of calf length : In the formula, The thigh proportion coefficient ranges from 0.23 to 0.27. The calf ratio coefficient ranges from 0.22 to 0.26, and the sum of the two values ​​satisfies the leg-to-body ratio range.

[0039] The subject performed a single-leg upright support movement on a ski simulator, entering a zero-speed state. At this moment, the spatial vector between the foot of the supporting leg and the waist tag is... It can be accurately obtained by a UWB positioning system under unobstructed conditions. Simultaneously, the pitch angle output by the support foot IMU... This reflects the degree of inclination of the leg relative to the direction of gravity. Using the geometric relationship of the right-angled triangle projected by the supporting leg, an equation is established in the sagittal plane: In the formula, and These are the angles between the thigh and lower leg relative to the direction of gravity, which can be determined by the IMU attachment position or attitude recursion.

[0040] In a single-leg standing support posture, because the knee joint is approximately in a straight and locked state, there is The above relationship can be simplified to: By combining the measured total length constraint with the initial estimated ratio, the final thigh length calibration value is calculated. and lower leg length calibration value The data is stored in the system for subsequent forward kinematics derivation. S302: Introducing muscle tension line constraints to limit the range of joint angles; during the integral solution of knee and hip joint angles, the calculated knee flexion-extension angles... and hip flexion-extension angle Apply hard boundary constraints based on anatomical tension lines: In the formula, The minimum physiological flexion angle of the knee joint is determined based on the individual lower limb muscle tension line anatomical characteristics of the subject. Negative or smaller values ​​represent hyperextension limit. The maximum physiological flexion angle of the knee joint is the limit angle when the heel is close to the buttocks. For a moment The hip flexion-extension angle, which is the angle between the longitudinal axis of the trunk and the longitudinal axis of the thigh in the sagittal plane; This is the minimum physiological flexion angle of the hip joint, usually a negative value, representing the extension limit; The maximum physiological flexion angle of the hip joint is the limit at which the thigh is raised forward and brought close to the abdomen. This constraint, as a customized cost function in the factor graph, imposes an exponential penalty on state variables that exceed the boundary, avoiding the calculation of anti-joint poses that violate human physiological structure due to IMU drift or NLOS outliers in UWB; this step is not a subjective rule of the human body, but an objective model constraint based on human anatomical parameters. S303: Calculate the virtual physical quantity of lumbar displacement; based on positive kinematics, using the constrained joint angles, derive the lumbar label point relative to the support base point. Theoretical spatial displacement vector This displacement vector contains the motion details sensed by the IMU and is unaffected by UWB signal occlusion.

[0041] Step S4: Multi-source factor graph fusion localization and spatiotemporal consistency conflict detection This step integrates the raw UWB ranging values, IMU pre-integrations, zero-velocity observations, and limb kinematic constraints into a factor graph framework for nonlinear optimization.

[0042] S401: Construct the state variable nodes to be optimized; the system's state vector. Location including waist label ,speed Posture Quaternions And the zero-bias IMU accelerometers and gyroscopes of the left and right feet; S402: Add each factor to the factor graph: UWB ranging factor: The factor is constructed using preprocessed ranging values ​​only when the conditions of stable ground contact and no deep occlusion are met; for each base station Its residual is: In the formula, For the reason about the Ranging residual scalar of a UWB base station; For waist label at all times The position vector, For the first The location vectors of each base station. This refers to the UWB ranging value in step S1; here, a kernel function based on kinematic conflict is introduced. When the Mahalanobis distance exceeds the threshold, the kernel function automatically reduces the weights to suppress non-line-of-sight errors. The mathematical expression for the kernel function is: In the formula, The residual vector of the UWB ranging factor; The covariance matrix of the UWB ranging factor; The Mahalanobis distance; Threshold; threshold The determination method is as follows: the Mahalanobis distance is regarded as a chi-square distribution variable, and its degree of freedom is equal to the upper 95th percentile of the number of effective base stations of the UWB ranging factor at the current time as the threshold. In this way, abnormal ranging values ​​that exceed the sensor noise level are statistically judged as kinematic conflicts and their weights are reduced in the form of an inverse proportional function. IMU pre-integration factor: Pre-integrates the IMU gyroscope and accelerometer measurements between two consecutive frames to construct a relative pose change factor, and connects the residuals to the state nodes at adjacent time points. and ; Zero-rate correction factor: generated in embedding step S2 ; Kinematic chain consistency factor: The muscle dynamics predicted displacement obtained in step S3 As a high-confidence relative displacement constraint, it is added between adjacent waist state nodes, and its residual term is: In the formula, This is the three-dimensional residual vector of the muscle kinetic chain consistency factor; For at any time The most recent zero-speed reference moment before ground contact Below, the optimized 3D position vector of the waist label; The timestamp of the last stable foot contact with the ground at zero speed is a specific anchor point moment in the gait cycle; S403: Perform nonlinear least squares optimization on the entire factor graph, and iteratively solve for the optimal state vector. .

[0043] Step S5: Anomaly detection based on inverse dynamic signature of projection residuals After obtaining the high-precision fused trajectory in step S4, instead of directly using the absolute coordinates of the trajectory points, anomaly detection is performed in the projection space.

[0044] S501: Project the spatial trajectory of the waist tag calculated in step S4 onto the local coordinate system formed by the propulsion direction of the ski carpet in the ski simulator and its normal vector; define the unit vector of the propulsion direction as... The normal vector is Extract the tangential displacement at each moment. and normal displacement : S502: Calculate the full kinematic signature in the projected coordinate system—that is... and Phase diagram closure and in step S3 The residual statistics of the final optimization result; the coupled analysis of this two-dimensional physical quantity (projection plane displacement and internal dynamic residual) constitutes the anomaly judgment rule of this method; when the phase diagram trajectory exceeds the preset normal rotation envelope range and the internal dynamic residual is in a non-stationary divergent state, the time segment is recorded as an attitude loss of control anomaly event, and the corresponding timestamp and spatial region are output; this method completely eliminates the interference of local positioning "flying points" caused by non-line-of-sight occlusion on the final anomaly judgment, because these "flying points" are filtered out at the bottom layer by the kernel function and kinematic chain constraints in step S4; Based on the coupling characteristics of projected displacement and spatial displacement residuals, motion pattern anomalies are identified, including: Calculate the phase diagram closure, the phase diagram is determined by tangential displacement. and normal displacement The calculation method for phase diagram closure is as follows: extract the trajectory line segment of a complete action cycle in the phase diagram, calculate the Euclidean distance between the first and last endpoints of the line segment, and use this distance as the closure index. A preset envelope is established by collecting data from several normal motion patterns. - Phase diagram trajectory, calculated by dividing the phase plane into angular bins and statistically averaging the radial distances. and standard deviation ,by It constitutes the normal motor envelope; The non-stationarity test is performed on the residual sequence of the motion chain consistency factor. The statistical test method for non-stationary divergence is as follows: the augmented Dickey-Fowler test is performed on the residual sequence. If the test statistic is greater than the critical value, the sequence is determined to be non-stationary. At the same time, the sliding window variance of the residual sequence is calculated. When the sliding variance increases continuously and exceeds a predetermined multiple of the initial stage variance, the residual is determined to be non-stationary divergence. An abnormal event is defined as when the phase diagram closure exceeds the preset closure threshold, the phase diagram trajectory exceeds the normal motion envelope, and the residual is simultaneously determined to be non-stationary divergence.

[0045] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for detecting motion trajectory anomalies by fusing multi-source positioning data, characterized in that, The detection method includes the following steps: The system collects raw ranging values ​​from the ultra-wideband positioning system, acceleration and angular velocity data from the inertial measurement unit, and plantar pressure sensor data, and synchronizes them in time. The plantar pressure sensor data is used to identify the zero-velocity state upon ground contact and to extract the zero-velocity virtual observation constraint. Based on the lower limb bone length calibration value and the angular velocity data of the inertial measurement unit, a limb kinematic chain model constrained by the joint angle range is constructed, and the spatial displacement of the waist relative to the supporting foot is derived. The original ranging value, the pre-integrated quantity of the inertial measurement unit, the zero-velocity virtual observation constraint, and the spatial displacement are incorporated into the factor graph for nonlinear optimization to solve the waist positioning trajectory; The waist positioning trajectory is projected onto a local coordinate system formed by the propulsion direction of the moving vehicle and the normal vector to obtain tangential displacement and normal displacement. Based on the phase diagram closure degree of the tangential displacement and normal displacement, as well as the residual statistical characteristics between the spatial displacement and the displacement increment obtained by factor diagram optimization, the motion mode is determined to be abnormal.

2. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 1, characterized in that: The identification of the zero-speed state upon ground contact includes: processing the plantar pressure sensing data using a continuous sampling window; when the proportion of ground contact marker frames within the window to the total number of frames is greater than 0.9, and the pressure value of the first frame of the window increases by more than 30% compared to the pressure value of the previous frame, the corresponding foot is determined to be in a zero-speed state, and a speed observation residual is established. In the formula, For a moment The zero-velocity virtual observation residual vector; For a moment feet The instantaneous velocity vector to be estimated; It is a three-dimensional zero vector; For physical time variables.

3. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 1, characterized in that: The lower limb bone length calibration value is obtained as follows: With the subject maintaining a zero-speed, single-leg upright support position, the pitch angle of the supporting leg is acquired by the inertial measurement unit. Combined with pre-inputted height parameters, initial estimates of the thigh and lower leg lengths are calculated using a human leg-to-body ratio coefficient of 0.45 to 0.

53. The thigh length calibration value is then obtained by solving the geometric relationship of the limb projection in the support position. and lower leg length calibration value ; The joint angle range constraints are as follows: the range of motion of the knee joint flexion and extension angle is limited to 0° to 130°, and the range of motion of the hip joint flexion and extension angle is limited to -15° to 120°. The range of motion is a hard boundary preset based on the physiological limit angle of the human lower limb skeletal structure in the sagittal plane without dislocation or fracture. When deriving the spatial displacement of the waist relative to the supporting foot, a penalty is imposed on the state variable that exceeds the range of motion.

4. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 3, characterized in that: The penalty is implemented by constructing a cost function, which is specifically in the form of: in the factor graph optimization objective function, the knee joint angle is... Introducing penalty items Regarding the hip joint angle Introducing penalty items ,in This is a function of the square of the one-sided excess, taking a value of zero when the angle is within the activity range and a value of the square of the excess when it exceeds the boundary; penalty coefficient. and The method for determining the penalty coefficient is as follows: based on the angular velocity noise variance of the inertial measurement unit and the passive resistance torque characteristics of human joints at extreme angles, the penalty coefficient is set so that the increment of the penalty term generated when the angle exceeds the constraint boundary by 1° is on the same order of magnitude as the typical variance value of the zero-velocity observation residual, thereby ensuring the reasonable weight of physiological constraints in optimization.

5. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 1, characterized in that: The factor graph includes: UWB ranging factor, whose residual form is... In the formula, For the reason about the Ranging residual scalar of a UWB base station; For waist label at all times The position vector, For the first The location vectors of each base station. The UWB ranging value in step S1; And the kinematic chain consistency factor, whose residual is the difference between the spatial displacement and the position increment of the waist state node at the current moment relative to the waist state node at the last zero-velocity update moment.

6. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 5, characterized in that: The UWB ranging factor incorporates a kernel function based on kinematic conflict for weight reduction processing. The mathematical expression of the kernel function is as follows: In the formula, The residual vector of the UWB ranging factor; The covariance matrix of the UWB ranging factor; The Mahalanobis distance; Threshold; threshold The method for determining the threshold is as follows: the Mahalanobis distance is regarded as a chi-square distribution variable, and its degree of freedom is equal to the upper 95th percentile of the number of effective base stations of the UWB ranging factor at the current time. In this way, abnormal ranging values ​​that exceed the sensor noise level are statistically judged as kinematic conflicts and their weights are reduced in the form of an inverse proportional function.

7. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 1, characterized in that: The process of solving the waist positioning trajectory includes: performing nonlinear least squares optimization on the factor map, and iteratively estimating the position, velocity, attitude of the waist tag, and the zero bias parameters of the inertial measurement unit.

8. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 1, characterized in that: The local coordinate system uses the propulsion direction of the moving vehicle and its normal vector as coordinate axes; the waist positioning trajectory is projected onto the local coordinate system to obtain the time series of tangential displacement and normal displacement.

9. The method for detecting motion trajectory anomalies by fusing multi-source positioning data according to claim 1, characterized in that: The determination of motion pattern anomalies based on the coupling characteristics of projected displacement and spatial displacement residual includes: Calculate the phase diagram closure, the phase diagram being determined by the tangential displacement. and normal displacement The phase diagram closure is calculated as follows: extract the trajectory line segment of a complete action cycle in the phase diagram, calculate the Euclidean distance between the first and last endpoints of the line segment, and use this distance as the closure index. The preset envelope is established by collecting data from several normal motion patterns. - Phase diagram trajectory, calculated by dividing the phase plane into angular bins and statistically averaging the radial distances. and standard deviation ,by It constitutes the normal motor envelope; The residual sequence of the motion chain consistency factor is subjected to a non-stationarity test. The statistical test method for non-stationary divergence is as follows: the augmented Dickey-Fowler test is performed on the residual sequence. If the test statistic is greater than the critical value, the sequence is determined to be non-stationary. At the same time, the sliding window variance of the residual sequence is calculated. When the sliding variance increases continuously and exceeds a predetermined multiple of the variance in the initial stage, the residual is determined to be non-stationary diverging. An abnormal event is determined when the phase diagram closure exceeds a preset closure threshold, the phase diagram trajectory exceeds the normal motion envelope, and the residual is simultaneously determined to be non-stationary divergence.

Images