Geometric positioning processing method based on linear array push-broom imaging satellite

Abstract

The invention belongs to the technical field of geometric positioning of satellites, and particularly discloses a geometric positioning processing method based on a linear array push-broom imaging satellite, and the method comprises the steps: constructing an offset matrix and a directional angle polynomial model, the directional angle polynomial model is used for simulating the change process of a directional angle along with the detection element image column, and the directional angle is the angle of the imaging light of the detection element relative to the direction perpendicular to the track and the direction along the track; and based on the offset matrix and the directional angle polynomial model, a rigorous geometric positioning model is adjusted, a joint calibration model is obtained, and the rigorous geometric positioning model is used for representing the geometric relationship between the image point and the ground point. According to the invention, the internal and external orientation elements are calibrated, so that the adjusted rigorous geometric positioning model can more accurately represent the geometric relationship between the image point and the ground point, and the geometric positioning precision based on the linear array push-broom imaging satellite is improved.

Description

Technical Field

[0001] The present application belongs to the technical field of satellite geometric positioning, and more specifically, to a geometric positioning processing method based on a linear array push-broom imaging satellite. Background Art

[0002] Satellite in-orbit geometric calibration is a key step in fully utilizing the geometric performance of satellites. The geometric positioning model can provide geographic location information for satellite images and is the key to subsequent measurement, positioning and other applications of remote sensing images. A general and rigorous geometric positioning model can be constructed for linear array push-broom imaging satellites, through which the geometric relationship between image points and real ground points can be established. However, the construction process of the positioning model involves the combination of many internal and external orientation element information. External orientation element errors and internal orientation element errors will affect the positioning accuracy of the geometric positioning model. How to improve the geometric positioning accuracy based on linear array push-broom imaging satellites is a technical problem that needs to be solved urgently in this field. Summary of the invention

[0003] In view of the defects of the prior art, the purpose of this application is to improve the geometric positioning accuracy of linear array push-broom imaging satellites.

[0004] To achieve the above objectives, in a first aspect, the present application provides a geometric positioning processing method based on a linear array push-broom imaging satellite, the method comprising: Construct a bias matrix and a pointing angle polynomial model. The bias matrix is ​​used to correct the pointing direction of the imaging light of the detector. The pointing angle polynomial model is used to simulate the process in which the pointing angle changes with the detector image sequence. The pointing angle is the angle of the imaging light of the detector relative to the direction perpendicular to the track and along the track. Based on the bias matrix and the pointing angle polynomial model, the rigorous geometric positioning model is adjusted to obtain the joint calibration model. The rigorous geometric positioning model is used to characterize the geometric relationship between image points and ground points.

[0005] In a second aspect, the present application provides a geometric positioning processing device based on a linear array push-broom imaging satellite, comprising: A construction module is used to construct a bias matrix and a pointing angle polynomial model. The bias matrix is ​​used to correct the pointing direction of the imaging light of the detector. The pointing angle polynomial model is used to simulate the process in which the pointing angle changes with the detector image sequence. The pointing angle is the angle of the imaging light of the detector relative to the direction perpendicular to the track and the direction along the track. The acquisition module is used to adjust the strict geometric positioning model based on the bias matrix and the pointing angle polynomial model to obtain the joint calibration model. The strict geometric positioning model is used to characterize the geometric relationship between the image points and the ground points.

[0006] In a third aspect, the present application provides an electronic device comprising: at least one memory for storing programs; and at least one processor for executing the programs stored in the memory. When the program stored in the memory is executed, the processor is used to execute the method described in the first aspect or any possible implementation of the first aspect.

[0007] In a fourth aspect, the present application provides a computer-readable storage medium, which stores a computer program. When the computer program runs on a processor, the processor executes the method described in the first aspect or any possible implementation of the first aspect.

[0008] In a fifth aspect, the present application provides a computer program product. When the computer program product runs on a processor, the processor executes the method described in the first aspect or any possible implementation of the first aspect.

[0009] In general, the above technical solutions conceived by this application have the following beneficial effects compared with the prior art: (1) By adjusting the rigorous geometric positioning model based on the bias matrix and the pointing angle polynomial model, the adjusted rigorous geometric positioning model can be used as a joint calibration model. The bias matrix in the joint calibration model can be used to correct the imaging light direction of the detector to calibrate the external orientation elements. The pointing angle polynomial model in the joint calibration model can be used to simulate the process in which the pointing angle changes with the detector image sequence to calibrate the internal orientation elements. By calibrating the internal and external orientation elements, the adjusted rigorous geometric positioning model can more accurately represent the geometric relationship between the image point and the ground point, thereby improving the geometric positioning accuracy of the linear array push-broom imaging satellite.

[0010] (2) By constructing a corresponding geometric calibration model based on the error sources and error characteristics of positioning accuracy, and using control point data to obtain The correction amount and the correction of the geometric positioning model can effectively improve the geometric positioning accuracy of the image. While improving the positioning accuracy, the consistency accuracy between multi-spectral images is taken into account, and the high-precision same-name points between the panchromatic and multispectral data are used for joint calibration to obtain The correction amount is used to modify the geometric positioning model to ensure the positioning consistency of the geometric positioning models of images in each band.

[0011] (3) For satellites with obvious attitude jitter, it is difficult to simulate the attitude models such as polynomial models accurately. However, the on-board attitude sampling frequency is high, and the penalized spline function model has a better fitting effect. It can accurately simulate the on-board attitude and eliminate the intersection error of the same-name points with obvious trends between spectral segments, effectively ensuring the positioning accuracy of the strict geometric positioning model. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 It is a flow chart of a geometric positioning processing method based on a linear array push-broom imaging satellite provided in an embodiment of the present application; Figure 2 is a schematic diagram of linear array push-scan imaging provided by an embodiment of the present application; Figure 3 is a schematic diagram of CCD rotation error provided in an embodiment of the present application; Figure 4 is a schematic diagram of the camera detector position provided in an embodiment of the present application; Figure 5 This is a schematic diagram of satellite attitude data information provided by an embodiment of the present application; Figure 6 It is a scatter plot of the dynamic change of satellite attitude angle provided in an embodiment of the present application; Figure 7 It is a polynomial model posture fitting diagram provided in an embodiment of the present application; Figure 8 It is a schematic diagram of the positioning error of the same-name points provided in the embodiment of the present application. DETAILED DESCRIPTION

