CN114781432A - A Displacement Solution Method Based on Multi-source Information Fusion and Detrended Fluctuation Analysis - Google Patents

Abstract

The invention provides a displacement resolving method based on multi-source information fusion and trend-removing fluctuation analysis, which can be applied to the fields of industrial automation, building bridges, earthquake monitoring, rehabilitation and nursing and the like. The method includes the steps that acceleration and angle data collected by a three-axis attitude sensor are utilized, space pose transformation is utilized to realize data fusion, the influence of gravity acceleration components in the acceleration data is eliminated, linear interpolation functions, a trend-removing fluctuation analysis method and a Kalman filtering algorithm are utilized to obtain acceleration data with equal time intervals, trend terms and noise in the acceleration data are eliminated, and finally a stable numerical integration method is utilized to calculate displacement. According to the method, on the premise that an external displacement sensor is not adopted, the actual displacement is obtained by resolving through fusing self acceleration and angle attitude information, the accurate solution of the spatial displacement is realized, and the defects that the displacement sensor can only measure one-dimensional displacement and the installation position space is limited are overcome.

Description

A Displacement Solution Method Based on Multi-source Information Fusion and Detrended Fluctuation Analysis

technical field

The invention relates to the technical field of sensors and signal processing, in particular to a displacement solution method based on multi-source information fusion and detrended fluctuation analysis.

Background technique

With the continuous development of sensor technology, various sensors are continuously used in various fields. Among them, the application range of displacement sensors is quite extensive, and they are often used in industrial automation and building bridges. The displacement sensor outputs electrical signals of different sizes according to the magnitude of the displacement, and then solves the magnitude of the displacement. It is mostly used for the length, distance, vibration, speed of equipment such as injection molding machines, hydraulic presses, hardware machinery, steel mill rolling rod adjustment, and shield machines. , azimuth and other physical quantities measurement. However, due to the large size of displacement sensors, it is not easy to install and use in smaller equipment, the measurement distance is easily limited, and the two-dimensional and three-dimensional motion states cannot be measured. Constraints on the state of the environment. In addition to directly obtaining the displacement signal through the displacement sensor, the acceleration signal can also be obtained through an acceleration sensor with a wider application field, and the displacement signal can be obtained by integrating the obtained acceleration signal. In the existing methods of obtaining displacement by integrating acceleration, the acceleration sensor is mostly fixed on the test platform, and the obtained acceleration data is used for analysis and processing. For acceleration sensors with vertical axis tilt, most of the existing methods use a high-pass filter to adjust the baseline, which not only removes the baseline error, but also eliminates the low-frequency content in the signal including residual displacement, so that the integral result is smaller than the true value as a whole. Therefore, in view of the characteristics of the acceleration sensor and numerical integration, the present invention provides a stable numerical integration method based on multi-source information fusion and detrended fluctuation analysis, which can be separated from the test bench to realize the movement and portability of the sensor, and is used to realize the acceleration sensor one. The displacement solution of 2D, 2D and 3D motion is the key to further solve the equipment displacement information under the condition of limited position space.

SUMMARY OF THE INVENTION

In order to solve the problem of realizing one-dimensional, two-dimensional and three-dimensional displacement solutions according to the time information, acceleration information and angle information output by the acceleration sensor when only the acceleration sensor is used, the present invention proposes a method combining information fusion and detrending fluctuation. The analytical displacement solution method provides a new idea and method for the displacement solution of equipment that is not suitable for using displacement sensors due to the limitation of position space.

A displacement calculation method based on multi-source information fusion and detrended fluctuation analysis proposed by the present invention is specifically:

Step 1: Detect the obtained attitude data, perform data preprocessing, and delete multiple duplicate data at the same time.

Step 2: Perform data fusion according to the obtained acceleration and spatial attitude information to eliminate the gravity acceleration component in the obtained acceleration.

Step 3: Perform interpolation processing according to the time data and the acceleration data to obtain acceleration data at equal time intervals.

Step 4: Eliminate the trend component in the acceleration according to the detrended fluctuation analysis technique.

