CN106595672B - Pulsar Arrival Time Estimation Method and System Based on Anti-noise Fast Compressed Sensing - Google Patents

Abstract

The present invention provides a kind of pulsar arrival time estimation method and system based on the perception of anti-noise Fast Compression, building including redundant dictionary, based on hadamard matrix, the observation energy for defining some row vector is the product variation range of the row vector and each column vector of dictionary, according to preset observation energy threshold, if the observation energy of some row vector is greater than thresholding, which is selected as a line of observing matrix；Match tracing is realized according to the correlation of measurement vector and each dictionary column vector.Compared with the compressed sensing based pulsar arrival time estimation method of tradition, technical solution of the present invention has stronger robustness to noise, and measurement vector number is less, and is convenient for hardware realization, has important application value.

Description

Pulsar arrival time estimation method and system based on the perception of anti-noise Fast Compression

Technical field

The invention belongs to Spacecraft Autonomous Navigation field, in particular to a kind of compressed sensing based pulsar arrival time Estimation method and system.

Background technique

X-ray pulsar navigation is a kind of emerging Spacecraft Autonomous Navigation method.X-ray pulsar navigation is to utilize arteries and veins Star radiation signal is rushed to navigate.The external stable radiation of X-ray pulsar, periodically pulsing signal.These signals are pressed It is accumulated according to the period, can get stable pulse profile.The profile is compared with calibration pulse profile, arteries and veins can be obtained The fractional part of (TOA, time-of-arrival) is rushed arrival time, and its number of cycles part can then be predicted by spacecraft It is estimated position.

Pulsar arrival time is the basis that pulsar navigation system can work normally.Currently, pulsar arrival time is estimated Meter is the research hotspot in pulsar navigation field.In recent years, existing scholar starts compressed sensing being applied to pulsar signal Processing.Such as: Su Zhe et al. delivered on " Chinese science: physics mechanics astronomy " " compressed sensing based in 2011 Pulsar profile developing algorithm "；Shen Lirong et al. has delivered " A robust compressed in 2016 on " Optik " Sensing based method for X-ray Pulsar Profile construction, robust compression sensing method exist Application in pulsar signal reconstruct ".Both methods focuses mainly on pulsar signal reconstruct.In fact, pulsar signal is gone It makes an uproar, the final purpose of reconstruct is all to obtain high-precision pulse star arrival time.Li Shengliang et al. was sent out on " Optik " in 2014 Table " Fleet algorithm for X-ray pulsar profile construction and TOA solution Based on compressed sensing, compressed sensing based pulsar profile and TOA fast algorithm "；Yu Hang et al. in " A sparse representation- has been delivered on " aerospace science and technology " within 2015 based optimization algorithm for measuring the time delay of pulsar Integrated pulse profile, using the pulsar accumulated pulse profile time delay measurement method of rarefaction representation ".This Two articles estimate compressed sensing applied to pulsar arrival time, achieve preferable effect.But in measurement vector All be in selection it is random, certain measurement vectors possibly can not play a role or act on very little.In this way, a large amount of useless observation arrows Amount is selected, and useful measurement vector is then missed.

The application in other field is perceived compared to conventional compression, the compression in current PRF star arrival time estimation method Perception has following feature:

(1) each element in dictionary is not orthogonal, redundancy.Element in dictionary possesses out of phase Same pulse profile.

(2) pulsar contour signal is considered as 1 rank sparse signal.

For following condition need to be met convenient for the realization pulsar arrival time estimation method on spacecraft:

(1) amplitude and phase are two parameters of pulsar profile.Since pulsar radiation intensity is unstable, the arteries and veins of acquisition It is also unknown for rushing star-wheel exterior feature amplitude.Pulsar TOA estimation emphasis is to obtain phase.Therefore, compressed sensing based pulsar TOA algorithm for estimating should not be interfered by amplitude.

(2) for the ease of hardware realization, observing matrix coefficient need to be { 1, -1 }.

(3) in order to reduce calculation amount, measurement vector number should be reduced to the greatest extent.

But there has been no the appearance of suitable technical solution at present.

Summary of the invention

