CN104142624A - Time synchronization method and system based on waveform matching - Google Patents

Abstract

The present invention proposes a time synchronization system and method based on waveform matching. The present invention uses a set of data containing reference time to synchronize accurate time information for a set of data without time stamps or with wrong time stamps through the principle of correlation. . In the present invention, the data to be matched is segmented first, and the correlation between each segment and different sub-segments of the reference data is solved. When the correlation coefficient is the largest, the time difference between the data of this segment and the reference data is calculated, thereby obtaining the time difference of all segments Then, the time difference sequence is used as the observation value, the time difference and time scale drift factor are the parameters to be estimated, and the least squares adjustment is used to solve the parameters to be estimated; finally, the time difference and time scale drift factor are used to synchronize the precise time for the data to be matched. The present invention can be widely applied to time synchronization during multi-sensor data fusion, and can also be used to detect the correctness of time synchronization by other methods.

Description

A time synchronization method and system based on waveform matching

technical field

The invention relates to the field of sensor applications, in particular to a time synchronization method and system based on waveform matching.

Background technique

At present, the research on multi-sensor data fusion technology is in the ascendant, and some mature technologies have been applied to practical projects and achieved good results. The application of multi-sensor data fusion technology can not only improve the accuracy and reliability of the system, but also improve the measurement range of the system, increase the reliability of the system, and shorten the reaction time of the system. However, multi-sensor data fusion is a complex uncertain information processing process. The prerequisite for fusion is that the information obtained from each sensor must be a description of the same target at the same time. This includes two aspects. First, it is necessary to ensure that the information obtained by each sensor is a description of the same parameter of the same target. Second, it is necessary to ensure that the time information of the data to be fused is synchronized. In a dynamic work environment, time synchronization issues are particularly prominent.

In different engineering practices, problem-specific time synchronization methods have emerged. For example, the time synchronization algorithm using the curve fitting method, the time synchronization algorithm using the sequence method to solve the time-varying observation, and the smoothing filter algorithm to synchronize the time synchronization of the measurement data between the sensors. Before these algorithms are applied, each sensor is time-stamped with the same time stamp at the beginning of the measurement, and then clocked according to its respective sampling rate. However, for the entire measurement process, each sensor does not have the same time scale, or the frequency scale drift of the processor itself due to the influence of temperature characteristics and other factors, the above algorithms cannot correctly realize the time between multiple sensors Synchronize.

A typical example of multi-sensor data fusion is the integrated navigation system composed of GPS (Global Positioning System, Global Positioning System) and INS (Inertial Navigation System, Inertial Navigation System). Among them, GPS is equipped with a high-precision atomic clock to provide accurate time, while INS can only count time through the set frequency, and the INS frequency standard often drifts. GPS/INS integrated navigation system usually adopts Kalman filter for data fusion, only when GPS and INS subsystem data are at the same time point, integrated navigation has practical significance. At present, the most common and effective method is to use the GPS second pulse signal as the time synchronization reference of the integrated navigation system, and the integrated structure design realizes time synchronization. However, when a new IMU (Inertial Measurement Unit, Inertial Measurement Unit) needs to be tested for integrated navigation, although the above method has high precision, it has high complexity, high cost, and requires hardware and software support, so it is not applicable. In the actual test process, a set of high-precision integrated navigation system is often used as a reference system, and the INS has been marked with high-precision GPS time scale.

The INS time scale in the integrated navigation system can be used as a reference time system. It only needs to design an appropriate algorithm to realize the time synchronization between the IMU to be tested and the reference system. This algorithm is not only suitable for time synchronization between inertial sensors, It can also be applied to time synchronization between other sensors with a reference time scale.

Contents of the invention

Aiming at the above problems, the present invention designs a time synchronization method and system based on waveform matching; this method and system have no restrictions on whether the original time systems of different sensor data are the same, and the operation is simple, the operation is fast, and the cost is low;

The technical solution of the present invention is: a time synchronization method based on waveform matching, comprising the following steps:

Step 1: Read the reference data R ₀ and the data to be matched T ₀ ;

Step 2: Unify the data processing format of the reference data R ₀ and the data to be matched T ₀ , and extract the converted reference data R ₁ and the data to be matched T ₁ ;

Step 3: Downsample R ₁ and T ₁ respectively, solve the cross-correlation sequence for the down-sampled reference data R ₂ and the data to be matched T ₂ , take the interval corresponding to the maximum value of the cross-correlation sequence, and find the current coarse time difference int_bias;

