A multi-view on-orbit geometric calibration method for agile imaging satellite

Abstract

A multi-view on-orbit geometric calibration system for agile imaging satellite: in view of the problem that imaging parameters change during the launch and on-orbit operation of agile imaging satellite, affecting the high-precision positioning of images, a rigorous imaging geometric model is constructed according to the collinear imaging principle, and the on-orbit geometric calibration of agile satellite imaging parameters is carried out through the internal and external step-by-step calibration strategy. The change of imaging angle of agile satellite causes the change of image imaging height, and the influence of atmospheric refraction is intensified, so that the imaging image calibration coefficient of the nadir point of the agile satellite is directly used for the geometric calibration of the inclined imaging image, and the uncontrolled positioning precision is reduced. By establishing the relationship between the change of imaging angle and the compensation of external calibration coefficient, the external calibration coefficient of the inclined imaging image is compensated, and the problem of the decline of uncontrolled positioning precision caused by the change of imaging height and the intensified influence of atmospheric refraction on the inclined imaging image is relieved. The present application can effectively realize the on-orbit geometric calibration of agile imaging satellite, and provide support for the subsequent high-precision geometric processing of images.

Description

TECHNICAL FIELD

[0001] The application relates to a multi-view on-orbit geometric calibration method for agile imaging satellites, and belongs to the technical field of photogrammetry and remote sensing. BACKGROUND

[0002] High-precision positioning of satellite remote sensing images is the basis for satellite image geometric processing, the basis for surface feature information quantification, and a key indicator for realizing satellite image mapping applications. To achieve high-precision positioning of satellite images, it is necessary to accurately obtain various imaging parameters of satellite sensors, and then restore a rigorous geometric imaging model of the satellite sensor and the surface features at the imaging time, and accurately establish the mapping relationship between the satellite image and the surface features.

[0003] Due to the device wear caused by long-term use of the satellite sensor and the influence of the operating environment and various disturbance forces during on-orbit operation, the imaging parameters of the agile satellite deviate from the initial values of the imaging parameters calibrated in the laboratory before launch, and finally result in reduced image uncontrolled positioning accuracy. Through on-orbit geometric calibration of systematic external errors and internal errors in the imaging system, a rigorous mathematical relationship between the ground point image coordinates and the ground coordinates thereof can be accurately established, and the image positioning accuracy can be improved.

[0004] With the improvement of the hardware capability of the sensor load, the satellite maneuverability is continuously improved. The agile imaging satellite can be multi-angle maneuvered in the rolling axis and the pitching axis. The complex imaging mode of the single camera agile satellite supports multiple observation modes at the same time, and puts forward new requirements for the uncontrolled positioning technology of the agile imaging satellite. The conventional on-orbit geometric calibration is based on the low imaging angle image calibration of the sub-satellite point, and then the calibration coefficients of the sub-satellite point image are obtained. The low-angle satellite imaging image can be well calibrated based on the calibration coefficients of the sub-satellite point image. However, the change of the imaging elevation of the satellite imaging angle and the influence of the atmospheric refraction cannot be ignored, so that the calibration coefficients of the sub-satellite point image at the low imaging angle cannot be directly used for the parameter calibration of the satellite large-angle imaging image, which is manifested as the greater the change of the imaging angle of the image compared with the sub-satellite point, the more the uncontrolled positioning accuracy of the image is reduced.

[0005] The existing strict imaging model based on the physical model has a clear meaning, however, these parameters on the long-focus, linear array CCD satellite camera show strong correlation, and direct solving easily causes unstable solving; moreover, the geometric calibration coefficients of the agile satellite based on the sub-satellite point image are directly used for the large-angle imaging image, and due to the change of the imaging height and the influence of the atmospheric refraction, the uncontrolled positioning accuracy is reduced. SUMMARY

[0006] The problem to be solved by the present application is that part of the imaging parameters of the linear array push-broom agile imaging satellite changes during on-orbit operation, the satellite positioning accuracy decreases, and the accuracy of the no-control positioning of the nadir image of the inclined imaging satellite decreases when using the calibration coefficient of the star image, and a multi-view on-orbit geometric calibration method for agile imaging satellites is proposed. On the one hand, the systematic external error and internal error in the rigorous geometric imaging model are calibrated to obtain the geometric calibration coefficient of the nadir image; and then the decline in the no-control positioning accuracy of the inclined image caused by the change in the imaging height is compensated based on the calibration coefficient of the nadir image, so as to realize the on-orbit high-precision geometric calibration of the linear array push-broom agile satellite.