The invention proposes a kind of pulsar arrival time estimation technique schemes based on the perception of anti-noise Fast Compression, it is intended to More accurate pulse arrival time estimation is rapidly provided for navigation system.

To achieve the above object of the invention, the technical scheme is that

A kind of pulsar arrival time estimation method based on the perception of anti-noise Fast Compression, includes the following steps,

Step 1, the building of redundant dictionary, realization is as follows,

If the redundant dictionary of pulsar accumulation profile is expressed as,

Wherein, n is serial number, and n={ 0,1,2 ... N-1 }, N are pulsion phase seat interval number, and redundant dictionary ψ is N × N's Matrix；

It is the column vector of N × 1 for n-th of column vector in dictionary, expression formula is as follows,

Wherein, s is pulsar nominal contour, and i is variable, i={ 0,1,2 ... N-1 },ForI-th of element, Mod () indicates modulus；

Step 2, observing matrix is set, realization is as follows,

If H is N rank hadamard matrix,

H=[h₀；h₁；h₂；…h_n；…h_N-1]

Wherein, h_nIt is the row vector in hadamard matrix, is 1*N, n={ 0,1,2 ... N-1 }；

Define some row vector h_nObservation energy e_nFor the product variation range of the row vector and each column vector of dictionary, e_n Expression formula is as follows,

e_n=max (h_n·Ψ)-min(h_n·Ψ)

Wherein, max (h_nΨ) indicate row vector h_nWith the maximum value of the product of each column vector of dictionary, min (h_nΨ) then Indicate row vector h_nWith the minimum value of the product of each column vector of dictionary；

According to preset observation energy threshold T, if some row vector h_nObservation energy e_nGreater than T, which is selected as A line of observing matrix, otherwise gives up；

It is the matrix of m × N if observing matrix is expressed as Φ, m is the line number being selected in hadamard matrix H；

Step 3, the acquisition of measurement vector, realization is as follows,

Y=Φ x

Wherein, y is measurement vector, is the vector of m × 1；X is that pulsar accumulates profile, is 1 × N vector；

Step 4, match tracing, realization is as follows,

First find out measurement vector y and each dictionary column vectorCorrelation, be denoted asN=0,1,2 ... N- 1 }, subscript T indicates transposition；

Serial number n corresponding with the dictionary vector of observation correlation maximum_τExpression formula it is as follows,

Pulsar arrival time estimated valueExpression formula is

Moreover, the expression formula that pulsar accumulates i-th of phason interval in profile x is as follows in step 3,

X (i)=Poisson (At λ_b/N+A·t·λ_s·s(mod(i-τ·N/p,N)))

Wherein, Poisson is Poisson distribution random signal generator, λ_bFor background noise photon flux density, λ_sFor pulse Star radiated photons flux density, A are X-ray detector effective area, and t is pulsar radiation signal observation time, and τ prolongs for the time Late, p is the pulsar radiation signal period.

It is including following the present invention also provides a kind of pulsar arrival time estimating system based on the perception of anti-noise Fast Compression Module, the first module, for the building of redundant dictionary, realization is as follows,

If the redundant dictionary of pulsar accumulation profile is expressed as,

Wherein, n is serial number, and n={ 0,1,2 ... N-1 }, N are pulsion phase seat interval number, and redundant dictionary ψ is N × N's Matrix；

It is the column vector of N × 1 for n-th of column vector in dictionary, expression formula is as follows,

Wherein, s is pulsar nominal contour, and i is variable, i={ 0,1,2 ... N-1 },ForI-th of element, Mod () indicates modulus；

Second module, for setting observing matrix, realization is as follows,

If H is N rank hadamard matrix,

H=[h₀；h₁；h₂；…h_n；…h_N-1]

Wherein, h_nIt is the row vector in hadamard matrix, is 1*N, n={ 0,1,2 ... N-1 }；

Define some row vector h_nObservation energy e_nFor the product variation range of the row vector and each column vector of dictionary, e_n Expression formula is as follows,

e_n=max (h_n·Ψ)-min(h_n·Ψ)

Wherein, max (h_nΨ) indicate row vector h_nWith the maximum value of the product of each column vector of dictionary, min (h_nΨ) then Indicate row vector h_nWith the minimum value of the product of each column vector of dictionary；

