Analysis of Kinematic Characteristics of Saanen Goat Spine under Multi-Slope

In order to improve the slope movement stability and flexibility of quadruped robot, a theoretical design method of a flexible spine of a robot that was based on bionics was proposed. The kinematic characteristics of the spine were analyzed under different slopes with a Saanen goat as the research object. A Qualisys track manager (QTM) gait analysis system was used to obtain the trunk movement of goats under multiple slopes, and linear time normalization (LTN) was used to calibrate and match typical gait cycles to characterize the goat locomotion gait under slopes. Firstly, the spatial angle changes of cervical thoracic vertebrae, thoracolumbar vertebrae, and lumbar vertebrae were compared and analyzed under 0°, 5°, 10°, and 15° slopes, and it was found that the rigid and flexible coupling structure between the thoraco–lumbar vertebrae played an obvious role when moving on the slope. Moreover, with the increase in slope, the movement of the spine changed to the coupling movement of thoraco–lumbar coordination movement and a flexible swing of lumbar vertebrae. Then, the Gaussian mixture model (GMM) clustering algorithm was used to analyze the changes of the thoraco–lumbar vertebrae and lumbar vertebrae in different directions. Combined with anatomical knowledge, it was found that the motion of the thoraco–lumbar vertebrae and lumbar vertebrae in the goat was mainly manifested as a left–right swing in the coronal plane. Finally, on the basis of the analysis of the maximin and variation range of the thoraco–lumbar vertebrae and lumbar vertebrae in the coronal plane, it was found that the coupling motion of the thoraco–lumbar cooperative motion and flexible swing of the lumbar vertebrae at the slope of 10° had the most significant effect on the motion stability. SSE, R2, adjusted-R2, and RMSE were used as evaluation indexes, and the general equations of the spatial fitting curve of the goat spine were obtained by curve fitting of Matlab software. Finally, Origin software was used to obtain the optimal fitting spatial equations under eight movements of the goat spine with SSE and adjusted-R2 as indexes. The research will provide an idea for the bionic spine design with variable stiffness and multi-direction flexible bending, as well as a theoretical reference for the torso design of a bionic quadruped robot.


Introduction
At present, most quadruped robots adopt rigid body structure, which is difficult to meet the needs of flexibility and bending when moving [1,2]. The flexible spine structure of tetrapods provides auxiliary power when they move so that they can show a variety of dexterous movement postures, demonstrating strong flexibility, stability, and environmental adaptability [3,4]. Hildebrand [5] studied the running process of cheetahs and horses, found that the torso of cheetahs was more curled and extended than that of horses, as well as the fact that the spine movement mainly acted on the legs, which could effectively increase the leg extension. They discussed the influence of leg movement on the running speed under the action of the spine. Hackert [6] observed the spinal movement of the pika. It was found that the spinal movement kept the body center of mass in a smooth motion trajectory that could effectively maintain the body balance and increase motion stability of the pika. Schilling et al. [7] used high-speed cineradiography to study the spinal intervertebral joint movement of mammals during running and semi-jumping, and it was found that spinal movement had an important impact on animal motor performance.
Therefore, the design of multidirectional flexible bending biomimetic spine structure, which is based on the motion characteristics of a tetrapod spine, is an important research direction in order to improve the kinematic capability of quadruped robots in an unstructured environment [8,9]. Zhang et al. [10] studied the spine structure and trunk movement characteristics of cheetahs and designed the passive compliant spine of a quadruped robot with adjustable stiffness on the basis of the bionics principle. It was found that when the spine stiffness was 681 N/m, the stable motion speed of the designed robot was the highest, which was 0.92 m/s. Zhao [11] designed and optimized the biomimetic trunk mechanism on the basis of the study of kinematic characteristics of goat torso and simulated it. It was found that the combination of goat kinematic characteristic parameters and robot motion control could effectively improve the robot motion effect. Wei [12] built a quadruped dynamic model to describe Transverse Gallop gait and Rotary Gallop gait. Through a large number of simulation experiments, it was found that the Transverse Gallop gait was designed for efficient, long-distance running, while the Rotary Gallop gait was designed for high-speed running in small-to-medium quadruped mobility systems. On the basis of Rotary Gallop gait, SQbot, a quadruped robot with bionic spinal joints, was developed. The overall trunk motion amplitude of the quadruped robot at higher moving speed was significantly reduced compared with that at lower moving speed. Li [13] studied the hopping gait of a quadruped robot with spinal joints and analyzed its gait parameters. Finally, a simplified model of a quadruped robot with spinal joints was established according to the bionics principle and the skeleton morphology of a cheetah. Physical validation found that the impedance coefficient of the robot with the best spine stiffness was 40% lower than that of the robot without the flexible spine. Borisova et al. [14] studied a method to match the dynamic model parameters of a cheetah with a running robot with a flexible spine through video tracking analysis of cheetah movement, and it was found that a flexible spine is the main feature of an efficient running robot.
Compared with other tetrapods, goats have smaller bodies, solid and powerful limb bones, and agile and light gait, and as well as this, they can adapt to different rugged terrain and road surface changes, showing strong movement stability [15][16][17]. In order to explore the structural characteristics of the goat spine, Van Engelen et al. [18] used a structural vibration test method to study the mechanical properties of goat spine motion segments. It was found that there were obvious formants in flexion and extension, lateral bending, and axial rotation of all goat lumbar spine motion segments, which verified the feasibility of the vibration test to evaluate the modal parameters of lumbar spine motion segments. Zhao [11] studied the trajectory changes of the trunk and the sensitivity of different parts by applying different load pressures on the trunk bones of goats. The difference in living environment between goats and sheep leads to the difference of their locomotor characteristics, but their spinal structures are the same, and thus the study of sheep static structural characteristics is of significance. Valentin et al. [19] found that the range of movement of the spine of sheep was smaller than that of humans, no matter what posture they adopted in comparison with the range of movement of the spine of people of different ages and sheep. DeVries et al. [20] tested the flexion and extension, lateral bending, and axial rotation of the cervical vertebrae of adult sheep (C2-C7) by a test device. They found the main behaviors and coupling behaviors of the C2-C7 segments of the cervical vertebrae of sheep. McDonald et al. [21] implanted steel beads into the vertebrae of sheep specimens, manually stretching neck extension and axial rotation. The three-dimensional position and direction of each cervical spine were determined using biplanar X-ray images. They verified the accuracy of modelbased tracking measurement of three-dimensional dynamic cervical spine kinematics. At present, studies on the biological characteristics of goat spine are mainly focused on the analysis of the structural characteristics and mechanical properties of the dissected spine using specific methods, as well as the motion characteristics of goat spine under different states using a high-speed camera. However, there are few studies on spatial motion characteristics of goat spine through intra-anatomical spine motion under multiple slopes. Moreover, few scholars have established spinal spatial curves through spinal motion characteristics.
To explore the dynamic motion characteristics of goat spine and construct a mathematical model of a spine curve, the QTM gait analysis system was used to obtain the trunk movement of goats under multiple slopes. Combined with typical gait cycle, the spatial dynamic motion characteristics of goat spine were analyzed. The Gaussian mixture clustering model was used to analyze the influence of the motion characteristics of the spine in the anatomical plane on its spatial motion characteristics. Moreover, the slope with a significant effect was analyzed according to the data variation range of the spine in the anatomical plane under multiple slopes. The research will provide an idea for the bionic spine design with variable stiffness levels and multi-direction flexible bending, as well as providing a theoretical reference for the torso design of a bionic quadruped robot.