[0007] The technical scheme of the present application is a multi-view on-orbit geometric calibration method for agile imaging satellites, which comprises:

[0008] Step 1, select the agile satellite optical image to be calibrated and the DOM image and DEM image of the same coverage area thereof;

[0009] Step 2, based on the matching algorithm, match the selected agile satellite optical image with the DOM image and the DEM image to obtain dense control points;

[0010] Step 3, use the imaging parameters at the dense control points to construct a geometric calibration model of the linear array push-broom agile satellite based on the pointing angle of the probe element, use a step-by-step calibration method to solve the internal and external calibration parameters, and obtain the geometric calibration result of the agile satellite camera;

[0011] Step 4, when only considering the change in the imaging pitch angle, based on the results obtained in steps 1-3, calculate the external calibration coefficient and the internal calibration coefficient of the agile satellite nadir image, and select agile satellite imaging images of different pitch angles; based on the internal calibration coefficient of the nadir point, externally calibrate the imaging images of different imaging angles to obtain the external calibration coefficient of the imaging images of different imaging angles;

[0012] Step 5, according to the external calibration coefficient of the imaging images of different angles obtained in step 4, calculate the change in the pitch angle and the change in the external calibration coefficient of the imaging images of different angles compared with the nadir image;

[0013] Step 6, according to the change in the pitch angle and the change in the external calibration coefficient of the imaging images of different angles compared with the nadir image obtained in step 5, establish a relationship between the change in the pitch angle and the change in the external calibration coefficient according to a polynomial fitting, and solve the fitting coefficient of the roll angle direction and the pitch angle direction external calibration correction according to the least square principle;

[0014] When only considering the imaging angle of agile satellite changing in roll direction, steps 4-6 are performed to obtain the relationship between the imaging roll angle change and the imaging exterior calibration coefficient change, and then the relationship between the roll angle change and the roll angle exterior calibration coefficient correction and the pitch angle exterior calibration coefficient correction when the roll angle changes is calculated to obtain the fitting coefficients of the exterior calibration correction in the roll direction and the pitch direction.

[0015] The obtained dense control points should be uniformly distributed on the agile satellite images.

[0016] Using the imaging parameters at the dense control points, a geometric calibration model of the linear array push-broom agile satellite based on the pointing angle of the probe element is constructed, and the interior and exterior calibration parameters are calculated using the step-by-step calibration method to obtain the geometric calibration results of the agile satellite camera, including:

[0017] Step 3.1, obtaining the initial interior parameters of the attitude, orbit, orbit elements and camera probe element pointing angle at the dense control points, and constructing a geometric calibration model according to the parameters;

[0018] Step 3.2, common calibration of each CCD in the exterior calibration, establishing the unit vector relationship between the object vector and the image vector, and constructing an exterior calibration error equation to calculate the exterior calibration coefficient;

[0019] Step 3.3, separate calibration of each CCD in the interior calibration, reading the calculated exterior calibration coefficient, constructing an interior calibration error equation, and calculating the polynomial fitting coefficients of the interior parameters;

[0020] Step 3.4, calculating the residual at the dense control points obtained in step 2 by using the exterior calibration coefficient and the interior calibration coefficient obtained in steps 3.2 and 3.3, if the residual is less than a preset residual threshold, completing the calibration calculation, and obtaining the final geometric calibration results of the agile satellite camera, otherwise returning to step 3.2.

[0021] The geometric calibration model constructed in step 3.1 has the following specific form:

[0022]

[0023]