Step 4: According to the obtained rough time difference int_bias, make the time period of the reference data R ₂ cover the time period of the data to be matched T ₂ , and each exceed a certain length of the data to be matched T ₂ , intercept the reference data R ₂ and the data to be matched T ₂ ; enter step 5;

Step 5: Unify the sampling rates of R ₂ and T ₂ to obtain the reference data R ₃ and the data to be matched T ₃ of the same sampling rate;

Step 6: Calculate the optimal time difference bias and time scale drift factor k between R ₃ and T ₃ .

Step 7: According to the best combination of bias and k, update the time T ₀ of the data to be matched, and store it;

Step 8: End.

The formula for the data processing format of the unified reference data _R0 and the data to be matched _T0 in the step 2 is as follows:

v=Δs/Δt _s

Among them, Δs is the incremental value measured by the gyro or accelerometer between the two samples before and after, Δt _s is the time interval between the two samples before and after, and v is the rate value corresponding to Δs;

The converted reference data _R1 and the data to be matched _T1 are composed of converted time series and one-axis data for waveform matching;

Described step 3 comprises the following steps:

Step 3.1: Downsample _R1 and _T1 respectively to obtain the downsampled reference data _R2 and the data to be matched _T2 ;

Step 3.2: After the downsampling, the shorter sequence of R ₂ and T ₂ is zero-filled until the two lengths are equal;

