A high-precision real-time single-point speed measurement method, computer equipment and storage medium based on high-frequency filtering fitting

Abstract

The present application relates to a high-precision real-time single-point speed measurement method, computer equipment and storage medium based on high-frequency filtering fitting, which includes obtaining high-frequency GNSS data and performing data preprocessing, detecting cycle slips of high-frequency GNSS data, eliminating satellites and their carrier phase observations that have cycle slips between adjacent epochs, constructing a Kalman_TDCP speed measurement model, performing Kalman_TDCP high-frequency speed measurement on the high-frequency GNSS data that has undergone data preprocessing, and based on the Kalman_TDCP speed measurement model, after obtaining the speed measurement results, storing the associated information of the residuals after the speed measurement results are verified; and performing polynomial fitting using the determined best fitting order and the best fitting time window. The present application has the effect of improving the speed measurement accuracy and real-time performance of high-precision single-point speed measurement and expanding its scope of application in different scenarios.

Description

Technical Field

[0001] The present application relates to the technical field of high-precision real-time speed measurement, and in particular to a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting, a computer device and a storage medium. Background Art

[0002] In recent years, the demand for high-precision real-time dynamic speed measurement has grown rapidly. Due to the characteristics of continuity, all-weather and global coverage of the Global Navigation Satellite System, GNSS, it can not only be used for navigation, positioning and timing, but also provides new technologies and methods for high-precision speed measurement, including differential speed measurement and single-point speed measurement.

[0003] There are two common single-point velocity measurement methods: raw Doppler observation velocity measurement and inter-epoch carrier differential velocity measurement (TDCP).

[0004] The speed measurement based on raw Doppler observations has good real-time performance, but its speed measurement accuracy is affected by the accuracy of the raw Doppler observations. The accuracy of the raw Doppler observations depends on the receiver hardware and cannot be improved by algorithms. It is also easily affected by complex environments and interferences, resulting in gross errors that cause the speed measurement to deviate seriously from the true value.

[0005] When using the TDCP method for dynamic speed measurement, there is an unavoidable defect. As the speed of the moving carrier changes dynamically, since its estimated speed is approximately equal to the average speed between the previous and next epochs, it will be different from the instantaneous speed of the moving carrier at the current moment, that is, there is a speed measurement delay, especially under dynamic conditions such as turning and rapid speed changes, the speed measurement delay is particularly significant.

[0006] Therefore, the current single-point speed measurement method has problems of low test accuracy and speed measurement delay, and its application scope is narrow, which needs to be improved. Summary of the invention

[0007] In order to improve the speed measurement accuracy and real-time performance of high-precision single-point speed measurement and expand its application scope, the present application provides a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting.

[0008] The above-mentioned invention objective of the present application is achieved through the following technical solutions:

[0009] A high-precision real-time single-point speed measurement method based on high-frequency filtering fitting comprises the following steps:

[0010] Acquire high-frequency GNSS data and perform data preprocessing; data preprocessing includes cycle slip detection of high-frequency GNSS data, eliminating satellites and their carrier phase observations that have cycle slips between adjacent epochs; then select a common satellite set between epochs to obtain a single-difference observation equation set of carrier phase between epochs that meets preset conditions; perform pre-verification gross error elimination based on the carrier-to-noise ratio (SNR) information of the observed quantity and the velocity prediction information of the previous beat;

[0011] Construct a Kalman_TDCP speed measurement model, perform Kalman_TDCP high-frequency speed measurement on the high-frequency GNSS data that has undergone data preprocessing, and obtain the average speed v between adjacent epochs based on the Kalman_TDCP speed measurement model, and use it as the instantaneous speed at the middle moment of the previous and next epochs, and store the instantaneous speed and corresponding time series information; after obtaining the speed measurement result, store the associated information of the post-verification residual of the speed measurement result at the same time, and the associated information includes the square error, square error and SSE of the speed measurement result, satellite geometry factor DOP, and reference observation quantity;

[0012] Performing scene adaptive fitting on the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results, including: calculating the dynamic state parameters of the current carrier motion based on the high-frequency speed measurement results, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; eliminating and weighting each speed measurement result through associated information to obtain the determined optimal fitting order and optimal fitting time window;

[0013] The polynomial fitting is performed using the determined best fitting order and best fitting time window.

[0014] In a preferred example of the present application, high-frequency GNSS data is acquired and data preprocessing is performed; the data preprocessing includes cycle slip detection of the high-frequency GNSS data, and elimination of satellites and carrier phase observations that have cycle slips between adjacent epochs; then a common satellite set between epochs is selected to obtain a carrier phase single-difference observation equation group between epochs that meets preset conditions; and the steps of pre-verification gross error elimination based on the carrier-to-noise ratio SNR information of the observed quantity and the speed prediction information of the previous beat include the following steps:

[0015] Receiver r receives multi-frequency signals from n satellites. The carrier phase observation equation received by receiver r at time k and k+1 respectively on the corresponding frequency of the satellite of system S is expressed as equation (1), where system S includes GPS, Galileo or BDS system;

[0016] (1)

[0017] In formula (1), the receiver r and the satellite In the ECEF framework, the positions at time k and k+1 are , , and , The code for the observation satellite. for One of them;

