阅读说明：本技术 基于wmad的gnss卫星钟差数据预处理方法 (GNSS satellite clock deviation data preprocessing method based on WMAD ) 是由 王旭 柴洪洲 王昶 李小舟 于 2019-08-15 设计创作，主要内容包括：本发明的目的是提供一种基于WMAD的GNSS卫星钟差数据预处理方法。本发明有效解决了经典粗差探测方法不能有效剔除钟差中量级较小的误差问题,本发明的方法具体是首先对GNSS卫星钟差进行一次差计算,得到钟差一次差分数据；然后对钟差一次差分数据进行小波分解,得到分解后小波系数；再利用MAD方法剔除较大的小波系数；计算各层阈值,进行阈值化处理；最后小波重构得到预处理后的钟差一次差分数据。本发明为卫星钟差数据预处理提供一种有效方法,提高钟差数据的质量。(The object of the present invention is to provide a kind of GNSS satellite clock deviation data preprocessing method based on WMAD.The present invention, which efficiently solves classical Detection of Gross Errors method, cannot effectively reject the lesser error problem of clock deviation middleweight, and method of the invention is specifically to carry out first difference calculating to GNSS satellite clock deviation first, obtain clock deviation first difference data；Then wavelet decomposition, wavelet coefficient after being decomposed are carried out to clock deviation first difference data；MAD method is recycled to reject biggish wavelet coefficient；Each layer threshold value is calculated, thresholding processing is carried out；Last wavelet reconstruction obtains pretreated clock deviation first difference data.The present invention provides a kind of effective ways for satellite clock correction data prediction, improves the quality of clock deviation data.)

1. the GNSS satellite clock deviation data preprocessing method based on WMAD, which comprises the following steps:

A: first difference calculating is carried out to GNSS satellite clock deviation, obtains clock deviation first difference data；

B: wavelet decomposition, wavelet coefficient after being decomposed are carried out to clock deviation first difference data；

C: biggish wavelet coefficient is rejected using MAD method；

D: calculating each layer threshold value, carries out thresholding processing；

E: wavelet reconstruction obtains pretreated clock deviation first difference data.

2. the GNSS satellite clock deviation data preprocessing method based on WMAD as described in claim 1, which is characterized in that described Step A the following steps are included:

A: the clock deviation value that L={ l (i), i=1,2 ..., n } is one group of difference epoch-making moment is set first；

B: satellite clock correction first difference data are calculated, i.e.,

Δ L (j)=l (i+1)-l (i), j=1,2 ... n-1

Wherein, Δ L (j) is clock deviation first difference data.

3. the GNSS satellite clock deviation data preprocessing method based on WMAD as described in claim 1, which is characterized in that described Step B the following steps are included:

C: wavelet function is determined；

It is utilized respectively different wavelet functions and wavelet decomposition is carried out to clock deviation first difference data, the wavelet coefficient after being decomposed, I.e.

Wherein, WT_f(a, b) is wavelet conversion coefficient, and f (g) is satellite clock correction first difference data, ψ_a,bIt (g) is wavelet mother function, A is contraction-expansion factor, and b shift factor, wherein g=0,1,2 ... T-1, T are satellite clock correction first difference data amount check；

The energy value of wavelet coefficient after calculating wavelet transformation, i.e.,

Wherein, E_sIt (W) is the energy value for extracting signal, N is the number of wavelet coefficient, and W (s, t) is wavelet coefficient, and s indicates scale Parameter, t indicate signal amplitude；The wavelet function for obtaining maximum energy value is the optimal wavelet function of selection；

D: the wavelet decomposition number of plies is determined；

According to the corresponding relationship of scale and frequency, it may be assumed that

Wherein f_s ^psIt is the corresponding pseudo frequency ps of scale s (pseudo frequency), f_cIt is the corresponding centre frequency of small echo, Δ t It is using the period；Pseudo frequency should cover useful signal f_sigWhole frequency ranges, it may be assumed that

Therefore minimum frequency (the minf of useful signal_sig) it should be the pseudo frequency of L not less than scale s, that is:

The Decomposition order of wavelet transform is that two are carried out to scale s into discrete sampling, therefore the corresponding Decomposition order j of L are as follows:

2^j-1≤L≤2^j

E: wavelet decomposition is carried out to clock deviation data and obtains wavelet coefficient.

4. the GNSS satellite clock deviation data preprocessing method based on WMAD as described in claim 1, which is characterized in that described Step C the following steps are included:

F: and then to every layer of wavelet coefficient WT_fMiddle number (MED) m of (a, b) and wavelet coefficient sequence adds the n of median (MAD) Times the sum of compare, even wavelet coefficient meet

|WT_f(a, b) | > (m+n × MAD)

It is exactly larger exceptional value that it is corresponding, which to be considered as the wavelet coefficient, while rejecting the wavelet coefficient.

5. the GNSS satellite clock deviation data preprocessing method based on WMAD as described in claim 1, which is characterized in that described Step D the following steps are included:

G: every layer of generic threshold value is calculated, it may be assumed that

Wherein: λ is generic threshold value, and σ is the standard deviation of noise signal, and N is the number of wavelet coefficient；

H: soft-threshold processing, soft-threshold function expression formula are carried out to high frequency coefficient are as follows:

Wherein WT_f(a, b) is wavelet coefficient,To estimate wavelet coefficient, λ is threshold value.

6. the GNSS satellite clock deviation data preprocessing method based on WMAD as described in claim 1, which is characterized in that described Step E the following steps are included:

I: wavelet reconstruction restores the clock deviation first difference data after denoising, the calculation formula of reconstruct are as follows:

Wherein, Sf (a, b) is low frequency coefficient, WT_f(a, b) is wavelet coefficient,WithIt corresponds respectively to reconstruct low Bandpass filter and reconstruct high-pass filter.

Technical field

The present invention relates to a kind of GNSS satellite clock deviation data processing methods, belong to geodesic survey and navigation surveying and mapping technology neck Domain.

Background technique

GNSS satellite atomic clock will appear the operations such as frequency modulation, clock switching in long-term operation, while also suffer from various The influence of uncertain factor, clock deviation data inevitably contain some rough errors, if carried out with the clock deviation data containing rough error Data analysis and Modeling and Prediction can largely effect on the quality of data analysis and Modeling and Prediction.Therefore, domestic and foreign scholars are in satellite clock A large amount of research has been carried out in terms of poor elimination of rough difference, the rough error in clock deviation data has been handled, common method has median (Median Absolute Deviation, MAD) method, Robust filter method and Bayesian method etc..Wherein, the side MAD Method has many advantages, such as that principle is simple, computational efficiency is high, is a kind of method generallyd use in current clock deviation data prediction.

The rough error contained in clock deviation data can be effectively rejected by the above method, but in actual clock deviation treatment process In, discovery, which has been rejected, can often have the lesser error of magnitude (exceptional value) in the clock deviation data of rough error, and the Detection of Gross Errors of classics Method these exceptional values are difficult effectively to be rejected, and these exceptional values also can to influence clock deviation data analysis and build The quality of mould forecast, so being carried out to the lesser error of clock deviation data middleweight rationally effective before using clock deviation data Pre-process the step that is very important.

Summary of the invention

The purpose of the present invention is to propose to a kind of the GNSS satellite clock deviation data preprocessing method based on WMAD, present invention analysis Wavelet analysis and median method correspond to meaning, and two methods fusion is referred to GNSS satellite clock deviation data prediction mistake Cheng Zhong improves the quality of GNSS satellite clock deviation data prediction.

To achieve the goals above, the present invention adopts the following technical solutions:

GNSS satellite clock deviation data preprocessing method based on WMAD, comprising the following steps:

A: first difference calculating is carried out to GNSS satellite clock deviation, obtains clock deviation first difference data；

B: wavelet decomposition, wavelet coefficient after being decomposed are carried out to GNSS satellite clock deviation first difference data；

C: biggish wavelet coefficient is rejected using MAD method；

D: calculating each layer threshold value, carries out thresholding processing；

E: wavelet reconstruction obtains pretreated clock deviation first difference data.

Further, the step A the following steps are included:

A: the clock deviation value that L={ l (i), i=1,2 ..., n } is one group of difference epoch-making moment is set first；

B: satellite clock correction first difference data are calculated, i.e.,

Δ L (j)=l (i+1)-l (i), j=1,2 ... n-1

Wherein, Δ L (j) is clock deviation first difference data.

Further, the step B the following steps are included:

C: wavelet function is determined；

It is utilized respectively different wavelet functions and wavelet decomposition is carried out to clock deviation first difference data, the wavelet systems after being decomposed Number, i.e.,

Wherein, WT_f(a, b) is wavelet conversion coefficient, and f (g) is satellite clock correction first difference data, ψ_a,b(g) female for small echo Function, a are contraction-expansion factor, b shift factor (g=0,1,2 ... T-1, T are satellite clock correction first difference data amount check)；

The energy value of wavelet coefficient after calculating wavelet transformation, i.e.,

Wherein, E_sIt (W) is the energy value for extracting signal, N is the number of wavelet coefficient, and W (s, t) is wavelet coefficient, and s is indicated Scale parameter, t indicate signal amplitude；The wavelet function for obtaining maximum energy value is the optimal wavelet function of selection；

D: the wavelet decomposition number of plies is determined；

According to the corresponding relationship of scale and frequency, it may be assumed that

Wherein f_s ^psIt is the corresponding pseudo frequency ps of scale s (pseudo frequency), f_cIt is the corresponding center frequency of small echo Rate, Δ t are using the period；Pseudo frequency should cover useful signal f_sigWhole frequency ranges, it may be assumed that

Therefore minimum frequency (the minf of useful signal_sig) it should be the pseudo frequency of L not less than scale s, that is:

The Decomposition order of wavelet transform is that two are carried out to scale s into discrete sampling, therefore the corresponding Decomposition order j of L Are as follows:

2^j-1≤L≤2^j

E: wavelet decomposition is carried out to clock deviation data and obtains wavelet coefficient.

Further, the step C the following steps are included:

F: and then to every layer of wavelet coefficient WT_fMiddle number (MED) m of (a, b) and wavelet coefficient sequence adds median (MAD) the sum of n times compares, and even wavelet coefficient meets

|WT_f(a, b) | > (m+n × MAD)

It is exactly larger exceptional value that it is corresponding, which to be considered as the wavelet coefficient, while rejecting the wavelet coefficient.

Further, the step D the following steps are included:

G: every layer of generic threshold value is calculated, it may be assumed that

Wherein.λ is generic threshold value, and σ is the standard deviation of noise signal, and N is the number of wavelet coefficient.

H: soft-threshold processing, soft-threshold function expression formula are carried out to high frequency coefficient are as follows:

Wherein WT_f(a, b) is wavelet coefficient,To estimate wavelet coefficient, λ is threshold value.

Further, the step E the following steps are included:

I: wavelet reconstruction restores the clock deviation first difference data after denoising, the calculation formula of reconstruct are as follows:

Wherein, Sf (a, b) is low frequency coefficient, WT_f(a, b) is wavelet coefficient,WithCorrespond respectively to weight Structure low-pass filter and reconstruct high-pass filter.

Compared with prior art, the beneficial effects of the present invention are: the present invention proposes a kind of GNSS clock deviation based on WMAD The method of data prediction, wavelet analysis (Wavelet) can effectively remove number due to having the characteristics that multiresolution analysis The smaller exceptional value as caused by white noise in, and median method can effectively reject number by adjusting the multiple n value of MAD The larger exceptional value in.The present invention analyzes wavelet analysis and median method corresponds to meaning, and (i.e. by two methods fusion WMAD it) refers in GNSS satellite clock deviation process of data preprocessing, improves the quality of GNSS satellite clock deviation data prediction.It is first First, multilevel wavelet decomposition is carried out to clock deviation data using wavelet function, the wavelet coefficient after being decomposed.Secondly, utilizing the side MAD Method rejects biggish wavelet coefficient, and the wavelet coefficient that recycles that treated calculates the threshold value of each layer and carries out thresholding processing, Finally, by treated, wavelet coefficient carries out clock deviation data after wavelet reconstruction is pre-processed.Therefore, the present invention combines two kinds of sides The advantages of method proposes a kind of clock deviation preprocess method of wavelet analysis in conjunction with median method, pre-processes to clock deviation data, To improve the quality of clock deviation data.

Detailed description of the invention

Fig. 1 is the GNSS satellite clock deviation data preprocessing method flow chart of the present invention based on WMAD；

Fig. 2 is PRN02 clock deviation first difference component in specific embodiment；

Fig. 3 is a differential data display figure after MAD pretreatment；

Fig. 4 is the pretreated first difference data display figure of wavelet threshold algorithms；

Fig. 5 is a differential data display figure after WMAD pretreatment；

Fig. 6 is that MAD pretreatment, wavelet threshold algorithms and WMAD of the present invention pre-process three kinds of pretreated forecast of method Error Graph.

Specific embodiment

