CN106097264A - Based on dual-tree complex wavelet and morphologic satellite telemetering data filtering method - Google Patents

Abstract

The invention discloses a kind of based on dual-tree complex wavelet and morphologic satellite telemetering data filtering method, for solving the technical problem of existing satellite telemetering data filtering method poor practicability.Technical scheme is to decompose first with dual-tree complex wavelet telemetry original to satellite, utilizes a kind of semi-soft threshold filter function to be filtered initial telemetry, effectively filters out the noise in data；Next utilizes complex morphological filtering method to be filtered primary data, and data filtering retains wherein abnormal data simultaneously；Finally ask for the difference of two kinds of filter result and carry out comprehensively, obtaining final filter result to two kinds of filter result by setting threshold value.Test result indicate that the noise in telemetry can be effectively filtered out and retain abnormal data by the inventive method, practical.

Description

Based on dual-tree complex wavelet and morphologic satellite telemetering data filtering method

Technical field

The present invention relates to a kind of satellite telemetering data filtering method, particularly to one based on dual-tree complex wavelet and morphology Satellite telemetering data filtering method.

Background technology

Document 1 " based on improving the dual-tree complex wavelet electro-ocular signal denoising method of threshold value, observation and control technology, 2015, Vol34 (8), p16-p18 " a kind of new threshold filter method based on dual-tree complex wavelet is disclosed, the method improve uniform threshold and excessively grip The problem killing wavelet coefficient, it is possible to suppressing noise the most preferably stick signal details, described in document, method is to electro-ocular signal Having better effects, method application signal is limited, and when there is significantly ANOMALOUS VARIATIONS in signal, calculated threshold value is to whole in signal Bulk noise filtration result declines.

Document 2 " electrocardiosignal based on shape filtering removes baseline drift method, Acta Physica Sinica, and 2014, Vol63 (9), P098701-1～p098701-6 " disclose a kind of ECG's data compression method based on two-stage morphological filter.The method is transported Electrocardiosignal feature is combined, it is proposed that based on morphologic ECG baseline drift minimizing technology, pass through with morphology theory The structural element using difformity and size carries out opening and closing, make and break cascading computing to signal, through checking the method energy Enough electrocardiosignal feature is preferably kept, improve Signal-to-Noise, reduce mean square deviation and also eliminate making an uproar of baseline drift Sound.Described in document, method is set by experiment, has better effects to electrocardiosignal, but uses two-stage shape for other types signal State filtering method can cause mutation amplitude characteristic in initial data to decay, unfavorable to abnormality detection.

Summary of the invention

In order to overcome the deficiency of existing satellite telemetering data filtering method poor practicability, the present invention provides a kind of based on double trees Phase information and morphologic satellite telemetering data filtering method.The method is first with dual-tree complex wavelet remote measurement original to satellite number According to decomposing, utilize a kind of semi-soft threshold filter function that initial telemetry is filtered, the noise in data is carried out Effectively filter out；Next utilizes complex morphological filtering method to be filtered primary data, retains data filtering the most different simultaneously Regular data；Finally ask for the difference of two kinds of filter result and carry out comprehensively, obtaining to two kinds of filter result by setting threshold value Whole filter result.Test result indicate that the noise in telemetry can be effectively filtered out and to abnormal number by the inventive method According to retaining, practical.

The technical solution adopted for the present invention to solve the technical problems: a kind of based on dual-tree complex wavelet and morphologic satellite Telemetry filtering method, is characterized in comprising the following steps:

Step one, utilize dual-tree complex wavelet that original telemetry X of satellite of N a length of when collecting is decomposed into M shell, obtain The individual-layer data x of different frequency bands_j,i(i=1,2 ..., N；J=1,2 ..., M), use semisoft shrinkage function T_j(x_j,i,α_j) logarithm According to being filtered, concrete formula is as follows:

$T_j\left(x_{j,i}{\mathrm{\α}}_j\right)=\left\{\begin{array}{c}{x}_{j,i}-sign\left(x_{j,i}\right)*\frac{{\mathrm{th}}_j}{{\mathrm{\α}}_j*\left|x_{j,i}\right|},\left|x_{j,i}\right|>{\mathrm{th}}_j\\ sign\left(x_{j,i}\right)*\frac{\left|x_{j,i}\right|}{{\mathrm{\α}}_j*{\mathrm{th}}_j},\left|x_{j,i}\right|\≤{\mathrm{th}}_j\end{array}\right.---\left(1\right)$

In formula, x_jFor input value, j is dual-tree complex wavelet decomposition scale, sign for taking sign function, α_jFor threshold value regulation it is Number, value is the biggest to be filtered closer to hard-threshold, is suitable for processing noise data, otherwise then closer to soft-threshold de-noising, to critical zone Wavelet coefficient retain, keep the singular data in legacy data, the therefore size of penalty coefficient and decomposition scale intermediate frequency Rate is inversely proportional to.th_jFor adaptive thresholding value function, its computing formula is:

${\mathrm{th}}_j={\mathrm{\σ}}_j\sqrt{2*\ln (N/2^j)}---\left(2\right)$

In formula, σ_jFor j layer x_jVariance, N is data length.Formula (2) is brought into formula (1) and obtains semi-soft threshold filter public affairs Formula is:

$T_j(x_{j,i},{\mathrm{\α}}_j,{\mathrm{\σ}}_j,N)=\left\{\begin{array}{c}{x}_{j,i}-sign\left(x_{j,i}\right)*\frac{{\mathrm{\σ}}_j\sqrt{2*\ln (N/2^j)}}{{\mathrm{\α}}_j*\left|x_{j,i}\right|},\left|x_{j,i}\right|>{\mathrm{th}}_j\\ sign\left(x_{j,i}\right)*\frac{{\mathrm{\σ}}_j\sqrt{2*\ln (N/2^j)}}{{\mathrm{\α}}_j*{\mathrm{th}}_j},\left|x_{j,i}\right|\≤{\mathrm{th}}_j\end{array}\right.---\left(3\right)$

Filtered data are reconstructed and obtain dual-tree complex wavelet semi-soft threshold filter result XD_i(i=1,2 ..., N).

Initial telemetry X is filtered by step 2, employing morphological method, and step is as follows:

The data shape choice structure element kept the most as required, and determine structural element parameter, selected structural elements Element should be with primary data form as reference, and parameter selects to be inversely proportional to sample frequency, or selects knot according to final filter effect Constitutive element；

2. determining that Morphological Filtering Algorithm, shape filtering rudimentary algorithm include erosion operation and dilation operation, concrete formula is such as Under:

Erosion operation Θ formula is:

$X\mathrm{\Θ}B=\{x_i:B+x_i\⋐X\}---\left(4\right)$

In formula, X is input data, x_iFor data value, B is structural element, utilizes structural element B sliding in input data X Dynamic, the data value x to each position_iCarrying out computing, structural element B has three kinds of results:

B Yu X maximal correlation；

B with X is unrelated；

B [x] ∩ X Yu B [x] ∩ X^cNot being the most empty, B is relevant to X section.

Dilation operation is the inverse operation of erosion operation, therefore dilation operationFormula is:

$X\⊕B={\[X^c\mathrm{\Θ}(-B)\]}^c---\left(5\right)$

Algorithmic procedure is identical with erosion operation with acquirement result.

In order to retain mutation data characteristics while satellite telemetering data is filtered, on the basis of basal morphological budget On, use open-close complex morphological filtering algorithm that data are processed, concrete formula is:

OC (X)=(X ο B) B (6)

It is filtered adopting by original telemetry, obtains based on complex morphological filter result XS_i(i=1,2 ..., N).

Step 3, ask for two kinds of filter result XD_iWith XS_iDifference e_i(i=1,2 ... N), by setting difference detection threshold Value th_eTo e_iDetect, by e_iIn less than threshold value th_eData acquisition dual-tree complex wavelet filter result XD_i, by e_iMore than detection Threshold value th_eData bit use shape filtering result XS_i, the filter result obtained comprehensively is obtained final filtering dataConcrete formula is:

$\left\{\begin{array}{cc}\widetilde{X_i}={\mathrm{XD}}_i,& if{e}_i<{\mathrm{th}}_e\\ \widetilde{X_i}={\mathrm{XS}}_i,& if{e}_i\&GreaterEqual;{\mathrm{th}}_e\end{array}\right.$

Threshold value th_eAccording to e_iMean square deviation is set.

The invention has the beneficial effects as follows: the inventive method is carried out first with dual-tree complex wavelet telemetry original to satellite Decompose, utilize a kind of semi-soft threshold filter function that initial telemetry is filtered, the noise in data is effectively filtered Remove；Next utilizes complex morphological filtering method to be filtered primary data, and data filtering retains wherein abnormal data simultaneously； Finally ask for the difference of two kinds of filter result and carry out comprehensively, finally being filtered knot to two kinds of filter result by setting threshold value Really.Test result indicate that the noise in telemetry can be effectively filtered out and protect abnormal data by the inventive method Stay, practical.

With detailed description of the invention, the present invention is elaborated below in conjunction with the accompanying drawings.

Accompanying drawing explanation

Fig. 1 is that the present invention is based on dual-tree complex wavelet and the flow chart of morphologic satellite telemetering data filtering method.

Fig. 2 is the initial telemetry of satellite.

Fig. 3 is to use dual-tree complex wavelet to satellite telemetering data decomposition result.

Fig. 4 is semi-soft threshold filter result.

Fig. 5 is to use OC shape filtering result.

Fig. 6 is filter result of the present invention.

Fig. 7 is " db4 " small echo penalty threshold filter result.

Detailed description of the invention

With reference to Fig. 1-7.The present invention is based on dual-tree complex wavelet and morphologic satellite telemetering data filtering method concrete steps As follows:

1, even numbers Phase information is used to be reconstructed after data being decomposed and utilizing semi-soft threshold filter.

Satellite telemetry primary data utilizes following formula to be simulated:

X=0.08*t+10*sin* (π t/15)+noise (1)

In formula, time t/length is 1000, and noise is the random number that (-1,1) is interval, adds 5 amplitude mutations in data Point.

Utilize dual-tree complex wavelet that the initial data X that length N is 1000 is decomposed into 4 layers, obtain decomposition data x_i,j(i=1, 2,…,1000；J=1,2,3,4,5).By data, calculate each variances sigma organizing data after decomposing_j, calculate the adaptive of each layer data Answer threshold function table th_j, computing formula is:

${\mathrm{th}}_j={\mathrm{\&sigma;}}_j\sqrt{2*\ln (N/2^j)}---\left(2\right)$

Obtain semi-soft threshold filter function T_jFor:

$T_j(x_{j,i},{\mathrm{\&alpha;}}_j,{\mathrm{\&sigma;}}_j,N)=\left\{\begin{array}{c}{x}_{j,i}-sign\left(x_{j,i}\right)*\frac{{\mathrm{\&sigma;}}_j\sqrt{2*\ln (N/2^j)}}{{\mathrm{\&alpha;}}_j*\left|x_{j,i}\right|},\left|x_{j,i}\right|>{\mathrm{th}}_j\\ sign\left(x_{j,i}\right)*\frac{{\mathrm{\&sigma;}}_j\sqrt{2*\ln (N/2^j)}}{{\mathrm{\&alpha;}}_j*{\mathrm{th}}_j},\left|x_{j,i}\right|\&le;{\mathrm{th}}_j\end{array}\right.---\left(3\right)$

α in formula_jFor each layer threshold value adjustment factor, it is gradually lowered owing to even numbers Phase information decomposes the frequency ranges of data obtained, because of This penalty coefficient value is inversely proportional to Decomposition order.

Semi-soft threshold filter function is utilized to be filtered initial telemetry then filtered data being reconstructed, Obtain filter result XD_i。

2, use complex morphological method that initial telemetry X is filtered.

Owing to initial telemetry main body is sinusoidal pattern, therefore use sinusoidal structured element, calculate primary data differential, root Determine that data structure element is 2 according to differential value.

Using open-close shape filtering method to be filtered initial telemetry X, computing formula is:

OC (X)=(X ο B) B (4)

In formula, structural element B is sinusoidal structured, and opening operation ο calculates and closed operation computing formula is:

$X\mathrm{\&Theta;}B=\{x_i:B+x_i\&Subset;X\}---\left(5\right)$

$X\&CirclePlus;B={\&lsqb;X^c\mathrm{\&Theta;}(-B)\&rsqb;}^c---\left(6\right)$

Utilizing complex morphological to be filtered satellite primary data, obtaining filter result is XS_i。

3, use difference detection threshold value that two kinds of filter result are carried out comprehensively.

Calculate two kinds of filter result difference e_i, the most also refer to e_iMean set difference detection threshold value th_e, This example sets th_eIt is 0.9.

e_i=| XD_i-XS_i| (i=1,2 ..., 1000) (7)

Utilize threshold value th_eThe difference obtained is carried out detection and obtains final filtering dataComputing formula is:

$\left\{\begin{array}{cc}\widetilde{X_i}={\mathrm{XD}}_i,& if{e}_i<{\mathrm{th}}_e\\ \widetilde{X_i}={\mathrm{XS}}_i,& if{e}_i\&GreaterEqual;{\mathrm{th}}_e\end{array}\right.---\left(8\right)$

Work as difference e_iLess than setting threshold value th_eTime, use dual-tree complex wavelet filter result XD_i, work as difference e_iMore than or equal to setting Determine threshold value th_eTime, use complex morphological filter result XS_i, obtain the final filter result of satellite telemetering data by threshold test

For control methods effectiveness, " db4 " small echo penalty threshold filter is used to be filtered obtaining to primary data Filter resultUtilize signal ratio (SNR) and mean square deviation (RMS) that involved 4 kinds of filtering methods are contrasted, be shown in Table 1.

Table 1 present document relates to method filter result

Find that shape filtering result is optimum to primary data filter effect by contrast, the inventive method filter result index It is between dual-tree complex wavelet filtering method and morphologic filtering method result, but noise data is effectively being pressed down by this method While system, it is effectively retained primary data exceedes the mutation data setting detection threshold value.It addition, contrast additive method The filter result of mutation data, this method Gibbs' effect to producing during wavelet filtering is inhibited.

Claims (1)

1. one kind based on dual-tree complex wavelet and morphologic satellite telemetering data filtering method, it is characterised in that include following step Rapid:

Step one, utilize dual-tree complex wavelet that original telemetry X of satellite of N a length of when collecting is decomposed into M shell, obtain difference The individual-layer data x of frequency band_j,i(i=1,2 ..., N；J=1,2 ..., M), use semisoft shrinkage function T_j(x_j,i,α_j) data are entered Row filtering, concrete formula is as follows:

$T_j(x_{j,i},{\mathrm{\&alpha;}}_j)=\left\{\begin{array}{c}{x}_{j,i}-sign\left(x_{j,i}\right)*\frac{{\mathrm{th}}_j}{{\mathrm{\&alpha;}}_j*\left|x_{j,i}\right|},\left|x_{j,i}\right|>{\mathrm{th}}_j\\ sign\left(x_{j,i}\right)*\frac{\left|x_{j,i}\right|}{{\mathrm{\&alpha;}}_j*{\mathrm{th}}_j},\left|x_{j,i}\right|\&le;{\mathrm{th}}_j\end{array}\right.---\left(1\right)$

In formula, x_jFor input value, j is dual-tree complex wavelet decomposition scale, sign for taking sign function, α_jFor threshold value adjustment factor, take It is worth the biggest closer to hard-threshold filtering, is suitable for processing noise data, otherwise then closer to soft-threshold de-noising, the small echo to critical zone Coefficient retains, and keeps the singular data in legacy data, and therefore the size of penalty coefficient becomes anti-with decomposition scale medium frequency Ratio；th_jFor adaptive thresholding value function, its computing formula is:

${\mathrm{th}}_j={\mathrm{\&sigma;}}_j\sqrt{2*\ln (N/2^j)}---\left(2\right)$

In formula, σ_jFor j layer x_jVariance, N is data length；Bring formula (2) into formula (1) and obtain semi-soft threshold filter formula and be:

$T_j(x_{j,i},{\mathrm{\&alpha;}}_j,{\mathrm{\&sigma;}}_j,N)=\left\{\begin{array}{c}{x}_{j,i}-sign\left(x_{j,i}\right)*\frac{{\mathrm{\&sigma;}}_j\sqrt{2*ln(N/2^j)}}{{\mathrm{\&alpha;}}_j*\left|x_{j,i}\right|},\left|x_{j,i}\right|>{\mathrm{th}}_j\\ sign\left(x_{j,i}\right)*\frac{{\mathrm{\&sigma;}}_j\sqrt{2*ln(N/2^j)}}{{\mathrm{\&alpha;}}_j*{\mathrm{th}}_j},\left|x_{j,i}\right|\&le;{\mathrm{th}}_j\end{array}\right.---\left(3\right)$

Filtered data are reconstructed and obtain dual-tree complex wavelet semi-soft threshold filter result XD_i(i=1,2 ..., N)；

Initial telemetry X is filtered by step 2, employing morphological method, and step is as follows:

The data shape choice structure element kept the most as required, and determine structural element parameter, selected structural element should With primary data form as reference, parameter selects to be inversely proportional to sample frequency, or according to final filter effect choice structure unit Element；

2. determining that Morphological Filtering Algorithm, shape filtering rudimentary algorithm include erosion operation and dilation operation, concrete formula is as follows:

Erosion operation Θ formula is:

$X\mathrm{\&Theta;}B=\{x_i:B+x_i\&Subset;X\}---\left(4\right)$

In formula, X is input data, x_iFor data value, B is structural element, utilizes structural element B to slide in input data X, right The data value x of each position_iCarrying out computing, structural element B has three kinds of results:

B Yu X maximal correlation；

B with X is unrelated；

B [x] ∩ X Yu B [x] ∩ X^cNot being the most empty, B is relevant to X section；

Dilation operation is the inverse operation of erosion operation, therefore dilation operationFormula is:

$X\&CirclePlus;B={\&lsqb;X^c\mathrm{\&Theta;}\left(-B\right)\&rsqb;}^c---\left(5\right)$

Algorithmic procedure is identical with erosion operation with acquirement result；

In order to retain mutation data characteristics while satellite telemetering data is filtered, on the basis of basal morphological budget, adopt Processing data with open-close complex morphological filtering algorithm, concrete formula is:

OC (X)=(X o B) B (6)

It is filtered adopting by original telemetry, obtains based on complex morphological filter result XS_i(i=1,2 ..., N)；

Step 3, ask for two kinds of filter result XD_iWith XS_iDifference e_i(i=1,2 ... N), by setting difference detection threshold value th_e To e_iDetect, by e_iIn less than threshold value th_eData acquisition dual-tree complex wavelet filter result XD_i, by e_iMore than detection threshold value th_eData bit use shape filtering result XS_i, the filter result obtained comprehensively is obtained final filtering dataTool Body formula is:

$\left\{\begin{array}{ccc}\widetilde{X_i}={\mathrm{XD}}_i,& if& e_i<{\mathrm{th}}_e\\ \widetilde{X_i}={\mathrm{XS}}_i,& if& e_i\&GreaterEqual;{\mathrm{th}}_e\end{array}\right.$

Threshold value th_eAccording to e_iMean square deviation is set.