[0013] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0014] In the embodiments of the present application, words such as "exemplary" or "for example" are used to indicate examples, illustrations or descriptions. Any embodiment or design described as "exemplary" or "for example" in the embodiments of the present application should not be interpreted as being more preferred or more advantageous than other embodiments or designs. Specifically, the use of words such as "exemplary" or "for example" is intended to present related concepts in a specific way.

[0015] In the description of the embodiments of the present application, unless otherwise specified, "multiple" means two or more than two. For example, multiple processing units refer to two or more processing units, etc.; multiple elements refer to two or more elements, etc.

[0016] The embodiments of the present application are described below in conjunction with the drawings in the embodiments of the present application.

[0017] Figure 1 is a flow chart of a geometric positioning processing method based on a linear array push-broom imaging satellite provided in an embodiment of the present application, such as Figure 1 As shown, the method includes the following steps S101 and S102.

[0018] Step S101, constructing a bias matrix and a pointing angle polynomial model, wherein the bias matrix is ​​used to correct the pointing direction of the imaging light of the detector, and the pointing angle polynomial model is used to simulate the process in which the pointing angle changes with the detector image sequence, and the pointing angle is the angle of the imaging light of the detector relative to the direction perpendicular to the track and the direction along the track; Step S102, based on the bias matrix and the pointing angle polynomial model, the strict geometric positioning model is adjusted to obtain a joint calibration model, where the strict geometric positioning model is used to characterize the geometric relationship between the image points and the ground points.

[0019] It can be understood that by adjusting the rigorous geometric positioning model based on the bias matrix and the pointing angle polynomial model, the adjusted rigorous geometric positioning model can be used as a joint calibration model. The bias matrix in the joint calibration model can be used to correct the imaging light direction of the detector to calibrate the external orientation elements. The pointing angle polynomial model in the joint calibration model can be used to simulate the process in which the pointing angle changes with the detector image series to calibrate the internal orientation elements. By calibrating the internal and external orientation elements, the adjusted rigorous geometric positioning model can more accurately characterize the geometric relationship between the image points and the ground points, thereby improving the geometric positioning accuracy of the linear array push-broom imaging satellite.

[0020] The geometric positioning processing method based on the linear array push-broom imaging satellite provided in the present application is exemplarily described below with reference to several examples.

[0021] (1) Construction of positioning model for linear array push-broom imaging camera and analysis of positioning error.

[0022] Linear array push-broom imaging satellites use linear array push-broom imaging cameras. Linear array push-broom imaging cameras are devices that use linear array image sensors for imaging. They capture image information by scanning line by line. When this camera is imaging, the photosensitive elements on the sensor expose the scene in sequence, converting the light signal into an electrical signal to form one-dimensional image data. As the camera or the scanned object moves, the accumulated data line by line is eventually combined into a complete two-dimensional image.

[0023] The geometric positioning model can provide geographic location information for satellite images and is the key to subsequent measurement and positioning of remote sensing images. The function of the geometric positioning model is to establish a mapping relationship between the coordinates of the image points on the image and the real geographic coordinates. During the imaging process, each pixel on the sensor can be simplified into a series of points to represent it. Therefore, according to the structural characteristics of different sensors, geometric positioning models corresponding to different satellite images can be established.

[0024] As the most accurate model in the positioning model, the rigorous geometric model is constructed strictly according to the geometric process of satellite imaging. In this process, it is necessary to obtain the real parameter information of the sensor devices when the satellite is in orbit, and the real-time orbit, attitude and other information of the satellite are needed to establish the connection between the coordinate systems.

[0025] (1-1) Construction of geometric positioning model of linear array push-scan imaging camera.

[0026] (1-1-1) Principles of rigorous (imaging) geometric models; For satellites using linear array push-broom imaging mode, the received electromagnetic wave signals of ground objects are converted into electrical signals through photosensors and recorded and stored. The electrical signals are converted into digital signals and transmitted to the ground receiving station through the digital transmission system. The linear array push-broom imaging camera images line by line as the satellite moves. The imaging diagram is shown in the figure below. Figure 2 shown.

[0027] During the satellite's in-orbit operation, the satellite remote sensing image is stitched together by line images as the satellite moves. At the moment of imaging each line of the image, the three points of the imaging payload CCD (Charge-Coupled Device), the optical perspective center and the target are collinear. Therefore, based on the collinearity equation, a rigorous geometric model can be established using the satellite sensor's internal orientation elements and the actual external orientation elements measured by the satellite platform during in-orbit operation to achieve precise positioning of the remote sensing image.

[0028] (1-1-2) Construction of rigorous geometric model; The general rigorous geometric positioning model is as follows: (1); In the formula, The transformation matrix from the camera coordinate system to the body coordinate system is obtained by measurement before the satellite is launched. represents the camera principal distance, is the CCD array's field deviation angle along the track, The position corresponding to the main visual axis (the vertical point through the principal point perpendicular to the CCD linear array), Indicates the corresponding position of each pixel in the CCD array, subscript Represents the index of the pixel. To find the size of the element, is the offset of the camera coordinate origin in the body coordinate system, It is the transformation matrix from the body coordinate system to the geocentric inertial coordinate system, where the J2000 (Julian 2000) coordinate system is the internationally commonly used geocentric inertial coordinate system. The attitude quaternion data measured by the satellite using star sensors in this application is the attitude in the satellite body coordinate system relative to the attitude in the J2000 coordinate system. The matrix can be determined using the quaternions parsed from the satellite auxiliary data. is the transformation matrix from the J2000 coordinate system to the WGS84 (World Geodetic System 1984) coordinate system at time t; is the coordinate position of the center origin of the GPS (Global Positioning System) module in the main body coordinate system, is the position vector of the phase center of the GPS module in the WGS84 coordinate system at time t; It indicates the coordinate value of the pixel point on the satellite image in the WGS84 coordinate system, and m is the scale factor.

[0029] The origin of the camera coordinate system is located at the camera projection center. The Z axis of the camera coordinate system is the camera principal optical axis and points to the positive direction of the focal plane. The Y axis of the camera coordinate system is parallel to the CCD array direction, and the X axis of the camera coordinate system roughly points to the satellite flight direction. The three-axis pointing meets the right-hand coordinate system rule.

[0030] The satellite body coordinate system is a coordinate system fixed to the satellite. Usually, the satellite center of mass is taken as the origin O, and the three main inertia axes of the satellite are taken as the XYZ axes. Among them, the OZ axis is positive when it points from the center of mass to the ground, the OX axis is positive when it points to the direction of the satellite's flight, and the OY axis is determined by the right-hand coordinate system rule.

[0031] This formula is a universal and rigorous geometric positioning model for linear array push-broom imaging satellites. Through this model, the geometric relationship between image points and real ground points can be established, which is the basis for geometric processing of remote sensing data.

[0032] (1-2) Positioning model error analysis; The construction process of the positioning model involves the combination of many internal and external orientation element information. Analyzing the error sources of the entire imaging process plays a key role in improving the positioning accuracy of the image and the consistency accuracy between multiple spectral segments. This is also the basis for subsequent optimization of the positioning model. By analyzing the process of building a geometric model, the error can be subdivided into internal orientation element error, external orientation element error, etc.

[0033] (1-2-1) Error analysis of external orientation elements; In view of the characteristics of linear array push-broom imaging satellite imaging, the external orientation elements refer to the three straight line elements used to describe the spatial coordinate values ​​of the photographic center at the moment of imaging a line of images and the three angular elements used to express the spatial attitude. The position of the photographic center can be indirectly obtained through the GPS module carried on the satellite, while the attitude of the satellite imaging can be obtained using the star sensor measurement data.

[0034] The orbital position error, i.e. the position error of the photographic center, can be decomposed into three directions: along the track, perpendicular to the track, and radial. According to the characteristics of the linear array push-scan imaging mode, the along-track error and the perpendicular to the track error are translation errors, while the image point offset caused by the radial error is a translation error and a proportional error. Based on the onboard high-precision orbit determination processing system, the proportional error caused by the radial error is small, so these three types of errors can be equivalent to the translation error.

[0035] The on-board attitude system error can also be divided into roll angle error, pitch angle error, and yaw angle error according to the rotation direction. The roll angle error is similar to the pitch angle error and can be equivalent to the translation error. The error caused by the yaw angle can be equivalent to the CCD rotation error.

[0036] During the process from satellite launch to on-orbit operation, camera installation, star sensors, etc. will produce errors due to the influence of the external environment, which is the main source of equipment installation errors. The translation error caused by displacement in this process can be equivalent to orbit deviation error, and the installation angle error can be equivalent to attitude measurement error.

[0037] (1-2-2) Error analysis of internal orientation elements; The internal orientation element errors mainly include the principal point offset error, the principal distance error of the satellite camera, the detector size error, the CCD rotation error, the error caused by the lens distortion, etc. The internal orientation element errors reflect the degree of internal distortion of the satellite.

[0038] The principal point offset error can be equivalent to the translation error along the track and perpendicular to the track. The principal distance error will cause the scaling of the image point coordinates, which can be equivalent to the proportional error along the track and perpendicular to the track. The detector size error causes the image point offset, which can be equivalent to the proportional error.

[0039] like Figure 3 As shown in the figure, the satellite is affected by external forces during its orbital operation, causing the CCD to be placed in a non-perpendicular position to the satellite's flight direction. At this time, the offsets in the x and y directions are: (2); Where: is the rotation angle, is the center of rotation, Indicates the y coordinate of the pixel.

[0040] The errors caused by the deformation of the camera lens surface or the displacement of its position are collectively referred to as errors caused by lens distortion. This distortion is mainly divided into radial distortion and eccentric distortion. Radial distortion is caused by the surface error of the lens, and eccentric distortion refers to the distortion caused by the lenses not being kept on a horizontal plane. The errors in the two directions caused by radial distortion are: (3); in: represents radial distortion, , and are polynomial fixed coefficients, are the variables in the polynomial, , like point Coordinates in the camera coordinate system , Represents the principal distance of the image point (camera principal distance).

[0041] The error caused by eccentricity distortion is: (4); in, , and are the coefficients of the equation to be solved.

[0042] (2) Calibration of internal and external orientation elements.

[0043] (2-1) External orientation element calibration; The purpose of external orientation element calibration is to eliminate the systematic errors in orbit, attitude and other data. The orbit measurement error can be treated equivalently with the attitude measurement error, and the load installation error is equivalent to the attitude measurement error. Therefore, in the process of external orientation element calibration, these three can be unified, and the bias matrix can be introduced to correct the imaging light in the imaging process so that the light points to the real position, so as to achieve the purpose of eliminating systematic errors. The bias matrix is ​​an orthogonal rotation matrix composed of the angle between the actual observation direction of the pixel in the body coordinate system and the theoretical observation direction. represents the bias matrix, Represents the rotation angle of the imaging light around the x-axis, y-axis, and z-axis of the body coordinate system , , , As shown below: (5); In order to obtain the true direction of the light in the body coordinate system, in actual application, remote sensing satellite images and control image data are matched to obtain control points, and the geographic coordinates of the control point data are used to inversely calculate the true direction of the light in the body coordinate system. The deviation between the light direction of the image point in the body coordinate system and the light line inversely calculated from the control point is calculated, and the light direction is corrected by the bias matrix. According to formula (1), the corrected form after adding the bias matrix is: (6); Formula (6) is the external orientation element compensation model.