According to preset observation energy threshold T, if some row vector h_nObservation energy e_nGreater than T, which is selected as A line of observing matrix, otherwise gives up；

It is the matrix of m × N if observing matrix is expressed as Φ, m is the line number being selected in hadamard matrix H；

Third module, for the acquisition of measurement vector, realization is as follows,

Y=Φ x

Wherein, y is measurement vector, is the vector of m × 1；X is that pulsar accumulates profile, is 1 × N vector；

4th module is used for match tracing, and realization is as follows,

First find out measurement vector y and each dictionary column vectorCorrelation, be denoted asN=0,1,2 ... N- 1 }, subscript T indicates transposition；

Serial number n corresponding with the dictionary vector of observation correlation maximum_τExpression formula it is as follows,

Pulsar arrival time estimated valueExpression formula is

Moreover, the expression formula that pulsar accumulates i-th of phason interval in profile x is as follows in third module,

X (i)=Poisson (At λ_b/N+A·t·λ_s·s(mod(i-τ·N/p,N)))

Wherein, Poisson is Poisson distribution random signal generator, λ_bFor background noise photon flux density, λ_sFor pulse Star radiated photons flux density, A are X-ray detector effective area, and t is pulsar radiation signal observation time, and τ prolongs for the time Late, p is the pulsar radiation signal period.

The present invention has studied under a kind of Low SNR, and the pulsar arrival time based on the perception of anti-noise Fast Compression is estimated Technical solution is counted, in compressed sensing, the pulse profile of out of phase constitutes the sparse dictionary of redundancy；Using hadamard matrix as base Plinth selects a small amount of Hadamard vectorial structure observing matrix to observe energy as criterion, and the observing matrix line number is less, anti-noise energy Power is strong；Finally, estimating pulse arrival time using match tracing.Estimate with tradition compressed sensing based pulsar arrival time Method is compared, which has stronger robustness to noise, and measurement vector number is less, and is convenient for hardware realization, In first place in the world, there is important market value.

Detailed description of the invention

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

Fig. 2 is the PSR B1937+21 calibration pulse profile schematic diagram of the embodiment of the present invention；

Fig. 3 is the observation energy diagram of the embodiment of the present invention；

Fig. 4 is the dependency diagram of the embodiment of the present invention.

Specific embodiment

Technical solution of the present invention is specifically described below in conjunction with drawings and examples.

In the prior art, a large amount of useless measurement vectors are selected, and useful measurement vector is then missed.For this Problem, invention defines observation energy, and useful measurement vector are selected according to this, and constitute observing matrix.

In embodiment, the invention will be further described by taking X-ray pulsar PSR B1937+21 as an example.Pulse profile number According to can be downloaded from European pulsar website (The European Pulsar Network).It will be to mark after the data normalization Quasi- pulse profile.Fig. 2 gives the calibration pulse profile of PSR B1937+21.

Embodiment provide process the following steps are included:

Step 1: the building of redundant dictionary.

The redundant dictionary of pulsar PSR B1937+21 accumulation profile may be expressed as:

Wherein, n is serial number, and n={ 0,1,2 ... N-1 }, N are pulsion phase seat interval number, in the present embodiment, N= 1024.Redundant dictionary ψ is the matrix of N × N.

It is the column vector of N × 1, expression formula for n-th of column vector in dictionary are as follows:

Wherein, s is pulsar nominal contour, and for waveform as shown in Fig. 2, waveform abscissa is phase, ordinate is photon Density.I is variable, i={ 0,1,2 ... N-1 },ForI-th of element, i.e. the element of n-th the i-th row of column in dictionary ψ. Mod () indicates modulus.

Step 2: the setting of observing matrix.

For the ease of hardware realization, observing matrix coefficient need to be { 1, -1 }.Hadamard (Hadamard) matrix be by+1 and- The orthogonal square matrix that 1 element is constituted.The present invention proposes based on hadamard matrix, designs observing matrix.If H is N rank Hadamard Matrix.In the present embodiment, H is 1024 rank hadamard matrixs.

H=[h₀；h₁；h₂；…h_n；…h_N-1] (3)

