CN104457745A - Method for estimating X-ray pulsar phase based on loop circuit tracking - Google Patents

Abstract

The invention discloses a method for estimating the X-ray pulsar phase based on loop circuit tracking. The method comprises the following steps: (1) a pulsar is determined and observed, and the pulse period P0 and the standard outline s(n) of the pulsar under the static condition are recorded; (2) within the observation time quantum, the arrival time of all the observed pulsar quantum photons is recorded; (3) the quantum photon arrival time series within the first tau s duration are superimposed to generate a measured pulse outline c(n) according to the pulse period P0 under the static condition; (4) the phase difference estimated value phi- of s(n) and c(n) is evaluated; (5) a one-order filter and a three-order filter are selected and combined to perform filtration on the phase difference estimated value phi- to obtain the frequency and phase feedback update information delta f and delta phi; (6) the delta f is utilized to correct the pulse period P0 under the static condition to obtain the corrected pulse period P- and obtain a new measured pulse outline c-(n) within the next tau s duration; (7) the delta phi is utilized to correct the standard outline s(n) to obtain the corrected standard pulse outline s-(n); (8) according to the s-(n) and the c-(n), the new phase difference estimated value phi-^ is obtained.

Description

Based on the X-ray pulsar phase estimation method of loop tracks

Technical field

The invention belongs to Spacecraft Autonomous Navigation field, relate to a kind of dynamic effect modification method of X-ray pulsar navigation signal processing procedure, be applicable to near-earth orbit, deep space and interplanetary flight spacecraft, and without the high-precision independent navigation and vehicle controL of dense atmosphere celestial body lander and surperficial stroller thereof.

Background technology

X-ray pulsar is a kind of natural beacon, the X-ray pulse signal of energy stable radiation.Utilizing X-ray pulsar to navigate is a kind of novel astronomical navigation method, is independent navigation truly.Observed by the timing of pulse signals, high-precision ranging information and temporal information can be obtained, by the imaging of paired pulses star, also can obtain attitude information.Pulsar navigation is applicable to whole solar autonomous navigation technology, the accommodation of carrying out navigation algorithm is only needed for task different phase, navigation accuracy can not reduce with the increase of flying distance, and this technology will be the most effective approach realizing seamless navigation and precise hard_drawn tuhes in survey of deep space and interplanetary mission future.

The ultimate principle of X-ray pulsar navigation is the impulse phase difference by the impulse phase and reference point (solar system barycenter) position measuring synchronization spacecraft, by certain navigation algorithm, obtain the position of observation moment spacecraft relative datum point.As can be seen here, pulsar phase estimation method is one of key link of X-ray pulsar navigation, only carries out accurate impulse phase measurement, just can carry out follow-up spacecraft navigator fix.The measuring method of impulse phase utilizes to observe the information obtained compare with the full sized pules profile recorded in advance, to obtain the phase differential of the two.Be photon time of arrival but not direct pulse information due to pulsar navigation X-ray detector record, therefore need to carry out corresponding pulsar navigation signal transacting link, to realize the conversion estimated to impulse phase time of arrival by photon.But for spacecraft, due to the Doppler effect and relativistic effect that cause of himself moving, when carrying out pulsar navigation signal transacting, the recurrence interval constantly changes, when making to adopt full sized pules profile to carry out phase estimation, there will be the situation that estimated accuracy declines.

Summary of the invention

The technical matters that the present invention solves is: overcome existing pulsar navigation signal transacting link utilize fixing full sized pules profile to realize estimated accuracy that impulse phase estimates to cause is degenerated, provide a kind of method utilizing feedback control loop acquisition pulse profile mechanical periodicity update information to realize dynamic spacecraft impulse phase accurately to estimate.

Technical solution of the present invention is: based on the X-ray pulsar phase estimation method of loop tracks, comprise the steps:

(1) determine to observe pulsar, under record quiescent conditions, observe the recurrence interval P of pulsar ₀with full sized pules profile s (n); Wherein full sized pules profile is that the mathematical expectation of pulsar photon arriving amt in one-period is at the n-th sampled point the sampled value at place, N is total number of sample points, n=1,2 ..., N;

(2) at observation time section τ ₀to τ ₀+ τ _obsin, record the time of arrival of all observation pulsar photons, wherein τ ₀for observation time starting point, τ _obsfor observing total duration;