[0044] (2-2) Internal orientation element calibration; Satellites can use three CCDs to stitch together for imaging. The internal orientation errors introduced by multiple CCDs mainly include the translation and rotation characteristic errors caused by the position deviation of the CCD device, as well as internal orientation element errors such as principal point offset error, principal distance error, and detector size error. These errors are highly correlated, and some errors cannot be modeled and eliminated. Therefore, it is difficult to establish a distortion model through the distortion caused by various errors to realize internal orientation element calibration.

[0045] In fact, the combination of various errors caused by the internal orientation elements can be regarded as the influence on the real light direction of each imaging probe, which causes the deviation of the real position of the image point. In practical applications, high-precision points with the same name as the control image and the satellite can be obtained as ground control points, and the control points can be used to restore the real direction of the pixel.

[0046] The relationship between the camera detector and the imaging center in the camera coordinate system is as follows: Figure 4 As shown, the imaging light of the detector According to the decomposition of perpendicular to the track direction and along the track direction, the angles of the imaging light relative to these two directions are respectively , then the image point Coordinates in the camera coordinate system It can be expressed by these two angles as: , Represents the principal distance of the image point (camera principal distance). At the same time, due to the existence of gross error points in the matching, the use of the pointing angle model to calibrate each pixel cannot take into account the local smoothness characteristics, which makes the detector pointing angle obtained by using the gross error point solution cause a large error. At the same time, the number of control points required to independently solve the pointing angle for each pixel is too large. Therefore, a polynomial is used to approximate the process of the pointing angle changing with the detector image sequence (the detector image sequence can be represented by the x-coordinate of the image point). Under normal circumstances, when the highest degree of the polynomial is 5, it can fully simulate the errors of various internal orientation elements. The pointing angle simulated by the polynomial method can be expressed as: (7); in, Indicates the x coordinate of the image point; At this point, the geometric model can be rewritten as: (8); Combining the above-mentioned internal and external orientation element models, the satellite joint calibration model is: (9).

[0047] (3) Analysis of satellite imaging geometric parameters and model construction; The construction of geometric positioning models for different satellites requires knowing the attitude and orbit information of the satellite during imaging. The imaging mode of linear array push-scan imaging satellites requires obtaining the time, attitude, and orbit corresponding to each line of images. Since the measurement frequency of the attitude and orbit measurement equipment on the satellite is limited, the orbit, attitude and other data obtained on the satellite are discrete, and generally, it is impossible to directly obtain the data information corresponding to each line. Therefore, it is necessary to perform error analysis on the attitude, orbit, and time characteristics of different satellites themselves, and use the most appropriate interpolation model according to the data characteristics, so that the model can effectively restore the external orientation element information of the satellite during imaging, so as to improve the positioning accuracy of the rigorous geometric model.

[0048] This application analyzes the corresponding imaging geometric parameter characteristics of satellite orbit, flight time, and attitude auxiliary data information, and constructs corresponding orbit, flight time, and attitude models.

[0049] (3-1) Track characteristics analysis and model construction In order to analyze the real situation of the data more intuitively, this application uses the real downlink data from the satellite for analysis. The satellite uses the global positioning system to record the position and speed information of the satellite at a certain moment in the WGS-84 system.

[0050] According to the satellite auxiliary data format, the corresponding orbital data information can be parsed. The original orbital data of the calibration scene image can be analyzed with time as the variable to obtain the scatter plot of the X, Y, and Z positions in the original satellite data with time as the variable under the WGS-84 system, and the X velocity under the WGS-84 system. , Y speed , Z speed Scatter plot with time as the variable. The orbital position of the satellite changes relatively steadily over time, without obvious fluctuations, and the satellite speed also tends to be stable over time. This shows that the orbital position of the satellite is relatively stable during its in-orbit operation. According to the smooth characteristics of its scatter plot, a polynomial model can be used to fit the orbital data, and the orbital information at any time can be obtained through this polynomial model.

[0051] The orbital polynomial fitting model is shown below: (10); Where: represents the polynomial coefficients, Represents the imaging moment.

[0052] (3-2) Analysis of travel time characteristics and model construction; The satellite records the time information on board, and then interpolates the time on board and transmits it to the ground.

[0053] The onboard time data obtained by analyzing the original onboard data is line-by-line time, so the image behavior independent variable can be used to analyze the change of each line imaging time with the line number.

[0054] The onboard imaging time shows a linear trend with the row number, and the row integration time is within 0.000148 seconds-0.000151 seconds. The row integration time is stable, and the maximum integration time jump is less than 0.000003 seconds. The linear interpolation method can be used to model the row time.

[0055] (3-3) Posture characteristic analysis and model construction.

[0056] (3-3-1) Posture characteristics analysis; The attitude of the satellite defines the rotation matrix between the Earth-fixed coordinate system and the orbital coordinate system. After the extended Kalman filter, it is expressed as the quaternion q in the J2000 coordinate system. During the on-orbit operation, the satellite calculates the quaternion value through the star sensor installed on the satellite and transmits it to the auxiliary data, so the initial auxiliary data is parsed to obtain the quaternion. The specific information is as follows Figure 5 However, in fact, quaternions have no geometric meaning, so it is impossible to analyze the satellite's on-board attitude by only observing the changes in quaternions. In order to better analyze the changes in the on-board attitude, the quaternions are converted into Euler angles with geometric meaning, and the changes in the on-board attitude are analyzed by observing the Euler angles.

