CN102901519A - optical push-broom satellite in-orbit stepwise geometric calibration method based on probe element direction angle - Google Patents

Abstract

The invention discloses an optical push-broom satellite in-orbit stepwise geometric calibration method based on a probe element direction angle. A probe element direction angle-based in-orbit geometric calibration model and a stepwise calibration resolving method are included. According to the method, a probe element direction angle model is used as a camera calibration model; calibration parameters comprise internal calibration parameters and external calibration parameters; and the external calibration parameters are resolved first and then the internal calibration parameters are resolved by using a stepwise calibration method, so that the condition that a reliable calibration result cannot be obtained due to strong correlation existing between the calibration parameters is avoided. By the effective optical push-broom satellite in-orbit geometric calibration method, in-orbit geometric calibration can be effectively performed on an optical push-broom satellite, and the geometric quality of image products is improved. Practice proves that the calibration method is feasible and effective, and the calibration result is stable, reliable and high in precision.

Description

A kind of pushing away based on spy unit sensing angle optics swept satellite in rail substep geometric calibration method

Technical field

The invention belongs to the remotely sensing image geometric process field, relate to a kind of pushing away based on spy unit sensing angle optics and sweep satellite in rail substep geometric calibration method.

Background technology

The key link of giving full play to high-resolution optical remote sensing satellite geometry performance at the rail geometric calibration, significant to improving how much quality of satellite image product.Owing to being subjected to the impact of the various factorss such as the change of release, imaging circumstances of stress in the satellite launch process and device aging, how much imaging parameters of satellite image change, the Laboratory Calibration result is no longer applicable, must be by again obtain the exact value of these parameters at the rail geometric calibration, otherwise, the geometry difficult quality guarantee of satellite image product can't satisfy the requirement of subsequent treatment and application, and its using value is reduced greatly.

At present, carry out a large amount of research and put into practice work for the geometric calibration of close shot camera, aerial camera, then studied lessly for the geometric calibration of optics linear array push-broom type satellite, not yet formed ripe theory and method.Because the high rail of spaceborne line-scan digital camera, narrow field angle and linear array push such as sweep at the special imaging characteristic, so that have strong correlation between how much imaging parameters, if adopt calibration model and the method for close shot camera, aerial camera, stable, reliable, high-precision the calibration results be must be difficult to obtain, therefore must suitable calibration model and method be made up for the characteristics of spaceborne line-scan digital camera.

Summary of the invention

Problem to be solved by this invention is for optics linear array push-broom type satellite, to provide a kind of effectively in rail geometric calibration method.

Technical scheme of the present invention is that a kind of pushing away based on spy unit sensing angle optics swept satellite in rail substep geometric calibration method, may further comprise the steps:

Step 1 is utilized the reference mark data, carries out the reference mark at image to be calibrated and measures, and obtains the reference mark measurement information.

Step 2 utilizes auxiliary data and Laboratory Calibration parameter to make up optics linear array push-broom type satellite based on the geometric calibration model of visiting unit sensing angle;

Step 3 utilizes the substep calibrating method to resolve inside and outside calibration parameter, obtains the calibration results; Implementation is, at first resolves outer calibration parameter based on least square adjustment, recovers camera coordinates and ties up to attitude in the space; Then on this basis, resolve the internal calibration parameter based on least square adjustment, each visits the sensing angle of unit under camera coordinates system to determine CCD.

And the described optics linear array of step 2 push-broom type satellite is as follows based on the geometric calibration model of visiting unit sensing angle,

$\left(\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ \tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\\ 1\end{array}\right)={\mathrm{\λR}}_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw}){(R_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}})---\left(1\right)$