Materials and Methods
The Saanen goat was taken as the research object. A young male goat with moderate body size and good body condition was selected for the study. Before the experiment, the weight of the goat was measured at 19.50 kg by a scale (accuracy 0.01 kg), and the goat was scanned by CT tomography. During the scanning, the voltage was set to 140 kV, the current was 288 mA, and the layer thickness was 0.6 mm. The goat was anesthetized and asleep. The Dicom format file of the goat trunk was generated by CT scanning, and the three-dimensional structure of spine was reconstructed by Mimics software (Version 19.0, Materialize, Leuven, Belgium).
The spine of the selected subject included 7 cervical vertebrae (purple part), 13 thoracic vertebrae (green part), and 6 lumbar vertebrae (blue part) without obvious lesions. The seventh cervical vertebra showed obvious movement, and therefore the cervical markers were C at the first cervical vertebra and T at the seventh cervical vertebra. The thoracic vertebrae are connected to the ribs, which to some extent limit the excessive bending of the thoracic vertebrae, and the stability of the thoracic vertebrae is significantly higher than that of other spinal regions. Therefore, the markers for thoracic vertebrae were L1 at the 13th thoracic vertebrae, and L2 at the third lumbar vertebrae and L3 at the sixth lumbar vertebrae were selected as markers for lumbar vertebrae. Hair on the spine of the test subject was removed before labeling. Reflective balls were pasted at T (the most obvious position at the junction between cervical vertebrae and thoracic vertebrae), L1 (the thoracic vertebrae connected with the last rib of the goat), L2 (the lumbar vertebrae connected with the hip bone of the goat), and L2 (the middle position of L1 and L3). In addition, reflective balls were pasted on the limbs to obtain the locomotion gait of the goat, and LFn (left front), LHn (left hind), RFn (right front), and RHn (right hind) were labeled. The attachment position of the reflective balls is shown in Figure 1. In this paper, the change of angle CTL1 represented the motion of cervical vertebrae relative to thoracic vertebrae of goats, the change of angle TL1L2 represented the motion of thoracic vertebrae relative to lumbar vertebrae of goats, and the change of angle L1L2L3 represented the internal motion of lumbar vertebrae of goats. They were labeled as θ 1 , θ 2 , and θ 3 , respectively, for subsequent analysis. vertebrae of goats, the change of angle TL1L2 represented the motion of thoracic vertebrae relative to lumbar vertebrae of goats, and the change of angle L1L2L3 represented the internal motion of lumbar vertebrae of goats. They were labeled as θ1, θ2, and θ3, respectively, for subsequent analysis.  A test platform consisting of a 6-DOF motion platform, a motion capture system (Qualisys Miqus M3 model, displacement error ≤0.1 mm, frame rate of 340 fps, system delay ≤5 ms), and a data acquisition and analysis system was used for the test. This is shown in Figure 2. Before the experiment, the platform was calibrated, and the basic space coordinate system was established. The positive direction of the X-axis was the movement direction of the goat, the positive direction of the Y-axis was the horizontal direction of the left side of the advance, and the positive direction of the Z-axis was the vertical geodetic direction. The sampling frequency of the motion capture system was set to 100 Hz, and the frame rate was set to 100 fps. The angle between the platform and the Z-axis is α. The slope test was carried out by changing the angle value of α, and the inclination angles α were 0°, 5°, 10°, and 15°. During the experiment, the goat's horns were tied with a rope and pulled by the experimenter at the back end. The rope was relaxed to ensure that the goat head was in a state of free movement, and the goat was lured forward by food induction. A test platform consisting of a 6-DOF motion platform, a motion capture system (Qualisys Miqus M3 model, displacement error ≤ 0.1 mm, frame rate of 340 fps, system delay ≤ 5 ms), and a data acquisition and analysis system was used for the test. This is shown in Figure 2. Before the experiment, the platform was calibrated, and the basic space coordinate system was established. The positive direction of the X-axis was the movement direction of the goat, the positive direction of the Y-axis was the horizontal direction of the left side of the advance, and the positive direction of the Z-axis was the vertical geodetic direction. The sampling frequency of the motion capture system was set to 100 Hz, and the frame rate was set to 100 fps. The angle between the platform and the Z-axis is α. The slope test was carried out by changing the angle value of α, and the inclination angles α were 0 • , 5 • , 10 • , and 15 • . During the experiment, the goat's horns were tied with a rope and pulled by the experimenter at the back end. The rope was relaxed to ensure that the goat head was in a state of free movement, and the goat was lured forward by food induction. Biomimetics 2022, 7, x FOR PEER REVIEW 5 of 21 Figure 2. The motion test platform consists of ① the motion capture system, ② 6-DoF motion platform, and ③ the three-dimensional biomechanical force measuring table.
It was found that goat mainly adopted a "three-phase supported, two-phase supported, three-phase supported" gait when moving on the slope. That is, when starting, one front leg swings and the other three legs support. Before the swing phase falls to the ground and becomes the supported phase, the other leg is raised and swings. In this study, the leg lifting sequence of "left front-right hind-right front-left hind" was analyzed as a typical gait cycle, and the locomotion gait is shown in Figure 3, which consisted of eight movements (①-⑨, with the ⑨ movement being the same as the ① movement). Goat locomotion gait on the slopes consisting of eight movements, which are ① the left front leg swung monophasically; ② the left front leg and the right hind leg swung biphasically; ③ the right hind leg swung monophasically; ④ the right hind leg and the right front leg swung biphasically; ⑤ the right front leg swung monophasically; ⑥ the right front leg and the left hind leg swung biphasically; ⑦ the left hind leg swung monophasically; ⑧ the left hind leg and the left front leg swung biphasically.
In order to accurately screen the gait cycle, the movement track of the limb markers was observed. The movement track of the third limb marker along the positive direction of X-, Y-, and Z-axes under the slope of 0° is shown in Figure 4. It was found that goat mainly adopted a "three-phase supported, two-phase supported, three-phase supported" gait when moving on the slope. That is, when starting, one front leg swings and the other three legs support. Before the swing phase falls to the ground and becomes the supported phase, the other leg is raised and swings. In this study, the leg lifting sequence of "left front-right hind-right front-left hind" was analyzed as a typical gait cycle, and the locomotion gait is shown in Figure 3, which consisted of eight movements ( 1 -9 , with the 9 movement being the same as the 1 movement). It was found that goat mainly adopted a "three-phase supported, two-phase supported, three-phase supported" gait when moving on the slope. That is, when starting, one front leg swings and the other three legs support. Before the swing phase falls to the ground and becomes the supported phase, the other leg is raised and swings. In this study, the leg lifting sequence of "left front-right hind-right front-left hind" was analyzed as a typical gait cycle, and the locomotion gait is shown in Figure 3, which consisted of eight movements (①-⑨, with the ⑨ movement being the same as the ① movement). Goat locomotion gait on the slopes consisting of eight movements, which are ① the left front leg swung monophasically; ② the left front leg and the right hind leg swung biphasically; ③ the right hind leg swung monophasically; ④ the right hind leg and the right front leg swung biphasically; ⑤ the right front leg swung monophasically; ⑥ the right front leg and the left hind leg swung biphasically; ⑦ the left hind leg swung monophasically; ⑧ the left hind leg and the left front leg swung biphasically.
In order to accurately screen the gait cycle, the movement track of the limb markers was observed. The movement track of the third limb marker along the positive direction of X-, Y-, and Z-axes under the slope of 0° is shown in Figure 4.  It is easy to see that the X and Z direction trajectories of the marker points had obvious periodicity, so when the leg was considered to be in the supported posture and it was the supported phase.
In order to study the kinematic characteristics of the goat spine at different slopes, the key roles of each part of the spine were explored through the spatial curvature of the spine and the maximin and variation range of θ1, θ2, and θ3. Then Gaussian mixture model (GMM) clustering method was used to determine the main direction of spinal curvature. Moreover, the main acting slope of the goat spine was analyzed, combined with the maximin and variation range of θ2 and θ3 in the main bending direction. Finally, the spine curve of eight movements under the main action slope was fitted. It is easy to see that the X and Z direction trajectories of the marker points had obvious periodicity, so when the leg was considered to be in the supported posture and it was the supported phase.
In order to study the kinematic characteristics of the goat spine at different slopes, the key roles of each part of the spine were explored through the spatial curvature of the spine and the maximin and variation range of θ 1 , θ 2 , and θ 3 . Then Gaussian mixture model (GMM) clustering method was used to determine the main direction of spinal curvature. Moreover, the main acting slope of the goat spine was analyzed, combined with the maximin and variation range of θ 2 and θ 3 in the main bending direction. Finally, the spine curve of eight movements under the main action slope was fitted.

