A multi-probe star camera internal and external parameter joint calibration method

Abstract

The present application relates to remote sensing satellite ground data processing technical field, especially to a kind of multi-probe star camera internal and external parameter joint calibration method, the method, including constructing initial light direction model, and solving initial attitude parameter;Based on initial attitude parameter, the theoretical image coordinate of star control point is solved, and the relative installation relationship initial value between multi-probe star camera is determined after solving attitude parameter using star point data after gross error elimination;Based on the angular distance constraint between the star points in single frame image, a single-probe internal parameter calibration model and a joint calibration model with internal parameters and relative installation relationship as unknowns are constructed;According to the weighting model, the observation equation is weighted, the internal and external parameters are updated based on the solution, and the star camera accurate internal and external parameters are obtained by iterative optimization until the convergence condition is met.The present application effectively suppresses the influence of observation error on parameter solution, and further improves the calibration accuracy and the reliability of parameter estimation.

Description

Technical Field

[0001] This invention relates to the field of remote sensing satellite ground data processing technology, and in particular to a method for joint calibration of intrinsic and extrinsic parameters of multi-probe satellite cameras. Background Technology

[0002] To achieve high-precision geometric positioning, high-resolution optical remote sensing satellites generally carry multi-probe satellite cameras as attitude measurement payloads. These cameras can acquire optical axis measurement data through simultaneous imaging by multiple probes, calculate the satellite's arcsecond-level attitude parameters in the reference coordinate system, and effectively avoid attitude errors introduced by time synchronization interpolation calculations. This has become the preferred solution for attitude payloads of current high-resolution optical remote sensing satellites.

[0003] The impact forces during satellite launch and the changes in the thermal environment during orbital operation can cause the intrinsic parameters of the satellite camera, such as the lens and focal plane, to deviate from the ground design values, with a maximum deviation of up to 10 pixels. At the same time, the relative installation extrinsic parameters between the probes of the multi-probe satellite camera will also undergo angular grade changes, directly leading to a 5 to 20-fold decrease in the accuracy of satellite attitude calculation, which seriously affects the geometric positioning accuracy of optical satellite imagery. Therefore, it is urgent to perform high-precision calibration of the intrinsic and extrinsic parameters of the on-orbit multi-probe satellite camera to compensate for the accuracy loss caused by parameter deviations. Existing multi-probe satellite camera calibration methods have several technical shortcomings. First, calibrating each satellite camera as an independent payload fails to consider the imaging constraints between multiple probes, making it impossible to achieve collaborative optimization of multi-probe parameters. Second, the step-by-step processing strategy of calibrating internal parameters first and then external parameters relative to the installation results in residuals from internal parameter calibration being directly transferred to the external parameter calibration process, leading to reduced external parameter calibration accuracy and causing inaccurate multi-probe references and decreased fusion attitude accuracy. Third, the gross error removal accuracy is low, and the parameter solution does not consider coupling. Solving multiple parameters at once is prone to divergence, making it difficult to achieve high-precision and stable solutions for camera internal distortion and relative installation parameters, thus failing to meet the high-precision positioning requirements of high-resolution optical remote sensing satellites. At present, a joint calibration method for the internal and external parameters of multi-probe satellite cameras is needed. Summary of the Invention

[0004] To address the technical problems of insufficient consideration of parameter coupling, low accuracy of gross error elimination, and difficulty in calibrating relative installation relationships in existing multi-probe satellite camera calibration, this invention provides a method for joint calibration of internal and external parameters of multi-probe satellite cameras.

[0005] Firstly, the present invention provides a method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera, which adopts the following technical solution:

[0006] A method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera includes:

[0007] An initial ray pointing model is constructed based on the theoretical parameters of the star camera, and the initial attitude parameters are calculated using the image-side coordinates of star points and the celestial coordinates of stars in the sequence star map.

[0008] The theoretical image-side coordinates of the stellar control points are calculated based on the initial attitude parameters, and gross errors are eliminated by comparison and analysis with the coordinates extracted from the star points.

[0009] The attitude parameters were recalculated using the star point data after gross error removal, and the initial values ​​of the relative installation relationship between the multi-probe star cameras were determined.

[0010] Based on the angular distance constraint between star points within a single frame image, a single-probe intrinsic parameter calibration model is constructed.

[0011] Based on the angular distance constraint of the star points simultaneously imaged by multiple probe star cameras, a joint calibration model is constructed with internal parameters and relative installation relationships as unknowns;

[0012] The observation equations are weighted according to the weighting model, and the model parameters are solved by a step-by-step, step-by-step strategy.

[0013] The intrinsic and extrinsic parameters are updated based on the solution results, and the precise intrinsic and extrinsic parameters of the star camera are obtained through iterative optimization until the convergence condition is met.

[0014] Furthermore, the construction of the initial ray pointing model based on the theoretical parameters of the star camera includes, based on the star camera's focal length... Principal point coordinates and theoretical value of pixel size The three-dimensional coordinates of each image point in the star map image in the image space coordinate system are calculated, and the three-dimensional coordinates are converted into ray directions. This is then used to calculate the grid points in the star map image. Based on the ray pointing and its corresponding ray direction, an internal ray pointing model of the star camera is constructed, and the model parameters are solved using the least squares method. The ray pointing model is as follows:

[0015] ,

[0016] in, Let i be the image coordinates of the i-th grid point. and These represent the ray pointing components along the X and Y directions of the grid point in the image space coordinate system. , Let m and n be the parameters to be determined for the ray pointing model, where m and n are non-negative integers and satisfy... .

[0017] Furthermore, the calculation of the initial attitude parameters includes establishing a collinearity condition equation between the ray pointing of star points and the celestial coordinates of stars based on the ray pointing model; constructing a set of error equations using the image-side coordinates of multiple star points in the sequence star image and the right ascension and declination coordinates of their corresponding stars; and obtaining the initial attitude parameters of each star image image through iterative solution. The collinearity condition equation is as follows:

[0018] ,