$\left\{\begin{array}{c}{\mathrm{\ψ}}_x\left(s\right)={\mathrm{ax}}_0+{\mathrm{ax}}_1\×s+{\mathrm{ax}}_2\×s^2+{\mathrm{ax}}_3\×s^3\\ {\mathrm{\ψ}}_y\left(s\right)={\mathrm{ay}}_0+{\mathrm{ay}}_1\×s+{\mathrm{ay}}_2\×s^2+{\mathrm{ay}}_3\×s^3\end{array}\right.$

In the formula, (X _g, Y _g, Z _g) and (X _Gps, Y _Gps, Z _Gps) represent respectively object space point and the coordinate of gps antenna phase center under the WGS84 coordinate system that picture point is corresponding;

Represent respectively the rotation matrix that rotation matrix that the WGS84 coordinate is tied to the J2000 coordinate system, rotation matrix, satellite body coordinate that the J2000 coordinate is tied to the satellite body coordinate system are tied to camera coordinates system; (B _X, B _Y, B _Z) _BodyEccentric vector the coordinate under satellite body coordinate system of representative from the sensor projection centre to the gps antenna phase center; (t) the expression parameter current is a time dependent amount; (ψ _x(s), ψ _y(s)) represent the sensing angle of the first s of spy under camera coordinates system, the s representative is visited first number;

In above geometric calibration model, parameter to be calibrated is divided into outer calibration parameter X _EWith the internal calibration parameter X _I, outer calibration parameter Pitch, roll, yaw are respectively pitching, rolling and yaw direction angle; The internal calibration parameter X _I=(ax ₀, ax ₁, ax ₂, ax ₃, ay ₀, ay ₁, ay ₂, ay ₃), ax ₀, ax ₁, ax ₂, ax ₃, ay ₀, ay ₁, ay ₂, ay ₃For internal calibration is visited the coefficient that unit points to angle model.

And step 3 comprises following substep,

Step 3.1 is located at and has measured K equally distributed Ground Nuclear Magnetic Resonance reference mark on the image to be calibrated as orientation point, is designated as RP _i, i=1,2,3 ..., K; The object space point that each orientation point is corresponding and picture side's point are designated as respectively RPG _iAnd RPM _i, object space point RPG _iThe WGS84 geocentric rectangular coordinate be (X _i, Y _i, Z _i), as side point RPM _iImage coordinate be (s _i, l _i);

Step 3.2 makes in the formula (1):

$\left\{\begin{array}{c}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]\left(\begin{array}{c}{R}_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}{\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}}\end{array}\right)\\ R_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw})=\left[\begin{array}{c}{a}_1,b_1,c_1\\ a_2,b_2,c_2\\ a_3,b_3,c_3\end{array}\right]\end{array}\right.$

Formula (1) is converted into formula (2):