[0024] In the formula, (φ x ,φ y ) is the pointing angle of the imaging pixel corresponding to the ground control point in the camera coordinate system; λ is the scale factor; is the rotation matrix from the camera coordinate system to the satellite body coordinate system; is the rotation matrix from the satellite body coordinate system to the satellite orbit coordinate system; is the rotation matrix from the satellite orbit coordinate system to the J2000 coordinate system; is the rotation matrix from the J2000 coordinate system to the WGS84 coordinate system; RU For the generalized compensation matrix, ; (X G ,Y G Z G (X) represents the spatial rectangular coordinates of the control point in the WGS coordinate system; S ,Y S Z S (m0, m1, m2, m3) represents the spatial rectangular coordinates of the GNSS antenna phase center in the WGS84 coordinate system; (m0, m1, m2, m3) represents the fitting coefficients of the intrinsic parameters along the track direction; (k0, k1, k2, k3) represents the fitting coefficients of the intrinsic parameters perpendicular to the track direction; and S represents the column of the image pixel corresponding to the control point.

[0025] In step 3.2, the unit vector relationship between the object vector and the image vector is established, and the external calibration error equation is constructed to calculate the external calibration coefficients, as follows:

[0026]

[0027] Where (X,Y,Z) are the object-space unit vectors, (x,y,z) are the phase-space unit vectors, and the scaling factor λ is 1; for equation (3), the following calculations are performed:

[0028]

[0029]

[0030]

[0031] in( ω and κ are the roll angle, pitch angle, and yaw angle, respectively; from equation (3) in Performing Taylor expansion at this point yields:

[0032]

[0033] For the nth ground control point (X) n ,Y n Z n ,x n ,y n According to equation (7), we have:

[0034]

[0035] Solve for the external scaling coefficients using the least squares adjustment principle.

[0036]

[0037] Among them Place P is a weight matrix of size 2n×2n.

[0038] In step 3.3, each CCD is calibrated individually for internal calibration. The calculated external calibration coefficients are read, the internal calibration error equation is constructed, and the polynomial fitting coefficients of the internal parameters are calculated, as follows:

[0039]

[0040] Where (X) I ,Y I Z I ) is the object-space unit vector. According to equation (10), the orbital intrinsic parameter φ of the probe position corresponding to the nth control point on a single CCD can be solved. xn With the internal parameter φ of the vertical rail yn for:

[0041]

[0042] Based on equations (2) and (11), for the nth control point on the ground:

[0043]

[0044] In formula (12) m = [m0, m1, m2, m3] T k = [k0, k1, k2, k3] T The internal calibration coefficients of the CCD along the track direction and perpendicular to the track direction are calculated sequentially based on the least squares adjustment, where the coefficients along the track direction are:

[0045]

[0046] Where in m0 = [m 00 ,m 01 ,m 02 ,m 03 ] T Place U is a weight matrix of size n×n; the internal scaling coefficients in the perpendicular direction are:

[0047]

[0048] In equation (14), k0 = [k 00 ,k 01 ,k 02 ,k 03 ] T Place T is a weight matrix of size n×n.

[0049] The calculation of the pitch angle change between the image at different angles and the nadir image includes: considering only the pitch angle change between the tilt image and the nadir image, and setting the nadir roll angle external calibration coefficient as follows: The external calibration coefficient for the nadir pitch angle is ω0, the external calibration coefficient for the nadir yaw angle is κ0, the nadir pitch angle is α0, and the pitch angle of the tilt imaging image is α. n Then the change in pitch angle for:

[0050]

[0051] Calculate the change in external calibration coefficients between images at different angles and nadir images, including: considering only the atmospheric refraction effect caused by changes in imaging pitch angle, and compensating for changes in imaging altitude by external calibration coefficients in both pitch and roll directions; then, calculate the change in external calibration coefficients of the tilt image compared to the nadir image in the roll angle direction. The change in pitch angle directional coefficient Δω 0n The calculation is as follows:

[0052]

[0053] in, With ω n Let the image imaging angle be α n The external calibration coefficients for the roll angle and pitch angle at that time.

[0054] Establish the relationship between the pitch angle change and the external calibration coefficient change. Based on the least squares principle, calculate the fitting coefficients of the external calibration correction in the roll angle direction and the pitch angle direction, including:

[0055] For imaging images at different tilt angles, the external calibration coefficient correction and angle change calculated based on the in-satellite calibration coefficients are as follows: Based on the polynomial fitting coefficients, the relationship between the angle change and the correction amount of the external calibration coefficients is established as follows:

[0056]

[0057] Based on the least squares algorithm, the fitting coefficients for the external calibration correction in the roll angle direction are obtained as r = [r0, r1, r2, r3]. T The fitting coefficients for the external calibration correction in the pitch direction are q = [q0, q1, q2, q3]. T .

[0058] Compared with the prior art, the present invention has the following advantages:

[0059] 1) According to the optical camera imaging geometry, the agile imaging satellite on-orbit geometric calibration model based on the probe element pointing angle is constructed. In the external calibration, the generalized compensation matrix is introduced, which can equivalent transform the external system line element error of the linear array push-broom agile camera into the angle element error, effectively eliminating the attitude, orbit measurement error, GPS eccentric error and camera installation error, etc. The internal calibration uses the probe element pointing angle model to transform the complex physical parameters of the camera into the polynomial coefficients of the probe element pointing angle model, and the control points are arranged in the along-track direction and the cross-track direction to quickly solve the internal calibration coefficients. The calibration method of external calibration first and internal calibration later is also conducive to reducing the influence of the correlation of complex imaging parameters, and realizing the on-orbit geometric calibration task of the linear array push-broom agile satellite.