[0019] in, Let be the rotation matrix of the k-th star image from the J2000 inertial coordinate system to the star camera coordinate system, denoted by the pitch angle. Roll angle and yaw angle Sure, Let J be the right ascension and declination coordinates of the star corresponding to the j-th star point in the k-th star map. These are the image-side coordinates of the star point. For vector normalization operations, , and These are the error components in three directions.

[0020] Furthermore, the gross error removal through comparison analysis with the extracted star point coordinates includes calculating the three-dimensional orientation of the stellar control point in the image space coordinate system based on the initial attitude parameters, solving its theoretical image-side coordinates using the image coordinate inversion model, constructing a two-factor interaction model based on the extracted centroid coordinates and theoretical image-side coordinates of the star point, calculating the theoretical predicted value and its residual of the inversion coordinates, and removing gross errors based on statistical thresholds. The two-factor interaction model is expressed as follows:

[0021] ,

[0022] in, Extract the coordinates of the centroid of the j-th star point. These are the theoretical image-side coordinates obtained by solving the image coordinate inversion model based on the initial attitude parameters. These are the theoretically predicted values ​​of the inversion coordinates calculated based on the two-factor interaction model. , , , and , , , The coefficients to be determined in the two-factor interaction model are obtained by fitting the extracted coordinates and inverted coordinates of the star points.

[0023] Furthermore, determining the initial values ​​of the relative installation relationship between the multi-probe star cameras includes obtaining the attitude parameters calculated by star camera probes A and B at the same imaging time T based on the stellar control coordinates after gross error removal, the star centroid coordinates, and the initial star camera ray pointing model. , Based on the attitude parameters of the image sequence, the rotation matrices of probe A and probe B in the inertial coordinate system can be obtained, and thus the relative positional relationship parameters of probe A and probe B can be obtained. :

[0024] ,

[0025] in, Let be the imaging time of the i-th sequence star map. For rotation matrix to Euler angle function, These are the rotation matrices constructed by probe A and probe B at the i-th imaging moment based on the refined attitude parameters.

[0026] Furthermore, the construction of the single-probe intrinsic parameter calibration model includes selecting any two star points to form a star point pair for a single frame star map image of a single-probe star camera, and calculating the inter-star angular distance in the celestial coordinate system based on the right ascension and declination coordinates of the corresponding stars.

[0027] Based on the image-side coordinates of star points and the ray pointing model of star camera, the ray pointing of star points in the camera coordinate system is calculated and normalized to a unit vector.

[0028] Based on the principle of invariance of inter-satellite angular distance, an error equation is established with the camera's intrinsic parameters as unknowns. A system of intrinsic parameter correction equations is then constructed through linearization. The error equation is expressed as follows:

[0029] ,

[0030] in, Let be the vector angle between star point i and star point j in the celestial coordinate system, calculated from the right ascension and declination coordinates of the corresponding stars. , Let be the normalized unit vectors of the ray pointing from star points i and j in the camera coordinate system, respectively. This is the inter-satellite angular distance error term.

[0031] Furthermore, the construction of the joint calibration model with internal parameters and relative installation relationships as unknowns includes:

[0032] For two probes that image simultaneously in a multi-probe star camera, select corresponding star points located in the star images of the two probes to form star point pairs, and calculate the inter-star angular distance in the celestial coordinate system based on the right ascension and declination coordinates of their corresponding stars.

[0033] Based on the image-side coordinates of the star points and the ray pointing models of each probe, the ray pointing of the star points in their respective camera coordinate systems is calculated and normalized to a unit vector.

[0034] By introducing the relative installation relationship between the two probes, a joint error equation is established with the intrinsic parameters of the two probes and the relative installation relationship as unknowns. A joint calibration equation system is then constructed through linearization. The expression for the joint error equation is as follows:

[0035] ,

[0036] in, The vector angle between star point I in the star image of probe A and star point J in the star image of probe B, in the celestial coordinate system, is calculated from the right ascension and declination coordinates of the corresponding stars. , Let I and J be the unit vectors pointing towards the light rays from star points I and J, respectively. Let the rotation matrix represent the relative mounting relationship between probe A and probe B, expressed by Euler angles. Sure, This is the inter-satellite angular distance error term.

[0037] Furthermore, the weighting of the observation equations according to the weighting model includes determining the variance of the dot product observations based on the inter-star angle equation and the error propagation theorem; setting the weights of each observation equation according to the size of the inter-star angle based on the relationship between the variance and the inter-star angle; and constructing a weight matrix. The formula for calculating the weights of the observation equations is as follows:

[0038] ,

[0039] in, The variance of the dot product observations, The standard deviation of the error. The actual angle between the stars. These are the weights for the corresponding observation equations.

[0040] Furthermore, the step-by-step, step-by-step strategy includes grouping the intrinsic parameters of probe A, the intrinsic parameters of probe B, and the relative installation parameters between the probes. By fixing the intrinsic parameters of probe A and probe B, the relative installation parameter correction between probe A and probe B is solved using weighted least squares. Then, by fixing the relative installation parameters, the lower-order and higher-order terms of the intrinsic parameters of probe A and probe B are solved sequentially. The formula for calculating the relative installation parameter correction is:

[0041] ,

[0042] in, For the relative installation parameter correction vector, , and These are the corrections for the relative Euler angles between probe A and probe B. This is a coefficient matrix of relative installation parameters. This is the corresponding weight matrix.

[0043] Furthermore, obtaining the precise intrinsic and extrinsic parameters of the star camera includes superimposing the calculated intrinsic parameter corrections and relative installation parameter corrections onto the current values ​​of the corresponding parameters, updating the star camera's intrinsic parameters and relative installation relationships, calculating the residuals of star point image-side ray pointing and inter-star angles based on the updated parameters, calculating the standard error, using a multiple of the standard error as a threshold for gross error removal, and then determining the convergence condition of the star point residuals. The expression for the standard error is:

[0044] ,