$\left\{\begin{array}{c}F(X_E,X_I)=\frac{a_1\mathrm{Ux}+b_1\mathrm{Uy}+c_1\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ G(X_E,X_I)=\frac{a_2\mathrm{Ux}+b_2\mathrm{Uy}+c_2\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\end{array}\right.---\left(2\right)$

In the following formula, vector $\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]$ Be object space vector U, vector the coordinate under body coordinate system of representative from the camera projection centre to object space point; a ₁, b ₁, c ₁a ₂, b ₂, c ₂a ₃, b ₃, c ₃Represent respectively camera 9 elements of matrix are installed; F (X _E, X _I), G (X _E, X _I) be respectively along rail and point to angle residual error and vertical rail sensing angle residual error;

Step 3.3 is externally calibrated parameter X E and internal calibration parameter X I initialize

Step 3.4 is considered as true value with the currency of internal calibration parameter X I, outer calibration parameter X E is considered as unknown parameter to be asked, with the currency of internal calibration parameter X I and outer calibration parameter X E

Substitution formula (2) to each orientation point RPi, is carried out linearization process to formula (2), sets up error equation (3),

V _i=A _iX-L _i, weigh and be P _i(3)

Wherein $A_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_E}\\ \frac{\∂G_i}{\∂X_E}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂\mathrm{pitch}}\frac{\∂F_i}{\∂\mathrm{roll}}\frac{\∂F_i}{\∂\mathrm{yaw}}\\ \frac{\∂G_i}{\∂\mathrm{pitch}}\frac{\∂G_i}{\∂\mathrm{roll}}\frac{\∂G_i}{\∂\mathrm{yaw}}\end{array}\right]$ $X=d{X}_E=\left[\begin{array}{c}\mathrm{dpitch}\\ \mathrm{droll}\\ \mathrm{dyaw}\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

In the formula, L _iTo utilize inside and outside calibration parameter currency

The error vector that substitution formula (2) calculates; A _iIt is the matrix of coefficients of error equation; The outer calibration of X representative parameter correction dX _E=(dpitch, droll, dyaw), d represents the correction symbol; P _iCurrent orientation point RP _iPower corresponding to picture point measurement accuracy; F _iAnd G _iBe respectively along rail and point to angle residual error F (X _E, X _I), the rail that hangs down points to angle residual error G (X _E, X _I) function model, obtain corresponding error equation behind the differential;

By formula (4) computing method equation coefficient matrix,

$A^T\mathrm{PA}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{A}_i$

（4）

$A^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{L}_i$

In the following formula, matrix $A=\left[\begin{array}{c}{A}_1\\ ...\\ A_i\\ ...\\ A_K\end{array}\right],$ Matrix

Matrix $L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

Utilize least square adjustment to calculate X, suc as formula (5),

X＝(A ^TPA) ^-1(A ^TPL) （5）

Utilize formula (6) to upgrade the currency of outer calibration parameter X E, then return execution in step 4 iterative computation, enter step 3.4 after the iteration stopping;

$X_E=X_E^o+X---\left(6\right)$

Step 3.5 is resolved the internal calibration parameter, and the currency of step 3.4 gained being calibrated parameter X E outward is considered as true value, and internal calibration parameter X I then is considered as unknown parameter to be asked, with the currency of internal calibration parameter X I and outer calibration parameter X E

Substitution formula (2) to each orientation point RPi, carries out linearization process to formula (2), sets up error equation (7),

V _i=B _iY-L _iPower is P _i(7)

Wherein $B_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_I}\\ \frac{\∂G_i}{\∂X_I}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂{\mathrm{ax}}_0}\frac{\∂F_i}{\∂{\mathrm{ax}}_1}\frac{\∂F_i}{\∂{\mathrm{ax}}_2}\frac{\∂F_i}{\∂{\mathrm{ax}}_3}\frac{\∂F_i}{\∂{\mathrm{ay}}_0}\frac{\∂F_i}{\∂{\mathrm{ay}}_1}\frac{\∂F_i}{\∂a{y}_2}\frac{\∂F_i}{\∂{\mathrm{ay}}_3}\\ \frac{\∂G_i}{\∂{\mathrm{ax}}_0}\frac{\∂G_i}{\∂{\mathrm{ax}}_1}\frac{\∂G_i}{\∂{\mathrm{ax}}_2}\frac{\∂G_i}{\∂{\mathrm{ax}}_3}\frac{\∂G_i}{\∂{\mathrm{ay}}_0}\frac{\∂G_i}{\∂{\mathrm{ay}}_1}\frac{\∂G_i}{\∂{\mathrm{ay}}_2}\frac{\∂G_i}{\∂{\mathrm{ay}}_3}\end{array}\right]$ $Y={\mathrm{dX}}_I=\left[\begin{array}{c}{\mathrm{dax}}_0\\ {\mathrm{dax}}_1\\ {\mathrm{dax}}_2\\ {\mathrm{dax}}_3\\ {\mathrm{day}}_0\\ {\mathrm{day}}_1\\ {\mathrm{day}}_2\\ {\mathrm{day}}_3\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

In the formula, L _iTo utilize inside and outside calibration parameter currency

The error vector that substitution formula (2) calculates; B _iIt is the matrix of coefficients of error equation; Y represents internal calibration parameter correction dX _I, d represents the correction symbol; P _iCurrent orientation point RP _iPower corresponding to picture point measurement accuracy; F _iAnd G _iBe respectively along rail and point to angle residual error F (X _E, X _I), the rail that hangs down points to angle residual error G (X _E, X _I) function model, obtain corresponding error equation behind the differential;

By formula (8) computing method equation coefficient matrix;

$B^T\mathrm{PB}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{B}_i$

（8）

$B^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{L}_i$

In the following formula, $B=\left[\begin{array}{c}{B}_1\\ ...\\ B_i\\ ...\\ B_K\end{array}\right],$

$L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

Utilize least square adjustment to calculate Y, suc as formula (9);

Y＝(B ^TPB) ^-1(B ^TPL) （9）

Utilize formula (10) to upgrade the currency of internal calibration parameter X I, then return execution in step 5 iterative computation, enter step 3.6 after the iteration stopping;

$X_I=X_I^o+Y---\left(10\right)$

Step 3.6 is calibrated the currency X of parameter outward according to step 3.4 and step 3.5 gained _EWith the internal calibration parameter X _ICurrency, upgrade the camera parameter file.

The invention has the advantages that: 1. what determine based on the calibration model of visiting unit and point to the angle no longer is that each visits unit position in viewing field of camera, but its sensing angle under camera coordinates system, can simplify the internal calibration model on the one hand, avoid between the inner various distortion errors of camera because the overparameterization problem that strong correlation causes; Simultaneously, because outer calibration parameter recovery camera coordinates ties up to the attitude (sensing of three axles) in the space, the internal calibration parameter is determined each spy unit sensing angle under camera coordinates system, outer calibration parameter provides reference data for the internal calibration parameter, the internal calibration parameter is calibrated outside under the reference data that parameter provides and is determined, inside and outside calibration parameter is determined respectively to visit the absolute sensing of unit in the space jointly, can avoid like this relativity problem between the inside and outside calibration parameter, the inside and outside calibration of raising parameter calculation stability; 2. adopt the substep calibrating method, resolving first outer calibration parameter recovers camera coordinates and ties up to attitude in the space, determine to provide reference data for the internal calibration parameter, then resolve on this basis the internal calibration parameter and determine each spy unit sensing angle under camera coordinates system (reference data), can avoid calibrating parameter and resolve simultaneously the equation ill-conditioning problem that causes owing to strong correlation, be conducive to obtain stable, reliable, high-precision the calibration results.

Description of drawings

Fig. 1 is the schematic flow sheet of the embodiment of the invention;

Fig. 2 points to the angle schematic diagram for visiting unit.

Embodiment

Describe the specific embodiment of the invention in detail below in conjunction with drawings and Examples.Referring to Fig. 1, the flow process of embodiment can be divided into three steps, and concrete grammar, formula and flow process that each step is implemented are as follows:

1. utilize the reference mark data, measurement control dot information on image to be calibrated.In order to guarantee the calculation accuracy of the calibration results, provide suggestion for the reference mark that measures in quantity and distribution: in image to be calibrated, the reference mark that measures is along being distributed on the rail direction in the shorter zone as far as possible, and the rail direction of hanging down is then answered the whole CCD scope of uniform fold.The number of control points aspect is many as far as possible under the prerequisite of reasonable cost.

2. utilize auxiliary data and the Laboratory Calibration parameters such as appearance rail, imaging time, make up optics linear array push-broom type satellite based on visiting geometric calibration model that unit points to the angle suc as formula (1):

$\left(\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ \tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\\ 1\end{array}\right)={\mathrm{\λR}}_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw}){(R_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}})---\left(1\right)$