Wherein, h_nIt is the row vector in hadamard matrix, is 1*N, n={ 0,1,2 ... N-1 }.

Define some row vector h_nObservation energy e_nFor the product variation range of the row vector and each column vector of dictionary, e_n Expression formula is as follows:

e_n=max (h_n·Ψ)-min(h_n·Ψ) (4)

Row vector h_nH is denoted as with each column vector product of dictionary ψ_nΨ specifically includes h_nWith N number of column vector(n=0, 1,2 ... N-1 }) product is taken respectively, N number of element, max (h is obtained_nIt Ψ) indicates the maximum value in these elements, swears at once Measure h_nWith the maximum value of the product of each column vector of dictionary；min(h_nΨ) then indicate minimum value, i.e. expression row vector h_nWith dictionary The minimum value of the product of each column vector.The two difference embodies amplitude of variation.If the value does not change or amplitude of variation very little, The element in dictionary cannot be distinguished according to the observation, i.e. the value restores without effect signal.Therefore, observation energy is bigger, h_n It can more play a role in signal recovery, it more can anti-noise.Corresponding observation energy such as Fig. 3 of each row vector in hadamard matrix It is shown.

Definition observation energy threshold T extracts the Hadamard row vector of high observation energy according to thresholding.If observation energy is greater than The thresholding, h_nThe a line that may be selected as observing matrix, otherwise gives up.It is the matrix of m × N if observing matrix is represented by Φ, m is The selected line number of hadamard matrix.In general, thresholding T is smaller, line number m is more, and precision is higher, but calculation amount increases；Instead ?.

In the present embodiment, if setting thresholding T as 1.9, m 6, that is, the 128th of hadamard matrix has been selected, 192, 256,320,384,448 rows；If setting thresholding T as 1.5, m 14；If setting thresholding T as 0.5, m 62；If setting door It is limited to 0.2, then m is 116；If setting door T is limited to 0, m 1024.In the present embodiment, by emulation experiment, in calculation amount It trades off between pulsar arrival time estimated accuracy, threshold sets 0.2.When it is implemented, those skilled in the art can be certainly Row presets the value of observation energy threshold T as needed.

Step 3: the acquisition of measurement vector.

The calculating formula of measurement vector is as follows:

Y=Φ x (5)

Wherein, y is measurement vector, is the vector of m × 1.X is that pulsar accumulates profile, is 1 × N vector.

Since the arrival time of photon obeys Poisson distribution, the expression formula at i-th of phason interval is as follows in x:

X (i)=Poisson (At λ_b/N+A·t·λ_s·s(mod(i-τ·N/p,N))) (6)

Wherein, Poisson is Poisson distribution random signal generator, λ_bFor background noise photon flux density, λ_sFor pulse Star radiated photons flux density, A are X-ray detector effective area, and t is pulsar radiation signal observation time, and τ prolongs for the time Late, p is the pulsar radiation signal period.

In the present embodiment, λ_bFor 0.005ph/cm²/ s, λ_sFor 4.99*10^-5ph/cm²/ s, A 384cm², t is 1000s, p 1.5578ms, τ 0.38945ms.

Step 4: the realization of match tracing.

Since the sparse order of signal is 1, matching pursuit algorithm of the present invention is not necessarily to iteration.Measurement vector y and each can first be found out Dictionary column vectorCorrelation, be denoted asN={ 0,1,2 ... N-1 }, subscript T indicate transposition.Different dictionary arrows It measures as shown in Figure 4 with the correlation of observation.Serial number n corresponding with the dictionary vector of observation correlation maximum_τIt is as required.It is comprehensive On, n_τExpression formula it is as follows:

Since matching process is unrelated with signal amplitude value, so, pulsar profile amplitude will not generate shadow to phase estimation It rings.

Pulsar arrival time estimated valueExpression formula is as follows:

To sum up, this method meets three conditions applied to spacecraft, has preferable performance.

500 Monte Carlo simulations are carried out to technical solution of the embodiment of the present invention, error 1.4669us is corresponding Position error is 440m.