[0045] in, This represents the mean square error of the inter-satellite angle residual. Let be the angle residual between star point i and star point j, and n be the number of star points in the current star map.

[0046] In summary, the present invention has the following beneficial technical effects:

[0047] 1. This invention constructs a unified joint calibration model for the internal and external parameters of a multi-probe satellite camera, achieving coordinated optimization of internal parameters, attitude parameters, and the relative installation relationship of multiple probes. This overcomes the technical defects of traditional methods that calibrate internal and external parameters separately, resulting in insufficient consideration of parameter coupling and limited calibration accuracy. It significantly improves the overall calibration accuracy and parameter consistency of the multi-probe satellite camera.

[0048] 2. This invention employs a two-factor interactive model for gross error removal and combines iterative optimization with angular distance constraints to achieve dual detection and dynamic removal of gross errors. Compared with the traditional single threshold method, it effectively improves the gross error identification rate, avoids interference from abnormal observations on parameter calculation, and significantly enhances the reliability of calibration results and the robustness of data processing.

[0049] 3. This invention utilizes the principle of inter-satellite angular distance invariance to construct a calibration model, enabling direct calculation of camera intrinsic parameters and relative installation relationships without relying on external high-precision attitude reference data. This reduces the dependence of the calibration process on external conditions. Furthermore, through a step-by-step solution strategy, it achieves step-by-step decoupling and progressive refinement of parameters, effectively solving the problem of easy divergence in simultaneous multi-parameter solutions and significantly improving the numerical stability of parameter solutions.

[0050] 4. This invention achieves optimal utilization of observation information by setting the weights of the observation equations according to the weighting model and performing adaptive weighting based on the size of the star point angle. This allows precise observations with smaller angles to receive larger weights, while noisy observations with larger angles receive smaller weights. This effectively suppresses the impact of observation errors on parameter calculation and further improves calibration accuracy and the reliability of parameter estimation.

[0051] 5. This invention achieves a gradual improvement in calibration accuracy and progressive optimization of intrinsic and extrinsic parameters by iteratively performing parameter calculation and gross error elimination until the convergence condition is met. Ultimately, it brings the star point residuals to sub-pixel level accuracy, providing effective technical support for high-precision geometric positioning of high-resolution optical remote sensing satellite images. It has significant application value and engineering practicality. Attached Figure Description

[0052] Figure 1 This is a schematic diagram of the overall process of a method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to an embodiment of the present invention.

[0053] Figure 2 This is an architecture diagram of a method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to an embodiment of the present invention.

[0054] Figure 3 This is an implementation effect diagram of a method for joint calibration of internal and external parameters of a multi-probe satellite camera according to an embodiment of the present invention. Detailed Implementation

[0055] The present invention will be further described in detail below with reference to the accompanying drawings.

[0056] Example 1: Refer to Figure 1 This embodiment of a method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera includes:

[0057] S1. Construct an initial ray pointing model based on the theoretical parameters of the star camera, and calculate the initial attitude parameters using the image-side coordinates of star points and the celestial coordinates of stars in the sequence star map;

[0058] S2. Calculate the theoretical image-side coordinates of the stellar control points based on the initial attitude parameters, and perform gross error elimination by comparing and analyzing the coordinates extracted from the star points;

[0059] S3. Recalculate the attitude parameters using the star point data after gross error removal, and determine the initial values ​​of the relative installation relationship between the multi-probe star cameras;

[0060] S4. Based on the angular distance constraints between star points within a single frame image, construct a single-probe intrinsic parameter calibration model;

[0061] S5. Based on the angular distance constraint of the star points simultaneously imaged by the multi-probe star cameras, a joint calibration model is constructed with the internal parameters and relative installation relationship as unknowns;

[0062] S6. Weight the observation equations according to the weighting model, and solve the model parameters using a step-by-step, step-by-step strategy;

[0063] S7. Update the intrinsic and extrinsic parameters based on the solution results, and iterate until the convergence condition is met to obtain the accurate intrinsic and extrinsic parameters of the star camera.

[0064] Specifically, a method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera includes the following:

[0065] like Figure 1 , Figure 2 As shown, step S1 aims to establish the initial imaging geometric model of the star camera and provide initial attitude values ​​for joint calibration. Specifically, it includes two steps: constructing the ray pointing model and calculating the initial attitude parameters.

[0066] In the ray pointing model construction stage, the three-dimensional coordinates in the image space coordinate system are first calculated based on the theoretical parameters of the star camera, and then based on the star camera's focal length f and the coordinates of the principal point. Theoretical value of pixel size For any point in the star map image, based on its image coordinates Its three-dimensional coordinates in the image space coordinate system can be calculated:

[0067] ,

[0068] Where X, Y, and Z are the coordinate components of the point along the three coordinate axes in the image space coordinate system, respectively. The negative sign indicates that the image point is located behind the focal plane. The above three-dimensional coordinates are then converted into ray directions. The ray direction describes the direction of a ray originating from the projection center and passing through the image point in the image space coordinate system. Its calculation formula is as follows:

[0069] ,

[0070] in, , These are the X and Y components of the ray direction, respectively. This transformation eliminates the influence of focal length, making the ray direction a dimensionless quantity related to the image point position, utilizing the grid points in the star map image. and the direction of its light A ray pointing model for the internal workings of the star camera was constructed. This model uses a cubic polynomial to describe the nonlinear mapping between image point coordinates and ray pointing, in order to compensate for lens distortion and focal plane distortion.

[0071] ,

[0072] in, Let i be the image coordinates of the i-th grid point. and These represent the ray pointing components along the X and Y directions of the grid point in the image space coordinate system. , Let m and n be the parameters to be determined for the ray pointing model, where m and n are non-negative integers and satisfy... This constraint ensures that the highest order of the polynomial is cubic, and includes 10 parameters: constant term, linear term, quadratic term, and cubic term. Specifically, it can be broken down as follows:

[0073] ,

