CN109725290A

CN109725290A - A kind of error extracting method, device, electronic equipment and readable storage medium storing program for executing

Info

Publication number: CN109725290A
Application number: CN201910096259.3A
Authority: CN
Inventors: 邓中亮; 韩可; 汤灿阳; 胡恩文; 朱棣; 唐诗浩; 刘延旭; 綦航; 赵鹤
Original assignee: Beijing University of Posts and Telecommunications
Current assignee: Beijing University of Posts and Telecommunications
Priority date: 2019-01-31
Filing date: 2019-01-31
Publication date: 2019-05-07
Anticipated expiration: 2039-01-31
Also published as: CN109725290B

Abstract

The embodiment of the invention provides a kind of error extracting method, device, electronic equipment and readable storage medium storing program for executing, applied to wireless location technology field, the described method includes: obtaining the reception signal of base station and observation station in each day preset time period respectively, difference processing is carried out to the pseudorange and carrier phase that receive signal, according to the position coordinates of obtained double difference pseudorange, double difference carrier phase and base station, observation station is obtained in the difference coordinate sequence in each day；Difference coordinate sequence is handled by two-dimensional movement weighted mean method, obtain smooth coordinate sequence, by Wavelet Packet Algorithm to it is any with reference to day and with this with reference to the maximum neighbouring day of day cross-correlation coefficient smooth coordinate sequence carry out M layers decompose, judge whether the destination layer J in each decomposition layer in the cross-correlation coefficient of low frequency part where maximum value is less than M；If so, the low frequency part extracted in destination layer J carries out wavelet reconstruction, Multipath Errors are obtained.The present invention improves the accuracy of error extraction.

Description

A kind of error extracting method, device, electronic equipment and readable storage medium storing program for executing

Technical field

The present invention relates to wireless location technology field, more particularly to a kind of error extracting method, device, electronic equipment and Readable storage medium storing program for executing.

Background technique

In short baseline high accuracy positioning, can be eliminated or be weakened by differential technique ionosphere delay, tropospheric delay, Satellite orbital error, receiver clock-offsets and satellite clock correction equal error, but multipath effect does not have correlation, nothing at baseline both ends Method is eliminated by difference, therefore Multipath Errors become the main error source of high accuracy positioning.

In order to detect, inhibit and eliminate Multipath Errors, related scholar designs from satellite-signal, the design of receiving antenna with Addressing, Digital Signal Processing and location navigation calculating etc. propose a variety of strategies.Wherein, location navigation is calculated without penetrating Any change is done to the processing for receiving signal in frequency front end and baseband digital signal processing module, only need to be to data after survey at Reason.Data processing method includes elevation angle weighting technique, Extended Kalman filter technology, utilizes signal-to-noise ratio skill after the survey proposed Art, empirical mode decomposition technology, Vondark filtering technique and wavelet technique etc..

Wherein, wavelet threshold denoising technology has a kind of " collection to determining signal using the mutative scale characteristic in wavelet transformation In " ability.If the energy of a signal concentrates on wavelet transformed domain minority wavelet coefficient, their value is naturally larger than In a large amount of signals of wavelet transformed domain self-energy dispersion and the wavelet coefficient of noise, the method for at this moment available threshold denoising.It is given One threshold value δ, wavelet coefficient of all absolute values less than δ turn to " noise ", their value is replaced with 0；And it is more than the small of threshold value δ The numerical value of wave system number, value again again after being contracted by.After wavelet reconstruction, required signal is obtained.But become using small echo The multiresolution analysis changed is that the multipath for decomposing, therefore extracting to the low frequency part of signal only has the more of low frequency part Path, and the Multipath Errors of signal high frequency section do not extract.In the Multipath Errors using Wavelet transformation to signal When carrying out multiresolution analysis, the influence of observation noise is had ignored, the effect is unsatisfactory so that signal is filtered.Therefore, existing Wavelet technique in, the accuracy of the Multipath Errors of extraction is lower, and filtered signal accuracy is lower.

Summary of the invention

The embodiment of the present invention is designed to provide a kind of error extracting method, device, electronic equipment and readable storage medium Matter, to improve the accuracy of error extraction.Specific technical solution is as follows:

The embodiment of the invention provides a kind of error extracting methods, which comprises

The reception signal of base station and observation station in each day preset time period is obtained respectively, for the base station and institute The daily reception signal in observation station is stated, the reception of each epoch in the preset time period to the base station and the observation station The pseudorange and carrier phase of signal carry out difference processing, obtain double difference pseudorange and double difference carrier phase；

Position resolving is carried out according to the position coordinates of the base station, the double difference pseudorange in each day and double difference carrier phase, is obtained Difference coordinate sequence to the observation station in each day；

Difference coordinate sequence by two-dimensional movement weighted mean method to the observation station in each day is handled, and obtains institute Observation station is stated in the smooth coordinate sequence in each day, the two-dimensional movement weighted mean method is based on time domain space and Value space carries out Moving weighted average；

For it is any refer to day, calculate the observation station this with reference to day and other adjacent to day smooth coordinate sequence it is mutual Relationship number chooses the smooth coordinate sequence that cross-correlation coefficient maximum neighbouring day in day is referred to reference to day and with this；

The smooth seat in the maximum neighbouring day of day cross-correlation coefficient is referred to this with reference to day and with this by Wavelet Packet Algorithm It marks sequence and carries out M layers of decomposition, calculate the cross-correlation coefficient of low frequency part in each decomposition layer, M is pre-set Decomposition order, is sentenced Whether disconnected destination layer J is less than M, and the destination layer J is maximum value institute in the cross-correlation coefficient of low frequency part in each decomposition layer The number of plies；

If so, the low frequency part extracted in the destination layer J carries out wavelet reconstruction, Multipath Errors are obtained；

If not, the value of M is added 1, return described mutual with reference to day with reference to day and with this to this by Wavelet Packet Algorithm The smooth coordinate sequence in the maximum neighbouring day of relationship number carries out the step of M layers of decomposition, until the destination layer J is less than M.

Optionally, the difference coordinate sequence by two-dimensional movement weighted mean method to the observation station in each day carries out Processing, obtains the observation station in the smooth coordinate sequence in each day, comprising:

If the number of the sample epoch coordinate is I, determine that the rolling average cycle T of the time domain space is no more than I Odd number；

If i-th of sample epoch coordinate is x_i, Indicate Value space, i is the integer of 1~I, according to formula:Calculate the average value of each sample epoch coordinate

According to formula:Calculate the standard deviation sigma of each sample epoch coordinate_I；

By the standard deviation sigma_IAs the threshold θ for choosing Value space epoch coordinate_I, then θ_I=σ_I,

If being smoothed to r-th of sample epoch coordinate, r is the integer of 1~I；

According to formula:Calculate each sample epoch coordinate and current sample epoch in Value space Coordinate x_rStandard deviation sigma_r,

By standard deviation sigma_rNo more than threshold θ_ISample epoch coordinate as fraction epoch coordinate x_k, the fraction epoch seat Mark x_kFor determining in r-th of smoothed out value of sample epoch coordinate of Value space, fraction epoch coordinate x_kNumber be r-th Rolling average period N of the sample epoch coordinate in Value space_r, N_r≤I；For r-th of sample epoch coordinate, if α_jWhen expression The weighting coefficient of j-th of domain space delay epoch coordinate, j areInteger, [] indicate round numbers,