[0057] For satellites with high-frequency attitude jitter (high-frequency attitude jitter satellites), the changes in the on-board attitude can also be analyzed by observing the Euler angles. The converted pitch, roll, and yaw angles are analyzed with the sampling time as the independent variable, and the following is obtained: Figure 6 Results shown.

[0058] Figure 6 In the a subgraph, b subgraph, c subgraph Figure 3The following figures show the scatter plots of pitch, roll, and yaw angles over time. All three figures have obvious overall upward or downward trends, and the changes of the three angles do not show smooth characteristics, but have obvious fluctuations, especially the roll angle and yaw angle, the attitude jitter is very obvious. The yaw angle has the highest oscillation frequency of the three angles, which is about 0.25Hz.

[0059] Commonly used posture models include polynomial models based on Euler angles. Figure 7 The figure shows the result of the polynomial model simulating the original data. Figure 7 It can be seen that the fitting effect of the third-order polynomial and the fourth-order polynomial on the original data is not ideal, and does not reflect the more obvious oscillation effect in the original data. Therefore, it is not advisable to use polynomials as the basis of the posture model.

[0060] analyze Figure 6 It can be found that the oscillation frequency of the yaw angle is the highest, which is about 0.25Hz, while the sampling frequency of the on-board attitude is 4Hz. If the correct attitude model is selected, the oscillation process of the attitude can be simulated more accurately. Therefore, this application proposes a penalty spline function model as an attitude model to simulate the on-board attitude for the attitude error characteristics of high-frequency attitude jitter satellites.

[0061] (3-3-2) Penalized spline function model; Common posture models cannot accurately simulate posture data with high posture oscillation frequency and large oscillation amplitude. In order to better solve this problem, the spline function model can be used to solve related problems.

[0062] The spline function is composed of polynomials on some subspaces with continuity conditions. Suppose there are n+1 points , these points satisfy , these points are called nodes. , specify an integer , on these nodes Subspline function The following two conditions are met: a) In each interval superior, Is a number polynomial of ; b) In superior have Continuous derivatives of order.

[0063] Common spline functions include cubic splines. Assume that there is a point , the cubic spline function is actually A spline function passes through these points and satisfies the above conditions.

[0064] According to the conditions of the cubic spline function, the following conditions can be obtained: (11); From the definition of the spline function, we know that in each small interval There are 4 coefficients to be determined, and there are n intervals in total, so a total of 4n coefficients need to be determined. According to the continuity condition of the second-order derivative, 3n-3 conditions can be determined, and the cubic spline function can determine n+1 conditions through all points, a total of 4n-2 conditions, so 2 conditions are needed to determine the spline function. . Usually in the interval Two conditions are added to the boundary to calculate all the unknown coefficient values. However, due to the noise in the measurement data, the jitter frequency of the posture is high, and the use of cubic spline function to simulate the posture data will result in overfitting.

[0065] Considering the existence of noise, we hope to use a loss function to better balance the relationship between noise and posture fitting accuracy. The form of the loss function is as follows: (12); In the formula represents the basis matrix of the spline function, is the parameter matrix of the spline function model, represents the undetermined coefficients of the spline function, a represents the posture record value, is the penalty factor. is a symmetric penalty positive definite matrix. The penalty method of Eilers and Marx is used to determine the symmetric penalty positive definite matrix by calculating the second-order derivative or difference of the basis function. , for example, if a spline function is used as the basis function, it can be constructed by the second-order derivative of the spline At the same time, we select B-spline functions as the basis functions of spline functions to determine , and the penalty factor It plays an important role in balancing noise and attitude fitting accuracy.

[0066] Wahba proposed a generalized cross-validation method to estimate , The estimated value is The minimum value of is defined as: (13); in, It is a (n+1)×(n+1) matrix. is defined as follows: (14); This method can be used to determine the penalty factor, thereby completing the construction of the penalty spline function.

[0067] It is understandable that the above example analyzes the change of attitude of high-frequency attitude jitter satellites during in-orbit operation. According to the on-board data, it is found that the satellite attitude jitter is obvious. Considering the relatively high stability of the on-board attitude measurement components, the random errors caused by the measurement are eliminated. For attitudes with more obvious jitter, it is difficult to simulate more accurately with commonly used attitude models such as polynomial models, but the on-board attitude sampling frequency is high, and the on-board attitude can be simulated by an accurate model. This application uses a penalty spline function as the basis for establishing an attitude model. The real data uploaded from the satellite are compared with the penalty spline model using different attitude models, and the accuracy of the attitude model is reflected according to the multi-spectral registration accuracy. Experiments have shown that the penalty spline function model has a good fitting effect, which can eliminate the intersection error of the same-name points with obvious trends between spectral bands, and ensure the positioning accuracy of this satellite product data.

[0068] (4) Solving the joint calibration model; In the joint calibration process, the matching data obtained include the control points matched by the panchromatic band image and the multispectral image with the control image (reference image used in geographic information system (GIS) or image registration); the same-name points matched between the multispectral and panchromatic; the same-name points between different panchromatic CCDs; the same-name points between different multispectral CCDs. Panchromatic band images (images captured within a wide band, usually covering the entire range of visible light) have higher resolution and can match high-precision control points, so the calibration effect is more ideal than multispectral images. However, the imaging time interval between multispectral and panchromatic band images is short, and high-precision same-name points can be matched between multispectral and panchromatic. At this time, adjusting the internal orientation elements of the multispectral through the same-name points between the panchromatic and multispectral will have a more obvious effect. For example, the satellite can generate five images of the same area in a short time interval, including panchromatic band images, red band images B1, green band images B2, blue band images B3, and infrared band images B4.