[0074] The least squares method is used to solve for the above model parameters. By minimizing the sum of squared residuals between the observed values ​​of the ray pointing at the grid points and the calculated values ​​of the model, the optimal parameter estimate is obtained, thereby establishing the initial ray pointing model of the star camera.

[0075] In the initial attitude parameter calculation stage, based on the aforementioned ray pointing model, a collinearity condition equation is established between the ray pointing of a star point and the celestial coordinates of the star. For the j-th star point in the k-th star chart, its ray pointing and the celestial coordinates of the star should satisfy the collinearity condition:

[0076] ,

[0077] in, Let be the rotation matrix of the k-th star image from the J2000 inertial coordinate system to the star camera coordinate system, denoted by the pitch angle. Roll angle and yaw angle Sure, The right ascension and declination coordinates of the star corresponding to the j-th star point in the k-th star map are obtained from star catalog data. The image-side coordinates of the star point are obtained from the star image using a star centroid extraction algorithm. To represent the vector normalization operation, , and Taking the error components in three directions as examples, the above error equation is linearized to obtain the error equation concerning the attitude parameter corrections:

[0078] ,

[0079] in, , , The partial derivatives of the error components with respect to the attitude parameters form the coefficient matrix of the error equation. , and These are the correction values ​​for the corresponding attitude parameters. , and This is the initial value of the error component calculated using the approximate values ​​of the current attitude parameters.

[0080] The above error equation system is constructed using M star points from the k-th star image, forming an overdetermined system of 3M equations and 3 unknowns. Attitude parameter corrections are obtained through iterative least squares solving, and the approximate attitude parameter values ​​are updated. This iteration is repeated until the corrections are less than a threshold, thus obtaining the initial attitude parameters for the k-th star image. .

[0081] S2. Geometric inversion of stellar control points is performed using the initial attitude parameters, and gross errors in star point extraction are identified and eliminated through a two-factor interaction model. This includes three steps: theoretical image-space coordinate inversion, two-factor interaction model construction, and statistical threshold gross error elimination. In the theoretical image-space coordinate inversion step, based on the initial attitude parameters calculated in step 1, the three-dimensional orientation of the stellar control points in the image space coordinate system is calculated. For the j-th stellar control point in the k-th star image, its three-dimensional orientation in the image space coordinate system is calculated by the following formula:

[0082] ,

[0083] in, , and These represent the pointing components of the star's control point along the three coordinate axes in the image space coordinate system. The rotation matrix is ​​constructed based on the initial attitude parameters. Here are the right ascension and declination coordinates of the star corresponding to the control point of that star. Using the image coordinate inversion model, the above image space orientation is converted into theoretical image-side coordinates. According to the principle of central projection imaging, the image point coordinates and the image space orientation satisfy a proportional relationship. Therefore, the inversion formula for the theoretical image-side coordinates can be obtained:

[0084] , ,

[0085] in, , Let be the theoretical image row and column coordinates of the j-th star control point, respectively. Focal length For pixel size, and For the principal point coordinates, this inversion process establishes a mapping from celestial coordinates to image coordinates.

[0086] In the two-factor interaction model construction phase, a two-factor interaction model is constructed based on the centroid-extracted coordinates of the star points and the theoretically inverted coordinates to describe the systematic deviation between the two. (Star point centroid-extracted coordinates) The coordinates are obtained from star map images using image processing algorithms, and theoretically inverted coordinates are used. The geometric inversion calculations described above show that, due to the existence of initial attitude error and extraction error, there is a nonlinear deviation between the two. This invention uses a two-factor interaction model to model this deviation.

[0087] ,

[0088] in, Extract the coordinates of the centroid of the j-th star point. These are the theoretical image-side coordinates obtained by solving the image coordinate inversion model based on the initial attitude parameters. These are the theoretically predicted values ​​of the inversion coordinates calculated based on the two-factor interaction model. , , , and , , , The coefficients to be determined for the two-factor interaction model are as follows: , For constant terms, , and , The coefficients of the linear term, , These are the interaction term coefficients, used to capture the coupling effects between row and column coordinates. Using the extracted and inverted coordinates of the M star points in the k-th star image, the above 8 model coefficients are solved through least-squares fitting to establish the mapping relationship between the extracted and inverted coordinates.

[0089] In the statistical threshold outlier removal stage, the residuals of each star point are calculated based on a two-factor interactive model, and statistical thresholds are used to identify outliers. For the j-th star control point, its actual inversion coordinates are compared with the model predictions to calculate the residuals. :

[0090] , ,

[0091] in, , Given the residuals in the row and column directions respectively, calculate the statistical standard deviation of all star point residuals:

[0092] , ,

[0093] in, , These are the standard deviations of the row and column residuals, respectively, with 3... 3 As a threshold, residuals with an absolute value greater than 3 are considered. Or 3 The star points were identified as gross errors and discarded. The criteria are based on the normal distribution assumption and can effectively identify outliers that deviate significantly from the normal distribution. By combining the above two-factor interaction model with statistical thresholds, dual detection and removal of gross errors in star point extraction are achieved, resulting in high-quality star point data after removing gross errors.

[0094] S3. The attitude parameters are refined again using the high-quality star point data after gross error removal, and the initial values ​​of the relative installation relationships between the multi-probe star cameras are calculated. In the attitude parameter recalculation stage, based on the star control coordinates, star centroid coordinates, and the initial star camera ray pointing model established in step 1 after gross error removal, the attitude parameters of the sequence of star images are recalculated. This process uses the same collinearity condition equation and iterative least squares solution method as in step 1, but because gross errors have been removed from the input data, the accuracy of the obtained attitude parameters is significantly improved. For the same imaging time T, attitude calculations are performed on the star imagery of star camera probes A and B respectively to obtain the attitude parameters of probe A. and the attitude parameters of probe B This synchronous imaging feature avoids attitude errors introduced by time interpolation.