Results
In order to describe the motion affected by the speed of goat movement, three groups of data were selected for analysis at different slopes to illustrate the general motion characteristics of the goat spine. For the same kind of movement, the goat movement speed and the step length were different, resulting in different gait cycle time. Linear time normalization (LTN) was used for sequence calibration and matching of gait cycle. The normalization mapped the data to the interval [0, 1] without changing the data features, so as to facilitate the comparison and analysis of spine movement during gait cycle. The normalization formula is as follows: where Frame min is the minimum frame in the gait cycle and Frame max is the maximum frame in the gait cycle. Figures 5-8 were drawn with the normalized gait cycle from 1 movement to 8 movement as the horizontal axis and the change of θ 1 , θ 2 , and θ 3 as the vertical axis.

Analysis of the Motion of Each Part of Spine
Combined with the eight movements under the 0 • slope as shown in Figure 5d (see Figure 3 for the explanation of movements), the general changes of θ 1 , θ 2 , and θ 3 as shown in Figure 5a-c were analyzed. It was found that the overall change trend of all angles was the same. When the front limb swung, the angle θ 1 decreased, and when the hind limb swung, the angle θ 1 increased. When the front limb was supported, the angle θ 1 reached the peak, and when the hind limb was supported, the angle θ 1 reached the minimum. When the right front leg swung, the angle of θ 2 reached the peak value of 170 • . When the motion state of the goat ipsilateral leg switched, the angle change rate of θ 2 decreased. The angle change rate of θ 2 increased when the motion state of the corresponding leg switched. When the left front leg swung monophasically, the angle of θ 3 increased slowly. When the left front leg and right hind leg swung, the angle change rate of θ 3 increased. When the right front leg swung monophasically, θ 3 reached the peak and then decreased, and the peak value was about 170 • . When the left hind leg swung monophasically, the angle θ 3 decreased.
Combined with the eight movements under the 5 • slope as shown in Figure 6d (see Figure 3 for the explanation of movements), the general changes of θ 1 , θ 2 , and θ 3 as shown in Figure 6a-c were analyzed. It was found that when the left front leg swung, the angle of θ 1 changed about 5 • , the angle of θ 2 decreased, and the angle of θ 3 decreased. When the left front leg and right hind leg swung, the angle of θ 1 increased, the angle of θ 2 decreased, and the angle of θ 3 changed about 5 • . When the right hind leg swung monophasically, the angle of θ 1 changed about 3 • , the angle of θ 2 changed about 3 • , and the angle of θ 3 changed about 4 • . When the right front leg and right hind leg swung on the same side, the angle of θ 1 decreased, the angle of θ 2 changed about 3 • , and the angle of θ 3 changed about 5 • . When the right front leg swung monophasically, the angle of θ 1 changed about 5 • , the angle of θ 2 increased, and the angle of θ 3 changed about 5 • . When the right front leg and left hind leg swung on the different side, the angle of θ 1 and θ 2 increased, while the angle of θ 3 decreased. When the left hind leg swung monophasically, the angle of θ 1 , θ 2 , and θ 3 all decreased.
Combined with the eight movements under 10 • slope as shown in Figure 7d (see Figure 3 for the explanation of movements), the general changes of θ 1 , θ 2 , and θ 3 as shown in Figure 7a-c were analyzed. It was found that the angle of θ 1 and θ 3 decreased, and the angle of θ 2 changed about 10 • when the left front leg swung. When the left front leg and right hind leg swung, the angle of θ 1 and θ 3 changed about 5 • , and the angle of θ 3 increased. When the right hind leg swung, the angle of θ 1 changed about 2 • , the angle of θ 2 changed about 5 • , and the angle of θ 3 increased. When the right front and right hind legs swung on the same side, the angle of θ 1 changed about 2 • , the angle of θ 2 decreased, and the angle of θ 3 changed about 10 • . When the right front leg swung monophasically, the angle of θ 1 changed about 2 • , and the angle of θ 2 and θ 3 increased. When the right front leg and left hind leg swung on the different side, the angle of θ 1 decreased, the angle of θ 2 increased, and the angle of θ 3 increased after decreasing. When the left hind leg swung monophasically, the angle of θ 1 decreased, the angle of θ 2 basically stayed the same, and the angle of θ 3 decreased by about 5 • . Combined with the eight movements under the 5° slope as shown in Figure 6d (see Figure 3 for the explanation of movements), the general changes of θ1, θ2, and θ3 as shown in Figure 6a-c were analyzed. It was found that when the left front leg swung, the angle of θ1 changed about 5°, the angle of θ2 decreased, and the angle of θ3 decreased. When the left front leg and right hind leg swung, the angle of θ1 increased, the angle of θ2 decreased, and the angle of θ3 changed about 5°. When the right hind leg swung monophasically, the angle of θ1 changed about 3°, the angle of θ2 changed about 3°, and the angle of θ3 changed about 4°. When the right front leg and right hind leg swung on the same side, the angle of θ1 decreased, the angle of θ2 changed about 3°, and the angle of θ3 changed about 5°. When the right front leg swung monophasically, the angle of θ1 changed about 5°, the angle of θ2 increased, and the angle of θ3 changed about 5°. When the right front leg and left hind leg swung on the different side, the angle of θ1 and θ2 increased, while the angle of θ3 decreased. When the left hind leg swung monophasically, the angle of θ1, θ2, and θ3 all decreased. Combined with the eight movements under the 15 • slope as shown in Figure 8d (see Figure 3 for the explanation of movements), the general changes of θ 1 , θ 2 , and θ 3 as shown in Figure 8a-c were analyzed. It was found that the angle of θ 1 decreased, the angle of θ 2 increased, and the angle of θ 3 changed to about 155 • when the left front leg swung. When the left front leg and right hind leg swung on the different side, the angle of θ 1 increased, the angle of θ 2 changed about 3 • , and the angle of θ 3 decreased. When the right hind leg swung, the angle of θ 1 changed about 5 • , the angle of θ 2 decreased, and the angle of θ 3 rose about 1 • . When the right front leg and right hind leg swung on the same side and the right front leg swung monophasically, the angle of θ 1 decreased, and the angle of θ 2 and θ 3 increased. When the right front leg and left hind leg swung on the different side, the angle of θ 1 increased, the angle of θ 2 changed about 4 • , and the angle of θ 3 changed about 5 • . When the left hind leg swung monophasically, the angle of θ 1 decreased, the angle of θ 2 increased, and the angle of θ 3 increased after it decreased. Biomimetics 2022, 7, x FOR PEER REVIEW 9 of 21 Combined with the eight movements under 10° slope as shown in Figure 7d (see Figure 3 for the explanation of movements), the general changes of θ1, θ2, and θ3 as shown in Figure 7a-c were analyzed. It was found that the angle of θ1 and θ3 decreased, and the angle of θ2 changed about 10° when the left front leg swung. When the left front leg and right hind leg swung, the angle of θ1 and θ3 changed about 5°, and the angle of θ3 increased. When the right hind leg swung, the angle of θ1 changed about 2°, the angle of θ2 changed about 5°, and the angle of θ3 increased. When the right front and right hind legs swung on the same side, the angle of θ1 changed about 2°, the angle of θ2 decreased, and the angle of θ3 changed about 10°. When the right front leg swung monophasically, the angle of θ1 changed about 2°, and the angle of θ2 and θ3 increased. When the right front leg and left hind leg swung on the different side, the angle of θ1 decreased, the angle of θ2 increased, and the angle of θ3 increased after decreasing. When the left hind leg swung Through the above analysis, it was found that the leg movement of goats was closely related to the spine movement, and the spine movement showed different phenomena at different slopes. Therefore, the spine played a significant role in the movement process of goats under slopes.

Analysis of the Maximin and Variation
Range of θ 1 , θ 2 , and θ 3 The analysis of the maximin and range of the change of each key angle of the spine under different slopes is an important method in terms of studying the action of each part of the spine. The motion of the cervical thoracic vertebrae, thoraco-lumbar vertebrae, and lumbar vertebrae were analyzed by the maximin and variation range of the key spinal angles (θ 1 , θ 2 , and θ 3 ). The maximin and variation range of them are shown in Figure 9. monophasically, the angle of θ1 decreased, the angle of θ2 basically stayed the same, and the angle of θ3 decreased by about 5°. Combined with the eight movements under the 15° slope as shown in Figure 8d (see Figure 3 for the explanation of movements), the general changes of θ1, θ2, and θ3 as shown in Figure 8a-c were analyzed. It was found that the angle of θ1 decreased, the angle of θ2 increased, and the angle of θ3 changed to about 155° when the left front leg swung. When the left front leg and right hind leg swung on the different side, the angle of θ1 increased, the angle of θ2 changed about 3°, and the angle of θ3 decreased. When the right hind leg swung, the angle of θ1 changed about 5°, the angle of θ2 decreased, and the angle of θ3 rose about 1°. When the right front leg and right hind leg swung on the same side and the right front leg swung monophasically, the angle of θ1 decreased, and the angle of θ2 and θ3 increased. When the right front leg and left hind leg swung on the different side, the angle of θ1 increased, the angle of θ2 changed about 4°, and the angle of θ3 changed about 5°. When the left hind leg swung monophasically, the angle of θ1 decreased, the angle of θ2 increased, and the angle of θ3 increased after it decreased. According to Figure 9a, the angle of θ 1 decreased first and then increased with the increase in the slope. At the same slope, θ 2 was larger than θ 3 , indicating that the role of thoracolumbar coordination was more significant in the process of goat movement. Since the thoracic vertebrae has better stability than the lumbar vertebrae to some extent, the motion of the thoracolumbar vertebrae could be regarded as the motion of rigid and flexible coupling structure. Through the analysis shown Figure 9b, it was found that the variation range of θ 2 was larger than that of θ 3 under the slope of 0 • . Conversely, the variation range of θ 2 was smaller than that of θ 3 under the slopes of 5 • , 10 • , and 15 • , indicating that with the increase in slope, the regulation function of the lumbar spine was increasingly more important. The main movement of the spine changed to the coupling movement between the thoracolumbar coordination movement and the flexible swing of the lumbar vertebrae. The range of θ 1 gradually increased with the increase in slope, indicating that the head had an important regulating effect on the movement stability of the goat under large slopes. This was consistent with the team's previous research results [16]. Through the above analysis, it was found that the leg movement of goats was closely related to the spine movement, and the spine movement showed different phenomena at different slopes. Therefore, the spine played a significant role in the movement process of goats under slopes.

