Disclosure of Invention
The invention aims to overcome the defects of the prior art and provide a multimode ultra-fast and low-power-consumption positioning receiver system and a method. The system provides millisecond ultra-fast (TTFF), intermittent and continuous, low-power, autonomous and reliable location services without receiver location information and without real-time wireless network connection.
The purpose of the invention is realized by the following technical scheme: a multi-mode, ultra-fast, low power consumption positioning receiver system comprising:
radio frequency front end device: the intermediate frequency signal sampling value of the GNSS satellite used for obtaining I/Q;
a storage device: the sampling device is used for storing the intermediate frequency signal sampling value;
a parameter acquisition device: the method comprises the steps of obtaining parameters including satellite ephemeris parameters, atmospheric model calibration parameters and prior information parameters;
the clock acquisition device: a time of arrival, TOA, of a receiver for obtaining the I/Q samples;
the measured value calculating means: the device is used for calculating the pseudo range and the Doppler measurement value of the multimode satellite system;
the PVT computing device: the device is used for calculating the accurate receiver position, the accurate satellite transmitting time TOT and the accurate receiver clock error under the conditions of no prior information and prior information, calculating the prior information parameters and storing the calculated prior information parameters into the parameter acquisition device.
The measurement value calculation device comprises:
a capture configuration parameter acquisition module: the acquisition configuration parameters comprise a satellite type, a satellite number, an acquisition mode, a code phase and initial Doppler;
the satellite judging and configuring module: the satellite acquisition module is used for acquiring the acquisition configuration parameters of the satellite, judging whether the satellite is configured or not according to the acquisition configuration parameters, entering the satellite fine acquisition module for processing if the satellite is configured, and entering the differential rapid acquisition satellite detection module for processing if the satellite is not configured;
the differential rapid acquisition satellite detection module comprises: the method is used for rapidly detecting the visible satellite of the current GNSS and obtaining a coarse code phase when the current prior information is unreliable and blind acquisition is needed;
a satellite coarse acquisition module: for obtaining a coarse doppler frequency and code phase;
a satellite fine acquisition module: further searching the vicinity of the rough Doppler frequency and the code phase to obtain the precise Doppler frequency and the precise code phase based on a rough acquisition module or acquisition configuration information;
the multi-phase code identification compensation module: the method is used for identifying multi-path codes of multi-path correlation peak values before and after the maximum correlation peak obtained by the satellite fine acquisition module, so as to obtain a compensation value of the current code phase and obtain a fraction pseudo range.
The differential rapid acquisition satellite detection module comprises:
a mixing packing unit: the device is used for mixing the stored I/Q original intermediate frequency data with a local fundamental frequency carrier signal, and then packaging the code phase according to Chip correlation integral;
a conjugate difference calculation unit: the device is used for carrying out forward and backward conjugate difference calculation on the packed integration result according to a Chip as a unit;
a coherent accumulation unit: the coherent accumulation is carried out according to the code phase value corresponding to the code period;
a new pseudo code acquisition unit: the method comprises the steps of carrying out conjugate difference on C/A codes of all generated satellites to obtain new pseudo codes of the satellites;
a cyclic correlation unit: the device is used for circularly correlating the coherent accumulation value obtained by the coherent accumulation unit with the new pseudo codes of different satellites obtained by the new pseudo code acquisition unit according to a code period;
a first peak-to-average ratio determination unit: judging whether the peak-to-average ratio after the corresponding satellite circular correlation exceeds a threshold: if the value is larger than the threshold, the acquisition is successful, and the coarse code phase of the corresponding code is obtained.
The satellite coarse acquisition module and the satellite fine acquisition module adopt the same structure and comprise:
a Chip packing result acquisition unit: the Chip packing result is used for obtaining the frequency mixing packing unit;
a correlation calculation unit: the device is used for carrying out correlation calculation on the integrated result packed by the Chip and the satellite C/A code to be captured;
correlation result halving unit: the device is used for averaging the correlation result of N milliseconds into P parts for coherent integration to obtain P coherent integration I/Q values;
a matching coherent unit: the method is used for carrying out millisecond coincidence matching coherence on N milliseconds of halved P coherent integration I/Q values based on a millisecond matching list;
FFT and maximum value finding unit: the matching unit is used for carrying out P-point FFT operation on the matching result of the matching coherent unit and searching the maximum value of the current matching amplitude;
a matching result storage unit: the matching result with the maximum amplitude is selected for storage;
differential conjugate coherent integration unit: the method is used for carrying out M differential conjugate coherent integrations, wherein the difference mode at each time is as follows: the phase difference is carried out on the P point of the N millisecond at the nth time and the P point of the N millisecond at the N +1 time, and corresponding phase amplitude value accumulation is carried out after the difference is finished;
a second peak-to-average ratio determination unit: the method is used for judging the peak-to-average ratio of the search space and judging whether the peak-to-average ratio exceeds a threshold: if the threshold is exceeded, the corresponding acquired Doppler frequency and code phase are output. Here, the search space is a two-dimensional space of codes and doppler.
The multi-phase code discrimination compensation module comprises:
code phase difference calculation unit: the device is used for selecting a plurality of correlation measurement value amplitudes with maximum code phase left-right symmetry under the frequency corresponding to the maximum correlation peak value, and obtaining a code phase difference value after normalization; the code phase difference value is a compensation value;
fractional pseudorange calculation unit: and the satellite fine acquisition module is used for calculating the fractional pseudorange according to the phase difference value and the phase value obtained by the satellite fine acquisition module.
The PVT computing device comprises:
a parameter acquisition module: the system is used for acquiring satellite ephemeris parameters, atmospheric model calibration parameters, time intervals, a coarse time TOA, a fractional pseudorange, Doppler frequency and position information;
a priori information judgment module: judging whether the first positioning/power-off time is too long or short-time power-off/intermittent power-off positioning/continuous positioning: if the first positioning/power-off time is too long, no prior information is considered, and the PVT calculation module without the prior information is used for processing; if the short-time power failure/intermittent power failure positioning/continuous positioning exists, the prior information is considered to exist, and the PVT calculation module with the prior information is used for processing;
the PVT calculation module without prior information: the method is used for carrying out least square iterative solution under the condition of no prior information so as to obtain an accurate receiver position, an accurate TOT and a receiver clock error;
the PVT calculation module with prior information: the method is used for calculating the PVT by using the constructed accurate full pseudorange measurement value Kalman filtering under the condition of prior information to obtain the accurate receiver position, the accurate TOT and the receiver clock error.
The PVT calculation module without prior information comprises:
a coarse position calculation unit: based on the Doppler frequency measurement value, using the least square iteration supporting multimode to calculate PVT, and obtaining a rough receiver position, a corrected rough TOA and a receiver clock drift;
a first full pseudo-range calculation unit: constructing a full pseudorange measurement value by using the position of the coarse receiver, the corrected coarse TOA, ephemeris information and a fractional pseudorange and adopting a suppression TOT and a position error;
a precise position calculation unit: based on the full pseudorange measurement value, the PVT is iteratively calculated by using the multimode-supporting least square iteration, and the accurate receiver position, the accurate TOT and the receiver clock error are obtained.
The coarse position calculating unit includes:
a parameter acquisition subunit: the method comprises the steps of obtaining input values required by current calculation, wherein the input values comprise a rough TOA, a Doppler frequency measurement value and an ephemeris parameter;
a data calculation subunit: the system is used for calculating and obtaining the satellite position, the satellite speed, the satellite clock error and the satellite clock drift based on the rough TOA and the ephemeris parameters;
receiver initial position calculation subunit: the device is used for respectively calculating the average value of each direction of XYZ under an ECEF coordinate system based on the position of each multimode satellite obtained by the data calculation subunit; then converting the average value into a WGS84 coordinate system, setting the altitude to be 0, then converting into an ECEF coordinate system, and taking a calculation result as an initial value of the position of the receiver;
doppler positioning iteration subunit: the device is used for performing five-state parameter least square Doppler positioning iteration supporting multimode based on the receiver position initial value and the Doppler frequency measured value;
an iteration judgment subunit: and the method is used for judging whether the iteration converges or not, outputting the rough receiver position and the corrected TOA if the iteration converges, and exiting otherwise.
The first full pseudo-range calculation unit includes:
satellite position calculation subunit: the TOA correction module is used for calculating the satellite position corresponding to the moment by using the corrected TOA obtained by the rough positioning;
a predicted geometric distance calculation subunit: the system comprises a receiver, a satellite position estimation unit and a satellite position estimation unit, wherein the receiver is used for estimating the position of the satellite;
initial reference satellite integer millisecond calculation subunit: for selecting an initial reference satellite and calculating an integer number of milliseconds of the initial reference satellite; the initial reference satellite is the satellite with the highest elevation angle;
a satellite integer millisecond calculation subunit with the shortest geometric distance from the reference satellite: the satellite selection method is used for selecting the satellite with the shortest geometric distance to the reference satellite and calculating an integer of milliseconds;
a full pseudo-range calculation subunit: for constructing full pseudoranges for all satellites.
The precise position calculating unit includes:
an initial value calculating operator unit: the system comprises a receiver, a GPS clock error correction module, a BDS clock error correction module, a positioning module and a positioning module, wherein the receiver is used for calculating a rough position, obtaining a rough receiver position and a corrected receiver time TOA, and the GPS clock error and the BDS clock error are 0 and are used as iteration initial variables; taking the full pseudorange constructed by the first full pseudorange calculation unit as a pseudorange measurement value;
predicted pseudorange calculation subunit: a pseudorange for deriving a prediction from the coarse position and the satellite coordinates;
a state variable update value operator unit: the state variable updating device is used for calculating a state variable updating value according to the deviation of the predicted pseudo range and the constructed full pseudo range measurement value;
a state variable update subunit: the state variable updating module is used for updating the state variable according to the state variable updating value;
an iteration convergence judgment subunit: and the pseudo range prediction calculation subunit is used for judging whether a convergence condition is met or not, if so, exiting, and otherwise, returning to the pseudo range prediction calculation subunit for calculation.
The PVT calculation module with the prior information comprises:
a second full pseudo-range calculation unit: the method comprises the steps of generating a full pseudorange measurement value after obtaining accurate historical TOA time, receiver clock error, receiver clock drift and historical receiver position information;
a kalman accurate position calculating unit: for computing PVT by kalman filtering, using full pseudorange measurements.
The second full pseudo-range calculation unit includes:
receiver clock error estimation subunit: the clock correction method is used for estimating the current receiver clock error according to the historical receiver clock error and the historical receiver clock drift time interval;
satellite transmission time estimation subunit: estimating the time of transmission TOT of the satellite according to the estimated current receiver clock error, the historical accurate TOA and the locally estimated time interval;
the multimode satellite TOT time position and clock error calculating subunit: the clock correction method comprises the steps of calculating satellite positions and clock errors of satellites at the TOT time of a multimode satellite respectively;
the expected geometric distance calculation subunit: for obtaining an estimated predicted geometric distance between the satellite and the receiver by calculation;
satellite integer millisecond calculation subunit: for calculating an integer number of milliseconds for each satellite;
a full pseudo-range calculation subunit: for constructing full pseudoranges for all satellites.
The Kalman accurate position calculation unit comprises:
the discrete state equation and the measurement equation establish a subunit: the system is used for establishing a discrete state equation and a measurement equation;
a Kalman filtering prediction subunit: the method is used for predicting the state value of the current epoch by using the state equation of the system on the basis of the state estimation value of the last epoch;
a Kalman filtering updating subunit: for correcting the state prior estimate obtained by the prediction process using the actual measurement.
A multimode, ultra-fast, low-power consumption positioning receiver method, comprising the steps of:
s1: in the arrival time of the interrupt signal, the measured value calculating device locks the intermediate frequency I/Q data converted by the radio frequency front end device and gradually stores the intermediate frequency I/Q data into the storage device, and simultaneously locks the current time of the clock acquisition device as a rough TOA;
s2: the measurement value calculation device acquires configuration parameters related to acquisition and calculates pseudo range and Doppler measurement values of the multimode satellite system;
s3: the PVT calculating device acquires Doppler measured values, rough TOA, ephemeris and atmospheric correction parameters, judges whether prior information is reliable or not, and calculates accurate receiver position, satellite transmitting time TOT and receiver clock error under the conditions of no prior information and prior information; and calculating prior information parameters and storing the calculated prior information parameters into the parameter acquisition device.
The step S2 includes the following sub-steps:
s21: acquiring configuration parameters captured based on prior information, and configuring different capturing modes;
s22: judging whether the satellite is configured, if so, entering the step S25, otherwise, entering the step S23;
s23: differential rapid acquisition satellite detection: rapidly detecting a visible satellite of a current GNSS and obtaining a coarse code phase;
s24: satellite coarse acquisition: obtaining a coarse doppler frequency and a code phase;
s25: satellite fine acquisition: further searching the vicinity of the rough Doppler frequency and the code phase to obtain the precise Doppler frequency and the precise code phase based on a rough acquisition module or acquisition configuration information;
s26: multi-phase code discrimination compensation: and performing multi-path code identification on multi-path correlation peak values before and after the maximum correlation peak obtained by the satellite fine acquisition module so as to obtain a compensation value of the current code phase and obtain a fractional pseudo range.
The step S23 includes the following sub-steps:
s231: mixing the stored I/Q original intermediate frequency data with a local fundamental frequency carrier signal, and then packaging code phases according to Chip correlation integral;
s232: carrying out forward and backward conjugate differential calculation on the packed integration result by taking a Chip as a unit, and storing the result for later phase accumulation;
s233: coherent accumulation is carried out according to code phase values corresponding to the code period, and the code phase values are stored;
s234: carrying out forward and backward conjugate difference on the C/A codes generated by each satellite to obtain a new pseudo code of the satellite;
s235: circularly correlating the coherent accumulated values stored in the step S233 with the new pseudo codes of different satellites in the step S234 according to a code period;
s236: judging whether the peak-to-average ratio after the corresponding satellite circular correlation exceeds a threshold: if the value is larger than the threshold, the acquisition is successful, and a coarse phase value of the corresponding code is obtained.
The step S24 includes the following sub-steps:
s241: inputting the Chip correlation integral packaging result obtained in the step S231;
s242: performing correlation calculation on the Chip-packed integration result and the satellite C/A code to be captured;
s243: p halving the N millisecond correlation results and performing coherent integration to obtain P coherent integration I/Q values;
s244: based on the millisecond matching list, carrying out millisecond coincidence matching coherence on P coherent integration I/Q values of N milliseconds;
s245: performing P-point FFT operation on the matching result, and searching the maximum value of the current matching amplitude;
s246: after matching of the matching list is completed, selecting a matching result with the maximum amplitude value for storage;
s247: performing M differential conjugate coherent integrations, wherein the differential mode is as follows: the phase difference is carried out on the P point of the N milliseconds at the 1 st time and the P point of the N milliseconds at the 2 nd time, and the like; after the difference is finished, corresponding phase amplitude values are accumulated;
s248: judging the peak-to-average ratio of the search space, and judging whether the peak-to-average ratio exceeds a threshold: and if the threshold is exceeded, outputting the corresponding acquisition frequency and code phase.
The step S26 includes the following sub-steps:
s261: selecting a plurality of correlation measurement value amplitudes with the maximum code phase left-right symmetry under the frequency corresponding to the maximum correlation peak value, and obtaining a code phase difference value after normalization; the code phase difference value is a compensation value;
s262: and calculating a fraction pseudo range according to the phase difference value and the phase value obtained by the satellite fine acquisition module.
The step S3 includes the following sub-steps:
s31: acquiring satellite ephemeris parameters, atmospheric model calibration parameters, time intervals, coarse time TOA, fractional pseudoranges, Doppler frequency and position information;
s32: judging whether the first positioning/power-off time is too long or short-time power-off/intermittent power-off positioning/continuous positioning: if the first positioning/power-off time is too long, no prior information is considered, and the step S33 is entered; if the positioning is short-time power-off/intermittent power-off/continuous positioning, the prior information is considered to be available, and the step S34 is entered;
s33: under the condition of no prior information, performing least square iterative solution to obtain an accurate receiver position, an accurate TOT and a receiver clock error;
s34: under the condition of prior information, the constructed accurate full pseudorange measurement value Kalman filtering is used for calculating PVT, and the accurate receiver position, the accurate TOT and the receiver clock error are obtained.
The step S33 includes the following sub-steps:
s331: based on the Doppler frequency measurement value, using the least square iteration supporting multimode to calculate PVT, and obtaining a rough receiver position, a corrected rough TOA and a receiver clock drift;
s332: constructing a full pseudorange measurement value by using the position of the coarse receiver, the corrected coarse TOA, ephemeris information and a fractional pseudorange and adopting a suppression TOT and a position error;
s333: based on the full pseudorange measurement value, the PVT is iteratively calculated by using the multimode-supporting least square iteration, and the accurate receiver position, the accurate TOT and the receiver clock error are obtained.
The step S331 includes the following sub-steps:
s3311: obtaining data required by calculation, wherein the data comprises a rough TOA (receiver time of arrival), a Doppler frequency measurement value and an ephemeris parameter;
s3312: calculating to obtain the position and the speed of the satellite, the clock error and the clock drift of the satellite by using the rough time TOA of the receiver and the ephemeris parameters;
s3313: based on the positions of the multimode satellites obtained in step S3312, the average values in each XYZ direction are respectively obtained in the ECEF coordinate system; then converting the average value into a WGS84 coordinate system, setting the altitude to be 0, then converting into an ECEF coordinate system, and taking a calculation result as an initial value of the position of the receiver;
s3314: based on the initial value of the receiver position and the Doppler frequency measurement value, performing five-state parameter least square Doppler positioning iteration supporting multiple modes;
s3315: and judging whether the iteration converges, if so, outputting the rough receiver position and the corrected TOA, and otherwise, exiting.
The step S332 includes the following substeps:
s3321: calculating the satellite position corresponding to the moment by using the corrected TOA obtained by rough positioning;
s3322: calculating to obtain an estimated predicted geometric distance between the satellite and the receiver by using the satellite position and the rough receiver position;
s3323: based on selecting an initial reference satellite and calculating an integer number of milliseconds of the initial reference satellite; the initial reference satellite is the satellite with the highest elevation angle;
s3324: selecting a satellite with the shortest geometric distance from a reference satellite, and calculating an integer of milliseconds;
s3325: a full pseudo-range calculation subunit: for constructing full pseudoranges for all satellites.
The step S333 includes the following sub-steps:
s3331: taking the rough receiver position obtained by calculation in the step S331, the corrected receiver time TOA, the GPS clock error and the BDS clock error as 0 as iteration initial variables, and simultaneously taking the accurate full pseudorange constructed in the step S332 as a pseudorange measurement value;
s3332: obtaining a predicted pseudorange from the coarse position and the satellite coordinates;
s3333: calculating an updated value of the state variable according to the deviation of the predicted pseudo range and the constructed full pseudo range measured value;
s3334: updating the state variable according to the state variable update value obtained in the step S3333;
s3335: judging whether an iteration convergence condition is met: and exiting if the condition is met, and returning to the step S3332 if the condition is not met.
The step S34 includes the following sub-steps:
s341: constructing an accurate full pseudorange measurement value based on the historical position, the clock error, the clock drift and the time interval;
s342: and (4) calculating PVT by using Kalman filtering of the full pseudorange measurement value, and acquiring an accurate receiver position, an accurate TOT and a receiver clock error.
The step S341 includes the following sub-steps:
s3411: estimating the clock error of the current receiver according to the clock error of the historical receiver and the clock drift time interval of the historical receiver;
s3412: estimating the sending time of the satellite according to the estimated current receiver clock error, the historical accurate TOA and the locally estimated time interval;
s3413: respectively calculating the satellite position of the TOT moment of the multimode satellite and the clock error of the satellite;
s3414: calculating to obtain the estimated predicted geometric distance between the satellite and the receiver;
s3415: calculating an integer millisecond of each satellite;
s3416: full pseudoranges are constructed for all satellites.
The step S342 includes the following sub-steps:
s3421: establishing a discrete state equation and a measurement equation;
s3422: kalman filtering iteration, comprising the following substeps:
s34221: on the basis of the state estimation value of the previous epoch, predicting the state value of the current epoch by using a state equation of the system;
s34222: the actual measurements are used to correct the state prior estimates from the prediction process. .
The invention has the beneficial effects that:
1. compared with the traditional receiver, the I/Q value of the short millisecond intermediate frequency data sampling is utilized, so that the data processing amount is greatly reduced, and the power consumption is reduced.
2. The measurement calculation device provides a method for supporting low-power consumption and ultra-fast acquisition of pseudo-range and Doppler measurement values of the multimode satellite system. Therefore, the problems that the tracking convergence speed and the measurement precision of the traditional tracking receiver are contradictory, the applicable scenes are reduced due to insufficient satellite number of a single mode, and the positioning speed is reduced due to the increase of the satellite number dimension of a multi-model system are solved. The main characteristics and the improvement points are as follows:
(1) the method only needs to carry out once simple coherent integration, then carries out cyclic correlation of corresponding codes on the satellite to be checked, and does not need operations with high complexity, such as FFT (fast Fourier transform algorithm) and the like, thereby reducing the operation amount and greatly improving the detection speed. The method is particularly suitable for the situation that the number of multimode satellites is increased.
(2) The structure of the satellite coarse capture and the satellite fine capture is unified, various functional requirements are realized through flexible configuration, and hardware resources are saved.
(3) The satellite coarse acquisition is based on the coarse Doppler frequency and the code phase of the coarse acquisition for further searching, so that the search space is reduced, the calculation amount is reduced, and the speed is increased.
(4) By utilizing the millisecond matching method, the limit of the length of coherent integration and the precision of frequency FFT acquisition caused by the difference of multimode signal modulation modes is broken through. The method carries out matching based on a possible matching list, expands the coherent integration length, increases the incoherent times and is compatible with the characteristics of multimode signals.
(5) The input required by the multi-phase code identification is that under the same frequency corresponding to the maximum correlation peak in the fine capture, the front and back symmetrical multi-channel code phase correlation peak values corresponding to the maximum phase are directly obtained, extra calculation is not needed, and the calculation amount is reduced.
(6) It has the effect of suppressing multipath from the two discrimination formulas for discriminating code phase errors. Therefore, the effect of multipath inhibition is achieved on the accurate code phase obtained from capturing, the positioning accuracy is improved, and the applicable scene range is enlarged.
(7) The main calculation amount of the multi-phase code discrimination is to discriminate the code phase error, and only 3 times of subtraction is done in terms of its formula. Compared with the method of performing curve fitting by using a curve function and performing iterative operation to obtain the code phase error, most of the operation is saved.
(8) According to different application scenes and prior information, flexible and quick capture configuration is provided, so that the search space is reduced, the speed of obtaining a measured value is increased, and the time required for positioning is reduced. Further reducing the amount of calculation and reducing the power consumption.
3. The PVT computing device is applied to different scenes respectively without receiver position information and real-time connection of a wireless network, so that the calculated amount is reduced, the positioning speed and the positioning precision are improved, and the positioning reliability and the robustness are improved. Specifically, the method comprises the following steps:
(1) the five-state variable Doppler positioning and six-state variable precise positioning method supporting multimode is improved. Compared with a single mode, the number of satellites is increased, the positioning geometric factor is improved, and the positioning accuracy is improved. On the other hand, the use scene of positioning is also increased. For example, in an occluded scene, a single mode is very likely to result in an insufficient number of positioning satellites.
(2) Due to the requirement of the LS iterative algorithm for the initial values: the more accurate the initial value, the faster the iterative convergence speed, and the lower the probability that the result is locally optimal. Compared with the traditional method which uses the geocentric as an initial value, the method has the advantages that the average value of the satellite positions projected on the ground is used as the initial value in the rough position calculation, so that the iteration times are reduced, the convergence speed is increased, and the accuracy of the rough position is improved.
(3) And analyzing the characteristic analysis of the estimation error generated by the satellite by the rough position and the rough time error, and selecting the satellite with the highest elevation angle and the satellite with the shortest distance to carry out gradual full pseudorange construction. Therefore, the first full pseudorange calculation has the performance of restraining the TOT and the position error to improve the construction of the full pseudorange measurement value under the condition of no prior information. Thereby improving the positioning usability and reliability.
(4) Depending on the differences and characteristics of the usage scenarios, within minutes of the positioning interval. The first full pseudorange calculation is improved by performing TOT estimation using time interval, clock bias, fractional pseudorange, and other information, and then rounding. Thereby estimating an accurate TOT and thus constructing an accurate full pseudorange. Compared with the first full pseudo range calculation, the improved second full pseudo range calculation reduces the calculation amount and improves the calculation speed. And the construction of the full pseudo range satellites is independent, and the estimation errors of other satellites caused by the estimation error of one satellite are avoided.
(5) And improving the multi-mode supporting six-state variable least square positioning calculation during intermittent outage or short-time outage. And using Kalman filtering solution supporting multimode to solve the problem that multiple iterations are needed in the least square positioning solution supporting multimode six state variables. Each iteration requires recalculation of satellite position and velocity information due to TOT updates. And the Kalman filtering solution supporting multimode only needs one-time prediction and updating, thereby further reducing the operation and improving the speed. In addition, compared with a least square algorithm, the filtering algorithm essentially improves the positioning precision and robustness.
4. The invention only needs millisecond-level short millisecond intermediate frequency data for each positioning, and the radio frequency can be powered off after the intermediate frequency data is obtained. The whole positioning process only needs the positioning time of millisecond, so the continuous positioning with the updating rate of 1s can be supported. Corresponding to some low-power-consumption applications, positioning is carried out once every few seconds or minutes, the system only needs millisecond-level time, only the clock device works at other time, and other equipment can be completely powered off, so that the power consumption is reduced.
Detailed Description
The technical scheme of the invention is further described in detail by combining the attached drawings:
when using short millisecond intermediate frequency data, of the order of milliseconds of sample length, GNSS receivers are enabled to operate at very low power and compute position very quickly. It is known that such short if data makes it impossible for the receiver to decode the navigation data of the satellite using such short sample length if data, thereby obtaining time information of the satellite. In addition, at such short sample lengths the frequency data is much lower than the convergence time required to employ a tracking loop, requiring other accurate estimation techniques to obtain pseudorange, doppler measurements.
The present embodiments describe a system and method for determining position, velocity, and time (PVT) from signal samples of a Global Navigation Satellite System (GNSS).
FIG. 1 is a system block diagram of an embodiment of the present invention, which includes: (1) a GNSS sample acquisition apparatus (radio frequency front end) for obtaining in-phase/quadrature (I/Q) intermediate frequency signal samples of GNSS satellites; (2) a storage device for storing the intermediate frequency signal acquisition value; (3) the parameter acquisition device is used for acquiring satellite ephemeris data, satellite clock parameters, atmosphere correction parameters, relevant prior information configuration parameters and the like; (4) a clock acquisition means for obtaining receiver time of reception (TOA) of said I/Q samples; (5) a measurement value calculation device, which provides a device for supporting low-power consumption and ultra-fast acquisition of pseudo range and Doppler measurement values of a multimode satellite system; (6) a PVT computing device is applied to different scenes respectively without receiver position information and real-time connection to a wireless network, and provides PVT related information.
Wherein the measurement value calculation means respectively include: (1) acquiring capture configuration parameters; (2) differential rapid acquisition satellite detection is used for rapid detection of visible satellites; (3) the satellite coarse acquisition is used for obtaining a coarse code phase and a Doppler measured value of the satellite without meeting the positioning requirement; (4) and (4) satellite fine acquisition, and further searching near the coarse acquisition result to obtain more accurate code phase and Doppler frequency value. (5) And (3) multi-phase code discrimination compensation, namely performing code phase difference discrimination with multipath effect resistance on the code phase result obtained by fine acquisition to obtain a phase compensation value. Thereby obtaining more accurate fractional pseudorange measurements.
In addition, the PVT computing device is divided into two mainlines, which are respectively adapted to different application scenarios, and respectively include:
first, no prior information (first positioning, or power-off time too long); (1) and (4) rough position calculation, namely obtaining a rough receiver position and a corrected rough TOA and a receiver clock drift by using Doppler measurement values. (2) A first full pseudorange calculation, using the suppressed TOT and the position error to construct a full pseudorange measurement. (3) And (4) calculating the accurate position, and acquiring the accurate position of the receiver, the accurate TOT and the receiver clock error by using the full pseudo-range measured value.
Secondly, prior information (short-time power-off, intermittent power-off positioning and continuous positioning) exists; (1) the second full pseudorange calculation, an improved method relative to the first full pseudorange calculation, also aims to construct an accurate full pseudorange measurement. (2) Calculating a Kalman fine position: and (4) calculating PVT by using Kalman filtering of the full pseudorange measurement value, and acquiring an accurate receiver position, an accurate TOT and a receiver clock error.
Fig. 2 is a block diagram of a hardware system of the present invention, which shows a block diagram of a hardware system according to an embodiment of the present invention. Correspondingly, the system may comprise a radio frequency front end that down-converts the analog radio frequency signals RF of the GNSS antenna into digital intermediate frequency signals IF. The intermediate frequency signals are then stored in a RAM, which is connected to an FPGA/IC for processing the time information corresponding to the measurements and intermediate frequency I/Q of the GNSS satellites, the time information being obtained by the RTC/crystal oscillator. The satellite measurement values are sent to a digital signal processor DSP/microprocessor ARM for PVT related operations. And the read-only memory/erasable programmable read-only memory ROM is respectively used for storing the firmware, including a processing algorithm and parameter information such as a plurality of prior information, ephemeris information, satellite clock error and the like.
FIG. 3 is a flow chart of the method of the present invention, and schematically shows a data flow chart of the system. When starting, the clock device is locked by the interrupt signal to obtain the rough time of reception (TOA) of the I/Q data, and the intermediate frequency data is locked and stored in the storage device, and the parameter information of the parameter device is obtained. The acquisition related configuration is used for the measurement calculation device and the ephemeris, atmospheric correction parameters are used for the PVT calculation device. The measurement value calculation device utilizes the intermediate frequency data I/Q to quickly acquire the measurement value of the GNSS satellite. The measured values, the coarse time of reception (TOA), the ephemeris, and the atmospheric correction parameters are subjected to PVT correlation. After obtaining the position, speed and time information, the current relevant information is calculated, and the prior information parameters are calculated and stored in the parameter equipment for the next positioning.
Wherein, at the arrival time of the interrupt signal, the measured value calculating means locks the intermediate frequency I/Q and gradually stores it in the storage device as the input information of the post-capture measurement. And simultaneously lock the current time of the clock device as a coarse TOA. At the same time, the parameters of the parameter means, such as ephemeris parameters (satellite orbit and satellite clock parameters), satellite predicted orbit parameters, satellite clock calibration, atmospheric model calibration parameters (ionosphere, troposphere) are read.
Fig. 4 is a flowchart of a calculation method of the measured value calculation apparatus.
The availability, reliability and accuracy of GNSS receivers depend to a large extent on the accuracy of the measurements. Conventional GNSS receivers, after successful bit and subframe synchronization, also ensure correct decoding of ephemeris to provide various parameters of the satellite, pseudoranges, doppler and carrier phase measurements.
However, the tracking loop can only use a short integration time before bit synchronization, resulting in high measurement noise and long convergence time. Once the electricity is short midway, the lock needs to be re-entered again. The tracking loop is generally not sufficient because the length of the sampled data is much less than the convergence time and the accuracy of the initial doppler and code phase is too poor for positioning, which section details a multi-stage measurement scheme with different variable modes configured based on different a priori scenarios, improving reliable fractional pseudorange and doppler measurements. The present embodiment mainly comprises 5 parts: acquiring acquisition configuration parameters, detecting a differential rapid acquisition satellite, roughly acquiring the satellite, finely acquiring the satellite and identifying and compensating the multi-phase code. The function of the 5 sections is outlined as follows.
First, obtain and catch the configuration parameter
In the embodiment, the acquisition configuration parameters are mainly based on prior information, different acquisition modes are configured, the acquisition time is shortened, and the acquisition sensitivity is improved. The prior information comprises information such as estimated satellite Doppler range, satellite code phase range, clock drift and the like.
Two, differential type fast acquisition satellite detection
The differential rapid acquisition satellite detection leads to the increase of the satellite dimensional search space due to the increase of the number of multimode satellites. Invisible satellites are eliminated through rapid searching, and therefore the acquisition speed is improved.
Correspondingly, when the differential rapid acquisition satellite detection is used for current prior information unreliable and blind acquisition is needed, the current GNSS visible satellite is rapidly detected, the coarse code phase is obtained, the search range of the code phase acquired in the later period is reduced, the number of satellites required to be acquired in the later period for coarse acquisition and fine acquisition of the satellite is reduced, and therefore the acquisition speed is improved. The method only needs simple coherent integration and does not need high-complexity operations such as FFT (fast Fourier transform) and the like, thereby reducing the operation amount.
After I/Q intermediate frequency data are input through a radio frequency front end, an intermediate frequency data complex signal after chip-level coherent integration packaging is as follows:
wherein A is the signal amplitude, InIs the nth Chip coherent integration value, QnIs the nth Chip coherent integration value, x (n) is the satellite C/A code, D (n) is the satellite navigation message, fiUnknown, is the doppler frequency of the intermediate frequency signal after mixing with the local fundamental frequency.
Then the signal is differential by conjugate before and after:
in this embodiment, the CodeRate of the GPS L1 is 1.023MHz, and the CodeRate of the corresponding BDS is 2.046 MHz.
Since the receiver performs forward and backward conjugate difference operation on the code level coherent integration result, the newly generated Sdif(n) then a new x is generateddifX (n) x (n-m) signal, which can be considered as a new pseudo code.
New input signal S
dif(n) has the following characteristics: (1) navigation messages of the original intermediate frequency signal are eliminated. (2) For a short time (tens of milliseconds required for capture)
iCan be regarded as unchanged, then
Is a constant that does not change with time, so the doppler of the original received signal is also cancelled. (3) The conventional receiver doppler frequency fluctuation range is 10KHz, for GNSS systems,
wherein, for the feature (1), the navigation bit information is eliminated, then long-term coherent integration can be performed, thereby improving the signal acquisition sensitivity. For the characteristics (2) and (3), Doppler information is eliminated, signals only contain new pseudo codes of visible satellites and are modulated together, new pseudo codes can be constructed by directly utilizing each satellite, the coherent accumulation result is directly subjected to circular correlation accumulation, whether acquisition is successful or not is judged, and the coarse code phase of the corresponding satellite is obtained. By doing so, the code can be searched first, and then the frequency dimension can be searched, so that the searching times can be reduced, and the speed can be improved. And only once coherent integration accumulation is needed, other satellites carry out cyclic correlation acquisition satellites based on the new pseudo codes, and independent coherent integration is not needed to be carried out on each satellite, so that the operation amount is reduced, and the speed is increased.
Fig. 5 is a flowchart of differential fast acquisition satellite detection, which shows detailed steps of differential fast acquisition satellite detection.
1. Mixing the stored I/Q original intermediate frequency data with a local fundamental frequency carrier signal, and then packaging the code phase according to Chip correlation integral.
2. And carrying out forward and backward conjugate differential calculation on the packed integration result by taking a Chip as a unit, and storing the result for later phase accumulation.
3. And carrying out coherent accumulation according to the code period corresponding to the code phase value, and storing.
4. C/A codes are generated by each satellite, and front-back conjugate difference is carried out to obtain new pseudo codes of the satellites.
5. And (4) circularly correlating the coherent accumulated values stored in the step (3) with new pseudo codes of different satellites respectively according to the code period.
6. And judging whether the peak-to-average ratio after the corresponding satellite circular correlation exceeds a threshold. If the value is larger than the threshold, the acquisition is successful, and a coarse phase value of the corresponding code is obtained.
Among them, corresponding to the C/A generator for generating C/A in step 4, a better method can be to calculate and store the C/A sequence in the nonvolatile memory in advance, so that the C/A sequence is always available in the calculation process and no additional calculation is needed.
Coarse acquisition of satellite
The coarse acquisition of the satellite reduces the search space with larger frequency and code interval, improves the acquisition sensitivity by increasing the differential coherent/noncoherent times, and is determined by the signal strength and the required performance of the receiver. The doppler and code phase measurement accuracy at this stage is not required to achieve the accuracy required for positioning, and a rough estimate is sufficient. The rough measurement value is only used for searching of fine capture, and the fine capture search space is reduced, so that the capture speed is improved, and the operation amount is reduced.
Fig. 6 is a flow chart of the satellite coarse/fine acquisition, describing the details of the satellite coarse/fine acquisition implementation.
1. The Chip packing result of differential quick capture calculation is input, and the calculated value is multiplexed, so that no extra calculation processing is needed. The calculation is reduced.
2. And carrying out correlation calculation on the integrated result packed by the Chip and the C/A code of the satellite to be acquired.
3. And dividing the N millisecond correlation results into P parts and performing coherent integration to obtain P coherent integration I/Q values.
4. Dividing P coherent integration values of N milliseconds, and performing millisecond matching coherence based on a millisecond matching list.
5. And performing P-point FFT operation on the matching result, and searching the maximum value of the current matching amplitude.
6. And after matching of the matching list is completed, selecting the matching result with the maximum amplitude value for storage.
7. And carrying out M differential conjugate coherent integrations, wherein the difference mode is that the phase of the P point of the N milliseconds at the 1 st time and the phase of the P point of the N milliseconds at the 2 nd time are differentiated, and the like, and corresponding phase amplitude value accumulation is carried out after the differentiation is finished.
8. And judging the peak-to-average ratio of the search space, and judging whether the peak-to-average ratio exceeds a threshold. And if the threshold is exceeded, outputting the corresponding acquisition frequency and code phase.
For the No-GEO of BDS, there is a 20ms period NH code, GEO is a 2ms navigation bit, and GPS is a 20ms length navigation bit. The length of coherent integration is limited and the accuracy of frequency FFT acquisition is also limited. The method is matched based on the possible matching list, the coherent integration length is expanded, and the multimode signal characteristics are compatible.
Although matching adds to some extent computation, here only the first fix is affected. Once positioned, the corresponding matching information can be estimated through the prior information, so that the matching times are reduced.
According to the prior information, the capture configuration parameters can be flexibly configured, so that the optimal result can be achieved in the aspects of sensitivity, search speed and search precision. Corresponding to coarse acquisition, a coarse frequency and a coarse code phase value need to be determined quickly, and the speed is improved. Then N may be set to be smaller, P larger, and M larger. The N value determines the capture frequency interval, the P value determines the number of the search points of the captured FFT and determines the frequency search range, and the M value improves the capture sensitivity. Since the original if data used is in the order of milliseconds, it is limited. N and M need to be decided according to the actual data length.
In this embodiment, N is set to 1 in response to the requirement of coarse acquisition, and it is not necessary to perform matching operation, thereby increasing the coarse acquisition speed.
Four, fine acquisition of satellite
The satellite fine acquisition aims to obtain more accurate code phase and Doppler frequency values, and further search is carried out based on the rough Doppler frequency and the vicinity of the code phase of the rough acquisition, so that the calculation amount is reduced, and the speed is increased. The structure of the method reuses coarse capture, so that hardware resources are saved for a hardware implementation mode. The main differences are as follows:
1. the input packed data I/Q are different, the packing interval is based on the code phase precision requirement, and the limit precision is the intermediate frequency data sampling interval. The frequency of the additional mixing is not the fundamental frequency but the fundamental frequency + the coarse acquisition doppler frequency.
2. The search range of the code phase of the fine acquisition is the precision range of the code phase of the coarse acquisition.
3. And finishing the search of the fine acquisition frequency, which is the frequency precision range of the coarse acquisition.
4. The precision of the acquisition frequency is determined based on the following requirement of Doppler positioning precision, and the precision of the acquisition frequency is improved by setting a larger N value. The frequency and code phase of the fine acquisition search are in a narrow range, and even if matching acquisition is carried out, the acquisition efficiency also meets the requirement.
Discrimination compensation of five-phase and multi-phase codes
The multi-phase code identification compensation is to utilize the multi-path correlation peak values before and after the maximum correlation peak of the fine capture to identify the multi-path code, thereby obtaining the compensation value of the current code phase and further improving the precision of the code phase.
From the foregoing, the accuracy of code pseudorange measurements is a key indicator for GNSS receivers. The code pseudorange measurement value obtained corresponding to the traditional tracking needs to be obtained after a loop is stabilized, and long time is needed. And after the tracking loop is powered off, the tracking loop needs to be locked again stably, and discontinuity exists. The method is independent in measurement time at each time, has no front-back correlation, and can be obtained by only using intermediate frequency data of tens of milliseconds.
Since the fine acquisition is obtaining a more accurate code phase, the maximum accuracy can only reach half of the sampling frequency of the clock. For low power consumption products, reduced sampling frequency is required, so that the accuracy is not sufficient to meet the positioning requirement. Further optimization and compensation of the code phase accuracy is required.
The I/Q values for coherently integrating the correlation of the signal for acquisition or tracking are as follows:
similarly, the result of coherent integration on the Q-path is
Wherein A is the signal amplitude, TcohFor the coherent integration time duration, R (τ) is the code autocorrelation function, where τ is the code phase offset, ωeFor angular frequency differences, θeFor frequency-phase difference,. epsilonIIs path I noise, εQD (n) is the satellite navigation message.
The I/Q integration results of 3.5.1 and 3.5.2 are simplified and combined into complex form using the sinc function as:
wherein the content of the first and second substances,
is a frequency error, epsilon
IQIs complex signal noise.
The frequency and code phase of the maximum correlation peak are obtained by fine acquisition and the frequency and code are searched at equal intervals. And selecting the amplitude of a plurality of correlation measurement values with the maximum code phase left-right symmetry under the frequency corresponding to the maximum correlation peak value, and normalizing. Ideally, an autocorrelation function measurement is performed by coherently accumulating millisecond-level intermediate frequency data signals and locally different code phases, thereby determining that one signal and another signal satisfy a trigonometric function. The multi-way correlation curve is shown in fig. 7: the multipath capture phase correlation peak configuration diagram. The 5-way shown in the figure, but not limited to 5-way, can select different numbers according to the fine acquisition code phase interval and the resource requirement.
Analysis of the effect of formula 3.5.3 frequency, sinc (f)
eT
coh) And
the values are equal for relative magnitude scaling of the amplitudes of the multipath correlation peaks. The selected multi-path correlation peak value is the same frequency corresponding to the maximum correlation peak value, and the normalization processing is carried out on the multi-path correlation peak value. With equation 3.5.3, the effect of the frequency difference is eliminated. In addition, for the same satellite signal, the signal strength is equal in a short time (in milliseconds), and then the correlation peak amplitude is only related to the autocorrelation function, which is simplified as:
r(n)=R(τ)+εIQ3.5.4
two compensation values for obtaining an accurate code phase using a multi-path correlation peak, a front-back slope method and a double difference discriminator method will be described below.
Fig. 7 is a schematic diagram of the configuration of the multipath acquisition phase correlation peak R (τ). The code interval based on the fine capture is equal interval, and the interval d is shown as the configuration diagram of the multipath correlation peak.
Front-to-back slope method, as shown in the figure. P is the phase of the maximum correlation peak, E1 and L1 are a pair of narrow-spaced correlator values that are 2d apart, and E2 and L2 are another pair of wider correlator values that are 4d apart. Then the GNSS receiver can measure the slope a of the autocorrelation peak advanced path by using the two advanced correlation values E1 and E2EThe slope a of the self-correlation peak lag posterior is measured using the two lag correlation values L1 and L2LI.e. by
Then, based on the trigonometric characteristics of the autocorrelation curve, the difference in code phase can be calculated according to the following formula:
equation 3.5.6 has the calculation equations for discriminating the code phase error and eliminating the effects of multipath. If there is no multipath, the left and right slopes are equal and correspond to opposite, i.e., aE=-aL. The back-and-forth slope method, however, requires an assumption that the signal power of the direct wave is greater than the signal power of the reflected wave.
Very close to the front-back slope method, the double-difference discriminator method also uses two pairs of leading and lagging correlation values to discriminate the code phase difference after eliminating the multipath influence, i.e. the compensation value, and the code phase difference calculation formula is as follows:
the fact that the trigonometric function of the autocorrelation curve is that of a GNSS C/a code with infinite bandwidth is a trigonometric function. In practice, the if signal is obtained by rf circuit through a band-limiting filter with limited bandwidth to suppress out-of-band noise. There is a flat top problem at the top of the main peak of the autocorrelation function curve and to overcome this problem, the interval d of the correlation peaks is not easily too narrow. The interval d of the correlation peak is determined by the code phase interval of the fine capture, and the large value reduces the search space of the capture, reduces the computation amount and improves the capture speed. Therefore, the value of the interval d of the correlation peak needs to be determined in the discrimination accuracy of the double difference discriminator.
Calculating to obtain the code phase difference deltacpAnd the phase value CodePhase in the millisecond of fine capture to obtain the fraction pseudo range z(k)
z(k)=mod(CodePhase+δcp,1) 3.5.8
Based on the above discussion, the multi-phase code discrimination compensation has the following improvements and features:
1. the input required by the multi-phase code identification is that under the same frequency corresponding to the maximum correlation peak in the fine capture, the front and back symmetrical multi-channel code phase correlation peak values corresponding to the maximum phase are directly obtained, extra calculation is not needed, and the calculation amount is reduced.
2. It has the effect of suppressing multipath from the two discrimination formulas for discriminating code phase errors. Therefore, the effect of multipath inhibition is achieved on the accurate code phase obtained from capturing, the positioning accuracy is improved, and the applicable scene range is enlarged.
3. The main calculation amount of the multi-phase code discrimination is to discriminate the code phase error, and only 3 times of subtraction is done in terms of its formula. Compared with the method that the code phase error is obtained by using curve function curve fitting and carrying out iterative operation, most of operation is saved.
Fig. 8 is a flow chart of a PVT computing device.
Several factors that a GNSS receiver can locate are necessary: the number of positioning satellites is sufficient, ephemeris, accurate time of transmission (TOT) of the satellites, and full pseudoranges that meet the positioning accuracy.
Conventional GNSS receivers obtain accurate satellite transmission time (TOT) by resolving a time of week (TOW) converted word (HOW) given by the navigation bit information. Since an accurate TOT cannot be obtained without decoding of the navigation information. The PVT calculating device obtains a rough position and a corrected TOT time through Doppler least square iteration rough positioning of five state variables supporting multiple modes, and then estimates an accurate TOT by performing six-state least square positioning by using accurate fractional pseudorange. Thus, the PVT computing device does not decode the navigation bit information.
The PVT calculating device mainly has the functions of utilizing the fraction pseudo range (accurate millisecond inner code phase value) and Doppler measured value of the measured value calculating device, correction parameters such as ephemeris, ionosphere and troposphere obtained by the parameter device, rough TOA time obtained by the clock device, and carrying out Least Square (LS) iteration or Kalman filtering positioning calculation under the conditions of no prior position and no TOT, so as to obtain relevant information of the receiver PVT (position, speed, clock error and the like).
The PVT computing device is divided into two main lines which are respectively applied to different scenes, so that the calculated amount is reduced, and the positioning speed and the positioning precision are improved.
First, no prior information (first fix, or power off time too long)
1. Coarse position calculation: the doppler measurements are used to iteratively compute PVT using a Least Squares (LS) supporting multiple modes to obtain a coarse receiver position and a corrected coarse TOA, receiver clock drift.
2. First full pseudorange calculation: and constructing a full pseudorange measurement value by using the position of the coarse receiver, the corrected coarse TOA, the ephemeris information and the fractional pseudorange and adopting the suppression TOT and the position error. The constructed accurate full pseudorange measurements are used for accurate position location.
3. And (3) calculating the accurate position: and iteratively calculating PVT by using the full pseudo-range measurement value and using the Least Square (LS) iteration supporting multimode to obtain the accurate receiver position, the accurate TOT and the receiver clock error.
Second, there is prior information (short-time power-off, intermittent power-off location, continuous location)
1. Second full pseudorange calculation: constructing accurate full pseudorange measurements based on historical position, clock bias, clock drift and time interval
2. Calculating a Kalman fine position: and (4) calculating PVT by using Kalman filtering of the full pseudorange measurement value, and acquiring an accurate receiver position, an accurate TOT and a receiver clock error.
The PVT computing device is applied to different scenes respectively without receiver position information and real-time connection of a wireless network, so that the calculated amount is reduced, the positioning speed and the positioning precision are improved, and the positioning reliability and the robustness are improved.
(1) The five-state variable Doppler positioning and six-state variable precise positioning method supporting multimode is improved. Compared with a single mode, the number of satellites is increased, the positioning geometric factor is improved, and the positioning accuracy is improved. On the other hand, the use scene of positioning is also increased. For example, in an occluded scene, a single mode is very likely to result in an insufficient number of positioning satellites.
(2) Due to the requirement of the LS iterative algorithm for the initial values: the more accurate the initial value is, the faster the iterative convergence speed is, and the lower the possible line of local optimum of the result is. Compared with the traditional method which uses the geocentric as an initial value, the method has the advantages that the average value of the satellite positions projected on the ground is used as the initial value in the rough position calculation, so that the iteration times are reduced, the convergence speed is increased, and the accuracy of the rough position is improved.
(3) And analyzing the characteristic analysis of the estimation error generated by the satellite by the rough position and the rough time error, and selecting the satellite with the highest elevation angle and the satellite with the shortest distance to carry out gradual full pseudorange construction. Therefore, the first full pseudorange calculation has the performance of constructing a full pseudorange measurement value by inhibiting the TOT and the position error under the condition of no prior information. Thereby improving the positioning usability and reliability.
(4) Depending on the differences and characteristics of the usage scenarios, within minutes of the positioning interval. The first full pseudorange calculation is improved by performing TOT estimation using time interval, clock bias, fractional pseudorange, and other information, and then rounding. Thereby estimating an accurate TOT and thus constructing an accurate full pseudorange. Compared with the first full pseudo range calculation, the improved second full pseudo range calculation reduces the calculation amount and improves the calculation speed. And the construction of the full pseudo range satellites is independent, and the estimation errors of other satellites caused by the estimation error of one satellite are avoided.
(5) And improving the multi-mode supporting six-state variable least square positioning calculation during intermittent outage or short-time outage. And using Kalman filtering solution supporting multimode to solve the problem that multiple iterations are needed in the least square positioning solution supporting multimode six state variables. Each iteration requires recalculation of satellite position and velocity information due to TOT updates. And the Kalman filtering solution supporting multimode only needs one-time prediction and updating, thereby further reducing the operation and improving the speed. In addition, compared with a least square algorithm, the filtering algorithm essentially improves the positioning precision and robustness.
The five main steps with and without prior information are described in turn below.
Firstly, rough position calculation: iterative computation of PVT using LS using Doppler measurements
The rough position calculation is to calculate and obtain a rough receiver position, a corrected receiver Time (TOA) and a receiver clock drift by using a least square Doppler positioning solution supporting multimode by taking a Doppler measurement value of a measurement value calculation device, a rough receiving Time (TOA) obtained by transposing a clock and an ephemeris parameter obtained by a parameter device as input. The computed results are used to construct a full pseudorange construction for accurate position location.
Fig. 9 is a flow chart of rough position calculation, which describes the main functions and processes of the implementation process.
1. Obtaining the current input value needed for the calculation: coarse receiver Time (TOA), doppler measurements, ephemeris parameters.
2. The satellite position, velocity, and satellite clock error, clock drift, etc. are calculated using the coarse receiver Time (TOA) and ephemeris parameters. The conventional calculation method is not repeated here.
3. And (5) obtaining the positions of the multimode satellites based on the step 2, and respectively calculating the average value of each direction of XYZ under an ECEF coordinate system. The average is then transformed into the WGS84 coordinate system, the altitude is set to 0, and then transformed into the ECEF coordinate system, the result being calculated as the initial value of the receiver position. In general, the captured satellites are projected onto the ground, distributed around the receiver, and averaged to a value close to the true receiver position, which is closer to the true value than the conventional method using the geocentric position as the initial value. Because the initial value has influence on LS iteration, the local optimal problem is easy to occur when the initial value is too large in deviation, and the iteration result is wrong. This method thus increases the reliability and availability of positioning.
4. And (4) performing five-state parameter least square Doppler positioning iteration supporting multimode based on the initial value of the receiver position in the step (3). The iterative model and process will be described in detail later.
5. And judging whether the iteration converges, outputting the rough receiver position and the corrected TOA by convergence, and exiting otherwise.
Doppler positioning LS iteration for six state variables supporting multiple modes is described in detail below.
In the precise position calculation, a multi-mode least square positioning equation system is supported under the condition that a rough time error exists. See 4.3 for details.4. Here, the coarse pseudorange location equation set is directly introduced: deltaρ=HδX+ε
Wherein:
p is the measured pseudorange (the exact full pseudorange constructed),
are predicted pseudoranges.
For a priori state vector X ═ X
u,y
u,z
u,t
G,t
B,t
C]H is the geometric matrix of the precise position calculation and epsilon is the measurement error of the pseudorange.
The simultaneous differentiation of both sides of 4.3.4 over time yields:
to the left of equation 4.1.2 is the deviant vector of the a priori Doppler measurements
Wherein f is
dIs the measured doppler vector and is the doppler vector,
is the predicted doppler vector.
The first equation to the right of equation 4.1.2 is actually the classical set of equations that measure the speed of the receiver.
Wherein
Is the updated state of the speed of the receiver,
the receiver is updated for the clock drift of the GPS and the clock drift of the BDS respectively,
is the coarse time error rate of change.
The second equation to the right of equation 4.1.2 is a direct relationship between doppler measurements and position.
Based on 4.3.5 and e(k)Is the unit observation vector of the receiver to satellite k,
for H simplification of 4.3.5, the following equation:
the H is brought into the reaction kettle,
to simplify the calculation, the superscript and the corresponding subscript of the navigation system can be omitted, and 4.1.6 also requires calculation
And
and let x be the position of the satellite,
in order to be the velocity of the satellite,
as satellite acceleration, x
xyz_
uIn order to be the position of the receiver,
in order to float the clock of the satellite,
the drift velocity of the satellite clock, namely the acceleration of the satellite clock, is generally small and can be ignored. R | | | x-x
xyz_uAnd | is the geometric distance between the satellite and the receiver.
Since the satellite clock acceleration is generally small and negligible, 4.1.9 is simplified to:
combining 4.1.2, 4.1.3, 4.1.4, 4.1.6, 4.1.10, the multi-mode coarse time doppler positioning equation system can be obtained as follows:
there is a set of equations for instantaneous doppler measurements for 12 state variables (3 receiver positions, 3 receiver velocities, 2 navigation system receiver clock offsets, 2 navigation system receiver clock drifts, coarse timing error rate of change) in the set of equations.
Receiver drift pairs for 2 navigation systems assuming receiver stationarity neglecting coarse time error rateThe same receiver is combined into 1 receiver clock drift variable, the coefficient of the receiver clock difference of 2 navigation systems is 0, and the equation set variable is reduced to 3 receiver positions
Clock drift of receiver
And
same, let it be
Coarse time error
Then 5 unknowns, 5 multimode 5 satellites are needed to obtain the receiver position, and 5 are needed for single mode as well.
Then the 4.1.10 equation is further simplified and the elimination of the ignored variables results in the following equation set:
substitution into 4.1.7
And 4.1.9 obtaining
To obtain
In order to ensure that the water-soluble organic acid,
wherein: f. ofdAre satellite doppler frequency measurements. C is the constant of the speed of light (299792458 m/s). FreqRFIs the carrier frequency of the satellite. GPS has an L1 of 1575.42MHz and BDS has a B1 of 1561.098 MHz.
Based on the above derivation, the doppler positioning equation set is:
and solving an 4.1.16 equation, and iteratively calculating the state quantity of the receiver time, which is corrected by the receiver position, the receiver clock drift and the like through a weighted least square method. The weight W matrix is a diagonal matrix, and the peak-to-average ratio mapping weight based on fine capture is constructed, and the basic idea is that the larger the peak-to-average ratio is, the higher the weight is.
The calculation steps are as follows:
1. and setting the initial value of the position of the receiver in the mean coordinate projected by the satellite, wherein the initial value of the clock drift of the receiver is 0, and the rough receiving Time (TOA) obtained by the clock transposition is used as the rough time error initial value. The doppler measurements use the doppler measurements of the computing device.
2. The pseudo-range is calculated from the coarse position, satellite coordinates and satellite velocity in the second half of 4.1.15
3. From predicted pseudoranges
And Doppler measurement value f
dDeviation of (2)
Calculating state variable update values
Calculate the updated formula as
4. D, calculating the delta obtained in the step 3
χUpdate to a state variable
5. And judging whether the iteration convergence condition is met, and quitting if the iteration convergence condition is met. The jump in step 2 is not satisfied. The judgment condition of whether the iteration converges is
Time, indicates convergence.
Second, calculating a first full pseudorange: construction of full pseudorange measurements with suppressed TOT and position error
For GNSS receiver positioning, full pseudoranges are a necessary factor for positioning. The coarse position calculation only obtains one coarse position and a coarse corrected TOA, and the measurement calculation only obtains precise fractional pseudoranges within milliseconds. The full pseudorange needs to be constructed. How full pseudoranges are obtained based on the coarse corrected TOA, the coarse position, and the fractional pseudorange calculations is described below.
In a traditional GNSS receiver, accurate intra-week seconds are obtained by decoding navigation data after bit synchronization and frame synchronization, and then the accurate intra-week seconds are determined by adding the number of local tracking C/A codes on the basis of the intra-week seconds. And the C/a code repeats every 1 ms, so the fractional pseudoranges represent the fraction of a segment within 1 ms.
Let the system have N GPS satellites, M BDS satellites, and K ═ M + N. Constructed full pseudo range rho and measurement fraction pseudo range z(k)The relation of (A) is as follows:
we know the measured pseudoranges of a conventional receiver, expressed as:
wherein TOA is the TOA of the receiver signal relative to the receiver clock, TOT(k)Is the time of transmission TOT relative to the satellite clock, C is the speed of light constant.
If a conventional GNSS receiver does not have time information at the start, the TOA may be estimated based on the T0T value relative to the satellite plus some common offset (receiver clock error). This initial TOA will have an absolute error of about tens of milliseconds. If the receiver has been initialized, T0A is typically within 1 millisecond.
Then, the relationship between the measured pseudorange and the true signal travel distance is:
wherein the content of the first and second substances,
and
the geometrical range of the satellites and receivers of the GPS and BDS respectively,
and
satellite clock offsets, t, for GPS and BDS at TOT time, respectively
GAnd t
BRelative GPS and BDS clock offsets for the receiver, respectively, then the TOT can be estimated from 4.2.1 and 4.2.2, as follows:
full pseudoranges may be constructed based on the relationships of 4.2.1 and 4.2.4, with T0T constructed for each satellite based on 4.24 without decoding the satellite navigation data. T0T is constructed with respect to the TOA of an estimated GNSS signal, so it differs from the true T0T by an integer number of milliseconds, which does not achieve the accuracy of the clock used to calculate the satellite position, velocity, acceleration and ephemeris data. As we know, the common deviation has no influence on positioning, and the TOT accurate value is obtained by iteration calculation as a state variable in precise position calculation.
In the following, the subscripts of the corresponding systems of GPS and BDS are omitted for simplicity of the formula.
The principle of how to construct full pseudoranges with correct common bias is described below.
First, an integer millisecond value N of a selected reference satellite is calculated(0)So construct the full pseudorange as N(0)+z(0). (to distinguish the reference star from the star to be calculated, the reference star is labeled 0).
Wherein: r is
(0)Is the true geometric distance of the satellite and the receiver,
is the geometric distance, d, estimated from a priori coarse time (time of transmission) and a priori position
(0)Is that
The error present in (a), caused by a priori time and a priori position; t is t
(0)Is the clock error of the satellite clock (coarse time ephemeris acquisition), t
uIs the common deviation (clock difference between receiver and satellite), ε
(0)Is measurement error (including atmospheric error and thermal noise).
The main difficulty today is to assign an integer value to all measurements so that they have the same common deviation. Assuming this has been done successfully, there will be the following equation for satellite k:
now subtracting the two equations, one can see how N is constructed(k)So that b has the same value for each satellite.
Then:
the following two conclusions can be drawn from the above formula:
1. common deviation tuIs correctly eliminated.
2. Not knowing d(0)And d(k)It is the values of (c) that need to be found because their values depend on a priori position and time. However, if-d(k)+ε(k)+d(0)-ε(0)Less than 0.5ms, the correct N will be obtained(k)The calculation is as follows.
Wherein
z
(0)And z
(k)For the measured millisecond pseudoranges,
and delta
tIs obtained from ephemeris in the prior state.
Thus, no matter N is selected(0)Is what value of (1) is obtained as N(k)Implicit common deviation tuConsistent for all satellites.
The effect of TOT coarse time error and coarse position error on constructing full pseudoranges is analyzed as follows and an efficient method of suppressing each error is proposed.
1. Initial reference satellite selection (highest elevation satellite)
In calculating N(k)When it is indicated that if-d(k)+ε(k)+d(0)-ε(0)Is less than 0.5ms (150km), the correct N will be obtained(k). Then d(k)Is the geometric range estimation error of satellite k, which is mainly caused by a priori position and a priori time error. And due to measurement error epsilon(k)Much less than d(k)Therefore, the constraint is considered to be | d(k)-d(0)|<0.5ms。
The position and coarse time errors affect the satellites as follows:
1. the prior position error is 1-2 orders of magnitude better for the vertical component than the horizontal component
2. The prior coarse time error has a large effect on the satellite with the maximum range rate, namely, the satellite ascending and descending.
Thus, the initial reference satellite selects the satellite with the highest elevation angle. The effect of errors on both the a priori position and the a priori time is minimized.
2. Reference satellite selection (selecting the closest satellite to the historical reference satellite)
If the limits of the a priori errors are to be approached, then it is noted that our algorithm for computing integer values is effective because the common bias is kept consistent for each satellite. The reference satellite can be replaced arbitrarily but still maintain consistency against common errors. In this case, | d can be preferentially analyzed(k)-d(0)The smallest combination. Because the prior position error has similar influence on the satellites in the same area in the sky, and the prior time error also has similar influence on the satellites with similar distance change rates. Therefore, N of the second satellite can be calculated after the initial reference satellite is selected(k)The value is obtained. This second satellite may then be referred to as the reference satellite and N for a third satellite may be calculated(k)The value is obtained. Each time the satellite with the shortest distance to the reference satellite is selected. This way d is minimized in each step(k)-d(0)|。
In conjunction with the above analytical discussion, the specific steps of the full pseudorange calculation, as shown in fig. 10, are given as follows:
1. the satellite position at the corresponding time is calculated using the corrected TOA obtained from the coarse positioning.
2. Calculating an estimated geometric distance between the satellite and the receiver using the satellite position and the coarse receiver position
3. Selecting an initial reference satellite (the satellite with the highest elevation angle) based on an initial reference satellite selection principle, and calculating the whole millisecond number of the reference satellite
4. Selecting the satellite with the shortest geometric distance to the reference satellite to calculate an integer of milliseconds:
5. constructing full pseudoranges ρ for all satellites(k)=N(k)+z(k)。
Thirdly, calculating the accurate position: iterative computation of PVT using LS using full pseudorange measurements
The precise position calculation: and adopting the constructed full pseudorange measurement value, the ephemeris parameter and the corrected TOA as an input quantity and an initial value, and adopting the least square iteration supporting multimode to calculate the PVT so as to obtain the accurate receiver position, the accurate receiver receiving Time (TOA) and the clock error.
Compared with the conventional positioning equation, there is an accurate time of transmission (TOT), and here, the corrected TOT is obtained through the previous steps, which is actually not accurate. Calculating accurate position by using the constructed full pseudo-range measurement value and ephemeris parameter and introducing coarse time error deltatc。
How to introduce coarse timing errors and obtain accurate position is described further below.
It is known that the common bias affects the actual measurement values and affects each satellite by the same amount, without affecting the positioning. The TOT error influences the predicted value, and the influences are different in magnitude. To solve the coarse time error positioning problem, the pseudorange rate can be calculated by knowing the relative velocity and clock frequency of each satellite. Such a coarse timing error can be expressed as:
wherein, delta
toaRepresenting an updated value for the a priori coarse-time state quantity. t is t
txRepresenting the actual satellite transmission time.
Represents a pair t
txCoarse time estimation of (2).
Is the pseudorange rate. e.g. of the type
(k)Is the unit observation vector from the receiver to satellite k,
R
kis the geometric distance of the satellite to the receiver.
Is the satellite velocity vector.
Representing the satellite clock rate.
Are predicted pseudoranges.
Let state vector X be [ X ]u,yu,zu,tG,tB,tC]Wherein [ x ]u,yu,zu]As a receiver position variable, tGFor the receiver GPS clock error, tBFor the BDS clock difference of the receiver, tCA coarse time state quantity of the time TOA is received for the prior receiver.
In order to ensure that the water-soluble organic acid,
an update vector to a prior state.
Therefore, it is necessary to establish for each satellite
And
the relationship between them.
Middle acceptor
The only part of the influence is
And is
And
the relationship therebetween was found to be 4.3.1. Thus for each satellite
And X is:
in correspondence with the GPS, the GPS receiver is,
in correspondence with the BDS,
wherein epsilon
(k)For pseudorange measurement noise, p
(k)To measure pseudoranges (the exact full pseudoranges constructed),
are predicted pseudoranges.
Then 4.3.2 and 2.3.3 are taken together, for K available satellites, N for GPS and M for BDS, and K equals M + N. A system of matrix equations may be obtained.
δρ=HδX+ε 4.3.4
Wherein
The meaning of the variables of 4.3.5 is further explained in connection with the observation vector diagram of the receiver pointing to the satellite of fig. 11.
Wherein the subscript B represents the BDS system and G represents the GPS system. The following description omits the subscript description. R is the geometric distance of the corresponding system satellite to the receiver position. [ x ] of
u,y
u,z
u]Respectively, represent the position of the receiver in the ECEF coordinate system.
Are pseudorange rates representing the corresponding system satellites. e.g. of the type
(k)Is the unit observation vector from the receiver to satellite k,
representing corresponding system satellite clock rates
Resolving 4.3.4 equation, and iteratively calculating [ x ] by weighted least squares methodu,yu,zu]For receiver position, GPS clock error tGBDS clock difference tBAccurate receiver time tCAnd the value of the variable is equal. The weight W matrix is a diagonal matrix, and the peak-to-average ratio mapping weight based on fine capture is constructed, and the basic idea is that the larger the peak-to-average ratio is, the higher the weight is.
The calculation steps are as follows:
1. the coarse receiver position obtained by the coarse position calculation, the corrected receiver Time (TOA), the GPS clock offset and the BDS clock offset are 0 as iteration initial variables. And constructing accurate full pseudoranges as pseudorange measurements.
2. Pseudoranges predicted from coarse position and satellite coordinates
3. From predicted pseudoranges
Deviation delta from constructed full pseudorange measurement rho
ρCalculating the state variable update value
The calculation update formula is delta
X=(H
TWH)
-1H
TWδ
ρ
4. D, calculating the delta obtained in the step 3XUpdate to state variable X ═ Xu,yu,zu,tG,tB,tC]
5. And judging whether the iteration convergence condition is met, and quitting if the iteration convergence condition is met. The jump in step 2 is not satisfied.
Fourthly, calculating a second full pseudorange: constructing full pseudorange measurements based on historical position, clock bias, clock drift and time interval
Through the foregoing analysis, the first full pseudorange calculation is constructed based on the coarse TOA time, which performs a more complicated satellite selection process and calculates an estimation error of one satellite to cause an estimation error of other satellites in order to suppress the coarse time error and the coarse position error. In response to these problems, it is further improved.
When the position has been located, accurate historical TOA time, receiver clock error, receiver clock drift, and historical receiver position information are obtained. According to the differences and characteristics of the use scenarios, the error is generally less than 0.5ms caused by constructing the full pseudorange through the TOA, the receiver clock, and the receiver position within a few minutes of the positioning interval time. That is, receiver position motion over a short period of time, errors in the clock error estimate, and changes in pseudorange due to satellite position errors caused by clock error estimate errors directly affect the supported gap time interval corrections. Compared with the first full pseudo range calculation, the improved second full pseudo range calculation reduces the calculation amount and improves the calculation speed. And the construction of the full pseudo range satellites is independent, and the estimation errors of other satellites caused by the estimation error of one satellite are avoided.
As shown in fig. 12, the calculation step of the second full pseudorange calculation:
1. and estimating the current receiver clock difference by using the historical receiver clock difference and the historical receiver clock drift time interval T. Wherein:
in the case of the GPS system, it is,
in the case of the BDS system,
2. estimated current receiver clock error
The time interval T of the historical accurate TOA and local estimate estimates the time of transmission of the satellite, which estimates have errors on the order of milliseconds, but which affect the estimated satellite position and satellite clock error by far less than 0.5 ms.
In the case of the GPS system, it is,
in the case of the BDS system,
wherein. The GPS is N satellites, and N represents the nth satellite. Similarly, the BDS is M satellites, and M represents the M-th satellite.
And
historical satellite-to-receiver geometric distances.
3. Computing
And
the time of day satellite position and the clock error of the satellite.
4. Calculating an estimated predicted geometric distance between the satellite and the receiver
And
5. calculate the integer number of milliseconds for each satellite:
in the case of the GPS system, it is,
in the case of the BDS system,
6. constructing full pseudoranges ρ for all satellites(k)=N(k)+z(k)。
Fifthly, Kalman fine position filtering calculation: calculation of PVT using full pseudorange measurement Kalman filtering
And when intermittent power failure or short-time power failure occurs, accurate pseudorange and TOT can be obtained through second full pseudorange calculation. Positioning calculation can be carried out through a traditional multi-mode Least Square (LS) pair, and operations of recalculating satellite positions and speeds and the like in each iterative operation are reduced. However, since the least square positioning result does not relate the motion characteristics at different time instants, the least square result is usually rather rough and cluttered, thereby affecting the positioning accuracy and robustness.
The invention improves the positioning precision and robustness and reduces the operation by supporting the multi-mode Kalman filtering solution. The Kalman filtering solution supporting multi-mode only needs one prediction and update, thereby further reducing operation and improving speed.
The model and the updating method are described in detail as follows.
The basic idea of the kalman filter is to use two equations to represent the unknown state transition process and the measurement system input and output relationships, respectively, to relate the state value at a certain time to the current and previous time measurements. The essence of the kalman filtering problem is therefore to solve the problem of unknown state equations and measurement equations jointly in some optimal way. Discrete equations of state and measurement are established as in equations 4.5.1 and 4.5.2.
s(k)=F*s(k-1)+w(k-1) (4.5.1)
z(k)=G*s(k)+v(k) (4.5.2)
Where s (k) and z (k) are the state variables and measurement vectors at time k, respectively, F and G are the state transition matrix and measurement matrix, respectively, and w (k) and v (k) are the process noise and measurement error, respectively, at time k.
The iterative process of the kalman filter is shown in equations (4.5.3) to (4.5.8) and fig. 13.
The prediction process, also called time update process, uses the state equation of the system to predict the state value of the current (i.e. k-th) epoch on the basis of the state estimation value of the last (i.e. k-1) epoch.
Update procedure-measurement update procedure, also known as correction procedure, uses actual measured values to correct the state prior estimates obtained by the prediction procedure.
Wherein the content of the first and second substances,
and
respectively representing the predicted value and the estimated value of the state variable at the k-th moment,
and P
kCovariance matrices of prediction and estimation errors at time k, Q and R covariance matrices of w (k) and v (k), respectively, e
kIs a measurement vector z
kCorresponding innovation process, K
kIs the kalman gain at time k. It can be seen from the above solving process that the solution of the kalman filter is calculated recursively, and each updated estimate of the state is calculated from the previous estimate and new input data, so that only the previous estimate needs to be stored and real-time processing can be achieved.
wherein:
process noise covariance matrix:
wherein the content of the first and second substances,
wherein the content of the first and second substances,
is the maximum value of the variance of the acceleration in each direction,
t is a time interval.
System measurement matrix:
wherein:
wherein,
And
is the geometric distance, I, from the corresponding satellite to the receiver
GAnd I
BIs an identity matrix.