[0095] In the initial value determination stage of the relative installation relationship, based on the attitude parameter calculation results of the image sequence, the relative installation relationship between probe A and probe B is determined. According to the rigid body attitude transformation principle, if both probes observe the same inertial space at the same time, there is a definite relative transformation relationship between their attitude parameters. The attitude parameters of probe A and probe B are converted into rotation matrix form. The rotation matrix of probe A from the inertial coordinate system to the home system is denoted as... The rotation matrix of probe B from the inertial coordinate system to the home system is denoted as... Then the relative rotation matrix of probe B relative to probe A satisfy:

[0096] ,

[0097] This equation shows that the relative rotation matrix can be obtained by multiplying the rotation matrices of the two probes, where This is the transpose of the rotation matrix of probe A, i.e., its inverse matrix. For sequential star map images, multiple imaging times are acquired. The relative rotation matrix below The estimation accuracy of relative installation relationships is improved by averaging. The relative rotation matrix at each time point is converted into Euler angles, and then the Euler angles are arithmetically averaged to obtain the final relative positional relationship parameters. :

[0098] ,

[0099] in, Let be the imaging time of the i-th sequence star map. This is the function to convert a rotation matrix to Euler angles, which transforms a 3×3 rotation matrix into three Euler angles. , representing the rotation angles about different axes, , Let be the rotation matrices constructed by probe A and probe B at the i-th imaging moment based on refined attitude parameters, and M be the number of sequence star maps involved in the calculation. This averaging process effectively suppresses the influence of random errors, resulting in better stability of the initial values ​​of the relative installation relationship.

[0100] S4. Utilizing the principle of invariance of interstellar angular distance, a calibration model is constructed with only the camera's intrinsic parameters as unknowns, achieving effective separation of intrinsic and attitude parameters. In the interstellar angular distance calculation stage, for a single frame of star imagery from a single-probe star camera, any two star points are selected to form a star point pair, and the interstellar angular distance in the celestial coordinate system is calculated based on the right ascension and declination coordinates of their corresponding stars. For star point i and star point j, the right ascension and declination coordinates of their corresponding stars are respectively... and Then the vector angle between the two star points in the celestial coordinate system satisfy:

[0101] ,

[0102] This formula uses the principle of spherical trigonometry to convert the right ascension and declination coordinates of stars into unit vectors in the celestial coordinate system. Then, it calculates the dot product of the two vectors to obtain the cosine of the included angle, and the interstellar angular distance. It is an inherent quantity determined by the celestial coordinates of a star, does not change with the observation attitude, has angular invariance, and can be used as a reference constraint for calibration.

[0103] In the ray pointing normalization step, based on the image-side coordinates of the star point and the star camera ray pointing model, the ray pointing of the star point in the camera coordinate system is calculated and normalized to a unit vector. For star point i, its image-side coordinates are... Calculate the ray direction based on the current ray direction model parameters:

[0104] ,

[0105] in, , These are the ray pointing components of star point i along the X and Y directions in the camera coordinate system of probe A, respectively, calculated from the ray pointing model. The ray pointing is then normalized to a unit vector:

[0106] ,

[0107] in, As the normalization factor, This is the normalized unit vector with a magnitude of 1, representing only the direction of light rays. Similarly, performing the same operation on star point j yields the normalized unit vector. .

[0108] In the process of constructing the intrinsic parameter error equation, based on the principle of invariance of inter-satellite angular distance, an error equation is established with the camera's intrinsic parameters as unknowns. The invariance of inter-satellite angular distance indicates that the inter-satellite angle in the celestial coordinate system should be equal to the angle between the rays pointing from two stars in the camera coordinate system. Therefore, the error equation is established:

[0109] ,

[0110] in, Let be the vector angle between star point i and star point j in the celestial coordinate system, calculated from the right ascension and declination coordinates of the corresponding stars. , Let be the normalized unit vectors of the ray pointing from star points i and j in the camera coordinate system, respectively. The inter-satellite angular distance error term should ideally be zero. Using the ray pointing model parameters of the star camera probe A as unknowns, the above error equation is linearized to obtain the intrinsic parameter correction equation:

[0111] ,

[0112] in, These are the corrections for the intrinsic parameters of probe A, with a total of 20 unknowns. , The partial derivatives of the error term with respect to each intrinsic parameter constitute the coefficients of the error equation. This is the initial value for the error term calculated using the current approximation of the intrinsic parameters. Using star points on single frames of the star map sequence from probe A, pairwise combinations of star points are performed to obtain a large number of inter-star angle observations, constructing an overdetermined set of intrinsic parameter correction equations. Similarly, for probe B, based on the star point angles in each frame of its star map sequence, the same principle is used to construct the intrinsic parameter correction equations for probe B. This model only uses the camera's intrinsic parameters as unknowns and does not involve attitude parameters, achieving effective separation between intrinsic parameters and external attitude.

[0113] S5. Utilizing the angular distance constraint of star points simultaneously imaged by multiple probe star cameras, a joint calibration model is constructed with internal parameters and relative installation relationships as unknowns. In the cross-probe star point pair selection stage, for two probes simultaneously imaging in the multi-probe star camera system, corresponding star points located in the star imagery of each probe are selected to form star point pairs. Since probe A and probe B have simultaneous imaging capabilities, at the same imaging moment, both probes acquire local star images of the same night sky, which inevitably involves joint observations of the same star or cross-observations of adjacent stars. The star image from probe A is selected. Star point I and probe B star map The star points J in the image form a cross-probe star point pair. It is required that the celestial coordinates of the two star points are known and the included angle is appropriate to ensure the geometric strength of the angular distance observation.

[0114] In the interstellar angular distance calculation stage, the interstellar angular distance in the celestial coordinate system is calculated based on the right ascension and declination coordinates of the corresponding stars for the cross-probe star points. For star point I and star point J, the right ascension and declination coordinates of their corresponding stars are respectively... and Then the vector angle between the two star points in the celestial coordinate system satisfy:

[0115] ,

[0116] The inter-satellite angular distance is used as a reference true value to constrain the geometric relationship between the two probes.