[0069] (4-1) Control point equation (standard equation fitted by correct error-free control points); The joint calibration model uses the bias matrix and polynomial model to eliminate the errors of internal and external orientation elements. The equation constructed by matching the panchromatic and multispectral band images with the control image is: (15); Where: Represent the three-dimensional coordinates of the control points that match the panchromatic band, multispectral red band, and the control image. , represents the pointing angle polynomial model coefficient. The variables in this formula are and polynomial coefficients. Let: .

[0070] Formula (15) can be simplified as: (16); in: , (17); According to formula (17), after eliminating m, we can set: (18); After linearizing the above equation, we get: (19); In the formula Respectively represent The above process is the control point solution process in the joint calibration model.

[0071] (4-2) Small intersection angle connection point equation (In the field of remote sensing, especially in satellite image processing, "small intersection angle" usually refers to a small angle between the lines of sight (LOS) of two sensors).

[0072] like Figure 8 As shown in the figure, if two satellites A and B image the same ground point, the image points are located at Department, is the imaging position corresponding to satellite A, is the imaging position corresponding to satellite B. At this time, the imaging equation of the two points can be expressed as: (20); Subscript and Used to distinguish The corresponding image points and The corresponding image point, and Indicates the pointing angle in the x and y directions.

[0073] If the imaging geometric parameters of the two stars are completely consistent with the actual satellite conditions and the elevation of the ground object point S is known, then the geographic coordinates of the two points calculated according to the imaging model should be the same. (twenty one); Usually, DEM (Digital Elevation Model) is used to obtain elevation. Due to the limited accuracy of DEM, the error caused by elevation is difficult to avoid. Even if the imaging geometric parameters of the two stars are completely correct, the geographic coordinates of the two points are still not equal. The error caused by elevation error can be approximately expressed as: (twenty two); Where: and is the attitude angle of the two previous and subsequent images; is the elevation error. According to formula (22), Depend on as well as and Control, if possible and Infinitely close, that is, when the satellite has very close attitude angles during two shots, the intersection error of the same-name points caused by the elevation can be eliminated. At this time, the positioning error of the same-name points can directly reflect the internal and external orientation element errors of the two bands during imaging.

[0074] The satellite's panchromatic and multispectral imaging CCD elements can be installed on the same camera and the distance between CCDs of different bands is very small. When the satellite is push-scanning, the time interval between imaging the same object in different bands is short. It can be considered that the attitude angles of the two shots are very close. At this time, the intersection error of the same-name point caused by the elevation error can be effectively controlled.

[0075] Since the imaging attitude angles of two points in different bands are close, it is necessary to introduce the earth ellipsoid equation when solving the three-dimensional coordinates of the same-name points. There are two common methods for solving the three-dimensional coordinates of the ground points where the same-name points intersect using the ellipsoid equation: iterative transformation based on DEM assistance and small intersection angle stereo intersection method.

[0076] (4-2-1) DEM-assisted iterative transformation; According to the earth ellipsoid model, the three-dimensional coordinates of any point on the ellipsoid surface The following relationship should be satisfied: (twenty three); Where: and Represent the major and minor semi-axes of the ellipsoid respectively; Represents the height of the ground point. For the WGS84 ellipsoid, m, m. It can be obtained through DEM, thus establishing and , The relationship between.

[0077] According to formula (20), it can be determined , , and relationship, this time (twenty four); In the formula All are known values, and the first formula of formula (20) can be transformed into: (25); Substituting this formula into the ellipsoid equation, we get (26); The only unknown quantity in this formula is , so by solving the quadratic equation we can get . At this time, Bringing it back to formula (25), we can solve for a three-dimensional coordinate .make = As the initial value of the elevation, Bring it into DEM and interpolate to get the corresponding elevation value , replace the formula (32) with Updated to Calculate the new , again calculate new elevation values ​​in the DEM , repeat the previous process until twice The difference is less than a certain threshold, then the calculated That is the three-dimensional coordinate corresponding to the point.

[0078] Using the above method, we can calculate , Two points ( and Used to distinguish the geographic coordinates corresponding to two image points), and the average value is the geographic coordinates of the ground point corresponding to the two points.

[0079] (4-2-2) Small intersection angle stereo intersection method; According to the principle of stereo intersection, when multiple satellite images intersect the ground with imaging rays at different angles, the three-dimensional coordinates of the ground points can be calculated without relying on DEM.

[0080] However, since the time interval between panchromatic and multispectral imaging is extremely short and they are located on the same camera, the attitude angles of the two images are very close when imaging the same-name points, which makes the elevation solution in three-dimensional coordinates unstable. Therefore, the earth ellipsoid model needs to be introduced. Through the intersection model, the three-dimensional coordinates of the same-name points are solved using the adjustment method.

[0081] The first equation in formula (26) (the equation on the left side of the first plus sign in formula (26)) is transformed into: (27); Will , are all known values, let , then: (28); Substituting the above formula into formula (33), we get: (29); Eliminate according to the above formula Then establish the equation, let . and Only , , is an unknown number, and according to the ellipsoid equation: (30); There are only two independent variables in this equation. , ,right and After expanding it according to the Taylor series and retaining the first-order term, we get: (31); in: ; According to formula (30), it can be calculated that: (32); At this time, according to formula (31), the error equation can be constructed to solve the three-dimensional coordinates of the same-name points.

[0082] (4-2-3) Solving the small intersection angle connection point of the joint calibration model; Combined with formula (26), the data of the same-name points of panchromatic and multispectral matching, the same-name points between different panchromatic CCDs, and the same-name points between different multispectral CCDs are combined. The bias matrix and pointing angle polynomial model are also used to eliminate the errors of the internal and external orientation elements. The equation can be rewritten as: (33); In the formula Represents the same-name points matched between the panchromatic band and a certain band in the multispectral band. Since the imaging attitude angles between the same-name points of different bands are approximately the same, the same-name point intersection error caused by elevation will be controlled within a certain range. The same-name point intersection error of the two points reflects the errors in the external and internal orientation elements of the two images. Therefore, the geographic coordinates of the same-name points of the panchromatic band and each band of the multispectral band can be calculated as the initial values. The three-dimensional coordinates of the same-name points are not control data and have errors themselves. Therefore, in addition to solving And polynomial coefficients, we also need to solve the geographic coordinates of the same-name points The amount of correction.