$\left\{\begin{array}{c}{\mathrm{\ψ}}_x\left(s\right)={\mathrm{ax}}_0+{\mathrm{ax}}_1\×s+{\mathrm{ax}}_2\×s^2+{\mathrm{ax}}_3\×s^3\\ {\mathrm{\ψ}}_y\left(s\right)={\mathrm{ay}}_0+{\mathrm{ay}}_1\×s+{\mathrm{ay}}_2\×s^2+{\mathrm{ay}}_3\×s^3\end{array}\right.$

In the formula, (X _g, Y _g, Z _g) and (X _Gps, Y _Gps, Z _Gps) represent respectively object space point and the coordinate of gps antenna phase center under the WGS84 coordinate system that picture point is corresponding, wherein, the latter is obtained by the GPS that carries on the satellite.

Represent respectively the rotation matrix that rotation matrix that the WGS84 coordinate is tied to the J2000 coordinate system, rotation matrix, satellite body coordinate that the J2000 coordinate is tied to the satellite body coordinate system are tied to camera coordinates system; (B _X, B _Y, B _Z) _BodyEccentric vector the coordinate under satellite body coordinate system of representative from the sensor projection centre to the gps antenna phase center; Wherein,

, gyro quick by star obtains (B by integrated attitude determination _X, B _Y, B _Z) _BodyThen before satellite launch by the laboratory calibration, (t) expression parameter current be a time dependent amount; (ψ _x(s), ψ _y(s)) represent the sensing angle of the first s of spy under camera coordinates system, the s representative is visited first number, and as shown in Figure 2: camera coordinates is X ₁Y ₁Z ₁Initial point be designated as O ₁, visit first V _ImageThe sensing angle be (ψ _x, ψ _y); In this calibration model, parameter to be calibrated is divided into outer calibration parameter X _EWith the internal calibration parameter X _I, wherein, outer calibration parameter

(pitch, roll, yaw are respectively pitching, rolling and yaw direction angle) is used for compensation camera error of fixed angles, and the recovery camera coordinates ties up to the sensing in the space, is the definite reference data of resolving of internal calibration parameter; The internal calibration parameter X _I=(ax ₀, ax ₁, ax ₂, ax ₃, ay ₀, ay ₁, ay ₂, ay ₃) then be used for compensation because the picture point error that the inner various distortion of camera cause determines that CCD respectively visits the sensing angle of unit under camera coordinates system (reference data).Internal calibration parameter and outer calibration parameter are recovered CCD jointly, and each visits the absolute sensing of unit in the space.Be with conventional camera calibration model difference, what the calibration model based on visiting unit sensing angle that the present invention proposes was determined no longer is that each visits the position of unit in viewing field of camera, but its sensing angle under camera coordinates system, can avoid like this calibrating between the parameter can't obtain to stablize owing to there being strong correlation, reliable the calibration results.

3. utilize the substep calibrating method to resolve inside and outside calibration parameter, the inside and outside calibration of step solution parameter, outer calibration parameter provides reference data for the internal calibration parameter, and the internal calibration parameter is calibrated outside under the reference data that parameter provides and is determined.Its principle is: at first resolve outer calibration parameter based on least square adjustment, the recovery camera coordinates ties up to the attitude in the space; Then on this basis, resolve the internal calibration parameter based on least square adjustment, each visits the sensing angle of unit under camera coordinates system to determine CCD.

Solution formula and the step of embodiment are as follows:

Step 3.1 is supposed to have measured K equally distributed Ground Nuclear Magnetic Resonance reference mark as orientation point at image to be calibrated, and is designated as RP _i, here, i=1,2,3 ..., K; The object space point that it is corresponding and picture side's point are designated as respectively RPG _iAnd RPM _i, object space point RPG _iThe WGS84 geocentric rectangular coordinate be (X _i, Y _i, Z _i), as side point RPM _iImage coordinate be (s _i, l _i);

Step 3.2 makes in the formula (1):

$\left\{\begin{array}{c}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]\left(\begin{array}{c}{R}_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}{\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}}\end{array}\right)\\ R_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw})=\left[\begin{array}{c}{a}_1,b_1,c_1\\ a_2,b_2,c_2\\ a_3,b_3,c_3\end{array}\right]\end{array}\right.$