[0117] In the relative installation relationship introduction stage, based on the image-side coordinates of the star points and the ray pointing models of each probe, the ray pointing of the star points in their respective camera coordinate systems is calculated and normalized to a unit vector. For star point I in probe A, its image-side coordinates are... The ray direction of probe A is calculated and normalized based on the ray direction model of probe A:

[0118] ,

[0119] in, Let I be the unit vector of star point I in the camera coordinate system of probe A. Similarly, for star point J in probe B, its image-side coordinates are... The normalized unit vector is calculated based on the ray pointing model of probe B:

[0120] ,

[0121] Because the two probes are installed relative to each other, the unit vector in the camera coordinate system of probe B needs to be transformed to the camera coordinate system of probe A by a relative rotation matrix before it can be used in conjunction with the probe. Angle comparisons are performed. The relative installation relationship between probe A and probe B is introduced. Euler angles Constructing the relative rotation matrix This matrix describes the attitude transformation of the camera coordinate system of probe B relative to the camera coordinate system of probe A.

[0122] In the joint error equation construction stage, a joint error equation is established with the intrinsic parameters of the two probes and their relative installation relationship as unknowns. Based on the principle of invariance of inter-satellite angular distance, the inter-satellite angle in the celestial coordinate system should be equal to the angle between the rays pointing from the two stars after the relative installation relationship transformation. Therefore, the joint error equation is established:

[0123] ,

[0124] in, The vector angle between star point I in the star image of probe A and star point J in the star image of probe B, in the celestial coordinate system, is calculated from the right ascension and declination coordinates of the corresponding stars. , Let I and J be the unit vectors pointing towards the light rays from star points I and J, respectively. Let the rotation matrix represent the relative mounting relationship between probe A and probe B, expressed by Euler angles. Sure, The equation for the inter-satellite angular distance error term includes the ray pointing model parameters of probe A, the ray pointing model parameters of probe B, and the equation for the relative installation parameters between probe A and probe B. Linearizing the above equation yields:

[0125] ,

[0126] in, Here are the corrections to the intrinsic parameters of probe A. This is the correction value for the intrinsic parameters of probe B. This represents the correction for the relative installation parameters between probe A and probe B. The partial derivative terms constitute the coefficient matrix of the error equation. To calculate the initial value of the error term using the approximate values ​​of the current parameters, star points on the star map images simultaneously imaged by probes A and B are combined pairwise to obtain the inter-star angles. This allows for the construction of a system of equations for the corrections to the intrinsic parameters of probe A, the intrinsic parameters of probe B, and the relative installation parameters between probes A and B. For star map images imaged sequentially and simultaneously, a system of equations for simultaneous imaging of star map sequences can be constructed.

[0127] S6. In the weighting model construction stage, based on the inter-satellite angle equation and the error propagation theorem, the variance of the dot product observations is determined, and the weights of each observation equation are set based on the relationship between the variance and the inter-satellite angle. For inter-satellite angular distance observations, the dot product observations are... The error propagation characteristics are closely related to the size of the included angle. According to the error propagation theorem, when the angle observation has a standard deviation of... When the random error is , the variance of the dot product observations is:

[0128] ,

[0129] in, The variance of the dot product observations, The standard deviation of the angle error. Let be the actual angle between the star points; this formula shows that when the angle between the star points is small, When the variance approaches zero, the variance of the dot product observations is small, resulting in higher observation accuracy; when the included angle is close to 90 degrees... The variance reaches its maximum when the angle approaches 1; it decreases again when the angle approaches 180 degrees. To balance the differences in accuracy between observations at different angles, the weights of each observation equation are set according to the size of the star point angle, and the weights are inversely proportional to the square root of the variance.

[0130] ,

[0131] in, The weights of the observation equations formed by star points i and j are determined by a weighting strategy that assigns greater weights to precise observations with smaller included angles and less weights to noisy observations with larger included angles, effectively suppressing the impact of observation errors on parameter calculation. A weight matrix is ​​constructed based on the weights of each observation equation. This matrix is ​​a diagonal matrix, and its diagonal elements consist of the weight values ​​of each pair of star points, which are used for subsequent weighted least squares parameter calculation.

[0132] In the step-by-step parameter calculation process, the intrinsic parameters of probe A, the intrinsic parameters of probe B, and the relative installation parameters between probes are grouped. By fixing the known parameter groups and releasing the parameter groups to be calculated step by step, various parameters are solved sequentially, avoiding strong coupling and divergence problems caused by solving multiple parameters simultaneously. The specific calculation process consists of four steps: First, the current intrinsic parameters of probe A and probe B are fixed as known quantities. The relative installation parameter correction between probe A and probe B is solved using weighted least squares, and the formula is as follows:

[0133] ,

[0134] in, For the relative installation parameter correction vector, , and These are the corrections for the relative Euler angles between probe A and probe B. The coefficient matrix for relative installation parameters is composed of the partial derivatives of the joint error equation with respect to the relative installation parameters. Given the corresponding weight matrix, the obtained relative installation parameters are fixed as known quantities. Then, the lower-order terms of the intrinsic parameters of probe A and probe B are solved using weighted least squares, i.e., the coefficients of the first-order terms of the intrinsic parameters of probe A. Probe B internal parameters The solution formula is:

[0135] ,

[0136] in, This is the coefficient matrix of the low-order intrinsic parameters. Next, with the relative installation parameters and low-order intrinsic parameters fixed as known quantities, the second-order terms of the intrinsic parameters of probe A and probe B are solved, i.e. and The solution formula is:

[0137] ,

[0138] in, Given the coefficient matrix of the second-order intrinsic parameters, and finally, fixing the relative installation parameters, lower-order and second-order intrinsic parameters as known quantities, solve for the third-order terms of the intrinsic parameters of probe A and probe B, i.e., the intrinsic parameters of probe A. and probe B intrinsic parameters The solution formula is:

[0139] ,