(3) for first τ _sobservation pulsar photon sequence τ time of arrival in duration _k∈ [τ ₀, τ ₀+ τ _s], according to quiescent conditions lower recurrence interval P ₀fold, generate ranging pulse profile c (n),

c(n)＝αs(n+φ×N)+β+ω(n)

Wherein α is scale factor, and β is shift factor, and ω (n) is measurement noises, and φ is the phase differential between ranging pulse profile and full sized pules profile;

(4) according to full sized pules profile s (n) and ranging pulse profile c (n), the estimated value of the phase difference of the two is obtained

(5) select single order and three rank filter banks, aim at the phase difference estimation value of pulse profile s (n) and ranging pulse profile c (n) carry out filtering, obtain frequency feedback update information Δ f and phase feedback update information Δ φ; Wherein filter bank is expressed as:

$\mathrm{\Δf}=((a_{31}\widetilde{\mathrm{\φ}}+\∫(a_{32}\widetilde{\mathrm{\φ}\∫}+a_{33}\widetilde{\mathrm{\φ}}))+a_{11}\widetilde{\mathrm{\φ}})/2$

Δφ＝∫Δf

Wherein, Δ φ, Δ f are respectively the feedback modifiers information of phase place and frequency, a ₃₁, a ₃₂, a ₃₃represent each integral element amplification coefficient of three rank wave filters respectively, a ₁₁represent firstorder filter amplification coefficient; Integral operation in described filter bank adopts square wave digital integrator to realize;

(6) the recurrence interval P under utilizing frequency feedback update information Δ f to revise quiescent conditions ₀, obtain the revised recurrence interval and utilize to next τ _sin duration, observation pulsar photon sequence time of arrival folds, and obtains next τ _sin duration, observation pulsar photon sequence time of arrival, obtains new ranging pulse profile thus

(7) utilize phase feedback update information Δ φ to revise full sized pules profile s (n), obtain revised full sized pules profile $\hat{s}\left(n\right)=s\left(\mathrm{mod}(n+\frac{\mathrm{\Δ\φ}\×N}{\hat{P}},N)\right);$

(8) basis with obtain new phase difference estimation value

(9) for observation time section τ ₀to τ ₀+ τ _obsinterior each process step-length τ _s, repeat step (3) ~ step (8), until all observation pulsar photon series processing time of arrival is complete, obtain final phase difference estimation value.

The present invention's advantage is compared with prior art:

(1) the inventive method can realize the stable observation of paired pulses star under spacecraft current intelligence.Due to the motion of spacecraft self, the pulse signal received will introduce Doppler frequency change, observe the recurrence interval obtained inconsistent under the recurrence interval that actual observation is obtained and quiescent conditions.The inventive method utilizes loop filtering feedback compensation structure, can before carrying out profile and be folding, and corrected impulse periodic quantity, ensures that profile folds the correctness of result.

(2) the present invention can improve the precision that impulse phase is measured.When carrying out impulse phase and measuring, due to the uncertainty of photon time of arrival, impulse phase measurement result is caused to there is randomness.The inventive method utilizes phase feedback track loop to produce filter effect, and the stochastic error comprised in paired pulses phase measurement carries out filtering, improves the measuring accuracy of impulse phase.

Accompanying drawing explanation

Fig. 1 is the FB(flow block) of the inventive method;

Fig. 2 is the Laplace transformation schematic diagram of wave filter corresponding in step 4 of the present invention;

Fig. 3 is the change schematic diagram adopting phase measurement error after phase-tracking method correction of the present invention in the embodiment of the present invention.

Embodiment

As shown in Figure 1, be the process flow diagram of the inventive method, key step is as follows:

Step one: determine observation pulsar and its recurrence interval P under recording quiescent conditions ₀with nominal contour s (n), wherein nominal contour is that the mathematical expectation of pulsar photon arriving amt in one-period is at the n-th sampled point the sampled value at place, N is total number of sample points, n=1,2 ..., N.

Assuming that observation time starting point is τ ₀, process step-length τ _s, observe total duration τ _obs.(τ in observation time section ₀to τ ₀+ τ _obs), record all pulsar photons time of arrival.For a kth photon, recording its time of arrival is τ _k.