Analysis of the Maximin and Variation Range of θ1, θ2, and θ3
The analysis of the maximin and range of the change of each key angle of the spine under different slopes is an important method in terms of studying the action of each part of the spine. The motion of the cervical thoracic vertebrae, thoraco-lumbar vertebrae, and lumbar vertebrae were analyzed by the maximin and variation range of the key spinal angles (θ1, θ2, and θ3). The maximin and variation range of them are shown in Figure 9. According to Figure 9a, the angle of θ1 decreased first and then increased with the increase in the slope. At the same slope, θ2 was larger than θ3, indicating that the role of thoracolumbar coordination was more significant in the process of goat movement. Since the thoracic vertebrae has better stability than the lumbar vertebrae to some extent, the motion of the thoracolumbar vertebrae could be regarded as the motion of rigid and flexible coupling structure. Through the analysis shown Figure 9b, it was found that the var-

Analysis of Motion of Thoracic and Lumbar Vertebrae under Different Slopes
The goat advanced along the positive direction of the X-axis. In order to further describe the motion of thoracolumbar vertebrae and lumbar vertebrae, the motion of θ 2 and θ 3 angles under different slopes was studied in combination with anatomical orientation terms. In the motion coordinate system, the motion of the spine was mainly represented by the motion (up and down fluctuation) in the XZ plane, which is the sagittal plane, and the motion (left and right swing) in the XY plane, which is the coronal plane. The Gaussian mixture model (GMM) is a common data clustering method that combines similar data together according to the similarity [22,23]. On the basis of the GMM method and expected maximum (EM) algorithm, the effects of covariance matrix types and initial conditions on the clustering results of Gaussian mixture models was studied. The optimal combination of two covariance matrices of diagonal and full and the conditions of identical and non-identical covariance matrices was explored.
Firstly, the clustering model between the data of θ 2 in the XZ and XY planes and the spatial data was established to analyze the influence of the motion of the thoracic-lumbar vertebra coupling structure in the two planes. The clustering results of θ 2 data in XZ and XY planes are shown in Figure 10. (g) (h) Figure 10. Clustering results of θ2 data in the XZ and XY planes: (a) angle data of θ2 in the XZ plane, when Sigmal was diagonal, and shared covariance was true; (b) angle data of θ2 in the XZ plane, when Sigmal was diagonal, and shared covariance was false; (c) angle data of θ2 in the XZ plane, when Sigmal was full, and shared covariance was true; (d) angle data of θ2 in the XZ plane, when Sigmal was full, and shared covariance was false; (e) angle data of θ2 in the XY plane, when Sigmal was diagonal, and shared covariance was true; (f) angle data of θ2 in the XY plane, when Sigmal was diagonal, and shared covariance was false; (g) angle data of θ2 in the XY plane, when Sigmal was full, and shared covariance was true; (h) angle data of θ2 in the XY plane, when Sigmal was full, and shared covariance was false.
By comparing (a) with (e), (b) with (f), and (c) and (g) in Figure 10, it was found that under the same covariance matrix and initial conditions, the data clustering effect of θ2 in the XY plane was obvious and could be clearly divided into four independent and complete regions. That was, the data could be accurately divided into the data under four slopes. According to this, the motion of θ2 in the XY plane had the most obvious corresponding relationship with the slope. In this case, the motion of the thoracolumbar spine under different slopes was mainly represented by the left and right swing in the coronal plane.
As shown in Figure 10e-g, it was found that when the covariance matrix was diagonal and different covariance matrices were selected, the four slope regions were most clearly delimited, and there was a certain distance between the regions. This indicated that the GMM had the best clustering effect under this condition, and the angle data under different slopes accurately corresponded to the corresponding slopes, which could be used to analyze the main direction of spinal curvature. Then, under the condition, clustering models between the angle data of θ3 in XZ and XY planes and spatial data were established in order to analyze the influence of lumbar curvature in the sagittal plane and in the coronal plane on spinal motion. The clustering results of θ3 in the XZ and XY planes are shown in Figure 11. Figure 10. Clustering results of θ 2 data in the XZ and XY planes: (a) angle data of θ 2 in the XZ plane, when Sigmal was diagonal, and shared covariance was true; (b) angle data of θ 2 in the XZ plane, when Sigmal was diagonal, and shared covariance was false; (c) angle data of θ 2 in the XZ plane, when Sigmal was full, and shared covariance was true; (d) angle data of θ 2 in the XZ plane, when Sigmal was full, and shared covariance was false; (e) angle data of θ 2 in the XY plane, when Sigmal was diagonal, and shared covariance was true; (f) angle data of θ 2 in the XY plane, when Sigmal was diagonal, and shared covariance was false; (g) angle data of θ 2 in the XY plane, when Sigmal was full, and shared covariance was true; (h) angle data of θ 2 in the XY plane, when Sigmal was full, and shared covariance was false.
By comparing (a) with (e), (b) with (f), and (c) and (g) in Figure 10, it was found that under the same covariance matrix and initial conditions, the data clustering effect of θ 2 in the XY plane was obvious and could be clearly divided into four independent and complete regions. That was, the data could be accurately divided into the data under four slopes. According to this, the motion of θ 2 in the XY plane had the most obvious corresponding relationship with the slope. In this case, the motion of the thoracolumbar spine under different slopes was mainly represented by the left and right swing in the coronal plane.
As shown in Figure 10e-g, it was found that when the covariance matrix was diagonal and different covariance matrices were selected, the four slope regions were most clearly delimited, and there was a certain distance between the regions. This indicated that the GMM had the best clustering effect under this condition, and the angle data under different slopes accurately corresponded to the corresponding slopes, which could be used to analyze the main direction of spinal curvature. Then, under the condition, clustering models between the angle data of θ 3 in XZ and XY planes and spatial data were established in order to analyze the influence of lumbar curvature in the sagittal plane and in the coronal plane on spinal motion. The clustering results of θ 3 in the XZ and XY planes are shown in Figure 11.  Figure 11. Clustering results of θ3 data in the XZ and XY planes: (a) angle data of θ3 in the XZ plane, when Sigmal was diagonal, and shared covariance was false; (b) angle data of θ3 in the XY plane, when Sigmal was diagonal, and shared covariance was false.
As can be seen from Figure 11, the data clustering effect of the θ3 angle in the XZ plane was not obvious and could not be clearly divided into four independent and complete regions. However, the clustering effect of the data in the XY plane was obvious, that is, the data could be accurately divided according to four slopes. Therefore, it was concluded that the change of θ3 in the XY plane was more significant, that is, the motion of goat lumbar vertebrae was mainly manifested as left-right swing in the coronal plane. Figure 11. Clustering results of θ 3 data in the XZ and XY planes: (a) angle data of θ 3 in the XZ plane, when Sigmal was diagonal, and shared covariance was false; (b) angle data of θ 3 in the XY plane, when Sigmal was diagonal, and shared covariance was false.
As can be seen from Figure 11, the data clustering effect of the θ 3 angle in the XZ plane was not obvious and could not be clearly divided into four independent and complete regions. However, the clustering effect of the data in the XY plane was obvious, that is, the data could be accurately divided according to four slopes. Therefore, it was concluded that the change of θ 3 in the XY plane was more significant, that is, the motion of goat lumbar vertebrae was mainly manifested as left-right swing in the coronal plane.