Formula (1) can be converted into formula (2):

$\left\{\begin{array}{c}F(X_E,X_I)=\frac{a_1\mathrm{Ux}+b_1\mathrm{Uy}+c_1\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ G(X_E,X_I)=\frac{a_2\mathrm{Ux}+b_2\mathrm{Uy}+c_2\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\end{array}\right.---\left(2\right)$

In the following formula, vector $\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]$ Be object space vector U, vector the coordinate under body coordinate system of representative from the camera projection centre to object space point; a ₁, b ₁, c ₁a ₂, b ₂, c ₂a ₃, b ₃, c ₃Represent respectively camera 9 elements of matrix are installed; F (X _E, X _I), G (X _E, X _I) be respectively along rail and point to angle residual error and vertical rail sensing angle residual error;

Step 3.3 is externally calibrated parameter X _E, the internal calibration parameter X _IInitialize

Here be the Laboratory Calibration value.

Step 3.4 is with the internal calibration parameter X _ICurrency be considered as " true value ", with outer calibration parameter X _EBe considered as unknown parameter to be asked.With their currency

Substitution formula (2) is to each orientation point RP _i, formula (2) is carried out linearization process, (subscript i representative utilizes orientation point RP to set up error equation (3) _iThe error equation of setting up):

Vi=A _iX-L _i, weigh and be P _i(3)

Wherein $A_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_E}\\ \frac{\∂G_i}{\∂X_E}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂\mathrm{pitch}}\frac{\∂F_i}{\∂\mathrm{roll}}\frac{\∂F_i}{\∂\mathrm{yaw}}\\ \frac{\∂G_i}{\∂\mathrm{pitch}}\frac{\∂G_i}{\∂\mathrm{roll}}\frac{\∂G_i}{\∂\mathrm{yaw}}\end{array}\right]$ $X=d{X}_E=\left[\begin{array}{c}\mathrm{dpitch}\\ \mathrm{droll}\\ \mathrm{dyaw}\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

In the formula, L _iTo utilize inside and outside calibration parameter currency

The error vector that substitution formula (2) calculates; A _iIt is the matrix of coefficients of error equation; The outer calibration of X representative parameter correction dX _E, i.e. (dpitch, droll, dyaw), d represents the correction symbol; P _iCurrent orientation point RP _iPower corresponding to picture point measurement accuracy; F _iAnd G _iBe respectively along rail and point to angle residual error F (X _E, X _I), the rail that hangs down points to angle residual error G (X _E, X _I) function model, obtain corresponding error equation behind the differential.

By formula (4) computing method equation coefficient matrix,

$A^T\mathrm{PA}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{A}_i$

（4）

$A^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{L}_i$

In the following formula, matrix $A=\left[\begin{array}{c}{A}_1\\ ...\\ A_i\\ ...\\ A_K\end{array}\right],$ Matrix

Matrix $L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

Utilize least square adjustment to calculate X, suc as formula (5);

X＝(A ^TPA) ^-1(A ^TPL) （5）

Utilize formula (6) to upgrade outer calibration parameter X _ECurrency, then return execution in step 4 iterative computation.3 parameters (dpitch, droll, dyaw) are all less than threshold value 10 in outer calibration parameter correction X ^-12The time, iteration stopping.

$X_E=X_E^o+X---\left(6\right)$

Step 3.5 is resolved the internal calibration parameter, and to resolve internal calibration parameter implementation consistent with step 4, and this moment, the currency with outer calibration parameter was considered as " true value ", and the internal calibration parameter then is considered as unknown parameter to be asked, with the currency of inside and outside calibration parameter

Substitution formula (2), at this moment

For step 3.4 gained is calibrated parameter X outward _ECurrency.The internal calibration parameter X _IWith outer calibration parameter X _ETo each orientation point RP _i, formula (2) is carried out linearization process, (subscript i representative utilizes orientation point RP to set up error equation (7) _iThe error equation of setting up): concrete solution formula and process are as follows:

V _i=B _iY-L _iPower is P _i(7)

Wherein $B_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_I}\\ \frac{\∂G_i}{\∂X_I}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂{\mathrm{ax}}_0}\frac{\∂F_i}{\∂{\mathrm{ax}}_1}\frac{\∂F_i}{\∂{\mathrm{ax}}_2}\frac{\∂F_i}{\∂{\mathrm{ax}}_3}\frac{\∂F_i}{\∂{\mathrm{ay}}_0}\frac{\∂F_i}{\∂{\mathrm{ay}}_1}\frac{\∂F_i}{\∂a{y}_2}\frac{\∂F_i}{\∂{\mathrm{ay}}_3}\\ \frac{\∂G_i}{\∂{\mathrm{ax}}_0}\frac{\∂G_i}{\∂{\mathrm{ax}}_1}\frac{\∂G_i}{\∂{\mathrm{ax}}_2}\frac{\∂G_i}{\∂{\mathrm{ax}}_3}\frac{\∂G_i}{\∂{\mathrm{ay}}_0}\frac{\∂G_i}{\∂{\mathrm{ay}}_1}\frac{\∂G_i}{\∂{\mathrm{ay}}_2}\frac{\∂G_i}{\∂{\mathrm{ay}}_3}\end{array}\right]$ $Y={\mathrm{dX}}_I=\left[\begin{array}{c}{\mathrm{dax}}_0\\ {\mathrm{dax}}_1\\ {\mathrm{dax}}_2\\ {\mathrm{dax}}_3\\ {\mathrm{day}}_0\\ {\mathrm{day}}_1\\ {\mathrm{day}}_2\\ {\mathrm{day}}_3\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

In the formula, L _iTo utilize inside and outside calibration parameter currency

The error vector that substitution formula (2) calculates; B _iIt is the matrix of coefficients of error equation; Y represents internal calibration parameter correction dX _I, d represents the correction symbol; P _iCurrent orientation point RP _iPower corresponding to picture point measurement accuracy; F _iAnd G _iBe respectively along rail and point to angle residual error F (X _E, X _I), the rail that hangs down points to angle residual error G (X _E, X _I) function model, obtain corresponding error equation behind the differential.

By formula (8) computing method equation coefficient matrix;

$B^T\mathrm{PB}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{B}_i$

（8）

$B^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{L}_i$

In the following formula, $B=\left[\begin{array}{c}{B}_1\\ ...\\ B_i\\ ...\\ B_K\end{array}\right],$

$L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

Utilize least square adjustment to calculate Y, suc as formula (9);

Y＝(B ^TPB) ^-1(B ^TPL) （9）

Utilize formula (10) to upgrade the internal calibration parameter X _ICurrency, then return execution in step 5 iterative computation.8 parameter (dax in internal calibration parameter correction Y ₀, dax ₁, dax ₂, dax ₃, day ₀, day ₁, day ₂, day ₃) all less than threshold value 10 ^-12The time, iteration stopping.

$X_I=X_I^o+Y---\left(10\right)$

Step 3.6, according to the calibration results that obtains, namely inside and outside calibration parameter currency upgrades the camera parameter file.

Specific embodiment described herein only is to the explanation for example of the present invention's spirit.Those skilled in the art can make various modifications or replenish or adopt similar mode to substitute described specific embodiment, but can't depart from spirit of the present invention or surmount the defined scope of appended claims.

Claims (3)

1. one kind is pointed to angle optics and pushes away and sweep satellite in rail substep geometric calibration method based on visiting unit, may further comprise the steps:

Step 1 is utilized the reference mark data, carries out the reference mark at image to be calibrated and measures, and obtains the reference mark measurement information.

Step 2 utilizes auxiliary data and Laboratory Calibration parameter to make up optics linear array push-broom type satellite based on the geometric calibration model of visiting unit sensing angle;

Step 3 utilizes the substep calibrating method to resolve inside and outside calibration parameter, obtains the calibration results; Implementation is, at first resolves outer calibration parameter based on least square adjustment, recovers camera coordinates and ties up to attitude in the space; Then on this basis, resolve the internal calibration parameter based on least square adjustment, each visits the sensing angle of unit under camera coordinates system to determine CCD.

2. push away based on the first sensing of spy angle optics as claimed in claim 1 and sweep star in rail substep geometric calibration method, it is characterized in that: the described optics linear array of step 2 push-broom type satellite is as follows based on the geometric calibration model of visiting first sensing angle,

$\left(\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ \tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\\ 1\end{array}\right)={\mathrm{\λR}}_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw}){(R_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}})---\left(1\right)$