β_kIndicate that the weighting coefficient of k-th of fraction epoch coordinate of Value space, k are 1~N_rInteger,

According to formula:Calculate each sample epoch in the rolling average cycle T of the time domain space Coordinate is relative to the average retardation epoch coordinate of r-th of sample epoch coordinateΔτ_jIt is opposite for j-th of sample epoch coordinate In r-th of sample epoch coordinate x_rDelay epoch,

According to formula:It calculates time domain space and respectively postpones the dissimilarity coefficient of epoch coordinate and be According to formula:It calculates j-th of time domain space and prolongs The weighting coefficient α of slow epoch coordinate_j；

According to formula:Calculate fraction epoch coordinate x in Value space_kAverage valueΔμ_kFor kth Difference of a fraction epoch coordinate relative to r-th of sample epoch coordinate；

According to formula:Calculate the dissimilarity coefficient of k-th of fraction epoch coordinate of Value space

According to formula:Calculate the weighting coefficient β of k-th of fraction epoch coordinate of Value space_k；

According to formula:

By r-th of sample epoch coordinate x_rRespectively in time domain space and Value space with x_rCentered on to carry out mobile weighting flat It handles, obtains the time domain space moving weighted average value of r-th of sample epoch coordinateAnd Value space rolling average adds Weight

According to formula:Calculate x_rAverage value under two-dimensional space

According to formula:Calculate the dissimilarity coefficient of time domain space

According to formula:Calculate the dissimilarity coefficient of Value space

According to formula:Calculate the weighting coefficient γ of time domain space_T, according to formula:Calculate the weighting coefficient of Value space

According to formula:It willWithIt is weighted processing, is obtained by smoothly locating Coordinate after reason

Optionally, the reception signal to the base station and the observation station each epoch in the preset time period Pseudorange and carrier phase carry out difference processing, obtain double difference pseudorange and double difference carrier phase, comprising:

According to formula:

Meter Observation station r and base station b is calculated in the double difference pseudorange of same epoch observation satellite n and satellite mAnd observation station r and Double difference carrier phase of the base station b in same epoch observation satellite n and satellite m

Wherein,For double difference operator, S indicates that GPS system, satellite n and satellite m are the satellite in the GPS system,Indicate observation station r and base station b in same epoch observation satellite n and satellite m satellite hub to receiver phase Geometric distance between the center of position, λ is carrier wavelength；It is observation station r and base station b in same epoch observation satellite The unknown integer ambiguity of the carrier phase of n and satellite m,Indicate that base station b and observation station r is seen in same epoch The pseudorange multipath error of satellite n and satellite m is surveyed,Indicate that observation station r and base station b is observed in same epoch The carrier phase Multipath Errors of satellite n and satellite m,Indicate that observation station r and base station b is observed in same epoch The pseudo-code of satellite n and satellite m,Indicate observation station r and base station b same epoch observation satellite n's and satellite m Carrier phase measurement noise.

Optionally, it is described by Wavelet Packet Algorithm to this with reference to day and with this with reference to day cross-correlation coefficient it is maximum neighbouring It smooth coordinate sequence carries out M layers of decomposition, comprising:

If this is x with reference to day and with the smooth coordinate sequence with reference to the maximum neighbouring day of day cross-correlation coefficient_i, according to Formula:

To x_iWAVELET PACKET DECOMPOSITION is carried out,

I is the integer of 1~I, and ω is the integer of 1~2M, and q ∈ I, q are with reference to epoch；

Wavelet packet when indicating undecomposed,Indicate the 2nd ω -1 wavelet packet coefficients on M layer,Indicate the ω wavelet packet coefficient on M-1 layer,Indicate the 2nd ω wavelet packet system on M layer Number, G are scaling function resolution filter, and H is wavelet function resolution filter,

The low frequency part extracted in the destination layer J carries out wavelet reconstruction, obtains Multipath Errors, comprising:

According to formula: To x_iCarry out wavelet package reconstruction；

Indicate the 2nd ω -1 wavelet packet coefficients on J+1 layer,Indicate the 2nd on J+1 layer ω wavelet packet coefficient, g are scaling function reconfigurable filter, and h is wavelet function reconfigurable filter.

The embodiment of the invention provides a kind of error extraction element, described device includes:

Double difference computing module, for obtaining the reception signal of base station and observation station in each day preset time period respectively, For the daily reception signal of the base station and the observation station, to the base station and the observation station when described default Between in section the pseudorange of the reception signal of each epoch and carrier phase carry out difference processing, obtain double difference pseudorange and double difference carrier wave phase Position；

Difference coordinate sequence computing module, for according to the position coordinates of the base station, the double difference pseudorange in each day and double Poor carrier phase carries out position resolving, obtains the observation station in the difference coordinate sequence in each day；

Smoothing module, for passing through difference coordinate sequence of the two-dimensional movement weighted mean method to the observation station in each day Column are handled, and obtain the observation station in the smooth coordinate sequence in each day, the two-dimensional movement weighted mean method is based on time domain Space and Value space carry out moving weighted average；

Smoothly coordinate sequence chooses module, and for referring to day for any, calculating the observation station, this refers to day and other The cross-correlation coefficient of the smooth coordinate sequence in neighbouring day chooses this with reference to day and with this with reference to the maximum neighbour of day cross-correlation coefficient The smooth coordinate sequence in nearly day；

Decomposing module, for referring to the maximum neighbour of day cross-correlation coefficient by Wavelet Packet Algorithm to this with reference to day and with this The smooth coordinate sequence in nearly day carries out M layers of decomposition, calculates the cross-correlation coefficient of low frequency part in each decomposition layer, and M is to preset Decomposition order；

Judgment module, for judging whether destination layer J is less than M, the destination layer J is low frequency part in each decomposition layer Cross-correlation coefficient in the number of plies where maximum value；

Reconstructed module, for when the judgment result of the judgment module is yes, extracting the low frequency portion in the destination layer J Divide and carry out wavelet reconstruction, obtains Multipath Errors；

When loop module for the judging result in the judgment module is no, the value of M is added 1, is returned described by small Wave packet algorithm carries out M layers points to this with reference to day and with the smooth coordinate sequence with reference to the maximum neighbouring day of day cross-correlation coefficient The step of solution, until the destination layer J is less than M.

Optionally, the smoothing module is specifically used for, if the number of the sample epoch coordinate be I, determine described in The rolling average cycle T of time domain space is the odd number no more than I；

According to formula:

According to formula:Calculate x_rAverage value under two-dimensional space

According to formula:Calculate the dissimilarity coefficient of Value space

Optionally, the double difference computing module is specifically used for according to formula:

Wherein,For double difference operator, S indicates that GPS system, satellite n and satellite m are the satellite in the GPS system,Indicate observation station r and base station b in same epoch observation satellite n and satellite m satellite hub to receiver phase Geometric distance between the center of position, λ is carrier wavelength,It is observation station r and base station b in same epoch observation satellite The unknown integer ambiguity of the carrier phase of n and satellite m,Indicate that base station b and observation station r is seen in same epoch The pseudorange multipath error of satellite n and satellite m is surveyed,Indicate that observation station r and base station b is observed in same epoch The carrier phase Multipath Errors of satellite n and satellite m,Indicate that observation station r and base station b is observed in same epoch The pseudo-code of satellite n and satellite m,Indicate observation station r and base station b same epoch observation satellite n's and satellite m Carrier phase measurement noise.

Optionally, the decomposing module, if maximum with reference to day cross-correlation coefficient with reference to day and with this specifically for this The smooth coordinate sequence in neighbouring day is x_i, according to formula:

To x_iWAVELET PACKET DECOMPOSITION is carried out,

The reconstructed module is specifically used for according to formula:

To x_iIt carries out Wavelet package reconstruction；

The embodiment of the invention provides a kind of electronic equipment, comprising: processor, communication interface, memory and communication bus, Wherein, the processor, the communication interface, the memory complete mutual communication by the communication bus；

The memory, for storing computer program；

The processor when for executing the program stored on the memory, realizes any of the above-described error The step of extracting method.

The embodiment of the invention provides a kind of computer readable storage medium, storage in the computer readable storage medium There is computer program, when the computer program is executed by processor, realizes the step of any of the above-described error extracting method Suddenly.

Error extracting method, device, electronic equipment and readable storage medium storing program for executing provided in an embodiment of the present invention, obtain base respectively Quasi- station and reception signal of the observation station in each day preset time period, it is right for the daily reception signal of base station and observation station The pseudorange of the reception signal of each epoch and carrier phase carry out difference processing, and root within a preset period of time for base station and observation station Position resolving is carried out according to the position coordinates of the double difference pseudorange in obtained each day, double difference carrier phase and base station, is observed It stands in the difference coordinate sequence in each day；Difference coordinate sequence is handled by two-dimensional movement weighted mean method, is observed It stands in the smooth coordinate sequence in each day, two-dimensional movement weighted mean method is based on time domain space and Value space carries out mobile weighting and puts down ?；M layers of decomposition are carried out by the Wavelet Packet Algorithm smooth coordinate sequence with reference to day and neighbouring day maximum to cross-correlation coefficient, The cross-correlation coefficient for calculating low frequency part in each decomposition layer, judges whether destination layer J is less than M, and destination layer J is low in each decomposition layer The number of plies in the cross-correlation coefficient of frequency part where maximum value；If so, the low frequency part extracted in destination layer J carries out small echo weight Structure obtains Multipath Errors；If not, the value of M is added 1, Wavelet Packet Algorithm is passed back through to the maximum reference of cross-correlation coefficient The step of smooth coordinate sequence in it and neighbouring day carries out M layers of decomposition, until destination layer J is less than M；It extracts in destination layer J Low frequency part carries out wavelet reconstruction, obtains Multipath Errors.It, can be with by two-dimensional movement weighted mean method in the embodiment of the present invention Weaken influence of the low frequency part observation noise to signal, and solve the end effect in signal smoothing processing, reaches and protect side denoising Purpose.Wavelet Packet Algorithm can be divided signal band at many levels, further be decomposed to high frequency section, and radio-frequency head is eliminated Observation noise in point, to improve the accuracy of Multipath Errors extraction.Certainly, it implements any of the products of the present invention or method It does not necessarily require achieving all the advantages described above at the same time.

Detailed description of the invention

In order to more clearly explain the embodiment of the invention or the technical proposal in the existing technology, to embodiment or will show below There is attached drawing needed in technical description to be briefly described, it should be apparent that, the accompanying drawings in the following description is only this Some embodiments of invention for those of ordinary skill in the art without creative efforts, can be with It obtains other drawings based on these drawings.

Fig. 1 is the flow chart of the error extracting method of the embodiment of the present invention；

Fig. 2 is the WAVELET PACKET DECOMPOSITION tree structure diagram in the embodiment of the present invention after 3 layers of decomposition；

Fig. 3 is to put down by wavelet threshold denoising algorithm, two-dimensional movement Weighted Average Algorithm and based on two-dimensional movement weighting The Wavelet Packet Algorithm handled is filtered the signal after denoising；

Fig. 4 is to put down by wavelet threshold denoising algorithm, two-dimensional movement Weighted Average Algorithm and based on two-dimensional movement weighting The Wavelet Packet Algorithm handled be filtered after residual signals；

Fig. 5 is the structure chart of the error extraction element of the embodiment of the present invention；

Fig. 6 is the structure chart of the electronic equipment of the embodiment of the present invention.

Specific embodiment

Following will be combined with the drawings in the embodiments of the present invention, and technical solution in the embodiment of the present invention carries out clear, complete Site preparation description, it is clear that described embodiments are only a part of the embodiments of the present invention, instead of all the embodiments.It is based on Embodiment in the present invention, it is obtained by those of ordinary skill in the art without making creative efforts every other Embodiment shall fall within the protection scope of the present invention.

When carrying out deformation observation using high-precision location technique, since observation station is generally fixed and ambient enviroment is protected substantially Hold constant, GPS (Global Positioning System, global positioning system) satellite and observation station and peripheral reflection face Geometric configuration there is repeatability, with a sidereal day (11 hours 58 points) for the period.There are multipath effect, The characteristic of its repeatability allows for closing on the coordinate sequence between day, and there are correlations.Based on this characteristic, the embodiment of the present invention is mentioned A kind of error extracting method, device, electronic equipment and readable storage medium storing program for executing have been supplied, to improve the accuracy of error extraction, thus Weaken multipath effect.

Error extracting method is provided for the embodiments of the invention first below to describe in detail.

It is the flow chart of the error extracting method of the embodiment of the present invention referring to Fig. 1, Fig. 1, comprising the following steps:

S101 obtains reception signal in each day preset time period of base station and observation station respectively, for base station and The daily reception signal in observation station, to the pseudorange and load of base station and observation station the reception signal of each epoch within a preset period of time Wave phase carries out difference processing, obtains double difference pseudorange and double difference carrier phase.

In the embodiment of the present invention, since observation station antenna will lead to multichannel from small variation with reception signal reflex identity distance Diameter effect changes a lot, and this variation shows as periodically changing, and will make the duplicate period hair of multipath effect Changing, therefore observation station changes in range as small as possible and tends to be constant.It is more when the incidence angle for reflecting signal is excessive The repeatability of path effects, which will be lower, even to disappear.Therefore, the incidence angle for reflecting signal is small as far as possible and less than 15 °.

The position coordinates of base station be it is determining, the position coordinates of observation station be it is unknown, base station and sight can be passed through The reception signal of survey station determines the position of observation station.It, can be with specifically, due to closing on the coordinate sequence between day there are correlation The reception signal of base station and observation station in each day preset time period is obtained, when preset time period can be arbitrary in one day Between section, such as 8:00-12:00, it is not limited here.For the reception letter of base station and observation station daily within a preset period of time Number, preset time period can be divided into multiple epoch, epoch is that some astronomical parameters are used as reference at the time of point in astronomy, that , difference processing can be carried out in the pseudorange and carrier phase of the reception signal of each epoch to base station and observation station, obtained double Poor pseudorange and double difference carrier phase.

In a kind of implementation of the invention, according to formula:

S102 carries out position resolving according to the position coordinates of base station, the double difference pseudorange in each day and double difference carrier phase, obtains Difference coordinate sequence to observation station in each day.

In the embodiment of the present invention, difference processing can be carried out in RTKLIB (positioning clearing software), eliminate or weaken institute Obtain ionosphere delay, tropospheric delay, satellite orbital error, receiver and the satellite clock correction equal error in difference coordinate sequence.

S103, the difference coordinate sequence by two-dimensional movement weighted mean method to observation station in each day are handled, are obtained Observation station is in the smooth coordinate sequence in each day, and two-dimensional movement weighted mean method is based on time domain space and Value space carries out mobile add Weight average.

Specifically, when observation station carries out static observation to base station (ambient enviroment does not change), Multipath Errors Variation between adjacent epoch tends to be steady, and due to the presence of observation noise, so that the error of coordinate between adjacent epoch fluctuates Aggravation.In order to weaken the observation noise in coordinate sequence obtained by Differential positioning, obtains accurate Multipath Errors and extract model, it can To be filtered pretreatment by two-dimensional movement weighted mean method, weaken observation noise.Two-dimensional movement weighted mean method refers to The moving weighted average method of two dimensions of time domain space and Value space, will hereafter describe to this method in detail.

S104, for it is any refer to day, calculating observation station this with reference to day and other adjacent to day smooth coordinate sequence it is mutual Related coefficient chooses the smooth coordinate sequence that cross-correlation coefficient maximum neighbouring day in day is referred to reference to day and with this.

In the embodiment of the present invention, multipath has strong correlation closing in the day same period, by two-dimensional movement plus The smooth coordinate sequence of observation point after weight average method smoothing processing still contains observation noise, but Multipath Errors are predominantly Position, so also there is very strong correlation.Cross-correlation coefficient ρ_abSize reflect the degree of correlation between two stochastic variables, phase Relationship number | ρ_ab| closer to 1, then A and B is closer to linear correlation, i.e., degree of correlation is bigger.The calculation formula of related coefficient can table It is shown as:

Wherein, A and B is stochastic variable, ρ_abFor the cross-correlation coefficient of variables A and B, For the covariance of variables A and B,For the side of variables A Difference,For the variance of variable B,For become A average value,For the average value of variable B.

In the embodiment of the present invention, A and B respectively refer to smooth coordinate sequence of the smooth coordinate sequence in day with other adjacent to day Column.But due to by remaining observation noise, the visual condition of satellite and observation the distance between station antenna and reflecting surface variation etc. The influence of factor, so that the cross-correlation coefficient of the smooth coordinate sequence closed between day becomes smaller.By between each day of observation point Smooth coordinate sequence carries out correlation analysis, chooses the best smooth coordinate sequence for closing on day of correlation and carries out Multipath Errors Extraction.

S105 is put down to this with reference to day and with this with reference to the maximum neighbouring day of day cross-correlation coefficient by Wavelet Packet Algorithm Sliding coordinate sequence carries out M layers of decomposition, calculates the cross-correlation coefficient of low frequency part in each decomposition layer, and M is pre-set decomposition layer Number.

Wherein, the smooth coordinate sequence of observation station is discrete series, the prior art mainly adopts the tower algorithm of Mallat into Row DMT modulation carries out the extractions of Multipath Errors.But since its scale is changed by binary system, high again and again Its frequency resolution of section is poor, and in low-frequency range, its temporal resolution is poor, smoothly sits so observation station cannot be analyzed well Mark the Multipath Errors of sequence high frequency section.The embodiment of the present invention is by Wavelet Packet Algorithm in the smooth coordinate sequence in observation station Multipath Errors further progress on the basis of two-dimensional movement weighted mean method is extracted.This method can carry out frequency band multi-level It divides, half is further dropped to the high frequency section for not having subdivision in dyadic wavelet transform and is divided, and can be according to analyzed signal Feature is allowed to match with signal spectrum.Based on this characteristic, the multipath of the high frequency section in the smooth coordinate sequence in observation station is missed Difference will be effectively extracted, and be will further impair multipath effect and carried out the influence in deformation observation using high accuracy positioning.

Pre-set Decomposition order M is related to the validity that Multipath Errors extract, and Decomposition order is excessive, it will loses Important information in smooth coordinate sequence, Decomposition order is very few, then cannot be complete by the Multipath Errors in smooth coordinate sequence It extracts.Meet relationship between Decomposition order M and the length I of smooth coordinate sequence: I=2^M+ l, l are remainder item.

Referring to fig. 2, Fig. 2 is the WAVELET PACKET DECOMPOSITION tree structure diagram in the embodiment of the present invention after 3 layers of decomposition, wherein A indicates low Frequency part, D indicate that high frequency section, the serial number at end indicate the number of plies of wavelet decomposition.Exploded relationship is S=AAA3+DAA3+ADA3 +DDA3+AAD3+DAD3+ADD3+DDD3。

S106 judges whether destination layer J is less than M, destination layer J be in each decomposition layer in the cross-correlation coefficient of low frequency part most The number of plies where big value.

The present invention can determine optimal Decomposition according to the selected size for closing on the cross-correlation coefficient that day corresponds to low frequency part Layer, i.e. destination layer, when the cross-correlation coefficient for determining corresponding decomposition layer low frequency part reaches maximum, which is that extraction is more The destination layer of tracking error.Cross-correlation coefficient shows that more greatly the accounting of Multipath Errors is bigger in the decomposition layer, but low frequency portion (Decomposition order is at the beginning as the increase of Decomposition order can show a convex parabolical variation tendency for the cross-correlation coefficient divided Increase can more filter out noise, but continue to increase, and be that the part of Multipath Errors is also filtered out originally), that is to say, that it looks for Maximum cross-correlation coefficient is had found to highest point.Assuming that pre-set Decomposition order is M, if having found maximum at M-1 layers Cross-correlation coefficient then can determine that the point is the most a little louder, but if maximum correlation coefficient appears in M layers, then not can determine that Whether maximum point, layer be compared as soon as only decomposing again, so circulation can determine that maximum cross-correlation coefficient and needs repeatedly The number of plies of decomposition.Therefore, judge whether destination layer J is less than M, if so, executing S107, otherwise, execute S108.

S107, the low frequency part extracted in destination layer J carry out wavelet reconstruction, obtain Multipath Errors.

In this step, after determining destination layer J, by the small of the high frequency section comprising influences such as random noises in destination layer J Wave packet coefficient zero setting carries out wavelet reconstruction to the low frequency part in destination layer J to get Multipath Errors are arrived.

The value of M is added 1, returns to S105 by S108, until destination layer J is less than M.