[0060] 2) Through the simulation of the imaging pitch and roll direction angle change and the external calibration coefficient compensation relationship, the polynomial coefficients are used to compensate the atmospheric refraction influence caused by the change of imaging angle, the change of imaging elevation and the change of imaging angle, and the imaging angle image external calibration coefficient is obtained, which is more suitable for the calibration coefficient of the angle imaging image, and the on-orbit geometric calibration task of the multi-angle imaging image is realized, and the uncontrolled positioning accuracy of the large-angle agile satellite image is effectively improved. BRIEF DESCRIPTION OF DRAWINGS

[0061] Figure 1 It is the on-orbit geometric calibration flowchart in the application.

[0062] Figure 2 It is the probe element pointing angle model schematic diagram involved in the application.

[0063] Figure 3 It is the inclination image external calibration compensation coefficient calculation flowchart in the application. DETAILED DESCRIPTION

[0064] A multi-view on-orbit geometric calibration method for agile imaging satellites, comprising the following steps:

[0065] Step 1, selecting the agile satellite optical image to be calibrated and the DOM image and DEM image covering the same area.

[0066] Step 2, realizing the image matching of the agile satellite optical image, the DOM image and the DEM image based on the matching algorithm, obtaining dense control points, and further screening the obtained control points, the obtained control point data should be uniformly distributed on the agile satellite image.

[0067] Step 3, constructing the linear array push-broom agile satellite geometric calibration model based on the probe element pointing angle using the imaging parameters at the control points, using the step-by-step calibration method to solve the internal and external calibration parameters, and obtaining the geometric calibration results of the agile satellite camera.

[0068] Step 3.1, according to the row and column where the control point is located, the attitude, the track, the track root number and the inner parameter of the initial probe element pointing angle at the control point are obtained by interpolation.

[0069] Step 3.2, common calibration of each piece of CCD, according to the principle of collinear imaging, the unit vectors of the object vector and the image vector are calculated to construct the external calibration error equation to calculate the external calibration coefficient.

[0070] Step 3.3, individual calibration of each piece of CCD, the external calibration coefficient calculated is read, the inner parameter at the control point is calculated from the image side, the inner parameter polynomial fitting coefficient is calculated by constructing the inner calibration error equation, and the inner orientation elements at the probe element on each piece of CCD and the inner orientation elements at the control point are recalculated.

[0071] Step 3.4, the residual at the control point is calculated through the external and inner calibration coefficients obtained in steps 3.1 and 3.2, if the residual is less than the preset residual threshold, the calibration calculation is completed, and the final geometric calibration result of the agile satellite camera is obtained, otherwise steps 3.2 and 3.3 are repeated.

[0072] In step 3, according to the imaging geometric relationship of the linear array push-broom agile satellite, the on-orbit geometric calibration model of the linear array push-broom agile satellite based on the probe element pointing angle is constructed, as shown in formula (1) and (2):

[0073]

[0074]