[0083] The same-name points matched in the CCD overlapping area of ​​each band image have the same characteristics as the same-name points between spectral bands, the imaging attitude angles are approximately the same, and the error equation construction is the same as that of the same-name points between spectral bands.

[0084] Transform formula (33) similarly to the control point equation, and let: (34); Where: , , … Represents the values ​​of the bias matrix.

[0085] After linearizing the above equation, we get: (35); In the formula Respectively represent This process is the process of solving the connection points of small intersection angles in the joint calibration model.

[0086] (4-3) Joint calibration parameter solution The error equation is constructed as follows: (36); : Error vector (or observation vector) that represents the difference between the observed values ​​and the values ​​predicted by the model. In remote sensing, this usually refers to the difference between the actual coordinates of the ground control points and the coordinates predicted by the geometric model. Each element of is an error term.

[0087] : The coefficient matrix (design matrix) contains the coefficients associated with the model parameters and is used to convert changes in the parameters into changes in the observed values. Each row of corresponds to an observation value, and each column corresponds to a model parameter. In the design matrix, different geometric relationships and physical models can be represented by different combinations of rows and columns.

[0088] : Parameter vector, which contains the model parameters that need to be estimated, such as the position of the sensor, attitude parameters, coordinates of the ground points, etc. Each element of is a model parameter.

[0089] : Observation vector (or constant term vector), containing the actual observation values, usually the coordinates of the ground control points. Each element of is an observation.

[0090] : Weight matrix, which is a diagonal matrix, the elements on the diagonal represent the precision or weight of each observation. In remote sensing data processing, different observations may have different reliability or precision, so their weights in the error equation are also different. Used in the least squares method to adjust for the effect of different observations on the solution.

[0091] In the formula The coefficient matrix is ​​expressed as: (37); in Represents the attitude angle coefficient matrix, which includes: (38); in Respectively represent the two equations established through each pair of observation data. According to formula (34), let: ; but The result of each item in the formula (39) is as follows:

[0092] Substituting the initial value into formula (39), the coefficient values ​​can be obtained. Since the CCDs of each panchromatic and multispectral band are located on the same camera, the same bias matrix is ​​used to correct the exterior orientation elements of each band.

[0093] In formula (37), Represent the panchromatic and multispectral band pointing angle coefficient matrices respectively. For example, The matrix is ​​specifically: (40); in Represent the two equations established by each pair of observation data. ; Substituting into formula (34) we get: (41); but After calculating each item of: , (42); The solution method is the same as same.

[0094] In formula (37), represents the connection point coefficient matrix, expressed as: (43); in represents the three-dimensional coordinates of the ground points corresponding to the nth pair of points with the same name. According to formula (43), let: Then we get: (44); According to the ellipsoid model, Can be , then the independent variable in the formula is only The coefficient matrix The value of is calculated as: (45); Where: ; in and Represent the major and minor semi-axes of the ellipsoid respectively; Represents the height of the ground point. For the WGS84 ellipsoid, m, m.

[0095] In formula (45), the coefficient matrix Corresponding adjustment parameters for: (46); in Represents the pointing angle model coefficient of the nth band, Represents the geographic coordinates corresponding to the nth pair of points with the same name.

[0096] In formula (45) represents the initial value calculated based on each observation value ( ), Represents the observation weight matrix.

[0097] According to the least squares principle, we can calculate: All the values ​​in the formula are known quantities. The initial values ​​obtained by calculation are continuously updated through iteration. When the value is less than a certain threshold, the result is the joint calibration parameter.

[0098] It is understandable that this application analyzes the error sources and characteristics of the internal and external orientation elements in the process of building the linear array push-scan imaging model, analyzes and establishes corresponding orbit, attitude, and travel time models based on the satellite's own orbit, attitude, and travel time characteristics, pays attention to the attitude jitter problem of the satellite due to the instability of the satellite platform, and proposes a method of establishing an attitude model based on a penalty spline function model to reduce the model error caused by insufficient model accuracy and ensure the positioning accuracy of a rigorous model. The corresponding geometric calibration model is constructed based on the error sources and error characteristics of the positioning accuracy, and the control data is used to improve the geometric positioning accuracy of the image. While improving the positioning accuracy, the consistency accuracy between multi-spectral images is taken into account. The high-precision same-name points between the panchromatic and multi-spectral data are used for joint calibration, and the geometric positioning model is corrected to ensure the positioning consistency of the geometric positioning model of each band image, providing precision support for the final landing application of the satellite.

[0099] The following is a description of the geometric positioning processing device based on the linear array push-broom imaging satellite provided in the present application. The geometric positioning processing device based on the linear array push-broom imaging satellite described below and the geometric positioning processing method based on the linear array push-broom imaging satellite described above can be referenced to each other.

[0100] The present application provides a geometric positioning processing device based on a linear array push-broom imaging satellite, comprising: a construction module and an acquisition module. Among them: A construction module is used to construct a bias matrix and a pointing angle polynomial model. The bias matrix is ​​used to correct the pointing direction of the imaging light of the detector. The pointing angle polynomial model is used to simulate the process in which the pointing angle changes with the detector image sequence. The pointing angle is the angle of the imaging light of the detector relative to the direction perpendicular to the track and the direction along the track. The acquisition module is used to adjust the strict geometric positioning model based on the bias matrix and the pointing angle polynomial model to obtain the joint calibration model. The strict geometric positioning model is used to characterize the geometric relationship between the image points and the ground points.