[0140] in, Let be the coefficient matrix of the third-order intrinsic parameters. This step-by-step, order-wise strategy effectively avoids strong coupling interference between higher-order and lower-order parameters by releasing the degrees of freedom of the parameters in stages. It ensures that the lower-order parameters first form a stable solution benchmark, providing reliable initial values ​​for solving the higher-order parameters, and finally achieving high-precision solution of the corrections for all-order parameters.

[0141] S7. Update the intrinsic and extrinsic parameters based on the solution results, and iterate until the convergence condition is met to obtain the accurate intrinsic and extrinsic parameters of the star camera.

[0142] Using the calibrated intrinsic and extrinsic parameters of the star camera and its relative installation relationship, gross errors are eliminated according to angular distance. Steps 4, 5, and 6 are iteratively executed until the error term is less than the threshold, thus obtaining the final accurate intrinsic and extrinsic parameters of the star camera. Using the solution result from step 6, the image-side ray pointing result of the current star point is calculated.

[0143] ,

[0144] Calculate the residual angle between the two star points:

[0145] ,

[0146] in, This represents the direction of the ray calculated from star point i using the current camera parameters. This indicates the direction of the ray calculated from star point j using the current camera parameters. This is the interstellar angle calculated from the right ascension coordinates of a star.

[0147] Given a star chart with n pairs of points, the error is:

[0148] ,

[0149] by As a threshold, outlier removal is performed on star points. Further, steps 4, 5, and 6 are repeated until the star point residuals are all less than [a certain value]. Threshold or less The iteration is terminated, and the final calibration result is output.

[0150] Example 2: The difference between this example and Example 1 is that this example provides a simulation experiment verification of a method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera;

[0151] To verify the effectiveness of the proposed method for joint calibration of intrinsic and extrinsic parameters of multi-probe satellite cameras, measured image data from satellite cameras A and B on an in-orbit satellite were used to conduct a comparative experiment with the traditional satellite camera calibration method based on a rigorous imaging model.

[0152] A multi-probe satellite camera system aboard a high-resolution optical remote sensing satellite was selected as the experimental object. This system consists of probe A and probe B and has simultaneous imaging capability. Using the sequence of satellite image data acquired in orbit, parameter calculations were performed using both the traditional rigorous imaging model calibration method and the method described in this embodiment. The traditional method employs a step-by-step calibration strategy for intrinsic and extrinsic parameters. First, the intrinsic parameters of each probe are calibrated independently, and then the relative installation relationships are calculated based on the calibrated intrinsic parameters.

[0153] Image-square residuals were used as the evaluation index for calibration accuracy. Image-square residuals are defined as the deviation between the extracted coordinates of the star centroid and the theoretical image-square coordinates retrieved based on the calibration parameters. A smaller image-square residual indicates that the calibration parameters more accurately describe the imaging geometry of the star camera, resulting in higher calibration accuracy. Experimental results are shown below. Figure 3 As shown in the figure, the horizontal axis represents the star point number, and the vertical axis represents the image-side residual value (unit: pixel). Comparative analysis shows that the image-side residual distribution of the traditional rigorous imaging model method is relatively dispersed, with many residual points greater than 0.6 pixels, and the maximum residual is close to 1 pixel; while the image-side residual of the method of this invention is significantly reduced, with most star point residuals less than 0.6 pixels, the residual distribution is more concentrated, and no obvious outliers are found.

[0154] Quantitative statistical results show that the image-side residual mean error of the traditional method is 0.3 pixels, while that of the method of this invention is 0.2 pixels. This result indicates that this embodiment effectively suppresses the influence of parameter coupling error and observation outliers by constructing a joint optimization model of internal and external parameters, using a two-factor interactive model for gross error elimination, and implementing a step-by-step solution strategy. This significantly improves the calibration accuracy of multi-probe satellite cameras and provides reliable technical support for high-precision geometric positioning of high-resolution optical remote sensing satellites.

[0155] The above are all preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Therefore, all equivalent changes made in accordance with the structure, shape and principle of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera, characterized in that, include: An initial ray pointing model is constructed based on the theoretical parameters of the star camera, and the initial attitude parameters are calculated using the image-side coordinates of star points and the celestial coordinates of stars in the sequence star map. The theoretical image-side coordinates of the stellar control points are calculated based on the initial attitude parameters, and gross errors are eliminated by comparison and analysis with the coordinates extracted from the star points. The attitude parameters were recalculated using the star point data after gross error removal, and the initial values ​​of the relative installation relationship between the multi-probe star cameras were determined. Based on the angular distance constraint between star points within a single frame image, a single-probe intrinsic parameter calibration model is constructed. Based on the angular distance constraint of the star points simultaneously imaged by multiple probe star cameras, a joint calibration model is constructed with internal parameters and relative installation relationships as unknowns; The observation equations are weighted according to the weighting model, and the model parameters are solved by a step-by-step, step-by-step strategy. The intrinsic and extrinsic parameters are updated based on the solution results, and the precise intrinsic and extrinsic parameters of the star camera are obtained through iterative optimization until the convergence condition is met.

2. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The initial ray pointing model is constructed based on the theoretical parameters of the star camera, including the focal length of the star camera. Principal point coordinates and theoretical value of pixel size The three-dimensional coordinates of each image point in the star map image in the image space coordinate system are calculated, and the three-dimensional coordinates are converted into ray directions. This is then used to calculate the grid points in the star map image. Based on the ray pointing and its corresponding ray direction, an internal ray pointing model of the star camera is constructed, and the model parameters are solved using the least squares method. The ray pointing model is as follows: ， in, Let i be the image coordinates of the i-th grid point. and These represent the ray pointing components along the X and Y directions of the grid point in the image space coordinate system. , Let m and n be the parameters to be determined for the ray pointing model, where m and n are non-negative integers and satisfy... .

3. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 2, characterized in that, The calculation of the initial attitude parameters includes establishing a collinearity condition equation between the ray pointing of star points and the celestial coordinates of stars based on the ray pointing model; constructing a set of error equations using the image-side coordinates of multiple star points in the sequence star image and the right ascension and declination coordinates of their corresponding stars; and obtaining the initial attitude parameters of each star image image through iterative solution. The collinearity condition equation is as follows: ， in, Let be the rotation matrix of the k-th star image from the J2000 inertial coordinate system to the star camera coordinate system, denoted by the pitch angle. Roll angle and yaw angle Sure, Let J be the right ascension and declination coordinates of the star corresponding to the j-th star point in the k-th star map. These are the image-side coordinates of the star point. This represents the vector normalization operation. , and These are the error components in three directions.

4. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The gross error removal process, which involves comparing the star point's extracted coordinates with the external coordinates, includes: calculating the three-dimensional orientation of the star control point in the image space coordinate system based on initial attitude parameters; solving for its theoretical image-side coordinates using an image coordinate inversion model; constructing a two-factor interactive model based on the extracted centroid coordinates and theoretical image-side coordinates of the star points; calculating the theoretical predicted values ​​and residuals of the inverted coordinates; and removing gross errors based on statistical thresholds. The two-factor interactive model is expressed as follows: ， in, Extract the coordinates of the centroid of the j-th star point. These are the theoretical image-side coordinates obtained by solving the image coordinate inversion model based on the initial attitude parameters. These are the theoretically predicted values ​​of the inversion coordinates calculated based on the two-factor interaction model. , , , and , , , The coefficients to be determined in the two-factor interaction model are obtained by fitting the extracted coordinates and inverted coordinates of the star points.

5. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The determination of the initial values ​​for the relative installation relationship between the multi-probe star cameras includes obtaining the attitude parameters calculated by star camera probe A and probe B at the same imaging time T based on the stellar control coordinates after gross error removal, the star centroid coordinates, and the initial star camera ray pointing model. , Based on the attitude parameters of the image sequence, the rotation matrices of probe A and probe B in the inertial coordinate system can be obtained, and thus the relative positional relationship parameters of probe A and probe B can be obtained. : ， in, Let be the imaging time of the i-th sequence star map. For rotation matrix to Euler angle function, These are the rotation matrices constructed by probe A and probe B at the i-th imaging moment based on the refined attitude parameters.

6. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The construction of the single-probe intrinsic parameter calibration model includes selecting any two star points to form a star point pair for a single frame star map image of a single-probe star camera, and calculating the inter-star angular distance in the celestial coordinate system based on the right ascension and declination coordinates of the corresponding stars. Based on the image-side coordinates of star points and the ray pointing model of star camera, the ray pointing of star points in the camera coordinate system is calculated and normalized to a unit vector. Based on the principle of invariance of inter-satellite angular distance, an error equation is established with the camera's intrinsic parameters as unknowns. A system of intrinsic parameter correction equations is then constructed through linearization. The error equation is expressed as follows: ， in, Let be the vector angle between star point i and star point j in the celestial coordinate system, calculated from the right ascension and declination coordinates of the corresponding stars. , Let be the normalized unit vectors of the ray pointing from star points i and j in the camera coordinate system, respectively. This is the inter-satellite angular distance error term.

7. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The construction of the joint calibration model with internal parameters and relative installation relationships as unknowns includes: For two probes that image simultaneously in a multi-probe star camera, select corresponding star points located in the star images of the two probes to form star point pairs, and calculate the inter-star angular distance in the celestial coordinate system based on the right ascension and declination coordinates of their corresponding stars. Based on the image-side coordinates of the star points and the ray pointing models of each probe, the ray pointing of the star points in their respective camera coordinate systems is calculated and normalized to a unit vector. By introducing the relative installation relationship between the two probes, a joint error equation is established with the intrinsic parameters of the two probes and the relative installation relationship as unknowns. A joint calibration equation system is then constructed through linearization. The expression for the joint error equation is as follows: ， in, The vector angle between star point I in the star image of probe A and star point J in the star image of probe B, in the celestial coordinate system, is calculated from the right ascension and declination coordinates of the corresponding stars. , Let I and J be the unit vectors pointing towards the light rays from star points I and J, respectively. Let the rotation matrix represent the relative mounting relationship between probe A and probe B, expressed by Euler angles. Sure, This is the inter-satellite angular distance error term.

8. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The weighting of the observation equations according to the weighting model includes determining the variance of the dot product observations based on the inter-star angle equation and the error propagation theorem; setting the weights of each observation equation according to the size of the inter-star angle based on the relationship between the variance and the inter-star angle; and constructing a weight matrix. The formula for calculating the weights of the observation equations is as follows: ， in, The variance of the dot product observations, The standard deviation of the error. The actual angle between the stars. These are the weights for the corresponding observation equations.

9. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 8, characterized in that, The step-by-step, step-by-step strategy includes grouping the intrinsic parameters of probe A, the intrinsic parameters of probe B, and the relative installation parameters between the probes. By fixing the intrinsic parameters of probe A and probe B, the relative installation parameter correction between probe A and probe B is solved using weighted least squares. Then, by fixing the relative installation parameters, the lower-order and higher-order terms of the intrinsic parameters of probe A and probe B are solved sequentially. The formula for calculating the relative installation parameter correction is: ， in, For the relative installation parameter correction vector, , and These are the corrections for the relative Euler angles between probe A and probe B. This is a coefficient matrix of relative installation parameters. This is the corresponding weight matrix.

10. The method for joint calibration of intrinsic and extrinsic parameters of a multi-probe satellite camera according to claim 1, characterized in that, The process of obtaining precise intrinsic and extrinsic parameters of the star camera includes superimposing the calculated intrinsic parameter corrections and relative installation parameter corrections onto the current values ​​of the corresponding parameters, updating the star camera's intrinsic parameters and relative installation relationships, calculating the residuals of star point image-side ray pointing and inter-star angles based on the updated parameters, calculating the standard error, using a multiple of the standard error as a threshold for gross error removal, and then determining the convergence condition of the star point residuals. The expression for the standard error is as follows: ， in, This represents the mean square error of the inter-satellite angle residual. Let be the angle residual between star point i and star point j, and n be the number of star points in the current star map.

Images