Step 3.3: Calculate the cross-correlation sequence of R ₂ and T ₂ ; the calculation formula of the cross-correlation sequence is as follows:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; R</span>}}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\underset{\mathrm{<span\; class="notranslate">\; K</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; K</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; \&Greater\; Equal;</span>\mathrm{<span\; class="notranslate">\; 0</span>}\\ {\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; R</span>}}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right.$

Among them, x and y are the sequences for correlation processing, the lengths of x and y are both N, K and m are used to represent the ordinal numbers of the values in the sequence x and y respectively, and the value of m is -(N-1) to + (N-1);

Step 3.4: Take the cross-correlation sequence The interval corresponding to the maximum value m, and the rough time difference int_bias at this time is calculated, then the rough time difference can be obtained using the following formula:

int_bias=R ₂ (t ₀ )-T ₂ (t ₀ )-T ₂ _interval×m

Among them, R ₂ (t ₀ ) is the start time of R ₂ , T ₂ (t ₀ ) is the start time of T ₂ , T ₂ _interval is the time length represented by the staggered interval, which is the sampling interval of T ₂ , then T ₂ _interval ×m is the staggered time of two columns of data;

In the described step 6.15, the following sub-steps are included:

Step 6.1: Calculate the sub-sampling interval T_inter;

Step 6.2: Using T_inter as the multiple factor, calculate the δs of the translation distance to the right of the data to be matched;

Step 6.3: Determine whether δs is a positive integer multiple of the sampling interval; if so, go to step S6.15, if not, go to step S6.4;

Step 6.4: Translating the data T ₃ to be matched by δs to obtain T ₄ ;

Step 6.5: Take Δt as the window length, divide T ₄ into finite n segments, and take the current segment data T ₄ (i), 1≤i≤n;

Step 6.6: Determine whether the data to be matched is extracted to the nth segment; if so, go to step S6.13; if not, go to step S6.7;

Step 6.7: Set the sliding window W ₁ , the window length is Δt, and the starting point is the starting time of the current segment;

Step 6.8: Determine whether the sliding distance of the window is greater than the set length σ _s ; if yes, go to step S6.11; otherwise, take W ₁ as the window and take the reference data R ₃ (W ₁ );

Step 6.9: Calculate the Pearson correlation coefficient between T ₄ (i) and R ₃ (W ₁ ); the Pearson correlation coefficient is calculated as follows:

$\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; =</span>\frac{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span>}{\sqrt{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&Center\; Dot;</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}}$

in is the mean value of the x series, is the mean value of the y sequence, N is the length of the sequence x and y, and r is the Pearson correlation coefficient between the sequence x and y;

Step 6.10: Determine whether the reference data has read the last bit, if so, go to step S6.11; otherwise, slide the window W ₁ backward for a sampling interval of the reference data, and then go to step S6.8;

Step 6.11: Calculate the time difference Bias(i) and the corresponding time time(i);

Step 6.12: Extract the next segment of data to be matched. Take the next segment T ₄ (i) and go to step 6.6;

Step 6.13: Use the least squares adjustment to find the time difference bias and time scale drift factor k; take bias and k as the parameters to be estimated, Bias(i) as the observed value, and the observation equation is as follows:

Bias(i)=k×[time(i)-t ₀ ]+bias

Among them, the initial value of bias, Bias0, is the median value of the Bias(i) sequence, and the corresponding time is t ₀ , and the initial value of k is 0;

Step 6.14: Determine whether δs is greater than a sampling interval, if so, go to step S6.15; otherwise, double δs by T_inter, go to step S6.3;

Step 6.15: Select the best combination of bias and k; use different sets of bias and k to obtain data to be matched after synchronization at different times;

In the step 6.15, if the data is uniaxial, select the set of bias and k with the largest covariance as the best combination of bias and k; if the data is multi-axial, select the most credible one-axis data for covariance Solve, choose the best combination of bias and k.

In described step 7, utilize following formula to calculate and update T ₀ time:

t'=(tt ₀ )×k+bias+t

Among them, t is the time column of T ₀ , t ₀ is the time corresponding to the initial value of the time difference selected in step 6.12, k is the time scale drift factor, bias is the time difference, and t′ is the time after t is updated.

A time synchronization system based on waveform matching, characterized in that it includes the following modules:

Reading module: used to read reference data R ₀ and data to be matched T ₀ ;

Data processing format module: used to unify the data processing format of the reference data R ₀ and the data T ₀ to be matched, and extract the converted reference data R ₁ and the data T ₁ to be matched;

Downsampling module: it is used to downsample _R1 and _T1 respectively, solve the cross-correlation sequence for the reference data _R2 after downsampling and the data _T2 to be matched, and obtain the corresponding interval when the maximum value of the cross-correlation sequence is obtained. The coarse time difference int_bias at this time;

Coarse time difference module: used to make the time period of the reference data _R2 cover the time period of the data to be matched _T2 according to the obtained rough time difference int_bias, and each of the front and rear exceeds a certain length of the data to be matched _T2 , intercepting the reference data _R2 and Data to be matched T ₂ ;

Sampling rate calculation module: used to unify the sampling rates of R ₂ and T ₂ to obtain reference data R ₃ and data to be matched T ₃ of the same sampling rate;

Calculation module: used to obtain the optimal time difference bias and time scale drift factor k of R ₃ and T ₃ ;

Update module: used to update the data to be matched T ₀ time according to the best bias and k combination, and store it.

Description of drawings

Fig. 1 is a schematic flow sheet of the present invention;

Fig. 2 is the schematic flow chart of step 6 of the present invention;

Figure 3-1 is a schematic diagram of the angular rate waveform in the vertical direction of the EPSON gyroscope;

Figure 3-2 is a schematic diagram of the angular rate waveform in the vertical direction of the FSAS gyroscope;

Figure 4-1 is a schematic diagram of the maximum correlation coefficient obtained by segmenting the three axes of the EPSON gyroscope;

Figure 4-2 is a schematic diagram of the maximum correlation coefficient obtained by segmenting the three axes of the EPSON accelerometer;

Figure 5-1 is a schematic diagram of the time scale difference sequence solved by the three axes of the EPSON gyroscope;

Figure 5-2 is a schematic diagram of the time scale difference sequence solved by the three axes of the EPSON accelerometer;

Fig. 6 is a schematic diagram of the time matching result in the vertical direction of the gyroscope;

Fig. 7 is a partial schematic diagram of the time matching result in the vertical direction of the gyroscope;

Fig. 8 is a schematic diagram of the system structure of the present invention.

Detailed ways

The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments.

The embodiment is described by taking data collected by an inertial sensor. The inertial sensor can collect multi-axis data. This embodiment takes multi-axis data processing as an example.

Data source in the embodiment: during on-vehicle test, carry the GPS/INS integrated navigation system of low-accuracy IMU EPSON and containing high-accuracy IMU FSAS simultaneously, wherein in the data that EPSON gathers with counter timing, FSAS passes through and GPS in the integrated navigation system Integrated design, has been marked with accurate GPS time.

Data format in the embodiment: EPSON sampling rate is 125Hz, the data contains 7 columns, the first column is the time column, and the 2nd to 4th columns are the three columns of angular rates collected by the gyroscope, which are forward, right and vertical respectively , the 5th to 7th columns are the three columns of acceleration collected by the accelerometer, which are forward, right and vertical; Columns of angular velocity increments are forward, right and vertical, respectively. Columns 5 to 7 are three columns of acceleration increments collected by the accelerometer, which are forward, right and vertical.

The present invention can synchronize accurate time information for data without a time mark or with an inaccurate time mark. The implementation steps are shown in Figure 1, and Figure 2 shows the steps for solving the optimal time difference and time scale drift factor in Figure 1.

S1: Read reference data R ₀ and data to be matched T ₀ . The FSAS data is used as the reference data R ₀ , and the EPSON data is used as the data to be matched T ₀ .

S2: Unify the data processing format of the reference data R ₀ and the data to be matched T ₀ , extract the converted reference data R ₁ and the data to be matched T ₁ (useful data column); judge the relationship between R ₀ and Whether the parameters represented by T ₀ are the same. If they are different, the conversion will be performed according to the relationship between the parameters, and finally the commonly used data format shall prevail. In this example, EPSON is the rate format, while FSAS is the incremental format, and the commonly used data processing format is the rate format. The formula for converting increments to rates is as follows:

v=Δs/Δt _s

Among them, Δs is the incremental value measured by the gyro or accelerometer between the two samples before and after, Δt _s is the time interval between the two samples before and after, and v is the velocity value corresponding to Δs.

Use the above formula to convert R ₀ , and the corresponding time is the mean value of the two sampling times before and after. If one of the axis data is used for waveform matching, the converted time series and current axis data are extracted to form new reference data R ₁ . In this example, waveform matching is performed on the data of the 6 axes, the final FSAS data R ₁ is the converted data of R ₀ , and the EPSON data T ₁ is T ₀ .

Figure 3-1 and Figure 3-2 show the vertical gyro data waveforms in the fourth column of R ₁ and T _1. It can be seen from the figure that the two waveforms are highly consistent, and other axes also have the same characteristics, so that Waveform matching is performed according to the correlation coefficient of the two data.

S3: Downsample R ₁ and T ₁ , solve the cross-correlation sequence for the down-sampled reference data R ₂ and the data to be matched T ₂ , take the interval corresponding to the maximum value of the cross-correlation sequence, and find the rough time difference at this time int_bias;

First, the matching data and reference data are down-sampled. The frequency after down-sampling is the common factor of the original frequency of the two sets of data. On the basis of meeting the accuracy of the coarse time difference, the smaller the better. In the embodiment, _R1 and _T1 are down-sampled to 5 Hz, and the cross-correlation sequence is solved for the FSAS data _R2 and EPSON data _T2 after the down-sampling, and the corresponding interval when the maximum value of the cross-correlation sequence is taken is obtained to obtain the rough Time difference int_bias. Specific steps are as follows:

S3.1: down-sampling R ₁ and T ₁ to obtain reference data R ₂ after down-sampling and data T ₂ to be matched;

S3.2: Pad zero after the shorter sequence of R ₂ and T ₂ until the two lengths are equal;

S3.3: Calculate the cross-correlation sequence of R ₂ and T ₂ . The general principle formula of cross-correlation sequence calculation is as follows:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; R</span>}}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\underset{\mathrm{<span\; class="notranslate">\; K</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; K</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; \&Greater\; Equal;</span>\mathrm{<span\; class="notranslate">\; 0</span>}\\ {\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; R</span>}}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right.$

Where x and y are sequences for correlation processing, the lengths of x and y are both N, the subscripts K and m are used to represent the ordinal numbers of the values in the sequences x and y, and m ranges from -(N-1) to + (N-1). is a cross-correlation sequence, When m is less than zero, use the formula of m greater than zero to calculate the cross-correlation sequence. When m takes different values from -(N-1) to +(N-1), the above formula can be used to calculate a cross-correlation value Finally, a sequence value of the cross-correlation function with a vector length of 2×N-1 is obtained.

Step 3.4: Take the cross-correlation sequence The corresponding interval when the maximum value is m, and calculate the rough time difference int_bias at this time. In this example, the fourth column of R ₂ and T ₂ after padding at the end and equal length is the sequence x and y that need to be correlated , using the above formula to calculate the cross-correlation sequence value The m corresponding to the maximum value is the staggered interval when the correlation between the two columns of data is maximum. Then the rough time difference can be obtained by using the following formula:

int_bias=R ₂ (t ₀ )-T ₂ (t ₀ )-T ₂ _interval×m

R ₂ (t ₀ ) is the start time of R ₂ , T ₂ (t ₀ ) is the start time of T ₂ , T ₂ _interval is the time length represented by the staggered interval, in this example it is the sampling interval of T ₂ , then T ₂ _interval×m is the time when two columns of data are staggered.

S4: According to the obtained rough time difference int_bias, make the time period of the reference data R ₂ cover the time period of the data to be matched T ₂ , and each exceed a certain length of the data to be matched T ₂ before and after, and intercept the reference data R ₂ and the data to be matched T 2 ₂ . The purpose of this step is to find the corresponding time period between the data to be matched T ₂ and the reference data R ₂ based on the rough time difference, so that the reference data R ₂ contains the data to be matched T ₂ , and each exceeds the data to be matched T ₂ by a certain length According to the principle, the reference data R ₂ and the data to be matched T ₂ are intercepted. In this embodiment, the obtained rough matching time accuracy is at the level of 1 second, so when fetching the FSAS data, the front and rear are respectively exceeded by 0.5 seconds to enter S5.

S5: Unify the sampling rate of R ₂ and T ₂ to obtain the reference data R ₃ and the data to be matched T ₃ with the same sampling rate; the sampling rate of FSAS is 200Hz, the sampling rate of EPSON is 125Hz, and finally perform linear interpolation on T ₀ to make it The sampling rate is 200Hz.

After all the above steps, the FSAS data is R ₃ , and the EPSON data is T ₃ .

S6: Obtain the optimal time difference bias between R ₃ and T ₃ and the scale drift factor k; the calculation of the time difference and time scale drift factor is the core of this algorithm, and it is the most computationally intensive and most complicated step.

S6.1: Calculate the sub-sampling interval T_inter, preferably between 1/10 and 1/2 of the sampling interval. The sub-sampling interval is defined as any value lower than the sampling interval. In the specific implementation, the user's accuracy requirements are taken as the criterion, and the appropriate value is selected. If T_inter is too small, the accuracy is improved, but the accuracy is not improved, and the amount of calculation is increased; if T_inter is too large, it cannot meet the requirement of improving accuracy. For example, in this example, the sampling interval between R ₃ and T ₃ is 0.005 second, and the sub-sampling interval T_inter is defined as 1/5 second of the sampling interval.

S6.2: Using T_inter as the multiple factor, calculate the δs of the translation distance to the right of the data to be matched. In the first cycle, δs is 0 times T_inter, that is, δs=0×T_inter=0.

S6.3: Determine whether δs is a positive integer multiple of the sampling interval. If yes, go to step S6.15; if not, go to step S6.4.

S6.4: Translate the EPSON data by δs. First, the time column of T ₃ is shifted to the right, and the translation distance is δs. After the shift, the time column is used as the interpolation point, and the data of other columns are linearly interpolated to obtain a new value T ₄ . For example, there are original times t ₁ and t ₂ , and the corresponding values are val ₁ and val ₂ , then for the new time t ₃ (t ₁ <t ₃ <t ₂ ), the corresponding value val ₃ is linearly interpolated by the following formula :

${\mathrm{<span\; class="notranslate">\; val</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; val</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span>\frac{{\mathrm{<span\; class="notranslate">\; val</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; val</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}$

S6.5: Take Δt as the window length, divide T ₄ into n limited segments, and take the current segment data T ₄ (i), 1≤i≤n. Within the window length of Δt, it should be able to contain rich dynamic information to restrain the influence of random noise. In this example, Δt is taken as 60 seconds, T ₄ is divided into 61 sections, and the first section T ₄ (1) is taken.

S6.6: Determine whether the data to be matched is extracted to the nth segment. If yes, go to step S6.13. If not, go to step S6.7.

S6.7: Set the sliding window W ₁ , the window length is Δt, and the starting point is the starting time of the current segment. The window length is the same as the window length in S6.5, which is 60 seconds. W ₁ shifts one sampling interval backward along the time axis each time to obtain reference data. The sliding distance of the window is the total time of backward shifting by W ₁ , and the calculation method is: W ₁ Shift backwards by the number of sample intervals multiplied by the sample interval time.

S6.8: Determine whether the sliding distance of the window is greater than the set length σ _s (σ _s depends on the rough matching time precision, which is about twice the rough matching time precision), if yes, go to step S6.11. Otherwise, take W ₁ as the window, and take the reference data R ₃ (W ₁ ). In this example, σ _s is 1 second.

S6.9: Calculate the Pearson correlation coefficient between T ₄ (i) and R ₃ (W ₁ ). The calculation principle of Pearson correlation coefficient is as follows:

$\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; =</span>\frac{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span>}{\sqrt{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&Center\; Dot;</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}}$

in is the mean value of the x series, is the mean value of the y sequence, N is the length of the sequence x and y, and r is the Pearson correlation coefficient between the sequence x and y. In this example, T ₄ (i) and R ₃ (W ₁ ) each contain six columns of data, and the correlation coefficients of each corresponding column of data are calculated independently according to the above formula, and each corresponding axis obtains a correlation coefficient, a total of six.

S6.10: Determine whether the last digit of the reference data has been read, if so, go to step S6.11. Otherwise, after sliding the window W ₁ backward by one sampling interval of the reference data, go to step S6.8.

S6.11: Calculate the time difference Bias(i) and the corresponding time time(i). For this segment of EPSON data, during the sliding process of W ₁ , the corresponding correlation coefficient series can be obtained between each corresponding column data of T ₄ (i) and R ₃ (W ₁ ). For one set of corresponding column data, when the correlation coefficient is the largest, the waveforms of the two columns of data match best. At this time, the maximum correlation coefficient max_corr, the number of sampling intervals W ₁ slides to the right is num. Then this axis can use the following formula to calculate a Bias(i), and the other axes are the same.

Bias(i)=R ₃ (t ₀ )-T ₄ (t ₀ )+R ₃ _interval×num

R ₃ (t ₀ ) is the start time of R ₃ , T ₄ (t ₀ ) is the start time of T ₄ , and R ₃ _interval is the time represented by W ₁ moving to the right by one sampling interval. This example is the sampling interval of R ₃ , Then R ₃ _interval×num is the total duration of W ₁ sliding to the right.

After calculating the Bias(i) of the six axes, take the max_corr of each axis as the weight, and obtain the unique Bias(i) of T ₄ (i) on a weighted average. The max_corr(i) of this section is the arithmetic mean of the max_corr of each axis. time(i) is the middle time of T ₄ (i).

The time difference between the data to be matched and the reference data in this period has been calculated.

S6.12: Extract the next segment of data to be matched. Take the next segment of T ₄ (i) and go to step 6.6.

As shown in Figure 4-1 and 4-2, the maximum correlation coefficient obtained for all segments of the six axes of EPSON. The larger the correlation coefficient, the higher the matching degree of the data waveform. It can also be seen from the figure that after the 16th segment, the correlation coefficient is close to 1, indicating that the waveform consistency of EPSON and FSAS data is high, and the best matching segment can be identified by this method.

As shown in Figures 5-1 and 5-2, for the time difference sequence calculated from the maximum correlation coefficient of each segment, select the period when the correlation coefficient is relatively large. It can be seen from the figure that this method can not only calculate the time difference of EPSON, but also find that the frequency standard of EPSON drifts.

S6.13: Calculate the time difference bias and time scale drift factor k by least squares adjustment. Taking bias and k as the parameters to be estimated, and Bias(i) as the observed value, the observation equation is as follows:

Bias(i)=k×[time(i)-t ₀ ]+bias

The initial value of bias, Bias0, is the median value of the Bias(i) sequence, and the corresponding time is t ₀ , and the initial value of k is 0.

The weight of Bias(i) is max_corr(i). In order to remove gross errors in the observed values, set a correlation coefficient threshold σ _r , in this example, σ _r is 0.9. Only when max_corr(i) is greater than 0.9, the observed value is substituted into the observation equation, and the least squares adjustment is performed to obtain a set of bias and k.

S6.14: Determine whether δs is greater than a sampling interval, if so, go to step S6.15. Otherwise, double δs by T_inter and go to step S6.3.

S6.15: Select the best combination of bias and k. Different sets of bias and k are used to obtain data to be matched after synchronization at different times. Calculate the covariance between the data to be matched and the reference data. If the data is single-axis, select the group of bias and k with the largest covariance as the best combination of bias and k; if the data is multi-axis, select the most credible one-axis data for covariance based on actual experience Solve, choose the best combination of bias and k.

S7: According to the best combination of bias and k, update the time T ₀ of the data to be matched, and store it. In this example, the most credible data is the fourth column. After time synchronization, the third group of data has the largest covariance. Use the following formula to calculate the update T ₀ time:

t'=(tt ₀ )×k+bias+t

Among them, t is the time column of T ₀ , t ₀ is the time corresponding to the initial value of the time difference selected in step 6.12, k is the time scale drift factor, bias is the time difference, and t′ is the time after t is updated.

Replace the time column in T ₀ with t', and re-store the data as a data file in the same format as T ₀ .

S8: End.

As shown in Fig. 6, after the EPSON and FSAS are time-synchronized using the algorithm of the present invention, the waveform diagram is marked at the same time in the vertical direction of the gyro. It can be seen from this figure that this algorithm can correctly perform time synchronization operations.

As shown in FIG. 7 , it is a partial detailed diagram of the waveform matching in a certain period of time in FIG. 6 . It can be seen from the figure that this algorithm can not only correctly synchronize the time between the data to be matched and the reference data, but also have high synchronization accuracy.

A time synchronization system based on waveform matching, characterized in that it includes the following modules:

Reading module: used to read reference data R ₀ and data to be matched T ₀ ;

Data processing format module: used to unify the data processing format of the reference data R ₀ and the data T ₀ to be matched, and extract the converted reference data R ₁ and the data T ₁ to be matched;

Downsampling module: it is used to downsample _R1 and _T1 respectively, solve the cross-correlation sequence for the reference data _R2 after downsampling and the data _T2 to be matched, and obtain the corresponding interval when the maximum value of the cross-correlation sequence is obtained. The coarse time difference int_bias at this time;

Coarse time difference module: used to make the time period of the reference data _R2 cover the time period of the data to be matched _T2 according to the obtained rough time difference int_bias, and each of the front and rear exceeds a certain length of the data to be matched _T2 , intercepting the reference data _R2 and Data to be matched T ₂ ;

Sampling rate calculation module: used to unify the sampling rates of R ₂ and T ₂ to obtain reference data R ₃ and data to be matched T ₃ of the same sampling rate;

Calculation module: used to obtain the optimal time difference bias and time scale drift factor k of R ₃ and T ₃ ;

Update module: used to update the data to be matched T ₀ time according to the best bias and k combination, and store it.

Claims (7)

1. the method for synchronizing time based on Waveform Matching, is characterized in that, comprises the steps:

Step 1: read reference data R ₀and data T to be matched ₀;

Step 2: unified reference data R ₀and data T to be matched ₀data processing form, extract the reference data R after conversion ₁and data T to be matched ₁;

Step 3: respectively to R ₁and T ₁down-sampled, to the reference data R after down-sampled ₂and data T to be matched ₂carry out simple crosscorrelation sequence and solve, corresponding interval while getting simple crosscorrelation sequence maximal value, obtains thick mistiming int_bias now;

Step 4: the thick mistiming int_bias according to trying to achieve, makes reference data R ₂time period cover data T to be matched ₂time period, and front and back respectively exceed data T to be matched ₂certain length, intercepting reference data R ₂and data T to be matched ₂, enter step 5;

Step 5: unified R ₂and T ₂sampling rate, obtains the reference data R of identical sampling rate ₃with data T to be matched ₃;

Step 6: try to achieve R ₃and T ₃the poor bias of optimal time and markers drift parameter k;

Step 7: according to best bias and k combination, upgrade data T to be matched ₀time, and storage;

Step 8: finish.

2. a kind of method for synchronizing time based on Waveform Matching according to claim 1, is characterized in that: unified reference data R in described step 2 ₀and data T to be matched ₀the formula of data processing form as follows:

v＝Δs/Δt _s

Wherein, Δ s is gyro or the measured increment size of accelerometer between former and later two samplings, Δ t _sfor the time interval between former and later two samplings, v is the rate value that Δ s is corresponding;

Reference data R after described conversion ₁and data T to be matched ₁that the time row after changing and an axis data of carrying out Waveform Matching form.

3. a kind of method for synchronizing time based on Waveform Matching according to claim 1, is characterized in that: in described step 3, comprise the steps:

Step 3.1: respectively to R ₁and T ₁down-sampled, obtain the reference data R after down-sampled ₂and data T to be matched ₂;

Step 3.2: to down-sampled rear R ₂and T ₂among shorter sequence trailing zero, until both are equal in length;

Step 3.3: calculate R ₂with T ₂simple crosscorrelation sequence; Simple crosscorrelation sequence computing formula is as follows:

${\hat{R}}_{\mathrm{xy}}\left(m\right)=\left\{\begin{array}{cc}\underset{K=0}{\overset{N-m-1}{\mathrm{\&Sigma;}}}{x}_{K+m}{y}_m& m\&GreaterEqual;0\\ {\hat{R}}_{\mathrm{xy}}(-m)& m<0\end{array}\right.$

Wherein, x and y are the sequence of carrying out correlativity processing, and x and y length are N, and K and m are respectively used to represent the ordinal number of the value in sequence x and y, m value is-(N-1) to+(N-1);

Step 3.4: get simple crosscorrelation sequence corresponding interval during maximal value m, obtains thick mistiming int_bias now, and the thick mistiming can utilize following formula to try to achieve:

int_bias＝R ₂(t ₀)-T ₂(t ₀)-T ₂_interval×m

Wherein, R ₂(t ₀) be R ₂start time, T ₂(t ₀) be T ₂start time, T ₂_ interval is the time span of the interval representative that staggers, is T ₂sampling interval, T ₂_ interval * m is the time that two column datas stagger.

4. a kind of method for synchronizing time based on Waveform Matching according to claim 1, is characterized in that: in described step 6.15, comprise following sub-step:

Step 6.1: calculate sub-sampling T_inter interval time;

Step 6.2: take T_inter as factor, calculate the data to be matched δ s of translation distance to the right;

Step 6.3: judge whether δ s is the positive integer times of sampling interval; If so, forward step S6.15 to, if not, forward step S6.4 to;

Step 6.4: by data data T to be matched ₃translation δ s obtains T ₄;

Step 6.5: take Δ t as length of window, by T ₄be divided into limited n section, get present segment data T ₄(i), 1≤i≤n;

Step 6.6: judge whether data to be matched extract n section; If so, forward step S6.13 to; If not, forward step S6.7 to;

Step 6.7: set moving window W ₁, window is long is Δ t, starting point is the present segment zero hour;

Step 6.8: judge whether window sliding distance is greater than preseting length σ _s; As if so, forwarded step S6.11 to;

Otherwise, with W ₁for window, get reference data R ₃(W ₁);

Step 6.9: calculate T ₄and R (i) ₃(W ₁) between Pearson's related coefficient; Pearson's Calculation of correlation factor as shown in the formula:

$r=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_i-\stackrel{\&OverBar;}{x})(y_i-\stackrel{\&OverBar;}{y})}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(x_i-\stackrel{\&OverBar;}{x})}^2\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2}}$

Wherein for x serial mean, for y serial mean, N is the length of sequence x and y, and r is the Pearson's related coefficient between sequence x and y;

Step 6.10: judge that whether reference data reads last position, if so, forwards step S6.11 to; Otherwise,

By window W ₁slide backward after the sampling interval of a reference data, forward step S6.8 to;

Step 6.11: calculate this section of mistiming Bias (i) and corresponding time time (i);

Step 6.12: extract next section of data to be matched.Take off one section of T ₄(i), forward step 6.6 to;

Step 6.13: with the poor bias of least square adjustment seeking time and markers drift parameter k; Using bias and k as solve for parameter, and Bias (i) is observed reading, and observation equation is as follows:

Bias(i)＝k×[time(i)-t ₀]+bias

Wherein the initial value Bias0 of bias is Bias (i) sequence intermediate value, and the corresponding time is t ₀, the initial value of k is 0;

Step 6.14: judge that whether δ s is greater than a sampling interval, if so, forwards step S6.15 to; Otherwise, by the δ s T_inter that doubles, forward step S6.3 to;

Step 6.15: select best bias and k combination; Utilize not bias and k on the same group, obtain the to be matched data of different time after synchronous.

5. a kind of method for synchronizing time based on Waveform Matching according to claim 4, is characterized in that: in described step 6.15, if data are single shaft, that group bias while selecting covariance maximum and k are as best bias and k combination; If data are multiaxis, select the most believable axis data to carry out covariance and solve, select best bias and k combination.

6. a kind of method for synchronizing time based on Waveform Matching according to claim 1, is characterized in that: in described step 7, utilize following formula to calculate and upgrade T ₀time:

t′＝(t-t ₀)×k+bias+t

Wherein, t is T ₀time row, t ₀for time corresponding to mistiming initial value of selecting in step 6.12, k is markers drift parameter, and bias is the mistiming, and t ' is the rear time for t upgrades.

7. the clock synchronization system based on Waveform Matching, is characterized in that, comprises as lower module:

Read module: for reading reference data R ₀and data T to be matched ₀;

Data processing format module: for unified reference data R ₀and data T to be matched ₀data processing form, extract the reference data R after conversion ₁and data T to be matched ₁;

Down-sampled module: for respectively to R ₁and T ₁down-sampled, to the reference data R after down-sampled ₂and data T to be matched ₂carry out simple crosscorrelation sequence and solve, corresponding interval while getting simple crosscorrelation sequence maximal value, obtains thick mistiming int_bias now;

Interception module: for according to the thick mistiming int_bias trying to achieve, make reference data R ₂time period cover data T to be matched ₂time period, and front and back respectively exceed data T to be matched ₂certain length, intercepting reference data R ₂and data T to be matched ₂;

Sampling rate computing module: for unified R ₂and T ₂sampling rate, obtains the reference data R of identical sampling rate ₃with data T to be matched ₃;

Computing module: for trying to achieve R ₃and T ₃the poor bias of optimal time and markers drift parameter k;

Update module: for according to best bias and k combination, upgrade data T to be matched ₀time, and storage.