When it is implemented, method provided by the present invention can realize automatic running process based on software technology, mould can also be used Block mode realizes corresponding system.The embodiment of the present invention accordingly provides a kind of pulsar arrival based on the perception of anti-noise Fast Compression Time Estimate system, comprises the following modules,

First module, for the building of redundant dictionary, realization is as follows,

If the redundant dictionary of pulsar accumulation profile is expressed as,

Wherein, n is serial number, and n={ 0,1,2 ... N-1 }, N are pulsion phase seat interval number, and redundant dictionary ψ is N × N's Matrix；

It is the column vector of N × 1 for n-th of column vector in dictionary, expression formula is as follows,

Wherein, s is pulsar nominal contour, and i is variable, i={ 0,1,2 ... N-1 },ForI-th of element, Mod () indicates modulus；

Second module, for setting observing matrix, realization is as follows,

If H is N rank hadamard matrix,

H=[h₀；h₁；h₂；…h_n；…h_N-1]

Wherein, h_nIt is the row vector in hadamard matrix, is 1*N, n={ 0,1,2 ... N-1 }；

Define some row vector h_nObservation energy e_nFor the product variation range of the row vector and each column vector of dictionary, e_n Expression formula is as follows,

e_n=max (h_n·Ψ)-min(h_n·Ψ)

Wherein, max (h_nΨ) indicate row vector h_nWith the maximum value of the product of each column vector of dictionary, min (h_nΨ) then Indicate row vector h_nWith the minimum value of the product of each column vector of dictionary；

According to preset observation energy threshold T, if some row vector h_nObservation energy e_nGreater than T, which is selected as A line of observing matrix, otherwise gives up；

It is the matrix of m × N if observing matrix is expressed as Φ, m is the line number being selected in hadamard matrix H；

Third module, for the acquisition of measurement vector, realization is as follows,

Y=Φ x

Wherein, y is measurement vector, is the vector of m × 1；X is that pulsar accumulates profile, is 1 × N vector；

4th module is used for match tracing, and realization is as follows,

First find out measurement vector y and each dictionary column vectorCorrelation, be denoted asN=0,1,2 ... N- 1 }, subscript T indicates transposition；

Serial number n corresponding with the dictionary vector of observation correlation maximum_τExpression formula it is as follows,

Pulsar arrival time estimated valueExpression formula is

Each module specific implementation can be found in corresponding steps, and it will not go into details by the present invention.

Specific embodiment described herein is only an example for the spirit of the invention.The neck of technology belonging to the present invention The technical staff in domain can make various modifications or additions to the described embodiments or replace by a similar method In generation, however, it does not deviate from the spirit of the invention or beyond the scope of the appended claims.

Claims (4)