Error extracting method provided in an embodiment of the present invention obtains base station and observation station in each day preset time period respectively Interior reception signal, for the daily reception signal of base station and observation station, within a preset period of time to base station and observation station The pseudorange and carrier phase of the reception signal of each epoch carry out difference processing, and according to the double difference pseudorange in obtained each day, double difference The position coordinates of carrier phase and base station carry out position resolving, obtain observation station in the difference coordinate sequence in each day；Pass through Two-dimensional movement weighted mean method handles difference coordinate sequence, obtains observation station in the smooth coordinate sequence in each day, two dimension Moving weighted average method is based on time domain space and Value space carries out moving weighted average；By Wavelet Packet Algorithm to cross correlation Number is maximum to carry out M layer with reference to the smooth coordinate sequences in day and neighbouring day and decomposes, and calculates the mutual of low frequency part in each decomposition layer Related coefficient, judges whether destination layer J is less than M, and destination layer J is maximum value in the cross-correlation coefficient of low frequency part in each decomposition layer The number of plies at place；If so, the low frequency part extracted in destination layer J carries out wavelet reconstruction, Multipath Errors are obtained；If not, The value of M is added 1, it is maximum to cross-correlation coefficient with reference to day and adjacent to the smooth coordinate sequence in day to pass back through Wavelet Packet Algorithm The step of carrying out M layers of decomposition, until destination layer J is less than M；The low frequency part extracted in destination layer J carries out wavelet reconstruction, obtains more Tracking error.In the embodiment of the present invention, low frequency part observation noise can be weakened to signal by two-dimensional movement weighted mean method Influence, and solve signal smoothing processing in end effect, achieve the purpose that protect side denoising.Wavelet Packet Algorithm can be by signal Frequency band is divided at many levels, is further decomposed to high frequency section, the observation noise in high frequency section is eliminated, to improve multichannel The accuracy that diameter error is extracted.

In a kind of implementation of the invention, in Fig. 1 embodiment S103, by two-dimensional movement weighted mean method to observation station Difference coordinate sequence in each day is handled, and obtains observation station in the smooth coordinate sequence in each day, comprising the following steps:

The first step determines that the rolling average cycle T of time domain space is no more than I's if the number of sample epoch coordinate is I Odd number.

Wherein, the rolling average cycle T of time domain space and the influence period of multipath effect are consistent, therefore are taken as being not more than The odd number of I.

Second step, if i-th of sample epoch coordinate is x_i, Indicate that Value space, i are the integer of 1~I, root According to formula:Calculate the average value of each sample epoch coordinate

By standard deviation sigma_IAs the threshold θ for choosing Value space epoch coordinate_I, then θ_I=σ_I,

Third step, if being smoothed to r-th of sample epoch coordinate, r is the integer of 1~I；

By standard deviation sigma_rNo more than threshold θ_ISample epoch coordinate as fraction epoch coordinate x_k, fraction epoch coordinate x_k For determining in r-th of smoothed out value of sample epoch coordinate of Value space, fraction epoch coordinate x_kNumber be r-th of sample Rolling average period N of the epoch coordinate in Value space_r, N_r≤I。

In the embodiment of the present invention, rolling average period of Value space according to sample epoch coordinate discrete distribution situation into Row value.Rolling average period N of r-th of sample epoch coordinate in Value space_r, that is to say, that different sample epoch sits Mark, the rolling average period of corresponding different Value space.

4th step, for r-th of sample epoch coordinate, if α_jIndicate the weighting of j-th of time domain space delay epoch coordinate Coefficient, j areInteger, [] indicate round numbers,

According to formula:Calculate each sample epoch coordinate in the rolling average cycle T of time domain space Average retardation epoch coordinate relative to r-th of sample epoch coordinate isΔτ_jIt is j-th of sample epoch coordinate relative to r A sample epoch coordinate x_rDelay epoch,

According to formula:It calculates time domain space and respectively postpones the dissimilarity coefficient of epoch coordinate and be According to formula:It calculates j-th of time domain space and prolongs The weighting coefficient α of slow epoch coordinate_j。

In time domain space, the epoch for calculating sample epoch coordinate and r-th of sample epoch coordinate postpones, and epoch delay is got over Small, the weight of the value of corresponding epoch is bigger.

5th step, according to formula:Calculate fraction epoch coordinate x in Value space_kAverage valueΔ μ_kDifference for k-th of fraction epoch coordinate relative to r-th of sample epoch coordinate；

According to formula:Calculate the weighting coefficient β of k-th of fraction epoch coordinate of Value space_k。

In Value space, the difference of all fraction epoch coordinates and r-th of sample epoch coordinate value is calculated, according to difference Size carries out weight distribution to corresponding each fraction epoch coordinate, and difference is smaller, illustrates that similarity is higher, and corresponding fraction epoch sits Target weight is also bigger.

6th step, according to formula:

7th step, according to formula:Calculate x_rAverage value under two-dimensional space

According to formula:Calculate the dissimilarity coefficient of Value space

It can solve the end effect in smoothing processing by above-mentioned two-dimensional movement weighted mean method, reach and protect side denoising Purpose.

In a kind of implementation of the invention, if this refers to the maximum neighbouring day of day cross-correlation coefficient with reference to day and with this Smooth coordinate sequence be x_i, according to formula:

To x_iWAVELET PACKET DECOMPOSITION is carried out,

Wavelet packet when indicating undecomposed,Indicate the 2nd ω -1 wavelet packet coefficients on M layer,Indicate the ω wavelet packet coefficient on M-1 layer,Indicate the 2nd ω wavelet packet system on M layer Number, G are scaling function resolution filter, and H is wavelet function resolution filter.

In the embodiment of the present invention, when carrying out WAVELET PACKET DECOMPOSITION, since there are different isolations, i.e. subspace in subspace There is different orthogonal basis.In theory analysis, the cost function of a sequence can be usually defined, from all small echos of wavelet library Finding in base makes the smallest base of cost function, and for a given vector, cost minimum is exactly most effective expression, the small echo Base is " optimal base ".In practical applications, wavelet basis is generally selected according to the purpose of signal analysis and experience.In GPS number It theoretically proposes to utilize Symlets small echo according to processing aspect, but mostly uses Daubechies small echo greatly in practical applications.

The low frequency part extracted in destination layer J carries out wavelet reconstruction, obtains Multipath Errors, comprising:

According to formula:To x_iCarry out wavelet packet weight Structure；

Specifically, by the wavelet packet coefficient zero setting of the high frequency section comprising influences such as random noises to get arrive coordinate residual error Multipath Errors in sequence.

Embodiment one

If analog signal model are as follows:

S by a cycle be 1200s cosine signal and two periods be respectively 900s and 300s sinusoidal signal and One Gaussian sequence e composition.The noise level for taking e is N (0,1.52), and data sampling rate is set as 1s, and number of samples is 3000。

The analog signal is passed through into wavelet threshold denoising algorithm, two-dimensional movement Weighted Average Algorithm respectively and based on two dimension The Wavelet Packet Algorithm of moving weighted average processing is filtered denoising, and wavelet function selects db8, and Decomposition order is 5 layers, denoising Signal is as shown in figure 3, residual signals are as shown in Figure 4.

The Wavelet Packet Algorithm based on two-dimensional movement weighted average processing of the embodiment of the present invention it can be seen from Fig. 3 and Fig. 4 Signal after denoising can preferably restore original signal, and slightly through the processed signals revivification of two-dimensional movement Weighted Average Algorithm It is better than wavelet threshold denoising algorithm.