Step 2: to first group of photon sequence τ time of arrival _k∈ [τ ₀, τ ₀+ τ _s], according to quiescent conditions lower recurrence interval P ₀, access time sequence mid point τ ₀+ τ _s/ 2 as folding reference point, and choosing total number of sample points is N, generates ranging pulse profile c (n), namely carries out pulse profile and folds.Due to the alternate position spike of spacecraft and navigational coordinate system initial point, measure between profile and nominal contour and there is phase shifts.Meanwhile, due to the existence of observation noise, measure the proportional zoom and displacement and measurement noises that also there is amplitude between profile and nominal contour.Therefore the relation measuring profile can be expressed as:

c(n)＝αs(n+φ×N)+β+ω(n)

Wherein, α is scale factor, and β is shift factor, and ω (n) is measurement noises, embodies observation noise to the impact of measuring profile; φ measures the phase differential between profile and nominal contour, and for pulsar navigation, other resolve link use.

Step 3: according to nominal contour s (n) and measurement profile c (n), obtain the estimated value of phase difference .

Step 4: utilize wave filter to phase difference estimation value carry out filtering.

Wave filter is made up of multiple integral.The present invention selects single order and three rank filter banks, and its formula can be written as:

$\mathrm{\Δf}=((a_{31}\widetilde{\mathrm{\φ}}+\∫(a_{32}\widetilde{\mathrm{\φ}\∫}+a_{33}\widetilde{\mathrm{\φ}}))+a_{11}\widetilde{\mathrm{\φ}})/2$

Δφ＝∫Δf

Wherein, Δ φ, Δ f are respectively the feedback modifiers information of phase place and frequency.A ₃₁, a ₃₂, a ₃₃represent each integral element amplification coefficient of three rank wave filters respectively, a ₁₁represent firstorder filter amplification coefficient. with represent single order and three rank filtering computation structures respectively, form complete filter function by feedback control loop, suppress the high fdrequency component in feedback modifiers information, reduce the noise of feedback modifiers information.According to feedback control loop mathematical feature, a can be obtained ₁₁, a ₃₁, a ₃₂and a ₃₃optimal value, as shown in Figure 2.

For carrying out integral and calculating under discrete conditions, square wave digital integrator can be adopted, namely by carrying out y (t)=y (t-Δ t)+Δ t × x (t) recursive operation, simulation y (t)=∫ x (t).

Step 5: the recurrence interval P0 under utilizing the frequency information Δ f of feedback to revise quiescent conditions.According to the relation of cycle and frequency, can obtain namely

$\hat{P}={(P_0^{-1}-\mathrm{\Δf})}^{-1}$

By above operation, according to upper one group of photon sequence time of arrival observation information, by the feedback compensation to the cycle, reduce the revised recurrence interval with the difference in true pulse cycle.Utilize the recurrence interval after upgrading next group photon sequence time of arrival is folded, and obtains next group ranging pulse profile c (n), improve the similarity of photon sequence time of arrival and true pulse profile simultaneously.

Step 6: utilize the phase information Δ φ of feedback to revise nominal contour s (n), that is:

$\hat{s}\left(n\right)=s\left(\mathrm{mod}(n+\frac{\mathrm{\Δ\φ}\×N}{\hat{P}},N)\right);$

The object of this correction carries out translation between primary standard profile s (n), makes its phase value approach next group ranging pulse profile c (n), thus make the phase difference estimation value of output maintain the state being tending towards 0, realize revising nominal contour to the tracking of c (n).Now phase feedback update information can be reflected as the phase differential of s (n) and c (n), i.e. phase measurement.

Arrive the change of probability because nominal contour s (n) characterizes photon in one-period, variable n variation range is 1,2 ..., N.When carrying out phase information correction, need to consider to realize ring shift by MOD function mod (), describe phase shifts to the impact in the whole cycle.

Step 7: according to the nominal contour after renewal measurement profile c (n) that current photon sequence time of arrival obtained with step 5 generates, obtains new phase difference estimation value

Step 8: repeat step 4 to seven, until all photon sequence time of arrival calculates complete.Also namely with τ _sfor process step-length, process total duration τ _obs.

Embodiment