Analysis of the Variation Range of θ 2 and θ 3 in the XY Plane
The motion of thoracic vertebrae and lumbar vertebrae was mainly manifested as left-right swing in the coronal plane, and the variation range of θ 2 and θ 3 reflected the influence of thoracic vertebrae and lumbar vertebrae on the body of goats at different slopes. The larger the range of variation, the greater the regulation effect of this part on the movement of goats. The variation range of θ 2 and θ 3 in the XY plane under different slopes is shown in Figure 12.
Biomimetics 2022, 7, x FOR PEER REVIEW 16 of Figure 12. The variation range of θ2 and θ3.

Construction of Spine Spatial Curve Models for Eight Movements at 10°
Matlab software was used to fit the spatial curve polynomial of the eight movemen of the goat spine under 10° slope. The spinal space curve was generated by combining th two-dimensional fitting curves in the coronal plane and sagittal plane. The data of C, L1, L2, and L3 points (x, y, z) were imported into Matlab software, and the curve of th spine in the coronal plane was fitted with x as the independent variable and y as the d pendent variable. Conversely, the curve of spine in sagittal plane was fitted with x as th independent variable and z as the dependent variable. Finally, the general equation spine space fitting curve was obtained after analyzing the fitting results. The fitting resul of different indices of x in the coronal plane and sagittal plane are shown in Tables 1 an 2, respectively.  It can be seen from Figure 12 that the variation range of θ 2 in the XY plane was essentially stable at 20-23 • , while the variation range of θ 3 in the XY plane varied greatly. As the slope increased to 10 • , the variation range of θ 3 in the XY plane increased gradually to 25 • , and the variation range decreased when the slope increased to 15 • . A comprehensive analysis of the variation range of θ 2 and θ 3 in the XY plane showed that the variation range of θ 2 and θ 3 in the XY plane was the largest at the slope of 10 • , which meant that the coupling motion of thoraco-lumbar cooperative motion and flexible swing of lumbar vertebrae at the slope of 10 • had the most significant regulating effect.

