CN111862332A - Method and system for correcting fitting error of satellite image general imaging model - Google Patents
Method and system for correcting fitting error of satellite image general imaging model Download PDFInfo
- Publication number
- CN111862332A CN111862332A CN202010748246.2A CN202010748246A CN111862332A CN 111862332 A CN111862332 A CN 111862332A CN 202010748246 A CN202010748246 A CN 202010748246A CN 111862332 A CN111862332 A CN 111862332A
- Authority
- CN
- China
- Prior art keywords
- grid
- imaging
- correction
- image
- image point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 173
- 238000000034 method Methods 0.000 title claims abstract description 30
- 238000012937 correction Methods 0.000 claims abstract description 139
- 238000005457 optimization Methods 0.000 claims abstract description 12
- 230000009466 transformation Effects 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 15
- 239000011159 matrix material Substances 0.000 claims description 11
- 238000000528 statistical test Methods 0.000 claims description 9
- 238000007689 inspection Methods 0.000 abstract description 3
- 238000012545 processing Methods 0.000 description 5
- 238000012892 rational function Methods 0.000 description 4
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000004590 computer program Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000009897 systematic effect Effects 0.000 description 3
- 206010044565 Tremor Diseases 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Software Systems (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Physics (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Image Processing (AREA)
Abstract
The invention provides a method and a system for correcting fitting errors of a general imaging model of satellite images, wherein a three-dimensional space grid based on a strict geometric imaging model is established; performing re-imaging correction grid establishment and image point grid coordinate correction, wherein the re-imaging correction grid establishment and the image point grid coordinate correction are performed, the re-imaging correction grid is set, the attitude polynomial corresponding to an ideal image is established, the re-imaging correction grid is set, the initial value is calculated by attitude transformation re-imaging, and the correction value corresponding to the image point grid is calculated by interpolation of the re-imaging correction grid; carrying out fitting residual statistical inspection and re-imaging correction grid optimization, wherein the fitting residual statistical inspection and re-imaging correction grid optimization are carried out, an object space three-dimensional space grid and a corrected image point grid are utilized, an RPC parameter is calculated by adopting a least square method, image point coordinates corresponding to the re-imaging correction grid are calculated by the RPC parameter, a fitting residual of a general imaging model is calculated by utilizing the re-imaging correction grid and the image point coordinates, the correction quantity of the re-imaging correction grid is calculated by utilizing the fitting residual corresponding to the grid, and the optimized re-imaging correction grid is realized.
Description
Technical Field
The invention belongs to the technical field of remote sensing, and particularly relates to compensation and correction of fitting errors of a satellite image rational polynomial universal imaging model.
Background
The satellite image imaging model expresses the mathematical relationship between the image point coordinates and the object space coordinates of the corresponding target point, and is the basis of the geometric processing of the satellite image. The satellite image geometric imaging Model can be divided into a strict geometric imaging Model (RSM) and a general Rational polynomial imaging Model (RFM). Due to the advantages of simple form, strong universality, independence on sensor structure and the like, the universal imaging model based on rational function has become the most widely used imaging model and is commonly used by satellite image providers, and almost all satellite images provide RPC (rational polymeric coeffients) parameters of the rational function model. Meanwhile, almost all satellite photogrammetry processing software supports a rational function model to perform satellite image geometric processing.
RPC parameters of the rational function model are generally obtained by fitting according to a strict geometric model of a satellite image by a terrain-independent method, and the establishment of the strict geometric imaging model according to the structure and imaging characteristics of a sensor is a premise of RPC modeling. The rigorous geometric imaging model is generally based on a collinear equation and is correspondingly expanded according to the structure and the imaging characteristics of the sensor. The satellite platform is generally subjected to tremor of different degrees under the influence of external space environment and internal mechanical operation, so that the satellite attitude has periodic oscillation phenomenon. The ultrahigh resolution sensor is generally carried on an agile satellite platform to acquire a strip image in an asynchronous push-broom imaging mode, and because the attitude of the sensor changes dynamically in the imaging process, the high-frequency attitude disturbance is more obvious. Meanwhile, with the improvement of resolution, the imaging integration time is shorter and shorter, so that the image deformation is more and more sensitive to the attitude disturbance of the platform. General imaging models based on rational polynomials are generally difficult to express local image deformation caused by high-frequency tremor, so that stereoscopic images have systematic vertical parallax (y-paralaxes), and furthermore, DSMs produced by the images have obvious fringe effect, and practice shows that the periodic elevation system error is not acceptable in application. How to effectively eliminate high-frequency disturbance in the sensor attitude and image deformation caused by the high-frequency disturbance to establish a high-precision general imaging model and calculate corresponding RPC parameters has important value for high-precision positioning and application of high-resolution satellite images.
In order to ensure higher RPC modeling accuracy, generally speaking, for a long strip image, dividing the long strip image into standard scenes (scenes) according to a scene division rule (such as World Reference System), and then fitting RPC parameters to the standard scenes. In the local image range of the standard shot, the RPC parameters can be considered to be equivalent to the strict model parameters under the condition that the image length is not large, the orbit and the posture are smooth, and the fitting error is small (0.01 pixel level). The standard shot mode is generally used for medium and high resolution satellite images (e.g., SPOT, IRS P5, Ziyuan-3, RapidEye). However, sub-level agile satellites such as QuickBird, WorldView, etc. generate standard level 1 video products by using a segmentation (segment) mode, and since the video length is large, in order to distribute and store, a stripe is divided into a plurality of sub-blocks (tiles) and each sub-block corresponds to a different RPC, which causes application inconvenience:
(1) corresponding image points of the same object space coordinate in the overlapped area are not completely the same;
(2) the requirement of image positioning on ground control points is increased;
(3) the three-dimensional model is complex in structure;
(4) the image coordinate conversion is troublesome, and particularly when the ground coordinate is projected to the image coordinate, it is necessary to determine which sub-block the ground point is located in by iterative calculation, and the corresponding RPC parameter is used.
In order to realize generalized and convenient photogrammetric processing, it is necessary to perform integral RPC modeling on the image of the long strip and generate RPC of the whole segmented image, so that a relatively strict geometric imaging model of the RFM has a large fitting error, and the positioning accuracy and the application of the general imaging model are affected.
Aiming at the problems, the invention provides an ideal re-imaging grid correction technical scheme for the fitting error of a general imaging model of satellite images.
Disclosure of Invention
The invention provides a re-imaging grid correction technical scheme for fitting errors of a satellite image general imaging model, which aims to effectively eliminate high-frequency disturbance in the attitude of a sensor and image deformation caused by the high-frequency disturbance, establish a high-precision general imaging model, calculate corresponding RPC parameters and simultaneously realize generalized and convenient photogrammetry processing.
The technical scheme provided by the invention provides a method for correcting the fitting error of a general imaging model of satellite images, which comprises the following steps,
step 1, establishing a three-dimensional space grid based on a strict geometric imaging model;
step 2, re-imaging correction grid establishment and image point grid coordinate correction are carried out, and the method comprises the following substeps,
step 2.1, fitting the attitude metadata by a cubic polynomial, and establishing an attitude polynomial corresponding to the ideal image;
step 2.2, setting a re-imaging correction grid according to the attitude metadata time interval and the image width;
step 2.3, re-imaging by utilizing posture transformation, and calculating an initial value of a re-imaging correction grid;
step 2.4, interpolating and calculating a correction value corresponding to the image point grid by using the re-imaging correction grid to generate a corrected image point grid;
step 3, carrying out fitting residual statistical test and re-imaging correction grid optimization, comprising the following substeps,
step 3.1, calculating RPC parameters by using an object space three-dimensional space grid and the corrected image point grid and adopting a least square method,
step 3.2, calculating the image point coordinates corresponding to the re-imaging correction grid by using the RPC parameters;
and 3.3, calculating the fitting residual error of the general imaging model by using the re-imaging correction grid and the image point coordinates.
And 3.4, calculating the correction quantity of the re-imaging correction grid according to the fitting residual errors corresponding to the grid points, and realizing the optimization of the re-imaging correction grid.
Furthermore, step 2.1 is implemented in such a way that,
given sensor pose metadataThe unit quaternion is converted to the corresponding euler angle omega,
kappa, then calculating a cubic fitting polynomial corresponding to the time series of the three corner elements,
wherein, ω (l),
kappa (l) is the Euler angle corresponding to any scan line l, l0The coefficient a of the cubic polynomial is the scan line corresponding to the image centerm，bm，cmM is more than or equal to 0 and less than or equal to 3 through least square adjustment calculation.
Furthermore, step 2.2 is implemented in such a way that,
according to the time interval of the attitude metadata, the grid number nr of the row direction of the re-imaging correction grid is calculated as follows,
nr＝floor(rows×TDI/dt)+1
wherein, floor () represents rounding down, rows is the total number of lines of the image, TDI is the integral time corresponding to each scanning line, and dt is the time interval corresponding to the attitude metadata;
the grid spacing of the re-imaging correction grid is calculated according to the following formula,
wherein ceil () represents rounding up, cols is the total column number of the image, dr is the row direction interval of the re-imaging correction grid, dc is the column direction interval of the re-imaging correction grid, and the grid number nc in the column direction adopts a preset value;
calculating the coordinates r of the starting point of the correction grid according to the following formula0，c0，
Dot coordinates(s) of grid pointsi2，lj2) The following calculation is carried out,
where i2, j2 is the grid index in the re-imaging correction grid.
Furthermore, the step 2.3 is implemented by performing the following calculations on all grid points of the re-imaging correction grid to obtain the correction values of the grid points,
let the rotation matrix of the sensor coordinate system corresponding to the scanning line l relative to the orbit reference coordinate system be Rs(l) Forming an instantaneous rotation matrix according to the Euler angles calculated by the cubic polynomial in the step 2.1
The initial value of the re-imaging correction grid is calculated according to the following image point coordinate re-imaging transformation formula,
wherein, (s, l) is the image point coordinate of the grid point on the original image, x0，y0F is an image internal orientation element, delta l, delta s are re-imaging correction quantity of grid points (s, l), and lambda is a proportionality coefficient.
Furthermore, step 2.4 is implemented in such a way that,
for any image point coordinate (s, l), using the grid parameters obtained in step 2.2, the index (i) in the re-imaging correction grid where the image point is located is calculated according to the following formula0，j0)，
And calculating the re-imaging correction quantity of the image point coordinates by bilinear interpolation according to the correction value corresponding to the grid.
Furthermore, step 3.2 is implemented in such a way that,
calculating the longitude and latitude of an object space corresponding to the grid point by a strict geometric imaging model according to the image point coordinates and the object grid elevation corresponding to the re-imaging correction grid, and generating an object space grid corresponding to the re-imaging correction grid;
then, using the RPC parameters obtained in step 3.1, calculating coordinates of the corresponding image points of the object space grid points on the ideal image, and recording as(s)1，l1) And obtaining an ideal image point grid.
Furthermore, step 3.3 is implemented in such a way that,
calculating a correction value corresponding to the image point coordinate by using the reimaging correction grid in a corresponding mode of the step 2.4 and through bilinear interpolation, calculating the image point coordinate residual error of the RFM imaging model according to the following formula,
wherein i, j is the index of the image point grid, k is the index of the object space grid elevation layering, (s, l) is the image point coordinate of the grid point, (s1，l1) The coordinates of the image points calculated by the RFM imaging model in step 3.2 are (delta s, delta l) correction values corresponding to the re-imaging correction grid, vs(i，j，k)，vl(i, j, k) is the fitting residual of the generic imaging model to the rigorous geometric model.
Furthermore, step 3.4 is implemented in such a way that,
and (4) performing statistical test on the fitting residual errors obtained in the step (3.3), and if the size of the residual errors exceeds a limit difference or the residual error distribution has systematicness, taking the average value of the residual errors of a plurality of different elevations corresponding to the grid points as the correction quantity of the fitting errors of the grid points to realize the optimization of the correction grid.
The invention also provides a system for correcting the fitting error of the satellite image general imaging model, which is used for realizing the method for correcting the fitting error of the satellite image general imaging model.
Compared with the prior art, the invention has the advantages and beneficial effects that:
the invention provides a re-imaging grid correction scheme of satellite image general rational polynomial imaging model fitting errors, a virtual image is established through a posture angle (Euler angle) of a cubic polynomial fitting sensor under a satellite orbit reference coordinate system, the conversion of image point coordinates between the virtual image and an original image is realized through a re-imaging algorithm based on posture transformation, the fitting error correction is concluded to be the re-imaging grid correction, the statistical test and the systematic correction are further carried out on the fitting errors of the corrected grid, the complete compensation and the correction of the general imaging model fitting strict geometric imaging model errors are realized, and the direct positioning accuracy of a general geometric model can be effectively improved.
Drawings
Fig. 1 is a flowchart of establishing a three-dimensional spatial grid based on a strict geometric imaging model according to an embodiment of the present invention.
FIG. 2 is a flow chart illustrating the process of re-imaging correction mesh creation and image point mesh coordinate correction according to an embodiment of the present invention.
FIG. 3 is a flowchart of fitting error statistical test and re-imaging correction grid optimization according to an embodiment of the present invention.
Detailed Description
The technical solution of the present invention is specifically described below with reference to the accompanying drawings and examples.
How to effectively eliminate high-frequency disturbance in the sensor attitude and image deformation caused by the high-frequency disturbance to establish a high-precision general imaging model and calculate corresponding RPC parameters has important value for high-precision positioning and application of high-resolution satellite images. Meanwhile, in order to realize generalized and convenient photogrammetry processing, it is necessary to perform integral RPC modeling on a satellite image with a long bandwidth breadth and generate uniform RPC parameters of the whole segmented image, so that a relatively strict geometric imaging model of the RFM has a large fitting error, and the positioning accuracy and the application of the general imaging model are influenced.
The invention provides an ideal re-imaging grid correction method for improving the geometric positioning precision of a satellite image by correcting and compensating the fitting error of a general rational polynomial imaging model of a high-resolution satellite. The method comprises the steps of establishing a virtual image through an attitude angle (Euler angle) of a cubic polynomial fitting sensor under a satellite orbit reference coordinate system, realizing the conversion of image point coordinates between the virtual image and an original image through a re-imaging algorithm based on attitude transformation, converting fitting error correction into re-imaging grid correction, further performing statistical inspection and systematic correction on the fitting error of the corrected grid, realizing complete compensation and correction of the fitting strict geometric imaging model error of a universal imaging model, and effectively improving the direct positioning precision of the universal geometric model.
The embodiment of the invention provides a method for correcting fitting errors of a general imaging model of satellite images, which comprises the following steps: 1. three-dimensional space grid establishment based on strict geometric imaging model
This step is a prerequisite for the subsequent steps and can be implemented with reference to the prior art. For ease of reference, an implementation description of the embodiments is provided. Referring to fig. 1, the three-dimensional spatial grid creation is divided into the following 5 sub-steps, which are described below.
1.1 establishing a strict geometric imaging model of the original satellite image by using the metadata.
The earth center rectangular coordinate system is selected as an object space coordinate system, the projection center is used as the origin of a satellite orbit reference coordinate system, the position vector of the earth center rectangular coordinate system under the object space coordinate system is used as a z-axis, the x-axis is positioned in a plane formed by a satellite flight motion vector and the z-axis and points to the direction of the motion vector, the y-axis is perpendicular to an xz plane, and the direction of the y-axis is determined by a right-hand spiral rule. The imaging geometry corresponding to any scan line in the image can be expressed as a collinear condition equation of the instantaneous projection center, the image point and the corresponding target point. The general form of the ideal geometric imaging model of the satellite image can be expressed as:
where (s, l) is the coordinates of the image point, (x)0，y0) Is the principal point-like coordinate, and f is the equivalent principal distance. [ X, Y, Z ]]Is the coordinate of the target point under the rectangular coordinate system of the earth center, [ X ]s(l)，Ys(l)，Zs(l)]Is scanning line l pairsAnd the corresponding projection center instantaneous position is in coordinates under a geocentric rectangular coordinate system. Ro(l) The instantaneous attitude matrix of the satellite orbit coordinate system corresponding to the scanning line l under the earth center rectangular coordinate system can be obtained by calculating the instantaneous position vector of the satellite and the corresponding velocity vector (the first derivative of the orbit equation). Rs(l) The attitude matrix of the sensor corresponding to the scanning line l in the orbit reference coordinate system expresses the rotation between the sensor coordinate system and the satellite orbit coordinate system, such as the pitch angles of the front-view camera and the back-view camera, the yaw angle of the side-view camera, and the like. λ is a scaling factor.
1.2 calculating the ground coverage range and the ground elevation range of the image according to the image point coordinates corresponding to the four corners of the image.
The ground elevation of a given target point is 0m, namely the point is positioned on the earth reference ellipsoid, the coordinates of the earth center rectangular coordinate system of the target point meet the rotating ellipsoid equation,
wherein a and b are the major and minor axes of the reference ellipsoid. The equation (1) is substituted into the formula to obtain a quadratic equation of a first order about the proportionality coefficient lambda, so that the proportionality coefficient lambda can be solved, and then the equation (1) is substituted to obtain the ground three-dimensional coordinates of the target point and convert the ground three-dimensional coordinates into longitude and latitude. And respectively calculating the ground longitude and latitude corresponding to the four corner points of the image to obtain the minimum and maximum longitude and latitude corresponding to the image coverage range.
According to the image coverage range, further counting the minimum and maximum elevation values of the grid points of the digital elevation model in the quadrilateral range formed by the minimum and maximum longitude and latitude, and recording the minimum and maximum elevation values as hmin，hmax。
1.3 determining the elevation range and the elevation layering of the three-dimensional space grid.
Given an elevation extension value deltah, e.g. 500m, the minimum maximum elevation obtained in step 1.2 is extended,
considering the actual range of the ground elevation, the value of the expanded elevation range is limited to Hmin≥-500m，Hmax10000m or less, thereby obtaining the minimum maximum height H after expansionmin，Hmax。
According to the elevation range determined above, n h elevation layers (nh is more than 3) are set at equal intervals, the elevation corresponding to each elevation layer is calculated according to the following formula,
Hi＝Hmin+k×(Hmax-Hmin)/nh，0≤k≤nh (4)
1.4 wherein k is the number of the layers. A grid is set on the original image at certain intervals, and the image point coordinates corresponding to the grid are calculated.
An image point grid is set on the 'original image' with the interval of lines and columns being dl and ds respectively, the image point coordinates corresponding to the grid points are calculated according to the following formula,
where i1, j1 is the grid index,(s)i1，lj1) And nl and ns are the grid numbers of the image point grids in the row and column directions respectively.
1.5 calculating object coordinates corresponding to the image point grid points through a strict geometric imaging model, and establishing a three-dimensional space grid.
According to the ephemeris metadata, satellite orbit ephemeris data given in Chebyshev polynomial fitting metadata are used, orbit polynomials in three directions under a geocentric rectangular coordinate system can be obtained respectively, the instantaneous projection center position corresponding to any scanning line can be obtained by utilizing polynomial calculation, and the corresponding speed vector can be obtained by utilizing the first derivative of the polynomial calculation.
In the embodiment, the unit quaternion is used for representing the sensor attitude metadata, and the unit quaternion corresponding to the instantaneous attitude at any time t can be obtained by calculation through a spherical linear interpolation (SLERP) method.
Qt＝(sin((1-t)×θ)×Q1+sin(t×θ)×Q2)/sin(θ) (6)
Wherein θ is a quaternion Q1And Q2The included angle between the two is calculated according to the following formula,
O＝arccos(Q1·Q2) (7)
therefore, the rotation matrix R of the corresponding instantaneous attitude of any scanning line l can be conveniently calculateds(l)。
When the elevation of the object space point is HkWhen the earth center three-dimensional rectangular coordinate is obtained, the corresponding earth center three-dimensional rectangular coordinate meets the following ellipsoid equation,
wherein a and b are the major and minor axes of the reference ellipsoid. Substituting equation (1) into the above formula can obtain a quadratic equation of a first order about the proportionality coefficient lambda, so as to solve the proportionality coefficient lambda, and then substituting equation (1) can calculate the ground three-dimensional coordinates of the target point and convert the ground three-dimensional coordinates into longitude and latitude. And repeating the calculation on all grid points to obtain the object space three-dimensional space grid corresponding to the image point grid.
2-time imaging correction grid establishment and image point grid coordinate correction
Referring to fig. 2, the present invention proposes to perform re-imaging correction grid creation and image point grid correction, including 4 sub-steps, which are respectively described as follows.
2.1 fitting the attitude metadata with a cubic polynomial to establish an attitude polynomial corresponding to the 'ideal image' with steady change of attitude.
First, a unit quaternion given by the sensor attitude metadata is converted into a corresponding euler angle omega,
kappa, then calculating a cubic fitting polynomial corresponding to the time series of the three corner elements,
wherein, ω (l)，
Kappa (l) is the Euler angle corresponding to any scan line l, l0The coefficient a of the cubic polynomial is the scan line corresponding to the image centerm，bm，cm(m is more than or equal to 0 and less than or equal to 3) can be obtained by least square adjustment calculation.
After the euler angle is calculated through the attitude polynomial, an attitude matrix corresponding to the ideal image can be further calculated, namely the 'ideal attitude matrix'.
And 2.2, setting a re-imaging correction grid at certain intervals according to the attitude metadata time interval and the image width. Wherein the re-imaging correction grid is a grid which can realize re-imaging correction through interpolation calculation.
According to the time interval of the attitude metadata, the number nr of grids in the row direction of the re-imaging correction grid is calculated according to the following formula,
nr＝floor(rows×TDI/dt)+1 (10)
wherein, floor () represents rounding down, rows is the total number of lines of the image, TDI is the integration time corresponding to each scan line, and dt is the time interval corresponding to the gesture metadata. The number nc of the grids in the column direction may be a predetermined value, and may be generally set to about 15, and may be appropriately adjusted according to the image width.
Further, the grid interval of the re-imaging correction grid is calculated according to the following formula,
where ceil () denotes rounding up, cols is the total number of columns of the image, dr is the column direction interval of the re-imaging correction grid, and dc is the column direction interval of the re-imaging correction grid.
Then calculating the coordinates r of the starting point of the correction grid according to the following formula0，c0，
Dot coordinates(s) of grid pointsi2，lj2) The following calculation is carried out,
where i2, j2 is the grid index, here the grid index in the re-imaging correction grid.
And 2.3, calculating an initial value of the re-imaging correction grid by using a posture transformation re-imaging method.
The initial value of the re-imaging correction grid is calculated according to the following formula of 'image point coordinate re-imaging transformation',
wherein, (s, l) is the image point coordinate of the grid point on the original image, x0，y0F is an in-image orientation element, and (x) is the same as (1)0，y0) Is the coordinate of the principal point, and f is the equivalent principal distance; Δ l, Δ s is the amount of re-imaging correction of the grid dots (s, l). Interpolating the unit quaternion corresponding to the scanning line l according to the spherical linear interpolation (SLERP) method in step 1.5, and further calculating to obtain a rotation matrix R of the corresponding sensor coordinate system relative to the orbit reference coordinate systems(l) Instantaneous rotation matrix formed by Euler angles calculated according to the cubic polynomial described in step 2.1
Repeating the calculation on all grid points of the heavy imaging correction grid to obtain the correction values of all grid points.
And 2.4, interpolating and calculating a correction value corresponding to the image point grid by using the re-imaging correction grid to generate the corrected image point grid.
For any image point coordinate (s, l), using the grid parameters given in step 2.2, the index (i) of the "re-imaging correction grid" in which the image point is located is calculated according to the following formula0，j0)，
And calculating the re-imaging correction quantity of the image point coordinates by bilinear interpolation according to the correction value corresponding to the grid. Repeating the calculation for all points of the image point grid to obtain the image point grid coordinates after re-imaging correction.
3 fitting residual error statistical test and re-imaging correction grid optimization
Referring to fig. 3, the fitting residual statistical test and re-imaging correction grid optimization includes 4 sub-steps, which are described below.
3.1, calculating RPC parameters by using the object space three-dimensional space grid and the corrected image point grid and adopting a least square method.
This step can be implemented with reference to the prior art. For ease of reference, the implementation of the embodiments is provided as follows:
a general expression of the satellite imagery generic imaging model RFM is as follows,
where (S, L) is the normalized image point coordinates, (E, N, H) is the normalized ground coordinates, Num (E, N, H) and Den (E, N, H) are cubic polynomials, Nums(E，N，H)、Dens(E，N，H)、Numl(E，N，H)、Denl(E, N, H) are respectively used for identifying polynomials in column direction and row direction, the order of each variable in the polynomials does not exceed 3, and the sum of the orders of all the variables does not exceed 3.
The normalized translation (offset) and scaling (scale) parameters are calculated using the maximum (max _ value) and minimum (min _ value) values of the image point coordinate grid and the virtual space grid, respectively, as follows,
then, the RPC parameters of the rational polynomial imaging model are calculated by least squares adjustment through the error equation set forth in formula (16).
And 3.2, calculating the coordinates of the image points corresponding to the re-imaging correction grid by using the RPC parameters.
And (3) calculating the longitude and latitude of an object corresponding to the grid point by using the image point coordinates and the object grid elevation corresponding to the re-imaging correction grid according to the mode in the step 1.7 through a strict geometric imaging model, and generating the object space grid corresponding to the re-imaging correction grid.
Then, using the RPC parameters obtained in step 3.1, calculating the coordinates of the corresponding image points of the object space grid points on the ideal image according to the formula (15), and recording the coordinates as(s)1，l1) And obtaining the ideal image point grid.
And 3.3, calculating the fitting residual error of the universal imaging model by using the re-imaging correction grid and the image point coordinates.
Calculating a correction value corresponding to the image point coordinates by using the reimaging correction grid and bilinear interpolation in the mode of the step 2.4, further calculating the image point coordinate residual error of the RFM imaging model according to the following formula,
wherein i, j is the index of the image point grid, k is the index of the object space grid elevation layering, (s, l) is the image point coordinate of the grid point, (s1，l1) The coordinates of the image points calculated by the RFM imaging model in step 3.2 are (delta s, delta l) correction values corresponding to the re-imaging correction grid, vs(i，j，k)，vl(i, j, k) is the fitting residual of the generic imaging model to the rigorous geometric model.
And 3.4, calculating the correction quantity of the re-imaging correction grid according to the fitting residual errors corresponding to the grid points.
And (3) performing statistical test on the fitting residual errors obtained in the step (3.3), if the size of the residual errors exceeds a limit difference or the residual error distribution has systematicness, taking the average value of the residual errors of nh different elevations corresponding to the grid points as the correction quantity of the fitting errors of the grid points, realizing the optimization of the correction grid, and optimizing the re-imaging correction grid according to the following formula.
Wherein Δ s (i, j) and Δ l (i, j) are correction amounts of the fitting error of the grid points.
In specific implementation, a person skilled in the art can implement the automatic operation process by using a computer software technology, and a system device for operating the method, such as a computer-readable storage medium storing a corresponding computer program according to the technical solution of the present invention and a computer device including a corresponding computer program for operating the corresponding computer program, should also be within the scope of the present invention.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.
Claims (9)
1. A method for correcting fitting errors of a general imaging model of satellite images is characterized by comprising the following steps: comprises the following steps of (a) carrying out,
step 1, establishing a three-dimensional space grid based on a strict geometric imaging model;
step 2, re-imaging correction grid establishment and image point grid coordinate correction are carried out, and the method comprises the following substeps,
step 2.1, fitting the attitude metadata by a cubic polynomial, and establishing an attitude polynomial corresponding to the ideal image;
step 2.2, setting a re-imaging correction grid according to the attitude metadata time interval and the image width;
step 2.3, re-imaging by utilizing posture transformation, and calculating an initial value of a re-imaging correction grid;
step 2.4, interpolating and calculating a correction value corresponding to the image point grid by using the re-imaging correction grid to generate a corrected image point grid;
step 3, carrying out fitting residual statistical test and re-imaging correction grid optimization, comprising the following substeps,
step 3.1, calculating RPC parameters by using an object space three-dimensional space grid and the corrected image point grid and adopting a least square method,
step 3.2, calculating the image point coordinates corresponding to the re-imaging correction grid by using the RPC parameters;
3.3, calculating the fitting residual error of the general imaging model by using the re-imaging correction grid and the image point coordinates;
and 3.4, calculating the correction quantity of the re-imaging correction grid according to the fitting residual errors corresponding to the grid points, and realizing the optimization of the re-imaging correction grid.
2. The method for correcting fitting errors of the satellite image general imaging model according to claim 1, wherein: the step 2.1 is implemented in such a way that,
the unit quaternion given by the sensor attitude metadata is converted into a corresponding euler angle omega,
kappa, then calculating a cubic fitting polynomial corresponding to the time series of the three corner elements,
wherein, ω (l),
kappa (l) is the Euler angle corresponding to any scan line l, l0The coefficient a of the cubic polynomial is the scan line corresponding to the image centerm,bm,cmM is more than or equal to 0 and less than or equal to 3 through least square adjustment calculation.
3. The method for correcting fitting errors of the satellite image general imaging model according to claim 2, wherein: the step 2.2 is implemented in such a way that,
according to the time interval of the attitude metadata, the grid number nr of the row direction of the re-imaging correction grid is calculated as follows,
nr＝floor(rows×TDI/dt)+1
wherein, floor () represents rounding down, rows is the total number of lines of the image, TDI is the integral time corresponding to each scanning line, and dt is the time interval corresponding to the attitude metadata;
the grid spacing of the re-imaging correction grid is calculated according to the following formula,
wherein ceil () represents rounding up, cols is the total column number of the image, dr is the row direction interval of the re-imaging correction grid, dc is the column direction interval of the re-imaging correction grid, and the grid number nc in the column direction adopts a preset value;
calculating the coordinates r of the starting point of the correction grid according to the following formula0,c0，
Dot coordinates(s) of grid pointsi2,lj2) The following calculation is carried out,
where i2, j2 is the grid index in the re-imaging correction grid.
4. The method for correcting fitting errors of the satellite image general imaging model according to claim 3, wherein: the step 2.3 is realized by calculating all grid points of the heavy imaging correction grid to obtain the correction values of the grid points,
let the rotation matrix of the sensor coordinate system corresponding to the scanning line l relative to the orbit reference coordinate system be Rs(l) Forming an instantaneous rotation matrix according to the Euler angles calculated by the cubic polynomial in the step 2.1
The initial value of the re-imaging correction grid is calculated according to the following image point coordinate re-imaging transformation formula,
wherein, (s, l) is the image point coordinate of the grid point on the original image, x0,y0F is an image internal orientation element, delta l, delta s are re-imaging correction quantity of grid points (s, l), and lambda is a proportionality coefficient.
5. The method for correcting fitting errors of the satellite image general imaging model according to claim 4, wherein: the step 2.4 is implemented in such a way that,
for any image point coordinate (s, l), using the grid parameters obtained in step 2.2, the index (i) in the re-imaging correction grid where the image point is located is calculated according to the following formula0,j0)，
And calculating the re-imaging correction quantity of the image point coordinates by bilinear interpolation according to the correction value corresponding to the grid.
6. The method for correcting fitting errors of the satellite image general imaging model according to claim 5, wherein: the step 3.2 is implemented in such a way that,
calculating the longitude and latitude of an object space corresponding to the grid point by a strict geometric imaging model according to the image point coordinates and the object grid elevation corresponding to the re-imaging correction grid, and generating an object space grid corresponding to the re-imaging correction grid;
then, using the RPC parameters obtained in step 3.1, calculating coordinates of the corresponding image points of the object space grid points on the ideal image, and recording as(s)1,l1) And obtaining an ideal image point grid.
7. The method for correcting fitting errors of the satellite image general imaging model according to claim 6, wherein: the step 3.3 is implemented in such a way that,
calculating a correction value corresponding to the image point coordinate by using the reimaging correction grid in a corresponding mode of the step 2.4 and through bilinear interpolation, calculating the image point coordinate residual error of the RFM imaging model according to the following formula,
wherein i, j is the index of the image point grid, k is the index of the object space grid elevation layering, (s, l) is the image point coordinate of the grid point, (s1,l1) The coordinates of the image points calculated by the RFM imaging model in step 3.2 are (delta s, delta l) correction values corresponding to the re-imaging correction grid, vs(i,j,k),vl(i, j, k) is the fitting residual of the generic imaging model to the rigorous geometric model.
8. The method for correcting fitting errors of the satellite image general imaging model according to claim 7, wherein: the implementation of step 3.4 is that,
and (4) performing statistical test on the fitting residual errors obtained in the step (3.3), and if the size of the residual errors exceeds a limit difference or the residual error distribution has systematicness, taking the average value of the residual errors of a plurality of different elevations corresponding to the grid points as the correction quantity of the fitting errors of the grid points to realize the optimization of the correction grid.
9. A correction system for fitting errors of a general imaging model of satellite images is characterized in that: the correction method for realizing the fitting error of the satellite image universal imaging model as claimed in claims 1 to 8.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010748246.2A CN111862332B (en) | 2020-07-30 | 2020-07-30 | Correction method and system for fitting errors of general imaging model of satellite image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010748246.2A CN111862332B (en) | 2020-07-30 | 2020-07-30 | Correction method and system for fitting errors of general imaging model of satellite image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111862332A true CN111862332A (en) | 2020-10-30 |
| CN111862332B CN111862332B (en) | 2024-06-18 |
Family
ID=72946414
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010748246.2A Active CN111862332B (en) | 2020-07-30 | 2020-07-30 | Correction method and system for fitting errors of general imaging model of satellite image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111862332B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113222871A (en) * | 2021-05-31 | 2021-08-06 | 哈尔滨工程大学 | Multi-view satellite image digital surface model fusion method |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101295747B1 (en) * | 2012-04-25 | 2013-08-20 | 서울시립대학교 산학협력단 | System of automatic geometric correction using rational polynomial cofficient and method thereof |
| CN103310487A (en) * | 2013-06-21 | 2013-09-18 | 中国科学院遥感与数字地球研究所 | Generating method for universal time variable based imaging geometric model |
| CN105091906A (en) * | 2015-06-30 | 2015-11-25 | 武汉大学 | High-resolution optical push-broom satellite steady-state reimaging sensor calibration method and system |
| CN107689064A (en) * | 2017-08-08 | 2018-02-13 | 武汉大学 | Take the strict geometry imaging model construction method of satellite optical of aberration correction into account |
| CN110500995A (en) * | 2019-07-12 | 2019-11-26 | 武汉大学 | The method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter |
-
2020
- 2020-07-30 CN CN202010748246.2A patent/CN111862332B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101295747B1 (en) * | 2012-04-25 | 2013-08-20 | 서울시립대학교 산학협력단 | System of automatic geometric correction using rational polynomial cofficient and method thereof |
| CN103310487A (en) * | 2013-06-21 | 2013-09-18 | 中国科学院遥感与数字地球研究所 | Generating method for universal time variable based imaging geometric model |
| CN105091906A (en) * | 2015-06-30 | 2015-11-25 | 武汉大学 | High-resolution optical push-broom satellite steady-state reimaging sensor calibration method and system |
| CN107689064A (en) * | 2017-08-08 | 2018-02-13 | 武汉大学 | Take the strict geometry imaging model construction method of satellite optical of aberration correction into account |
| CN110500995A (en) * | 2019-07-12 | 2019-11-26 | 武汉大学 | The method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter |
Non-Patent Citations (2)
| Title |
|---|
| WANG, TAOYANG 等: "Multi-Mode GF-3 Satellite Image Geometric Accuracy Verification Using the RPC Model", SENSORS * |
| 李海鸿;曹辉;施俊;: "高分辨率光学卫星影像高精度定位技术与实践", 地理空间信息, no. 05 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113222871A (en) * | 2021-05-31 | 2021-08-06 | 哈尔滨工程大学 | Multi-view satellite image digital surface model fusion method |
| CN113222871B (en) * | 2021-05-31 | 2022-05-20 | 哈尔滨工程大学 | Multi-view satellite image digital surface model fusion method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111862332B (en) | 2024-06-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109903352B (en) | Method for making large-area seamless orthoimage of satellite remote sensing image | |
| CN110500995B (en) | Method for establishing high-resolution satellite image equivalent geometric imaging model by using RPC parameters | |
| CN107705329B (en) | High-resolution optical satellite staring image registration method based on geometric constraint | |
| Oh et al. | A piecewise approach to epipolar resampling of pushbroom satellite images based on RPC | |
| JP3869814B2 (en) | Orthorectification processing method for satellite images | |
| CN107341778B (en) | SAR image orthorectification method based on satellite control point library and DEM | |
| CN111174753B (en) | Optical image and laser height measurement data adjustment method based on rational function model | |
| CN106570905B (en) | A kind of noncooperative target point cloud initial attitude verification method | |
| CN109709551B (en) | Area network plane adjustment method for satellite-borne synthetic aperture radar image | |
| CN110030978B (en) | Method and system for constructing geometric imaging model of full-link optical satellite | |
| CN111724465B (en) | Satellite image adjustment method and device based on plane constraint optimization virtual control point | |
| CN107040695B (en) | satellite-borne video image stabilization method and system based on RPC positioning model | |
| CN108226982B (en) | Single linear array satellite laser combined high-precision positioning processing method | |
| CN115439571A (en) | Method and device suitable for generating linear array push-broom satellite image epipolar image | |
| CN112270698A (en) | Non-rigid geometric registration method based on nearest curved surface | |
| CN110660099B (en) | Rational function model fitting method for remote sensing image processing based on neural network | |
| KR102159134B1 (en) | Method and system for generating real-time high resolution orthogonal map for non-survey using unmanned aerial vehicle | |
| CN112070891A (en) | Image area network adjustment method and system with digital ground model as three-dimensional control | |
| CN111508028A (en) | Autonomous in-orbit geometric calibration method and system for optical stereo mapping satellite camera | |
| CN103278140B (en) | Coordinate back calculation method for TDICCD (time delay and integration charge coupled devices) linear array push-sweep sensor | |
| CN111862332A (en) | Method and system for correcting fitting error of satellite image general imaging model | |
| CN104537614B (en) | CCD image orthorectification method for environment satellite I | |
| CN113610924B (en) | High-precision coordinate back-calculation method, equipment and medium for linear approximation considering terrain effect | |
| CN115326025A (en) | Binocular image measuring and predicting method for sea waves | |
| CN111044076B (en) | Geometric calibration method for high-resolution first-number B satellite based on reference base map |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20240519 Address after: Room 1510, Floor 1-21, Building 7, Contemporary Science and Technology Park (Huaxia Entrepreneurship Center), No. 20 Guannanyuan 1st Road, Donghu New Technology Development Zone, Wuhan City, Hubei Province, 430070 Applicant after: Wuhan duopuyun Technology Co.,Ltd. Country or region after: China Address before: 7th Floor, Beidou Building, No. 980 Gaoxin Avenue, Donghu New Technology Development Zone, Wuhan City, Hubei Province, 430014 Applicant before: Wuhan Lvzhu Linjia Technology Co.,Ltd. Country or region before: China |
|
| TA01 | Transfer of patent application right | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |