CN1939807A - Star sensor online aligning method based on weng model - Google Patents

Abstract

An on orbit calibration method based on Weng model for the star-sensitive device includes such steps as creating a posture transform array of star-sensitive device, creating the Weng imaging model, acquiring data, calculating parameters, and calibrating.

Description

Based on the star sensor of weng model at the rail alignment method

Technical field

The invention belongs to the aerospace measurement technology, relate to the improvement of star sensor at the rail alignment method.

Background technology

Star sensor is a kind of star observation that utilizes, and the aerospace measurement instrument of high-precision attitude information is provided for space vehicle.Star sensor the rail alignment method be star sensor research and use in a gordian technique.Usually before star sensor is launched, its inner parameter such as principal point, meetings such as focal length and distortion factor are demarcated on ground.Calibration method has night sky demarcation and the demarcation of starlight laboratory etc.But after the aircraft emission, because the impact during emission and the situation of change of working environment, all can be different from surface state as gravity, atmosphere and temperature etc., make the inner parameter of star sensor deviation occur, wearing out after the long-term use in addition also can influence operating accuracy.

Existingly mainly divide two classes at the rail alignment method, a class is the calibration that depends on attitude information; One class is the calibration according to angular distance invariance principle (not relying on attitude information) in the star.The calibration that depends on attitude information is to provide a known attitude by the outside, can calculate the ideal position of current visual field fixed star on the star sensor target surface according to this attitude, as the star sensor inner parameter, as principal point, when focal length and distortion factor change, the fixed star measuring position will with ideal position generation deviation.According to a series of deviates that obtain, can recalibrate inner parameter.Obvious defects of this method is owing to need the outside that an attitude is provided, thereby if outside attitude has error, then this error can be introduced in the calibration process, influences calibration accuracy.

Do not rely on attitude information the rail alignment method be according to star in angular distance orthogonal transformation invariance principle, detect the deviation between the angular distance observed reading and actual value in rail flight course culminant star, utilize basic function to come match angular distance deviation.The shortcoming that this method exists is that because the fit procedure employing is star interior angle relative value, 0 rank item of error is by cancellation; Simultaneously, basic function is selected the fitting precision influence bigger, and basic function is selected improper, and the unsettled situation of numerical value appears in the possibility of result.

U.S. JPL laboratory has and adopts radial base neural net to carry out star sensor to calibrate at rail in addition.This neural network groundwork also is a kind of method that relies on attitude, just utilizes the concurrent operation ability of neural network to handle the deviation of estimating between the star vectorial sum measurement star vector.This method shortcoming is that precision is not high, and the calculation resources consumption rate is bigger.

Summary of the invention

The objective of the invention is: at present star sensor in the problem that rail calibration exists, propose a kind of based on the weng model at the rail alignment method.This method is utilized the weng model in the machine vision, divides linear parameter and two steps of nonlinear parameter to estimate interior square element.This method distortion model is complete, can handle complicated lens of star sensor distortion situation.Complicated parameter optimization is finished and avoided to data handing in two steps, data convergence fast and stable.

Technical scheme of the present invention is:

At the rail alignment method, it is characterized in that the step of calibration is as follows based on the star sensor of weng model:

1, sets up star sensor attitude transition matrix;

1.1, set up system of celestial coordinates and star sensor system of axes; If O '-XnYnZn is a system of celestial coordinates, O-XYZ is the star sensor system of axes, and the attitude angle of star sensor is by right ascension α ₀, declination β ₀, and roll angle φ ₀Form α ₀Be the projection of Z axle on the XnYn face and the angle of Xn axle; β ₀Be Z axle and its angle between the projection on the XnYn face; φ ₀Angle for meridian plane and XY hand-deliver line and Y-axis;

1.2, set up star sensor attitude transition matrix; The switching process that is tied to the star sensor system of axes from celestial coordinates is: the first step, rotate around the Zn axle

The angle makes the Xn axle vertical with meridian plane; In second step, rotate around new Xn axle The angle makes the Zn axle consistent with the Z axle; In the 3rd step, rotate φ around new Zn axle ₀The angle makes system of celestial coordinates and star sensor system of axes overlap.Then the direction cosine matrix R of this rotation process is the attitude transition matrix, is expressed as:

In the formula, s represents sinusoidal function, and c represents cosine function;

R can be simplified shown as: $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right],$ r ₁～r ₉9 elements of expression R matrix;

2, set up star sensor weng imaging model;

If the starlight direction vector is: $w=\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[2\right]$

In the formula, w1, w2 and w3 represent the x of starlight vector under the star sensor system of axes, y, z coordinate figure; α, β are right ascension, the declination coordinate of fixed star;

2.1, use star sensor that starlight is carried out the radial distortion imaging; O-XYZ is the star sensor system of axes, and o-xy is 2 dimension target surface system of axess, and П is a target surface, and distance is a focal length between the Oo; The desirable perspective imaging point of starlight vector w is P ', because radial distortion takes place, actual imaging point is P;

2.2, set up the linear imaging model; Owing to have uncertain factor of proportionality between x and the y direction, use f _uAnd f _vThe focal length of representing x and y direction respectively, set up the linear imaging model and be:

$\frac{x^{\′}}{f_u}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}---\left[3\right]$

$\frac{y^{\′}}{f_v}=\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

In the formula, x ', y ' defines linear parameter vector d for the asterism position of perspective projection ₁=(f _uf _vα ₀β ₀_ ₀);

2.3, consider distortion effects, set up weng nonlinear distortion model formation and be:

$\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}=\hat{u}(g_1+g_3){\hat{u}}^2+g_4\hat{u}\hat{v}+g_1{\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)---\left[4\right]$

$\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}=\hat{v}+g_2{\hat{u}}^2+g_3\hat{u}\hat{v}+(g_2+g_4){\hat{v}}^2{+k}_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

In the formula, $\hat{u}=\frac{x}{f_u},$ $\hat{v}=\frac{y}{f_v},$ X, y are asterism actual measurement coordinate, and k1 is a coefficient of radial distortion, and g1～g4 is the aggregative formula coefficient of decentering distortion and thin prism distortion, so nonlinear distortion parameter has 5, and defined parameters vector d ₂=(k ₁g ₁g ₂g ₃g ₄);

3, data acquisition; According to the frame star chart that star sensor collects, whole asterism data acquisitions that definition collects are Ω ₀, the asterism number is n; With the principal point is the center, and diameter is that picture size 1/2 is done circle, and the data acquisition in the note circle is Ω ₁, the asterism number is n ₁

4, calibrate by calculation of parameter;

Calibration process be divided into two the step finish, according to centre data Ω ₁Be subjected to the less principle of distortion effects, the first step is utilized Ω ₁Construct linear computation model, calculate the initial value of parameter vector d1; The parameter of second step to d1 is optimized processing; In the 3rd step, according to the weng model, structure distortion factor linear equation is at global data Ω ₀Following calculating parameter vector d2; At last, whether judgment data result meets the demands, if do not satisfy then return the second stepping row iteration and calculate, if satisfied the output parameter vector finally separate d1 and d2; Concrete computation process is as follows:

4.1, carry out linear dimensions and estimate; Have according to formula [3]:

x′w ₁r ₇+x′w ₂r ₈+x′w ₃r ₉＝w ₁r ₁f _u+w ₂r ₃f _u+w ₃r ₃f _u [5]

y′w ₁r ₇+y′w ₂r ₈+y′w ₃r ₉＝w ₁r ₄f _v+w ₂r ₅f _v+w ₃r ₆f _v

Arrangement has:

$w_1=\frac{r_1{f}_u}{r_9}+w_2\frac{r_2{f}_u}{r_9}+w_3\frac{r_3{f}_u}{r_9}-{\mathrm{xw}}_1\frac{r_7}{r_9}-{\mathrm{xw}}_2\frac{r_8}{r_9}={\mathrm{xw}}_3---\left[6\right]$

$w_1=\frac{r_4{f}_v}{r_9}+w_2\frac{r_5{f}_v}{r_9}+w_3\frac{r_6{f}_v}{r_9}-{\mathrm{yw}}_1\frac{r_8}{r_9}-{\mathrm{yw}}_2\frac{r_8}{r_9}={\mathrm{yw}}_3$

Suppose $r_1^{\&prime;}=\frac{r_1{f}_u}{r_9},r_2^{\&prime;}=\frac{r_2{f}_u}{r_9},r_3^{\&prime;}=\frac{r_3{f}_u}{r_9},{\mathrm{<figure-callout\; id="4"\; label="r"\; filenames="A20061014048100172.png"\; state="\{\{state\}\}">r</figure-callout>}}_4^{\&prime;}=\frac{r_4{f}_v}{r_9},r_5^{\&prime;}=\frac{r_5{f}_u}{r_9},r_6^{\&prime;}=\frac{r_6{f}_v}{r_9},r_7^{\&prime;}=\frac{r_7}{r_9},r_8^{\&prime;}=\frac{r_8}{r_9},$ Totally 8 unknown numberes, defined parameters vector m=[r ₁＇ r ₂＇ r ₃＇ r ₄＇ r ₅＇ r ₆＇ r ₇＇ r ₈＇] ^T, can obtain equation for each asterism data:

[w ₁ w ₂ w ₃ 0 0 0 -xw ₁ -xw ₂]m＝xw ₃ [7]

[0 0 0 w ₁ w ₂ w ₃ -yw ₁ -yw ₂]m＝yw ₃

For data acquisition Ω ₁Interior asterism can obtain 2n ₁Individual equation constitutes " Am=b " form; A is 2n ₁The matrix that individual equation levoform constitutes, b is 2n ₁The vector that the right formula of individual equation constitutes; Available method of least square makes Solve 8 unknown number r ₁＇～r ₈＇; Be the orthogonal matrix characteristics according to R then, have:

$r_9=\sqrt{1+r_7^{\&prime;}+r_8^{\&prime;}},r_7=r_7^{\&prime;}\&CenterDot;r_9,r_8=r_8^{\&prime;}\&CenterDot;r_9$

$f_u=r_9\&CenterDot;\sqrt{r_1^{\&prime;2}+r_2^{\&prime;2}+r_3^{\&prime;2}}$

r ₁＝r ₁′·r ₉/f _u，r ₂＝r ₂′·r ₉/f _u，r ₃＝r ₃′·r ₉/f _u

$f_v=r_9\&CenterDot;\sqrt{{\mathrm{<figure-callout\; id="4"\; label="r"\; filenames="A20061014048100172.png"\; state="\{\{state\}\}">r</figure-callout>}}_4^{\&prime;2}+r_5^{\&prime;2}+r_6^{\&prime;2}}$

r ₄＝r ₄′·r ₉/f _v，r ₅＝r ₅′·r ₉/f _v，r ₆＝r ₆′·r ₉/f _v

Earlier hypothesis r9 is for just, obtain the attitude estimated matrix, select a certain distance center fixed star far away in the visual field to judge that whether r9 is really for just, if starlight vector coordinate after the imaging of attitude estimated matrix is consistent with meeting of measurement coordinate then, r9 then is described for just, otherwise for negative;

4.2, optimize parameter vector d1;

If the actual value of linear dimensions vectorial sum nonlinear parameter vector is d ₁, d ₂, estimated valve is

Optimum estimated valve should satisfy:

$\underset{d}{\min }\left|\right|P-f({\hat{d}}_1,{\hat{d}}_2)\left|\right|$

In the formula, P is point set Ω ₁The asterism coordinate; With separating of obtaining of formula [8] is initial value, adopts nonlinear optimization method that d1 is optimized; If:

x′＝f _x(d ₁， d ₂) [9]

y′＝f _y(d ₁， d ₂)

In the formula, d ₂Expression parameter vector d2 fixes; Because f _xAnd f _yBe nonlinear function, therefore adopt non-linear linearization method to come estimated parameter vector d1; Suppose

Be estimated valve, x, y are observed reading, Δ d ₁Be vectorial estimated bias;

$\mathrm{\&Delta;x}=\stackrel{\&OverBar;}{x}-\hat{x}\&ap;\mathrm{A\&Delta;}{d}_1---\left[10\right]$

$\mathrm{\&Delta;y}=\stackrel{\&OverBar;}{y}-\hat{y}\&ap;\mathrm{B\&Delta;}{d}_1$

In the formula, A, B are sensitive matrix, ask every partial derivative to obtain;

So obtain parameter vector updating value d ₁=d ₁+ Δ d ₁

4.3, calculate distortion factor vector d2;

After obtaining linear dimensions, next step finds the solution nonlinear parameter d ₂According to formula [4], construct linear estimation model:

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& {\hat{u}}^2+{\hat{v}}^2& 0& {\hat{u}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=d_u---\left[11\right]$

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& 0& {\hat{u}}^2+{\hat{u}}^2& {\hat{v}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=d_v$

In the formula:

$\mathrm{du}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}-\hat{u}$

$\mathrm{dv}=\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}-\hat{v}$

Constitute " Bd ₂=c " form, B is the matrix that 2n equation levoform constitutes, c is the vector that the right formula of 2n equation constitutes; Utilize method of least square to make

Obtain the estimated valve of distortion parameter vector d2, the data here are whole point set Ω ₀

After obtaining whole parameter d 1 and d2 of model, judge by contrast ideal coordinates value and measurement coordinate figure; If error is still bigger, then get back to step 1.4.2 and further optimize d1, carry out iterative computation; If error meets the demands, then stop iteration, the output result.

Advantage of the present invention is:

1, reduced the related influence of coupling of ambient parameter and inner parameter;

2, distortion model is complete, and good comformability is arranged;

3, find the solution parameter based on linear enclosed, calculated amount is moderate.

Description of drawings

Fig. 1 is the star sensor attitude angle scheme drawing among the present invention.

Fig. 2 is the star sensor starlight imaging scheme drawing among the present invention.

Fig. 3 is the parameter calculation procedure scheme drawing among the present invention.

Fig. 4 is a simulation stellar field scheme drawing.

The specific embodiment

Below the present invention is described in further details.Of the present inventionly it is characterized in that at the rail alignment method step of calibration is as follows:

1, sets up star sensor attitude transition matrix;

1.1, set up system of celestial coordinates and star sensor system of axes; Referring to Fig. 1, establishing O '-XnYnZn is system of celestial coordinates, and O-XYZ is the star sensor system of axes, and the attitude angle of star sensor is by right ascension α ₀, declination β ₀, and roll angle φ ₀Form α ₀Be the projection of Z axle on the XnYn face and the angle of Xn axle; β ₀Be Z axle and its angle between the projection on the XnYn face; φ ₀Angle for meridian plane and XY hand-deliver line and Y-axis;

1.2, set up star sensor attitude transition matrix; The switching process that is tied to the star sensor system of axes from celestial coordinates is: the first step, rotate around the Zn axle

The angle makes the Xn axle vertical with meridian plane; In second step, rotate around new Xn axle The angle makes the Zn axle consistent with the Z axle; In the 3rd step, rotate φ around new Zn axle ₀The angle makes system of celestial coordinates and star sensor system of axes overlap.Then the direction cosine matrix R of this rotation process is the attitude transition matrix, is expressed as:

In the formula, s represents sinusoidal function, and c represents cosine function;

R can be simplified shown as: $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right],$ r ₁～r ₉9 elements of expression R matrix;

2, set up star sensor weng imaging model;

If the starlight direction vector is: $w=\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\&beta;}\cos \mathrm{\&alpha;}\\ \cos \mathrm{\&beta;}\sin \mathrm{\&alpha;}\\ \sin \mathrm{\&beta;}\end{array}\right)---\left[2\right]$

In the formula, w1, w2 and w3 represent the x of starlight vector under the star sensor system of axes, y, z coordinate figure; α, β are right ascension, the declination coordinate of fixed star;

2.1, use star sensor that starlight is carried out the radial distortion imaging; Be illustrated in figure 2 as star sensor radial distortion imaging scheme drawing.O-XYZ is the star sensor system of axes, and o-xy is 2 dimension target surface system of axess, and П is a target surface, and distance is a focal length between the Oo; The desirable perspective imaging point of starlight vector w is P ', because radial distortion takes place, actual imaging point is P;

2.2, set up the linear imaging model; Owing to have uncertain factor of proportionality between x and the y direction, use f _uAnd f _vThe focal length of representing x and y direction respectively, set up the linear imaging model and be:

$\frac{x^{\&prime;}}{f_u}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}---\left[3\right]$

$\frac{y^{\&prime;}}{f_v}=\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

In the formula, x ', y ' defines linear parameter vector d for the asterism position of perspective projection ₁=(f _uf _vα ₀β ₀_ ₀);

2.3, consider distortion effects, set up weng nonlinear distortion model formation and be:

$\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}=\hat{u}(g_1+g_3){\hat{u}}^2+g_4\hat{u}\hat{v}+g_1{\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)---\left[4\right]$

$\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}=\hat{v}+g_2{\hat{u}}^2+g_3\hat{u}\hat{v}+(g_2+g_4){\hat{v}}^2{+k}_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

In the formula, $\hat{u}=\frac{x}{f_u},\hat{v}=\frac{y}{f_v},$ X, y are asterism actual measurement coordinate, and k1 is a coefficient of radial distortion, and g1～g4 is the aggregative formula coefficient of decentering distortion and thin prism distortion, so nonlinear distortion parameter has 5, and defined parameters vector d ₂=(k ₁g ₁g ₂g ₃g ₄);

3, data acquisition; According to the frame star chart that star sensor collects, whole asterism data acquisitions that definition collects are Ω ₀, the asterism number is n; With the principal point is the center, and diameter is that picture size 1/2 is done circle, and the data acquisition in the note circle is Ω ₁, the asterism number is n ₁

4, calibrate by calculation of parameter;

Fig. 3 is the parameter calculation procedure scheme drawing.Calibration process be divided into two the step finish, according to centre data Ω ₁Be subjected to the less principle of distortion effects, the first step is utilized Ω ₁Construct linear computation model, calculate the initial value of parameter vector d1; The parameter of second step to d1 is optimized processing; In the 3rd step, according to the weng model, structure distortion factor linear equation is at global data Ω ₀Following calculating parameter vector d2; At last, whether judgment data result meets the demands, if do not satisfy then return the second stepping row iteration and calculate, if satisfied the output parameter vector finally separate d1 and d2; Concrete computation process is as follows:

4.1, carry out linear dimensions and estimate; Have according to formula [3]:

x′w ₁r ₇+x′w ₂r ₈+x′w ₃r ₉＝w ₁r ₁f _u+w ₂r ₂f _u+w ₃r ₃f _u [5]

y′w ₁r ₇+y′w ₂r ₈+y′w ₃r ₉＝w ₁r ₄f _v+w ₂r ₅f _v+w ₃r ₆f _v

Arrangement has:

$w_1=\frac{r_1{f}_u}{r_9}+w_2\frac{r_2{f}_u}{r_9}+w_3\frac{r_3{f}_u}{r_9}-{\mathrm{xw}}_1\frac{r_7}{r_9}-{\mathrm{xw}}_2\frac{r_8}{r_9}={\mathrm{xw}}_3---\left[6\right]$

$w_1=\frac{r_4{f}_v}{r_9}+w_2\frac{r_5{f}_v}{r_9}+w_3\frac{r_6{f}_v}{r_9}-{\mathrm{yw}}_1\frac{r_8}{r_9}-{\mathrm{yw}}_2\frac{r_8}{r_9}={\mathrm{yw}}_3$

Suppose $r_1^{\&prime;}=\frac{r_1{f}_u}{r_9},r_2^{\&prime;}=\frac{r_2{f}_u}{r_9},r_3^{\&prime;}=\frac{r_3{f}_u}{r_9},{\mathrm{<figure-callout\; id="4"\; label="r"\; filenames="A20061014048100172.png"\; state="\{\{state\}\}">r</figure-callout>}}_4^{\&prime;}=\frac{r_4{f}_v}{r_9},r_5^{\&prime;}=\frac{r_5{f}_u}{r_9},r_6^{\&prime;}=\frac{r_6{f}_v}{r_9},r_7^{\&prime;}=\frac{r_7}{r_9},r_8^{\&prime;}=\frac{r_8}{r_9},$ Totally 8 unknown numberes, defined parameters vector m=[r ₁＇ r ₂＇ r ₃＇ r ₄＇ r ₅＇ r ₆＇ r ₇＇ r ₈＇] ^T, can obtain equation for each asterism data:

[w ₁ w ₂ w ₃ 0 0 0 -xw ₁ -xw ₂]m＝xw ₃ [7]

[0 0 0 w ₁ w ₂ w ₃ -yw ₁ -yw ₂]m＝yw ₃

For data acquisition Ω ₁Interior asterism can obtain 2n ₁Individual equation constitutes " Am=b " form; A is 2n ₁The matrix that individual equation levoform constitutes, c is 2n ₁The vector that the right formula of individual equation constitutes; Available method of least square makes

Solve 8 unknown number r ₁＇～r ₈＇; Be the orthogonal matrix characteristics according to R then, have:

$r_9=\sqrt{1+r_7^{\&prime;}+r_8^{\&prime;}},r_7=r_7^{\&prime;}\&CenterDot;r_9,r_8=r_8^{\&prime;}\&CenterDot;r_9$

$f_u=r_9\&CenterDot;\sqrt{r_1^{\&prime;2}+r_2^{\&prime;2}+r_3^{\&prime;2}}$

r ₁＝r ₁′·r ₉/f _u，r ₂＝r ₂′·r ₉/f _u，r ₃＝r ₃′·r ₉/f _u

$f_v=r_9\&CenterDot;\sqrt{{\mathrm{<figure-callout\; id="4"\; label="r"\; filenames="A20061014048100172.png"\; state="\{\{state\}\}">r</figure-callout>}}_4^{\&prime;2}+r_5^{\&prime;2}+r_6^{\&prime;2}}$

r ₄＝r ₄′·r ₉/f _v，r ₅＝r ₅′·r ₉/f _v，r ₆＝r ₆′·r ₉/f _v

Earlier hypothesis r9 is for just, obtain the attitude estimated matrix, select a certain distance center fixed star far away in the visual field to judge that whether r9 is really for just, if starlight vector coordinate after the imaging of attitude estimated matrix is consistent with meeting of measurement coordinate then, r9 then is described for just, otherwise for negative;

4.2, optimize parameter vector d1;

If the actual value of linear dimensions vectorial sum nonlinear parameter vector is d ₁, d ₂, estimated valve is

Optimum estimated valve should satisfy:

$\underset{d}{\min }\left|\right|P-f({\hat{d}}_1,{\hat{d}}_2)\left|\right|$

In the formula, P is point set Ω ₁The asterism coordinate; With separating of obtaining of formula [8] is initial value, adopts nonlinear optimization method that d1 is optimized; If:

$x^{\&prime;}=f_x(d_1,{\stackrel{\&OverBar;}{d}}_2)---\left[9\right]$

$y^{\&prime;}=f_y(d_1,{\stackrel{\&OverBar;}{d}}_2)$

In the formula, d ₂Expression parameter vector d2 fixes; Because f _xAnd f _yBe nonlinear function, therefore adopt non-linear linearization method to come estimated parameter vector d1; Suppose

Be estimated valve, x, y are observed reading, Δ d ₁Be vectorial estimated bias;

$\mathrm{\&Delta;x}=\stackrel{\&OverBar;}{x}-\hat{x}\&ap;\mathrm{A\&Delta;}{d}_1---\left[10\right]$

$\mathrm{\&Delta;y}=\stackrel{\&OverBar;}{y}-\hat{y}\&ap;\mathrm{B\&Delta;}{d}_1$

In the formula, A, B are sensitive matrix, ask every partial derivative to obtain;

So obtain parameter vector updating value d ₁=d ₁+ Δ d ₁