Construction of Spine Spatial Curve Models for Eight Movements at 10 •
Matlab software was used to fit the spatial curve polynomial of the eight movements of the goat spine under 10 • slope. The spinal space curve was generated by combining the two-dimensional fitting curves in the coronal plane and sagittal plane. The data of C, T, L1, L2, and L3 points (x, y, z) were imported into Matlab software, and the curve of the spine in the coronal plane was fitted with x as the independent variable and y as the dependent variable. Conversely, the curve of spine in sagittal plane was fitted with x as the independent variable and z as the dependent variable. Finally, the general equation of spine space fitting curve was obtained after analyzing the fitting results. The fitting results of different indices of x in the coronal plane and sagittal plane are shown in Tables 1 and 2, respectively. Because the fitting result of the spinal curve in the coronal plane was uncertain when the independent variable index was four, the fitting effect of the spinal curve was analyzed when the independent variable index increased from one to three. As can be seen from Table 1, with the increase in independent variable index, the SSE of the spine curve fitting in the coronal plane gradually decreased and the R 2 increased. In addition, the influence of the number of independent variables included in the equation should be considered. The Adjusted-R 2 of group 2 and 3 increased with the increase in independent variables. When n was 3, SSE were 159.53, 318.17, and 0.33; R 2 were 0.9760, 0.9881, and 1.0000, Adjusted-R 2 were 0.9041, 0.9524, and 0.9999; and RMSE were 12.631, 17.837, and 0.573 for groups 1, 2, and 3, respectively. It met the fitting requirements.
As can be seen from Table 2, with the increase in the independent variable index, SSE and RMSE decreased, while R 2 and Adjusted-R 2 constantly increased, indicating that with the increase in the index of x, the fitted curve was closer to the curve characteristics of the actual spine. When the independent variable index was 2 or 3, R 2 was above 0.9885, and Adjusted-R 2 was above 0.9770. When the index of x increased from 2 to 3, SSE decreased by 127.81, 55.23, and 33.44; R 2 increased by 0.0111, 0.0049, and 0.0067; Adjusted-R 2 increased by 0.0212, 0.0008, and 0.0124; and RMSE decreased by 5.876, 4.5178, and 3.8783 for groups 1, 2, and 3, respectively. However, with the increase in the index of independent variables, the complexity of the spinal space curve increased. When n was 3, SSE were 159.53, 318.17, and 0.33; R2 were 0.9760, 0.9881, and 1.0000; Adjusted-R2 were 0.9041, 0.9524, and 0.9999; and RMSE were 12.631, 17.837, and 0.573 for groups 1, 2, and 3, respectively. It met the fitting requirements. Therefore, the general equation of fitting curve of goat spine in coronal plane was y = f(x) = ax 3 + bx 2 + cx + d (4) the general equation of fitting curve in sagittal plane was the spatial fitting curve equations of goat spine could be expressed as The independent variable x, dependent variables y and z, parameter names, and functions were set through the function editor function of Origin software, and the setup function was compiled. The nonlinear fitting function was used to fit the spine curve, and the equation parameters A, B, C, D, E, F, and G were obtained. As three groups of data were selected for fitting, the curve equations with the best fitting accuracy were selected to represent the spine curves. The fitting results under the eight movements are shown in Table 3. As can be seen from Table 3, the Adjusted-R 2 of the second group of data was smaller than that of the first and third groups, while the SSE was larger, indicating that the fitting effect of the second group of data was the worst. A comprehensive comparison of the Adjusted-R 2 and SSE of the first group and the third group showed that the Adjusted-R 2 of the two groups were above 0.9888. The SSE of the first group of 1 -6 movements were higher than that of the third group, and the SSE of the first group of 7 and 8 movements were lower than that of the third group. Therefore, the third set of fitting equations was chosen to represent the spinal space curves. The specific curve equations of the eight movements were as follows: 1 The left front leg swung monophasically, 2 The left front leg and the right hind leg swung biphasically, 3 The right hind leg swung monophasically, 4 The right hind leg and the right front leg swung biphasically, 5 The right front leg swung monophasically, 6 The right front leg and the left hind leg swung biphasically, 7 The left hind leg swung monophasically, 8 The left hind leg and the left front leg swung biphasically,