[0018] The subscripts k and k+1 of each variable are the corresponding moments of each variable; For observation satellite The corresponding frequency The carrier phase observation value on , and the unit is m; From the receiver r to the observation satellite The Euclidean distance of is the speed of light; is the mobile station receiver clock error; For observation satellite The clock difference is the tropospheric delay in the phase observation at receiver r; For mobile stations Phase observation corresponding to frequency ionospheric delay; For frequency The corresponding wavelength; For receiver The undifferenced integer ambiguity of The receiver r is the frequency Carrier phase of the signal Carrier phase hardware delay; For observation satellite Frequency Carrier phase of the signal Carrier phase hardware delay; For noise.

[0019] In a preferred example of the present application, high-frequency GNSS data is obtained and data preprocessing is performed; the data preprocessing includes cycle slip detection of the high-frequency GNSS data, and elimination of satellites and carrier phase observations that have cycle slips between adjacent epochs; then a common satellite set between epochs is selected to obtain a carrier phase single-difference observation equation group between epochs that meets preset conditions; and a step of eliminating gross errors before verification based on the carrier-to-noise ratio SNR information of the observed quantity and the speed prediction information of the previous beat is also included:

[0020] The carrier phase observation value of a single GNSS receiver r is differentiated between epochs, and the change in the position of the receiver r between adjacent epochs and the change in the clock error of the receiver r between adjacent epochs are used as the parameters to be estimated. The inter-epoch differential carrier phase observation equation is established to solve the velocity, including:

[0021] The corresponding frequency of satellite p of a certain system S at receiver r at time k+1 The inter-epoch differential carrier phase observation equation on can be expressed as:

[0022] (2)

[0023] In formula (2), the subscript is the time k+1 corresponding to each variable, and ; For observation satellite frequency Single-difference observations between epochs on ; is the difference in the mobile station receiver clock between epochs k and k+1; For observation satellite The difference in satellite clock errors between epochs k and k+1; For the mobile station r to the satellite The difference in station-satellite distance between epochs k and k+1; Observation satellite for mobile station r The non-differenced integer ambiguity change between epochs k and k+1. If the carrier phase observation does not experience a cycle slip, this value is 0. For noise.

[0024] In a preferred example of the present application: after performing inter-epoch differentials on the carrier phase observation values ​​of a single GNSS receiver r, taking the change in the position of the receiver r between adjacent epochs and the change in the clock error of the receiver r between adjacent epochs as the parameters to be estimated, and establishing the inter-epoch differential carrier phase observation equation for velocity calculation, the following steps are further performed: If the carrier phase observation value without cycle slip is selected, after linearizing equation (2), it can be obtained:

[0025] (3)

[0026] In formula (3), represents the direction cosine matrix from receiver r to satellite p at time k+1, is the parameter to be estimated at time k+1, which is related to the speed There is a definite conversion formula:

[0027] (4)

[0028] If at time k+1, the receiver r receives m common carrier observations without cycle slips between epochs and n common satellites, then an m-column double-difference observation equation system between epochs can be formed. When both m and n are greater than or equal to 4, the solution can be obtained by Perform velocity calculation.

[0029] In a preferred example of the present application: the step of constructing a Kalman_TDCP speed measurement model, performing Kalman_TDCP high-frequency speed measurement on the high-frequency GNSS data after data preprocessing, obtaining the average speed v between adjacent epochs based on the Kalman_TDCP speed measurement model, and using it as the instantaneous speed at the middle moment of the previous and next epochs, storing the instantaneous speed and the corresponding timing information; after obtaining the speed measurement result, storing the associated information of the post-verification residual of the speed measurement result at the same time, the associated information including the square error, square error and SSE of the speed measurement result, satellite geometry factor DOP and reference observation quantity, including the steps of:

[0030] Based on the Kalman_TDCP model, the state equation and observation equation of the discrete-time Kalman filter at time k can be expressed as:

[0031] (5.1)

[0032] In formula (5.1), Represents the position change information at time k and the receiver clock error epoch difference The parameter vector to be estimated is composed of Predicted value, is the state transfer matrix from time k-1 to time k. High-frequency TDCP is used for velocity measurement. The position change between epochs of the model at time t-1 is equal to the position change between the epochs before and after time t. Therefore, in this model ;

[0033] are system process noise and observation noise respectively, which obey Gaussian distribution with mean 0; is the observation vector at time k, corresponding to the sum of the difference between the epoch carrier phase observations and the difference between the epoch satellite clock errors; H is the observation matrix, which can be constructed based on formula (3); E is the fourth-order unit matrix.

[0034] In a preferred example of the present application: the Kalman_TDCP speed measurement model is constructed, the high-frequency GNSS data that has undergone data preprocessing is subjected to Kalman_TDCP high-frequency speed measurement after Kalman filtering, and based on the Kalman_TDCP speed measurement model, the average speed v between adjacent epochs is obtained, and it is used as the instantaneous speed at the middle moment of the previous and next epochs, and the instantaneous speed and the corresponding timing information are stored; after obtaining the speed measurement result, the associated information of the post-verification residual of the speed measurement result is stored at the same time, and the associated information includes the error square, square error and SSE of the speed measurement result, the satellite geometry factor DOP and the reference observation quantity, and also includes the steps of:

[0035] After obtaining the state prior estimate through time update by formula (5.1), it is necessary to use the actual measurement value to correct the state prior estimate obtained by formula (5.1) to obtain the optimal estimate of the system state: the calculation formula is as follows:

[0036] (5.2)

[0037] In formula (5.2), the first two steps are generally called time updates, and the last three steps are measurement updates; is the covariance matrix of the system process noise; is the covariance matrix of the observation noise; is the prior estimate of the state The weight matrix of is the coefficient matrix, called the Kalman filter gain; is the optimal state estimate and its mean square error matrix; is the fourth-order identity matrix;

[0038] When the optimal state estimate is obtained After that, the Kalman_TDCP speed measurement result can be obtained through formula (4).

[0039] In a preferred example of the present application: the scene adaptive fitting of the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results includes: based on the high-frequency speed measurement results, calculating the dynamic state parameters of the current carrier motion, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; through the associated information, eliminating and weighting each speed measurement result to obtain the determined optimal fitting order and optimal fitting time window, including the steps of:

[0040] Through high-frequency Kalman_TDCP speed measurement, obtain arrive m speed measurement results and their associated information at the time, and m calculation samples are obtained ;

[0041] Assume that the model estimate is , then the calculation formulas for the total sum of squares SST, error sum of squares SSE and regression sum of squares SSR of the calculated sample are as follows:

[0042] (6)

[0043] Its coefficient of determination The calculation formula is as follows:

[0044] (7)

[0045] The value range is between 0 and 1;

[0046] Coefficient of determination Set a confidence threshold using the sum of squared errors (SSE) and , when the fitting order is iteratively fitted from low to high order, or the fitting time window is iteratively fitted from short to long, if or , the fitting order is considered to be optimal. If it is always greater than the confidence threshold, the minimum value of the coefficient of determination is selected. as the best fit order.

[0047] In a preferred example of the present application: the scene adaptive fitting of the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results includes: based on the high-frequency speed measurement results, calculating the dynamic state parameters of the current carrier movement, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; through the associated information, eliminating and weighting each speed measurement result to obtain the determined optimal fitting order and optimal fitting time window, and also includes the steps of:

[0048] Set the minimum fitting time window Restrictions to improve polynomial fitting stability;

[0049] A polynomial fitting of greater than 5th order is performed on the total samples to obtain the dynamic state parameters of the carrier motion. The dynamic state parameters include velocity range, acceleration and acceleration information, and the upper and lower limits of the fitting order are determined.

[0050] The above-mentioned second objective of the present application is achieved through the following technical solutions:

[0051] A computer device comprises a memory, a processor and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of the high-precision real-time single-point speed measurement method based on high-frequency filtering fitting are implemented.

[0052] The third objective of the present application is achieved through the following technical solutions:

[0053] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of the high-precision real-time single-point speed measurement method based on high-frequency filtering fitting.

[0054] In summary, the present application includes at least one of the following beneficial technical effects:

[0055] 1. Based on the TDCP model, comprehensive high-frequency solution is used, and the frequency of 10 Hz is generally selected. By utilizing the time and space correlation of observation errors, the carrier phase epoch difference can eliminate most observation errors such as satellite orbit error, carrier phase hardware delay, ionosphere and troposphere delay, etc., which greatly reduces the acquisition error of high-frequency GNSS data. In subsequent calculations, the residual part of these observation errors has a negligible impact on speed measurement. Furthermore, Kalman_TDCP high-frequency speed measurement and Kalman filtering smoothing can greatly reduce speed measurement delay, and filtering smoothing improves speed measurement accuracy and stability.

[0056] The high-frequency speed measurement results are adaptively fitted to obtain high-precision speed measurement results without delay. The polynomial fitting order and window are adaptively adjusted to perform optimal speed measurement fitting, so as to reduce speed measurement noise, improve speed measurement accuracy, eliminate speed measurement delay, and further obtain high-precision real-time speed measurement results. Only one receiver is needed, and the dynamic speed measurement results can meet both high-precision and real-time characteristics without the need for a base station. This provides a reliable solution for GNSS products that have high-precision real-time speed measurement requirements and have application scenarios with or without base station coverage.

[0057] 2. Through formula (5.1), the estimated value of the system state parameter to be estimated has the minimum mean square error (MMSE).

[0058] 3. Through high-frequency Kalman_TDCP speed measurement, m speed measurement results and their corresponding reliability indicators from t1 to t2 are obtained. The scene adaptive fitting model first calculates and configures the scene parameters based on these m sample data, and then performs weighted polynomial fitting after setting the scene parameters to obtain high-precision real-time speed results for the current epoch.

[0059] 4. In order to obtain the optimal speed measurement result of the current epoch, it is necessary to select an optimal polynomial model to approximate the data points. Too high a polynomial order may lead to overfitting, while too low an order may lead to underfitting. Similarly, too long a time window fitting length may make it difficult for the polynomial model to obtain a satisfactory fitting result, while too short a time window fitting length may cause the fitting result to be affected by measurement noise and gross errors, resulting in fitting deviation true value anomalies. Therefore, before fitting, it is necessary to select an appropriate polynomial fitting order and a suitable time window fitting length to ensure that the subsequent weighted polynomial fitting can solve the optimal fitting speed measurement result.

[0060] 5. The key point of scene parameter calculation and configuration is to determine the optimal fitting order and optimal fitting time window length of the current motion state of the receiver. Goodness of fit refers to the degree of fit of the regression line to the observed value. The parameter statistic that measures the goodness of fit in the fitting model is the goodness of fit determination coefficient. (Coefficient of Determination) and SSE (Residual Sum of Squares).

[0061] 6. Since the speed data of the fitting model is obtained based on high-frequency Kalman filter speed measurement, its noise is larger than the low-frequency TDCP speed measurement result, so it is necessary to set the minimum fitting time window Restrictions are set to ensure the stability of polynomial fitting. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 A schematic diagram of an embodiment of a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting in the present application;

[0063] Figure 2 The figure shows the eastward velocity measurement result and the true value result obtained by the instantaneous Doppler velocity measurement method;

[0064] Figure 3 The figure shows the velocity measurement result and the true value result of the north direction obtained by the instantaneous Doppler velocity measurement method;

[0065] Figure 4 The vertical velocity measurement result and the true value result diagram obtained by the instantaneous Doppler velocity measurement method;

[0066] Figure 5 The speed measurement result and the true value result diagram of the east direction obtained by the low-frequency TDCP speed measurement method;

[0067] Figure 6 The speed measurement result and the true value result diagram of the north direction obtained by the low-frequency TDCP speed measurement method;

[0068] Figure 7 The vertical velocity measurement result and the true value result diagram obtained by the low-frequency TDCP velocity measurement method;

[0069] Figure 8 This is a diagram of the eastward speed measurement result and the true value result obtained by the high-precision real-time speed measurement method based on high-frequency fitting in this application;

[0070] Fig. 9 This is a diagram of the north-bound speed measurement result and the true value result obtained by the high-precision real-time speed measurement method based on high-frequency fitting in this application;

[0071] Fig.10 This is a diagram of the vertical speed measurement results and true value results obtained by the high-precision real-time speed measurement method based on high-frequency fitting in this application. DETAILED DESCRIPTION

[0072] The following is combined with Figure 1-Figure 10 This application is described in further detail.

[0073] As one of the important parameters of a moving carrier, speed has important applications in vehicle navigation and positioning, drone autopilot, inertial navigation system calibration, precision agriculture seeding and fertilization, satellite orbit docking, and aerial gravity measurement. In recent years, the demand for high-precision real-time dynamic speed measurement has grown rapidly. Due to the characteristics of GNSS such as continuity, all-weather and global coverage, it can not only be used for navigation, positioning and timing, but also provides new technologies and methods for high-precision speed measurement. At present, using satellites to measure speed is one of the common methods and important technical means in the field of speed measurement.

[0074] Depending on the operating method and cost control, GNSS speed measurement can be divided into differential speed measurement and single-point speed measurement.

[0075] GNSS differential speed measurement is based on GNSS high-precision positioning results, and obtains speed through position differential. If this speed measurement method is to ensure the high accuracy of the speed measurement results, it needs to rely on the base station or the virtual reference station data on the network CORS side, and is limited by the distance between the base station and the mobile station. Its working range is quite limited. When the mobile carrier is in an environment without base station coverage or enters a complex environment and cannot stably maintain a high-precision RTK fixed solution, the differential speed measurement solution will have applicability issues. GNSS single-point speed measurement relies solely on the receiver on the mobile carrier to obtain the corresponding speed information. Common single-point speed measurement methods mainly include raw Doppler observation value speed measurement and inter-epoch carrier differential speed measurement TDCP.

[0076] The speed measurement based on raw Doppler observations has good real-time performance, but its speed measurement accuracy is affected by the accuracy of the raw Doppler observations. The accuracy of the raw Doppler observations depends on the receiver hardware and cannot be improved by algorithms. It is also easily affected by complex environments and interferences, resulting in gross errors that cause the speed measurement to deviate seriously from the true value.

[0077] TDCP is a method for accurately estimating the receiver position in adjacent epochs by differentiating the carrier phase positioning observation equations of the same visible satellite between adjacent epochs, selecting satellites that do not have cycle slips between epochs to construct the equation group, and then using this to accurately estimate the receiver position in adjacent epochs. The amount of change between , and then solve the average carrier velocity between the adjacent epochs The accuracy of its position change The accuracy can reach mm / s level under static conditions and cm / s level under dynamic open conditions.

[0078] However, when using the TDCP method for dynamic speed measurement, there is an unavoidable defect. As the speed of the moving carrier changes dynamically, the estimated speed is approximately equal to the average speed between the previous and next epochs, so it will be different from the instantaneous speed of the moving carrier at the current moment, that is, there is a speed measurement delay, especially under dynamic conditions such as turning and fast speed changes. The speed measurement delay is particularly significant. The larger the value, the greater the speed measurement delay, which means that when performing speed time alignment, its real-time and matching performance will be worse, and The larger the value, the higher the probability of cycle slips, the fewer the number of cycle slip-free observation equations of the common satellites in adjacent epochs, and the worse the velocity measurement stability may be. The shorter the time, the better the real-time performance, but the average speed The error noise will also increase accordingly, resulting in poor speed measurement accuracy. The average delay of 1 Hz TDCP speed measurement is about 0.5s, and the delay of 10 Hz TDCP speed measurement is about 0.05s, but at the same time, the noise of 10 Hz speed measurement is 10 times that of 1 Hz speed measurement. In addition, in complex environments with frequent cycle jumps and few satellites, the speed measurement accuracy is inevitably affected.

[0079] In one embodiment, referring to Figure 1 The present application discloses a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting, which specifically includes the following steps:

[0080] S10: Acquire high-frequency GNSS data and perform data preprocessing; data preprocessing includes cycle slip detection of high-frequency GNSS data, eliminating satellites and their carrier phase observations that have cycle slips between adjacent epochs; then select a common satellite set between epochs to obtain a carrier phase single-difference observation equation group between epochs that meets preset conditions; and perform pre-verification gross error elimination based on the carrier-to-noise ratio (SNR) information of the observed quantity and the velocity prediction information of the previous beat;