Step 5: Eliminate the noise in the acceleration according to the Kalman filter method.

Step 6: According to the stable numerical integration method, integrate the acceleration data to solve the displacement.

The beneficial effects of the present invention are:

(1) Using a smaller acceleration sensor, compared with the displacement sensor, it is less affected by the position space and easy to transplant, and can be combined with a micro-miniature wireless transmission module to realize wireless communication with the host computer.

(2) Combined with the fusion algorithm of acceleration and corresponding angle information, the influence of gravitational acceleration is eliminated from the aspect of pose transformation, and the influence of the change of gravitational acceleration component caused by the spatial attitude transformation of the acceleration sensor during the movement process is eliminated.

(3) Using detrended fluctuation analysis and Kalman filtering method, the drift and noise of the acceleration sensor during long-term work are eliminated, so that the result of solving the displacement during long-term work is more accurate.

(4) From the transfer function, the reasons for the drift of the traditional numerical integration method in the integration process are considered, and the improved integration algorithm becomes a stable system by modifying the integration parameters.

Description of drawings

Fig. 1 is the displacement solution flow chart of the present invention;

2 is a schematic diagram of the Cartesian coordinate system pose transformation of the acceleration sensor of the present invention;

Fig. 3 is the acceleration data and angle data collected when the acceleration sensor is used in the present invention to move along the X-axis;

Fig. 4 is the acceleration data and angle data collected when the acceleration sensor is used in the present invention to move along the Y-axis;

Fig. 5 is the acceleration data and angle data collected when the acceleration sensor is used in the present invention to move along the X axis and the Y axis with an angle of -61°;

Fig. 6 is the acceleration data and the angle data collected when the acceleration sensor is used in the present invention to move along the X axis and the Y axis angle of 53°;

7 is the acceleration data collected when the acceleration sensor is used in the present invention to move along the X axis and the acceleration data processed by the method of the present invention;

8 is the acceleration data collected when the acceleration sensor is used in the present invention to move along the Y axis and the acceleration data processed by the method of the present invention;

9 is the acceleration data collected when the acceleration sensor is used in the present invention to move along the X axis and the Y axis at an angle of -61° and the acceleration data processed by the method of the present invention;

10 is the acceleration data collected when the acceleration sensor is used to move along the X-axis and the Y-axis at an angle of 53° and the acceleration data processed by the method of the present invention;

11 is a comparison diagram of the integral displacement and the actual displacement obtained by the method of the present invention when the acceleration sensor moves along the X axis in the present invention;

12 is a comparison diagram of the integral displacement and the actual displacement obtained by the method of the present invention when the acceleration sensor moves along the Y axis in the present invention;

13 is a comparison diagram of the integral displacement and the actual displacement obtained after the calculation by the method of the present invention when the acceleration sensor moves along the X axis and the Y axis at an angle of -61° in the present invention;

14 is a comparison diagram of the integral displacement and the actual displacement obtained after calculation by the method of the present invention when the acceleration sensor moves along the X axis and the Y axis at an angle of 53° in the present invention;

15 is a three-dimensional diagram of the displacement trajectory obtained by the method of the present invention after the acceleration data collected when the acceleration sensor is used to move in four different directions in the present invention is calculated.

Detailed ways

The displacement calculation method based on information fusion and detrended fluctuation analysis of the present invention will be described in detail below with reference to the accompanying drawings and embodiments.

Step 1: Detect the time data in the obtained data samples, delete the duplicate data groups in the process of data transmission and reception, retain only a group of valid data at the same time, and process the time data as absolute time, taking the first group of data as Start data, time node is set to 0.

Step 2: Introduce the concept of rotation matrix in robot kinematics to eliminate the gravitational acceleration component contained in the original acceleration data. As shown in Figure 2, the acceleration sensor is designed with a Cartesian coordinate system. According to the right-hand rule, the positive direction of rotation around the axes X ₀ , Y ₀ , and Z ₀ can be determined, and the rotation angles are expressed as α, β, and γ, respectively.