1. a pulsar time-of-arrival estimation method based on anti-noise fast compressed sensing, is characterized in that: comprise the following steps, Step 1, the construction of redundant dictionary is implemented as follows, Let the redundant dictionary representation of the pulsar cumulative profile be, Among them, n is the serial number, n={0,1,2,...N-1}, N is the number of pulse phase sub-intervals, and the redundancy dictionary ψ is an N×N matrix; is the nth column vector in the dictionary, which is an N×1 column vector, and the expression is as follows, Among them, s is the standard profile of the pulsar, i is a variable, i={0,1,2,...N-1}, for The i-th element of , mod() means modulo; Step 2, set the observation matrix, the implementation is as follows, Let H be a Hadamard matrix of order N, H=[h ₀ ; h ₁ ; h ₂ ; ... h _n ; ... h _N-1 ] Among them, h _n is the row vector in the Hadamard matrix, which is 1*N, n={0,1,2,...N-1}; Define the observed energy e _n of a row vector h _n as the variation range of the product of the row vector and each column vector of the dictionary. The expression of e _n is as follows: e _n =max(h _n ·Ψ)-min(h _n ·Ψ) Wherein, max(h _n ·Ψ) represents the maximum value of the product of the row vector h _n and each column vector of the dictionary, and min(h _n ·Ψ) represents the minimum value of the product of the row vector h _n and each column vector of the dictionary; According to the preset observation energy threshold T, if the observation energy e _n of a row vector h _n is greater than T, the row vector is selected as a row of the observation matrix, otherwise it is discarded; Let the observation matrix be expressed as Φ, which is an m×N matrix, and m is the number of selected rows in the Hadamard matrix H; Step 3, the acquisition of the observation vector is implemented as follows, y=Φ·x Among them, y is the observation vector, which is an m×1 vector; x is the cumulative contour of the pulsar, which is a 1×N vector; Step 4, matching tracking, is implemented as follows, First find the observation vector y and each dictionary column vector correlation, denoted as n={0,1,2,...N-1}, the superscript T means transposition; The expression of the sequence number n _τ corresponding to the dictionary vector with the greatest correlation with the observation value is as follows, Pulsar Arrival Time Estimation The expression is p is the period of the pulsar radiation signal. 2. the pulsar time-of-arrival estimation method based on anti-noise fast compressed sensing according to claim 1, is characterized in that: in step 3, the expression of the i-th phase sub-interval in the pulsar accumulation profile x is as follows, x(i)=Poisson(A·t·λ _b /N+A·t·λ _s ·s(mod(i-τ·N/p,N))) where Poisson is a Poisson distributed random signal generator, λ _b is the background noise photon flux density, λ _s is the pulsar radiation photon flux density, A is the effective area of the X-ray detector, t is the pulsar radiation signal observation time, τ is the time delay, and p is the period of the pulsar radiation signal. 3. A pulsar time-of-arrival estimation system based on anti-noise fast compressed sensing, characterized in that: comprising the following modules, The first module, used for the construction of redundant dictionary, is implemented as follows, Let the redundant dictionary representation of the pulsar cumulative profile be, Among them, n is the serial number, n={0,1,2,...N-1}, N is the number of pulse phase sub-intervals, and the redundancy dictionary ψ is an N×N matrix; is the nth column vector in the dictionary, which is an N×1 column vector, and the expression is as follows, Among them, s is the standard profile of the pulsar, i is a variable, i={0,1,2,...N-1}, for The i-th element of , mod() means modulo; The second module, used to set the observation matrix, is implemented as follows, Let H be a Hadamard matrix of order N, H=[h ₀ ; h ₁ ; h ₂ ; ... h _n ; ... h _N-1 ] Among them, h _n is the row vector in the Hadamard matrix, which is 1*N, n={0,1,2,...N-1}; Define the observed energy e _n of a row vector h _n as the variation range of the product of the row vector and each column vector of the dictionary. The expression of e _n is as follows: e _n =max(h _n ·Ψ)-min(h _n ·Ψ) Wherein, max(h _n ·Ψ) represents the maximum value of the product of the row vector h _n and each column vector of the dictionary, and min(h _n ·Ψ) represents the minimum value of the product of the row vector h _n and each column vector of the dictionary; According to the preset observation energy threshold T, if the observation energy e _n of a row vector h _n is greater than T, the row vector is selected as a row of the observation matrix, otherwise it is discarded; Let the observation matrix be expressed as Φ, which is an m×N matrix, and m is the number of selected rows in the Hadamard matrix H; The third module, used to obtain the observation vector, is implemented as follows, y=Φ·x Among them, y is the observation vector, which is an m×1 vector; x is the cumulative outline of the pulsar, which is a 1×N vector; The fourth module, used for matching tracking, is implemented as follows, First find the observation vector y and each dictionary column vector correlation, denoted as n={0,1,2,...N-1}, superscript T means transposition; The expression of the sequence number n _τ corresponding to the dictionary vector with the greatest correlation with the observation value is as follows, Pulsar Arrival Time Estimation The expression is p is the period of the pulsar radiation signal. 4. The pulsar time-of-arrival estimation system based on anti-noise fast compressed sensing according to claim 3, is characterized in that: in the third module, the expression of the ith phase sub-interval in the pulsar accumulation profile x is as follows, x(i)=Poisson(A·t·λ _b /N+A·t·λ _s ·s(mod(i-τ·N/p,N))) where Poisson is a Poisson distributed random signal generator, λ _b is the background noise photon flux density, λ _s is the pulsar radiation photon flux density, A is the effective area of the X-ray detector, t is the pulsar radiation signal observation time, τ is the time delay, and p is the period of the pulsar radiation signal.