$\left\{\begin{array}{c}{\mathrm{\ψ}}_x\left(s\right)={\mathrm{ax}}_0+{\mathrm{ax}}_1\×s+{\mathrm{ax}}_2\×s^2+{\mathrm{ax}}_3\×s^3\\ {\mathrm{\ψ}}_y\left(s\right)={\mathrm{ay}}_0+{\mathrm{ay}}_1\×s+{\mathrm{ay}}_2\×s^2+{\mathrm{ay}}_3\×s^3\end{array}\right.$

In the formula, (X _g, Y _g, Z _g) and (X _Gps, Y _Gps, Z _Gps) represent respectively object space point and the coordinate of gps antenna phase center under the WGS84 coordinate system that picture point is corresponding;

Represent respectively the rotation matrix that rotation matrix that the WGS84 coordinate is tied to the J2000 coordinate system, rotation matrix, satellite body coordinate that the J2000 coordinate is tied to the satellite body coordinate system are tied to camera coordinates system; (B _X, B _Y, B _Z) _BoayEccentric vector the coordinate under satellite body coordinate system of representative from the sensor projection centre to the gps antenna phase center; (t) the expression parameter current is a time dependent amount; (ψ _x(s), ψ _y(s)) represent the sensing angle of the first s of spy under camera coordinates system, the s representative is visited first number;

In above geometric calibration model, parameter to be calibrated is divided into outer calibration parameter X _EWith the internal calibration parameter X _I, outer calibration parameter

Pitch, roll, yaw are respectively pitching, rolling and yaw direction angle; The internal calibration parameter X _I=(ax ₀, ax ₁, ax ₂, ax ₃, ay ₀, ay ₁, ay ₂, ay ₃), ax ₀, ax ₁, ax ₂, ax ₃, ay ₀, ay ₁, ay ₂, ay ₃For internal calibration is visited the coefficient that unit points to angle model.

3. push away based on spy unit sensing angle optics as claimed in claim 2 and sweep star in rail substep geometric calibration method, it is characterized in that: step 3 comprises following substep,

Step 3.1 is located at and has measured K equally distributed Ground Nuclear Magnetic Resonance reference mark on the image to be calibrated as orientation point, is designated as RP _i, i=1,2,3 ..., K; The object space point that each orientation point is corresponding and picture side's point are designated as respectively RPG _iAnd RPM _i, object space point RPG _iThe WGS84 geocentric rectangular coordinate be (X _i, Y _i, Z _i), as side point RPM _iImage coordinate be (s _i, l _i);

Step 3.2 makes in the formula (1):

$\left\{\begin{array}{c}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]\left(\begin{array}{c}{R}_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}{\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}}\end{array}\right)\\ R_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw})=\left[\begin{array}{c}{a}_1,b_1,c_1\\ a_2,b_2,c_2\\ a_3,b_3,c_3\end{array}\right]\end{array}\right.$

