CN102506876B - Self-contained navigation method for measurement of earth ultraviolet sensor - Google Patents

Abstract

The invention relates to a self-contained navigation method for measurement of an earth ultraviolet sensor, which comprises the following steps of: initializing: utilizing satellite injection information to supply a state initial value X~1 and a covariance initial value P~1 outlier of a satellite at a starting moment, accumulating geocentric vector information and earth apparent radius information obtained by the ultraviolet sensor for a specified period of time, and then utilizing an outlier eliminating algorithm to eliminate the outlier in the information; batch processing: performing a least square method on the geocentric vector information and the earth apparent radius information after the outlier is eliminated, and calculating a state estimation X`k and a covariance estimation Pk of the satellite at an ending moment of observation; unscented transformation: utilizing the unscented transformation to realize the time updating of the state and the covariance, obtaining a state forecast X~(k+1) and covariance forecast P~(k+1) of the satellite at the next ending moment of observation; and recursion: continuously accumulating the geocentric vector information and the earth apparent radius information for a period of time, and repeating the steps and performing continuous recursion and calculation, thereby continuously obtaining the estimation of the state of the satellite. According to the method, the precision and the stability of a navigation system are increased.

Description

The autonomous navigation method that a kind of earth ultraviolet sensors is measured

Technical field

The invention belongs to Spacecraft Autonomous Navigation field, relate to the autonomous navigation method that a kind of earth ultraviolet sensors is measured.

Background technology

Earth ultraviolet sensors can detect the image at whole earth edge, utilize a large amount of marginal points information to fit to the geometric center of picture, thereby obtain the measurement data of the earth's core vector earth apparent radius, simultaneously ultraviolet sensors can also responsive stellar radiation, and the attitude information in satellite relative inertness space is provided.Utilize the abundant metrical information of earth ultraviolet sensors can realize the independent navigation of earth satellite completely, thereby improve independence and the autonomous viability of satellite, be conducive to reduce the degree of dependence of satellite for ground system simultaneously.

The metrical information of earth ultraviolet sensors inevitably contains noise, can not directly calculate the position and speed of satellite according to geometric relationship, but will estimate by filtering algorithm the position and speed of satellite.Conventional representational filtering algorithm has Kalman filtering and least-squares estimation etc. at present.EKF is a kind of Recursive Filtering algorithm, if filtering initial error very greatly, needs the long period to restrain, and is subject to wild value simultaneously and disturbs.Therefore be not suitable for using Initial orbit determination.LSR filter algorithm to the susceptibility of initial error a little less than, can reject wild value processes simultaneously, but need to realize by state-transition matrix the time renewal of state and covariance, this process exists very large linearization error, has affected the further raising of navigation accuracy.

Summary of the invention

Technology of the present invention is dealt with problems and is: overcome the deficiencies in the prior art, the autonomous navigation method that provides a kind of ultraviolet sensors to measure, has not only ensured the numerical stability of filtering but also improved navigation accuracy.

Technical solution of the present invention is: the autonomous navigation method that a kind of ultraviolet sensors is measured, and step is as follows:

(1) initialization utilizes the satellite information of entering the orbit to provide satellite at initial time state initial value

with covariance initial value

(2) elimination of burst noise

The time of the earth's core Vector Message that ultraviolet sensors is obtained and one section of appointment of earth apparent radius information accumulative total, then utilize unruly-value rejecting algorithm that the open country value in above-mentioned information is removed;

(3) batch processing

The earth's core Vector Message after the elimination of burst noise that step (2) is obtained and earth apparent radius information exchange are crossed least square method, the state estimation of the calculating observation satellite finish time

estimate P with covariance _k, specific as follows:

${\hat{x}}_k={\tilde{x}}_k+P_k{\tilde{H}}_k^{T}W(z_k-{\tilde{z}}_k)---\left(11\right)$

$P_k={({\tilde{P}}_k+{\tilde{H}}_k^{T}W{\tilde{H}}_k)}^{-1}---\left(12\right)$

Wherein, subscript k represents cycle index;

it is status predication; W is Weighting Matrices, gets

σ ₀it is the standard deviation of observation noise;

it is covariance prediction;

for total observing matrix, ${\tilde{H}}_k=\left[\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="H"\; filenames="HSA00000633147200022.png"\; state="\{\{state\}\}">H</figure-callout>}}_1\mathrm{\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_1)\\ H_2\mathrm{\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_2)\\ .\\ .\\ .\\ H_m\mathrm{\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_m)\end{array}\right],$

h _ibe observation equation, m represents the earth's core Vector Message and earth apparent radius acquisition of information order, Φ (t ₀, t _i) (i=1,2 ..., m) be state-transition matrix, by the numerical integration differential equation

obtain, Φ (t ₀, t ₀)=I, t ₀be initial time, I is unit matrix; z _k=[z _k1z _k2z _km] ^tobservation vector, i.e. the earth's core Vector Message and earth apparent radius information;

it is the predicted value of observation vector;

(4) without mark conversion (unscented transformation, UT)

Utilize the time renewal that realizes state and covariance without mark conversion, obtain the status predication of next section of observation satellite finish time

predict with covariance

detailed process is as follows:

A. calculate weighting coefficient

with

$w_0^m=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}---\left(13\right)$

$w_0^c=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}+1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;}---\left(14\right)$

$w_i^m=w_i^c=\frac{1}{2(n+\mathrm{\&lambda;})},i=1,...,2n---\left(15\right)$

Wherein, λ=α ²(n+ κ)-n is a scalar; N is state variable dimension; The distribution that constant alpha control sigma is ordered, its span is 10 ^-4≤ α≤1, regulates α can reduce the higher order term impact of nonlinear equation, and the value of κ should be guaranteed matrix (n+ λ) P _kfor positive semidefinite matrix, β is a parameter information-related with state prior distribution;

B. produce sigma point

S _k＝{chol(P _K)} ^T (16)

$X_k=\left[\begin{array}{ccc}{\hat{x}}_k& {\hat{x}}_k+\mathrm{\&gamma;}{S}_k& {\hat{x}}_k-{\mathrm{\&gamma;S}}_k\end{array}\right]---\left(17\right)$

Wherein, chol () represents that Qiao Lisiji (cholesky) decomposes, and positive definite matrix X is had to chol (X)=R, meets R ^tr=X, and R is upper triangular matrix;

C. the time upgrades

$X_{k+1|k}=X_k+{\&Integral;}_{t_k}^{t_{k+1}}f\left(x\right)\mathrm{dt}---\left(18\right)$

${\tilde{x}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_i^m{X}_{i,k+1|k}---\left(19\right)$

${\tilde{P}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_i^c(X_{i,k+1|k}-{\tilde{x}}_{k+1}){(X_{i,k+1|k}-{\tilde{x}}_{k+1})}^T---\left(20\right)$

Wherein, f (x) is system equation, is to set up according to the principle of dynamics of satellite; X _{i, k+1|k}vectorial X _k+1|kin since 0 counting i component, i=0,1 ..., 2n;

(5) recursion

Continue accumulative total a period of time the earth's core Vector Message and earth apparent radius information, and repeat above-mentioned steps (2) to step (4), so constantly recursive operation goes down to obtain continuously the estimation to satellitosis.

κ=3-n in described step (4) a, β=2.

Principle of the present invention: earth ultraviolet sensors can detect the image at whole earth edge, utilizes a large amount of marginal points information to fit to the geometric center of picture, thereby obtains the measurement data of the earth's core vector earth apparent radius.When initial orbit determination, earth ultraviolet sensors Continuous Observation a period of time, accumulate the earth's core vector apparent radius corresponding to multiple moment and measured.Then, this batch of measurement rejected to wild value to be processed, reject the measurement of estimated bias, again by all the other measure disposable send into Navigation Filter carry out a Post Orbit determine, estimate satellite position and speed in J2000.0 geocentric equatorial polar coordinate, then to upgrade the state that provides next observation segmental arc be its covariance prediction the elapsed time.After the measurement and orbit determination of multiple segmental arcs, position and velocity estimation on star can converge to degree of precision.For avoiding linearization error effect, introduce without mark and convert to realize time renewal.Thought without mark conversion is: suppose state vector Gaussian distributed, utilize sigma point that cluster is chosen meticulously after nonlinear transformation, to ask for single order and the second moment of state, utilizing can average and the covariance of accurate Calculation after nonlinear system equation transmits without mark conversion, has finally improved navigation accuracy.

The present invention's advantage is compared with prior art: the present invention is different from legacy card Kalman Filtering recursion processing metrical information, by the measurement accumulative total in a period of time, first reject wherein wild value, then measure common certain moment satellite orbit of estimation with all the other, utilize the split-second precision renewal of carrying out state and covariance without mark conversion simultaneously, this method of the present invention has been got rid of the interference of the wild value of measurement to orbit determination, the numerical stability of filtering and the speed of convergence are improved, avoid the linearization error interference of time renewal process, finally realized the raising of navigation accuracy.

Brief description of the drawings

Fig. 1 is the realization flow figure of the inventive method;

Fig. 2 is the autonomous navigation method schematic diagram that ultraviolet sensors of the present invention is measured;

Fig. 3 is site error and the velocity error of traditional autonomous navigation method;

Fig. 4 is site error and the velocity error of the autonomous navigation method that provides of the present invention.

Embodiment

Earth ultraviolet sensors can detect the image at whole earth edge, utilizes a large amount of marginal points information to fit to the geometric center of picture, thereby obtains the measurement data of the earth's core vector earth apparent radius.When initial orbit determination, earth ultraviolet sensors Continuous Observation a period of time, accumulate the earth's core vector apparent radius corresponding to multiple moment and measured, first reject wherein wild value, then do least-squares estimation with all the other measurements, estimate satellite position and speed in J2000.0 geocentric equatorial polar coordinate.Known system equation and observation equation are as follows

$\stackrel{\&CenterDot;}{x}=f\left(x\right)+w---\left(21\right)$

z＝h(x)+v (22)

Wherein, x is satellitosis; Z is observed quantity; W is system noise; V measures noise.

It is as follows that the present invention proposes the concrete calculation procedure of autonomous navigation method that earth ultraviolet sensors measures:

(1) initialization

Utilize the satellite information etc. of entering the orbit to be given in initial time state initial value

with covariance initial value

(2) elimination of burst noise

The time of one section of appointment of the earth's core vector earth apparent radius information accumulative total that earth ultraviolet sensors is obtained, then utilize unruly-value rejecting algorithm that the open country value in above-mentioned information is removed;

(3) batch processing

The earth's core Vector Message after the elimination of burst noise that step (2) is obtained and earth apparent radius information exchange are crossed least square method, the state estimation of the calculating observation satellite finish time

estimate P with covariance _k, detailed process is as follows:

A) using observation equation to the local derviation of state x the observing matrix H as single measurement _i, as follows

$H_i=\frac{{\&PartialD;h}_i}{\&PartialD;x}{|}_{x={\tilde{x}}_k},(i=\mathrm{1,2},...,m)---\left(23\right)$

In formula, h _iit is the corresponding observation equation of measuring for the i time; it is status predication; M represents the earth's core Vector Message and earth apparent radius acquisition of information order.

B) below numerical integration, the differential equation is tried to achieve state-transition matrix Φ (t ₀, t _i)

$\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_i)=\frac{\&PartialD;\stackrel{\&CenterDot;}{x}\left(t_i\right)}{\&PartialD;x\left(t_i\right)}\frac{\&PartialD;x\left(t_i\right)}{\&PartialD;x\left(t_0\right)}=\mathrm{A\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_i),(i=\mathrm{1,2},...,m)---\left(24\right)$

In formula,

Φ (t ₀, t ₀)=I, t ₀be initial time, I is unit matrix.

C) calculate total observing matrix

$\tilde{H}=\left[\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="H"\; filenames="HSA00000633147200022.png"\; state="\{\{state\}\}">H</figure-callout>}}_1\mathrm{\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_1)\\ H_2\mathrm{\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_2)\\ .\\ .\\ .\\ H_w\mathrm{\&Phi;}({\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HSA00000633147200022.png,HSA00000633147200023.png"\; state="\{\{state\}\}">t</figure-callout>}}_0,t_w)\end{array}\right]---\left(25\right)$

D) estimated state and covariance

$P_k={({\tilde{P}}_k+{\tilde{H}}_k^{T}W{\tilde{H}}_k)}^{-1}---\left(26\right)$

${\hat{x}}_k={\hat{x}}_k+P_k{\tilde{H}}_k^{T}W(z_k-{\tilde{z}}_k)---\left(27\right)$

Wherein, subscript k represents cycle index; W is Weighting Matrices, gets

σ ₀it is the standard deviation of observation noise; P _kthat covariance is estimated;

it is covariance prediction;

for observing matrix;

it is state estimation; it is status predication; z _k=[z _k1z _k2z _km] ^tit is observation vector;

it is the predicted value of observation vector.

(4) without mark conversion (unscented transformation, UT), utilize the time renewal that realizes state and covariance without mark conversion, obtain next observation window status predication finish time predict with covariance

detailed process is as follows:

A) calculate weighting coefficient

with

$w_0^m=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}--\left(28\right)$

$w_0^c=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}+1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;}---\left(29\right)$

$w_i^m=w_i^c=\frac{1}{2(n+\mathrm{\&lambda;})},i=1,...,2n---\left(30\right)$

Wherein, λ=α ²(n+ κ)-n is a scalar; N is state variable dimension.The distribution that constant alpha control sigma is ordered, its span is 10 ^-4≤ α≤1, regulates α can reduce the higher order term impact of nonlinear equation.The value of κ should be guaranteed matrix (n+ λ) P _kfor positive semidefinite matrix, get κ=3-n.β is a parameter information-related with state prior distribution, gets β=2.

B) produce sigma point

S _k＝{chol(P _k)} ^T (31)

$X_k=\left[\begin{array}{ccc}{\hat{x}}_k& {\hat{x}}_k+\mathrm{\&gamma;}{S}_k& {\hat{x}}_k-{\mathrm{\&gamma;S}}_k\end{array}\right]---\left(32\right)$

Wherein, chol () represents that Qiao Lisiji (cholesky) decomposes, and positive definite matrix X is had to chol (X)=R, meets R ^tr=X, and R is upper triangular matrix;

C) time upgrades

$X_{k+1|k}=X_k+{\&Integral;}_{t_k}^{t_{k+1}}f\left(x\right)\mathrm{dt}---\left(33\right)$

${\tilde{x}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_i^m{X}_{i,k+1|k}---\left(34\right)$

${\tilde{P}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_i^c(X_{i,k+1|k}-{\tilde{x}}_{k+1}){(X_{i,k+1|k}-{\tilde{x}}_{k+1})}^T---\left(35\right)$

Wherein, f (x) is that system equation is the expression formula of right-hand member in formula (21), is to set up according to the principle of dynamics of satellite; X _{i, k+1|k}vectorial X _k+1|kin since 0 counting i component, i=0,1 ..., 2n.

(5) recursion

Continue accumulative total a period of time the earth's core Vector Message and earth apparent radius information, and repeat above-mentioned steps (2) to step (4), so constantly recursive operation goes down to obtain continuously the estimation to satellitosis.

The unruly-value rejecting algorithm of mentioning in step (2) belongs to mature technology in this area, has the several different methods such as 3 σ thresholding diagnostic methods to adopt, and non-the claims in the present invention content, therefore be not described in detail.

Fig. 2 is the schematic diagram of the autonomous navigation method that proposes of the present invention, wherein

with

(j=1,2 ...) be status predication and the covariance prediction of satellite while carrying out the j time navigation; z _ji(j=1,2, I=1,2 ..., the earth's core Vector Message and the earth apparent radius information of the rejecting wild value that m) obtain for earth ultraviolet sensors in j section observation time; BLS is illustrated in the least square method providing in step (3);

and P _jbe respectively the state estimation and the covariance estimation that after the j time navigation calculated, obtain satellite; UT is illustrated in converting without mark of providing in step (4).

Simulation example: true track utilizes STK software to produce, and has adopted altitude of the apogee 5000km, the elliptical orbit of perigee altitude 500km, kinetic model is got 50 × 50 rank gravitational fields, has considered the perturbations such as the sun, lunar gravitation, atmospherical drag and sun optical pressure.0.1 ° of earth ultraviolet sensors pointing accuracy, 0.1 ° of apparent radius precision; The attitude determination accuracy of star sensor combination gyro is 0.03 °, and navigation sensor is exported with the sampling interval of 30 seconds.Preliminary orbit error: site error 100km, velocity error 10m/s.

Simulation result as shown in Figure 3 and Figure 4.Fig. 3 is the navigation results obtaining that adopts traditional autonomous navigation method, and its positional precision is 20km, and velocity accuracy is 8m/s.Fig. 4 is the navigation results that adopts the autonomous navigation method that provides of the present invention to obtain, and its positional precision is 6km, and velocity accuracy is 2m/s, compares classic method precision and is more than doubled.

In instructions of the present invention, other content not being described in detail belongs to those skilled in the art's known technology.

Claims (2)

1. the autonomous navigation method that earth ultraviolet sensors is measured, is characterized in that step is as follows:

(1) initialization

Utilize the satellite information of entering the orbit to provide satellite at initial time state initial value

with covariance initial value

(2) elimination of burst noise

The time of the earth's core Vector Message that earth ultraviolet sensors is obtained and one section of appointment of earth apparent radius information accumulative total, then utilize unruly-value rejecting algorithm that the open country value in above-mentioned information is removed;

(3) batch processing

The earth's core Vector Message after the elimination of burst noise that step (2) is obtained and earth apparent radius information exchange are crossed least square method, the state estimation of the calculating observation satellite finish time estimate P with covariance _k, specific as follows:

${\hat{x}}_k={\tilde{x}}_k+P_k{\tilde{H}}_k^{T}W(z_k-{\tilde{z}}_k)---\left(1\right)$

$P_k={({\tilde{P}}_k+{\tilde{H}}_k^{T}W{\tilde{H}}_k)}^{-1}---\left(2\right)$

Wherein, subscript k represents cycle index;

it is status predication; W is Weighting Matrices, gets σ ₀it is the standard deviation of observation noise;

it is covariance prediction;

for total observing matrix, ${\tilde{H}}_k=\left[\begin{array}{c}{H}_1\mathrm{\&Phi;}(t_0,t_1)\\ H_2\mathrm{\&Phi;}(t_0,t_2)\\ .\\ .\\ .\\ H_m\mathrm{\&Phi;}(t_0,t_m)\end{array}\right],$

h _ibe observation equation, m represents the earth's core Vector Message and earth apparent radius acquisition of information order, Φ (t ₀, t _i) (i=1,2 ..., m) be state-transition matrix, by the numerical integration differential equation

obtain, Φ (t ₀, t ₀)=I, t ₀be initial time, I is unit matrix; z _k=[z _k1z _k2z _km] ^tobservation vector, i.e. the earth's core Vector Message and earth apparent radius information;

it is the predicted value of observation vector;

(4) convert without mark

Utilize the time renewal that realizes state and covariance without mark conversion, obtain the status predication of next section of observation satellite finish time

predict with covariance

detailed process is as follows:

A. calculate weighting coefficient

with

$w_0^m=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}---\left(3\right)$

$w_0^c=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}+1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;}---\left(4\right)$

$w_i^m=w_i^c=\frac{1}{2(n+\mathrm{\&lambda;})},i=1,...,2n---\left(5\right)$

Wherein, λ=α ²(n+ κ)-n is a scalar; N is state variable dimension; The distribution that constant alpha control sigma is ordered, its span is 10 ^-4≤ α≤1, regulates α can reduce the higher order term impact of nonlinear equation, and the value of κ should be guaranteed matrix (n+ λ) P _kfor positive semidefinite matrix, β is a parameter information-related with state prior distribution;

B. produce sigma point

S _k＝{chol(P _k)} ^T (6)

$X_k=\left[\begin{array}{ccc}{\hat{x}}_k& {\hat{x}}_k+\mathrm{\&gamma;}{S}_k& {\hat{x}}_k-{\mathrm{\&gamma;S}}_k\end{array}\right]---\left(7\right)$

Wherein, chol () represents Cholesky factorization, and positive definite matrix X is had to chol (X)=R, meets R ^tr=X, and R is upper triangular matrix;

C. the time upgrades

$X_{k+1|k}=X_k+{\&Integral;}_{t_k}^{t_{k+1}}f\left(x\right)\mathrm{dt}---\left(8\right)$

${\tilde{x}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_i^m{X}_{i,k+1|k}---\left(9\right)$

${\tilde{P}}_{k+1}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_i^c(X_{i,k+1|k}-{\tilde{x}}_{k+1}){(X_{i,k+1|k}-{\tilde{x}}_{k+1})}^T---\left(10\right)$

Wherein, f (x) is system equation, is to set up according to the principle of dynamics of satellite; X _{i, k+1|k}vectorial X _k+1|kin since 0 counting i component, i=0,1 ..., 2n;

(5) recursion

Continue accumulative total a period of time the earth's core Vector Message and earth apparent radius information, and repeat above-mentioned steps (2) to step (4), so constantly recursive operation goes down to obtain continuously the estimation to satellitosis.

2. the autonomous navigation method that earth ultraviolet sensors according to claim 1 is measured, is characterized in that: the κ=3-n in described step (4) a, β=2.

Images