Step one: determine observation pulsar and its recurrence interval P under recording quiescent conditions ₀with nominal contour s (n).The pulsar that this example chooses observation is PSR B0531+21, and the recurrence interval under its quiescent conditions is P ₀=33.085ms.Observation starting point is τ ₀=55746.0MJD, process step-length τ _s=1s, observation interval τ _obs=2s.In observation time section, collect whole photon time of arrival, be designated as τ _k(k=1,2,3 ...).Total number of sample points N=1024.

Step 2: to first group of photon sequence τ time of arrival _k∈ [0s, 1s], folds according to quiescent conditions lower recurrence interval P0, and folding reference point is taken as 0.5s.First ranging pulse profile c (n)=0 is defined, (n=1,2 ..., 1024), if a kth photon meets time of arrival:

$\left\{\frac{{\mathrm{\τ}}_k-{\mathrm{\τ}}_s}{P_0}\right\}\×N\∈(n-1,n]$

Wherein { } means and rounds up, layout when then thinking that this photon to correspond to n-th of folding cycle time of arrival, corresponding ranging pulse profile function c (n)=c (n)+1, successively this group photon sequence data time of arrival is processed, until all photon data processing time of arrival is complete.

Step 3: according to nominal contour s (n) and measurement profile c (n), adopts both cross correlation algorithm calculating phase differential that is:

$\widetilde{\mathrm{\φ}}=\arg \left(\underset{k}{\max }\left(\mathrm{cov}(c\left(n\right),s_k\left(n\right))\right)\right)\×\frac{P_0}{1024}$

Wherein s _kn () expression carries out k position ring shift right to 1024 n dimensional vector ns, namely

s _k(n)＝s(mod(n+k,k))

Function cov () represents cross-correlation calculation, realizes by vector dot product mode.The time location that the maximum ring shift right amount of auto-correlation cross correlation value is corresponding and nominal contour and the phase differential measured between profile.

Step 4: utilize wave filter to phase difference estimation value carry out filtering.

According to the Laplace transformation of track loop, one group of feasible solution of each parameter in filtering computing formula can be pushed away to obtain, that is:

$\left\{\begin{array}{c}{a}_{33}=0.0166\\ a_{32}=0.0715\\ a_{11}+a_{31}=1.4119\end{array}\right.$

For subsequent integration calculate easy, definition intermediate variable x ₁, x ₂.Because current time is integration initial value, the calculating of the feedback modifiers information of frequency and phase place can be written as

$x_1=0.0166\widetilde{\mathrm{\φ}}\×{\mathrm{\τ}}_s$

$x_2=x_1+0.0715\widetilde{\mathrm{\φ}}\×{\mathrm{\τ}}_s$

$\mathrm{\Δf}=(1.4119\widetilde{\mathrm{\φ}}+x_2)/2$

Δφ＝Δf×τ _s

Δ φ is also that this algorithm of current slot exports simultaneously, the impulse phase estimated value after level and smooth.

Step 5: the recurrence interval P under utilizing the frequency information Δ f of feedback to revise quiescent conditions ₀, namely

$\hat{P}={(P_0^{-1}-\mathrm{\Δf})}^{-1}$

Utilize the recurrence interval P after upgrading to next group photon sequence τ time of arrival _k∈ [1s, 2s] folds, and its method is similar to step 3, and difference is that the recurrence interval P after with renewal substitutes the recurrence interval P under quiescent conditions ₀participate in calculating, and obtain this group photon sequence c time of arrival (n).

Step 6: utilize the phase information Δ φ of feedback to revise nominal contour s (n), that is:

$\hat{s}\left(n\right)=s\left(\mathrm{mod}(n+\frac{\mathrm{\Δ\φ}\×N}{P},N)\right)$

Step 7: according to the nominal contour after renewal obtain with step 5, by current photon sequence τ time of arrival _kmeasurement profile c (n) that ∈ [1s, 2s] generates, carries out cross-correlation calculation, poor to obtain this section of photon sequence phase time of arrival.Its method is similar to step 3.

Step 8: repeat step 4 to seven, because photon observation information is finished using, now only need the operation carrying out step 4:

$x_1\left(t\right)=x_1(t-1)+0.0166\widetilde{\mathrm{\φ}}\×{\mathrm{\τ}}_s$

$x_2\left(t\right)=x_2(t-1)+x_1\left(t\right)+0.0715\widetilde{\mathrm{\φ}}\×{\mathrm{\τ}}_s$

$\mathrm{\Δf}=(1.4119\widetilde{\mathrm{\φ}}+x_2\left(t\right))/2$

Δφ＝Δφ(t-1)+Δf×τ _s

Δ φ is also that this algorithm of current slot exports simultaneously, the impulse phase estimated value after level and smooth.

By observation time τ _obsextend to 7200s, the pulsar x-ray photon time of arrival that simulation geo-synchronous orbit satellite is observed, carry out this method simultaneously and do not adopt this method, by means of only cross-correlation phase estimation method, obtain two kinds of method phase estimation error curves as shown in Fig. 3 (a) He Fig. 3 (b).Cross-correlation method is only adopted to estimate that the phase error root mean square obtained is 5.8631e-004 (rad); The phase error root mean square adopting this method to obtain is 4.0603e-004 (rad), precision improvement 30%.Visible the inventive method effectively can improve pulsar navigation impulse phase measuring accuracy.

The content be not described in detail in instructions of the present invention belongs to the known technology of those skilled in the art.

Claims (2)

1., based on the X-ray pulsar phase estimation method of loop tracks, it is characterized in that comprising the steps:

(1) determine to observe pulsar, under record quiescent conditions, observe the recurrence interval P of pulsar ₀with full sized pules profile s (n); Wherein full sized pules profile is that the mathematical expectation of pulsar photon arriving amt in one-period is at the n-th sampled point the sampled value at place, N is total number of sample points, n=1,2 ..., N;

(2) at observation time section τ ₀to τ ₀+ τ _obsin, record the time of arrival of all observation pulsar photons, wherein τ ₀for observation time starting point, τ _obsfor observing total duration;

(3) for first τ _sobservation pulsar photon sequence τ time of arrival in duration _k∈ [τ ₀, τ ₀+ τ _s], according to quiescent conditions lower recurrence interval P ₀fold, generate ranging pulse profile c (n),

c(n)＝αs(n+φ×N)+β+ω(n)

Wherein α is scale factor, and β is shift factor, and ω (n) is measurement noises, and φ is the phase differential between ranging pulse profile and full sized pules profile;

(4) according to full sized pules profile s (n) and ranging pulse profile c (n), the estimated value of the phase difference of the two is obtained

(5) select single order and three rank filter banks, aim at the phase difference estimation value of pulse profile s (n) and ranging pulse profile c (n) carry out filtering, obtain frequency feedback update information Δ f and phase feedback update information Δ φ; Wherein filter bank is expressed as:

$\mathrm{\Δf}=((a_{31}\widetilde{\mathrm{\φ}}+\∫(a_{32}+{\∫a}_{33}\widetilde{\mathrm{\φ}}))+a_{11}\widetilde{\mathrm{\φ}})/2$

Δφ＝∫Δf

Wherein, Δ φ, Δ f are respectively the feedback modifiers information of phase place and frequency, a ₃₁, a ₃₂, a ₃₃represent each integral element amplification coefficient of three rank wave filters respectively, a ₁₁represent firstorder filter amplification coefficient;

(6) the recurrence interval P under utilizing frequency feedback update information Δ f to revise quiescent conditions ₀, obtain the revised recurrence interval and utilize to next τ _sin duration, observation pulsar photon sequence time of arrival folds, and obtains next τ _sin duration, observation pulsar photon sequence time of arrival, obtains new ranging pulse profile thus

(7) utilize phase feedback update information Δ φ to revise full sized pules profile s (n), obtain revised full sized pules profile $\tilde{s}\left(n\right)=s\left(\mathrm{mod}(n+\frac{\mathrm{\Δ\φ}\×N}{\hat{P}},N)\right);$

(8) basis with obtain new phase difference estimation value

(9) for observation time section τ ₀to τ ₀+ τ _obsinterior each process step-length τ _s, repeat step (3) ~ step (8), until all observation pulsar photon series processing time of arrival is complete, obtain final phase difference estimation value.

2. the X-ray pulsar phase estimation method based on loop tracks according to claim 1, is characterized in that: the integral operation in described filter bank adopts square wave digital integrator to realize.