Formula (1) is converted into formula (2):

$\left\{\begin{array}{c}F(X_E,X_I)=\frac{a_1\mathrm{Ux}+b_1\mathrm{Uy}+c_1\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ G(X_E,X_I)=\frac{a_2\mathrm{Ux}+b_2\mathrm{Uy}+c_2\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\end{array}\right.---\left(2\right)$

In the following formula, vector $\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]$ Be object space vector U, vector the coordinate under body coordinate system of representative from the camera projection centre to object space point; a ₁, b ₁, c ₁a ₂, b ₂, c ₂a ₃, b ₃, c ₃Represent respectively camera 9 elements of matrix are installed; F (X _E, X _I), G (X _E, X _I) be respectively along rail and point to angle residual error and vertical rail sensing angle residual error;

Step 3.3 is externally calibrated parameter X E and internal calibration parameter X I initialize

Step 3.4 is with the internal calibration parameter X _ICurrency be considered as true value, with outer calibration parameter X _EBe considered as unknown parameter to be asked, with the internal calibration parameter X _IWith outer calibration parameter X _ECurrency

Substitution formula (2) is to each orientation point RP _i, formula (2) is carried out linearization process, set up error equation (3),

V _i=A _iX-L _i, weigh and be P _i(3)

Wherein $A_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_E}\\ \frac{\∂G_i}{\∂X_E}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂\mathrm{pitch}}\frac{\∂F_i}{\∂\mathrm{roll}}\frac{\∂F_i}{\∂\mathrm{yaw}}\\ \frac{\∂G_i}{\∂\mathrm{pitch}}\frac{\∂G_i}{\∂\mathrm{roll}}\frac{\∂G_i}{\∂\mathrm{yaw}}\end{array}\right]$ $X=d{X}_E=\left[\begin{array}{c}\mathrm{dpitch}\\ \mathrm{droll}\\ \mathrm{dyaw}\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

In the formula, L _iTo utilize inside and outside calibration parameter currency

The error vector that substitution formula (2) calculates; Ai is the matrix of coefficients of error equation; The outer calibration of X representative parameter correction dX _E=(dpitch, droll, dyaw), d represents the correction symbol; P _iCurrent orientation point RP _iPower corresponding to picture point measurement accuracy; F _iAnd G _iBe respectively along rail and point to angle residual error F (X _E, X _I), the rail that hangs down points to angle residual error G (X _E, X _I) function model, obtain corresponding error equation behind the differential;

By formula (4) computing method equation coefficient matrix,

$A^T\mathrm{PA}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{A}_i$

（4）

$A^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{L}_i$

In the following formula, matrix $A=\left[\begin{array}{c}{A}_1\\ ...\\ A_i\\ ...\\ A_K\end{array}\right],$ Matrix Matrix $L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

Utilize least square adjustment to calculate X, suc as formula (5),

X＝(A ^TPA) ^-1(A ^TPL) （5）

Utilize formula (6) to upgrade outer calibration parameter X _ECurrency, then return execution in step 4 iterative computation, enter step 3.4 after the iteration stopping;

$X_E=X_E^o+X---\left(6\right)$

Step 3.5 is resolved the internal calibration parameter, and step 3.4 gained is calibrated parameter X outward _ECurrency be considered as true value, and internal calibration parameter X _IThen be considered as unknown parameter to be asked, with the internal calibration parameter X _IWith outer calibration parameter X _ECurrency Substitution formula (2) is to each orientation point RP _i, formula (2) is carried out linearization process, set up error equation (7),

V _i=B _iY-L _iPower is P _i(7)

Wherein $B_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_I}\\ \frac{\∂G_i}{\∂X_I}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂{\mathrm{ax}}_0}\frac{\∂F_i}{\∂{\mathrm{ax}}_1}\frac{\∂F_i}{\∂{\mathrm{ax}}_2}\frac{\∂F_i}{\∂{\mathrm{ax}}_3}\frac{\∂F_i}{\∂{\mathrm{ay}}_0}\frac{\∂F_i}{\∂{\mathrm{ay}}_1}\frac{\∂F_i}{\∂a{y}_2}\frac{\∂F_i}{\∂{\mathrm{ay}}_3}\\ \frac{\∂G_i}{\∂{\mathrm{ax}}_0}\frac{\∂G_i}{\∂{\mathrm{ax}}_1}\frac{\∂G_i}{\∂{\mathrm{ax}}_2}\frac{\∂G_i}{\∂{\mathrm{ax}}_3}\frac{\∂G_i}{\∂{\mathrm{ay}}_0}\frac{\∂G_i}{\∂{\mathrm{ay}}_1}\frac{\∂G_i}{\∂{\mathrm{ay}}_2}\frac{\∂G_i}{\∂{\mathrm{ay}}_3}\end{array}\right]$ $Y={\mathrm{dX}}_I=\left[\begin{array}{c}{\mathrm{dax}}_0\\ {\mathrm{dax}}_1\\ {\mathrm{dax}}_2\\ {\mathrm{dax}}_3\\ {\mathrm{day}}_0\\ {\mathrm{day}}_1\\ {\mathrm{day}}_2\\ {\mathrm{day}}_3\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

In the formula, L _iTo utilize inside and outside calibration parameter currency

The error vector that substitution formula (2) calculates; Bi is the matrix of coefficients of error equation; Y represents internal calibration parameter correction dX _I, d represents the correction symbol; P _iCurrent orientation point RP _iPower corresponding to picture point measurement accuracy; F _iAnd G _iBe respectively along rail and point to angle residual error F (X _E, X _I), the rail that hangs down points to angle residual error G (X _E, X _I) function model, obtain corresponding error equation behind the differential;

By formula (8) computing method equation coefficient matrix;

$B^T\mathrm{PB}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{B}_i$

（8）

$B^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{L}_i$

In the following formula, $B=\left[\begin{array}{c}{B}_1\\ ...\\ B_i\\ ...\\ B_K\end{array}\right],$

$L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

Utilize least square adjustment to calculate Y, suc as formula (9);

Y＝(B ^TPB) ^-1(B ^TPL) （9）

Utilize formula (10) to upgrade the currency of internal calibration parameter X I, then return execution in step 5 iterative computation, enter step 3.6 after the iteration stopping;

$X_I=X_I^o+Y---\left(10\right)$

Step 3.6 is calibrated the currency XE of parameter and the currency of internal calibration parameter X I outward according to step 3.4 and step 3.5 gained, upgrades the camera parameter file.

Images