Below in conjunction with embodiment, technical scheme in the embodiment of the invention is clearly and completely described, it is clear that Described embodiments are only a part of the embodiments of the present invention, instead of all the embodiments.Based on the implementation in the present invention Example, every other embodiment obtained by those of ordinary skill in the art without making creative efforts belong to The scope of protection of the invention.

Referring to Fig.1, the GNSS satellite clock deviation data preprocessing method based on WMAD, comprising the following steps:

A: first difference calculating is carried out to satellite clock correction, obtains clock deviation first difference data；

B: wavelet decomposition, wavelet coefficient after being decomposed are carried out to GNSS satellite clock deviation first difference data；

C: biggish wavelet coefficient is rejected using MAD method；

D: calculating each layer threshold value, carries out thresholding processing；

E: wavelet reconstruction obtains pretreated clock deviation data.

Further, the step A the following steps are included:

A: the clock deviation value that L={ l (i), i=1,2 ..., n } is one group of difference epoch-making moment is set first；

B: satellite clock correction first difference data are calculated, i.e.,

Δ L (j)=l (i+1)-l (i), j=1,2 ... n-1

Wherein, Δ L (j) is clock deviation first difference data.

Further, the step B the following steps are included:

C: wavelet function is determined；

It is utilized respectively different wavelet functions and wavelet decomposition is carried out to clock deviation first difference data, the wavelet systems after being decomposed Number, i.e.,

Wherein, WT_f(a, b) is wavelet conversion coefficient, and f (g) is satellite clock correction first difference data, ψ_a,b(g) female for small echo Function, a are contraction-expansion factor, b shift factor (g=0,1,2 ... T-1, T are satellite clock correction first difference data amount check).

The energy value of wavelet coefficient after calculating wavelet transformation, i.e.,

Wherein, E_sIt (W) is the energy value for extracting signal, N is the quantity of wavelet coefficient, and W (s, t) is wavelet coefficient, and s is indicated Scale parameter, t indicate signal amplitude.The wavelet function for obtaining maximum energy value is the optimal wavelet function of selection

D: the wavelet decomposition number of plies is determined；

According to the corresponding relationship of scale and frequency, it may be assumed that

Wherein f_s ^psIt is the corresponding pseudo frequency of scale s (pseudo frequency), f_cIt is the corresponding centre frequency of small echo, Δ t is using the period.Pseudo frequency should cover useful signal f_sigWhole frequency ranges, it may be assumed that

Therefore minimum frequency (the minf of useful signal_sig) it should be the pseudo frequency of L not less than scale, that is:

The Decomposition order of wavelet transform is that two are carried out to scale s into discrete sampling, therefore the corresponding Decomposition order j of L Are as follows:

2^j-1≤L≤2^j

E: wavelet decomposition is carried out to clock deviation data and obtains wavelet coefficient.

Further, the step C the following steps are included:

F: and then to every layer of wavelet coefficient WT_fMiddle number (MED) m of (a, b) and wavelet coefficient sequence adds median (MAD) the sum of n times compares, and even high frequency coefficient meets

|WT_f(a, b) | > (m+n × MAD)

It is exactly larger exceptional value that it is corresponding, which to be considered as the wavelet coefficient, while rejecting the wavelet coefficient.

Further, the step D the following steps are included:

G: every layer of generic threshold value is calculated, it may be assumed that

Wherein.λ is generic threshold value, and σ is the standard deviation of noise signal, and N is the number of wavelet coefficient.

H: soft-threshold processing, soft-threshold function expression formula are carried out to high frequency coefficient are as follows:

Wherein WT_f(a, b) is wavelet coefficient,To estimate wavelet coefficient, λ is threshold value.

Further, the step E the following steps are included:

I: wavelet reconstruction restores the clock deviation first difference data after denoising, the calculation formula of reconstruct are as follows:

Wherein, Sf (a, b) is low frequency coefficient, WT_f(a, b) is wavelet coefficient,WithCorrespond respectively to weight Structure low-pass filter and reconstruct high-pass filter.

In order to verify the validity that the present invention proposes clock deviation preprocess method, the GPS that the present embodiment is provided using the website IGS The Clock Bias data at the interval 30s carry out tentative calculation analysis in the final precise ephemeris of system.With GPS week1584 first day The PRN02 satellite clock correction data instance (other satellite datas) in (on May 17th, 2010), first does clock deviation data once Difference, PRN02 clock deviation first difference are as shown in Figure 2；Then first difference data are being utilized respectively MAD, wavelet threshold algorithms And WMAD pretreatment strategy proposed by the present invention is pre-processed, and is built respectively by BP neural network with pretreated data Mould forecasts satellite clock correction.