The signal restored after the denoising of three of the above algorithm and the residual signals of generation are subjected to quantitative analysis, analysis indexes packet Include the S after denoising in signal_RMSAnd the correlation coefficient ρ of the signal and original signal after denoising.For three kinds of algorithms of scientific evaluation Performance, simulation analysis has been carried out under 4 kinds of different white Gaussian noise simulated environments respectively, as a result such as 3 institute of table 1, table 2 and table Show.Table 1 is the performance indicator of wavelet threshold denoising under different white Gaussian noise simulated environments, and table 2 is different white Gaussian noise moulds The performance indicator that two-dimensional movement weighted mean method denoises under near-ring border, table 3 are that two are based under different white Gaussian noise simulated environments Tie up the performance indicator of the wavelet packet denoising of moving weighted average processing.

Table 1

Noise level	0.5	1.5	2.5	3.5
					S_RMS	0.0982	0.1648	0.1732	0.2231
ρ	0.9886	0.9539	0.9524	0.9493

Table 2

Noise level	0.5	1.5	2.5	3.5
					S_RMS	0.0268	0.0831	0.1638	0.2285
ρ	0.9892	0.9802	0.9742	0.9698

Table 3

Noise level	0.5	1.5	2.5	3.5
					S_RMS	0.0087	0.0401	0.0722	0.0836
ρ	0.9998	0.9882	0.9832	0.9801

By table 1, table 2 and table 3 it is found that the related coefficient of above-mentioned three kinds of algorithms index under different noise levels is most of all It has been more than 0.95, has illustrated that they can preferably restore the information of original signal.It is small but with the continuous increase of noise level The S of wave Threshold Filter Algorithms and two-dimensional movement weighted average Processing Algorithm_RMSIt obviously increases, but wavelet threshold denoising is in strong noise It is increased more slow under level condition, and two-dimensional movement weighted average Processing Algorithm is relatively fast, but its S_RMSIn strong noise water Flat lower is so better than wavelet threshold denoising algorithm.

The S of Wavelet Packet Algorithm based on two-dimensional movement weighted average processing_RMSIt is all relatively optimal under different noise levels, And its related coefficient is both greater than 0.98.

In conclusion the Wavelet Packet Algorithm based on two-dimensional movement weighted average processing can be preferably to adding the signal after making an uproar It is filtered denoising, shows that the present invention can accurately extract Multipath Errors.

Corresponding to above method embodiment, the embodiment of the invention provides a kind of error extraction elements, are referring to Fig. 5, Fig. 5 The structure chart of the error extraction element of the embodiment of the present invention, comprising:

Double difference computing module 501, for obtaining the reception letter of base station and observation station in each day preset time period respectively Number, for the daily reception signal of base station and observation station, to base station and observation station, each epoch is connect within a preset period of time The pseudorange and carrier phase of the collection of letters number carry out difference processing, obtain double difference pseudorange and double difference carrier phase；

Difference coordinate sequence computing module 502, for according to the position coordinates of base station, the double difference pseudorange in each day and double difference Carrier phase carries out position resolving, obtains observation station in the difference coordinate sequence in each day；

Smoothing module 503, for passing through difference coordinate sequence of the two-dimensional movement weighted mean method to observation station in each day Column are handled, and obtain observation station in the smooth coordinate sequence in each day, two-dimensional movement weighted mean method is based on time domain space and value Domain space carries out moving weighted average；

Smooth coordinate sequence chooses module 504, for referring to day for any, calculating observation station this with reference to day and other are adjacent It is maximum neighbouring with reference to day cross-correlation coefficient with reference to day and with this to choose this for the cross-correlation coefficient of the smooth coordinate sequence in nearly day It smooth coordinate sequence；

Decomposing module 505, for referring to day cross-correlation coefficient maximum by Wavelet Packet Algorithm to this with reference to day and with this Neighbouring day smooth coordinate sequence carry out M layer decompose, calculate the cross-correlation coefficient of low frequency part in each decomposition layer, M is preparatory The Decomposition order of setting；

Judgment module 506, for judging whether destination layer J is less than M, destination layer J is the mutual of low frequency part in each decomposition layer The number of plies in related coefficient where maximum value；

Reconstructed module 507, for the judging result in judgment module be when, extract destination layer J in low frequency part into Row wavelet reconstruction, obtains Multipath Errors；

Loop module 508, for the judging result in judgment module be it is no when, the value of M is added 1, passes back through wavelet packet Algorithm carries out M layers of decomposition to this with reference to day and with the smooth coordinate sequence with reference to the maximum neighbouring day of day cross-correlation coefficient Step, until destination layer J is less than M.

Error extraction element provided in an embodiment of the present invention obtains base station and observation station in each day preset time period respectively Interior reception signal, for the daily reception signal of base station and observation station, within a preset period of time to base station and observation station The pseudorange and carrier phase of the reception signal of each epoch carry out difference processing, and according to the double difference pseudorange in obtained each day, double difference The position coordinates of carrier phase and base station carry out position resolving, obtain observation station in the difference coordinate sequence in each day；Pass through Two-dimensional movement weighted mean method handles difference coordinate sequence, obtains observation station in the smooth coordinate sequence in each day, two dimension Moving weighted average method is based on time domain space and Value space carries out moving weighted average；By Wavelet Packet Algorithm to cross correlation Number is maximum to carry out M layer with reference to the smooth coordinate sequences in day and neighbouring day and decomposes, and calculates the mutual of low frequency part in each decomposition layer Related coefficient, judges whether destination layer J is less than M, and destination layer J is maximum value in the cross-correlation coefficient of low frequency part in each decomposition layer The number of plies at place；If so, the low frequency part extracted in destination layer J carries out wavelet reconstruction, Multipath Errors are obtained；If not, The value of M is added 1, it is maximum to cross-correlation coefficient with reference to day and adjacent to the smooth coordinate sequence in day to pass back through Wavelet Packet Algorithm The step of carrying out M layers of decomposition, until destination layer J is less than M；The low frequency part extracted in destination layer J carries out wavelet reconstruction, obtains more Tracking error.In the embodiment of the present invention, low frequency part observation noise can be weakened to signal by two-dimensional movement weighted mean method Influence, and solve signal smoothing processing in end effect, achieve the purpose that protect side denoising.Wavelet Packet Algorithm can be by signal Frequency band is divided at many levels, is further decomposed to high frequency section, the observation noise in high frequency section is eliminated, to improve multichannel The accuracy that diameter error is extracted.

Optionally, smoothing module is specifically used for, if the number of sample epoch coordinate is I, determines the shifting of time domain space Dynamic T average period is the odd number no more than I；

By standard deviation sigma_rNo more than threshold θ_ISample epoch coordinate as fraction epoch coordinate x_k, fraction epoch coordinate x_k For determining in r-th of smoothed out value of sample epoch coordinate of Value space, fraction epoch coordinate x_kNumber be r-th of sample Rolling average period N of the epoch coordinate in Value space_r, N_r≤I；For r-th of sample epoch coordinate, if α_jIndicate that time domain is empty Between j-th of delay epoch coordinate weighting coefficient, j isInteger, [] indicate round numbers,

According to formula:

According to formula:Calculate x_rAverage value under two-dimensional space

According to formula:Calculate the dissimilarity coefficient of Value space

Optionally, double difference computing module is specifically used for according to formula:

Wherein,For double difference operator, S indicates global position system GPS system, and satellite n and satellite m are in the GPS system Satellite,Indicate observation station r and base station b in same epoch observation satellite n and satellite m satellite hub to connecing Geometric distance between receipts machine phase center, λ are carrier wavelength,It is seen for observation station r and base station b in same epoch The unknown integer ambiguity of carrier phase of satellite n and satellite m is surveyed,Indicate base station b and observation station r same The pseudorange multipath error of epoch observation satellite n and satellite m,Indicate that observation station r and base station b is gone through same The carrier phase Multipath Errors of first observation satellite n and satellite m,Indicate that observation station r and base station b is gone through same The pseudo-code of first observation satellite n and satellite m,Indicate observation station r and base station b in same epoch observation satellite n and The carrier phase measurement noise of satellite m.

Optionally, decomposing module, if maximum neighbouring with reference to day cross-correlation coefficient with reference to day and with this specifically for this It smooth coordinate sequence is x_i, according to formula:

To x_iWAVELET PACKET DECOMPOSITION is carried out,

Reconstructed module is specifically used for according to formula:

To x_iIt carries out Wavelet package reconstruction；Indicate the 2nd ω -1 wavelet packet coefficients on J+1 layer,It indicates on J+1 layer 2nd ω wavelet packet coefficient, g are scaling function reconfigurable filter, and h is wavelet function reconfigurable filter.

The embodiment of the invention also provides a kind of electronic equipment, are the electronic equipment of the embodiment of the present invention referring to Fig. 6, Fig. 6 Structure chart, comprising: processor 601, communication interface 602, memory 603 and communication bus 604, wherein processor 601, logical Letter interface 602, memory 603 complete mutual communication by communication bus 604；

Memory 603, for storing computer program；

Processor 601 when for executing the program stored on memory 603, realizes any of the above-described error extracting method The step of.

It should be noted that the communication bus 604 that above-mentioned electronic equipment is mentioned can be PCI (Peripheral Component Interconnect, Peripheral Component Interconnect standard) bus or EISA (Extended Industry Standard Architecture, expanding the industrial standard structure) bus etc..The communication bus 604 can be divided into address bus, data/address bus, Control bus etc..Only to be indicated with a thick line in Fig. 6, it is not intended that an only bus or a seed type convenient for indicating Bus.

Communication interface 602 is for the communication between above-mentioned electronic equipment and other equipment.

Memory 603 may include RAM (Random Access Memory, random access memory), also may include Nonvolatile memory (non-volatile memory), for example, at least a magnetic disk storage.Optionally, memory may be used also To be storage device that at least one is located remotely from aforementioned processor.

Above-mentioned processor 601 can be general processor, comprising: CPU (Central Processing Unit, center Processor), NP (Network Processor, network processing unit) etc.；It can also be DSP (Digital Signal Processing, digital signal processor), ASIC (Application Specific Integrated Circuit, it is dedicated Integrated circuit), FPGA (Field-Programmable Gate Array, field programmable gate array) or other are programmable Logical device, discrete gate or transistor logic, discrete hardware components.

The embodiment of the invention also provides a kind of computer readable storage medium, it is stored in computer readable storage medium Computer program, when computer program is executed by processor, the step of realizing any of the above-described error extracting method.

It should be noted that, in this document, relational terms such as first and second and the like are used merely to a reality Body or operation are distinguished with another entity or operation, are deposited without necessarily requiring or implying between these entities or operation In any actual relationship or order or sequence.Moreover, the terms "include", "comprise" or its any other variant are intended to Non-exclusive inclusion, so that the process, method, article or equipment including a series of elements is not only wanted including those Element, but also including other elements that are not explicitly listed, or further include for this process, method, article or equipment Intrinsic element.In the absence of more restrictions, the element limited by sentence "including a ...", it is not excluded that There is also other identical elements in process, method, article or equipment including the element.

Each embodiment in this specification is all made of relevant mode and describes, same and similar portion between each embodiment Dividing may refer to each other, and each embodiment focuses on the differences from other embodiments.Especially for device, For electronic equipment and readable storage medium storing program for executing embodiment, since it is substantially similar to the method embodiment, so the comparison of description is simple Single, the relevent part can refer to the partial explaination of embodiments of method.

The foregoing is merely illustrative of the preferred embodiments of the present invention, is not intended to limit the scope of the present invention.It is all Any modification, equivalent replacement, improvement and so within the spirit and principles in the present invention, are all contained in protection scope of the present invention It is interior.

Claims

1. a kind of error extracting method, which is characterized in that the described method includes:

The reception signal of base station and observation station in each day preset time period is obtained respectively, for the base station and the sight The daily reception signal of survey station, to the reception signal of the base station and the observation station each epoch in the preset time period Pseudorange and carrier phase carry out difference processing, obtain double difference pseudorange and double difference carrier phase；

Position resolving is carried out according to the position coordinates of the base station, the double difference pseudorange in each day and double difference carrier phase, obtains institute Observation station is stated in the difference coordinate sequence in each day；

Difference coordinate sequence by two-dimensional movement weighted mean method to the observation station in each day is handled, and obtains the sight Survey station is in the smooth coordinate sequence in each day, and the two-dimensional movement weighted mean method is based on time domain space and Value space is moved Weighted average；

For it is any refer to day, calculate the observation station this with reference to day and other adjacent to the smooth coordinate sequence in day cross correlation Number chooses the smooth coordinate sequence that cross-correlation coefficient maximum neighbouring day in day is referred to reference to day and with this；

The smooth coordinate sequence in the maximum neighbouring day of day cross-correlation coefficient is referred to this with reference to day and with this by Wavelet Packet Algorithm Column carry out M layers of decomposition, calculate the cross-correlation coefficient of low frequency part in each decomposition layer, and M is pre-set Decomposition order, judge mesh Whether mark layer J is less than M, and the destination layer J is where maximum value in the cross-correlation coefficient of low frequency part in each decomposition layer The number of plies；

If not, the value of M is added 1, returns to the Wavelet Packet Algorithm that passes through and day is referred to this and refers to day cross correlation with this The step of smooth coordinate sequence in the maximum neighbouring day of number carries out M layers of decomposition, until the destination layer J is less than M.

2. error extracting method according to claim 1, which is characterized in that described to pass through two-dimensional movement weighted mean method pair Difference coordinate sequence of the observation station in each day is handled, and obtains the observation station in the smooth coordinate sequence in each day, packet It includes:

If the number of the sample epoch coordinate is I, determine that the rolling average cycle T of the time domain space is the surprise no more than I Number；

According to formula:Calculate each sample epoch coordinate and current sample epoch coordinate x in Value space_r Standard deviation sigma_r,

By standard deviation sigma_rNo more than threshold θ_ISample epoch coordinate as fraction epoch coordinate x_k, the fraction epoch coordinate x_kWith In determining in r-th of smoothed out value of sample epoch coordinate of Value space, fraction epoch coordinate x_kNumber be r-th of sample go through Rolling average period N of first coordinate in Value space_r, N_r≤I；For r-th of sample epoch coordinate, if α_jIndicate time domain space The weighting coefficient of j-th of delay epoch coordinate, j areInteger, [] indicate round numbers,

According to formula:Calculate each sample epoch coordinate phase in the rolling average cycle T of the time domain space Average retardation epoch coordinate for r-th of sample epoch coordinate isΔτ_jIt is j-th of sample epoch coordinate relative to r-th Sample epoch coordinate x_rDelay epoch,