[0075] In the formula, (φ x ,φ y ) is the probe element pointing angle of the corresponding imaging element of the ground control point in the camera coordinate system; λ is the scale factor; is the rotation matrix from the camera coordinate system to the satellite body coordinate system; is the rotation matrix from the satellite body coordinate system to the satellite orbit coordinate system; is the rotation matrix from the satellite orbit coordinate system to the J2000 coordinate system; is the rotation matrix from the J2000 coordinate system to the WGS84 coordinate system; R U is the generalized compensation matrix; (X G ,Y G ,Z G ) is the space rectangular coordinate of the control point in the WGS coordinate system; (X S ,Y S ,Z S(m0, m1, m2, m3) represents the spatial rectangular coordinates of the GNSS antenna phase center in the WGS84 coordinate system; (m0, m1, m2, m3) represents the fitting coefficients of the intrinsic parameters along the track direction; (k0, k1, k2, k3) represents the fitting coefficients of the intrinsic parameters perpendicular to the track direction; and S represents the column of the image pixel corresponding to the control point.

[0076] In step 3.2, based on the principles of collinear imaging and vector equality, the relationship between the object-side unit vector and the image-side unit vector is established according to equation (3) as follows:

[0077]

[0078] Where (X,Y,Z) are the object-side unit vectors, (x,y,z) are the image-side unit vectors, and the scale factor λ is 1. Equation (3) has the following calculation:

[0079]

[0080]

[0081]

[0082] in( ω and κ represent the roll angle, pitch angle, and yaw angle, respectively. From equation (3)... Performing Taylor expansion at this point yields:

[0083]

[0084] For the nth ground control point (X) n ,Y n Z n ,x n ,y n According to equation (7), we have:

[0085]

[0086] Solve for the unknowns using the least squares adjustment principle. At least two control points should be used:

[0087]

[0088] Among them Place P is a weight matrix of size n×n.

[0089] Step 3.2: After establishing a rigorous imaging model based on external calibration, perform internal calibration on each CCD. Based on the principle of vector equality, we have:

[0090]

[0091] Where (X) I ,Y I Z I ) is the object-space unit vector. According to equation (10), the orbital intrinsic parameter φ of the probe position corresponding to the nth control point on a single CCD can be solved. xn With the internal parameter φ of the vertical rail yn for:

[0092]

[0093] Based on equations (2) and (11), for the nth control point on the ground:

[0094]

[0095] In formula (12) m = [m0, m1, m2, m3] T k = [k0, k1, k2, k3] T Based on least squares adjustment, the internal calibration coefficients of the CCD along the track direction and perpendicular to the track direction can be calculated sequentially, where the coefficients along the track direction are:

[0096]

[0097] Where in m0 = [m 00 ,m 01 ,m 02 ,m 03 ] T Place U is a weight matrix of size n×n; the internal scaling coefficients in the perpendicular direction are:

[0098]

[0099] In equation (14), k0 = [k 00 ,k 01 ,k 02 ,k 03 ] T Place T is a weight matrix of size n×n.

[0100] Step 4: Based on steps 1-3, the imaging geometric parameters for the nadir image calibration of the Agile satellite can be obtained. Considering only the change in image elevation angle, multiple Agile satellite images with different elevation angles are selected. Since the satellite installation relationship does not change within a certain period of time, external calibration is performed on the images with different imaging angles based on the nadir internal calibration coefficients.

[0101] Step 4.1: According to equations (1)-(9), perform external calibration on the tilt imaging image based on the nadir point internal calibration coefficients to obtain the roll angle external calibration coefficients of the image during tilt imaging. ω n and κ n .

[0102] Step 5, calculate the change of the pitch angle and the outer scale factor of the images of different angles compared with the nadir image;

[0103] In step 5, only the change of the pitch angle of the oblique angle imaging image compared with the nadir image is considered, and the outer scale factor of the nadir roll angle is set as The outer scale factor of the nadir pitch angle is ω0, the outer scale factor of the nadir yaw angle is κ0, the nadir pitch angle is α0, and the pitch angle of the oblique angle imaging image in step 4.1 is α n , then the change of the pitch angle is The calculation is as follows:

[0104]

[0105] In step 5, only the atmospheric refraction caused by the change of the imaging pitch angle and the change of the imaging height caused by the outer scale factors of the pitch and roll directions are considered, and then the change of the outer scale factor of the oblique angle imaging image in the roll angle direction compared with the outer scale factor of the nadir is and the change of the outer scale factor in the pitch angle direction is Δω 0n The calculation is as follows:

[0106]

[0107] Step 6, establish the relationship between the angle change and the outer scale factor correction according to the polynomial fitting.

[0108] In step 6, for different oblique angle imaging images, the outer scale factor correction and the angle change calculated based on the inner scale factor of the nadir are According to the polynomial fitting coefficients, the relationship between the angle change and the outer scale factor correction is as follows

[0109]

[0110] According to formulas (12)-(13), combined with the least square algorithm, the fitting coefficients of the outer scale correction in the roll angle direction are r = [r0, r1, r2, r3] T , and the fitting coefficients of the outer scale correction in the pitch angle direction are q = [q0, q1, q2, q3] T .

[0111] When the imaging angle of the agile satellite only changes in the roll angle direction, the relationship between the imaging roll angle change and the image exterior calibration coefficient correction amount can be obtained by referring to steps 4-6, and then the relationship between the roll angle change, the roll angle exterior calibration coefficient correction amount and the pitch angle exterior calibration coefficient correction amount when the roll angle changes can be calculated.

[0112] By calculating the changes in the pitch angle and the roll angle between the imaging image and the nadir imaging image, the exterior calibration coefficient of the imaging image can be compensated according to the exterior calibration correction amount fitting coefficient, and then the uncontrolled positioning accuracy of the imaging image can be improved, and the problem of the uncontrolled positioning accuracy reduction caused by the increased atmospheric refraction and the imaging height change can be reduced.

[0113] When the new agile satellite image data is calibrated, the exterior and interior calibration coefficients of the new data can be recalculated by using the calculated nadir interior and exterior calibration coefficients and the exterior calibration compensation coefficient, and the subsequent high-precision geometric product production of the image can be assisted.

[0114] The specific embodiments of the present application are further illustrated below in combination with the flowchart and the embodiments shown in Figure 1 As shown in Figure 1 , the detailed description of each step of the embodiment flow is as follows:

[0115] Step 1, image data preparation, selecting the camera agile satellite optical image to be calibrated and the DOM image and the DEM image of the same coverage area of the agile satellite, for subsequent matching or manual production of control points.

[0116] Step 2, obtaining control points based on matching algorithm, based on the selected base image and the agile satellite image to be matched, performing automatic matching or manual point picking of control points. In automatic matching, the image matching of the to-be-matched image, the DOM and the DEM data is realized based on feature point recognition and gray cross-correlation algorithm, dense control points are obtained, and the obtained control points are further screened by random sample consensus algorithm. The obtained control point data should be uniformly distributed on the agile satellite image.

[0117] Step 3, obtaining image imaging parameters, constructing a geometric calibration model based on the pointing angle of the probe element. First, read the attitude, orbit and time data of the corresponding satellite transmitted by the selected agile satellite. According to the positions of the probe elements on the CCD in the camera coordinate system as shown in Figure 2 , the vertical and along-track pointing angles of the probe elements can be calculated as follows:

[0118]

[0119] In formula (1), the probe element pointing angle of the imaging pixel corresponding to the ground control point in the camera coordinate system, x ccd and y ccdis the position of the CCD detector in the camera coordinate system, and f is the focal length of the camera.

[0120] The imaging parameters at the control points are obtained by interpolation. During the on-orbit period of the satellite, the recording frequency of the orbit data and the recording frequency of the attitude data do not match the imaging push-broom speed, so it is necessary to interpolate the orbit and the attitude for each row. According to the row in which the control point is located, the corresponding row is interpolated, and the corresponding attitude and orbit are further interpolated according to the row time to calculate the orbit root. According to the column in which the control point is located, the parameters in the control point are interpolated. Thus, the imaging parameters at each control point for on-orbit geometric calibration are obtained.

[0121] A geometric calibration model of the linear array push-broom agile satellite based on the pointing angle of the detector is constructed using the imaging parameters at the control points. The internal and external calibration parameters are solved by using a step-by-step calibration method to obtain the geometric calibration result of the camera of the agile satellite. The implementation is as follows. First, the external calibration parameters are solved by least squares adjustment to restore the attitude of the camera coordinate system in space. Then, the internal calibration parameters are solved by least squares adjustment to determine the pointing angle of each detector of each slice CCD in the camera coordinate system.

[0122] Step 3.1, the external calibration uses all the CCD slices in the camera to calibrate together. According to the collinear imaging principle, the unit vectors of the object vector and the image vector are calculated to construct the external calibration error equation to calculate the external calibration coefficient.

[0123] Step 3.2, the internal calibration calibrates each slice CCD of the camera separately. The external calibration coefficient calculated in step 3.1 is read, and the internal parameters at the control point are calculated from the image side to construct the internal calibration error equation to calculate the polynomial fitting coefficient of the internal parameters, and the internal orientation elements at the detectors of each slice CCD and the internal orientation elements at the control point are recalculated.

[0124] Step 3.3, the residual at the control point is calculated by the external and internal calibration coefficients obtained in steps 3.1 and 3.2. If the residual is less than a preset residual threshold, the calibration solution is completed, and the final geometric calibration result of the camera of the agile satellite is obtained. Otherwise, steps 3.2 and 3.3 are repeated.

[0125] Step 4, the nadir external calibration coefficient and the internal calibration coefficient of the agile satellite are obtained according to steps 1-3. Only the change of the imaging pitch angle of the image is considered, and the imaging images of the agile satellite at different pitch angles are selected. Since the satellite installation relationship does not change within a certain time, the internal parameters of the satellite tend to be stable. Therefore, the external calibration of the agile satellite imaging images at different imaging angles is performed based on the internal calibration coefficient of the nadir to obtain the external calibration coefficient of the agile satellite imaging at different imaging angles.

[0126] Step 5, the pitch angle change and the external calibration coefficient change of the images at different angles compared with the images at the nadir are calculated according to the external calibration coefficients of the images at different angles obtained in step 4.

[0127] Step 6, according to the pitch angle variation and the outer calibration coefficient variation of the images of different angles obtained in step 5 compared with the images of the subsatellite point, a relationship between the pitch angle variation and the outer calibration coefficient variation is established according to polynomial fitting. The least square method is used to solve the equation to obtain the fitting coefficients of the roll angle outer calibration correction and the fitting coefficients of the pitch angle outer calibration correction.

[0128] When the imaging angle of the agile satellite only changes in the roll angle direction, the relationship between the image imaging roll angle variation and the image outer calibration coefficient correction can be obtained by referring to steps 6-9, and then the relationship between the roll angle variation and the outer calibration coefficient correction when the roll angle changes can be calculated to obtain the fitting coefficients of the roll angle outer calibration correction and the fitting coefficients of the pitch angle outer calibration correction.

[0129] Through the calculation of the pitch angle and the roll angle variation between the inclined angle imaging image and the subsatellite point imaging image, the outer calibration coefficient of the inclined angle imaging image can be compensated according to the fitting coefficients of the roll angle and the pitch angle outer calibration correction, and then the uncontrolled positioning accuracy of the inclined angle imaging image is improved, and the problem of uncontrolled positioning accuracy reduction caused by the aggravation of atmospheric refraction and the change of imaging height is reduced.

[0130] When new agile satellite image data is calibrated, the outer and inner calibration coefficients of the new data can be calculated by using the calculated subsatellite inner and outer calibration coefficients and the outer calibration compensation coefficients, and the subsequent high-precision geometric product production of the image is assisted.

[0131] The contents not described in detail in the specification of the present application are the known technology of those skilled in the art.

Claims

1. A multi-view on-orbit geometric calibration method for agile imaging satellites, characterized in that, The method comprises the following steps: Step 1, selecting agile satellite optical images to be calibrated and DOM images and DEM images of the same coverage area of the agile satellite optical images; Step 2, matching the selected agile satellite optical images with the DOM images and the DEM images based on a matching algorithm to obtain dense control points; Step 3, constructing a geometric calibration model of the linear array push-broom agile satellite based on a pointing angle of a probe element by using imaging parameters at the dense control points, and solving internal and external calibration parameters by using a step-by-step calibration method to obtain geometric calibration results of the agile satellite camera; Step 4, when only considering changes in imaging pitch angles of the images, calculating external calibration coefficients and internal calibration coefficients of the agile satellite sub-satellite point images based on the results obtained in steps 1-3, and selecting agile satellite imaging images of different pitch angles; and performing external calibration on the imaging images of different imaging angles based on the internal calibration coefficients of the sub-satellite point to obtain external calibration coefficients of the imaging images of different imaging angles; Step 5, calculating changes in pitch angles and changes in external calibration coefficients of the imaging images of different angles compared with the sub-satellite point images according to the external calibration coefficients of the imaging images of different angles obtained in step 4; Step 6, establishing a relationship between the changes in pitch angles and the changes in external calibration coefficients according to the changes in pitch angles and the changes in external calibration coefficients of the imaging images of different angles compared with the sub-satellite point images obtained in step 5, and solving fitting coefficients of external calibration correction amounts in the roll angle direction and the pitch angle direction according to the least square principle; When only considering changes in imaging angles of the agile satellite in the roll angle direction, the relationship between the changes in the roll angle and the changes in the external calibration coefficients of the images is obtained by executing steps 4-6, and then the relationship between the changes in the roll angle and the correction amounts of the roll angle external calibration coefficients and the pitch angle external calibration coefficients when the roll angle changes can be calculated, and fitting coefficients of the external calibration correction amounts in the roll angle direction and the pitch angle direction are obtained.

2. The multi-view on-orbit geometric calibration method for agile imaging satellites according to claim 1, wherein, The dense control points obtained should be uniformly distributed on the agile satellite images.

3. The multi-view on-orbit geometric calibration method for agile imaging satellites according to claim 1, wherein, The geometric calibration model of the linear array push-broom agile satellite based on the pointing angle of the probe element is constructed by using the imaging parameters at the dense control points, and the internal and external calibration parameters are solved by using the step-by-step calibration method to obtain the geometric calibration results of the agile satellite camera, which comprises the following steps: Step 3.1, obtaining initial internal parameters of the attitude, the orbit, the orbit root number and the camera probe element pointing angle at the dense control points, and constructing a geometric calibration model according to the parameters; Step 3.2, performing common calibration of each piece of CCD to establish a unit vector relationship between the object vector and the image vector, and constructing an external calibration error equation to calculate the external calibration coefficients; Step 3.3, performing separate calibration of each piece of CCD to read the calculated external calibration coefficients, construct an internal calibration error equation, and calculate polynomial fitting coefficients of the internal parameters; Step 3.4, calculating the residual at the dense control points obtained in step 2 by using the external calibration coefficients and the internal calibration coefficients obtained in steps 3.2 and 3.3, and if the residual is less than a preset residual threshold, completing the calibration solution to obtain the final geometric calibration results of the agile satellite camera, otherwise returning to step 3.

2.

4. The multi-view on-orbit geometric calibration method for agile imaging satellites according to claim 3, wherein, The geometric calibration model constructed in step 3.1 has the following specific form: In the formula, (φ x ,φ y ) is the pointing angle of the corresponding imaging pixel of the ground control point in the camera coordinate system. λ is a scale factor; R is a rotation matrix from the camera coordinate system to the satellite body coordinate system; R is a rotation matrix from the satellite body coordinate system to the satellite orbit coordinate system; R is a rotation matrix from the satellite orbit coordinate system to the J2000 coordinate system; R is a rotation matrix from the J2000 coordinate system to the WGS84 coordinate system; U R is a generalized compensation matrix; G , Y G , Z G are the spatial rectangular coordinates of the control point in the WGS coordinate system; S , Y S , Z S are the spatial rectangular coordinates of the GNSS antenna phase center in the WGS84 coordinate system; (m0, m1, m2, m3) are the fitting coefficients of the along-track direction internal parameters; (k0, k1, k2, k3) are the fitting coefficients of the across-track direction internal parameters; S is the column in which the corresponding image point of the control point is imaged.

5. The multi-view on-orbit geometric calibration method for agile imaging satellites according to claim 1, wherein, The calculation of the change of the pitch angle of the different angle image compared with the sub-sky point image comprises: only considering the change of the pitch angle of the tilt angle imaging image compared with the sub-sky point image, setting the sub-sky point roll angle outer calibration coefficient as The sub-sky point pitch angle outer calibration coefficient is ω0, the sub-sky point yaw angle outer calibration coefficient is κ0, the sub-sky point pitch angle is α0, and the tilt angle imaging image pitch angle is α n The change of the pitch angle is is:

6. The multi-view on-orbit geometric calibration method for agile imaging satellite according to claim 5, characterized in that, The variation of the exterior calibration coefficient of the different angle image compared with the exterior calibration coefficient of the nadir image is calculated, including: only considering the atmospheric refraction caused by the change of the imaging pitch angle, the change of the imaging height is compensated by the exterior calibration coefficient of the pitch and roll directions, and the variation of the exterior calibration coefficient of the tilt imaging image compared with the exterior calibration coefficient of the nadir image is the variation of the exterior calibration coefficient in the roll angle direction and the variation of the coefficient in the pitch angle direction Δω 0n The calculation is as follows: wherein, with ω n the roll angle outer scaling factor and the pitch angle outer scaling factor for an image imaging angle of α n when.

7. The multi-view on-orbit geometric calibration method for agile imaging satellites according to claim 6, wherein, The relationship between the pitch angle change amount and the outer calibration coefficient change amount is established, and the fitting coefficients of the outer calibration correction amount of the roll angle direction and the pitch angle direction are obtained by solving according to the least square principle, including: For different imaging images with different tilt angles, the outer calibration coefficient correction and the angle change based on the inner calibration coefficient of the subsatellite point are is the tilt angle change; according to the polynomial fitting coefficient, the relationship between the angle change and the outer calibration coefficient correction is established as follows: According to the least square algorithm, the fitting coefficients of the roll angle direction outer scaling correction amount r = [r0, r1, r2, r3] are solved T , and the fitting coefficients of the pitch angle direction outer scaling correction amount q = [q0, q1, q2, q3] T .

Images