4.3, calculate distortion factor vector d2;

After obtaining linear dimensions, next step finds the solution nonlinear parameter d ₂According to formula [4], construct linear estimation model:

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& {\hat{u}}^2+{\hat{v}}^2& 0& {\hat{u}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=d_u---\left[11\right]$

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& 0& {\hat{u}}^2+{\hat{u}}^2& {\hat{v}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=d_v$

In the formula:

$\mathrm{du}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}-\hat{u}$

$\mathrm{dv}=\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}-\hat{v}$

Constitute " Bd ₂=c " form, B is the matrix that 2n equation levoform constitutes, c is the vector that the right formula of 2n equation constitutes; Utilize method of least square to make

Obtain the estimated valve of distortion parameter vector d2, the data here are whole point set Ω ₀

After obtaining whole parameter d 1 and d2 of model, judge by contrast ideal coordinates value and measurement coordinate figure; If error is still bigger, then get back to step 1.4.2 and further optimize d1, carry out iterative computation; If error meets the demands, then stop iteration, the output result.

Emulation and interpretation of result

The star sensor basic specification of emulation is:

Visual field: 12 * 12;

Pel array: 1024 * 1024;

Pixel Dimensions: 0.015mm * 0.015mm;

Focal length: 73.0703mm

Based on the weng distortion model at the rail alignment method, the star number in the visual field is required many.Can adopt the two-dimensional random distributed points to simulate the distribution of fixed star in the visual field, the benefit of doing like this is easily to produce the stellar field that certain star number requires, and can study the stability of star number and calibration parameter convergence.Here suppose to have in the stellar field 100 stars.

Fig. 4 is simulation stellar field scheme drawing, and the circle diameter is 9.2mm among the figure, promptly 28.3% of the target surface area, in the circle 28 stars are arranged.Suppose that attitude angle is (0,0,0), the even random scattering of target surface asterism coordinate.

Simulation parameter is:

f _u＝76.7238mm，f _v＝73.0703mm。

k1＝-0.3，g1＝0.02，g2＝-0.009，g3＝-0.02，g4＝0.009。

The asterism noise level is the Gaussian white noise of 0.05 pixel variance.

The enclosed result of calculation of the first step is:

f u/mm	fv/mm	α 0	β0	φ0
					75.957	72.223	0.015	0.140	-0.00116

The iterative computation process in the second, three step is:

At first the computation process of formula parameter vector d1 sees the following form:

Iterations	f u/mm	fv/mm	α 0 /arcsec	β 0/arcsec	φ 0/arcsec
						1	77.373	74.495	-67.2277	939.3451	239.1954
2	76.917	73.206	-15.8128	201.4925	473.6030
						3	76.754	73.081	-2.4443	41.4535	103.5382
4	76.727	73.061	-0.7597	8.0169	20.8955
						5	76.724	73.058	-0.5817	1.0852	3.6631
6	76.723	73.058	-0.5772	-0.3536	0.0849
						7	76.723	73.058	-0.5828	-0.6531	-0.6587

Attitude Eulerian angles deviation from last table as can be seen, the outside attitude error that obtains can satisfy job requirement fully less than 1 rad.

The computation process of parameter vector d2 sees the following form:

Iterations	k1	g1	g2c	g3	g4
						1	-1.1	0.035399	0.43	-0.12328	0.024027
2	-0.40897	0.027119	0.079602	-0.045578	0.016899
						3	-0.31141	0.022141	0.0089006	-0.025147	0.011051
4	-0.29845	0.020565	-0.005874	-0.020933	0.0098816

5	-0.29733	0.020158	-0.008952	-0.020081	0.0096507
						6	-0.29744	0.020059	-0.009594	-0.019909	0.0096042
7	-0.29753	0.020035	-0.009728	-0.019874	0.0095947

Final asterism coordinate estimated valve and deviation of measuring value root of mean square are x direction 0.0571, y direction 0.0575 pixel.Because noise level is 0.05 pixel, therefore final parameter can satisfy the calibration accuracy requirement.

Claims (1)

1, based on the star sensor of weng model at the rail alignment method, it is characterized in that the step of calibration is as follows:

1.1, set up star sensor attitude transition matrix;

1.1.1, set up system of celestial coordinates and star sensor system of axes; If O '-XnYnZn is a system of celestial coordinates, O-XYZ is the star sensor system of axes, and the attitude angle of star sensor is by right ascension α ₀, declination β ₀, and roll angle φ ₀Form α ₀Be the projection of Z axle on the XnYn face and the angle of Xn axle; β ₀Be Z axle and its angle between the projection on the XnYn face; φ ₀Angle for meridian plane and XY hand-deliver line and Y-axis;

1.1.2, set up star sensor attitude transition matrix; The switching process that is tied to the star sensor system of axes from celestial coordinates is: the first step, rotate around the Zn axle