[0101] It can be understood that the detailed functional implementation of each of the above-mentioned units / modules can be found in the introduction of the aforementioned method embodiment, and will not be repeated here.

[0102] It should be understood that the above-mentioned device is used to execute the method in the above-mentioned embodiment. The implementation principle and technical effect of the corresponding program module in the device are similar to those described in the above-mentioned method. The working process of the device can refer to the corresponding process in the above-mentioned method, which will not be repeated here.

[0103] Based on the method in the above embodiment, the embodiment of the present application provides an electronic device, which may include: a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus. The processor may call the logic instructions in the memory to execute the method in the above embodiment.

[0104] In addition, the logic instructions in the above-mentioned memory can be implemented in the form of software functional units and can be stored in a computer-readable storage medium when sold or used as an independent product. Based on this understanding, the technical solution of the present application, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product, which is stored in a storage medium and includes a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present application.

[0105] Based on the method in the above embodiment, an embodiment of the present application provides a computer-readable storage medium, which stores a computer program. When the computer program runs on a processor, the processor executes the method in the above embodiment.

[0106] Based on the method in the above embodiment, an embodiment of the present application provides a computer program product. When the computer program product runs on a processor, the processor executes the method in the above embodiment.

[0107] It should be understood that the various numerical numbers involved in the embodiments of the present application are only used for the convenience of description and are not used to limit the scope of the embodiments of the present application.

[0108] It will be easily understood by those skilled in the art that the above description is only a preferred embodiment of the present application and is not intended to limit the present application. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present application shall be included in the scope of protection of the present application.

Claims

1. A geometric positioning processing method based on a linear array push-broom imaging satellite, characterized in that: include: Construct a bias matrix and a pointing angle polynomial model. The bias matrix is ​​used to correct the pointing direction of the imaging light of the detector. The pointing angle polynomial model is used to simulate the process in which the pointing angle changes with the detector image sequence. The pointing angle is the angle of the imaging light of the detector relative to the direction perpendicular to the track and along the track. Based on the bias matrix and the pointing angle polynomial model, the rigorous geometric positioning model is adjusted to obtain the joint calibration model. The rigorous geometric positioning model is used to characterize the geometric relationship between image points and ground points.

2. The geometric positioning processing method based on linear array push-broom imaging satellite according to claim 1 is characterized in that: The step of adjusting the strict geometric positioning model based on the bias matrix and the pointing angle polynomial model to obtain the joint calibration model includes obtaining the joint calibration model through the following formula: ； in, Indicates the coordinate value of the pixel point on the satellite image corresponding to the WGS84 coordinate system. express The position vector of the phase center of the GPS module in the WGS84 coordinate system at this moment, represents the proportionality coefficient, for The conversion matrix from the J2000 coordinate system to the WGS84 coordinate system at this moment, for The transformation matrix from the body coordinate system to the geocentric inertial coordinate system at the moment, represents the bias matrix, Represents the transformation matrix from the camera coordinate system to the body coordinate system; represents the pointing angle polynomial model, Represents pixel points coordinate, and Indicates the coefficients of the angular polynomial model, with the subscript Represents the degree of the polynomial.

3. The geometric positioning processing method based on linear array push-broom imaging satellite according to claim 2 is characterized in that: The bias matrix is ​​specifically: ； in, represents the rotation angle of the imaging light of the detector around the x-axis of the body coordinate system, represents the rotation angle of the imaging light of the detector around the y-axis of the body coordinate system, Represents the rotation angle of the imaging light of the detector around the z-axis of the body coordinate system.

4. The geometric positioning processing method based on linear array push-broom imaging satellite according to claim 3 is characterized in that: Also includes: Bring the control points into the joint calibration model and construct the control point equation: ； By solving the control point equation, we can obtain The amount of correction; in, Represents the three-dimensional coordinates of control points, which are obtained by matching the panchromatic band image to the control image or the control points are obtained by matching the multispectral image to the control image.

5. The geometric positioning processing method based on linear array push-broom imaging satellite according to claim 4 is characterized in that: Also includes: Bring the same-name points into the joint calibration model and construct the same-name point equation: ； By solving the equation of the same name point, we can obtain The amount of correction; in, Represents the three-dimensional coordinates of keypoints, which are obtained by matching the panchromatic band image to a band in the multispectral image.

6. The geometric positioning processing method based on linear array push-broom imaging satellite according to any one of claims 1 to 5, characterized in that: Also includes: Based on the spline function and the following loss function, a penalty spline function model is constructed: ； The posture data is simulated through a penalized spline function model; Where W represents the parameter matrix of the spline function model, B represents the basis matrix of the spline function, c represents the undetermined coefficient of the spline function, a represents the posture record value, represents the penalty factor, and P represents the penalty positive definite matrix.

7. A geometric positioning processing device based on a linear array push-broom imaging satellite, characterized in that: include: A construction module is used to construct a bias matrix and a pointing angle polynomial model. The bias matrix is ​​used to correct the pointing direction of the imaging light of the detector. The pointing angle polynomial model is used to simulate the process in which the pointing angle changes with the detector image sequence. The pointing angle is the angle of the imaging light of the detector relative to the direction perpendicular to the track and the direction along the track. The acquisition module is used to adjust the strict geometric positioning model based on the bias matrix and the pointing angle polynomial model to obtain the joint calibration model. The strict geometric positioning model is used to characterize the geometric relationship between the image points and the ground points.

8. An electronic device, characterized in that: include: at least one memory for storing a computer program; At least one processor is used to execute the program stored in the memory. When the program stored in the memory is executed, the processor is used to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium storing a computer program, characterized in that: When the computer program runs on a processor, the processor is caused to execute the method according to any one of claims 1 to 6.

10. A computer program product, characterized in that When the computer program product runs on a processor, the processor is caused to execute the method according to any one of claims 1 to 6.

Images