Satellite atomic clock normal condition be it is more stable, the numerical value change of clock deviation data is relatively small between adjacent epoch, The clock deviation data variation that PRN02 satellite clock correction data are influenced two neighboring epoch by noise from Fig. 2 is larger, carries out one Occur many peak values after secondary difference, it, can significantly if carrying out Modeling and Prediction clock deviation data with the data containing a large amount of exceptional values The precision of forecast is influenced, so it is necessary for carrying out reasonable pretreatment to clock deviation first difference data before modeling.

Fig. 3, Fig. 4 and Fig. 5 are that MAD pretreatment, wavelet threshold algorithms and WMAD of the present invention pre-process three kinds of methods respectively To the pretreated effect picture of PRN02 satellite clock correction first difference, can be evident that from Fig. 3, MAD method is to rejecting Relatively large exceptional value effect is obvious, but relatively small exceptional value is but not achieved satisfied effect, this is because clock Variable quantity is smaller after difference progress first difference and numerical value becomes smaller, and uses the median of first difference data as threshold value, often relatively Lesser exceptional value cannot be found by threshold value, therefore cannot reject relatively small exceptional value well；It can be obvious from Fig. 4 Find out, wavelet threshold algorithms can be good at falling relatively small outlier processing by frequency domain processing, but cannot be quasi- The relatively large exceptional value of true rejecting, this is because the noise variance of every layer of calculating can be made to become since larger exceptional value must influence It is small, so that the threshold value calculated can become smaller, it cannot detect and some represent larger exceptional value and obtain wavelet coefficient；And it can be with from Fig. 5 Find out, WMAD preprocess method of the invention compensates for the deficiency of both above method, i.e., effectively rejects larger exceptional value energy again Inhibit lesser noise error well, preferably purifies data.

In order to preferably verify the validity of proposition method of the present invention, the pretreated difference data of three kinds of methods is led to It crosses BP neural network modeling to predict clock deviation, using a difference data after forecast with known clock deviation data for being added just The forecast clock deviation at required moment can be obtained, finally, using the corresponding precise clock correction value of IGS as benchmark, with root-mean-square error (RMS) With the precision of the result of mean error (Mean) evaluation forecast；It is evaluated using maximum, minimal error absolute value of the difference (Range) The stability that algorithm is forecast.Wherein root-mean-square error calculation formula are as follows:

In formula, t_iIt is IGS precise clock correction；For clock deviation predicted value.

It is pre- using the pretreated difference data of no progress and MAD in order to verify the validity of WMAD preprocess method 1440 epoch before processing, wavelet threshold algorithms and the WMAD of the present invention pretreatment pretreated difference data of three kinds of methods (12h) carries out BP neural network modeling, forecasts the clock deviation of next 12 hours of 1440 epoch, the effect of comparative analysis forecast Fruit.Table 1, which gives, does not carry out pre-processing and utilize the pretreated forecast statistics value of three kinds of methods, prediction error figure such as Fig. 6 It is shown.

1 forecast result statistical form of table/ns

Using without carrying out pretreated first difference data and pre-processing a difference using three kinds of methods in contrast table 1 Clock deviation forecast statistics value can be seen that from RMS the and Mean value of forecast result by the pretreated clock deviation forecast of three kinds of methods Precision is superior to not by the forecast precision of pretreatment clock deviation, and the Range value of result then can be seen that according to weather report, by pre- Data prediction stability that treated is also superior to clock deviation after not pre-processing and forecasts stability；Through MAD and wavelet threshold algorithms Clock deviation forecast precision is substantially suitable after two methods pretreatment, forecasts that MAD method is better than wavelet threshold algorithms in stability；And it passes through All greatly better than two kinds of pretreatments of front on the pretreated clock deviation forecast precision of WMAD proposed by the present invention and stability Method.

From Fig. 6 prediction error figure, it can also be seen that, the clock deviation error of preprocessed first difference value forecast is not maximum, and have compared with Big floating；The clock of not preprocessed first difference value forecast is compared with the pretreated clock deviation prediction error of wavelet threshold through MAD Mistake difference reduce, MAD pretreatment after prediction error float it is larger, although and wavelet threshold pretreatment after prediction error float compared with The small but relatively large error robustness of appearance is not high；It is not only smaller through the pretreated clock deviation error amount of WMAD, and not compared with Big floating has good stability, demonstrates the validity of the method for the present invention.