The angle makes the Xn axle vertical with meridian plane; In second step, rotate around new Xn axle

The angle makes the Zn axle consistent with the Z axle; In the 3rd step, rotate φ around new Zn axle ₀The angle makes system of celestial coordinates and star sensor system of axes overlap.Then the direction cosine matrix R of this rotation process is the attitude transition matrix, is expressed as:

In the formula, s represents sinusoidal function, and c represents cosine function;

R can be simplified shown as: $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right],$ r ₁～r ₉9 elements of expression R matrix;

1.2, set up star sensor weng imaging model;

If the starlight direction vector is: $w=\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\&beta;}\cos \mathrm{\&alpha;}\\ \cos \mathrm{\&beta;}\sin \mathrm{\&alpha;}\\ \sin \mathrm{\&beta;}\end{array}\right)---\left[2\right]$

In the formula, w1, w2 and w3 represent the x of starlight vector under the star sensor system of axes, y, z coordinate figure; α, β are right ascension, the declination coordinate of fixed star;

1.2.1, use star sensor that starlight is carried out the radial distortion imaging; O-XYZ is the star sensor system of axes, and o-xy is 2 dimension target surface system of axess, and ∏ is a target surface, and distance is a focal length between the Oo; The desirable perspective imaging point of starlight vector w is P ', because radial distortion takes place, actual imaging point is P;

1.2.2, set up the linear imaging model; Owing to have uncertain factor of proportionality between x and the y direction, use f _uAnd f _vThe focal length of representing x and y direction respectively, set up the linear imaging model and be:

$\frac{x^{\&prime;}}{f_u}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

$\frac{y^{\&prime;}}{f_v}=\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$ [3]

In the formula, x ', y ' defines linear parameter vector d for the asterism position of perspective projection ₁=(f _uf _vα ₀β ₀_ ₀);

1.2.3, consider distortion effects, set up weng nonlinear distortion model formation and be:

$\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}=\hat{u}+(g_1+g_3){\hat{u}}^2+g_4\hat{u}\hat{v}+g_1{\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

$\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}=\hat{v}+g_2{\hat{u}}^2+g_3\hat{u}\hat{v}+(g_2+g_4){\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$ [4]

In the formula, $\hat{u}=\frac{x}{f_u},\hat{v}=\frac{y}{f_v},$ X, y are asterism actual measurement coordinate, and k1 is a coefficient of radial distortion, and g1～g4 is the aggregative formula coefficient of decentering distortion and thin prism distortion, so nonlinear distortion parameter has 5, and defined parameters vector d ₂=(k ₁g ₁g ₂g ₃g ₄);

1.3, data acquisition; According to the frame star chart that star sensor collects, whole asterism data acquisitions that definition collects are Ω ₀, the asterism number is n; With the principal point is the center, and diameter is that picture size 1/2 is done circle, and the data acquisition in the note circle is Ω ₁, the asterism number is n ₁

1.4, calibrate by calculation of parameter;

Calibration process be divided into two the step finish, according to centre data Ω ₁Be subjected to the less principle of distortion effects, the first step is utilized Ω ₁Construct linear computation model, calculate the initial value of parameter vector d1; The parameter of second step to d1 is optimized processing; In the 3rd step, according to the weng model, structure distortion factor linear equation is at global data Ω ₀Following calculating parameter vector d2; At last, whether judgment data result meets the demands, if do not satisfy then return the second stepping row iteration and calculate, if satisfied the output parameter vector finally separate d1 and d2; Concrete computation process is as follows:

1.4.1, carry out linear dimensions and estimate; Have according to formula [3]:

x′w ₁r ₇+x′w ₂r ₈+x′w ₃r ₉＝w ₁r ₁f _u+w ₂r ₂f _u+w ₃r ₃f _u

y′w ₁r ₇+y′w ₂r ₈+y′w ₃r ₉＝w ₁r ₄f _v+w ₂r ₅f _v+w ₃r ₆f _v [5]

Arrangement has:

$\begin{array}{c}{w}_1\frac{r_1{f}_u}{r_9}+w_2\frac{r_2{f}_u}{r_9}+w_3\frac{r_3{f}_u}{r_9}-{\mathrm{xw}}_1\frac{r_7}{r_9}-{\mathrm{xw}}_2\frac{r_8}{r_9}=x_{w3}\end{array}$

$w_1\frac{r_4{f}_v}{r_9}+w_2\frac{r_5{f}_v}{r_9}+w_3\frac{r_6{f}_v}{r_9}-y{w}_1\frac{r_8}{r_9}-y{w}_2\frac{r_8}{r_9}={\mathrm{yw}}_3$ [6]

Suppose $r_1^{\&prime;}=\frac{r_1{f}_u}{r_9},r_2=\frac{r_2{f}_u}{r_9},r_3^{\&prime;}=\frac{r_3{f}_u}{r_9},r_4^{\&prime;}=\frac{r_4{f}_v}{r_9},r_5^{\&prime;}=\frac{r_5{f}_u}{r_9},r_6=\frac{r_6{f}_v}{r_9},r_7^{\&prime;}=\frac{r_7}{r_9},r_8^{\&prime;}=\frac{r_8}{r_9},$ Totally 8 unknown numberes, defined parameters vector m=[r ₁' r ₂' r ₃' r ₄' r ₅' r ₆' r ₇' r ₈'] ^T, can obtain equation for each asterism data:

[w ₁ w ₂ w ₃ 0 0 0 -xw ₁ -xw ₂]m＝xw ₃

[0 0 0 w ₁ w ₂ w ₃ -yw ₁ -yw ₂]m＝yw ₃[7]

For data acquisition Ω ₁Interior asterism can obtain 2n ₁Individual equation constitutes " Am=b " form; A is 2n ₁The matrix that individual equation levoform constitutes, b is 2n ₁The vector that the right formula of individual equation constitutes; Available method of least square makes

Solve 8 unknown number r ₁'～r ₈'; Be the orthogonal matrix characteristics according to R then, have:

$r_9=\sqrt{1+r_7^{\&prime;}+r_8^{\&prime;}},r_7=r_7^{\&prime;}\&CenterDot;r_9,r_8=r_8^{\&prime;}\&CenterDot;r_9$

$f_u=r_9\&CenterDot;\sqrt{r_1^{\&prime;2}+r_2^{\&prime;2}+r_3^{\&prime;2}}$

$r_1=r_1^{\&prime;}\&CenterDot;r_9/f_u,r_2=r_2^{\&prime;}\&CenterDot;r_9/f_u,r_3=r_3^{\&prime;}\&CenterDot;r_9/f_u$

$f_v=r_9\&CenterDot;\sqrt{r_4^{\&prime;2}+r_5^{\&prime;2}+r_6^{\&prime;2}}$

$r_4=r_4^{\&prime;}\&CenterDot;r_9/f_v,r_5=r_5^{\&prime;}\&CenterDot;r_9/f_v,r_6=r_6^{\&prime;}\&CenterDot;r_9/f_v$ [8]

Suppose r earlier ₉For just, obtain the attitude estimated matrix, select a certain distance center fixed star far away in the visual field to judge r then, whether really for just,, r is described then if starlight vector coordinate after the imaging of attitude estimated matrix is consistent with meeting of measurement coordinate ₉For just, otherwise for negative;

1.4.2, optimize parameter vector d1;

If the actual value of linear dimensions vectorial sum nonlinear parameter vector is d ₁, d ₂, estimated valve is Optimum estimated valve should satisfy:

$\underset{d}{\min }\left|\right|P-f({\hat{d}}_1,{\hat{d}}_2)\left|\right|$

In the formula, P is point set Ω ₁The asterism coordinate; With separating of obtaining of formula [8] is initial value, adopts nonlinear optimization method that d1 is optimized; If:

x′＝f _x(d ₁， d ₂)

y′＝f _y(d ₁， d ₂)[9]

In the formula, d ₂Expression parameter vector d2 fixes; Because f _xAnd f _yBe nonlinear function, therefore adopt non-linear linearization method to come estimated parameter vector d1; Suppose

Be estimated valve, x, y are observed reading, Δ d ₁Be vectorial estimated bias;

$\mathrm{\&Delta;x}=\stackrel{\&OverBar;}{x}-\hat{x}\&ap;\mathrm{A\&Delta;}{d}_1$

$\mathrm{\&Delta;y}=\stackrel{\&OverBar;}{y}-\hat{y}\&ap;\mathrm{B\&Delta;}{d}_1$ [10]

In the formula, A, B are sensitive matrix, ask every partial derivative to obtain;

So obtain parameter vector updating value d ₁=d ₁+ Δ d ₁

1.4.3, calculate distortion factor vector d2;

After obtaining linear dimensions, next step finds the solution nonlinear parameter d ₂According to formula [4], construct linear estimation model:

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& {\hat{u}}^2+{\hat{v}}^2& 0& {\hat{u}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=\mathrm{du}$

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& 0& {\hat{u}}^2+{\hat{v}}^2& {\hat{v}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=\mathrm{dv}$ [11]

In the formula: $\mathrm{du}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}-\hat{u}$

$\mathrm{dv}=\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}-\hat{v}$

Constitute " Bd ₂=c " form, B is the matrix that 2n equation levoform constitutes, b is the vector that the right formula of 2n equation constitutes; Utilize method of least square to make

Obtain the estimated valve of distortion parameter vector d2, the data here are whole point set Ω ₀

After obtaining whole parameter d 1 and d2 of model, judge by contrast ideal coordinates value and measurement coordinate figure; If error is still bigger, then get back to step 1.4.2 and further optimize d1, carry out iterative computation; If error meets the demands, then stop iteration, the output result.

Images