Rotate α, β, γ around the X ₀ , Y ₀ , Z ₀ axes respectively, and the obtained rotation matrices are:

From the initial coordinate system 0-X ₀ Y ₀ Z ₀ , first rotate α around the X ₀ axis, then rotate β around the Y ₀ axis, and finally rotate γ around the Z ₀ axis to obtain the real pose coordinate system 0-XYZ during the actual operation. . The relationship between the final rotation matrix R and the two coordinate systems are:

Acc _init = R·Acc _gra (5)

Among them, Acc _init is the gravitational acceleration component of the gravitational acceleration on each axis in the world coordinate system, and Acc _gra is the gravitational acceleration component of the gravitational acceleration on each axis in the actual pose coordinate system. Due to the influence of gravity, when the acceleration sensor is stationary in the world coordinate system, there is always an acceleration whose magnitude is g and whose direction coincides with the positive direction of the Z ₀ axis, that is, Acc _init = [00g] ^T . Then in the actual pose coordinate system, the gravitational acceleration component of the gravitational acceleration in each axis can be expressed as:

Then in the actual pose coordinate system, the acceleration of each axis after removing the gravitational acceleration component is:

Acc _act = Acc - Acc _gra = [Acc _x +gsinβAcc _y -gcosβsinαAcc _z -gcosαcosβ] ^T (7)

Wherein Acc=[Acc _x Acc _y Acc _z ] ^T , Acc _x , Acc _y , and Acc _z respectively represent the collected three-axis acceleration data.

Step 3. The time intervals of the acceleration data after eliminating the gravitational acceleration component may not be equal. Interpolate the known acceleration data and the time data to ensure that the acceleration data in the subsequent integration process have the same time interval and can be To a certain extent, the simplicity of the integration process and the accuracy of the integration results are guaranteed. Using time data and acceleration data for piecewise linear interpolation, the interpolation function φ(x) in each small interval [x _k , x _k+1 ] can be expressed as:

Among them, x _k (k=0,1,...,n) represents the time node, and f _k (k=0,1,...,n) represents the acceleration at the time of x _k . Using the interpolation basis function l(x) to represent the interpolation function over the entire time interval [t _s , _te ] is:

Among them, l _i (x) is the interpolation basis function of the i-th time node (i=0,1,...,n), which can be expressed as:

Among them, t _s and t _e represent the start time and end time of the entire time interval, respectively. The equal interval time and the corresponding acceleration data can be determined as required by the linear interpolation function.

Step 4. The data output by the accelerometer during long-term operation contains certain trend components. If the trend components are not filtered out during analysis, the actual data will gradually deviate from the true value with the increase of time, and will be in the subsequent integration process. cause greater deviations. Therefore, the detrended fluctuation analysis technique is used to remove the trend components in the acceleration data for the interpolated acceleration data. The detrended volatility analysis method process is as follows:

For a sequence of length N {φ _j ,j=1,2,...,N}, create a new sequence:

为原序列φ_j的均值。in,

is the mean value of the original sequence φ _j .

The new sequence y(m) is divided into non-overlapping equal-length sub-intervals of length s, and the sequence of length N is divided into N _s =int(N/s) sub-intervals. Because the sequence length N is not necessarily divisible by the subinterval length s, in order to ensure that the original sequence information is not lost, it can be divided again from the end of the sequence before the reverse, and a total of 2N _s subintervals are obtained.