Discussion
In this study, the QTM gait analysis system was used to analyze the motion characteristics of the spine at different slopes in combination with the typical gait cycle of the goat. Compared with the traditional spine biological characteristics analysis method, it can obtain the changes of goat spine under real movement state, which provides a reference for the subsequent study of goat movement mechanism.
(1) The change of CTL1 represented the motion of cervical vertebrae relative to thoracic vertebrae, the change of TL1L2 represented the motion of thoracic vertebrae relative to lumbar vertebrae, and the change of L1L2L3 represented the internal motion of lumbar vertebrae, which were labeled as θ 1 , θ 2 , and θ 3 , respectively, for the convenience of subsequent research. On the basis of LTN, the sequence calibration and matching were performed. Moreover, the "left front-right hind-right front-left hind" leg lifting sequence was used as a typical gait cycle to analyze the changes of θ 1 , θ 2 , and θ 3 under different slopes. It was found that the changes of θ 1 , θ 2 , and θ 3 of the eight movements had a general law, which indicated that the spine motion characteristics were significant in the motion process of the goat under different slopes. It had a regulating effect on the stability of the goat motion.
(2) The maximin and variation range of θ 1 , θ 2 , and θ 3 under 0 • , 5 • , 10 • , and 15 • slopes were analyzed and compared. It was found that the maximin of θ 2 was greater than that of θ 3 at different slopes, and the variation range of θ 2 changed from larger than that of θ 3 to smaller than that of θ 3 with the increase in slope. That is to say, lumbar vertebrae played an increasingly important role in regulation. Spinal movement under large slopes was mainly manifested as a coupled motion of thoraco-lumbar coordination motion and flexible swing of lumbar vertebrae. Moreover, the head played an important role in regulating the motion stability of the goat under large slopes.
(3) According to the similarity between the data of anatomic plane and the spatial data, the Gaussian mixture model clustering algorithm was used to analyze the spatial variation of θ 2 in thoracolumbar vertebrae and θ 3 in a lumbar intervertebral body. It was found that when the covariance matrix was diagonal and the covariance matrix was different, the clustering effect of angle variation data under different slopes was obvious. The motion of thoracolumbar vertebrae and lumbar vertebrae was mainly manifested as left-right swing in the coronal plane.
(4) According to the range of data variation, the regulation effect of each part of the spine on goat movement was studied, and the range of angle variation in the coronal plane of θ 2 and θ 3 was comprehensively analyzed. It was found that the variation range of θ 2 in the coronal plane was essentially stable at 20-23 • when the slope changed from 0 • to 15 • , and the variation range of θ 3 in the coronal plane was about 25 • at 10 • . The results showed that the coupling motion of thoracolumbar joint motion and flexible swing of lumbar vertebrae at 10 • had the most significant effect on the kinematic stability of the goat.
The previous research of the team found that goats have excellent adhesion performance and can move stably in the slope environment. At present, only the motion characteristics of the spine under 0 • , 5 • , 10 • , and 15 • slopes in combination with the structural characteristics of the goat spine were analyzed in this study. Moreover, combined with the motion characteristics of the spine, the optimal fitting space equations were established under the eight movements. However, the best fitting curve based on the adhesion coefficient judgment still needs to be determined to further design the multidirectional flexible bending bionic spine. Limited by the experimental conditions, only laboratory slope tests have been conducted to obtain the laboratory data of goat climbing movement. However, the characteristics of goat spine movement under a real slope environment and larger slopes need to be further studied. The excellent locomotor stability of the goat is the coordinated action of spine and limbs. The further interaction mechanism of the spine and limbs remains to be studied. In future research, we will use high-speed camera technology to study the movement state of goats in a real environment and analyze the coordination between the spine and limbs on the basis of image research.

Conflicts of Interest:
The authors declare no conflict of interest, and the funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.