According to formula:It calculates time domain space and respectively postpones the dissimilarity coefficient of epoch coordinate and bej∈T,

According to formula:Calculate the weighting coefficient α of j-th of time domain space delay epoch coordinate_j；

According to formula:Calculate fraction epoch coordinate x in Value space_kAverage valueΔμ_kIt is k-th point Weigh difference of the epoch coordinate relative to r-th of sample epoch coordinate；

According to formula:

By r-th of sample epoch coordinate x_rRespectively in time domain space and Value space with x_rCentered on carry out moving weighted average at Reason, obtains the time domain space moving weighted average value of r-th of sample epoch coordinateAnd Value space rolling average weighted value

According to formula:Calculate x_rAverage value under two-dimensional space

According to formula:Calculate the dissimilarity coefficient of Value space

According to formula:Calculate the weighting coefficient γ of time domain space_T, according to formula:Meter Calculate the weighting coefficient of Value space

According to formula:It willWithIt is weighted processing, is obtained after smoothing processing Coordinate

3. error extracting method according to claim 1, which is characterized in that described to the base station and the observation station In the preset time period each epoch reception signal pseudorange and carrier phase carry out difference processing, obtain double difference pseudorange and Double difference carrier phase, comprising:

According to formula:

It calculates and sees Double difference pseudorange of the survey station r and base station b in same epoch observation satellite n and satellite mAnd observation station r and benchmark Stand b same epoch observation satellite n and satellite m double difference carrier phase

Wherein,For double difference operator, S indicates global position system GPS system, and satellite n and satellite m are defending in the GPS system Star,Indicate observation station r and base station b in same epoch observation satellite n and satellite m satellite hub to receiver Geometric distance between phase center, λ are carrier wavelength；For observation station r and base station b, in same epoch, observation is defended The unknown integer ambiguity of the carrier phase of star n and satellite m,Indicate base station b and observation station r in same epoch The pseudorange multipath error of observation satellite n and satellite m,Indicate that observation station r and base station b is seen in same epoch The carrier phase Multipath Errors of satellite n and satellite m are surveyed,Indicate that observation station r and base station b is seen in same epoch The pseudo-code of satellite n and satellite m are surveyed,Indicate observation station r and base station b in same epoch observation satellite n and satellite m Carrier phase measurement noise.

4. error extracting method according to claim 1, which is characterized in that described to refer to day to this by Wavelet Packet Algorithm And M layers of decomposition are carried out with the smooth coordinate sequence with reference to the maximum neighbouring day of day cross-correlation coefficient, comprising:

To x_iWAVELET PACKET DECOMPOSITION is carried out,

According to formula:To x_i Carry out wavelet package reconstruction；

Indicate the 2nd ω -1 wavelet packet coefficients on J+1 layer,Indicate the 2nd ω on J+1 layer Wavelet packet coefficient, g are scaling function reconfigurable filter, and h is wavelet function reconfigurable filter.

5. a kind of error extraction element, which is characterized in that described device includes:

Double difference computing module, for obtaining the reception signal of base station and observation station in each day preset time period respectively, for The daily reception signal of the base station and the observation station, to the base station and the observation station in the preset time period The pseudorange and carrier phase of the reception signal of interior each epoch carry out difference processing, obtain double difference pseudorange and double difference carrier phase；

Difference coordinate sequence computing module, for being carried according to the position coordinates of the base station, the double difference pseudorange in each day and double difference Wave phase carries out position resolving, obtains the observation station in the difference coordinate sequence in each day；

Smoothing module, for by two-dimensional movement weighted mean method to the observation station each day difference coordinate sequence into Row processing obtains the observation station in the smooth coordinate sequence in each day, and the two-dimensional movement weighted mean method is based on time domain space Moving weighted average is carried out with Value space；

Smooth coordinate sequence chooses module, and for referring to day for any, calculating the observation station, this is neighbouring with other with reference to day The cross-correlation coefficient of it smooth coordinate sequence chooses this with reference to day and with this with reference to the maximum neighbouring day of day cross-correlation coefficient Smooth coordinate sequence；

Decomposing module, for referring to the maximum neighbouring day of day cross-correlation coefficient by Wavelet Packet Algorithm to this with reference to day and with this Smooth coordinate sequence carry out M layer decomposition, calculate the cross-correlation coefficient of low frequency part in each decomposition layer, M is pre-set point Solve the number of plies；

Judgment module, for judging whether destination layer J is less than M, the destination layer J is the mutual of low frequency part in each decomposition layer The number of plies in related coefficient where maximum value；

Reconstructed module, for when the judgment result of the judgment module is yes, extract low frequency part in the destination layer J into Row wavelet reconstruction, obtains Multipath Errors；

Loop module, for the judging result in the judgment module be it is no when, the value of M is added 1, passes through wavelet packet described in return Algorithm carries out M layers of decomposition to this with reference to day and with the smooth coordinate sequence with reference to the maximum neighbouring day of day cross-correlation coefficient Step, until the destination layer J is less than M.

6. error extraction element according to claim 5, which is characterized in that the smoothing module is specifically used for, if The number of the sample epoch coordinate is I, determines that the rolling average cycle T of the time domain space is the odd number no more than I；

According to formula:

According to formula:Calculate x_rAverage value under two-dimensional space

According to formula:Calculate the dissimilarity coefficient of Value space

7. error extraction element according to claim 5, which is characterized in that the double difference computing module is specifically used for root According to formula:

Wherein,For double difference operator, S indicates global position system GPS system, and satellite n and satellite m are defending in the GPS system Star,Indicate observation station r and base station b in same epoch observation satellite n and satellite m satellite hub to receiver Geometric distance between phase center, λ are carrier wavelength,For observation station r and base station b, in same epoch, observation is defended The unknown integer ambiguity of the carrier phase of star n and satellite m,Indicate base station b and observation station r in same epoch The pseudorange multipath error of observation satellite n and satellite m,Indicate that observation station r and base station b is seen in same epoch The carrier phase Multipath Errors of satellite n and satellite m are surveyed,Indicate that observation station r and base station b is seen in same epoch The pseudo-code of satellite n and satellite m are surveyed,Indicate observation station r and base station b in same epoch observation satellite n and satellite m Carrier phase measurement noise.

8. error extraction element according to claim 5, which is characterized in that the decomposing module, if being specifically used for the ginseng It examines day and refers to the smooth coordinate sequence in the maximum neighbouring day of day cross-correlation coefficient with this as x_i, according to formula:

To x_iWAVELET PACKET DECOMPOSITION is carried out,

The reconstructed module is specifically used for according to formula:

To x_iCarry out small echo Packet reconstruct；

9. a kind of electronic equipment characterized by comprising processor, communication interface, memory and communication bus, wherein described Processor, the communication interface, the memory complete mutual communication by the communication bus；

The memory, for storing computer program；

The processor when for executing the program stored on the memory, realizes any mistake of claim 1-4 The step of poor extracting method.

10. a kind of computer readable storage medium, which is characterized in that be stored with computer in the computer readable storage medium Program when the computer program is executed by processor, realizes the step of any error extracting method of claim 1-4 Suddenly.