Perform polynomial regression _fitting on the data of each sub-interval v (v=1, 2, . Eliminate the trend in each sub-interval and calculate its variance mean F ² (v,s) as:

Determine the q-order wave function for the full sequence:

Among them, q can take any non-zero real number. When q=0, equation (14) becomes:

In actual calculation, when performing polynomial regression fitting on the data of each sub-interval v (v=1, 2, . . . , 2N _s ), q can be taken as any non-zero real number.

Step 5: The original acceleration signal contains noise generated in the process of data acquisition, and Kalman filtering algorithm is used to process the noise signal. The algorithm has the advantages of small computation and good real-time performance. The estimated value of the future motion state can be continuously revised by using the actual motion parameters, and the estimation accuracy can be improved on the premise of taking into account the real-time performance and stability.

The state vector prediction equation can be expressed as:

The state vector covariance matrix prediction can be expressed as:

The Kalman gain matrix can be expressed as:

The state vector update equation can be expressed as:

The state vector covariance update equation can be expressed as:

表示预测状态，

表示状态估计，A表示状态转移矩阵，

表示预测误差协方差，P_k表示估计误差协方差，Q表示过程误差，K_k表示卡尔曼增益，H表示测量矩阵，R表示测量误差，a_k表示测量值，I表示单位矩阵。in,

represents the predicted state,

represents the state estimation, A represents the state transition matrix,

is the prediction error covariance, P _k is the estimation error covariance, Q is the process error, K _k is the Kalman gain, H is the measurement matrix, R is the measurement error, a _k is the measurement value, and I is the identity matrix.

Step 6: Use the trapezoidal method to integrate the acceleration data processed above. For the time continuous function x(t) whose sampling interval is T, the trapezoidal integration formula in the interval [0, KT] is:

y represents the value obtained by integrating the time-continuous function x(t) with the sampling interval T. The transfer function of the trapezoidal integral formula is:

From the transfer function described in formula (22), it can be known that its pole is located on the unit circle, that is, it is in a critically stable state, and it is easy to cause divergence during the integration process. The formulas of velocity v(j) and displacement s(j) obtained by integrating the jth (j=1,...,N) acceleration a(j) are expressed as:

v(j)=T[a(j)+a(j-1)]/2+ωv(j-1) (23)

s(j)=T[v(j)+v(j-1)]/2+ωs(j-1) (24)

where ω is the stabilization factor.

Example

The method of the present invention is used to perform stable numerical integral calculation on the data collected by the three-axis acceleration sensor as follows:

The test instrument with built-in three-axis accelerometer is used to perform four groups of cyclic reciprocating motions, which move along the X-axis, Y-axis and X-axis, and the included angles of the Y-axis are -61° and 53°. The test instrument sends the time through the wireless module. data, acceleration data, and angle data. Check the time data, and delete the duplicate data to obtain the acceleration data and angle data as shown in Figure 3-Figure 6. In order to verify the correctness of the method of the present invention, a laser displacement sensor is used to collect actual displacement data for reference.

For the obtained triaxial acceleration data and the corresponding triaxial angle data, according to formula (7), the triaxial acceleration data after removing the gravitational acceleration component is obtained. Since the interval of time data obtained by the acceleration sensor is not equal, linear interpolation is used to obtain three-axis acceleration interpolation data with a time interval of T=0.01s. The trend component in the interpolated data is then removed using a detrended fluctuation analysis technique, where q=15. Then the Kalman filter algorithm with process error Q=0.001 and measurement error R=0.01 is used for final processing of the acceleration data. Figure 7-Figure 10.

According to the obtained final acceleration signal, the stable numerical integration algorithm of formula (23) and formula (24) is used to solve the displacement, where the stability factor ω=0.996, which is compared with the displacement data obtained by the laser displacement sensor. The integral result and the actual displacement result are shown in the figure 11—Figure 14, the three-dimensional results of the integration results are shown in Figure 15.

The peak error, difference error, absolute error and mean square error are used to evaluate the error of the acceleration integration result using the method of the present invention.

The peak error err_peak represents the average value of the relative difference between the peak value of the integration result and the peak value of the true value:

The difference error err_diff represents the average value of the difference between the result of the integration and the true value:

The absolute error err_abs represents the ratio of the difference between the integral result and the true value to the true value:

The difference in mean square error MSE represents the degree of difference between the integration result and the true value:

Among them, s(t) represents the integration result, s _t (t) represents the true value, and the displacement data collected by the laser displacement sensor is the true value by default.

According to formula (25) - formula (28), the calculation error results are shown in Table 1.

Table 1 Error Results

It can be seen from Table 1 that, applying the method of the present invention to process the acceleration, the accuracy of the integration result is obviously improved, and the average of err_peak, err_diff, err_abs and MSE is increased by 91.68%, 63.41%, 72.40% and 86.83% respectively.

Claims (10)

1. a displacement solution method based on multi-source information fusion and detrended fluctuation analysis, is characterized in that, comprises the steps as follows: Step 1: Detect the obtained attitude data, perform data preprocessing, and delete multiple duplicate data at the same time; Step 2: Perform data fusion according to the obtained acceleration and spatial attitude information to eliminate the gravitational acceleration component in the obtained acceleration; Step 3: Perform interpolation processing according to time data and acceleration data to obtain acceleration data at equal time intervals; Step 4: Eliminate the trend component in the acceleration according to the detrended fluctuation analysis technology; Step 5. Eliminate the noise in the acceleration according to the Kalman filter method; Step 6: According to the stable numerical integration method, integrate the acceleration data to solve the displacement. 2. a kind of displacement solution method based on multi-source information fusion and detrended fluctuation analysis according to claim 1, is characterized in that: in step 1, be specifically: detect the time data in the data sample obtained, delete In the process of data transmission and reception, only a group of valid data at the same time is retained, and the time data is processed as absolute time. 3. a kind of displacement calculation method based on multi-source information fusion and detrended fluctuation analysis according to claim 1, is characterized in that: in step 2, be specifically: introduce the concept of rotation matrix in robot kinematics to eliminate original The gravitational acceleration component contained in the acceleration data; the acceleration sensor is designed with a Cartesian coordinate system, and the positive direction of rotation around the X ₀ , Y ₀ , and Z ₀ axes is determined according to the right-hand rule, and the rotation angles are respectively expressed as α, β, γ; Rotate α, β, γ around the X ₀ , Y ₀ , Z ₀ axes respectively, and the obtained rotation matrices are:

From the initial coordinate system 0-X ₀ Y ₀ Z ₀ , first rotate α around the X ₀ axis, then rotate β around the Y ₀ axis, and finally rotate γ around the Z ₀ axis to obtain the real pose coordinate system 0-XYZ during the actual operation. ; The relationship between the final rotation matrix R and the two coordinate systems are:

Acc _init = R·Acc _gra (5) Among them, Acc _init is the gravitational acceleration component of the gravitational acceleration on each axis in the world coordinate system, and Acc _gra is the gravitational acceleration component of the gravitational acceleration on each axis in the actual pose coordinate system. 4. a kind of displacement calculation method based on multi-source information fusion and detrended fluctuation analysis according to claim 3, is characterized in that: in step 2, also further comprises: because of the influence of gravity, the acceleration sensor is in the world coordinate When stationary in the system, there is always an acceleration whose size is g, and the direction coincides with the positive direction of the Z ₀ axis, that is, Acc _init = [00g] ^T ; then in the actual pose coordinate system, the gravitational acceleration component of the gravitational acceleration in each axis is expressed as :

Then in the actual pose coordinate system, the acceleration of each axis after removing the gravitational acceleration component is: Acc _act = Acc - Acc _gra = [Acc _x +g sinβ Acc _y -g cosβsinα Acc _z -g cosαcosβ] ^T (7) Wherein Acc=[Acc _x Acc _y Acc _z ] ^T , Acc _x , Acc _y , and Acc _z respectively represent the collected three-axis acceleration data. 5. a kind of displacement solution method based on multi-source information fusion and detrended fluctuation analysis according to claim 1, is characterized in that: in step 3, be specifically: the time interval of the acceleration data after eliminating the gravitational acceleration component It may not be equal. Interpolate the known acceleration data and time data to ensure that the acceleration data in the subsequent integration process have the same time interval, and to a certain extent, ensure the simplicity of the integration process and the accuracy of the integration results; Using time data and acceleration data for piecewise linear interpolation, the interpolation function φ( _{x )} is expressed as:

Among them, x _k (k=0,1,...,n) represents the time node, and f _k (k=0,1,...,n) represents the acceleration at the time of x _k ; the interpolation basis function l( x) represents that the interpolation function over the entire time interval [t _s , _te ] is:

Among them, l _i (x) is the interpolation basis function of the i-th time node (i=0,1,...,n), which is expressed as:

Among them, t _s and _te represent the start time and end time of the entire time interval, respectively; the equal interval time and the corresponding acceleration data are determined by the linear interpolation function as required. 6. a kind of displacement solution method based on multi-source information fusion and detrended fluctuation analysis according to claim 1, it is characterized in that: in step 4, be specifically: the data that the acceleration sensor outputs during long-term operation contains a certain trend component, and it is necessary to use the detrended fluctuation analysis technique to remove the trend component in the acceleration data for the interpolated acceleration data; the process of the detrended fluctuation analysis method is as follows: For a sequence of length N {φ _j ,j=1,2,...,N}, create a new sequence:

in,

is the mean value of the original sequence φ _j ; Divide the new sequence y(m) into non-overlapping equal-length sub-intervals of length s, and the sequence of length N is divided into N _s =int(N/s) sub-intervals; because the sequence length N is not necessarily determined by the length of the sub-intervals Divide by s, in order to ensure that the original sequence information is not lost, it is divided again from the end of the sequence before the reverse, and a total of 2N _s subintervals are obtained. 7. a kind of displacement solution method based on multi-source information fusion and detrended fluctuation analysis according to claim 6, is characterized in that: in step 4, also further comprises: to each sub-interval v (v=1, 2,...,2N _s ) data are subjected to polynomial regression fitting, and the obtained local trend function is a first-order, second-order or higher-order polynomial; eliminate the trend in each sub-interval, and calculate its variance mean F ² (v , s) is:

Determine the q-order wave function for the full sequence:

where q takes any non-zero real number. 8. a kind of displacement solution method based on multi-source information fusion and detrended fluctuation analysis according to claim 7, is characterized in that: in step 4, also further comprises: when q=0, formula (14) becomes:

In actual calculation, when performing polynomial regression fitting on the data of each sub-interval v (v=1, 2, . . . , 2N _s ), q takes any non-zero real number. 9. a kind of displacement solution method based on multi-source information fusion and detrended fluctuation analysis according to claim 1, is characterized in that: in step 5, is specifically: the original acceleration signal contains the data generated in the process of data acquisition. Noise, the Kalman filter algorithm is used to process the noise signal; this algorithm has the advantages of small calculation amount and good real-time performance; using the actual motion parameters, the estimated value of the future motion state is continuously revised, taking into account the real-time performance and stability. Improve estimation accuracy; The state vector prediction equation is expressed as:

The state vector covariance matrix prediction is expressed as:

The Kalman gain matrix is expressed as:

The state vector update equation is expressed as:

The state vector covariance update equation is expressed as:

in,

represents the predicted state,

represents the state estimation, A represents the state transition matrix,

is the prediction error covariance, P _k is the estimation error covariance, Q is the process error, K _k is the Kalman gain, H is the measurement matrix, R is the measurement error, a _k is the measurement value, and I is the identity matrix. 10. A displacement calculation method based on multi-source information fusion and detrended fluctuation analysis according to claim 1, characterized in that: in step 6, specifically: adopting a trapezoidal method to perform the above-mentioned processing on the acceleration data. Integral, for the time continuous function x(t) whose sampling interval is T, the trapezoidal integration formula in the interval [0, KT] is:

Among them, y represents the value obtained after integrating the time continuous function x(t) with the sampling interval T; the transfer function of the trapezoidal integration formula is:

Obtained from the transfer function described in formula (22), its pole is located on the unit circle, that is, it is in a critically stable state, and it is easy to cause divergence in the integration process; for the jth (j=1,...,N) acceleration a( j) The velocity v(j) and displacement s(j) obtained by integration are expressed as: v(j)=T[a(j)+a(j-1)]/2+ωv(j-1) (23) s(j)=T[v(j)+v(j-1)]/2+ωs(j-1) (24) where ω is the stabilization factor.

Images