[0081] In this embodiment, if it is the first beat, no pre-verification elimination is performed based on the speed prediction information.

[0082] Step S10 includes the following steps:

[0083] S11: Considering the error in the process of satellite signal broadcasting to the receiver, receiver r receives multi-frequency signals from n satellites. Then the carrier phase observation equation received by receiver r at time k and k+1 respectively on the corresponding frequency of the satellite of system S is expressed as equation (1), where system S includes GPS, Galileo or BDS system;

[0084] (1)

[0085] In formula (1), the receiver r and the satellite In the ECEF framework, the positions at time k and k+1 are , , and , The code for the observation satellite. for One of them;

[0086] The subscripts k and k+1 of each variable are the corresponding moments of each variable; For observation satellite The corresponding frequency The carrier phase observation value on , and the unit is m; From the receiver r to the observation satellite The Euclidean distance of is the speed of light; is the mobile station receiver clock error; For observation satellite The clock difference is the tropospheric delay in the phase observation at receiver r; For mobile stations Phase observation corresponding to frequency ionospheric delay; For frequency The corresponding wavelength; For receiver The undifferenced integer ambiguity of The receiver r is the frequency Carrier phase of the signal Carrier phase hardware delay; For observation satellite Frequency Carrier phase of the signal Carrier phase hardware delay; For noise.

[0087] S12: Perform inter-epoch difference on the carrier phase observation value of a single GNSS receiver r, take the change of the position of the inter-epoch receiver r between adjacent epochs and the change of the clock error of the inter-epoch receiver r between adjacent epochs as the parameters to be estimated, and establish the inter-epoch differential carrier phase observation equation for velocity solution, including:

[0088] The corresponding frequency of satellite p of a certain system S at receiver r at time k+1 The inter-epoch differential carrier phase observation equation on can be expressed as:

[0089] (2)

[0090] In formula (2), the subscripts of the variables are This means that each variable corresponds to time k+1, and ; The calculation formula is ; For observation satellite frequency Single-difference observations between epochs on ; is the difference in the mobile station receiver clock between epochs k and k+1; For observation satellite The difference in satellite clock errors between epochs k and k+1 can be directly calculated from the ephemeris parameters; For the mobile station r to the satellite The difference in station-satellite distance between epochs k and k+1; Observation satellite for mobile station r The change in the non-differenced integer ambiguity between epochs k and k+1 is also called the epoch single-differenced integer ambiguity. If the carrier phase observation does not experience a cycle slip, the value is 0. For noise.

[0091] S13: If the carrier phase observation without cycle slip is selected, the linearization of equation (2) yields:

[0092] (3)

[0093] In formula (3), represents the direction cosine matrix from receiver r to satellite p at time k+1, is the parameter to be estimated at time k+1, which is related to the speed There is a definite conversion formula:

[0094] (4)

[0095] If at time k+1, the receiver r receives m common carrier observations without cycle slips between epochs and n common satellites, then an m-column double-difference observation equation system between epochs can be formed. When both m and n are greater than or equal to 4, the solution can be obtained by Perform velocity calculation.

[0096] In this embodiment, by using the temporal and spatial correlation of observation errors, carrier phase epoch-to-epoch differential can eliminate most observation errors such as satellite orbit error, carrier phase hardware delay, ionosphere and troposphere delay, etc. In subsequent calculations, the residual part of these observation errors has a negligible impact on speed measurement. Therefore, by performing inter-epoch differential on the carrier phase observation value of a single GNSS receiver, the receiver position between epochs is calculated based on the position of the receiver in adjacent epochs. The amount of change between and the receiver clock error in adjacent epochs The amount of change between As the parameter to be estimated, the differential carrier phase observation equation between epochs can be established to solve the velocity.

[0097] S20: constructing a Kalman_TDCP speed measurement model, performing Kalman_TDCP high-frequency speed measurement on the high-frequency GNSS data that has undergone data preprocessing, and obtaining the average speed v between adjacent epochs based on the Kalman_TDCP speed measurement model, and using it as the instantaneous speed at the middle moment of the previous and next epochs, storing the instantaneous speed and the corresponding time series information; after obtaining the speed measurement result, storing the associated information of the post-verification residual of the speed measurement result at the same time, the associated information includes the square error, square error and SSE of the speed measurement result, the satellite geometry factor DOP and the reference observation quantity;

[0098] In this embodiment, the storage of the associated information provides a reference for subsequent scene adaptive fitting model weighting.

[0099] Step S20 includes the following steps:

[0100] S21: Based on the Kalman_TDCP model, the state equation and observation equation of the discrete-time Kalman filter at time k can be expressed as:

[0101] ; (5.1)

[0102] In formula (5.1), Represents the position change information at time k and the receiver clock error epoch difference The parameter vector to be estimated is composed of Predicted value, is the state transfer matrix from time k-1 to time k. High-frequency TDCP is used for velocity measurement. The position change between epochs of the model at time t-1 is equal to the position change between the epochs before and after time t. Therefore, in this model ;

[0103] are system process noise and observation noise respectively, which obey Gaussian distribution with mean 0; is the observation vector at time k, corresponding to the sum of the difference between the epoch carrier phase observations and the difference between the epoch satellite clock errors; H is the observation matrix, which can be constructed based on formula (3); E is the fourth-order unit matrix.

[0104] S22: After obtaining the state prior estimation value through time update by formula (5.1), it is necessary to use the actual measurement value to correct the state prior estimation value predicted by formula (5.1) to obtain the optimal estimation result of the system state: the calculation formula is as follows:

[0105] (5.2)

[0106] In formula (5.2), the first two steps are generally called time updates, and the last three steps are measurement updates; is the covariance matrix of the system process noise; is the covariance matrix of the observation noise; is the prior estimate of the state The weight matrix of is the coefficient matrix, called the Kalman filter gain; is the optimal state estimate and its mean square error matrix; is the fourth-order identity matrix;

[0107] When the optimal state estimate is obtained After that, the Kalman_TDCP speed measurement result can be obtained through formula (4).

[0108] S30: performing scene adaptive fitting on the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results, including: calculating the dynamic state parameters of the current carrier motion based on the high-frequency speed measurement results, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; eliminating and weighting each speed measurement result through associated information to obtain the determined optimal fitting order and optimal fitting time window;

[0109] In this embodiment, an adaptive fitting model is constructed to perform scene adaptive fitting on high-frequency speed measurement results. The calculation criteria include determining the upper and lower limits of the polynomial fitting order; and then performing polynomial fitting regression iterative calculations with orders from low to high and time windows from short to long.

[0110] Step S30 includes the following steps:

[0111] S31: Through high-frequency Kalman_TDCP speed measurement, obtain arrive m speed measurement results and their associated information at the time, and m calculation samples are obtained ;

[0112] Assume that the model estimate is , then the calculation formulas for the total sum of squares SST, error sum of squares SSE and regression sum of squares SSR of the calculated sample are as follows:

[0113] (6)

[0114] Its coefficient of determination The calculation formula is as follows:

[0115] (7)

[0116] The value range is between 0 and 1;

[0117] S32: Coefficient of determination Set a confidence threshold using the sum of squared errors (SSE) and , when the fitting order is iteratively fitted from low to high order, or the fitting time window is iteratively fitted from short to long, if or , the fitting order is considered to be optimal. If it is always greater than the confidence threshold, the minimum value of the coefficient of determination is selected. As the best fit order, that is, the calculation criterion.

[0118] S33: Set the minimum fitting time window Restrictions to improve polynomial fitting stability;

[0119] S34: Perform a polynomial fitting of greater than 5th order on the total samples to obtain the dynamic state parameters of the carrier motion, the dynamic state parameters including velocity range, acceleration and acceleration information, and determine the upper and lower limits of the fitting order.

[0120] In this embodiment, The value range is between 0 and 1. The closer it is to 1, the closer the sum of squared errors is to 0. The smaller the error, the more accurate the fitting. Similarly, the smaller the SSE is. The closer it is to 0, the more accurate the relative fitting is.

[0121] S40: Perform polynomial fitting using the determined optimal fitting order and optimal fitting time window.

[0122] As a comparative example, a set of dual-frequency multi-system vehicle-mounted urban environment GNSS test data with combined inertial navigation was used for experiments, and speed measurement was performed using instantaneous Doppler speed measurement, low-frequency TDCP speed measurement and the high-precision real-time speed measurement method based on high-frequency filtering fitting of the present application. The speed information output by the combined inertial device during the test was used as the true value for accuracy analysis and statistics.

[0123] Method 1 is an instantaneous Doppler velocity measurement method, method 2 is a low-frequency TDCP velocity measurement method, and method 3 is a high-precision real-time velocity measurement method based on high-frequency fitting in this application.

[0124] Reference Figure 2-Figure 10 , Figure 2-4 The velocity measurement results and true value results of the east, north and vertical directions obtained by the instantaneous Doppler velocity measurement method of method 1 are shown in Figure 2. Figure 5-7 The speed measurement results and true value results of the east, north and vertical directions obtained by the low-frequency TDCP speed measurement method of method 2 are shown in Figure 2. Figure 8-10The speed measurement results and true value results of the east, north and vertical directions obtained by the high-precision real-time speed measurement method based on high-frequency fitting in method 3 are shown in Figure 1. Figure 2-Figure 10 The GPS week second is a way of expressing time in the Global Positioning System (GPS), Ve represents the eastward speed, Vn represents the northward speed, Vu represents the vertical speed, and Method 1, Method 2, and Method 3 are respectively the speed measurements of the instantaneous Doppler speed measurement, the low-frequency TDCP speed measurement, and the high-precision real-time speed measurement method based on high-frequency filtering fitting of the present application, and the true value is the speed measurement output by the fiber optic inertial navigation.

[0125] The accuracy statistics of the above three methods are performed to obtain the following accuracy statistics result table:

[0126]

[0127] Depend on Figure 2-Figure 10 As can be seen from the table above, the instantaneous Doppler speed measurement method of method 1 has more gross errors and poor speed measurement accuracy; the low-frequency TDCP speed measurement method of method 2 calculates the average speed between epochs, which is affected by its speed measurement delay, and its speed measurement accuracy is also poor, but its gross error is significantly less than that of method 1; the high-precision real-time speed measurement method of method 3 based on high-frequency fitting has a speed measurement accuracy significantly higher than the other two methods, with high accuracy and good real-time performance.

[0128] It should be understood that the size of the serial numbers of the steps in the above embodiments does not mean the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present application.

[0129] In one embodiment, a computer device is provided, which may be a server. The computer device includes a processor, a memory, a network interface, and a database connected via a system bus. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program, and a database. The network interface of the computer device is used to communicate with an external terminal via a network connection. When the computer program is executed by the processor, a high-precision real-time single-point speed measurement method based on high-frequency filter fitting is implemented.

[0130] In one embodiment, a computer device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting are implemented.

[0131] In one embodiment, a computer-readable storage medium is provided, on which a computer program is stored. When the computer program is executed by a processor, the steps of a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting are implemented:

[0132] Those of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, storage, database or other media used in the embodiments provided in this application may include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM) or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in many forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link (Synchlink) DRAM (SLDRAM), memory bus (Rambus) direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0133] Those skilled in the art can clearly understand that for the convenience and simplicity of description, only the division of the above-mentioned functional units and modules is used as an example. In actual applications, the above-mentioned functions can be distributed and completed by different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above.

[0134] The embodiments described above are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features may be replaced by equivalents. Such modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application, and should all be included in the protection scope of the present application.

Claims

1. A high-precision real-time single-point speed measurement method based on high-frequency filtering fitting, characterized in that: Includes steps: Acquire high-frequency GNSS data and perform data preprocessing; data preprocessing includes cycle slip detection of high-frequency GNSS data, eliminating satellites and their carrier phase observations that have cycle slips between adjacent epochs; then select a common satellite set between epochs to obtain a single-difference observation equation set of carrier phase between epochs that meets preset conditions; perform pre-verification gross error elimination based on the carrier-to-noise ratio (SNR) information of the observed quantity and the velocity prediction information of the previous beat; Construct a Kalman_TDCP speed measurement model, perform Kalman_TDCP high-frequency speed measurement on the high-frequency GNSS data that has undergone data preprocessing, and obtain the average speed v between adjacent epochs based on the Kalman_TDCP speed measurement model, and use it as the instantaneous speed at the middle moment of the previous and next epochs, and store the instantaneous speed and corresponding time series information; after obtaining the speed measurement result, store the associated information of the post-verification residual of the speed measurement result at the same time, and the associated information includes the square error, square error and SSE of the speed measurement result, satellite geometry factor DOP, and reference observation quantity; Performing scene adaptive fitting on the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results, including: calculating the dynamic state parameters of the current carrier motion based on the high-frequency speed measurement results, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; Through the correlation information, each speed measurement result is eliminated and weighted to obtain the best fitting order and the best fitting time window; The polynomial fitting is performed using the determined best fitting order and best fitting time window.

2. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 1 is characterized in that: Acquire high-frequency GNSS data and perform data preprocessing; the data preprocessing includes cycle slip detection of high-frequency GNSS data, and eliminating satellites and carrier phase observations that have cycle slips between adjacent epochs; then select a common satellite set between epochs to obtain a carrier phase single-difference observation equation group between epochs that meets preset conditions; and perform a step of pre-verification gross error elimination based on the carrier-to-noise ratio SNR information of the observed quantity and the speed prediction information of the previous beat, including the following steps: Receiver r receives multi-frequency signals from n satellites, then receiver r receives satellite signals from system S at time k and k+1 respectively. Corresponding frequency The carrier phase observation equation on is expressed as equation (1), where the system S includes GPS, Galileo or BDS system; （1） In formula (1), the receiver r and the satellite In the ECEF framework, the positions at time k and k+1 are , , and , The code of the observation satellite. for One of them; The subscripts k and k+1 of each variable are the corresponding moments of each variable; For observation satellite The corresponding frequency The carrier phase observation value on , and the unit is m; From the receiver r to the observation satellite The Euclidean distance of is the speed of light; is the mobile station receiver clock error; For observation satellite The clock difference is the tropospheric delay in the phase observation at receiver r; For mobile stations Phase observation corresponding to frequency ionospheric delay; For frequency The corresponding wavelength; For receiver The undifferenced integer ambiguity of The receiver r is the frequency Carrier phase of the signal Carrier phase hardware delay; For observation satellite Frequency Carrier phase of the signal Carrier phase hardware delay; For noise.

3. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 2 is characterized in that: The high-frequency GNSS data is obtained and preprocessed; the data preprocessing includes detecting the cycle slip of the high-frequency GNSS data, eliminating the satellites and their carrier phase observations that have cycle slips between adjacent epochs; then selecting the common satellite set between epochs to obtain the carrier phase single-difference observation equation group between epochs that meets the preset conditions; the step of eliminating the gross error before verification based on the carrier-to-noise ratio SNR information of the observed quantity and the speed prediction information of the previous beat also includes the steps of: The carrier phase observation value of a single GNSS receiver r is differentiated between epochs, and the change in the position of the receiver r between adjacent epochs and the change in the clock error of the receiver r between adjacent epochs are used as the parameters to be estimated. The inter-epoch differential carrier phase observation equation is established to solve the velocity, including: The corresponding frequency of satellite p of a certain system S at receiver r at time k+1 The inter-epoch differential carrier phase observation equation on can be expressed as: (2) In formula (2), the subscript is the time k+1 corresponding to each variable, and ; For observation satellite frequency Single-difference observations between epochs on ; is the difference in the mobile station receiver clock between epochs k and k+1; For observation satellite The difference in satellite clock errors between epochs k and k+1; For the mobile station r to the satellite The difference in station-satellite distance between epochs k and k+1; Observation satellite for mobile station r The non-differenced integer ambiguity change between epochs k and k+1. If the carrier phase observation does not experience a cycle slip, this value is 0. For noise.

4. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 3 is characterized in that: After performing the inter-epoch difference of the carrier phase observation value of a single GNSS receiver r, taking the change of the position of the receiver r between adjacent epochs and the change of the clock error of the receiver r between adjacent epochs as the parameters to be estimated, and establishing the inter-epoch difference carrier phase observation equation for velocity calculation, the following steps are further performed: selecting the carrier phase observation value without cycle slip, and linearizing equation (2) to obtain: （3） In formula (3), represents the direction cosine matrix from receiver r to satellite p at time k+1, is the parameter to be estimated at time k+1, which is related to the speed There is a definite conversion formula: （4） If at time k+1, the receiver r receives m common carrier observations without cycle slips between epochs and n common satellites, then an m-column double-difference observation equation system between epochs can be formed. When both m and n are greater than or equal to 4, the solution can be obtained by Perform velocity calculation.

5. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 4 is characterized in that: The step of constructing a Kalman_TDCP speed measurement model, performing Kalman_TDCP high-frequency speed measurement on the high-frequency GNSS data after data preprocessing and after Kalman filtering, obtaining the average speed v between adjacent epochs based on the Kalman_TDCP speed measurement model, and using it as the instantaneous speed at the middle moment of the previous and next epochs, storing the instantaneous speed and corresponding time series information; after obtaining the speed measurement result, storing the associated information of the post-verification residual of the speed measurement result at the same time, the associated information including the error square, square error and SSE of the speed measurement result, satellite geometric structure factor DOP and reference observation quantity, comprises the steps of: Based on the Kalman_TDCP model, the state equation and observation equation of the discrete-time Kalman filter at time k can be expressed as: （5.1） In formula (5.1), Represents the position change information at time k and the receiver clock error epoch difference The parameter vector to be estimated is composed of Predicted value, is the state transfer matrix from time k-1 to time k. High-frequency TDCP is used for velocity measurement. The position change between epochs at time t-1 is equal to the position change between the epochs before and after time t. In this model, ; are system process noise and observation noise respectively, which obey Gaussian distribution with mean 0; is the observation vector at time k, corresponding to the sum of the difference between the epoch carrier phase observations and the difference between the epoch satellite clock errors; H is the observation matrix, which can be constructed based on formula (3); E is the fourth-order unit matrix.

6. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 4 is characterized in that: The Kalman_TDCP speed measurement model is constructed, and the high-frequency GNSS data that has undergone data preprocessing is subjected to Kalman_TDCP high-frequency speed measurement after Kalman filtering. Based on the Kalman_TDCP speed measurement model, the average speed v between adjacent epochs is obtained, and the average speed v is used as the instantaneous speed at the middle moment of the previous and next epochs, and the instantaneous speed and the corresponding time series information are stored; after obtaining the speed measurement result, the associated information of the post-verification residual of the speed measurement result is stored at the same time, and the associated information includes the error square, square error and SSE of the speed measurement result, the satellite geometric structure factor DOP and the reference observation quantity, and also includes the steps of: After obtaining the state prior estimate through time update by formula (5.1), it is necessary to use the actual measurement value to correct the state prior estimate obtained by formula (5.1) to obtain the optimal estimate of the system state: the calculation formula is as follows: （5.2） In formula (5.2), the first two steps are generally called time updates, and the last three steps are measurement updates; is the covariance matrix of the system process noise; is the covariance matrix of the observation noise; is the prior estimate of the state The weight matrix of is the coefficient matrix, called the Kalman filter gain; is the optimal state estimate and its mean square error matrix; is the fourth-order identity matrix; When the optimal state estimate is obtained After that, the Kalman_TDCP speed measurement result can be obtained through formula (4).

7. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 6 is characterized in that: The method of performing scene adaptive fitting on the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results includes: calculating the dynamic state parameters of the current carrier motion based on the high-frequency speed measurement results, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; and eliminating and weighting each speed measurement result through associated information to obtain the determined optimal fitting order and optimal fitting time window, including the steps of: Through high-frequency Kalman_TDCP speed measurement, obtain arrive m speed measurement results and their associated information at the time, and m calculation samples are obtained ; Assume that the model estimate is , then the calculation formulas for the total sum of squares SST, error sum of squares SSE and regression sum of squares SSR of the calculated sample are as follows: （6） Its coefficient of determination The calculation formula is as follows: （7） The value range is between 0 and 1; Coefficient of determination Set a confidence threshold using the sum of squared errors (SSE) and , when the fitting order is iteratively fitted from low to high order, or the fitting time window is iteratively fitted from short to long, if or , the fitting order is considered to be optimal. If it is always greater than the confidence threshold, the minimum value of the coefficient of determination is selected. as the best fit order.

8. The high-precision real-time single-point speed measurement method based on high-frequency filtering fitting according to claim 7 is characterized in that: The method of performing scene adaptive fitting on the high-frequency speed measurement results to obtain high-precision and delay-free speed measurement results includes: calculating the dynamic state parameters of the current carrier motion based on the high-frequency speed measurement results, determining the calculation criteria of the polynomial fitting, and determining the optimal polynomial fitting order and the optimal time window fitting length based on the calculation criteria; eliminating and weighting each speed measurement result through associated information to obtain the determined optimal fitting order and optimal fitting time window, and also includes the steps of: Set the minimum fitting time window Restrictions to improve polynomial fitting stability; A polynomial fitting of greater than 5th order is performed on the total samples to obtain the dynamic state parameters of the carrier motion. The dynamic state parameters include velocity range, acceleration and jerk information, and the upper and lower limits of the fitting order are determined.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting as described in any one of claims 1 to 8 are implemented.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the steps of a high-precision real-time single-point speed measurement method based on high-frequency filtering fitting as described in any one of claims 1 to 8 are implemented.

Images