Method for repairing RPC of remote-sensing image

Abstract

Disclosed in the present invention is a method for repairing an RPC of a remote-sensing image. The method comprises: S1, acquiring information such as an image to be repaired and a reference image; S2, performing feature extraction and matching on the image to be repaired and the reference image by using an SIFT algorithm, so as to obtain a matching point pair and corresponding first pixel coordinates; S3, acquiring longitude and latitude coordinates; S4, acquiring geodetic coordinates and a control point; S5, constructing an RFM model, performing RPC positive calculation on an original RPC parameter of the image to be repaired and the geodetic coordinates, so as to obtain corresponding second pixel coordinates of the geodetic coordinates on a remote-sensing image, and separately calculating a mean value of error distances between corresponding pixel coordinates; S6, introducing a constraint condition by using an L1 regularization optimization method, so as to obtain an optimized RFM model; and S7, performing determination according to the mean value of the error distances and a preset threshold value, and selecting different repair policies to recalculate or correct the RPC parameter. In the present invention, an RFM model is optimized by using an L1 regularization method, and different repair policies are selected to perform correction or estimation on an RPC parameter, thereby improving the geometric positioning precision of a satellite image.

Description

A remote sensing image RPC repair method Technical Field

[0001] The present invention relates to the field of remote sensing processing, and in particular to a remote sensing image RPC repair method. Background Art

[0002] Geometric imaging models for high-resolution satellite imagery are generally divided into two categories: rigorous imaging models and general imaging models. All parameters in rigorous imaging models have clear geometric meanings, resulting in the highest positioning accuracy. However, establishing rigorous imaging models requires the use of satellite orbital positions, attitude observations, and sensor structural parameters, which lacks universality and is therefore rarely used as a form of satellite imagery data release. As a mathematical analogue of the rigorous imaging model, the general imaging model (RFM) based on cubic rational polynomials offers advantages such as sensor parameter confidentiality, strong versatility, a simple form, and fast and convenient computation. Consequently, it has been widely studied and applied. Currently, the RPC (Rational Polynomial Coefficient) parameter of the RFM model has become a standard for high-resolution satellite imagery. Since the general model contains 78 rational polynomial coefficients (RPCs), it requires a sufficient and well-distributed array of control points for iterative solution and calculation.

[0003] Satellite images with RPC parameters have typically undergone geometric preprocessing, including on-orbit geometric calibration and sensor system error correction, resulting in a certain degree of direct ground positioning accuracy. These images can be used for subsequent processing such as geometric positioning, orthophoto rectification, and stereo mapping. However, RFM still has certain flaws in practical applications, primarily due to inherent RPC issues and external influences. First, RPC suffers from inherent problems such as over-parameterization, strong correlation between parameters, and residual errors from fitting a rigorous imaging model. Second, due to technical limitations such as sensor technology, many satellite platforms lack stability and have limited instantaneous attitude measurement accuracy. These external factors lead to significant systematic errors in satellite imagery and their RPC parameters, such as errors caused by high-frequency vibrations of the satellite platform.

[0004] These problems will directly lead to the difficulty in eliminating the vertical parallax of satellite images in applications such as stereo mapping and image matching, making stereo observation difficult and reducing positioning accuracy.

[0005] Summary of the Invention

[0006] Based on the above technical problems, the present invention provides a remote sensing image RPC repair method, which uses the L1 regularization method to optimize the RFM model, adopts different RPC parameter repair strategies, corrects or estimates the original RPC parameters, and controls the complexity of the RFM model, reduces overfitting, eliminates the problems of reduced geometric positioning accuracy caused by model errors in RFM, and thus improves the geometric positioning accuracy of high-resolution satellite images.

[0007] To achieve the above technical objectives, the present invention provides a remote sensing image RPC repair method, comprising the following steps:

[0008] S1 obtains the image to be restored, reference image, geographic information data and DEM data;

[0009] S2 uses the SIFT algorithm to extract and match features between the image to be repaired and the reference image to obtain matching point pairs and the corresponding first pixel coordinates;

[0010] S3 performs geographic coordinate conversion on the first pixel coordinate of the matching point pair according to the geographic information data to obtain the latitude and longitude coordinates of the matching point pair;

[0011] S4 interpolates the longitude and latitude coordinates of the matching point pairs according to the DEM data to obtain geodetic coordinates, and uses the matching point pairs corresponding to the geodetic coordinates as control points;

[0012] S5 builds an RFM model, performs RPC forward calculation on the original RPC parameters of the image to be repaired and the geodetic coordinates, obtains the second pixel coordinates corresponding to the geodetic coordinates on the remote sensing image, and calculates the mean of the error distances of the corresponding first pixel coordinates and the second pixel coordinates respectively;

[0013] S6 uses the L1 regularization optimization method based on the least squares adjustment principle and introduces the set linear constraints to obtain the optimized RFM model;

[0014] S7 makes a judgment based on the mean value of the error distance and the preset threshold, selects different repair strategies, and recalculates or corrects the RPC parameters.

[0015] Specifically, S51 constructs the RFM model:

[0016] Among them, P i (i=1,2,...,4) is a third-order polynomial, (l,s) is the normalized second pixel coordinate, (X,Y,Z) is the normalized geodetic coordinate of the control point, and all values ​​are between [-1,1];

[0017] S52 expands (l, s) into two linear equations according to the Taylor series: P1(X, Y, Z)-lP2(X, Y, Z)=0 P3(X, Y, Z)-sP4(X, Y, Z)=0;

[0018] S53 uses the geodetic coordinates of the control points to construct the two linear equations into a system of equations and solves the equations: y = Ax0 + e

[0019] Where y∈R m×1 is the observation vector, A∈R m×n is the coefficient matrix and e∈R m×1 is the residual vector, x0∈R n×1 is the RPCs vector; m = 2k, where k is the number of geodetic coordinates of the control points, and n is 78;

[0020] Solve the system of equations: x0=N -1 A T y

[0021] Where N = A T A is a normal matrix, x0∈R n×1 is the RPCs vector;

[0022] Based on the RFM model, S54 performs RPC forward calculation on the original RPC parameters of the remote sensing image and the geodetic coordinates of the control point to obtain the second pixel coordinates of the geodetic coordinates of the control point on the remote sensing image.

[0023] S55 calculates the average of the error distances between the first pixel coordinate and the second pixel coordinate.

[0024] Specifically, the preset linear constraint condition in step S6 is:

[0025] Where ||x||1 is the sum of the absolute values ​​of the vector x, and γ is the regularization parameter.

[0026] Specifically, the different repair strategies selected in step S7 are:

[0027] When the mean value of the error distance meets the first judgment condition, the first restoration strategy is selected, and the RPC is recalculated according to the original RPC parameters of the image to be restored and the optimized RFM model to obtain the restored RPC parameters;

[0028] When the mean value of the error distance meets the second judgment condition, the second restoration strategy is selected, and RPC correction is performed according to the geodetic coordinates, the original RPC parameters of the image to be restored, and the optimized RFM model to obtain the restored RPC parameters.

[0029] Specifically, the first repair strategy in step S7 includes:

[0030] In the optimized RFM model, the preset linear constraints are updated as follows: min x,e ||x||1+||e||1

[0031] Under the updated linear constraints, according to the geodetic coordinates of the control points, the equation system is constructed as follows: y=Ax+e

[0032] Solve the system of equations: x = N -1 A T y

[0033] Where x is the repaired RPC parameter.

[0034] Specifically, the second repair strategy in step S7 includes:

[0035] Use the original RPC parameter as the initial value x initial ;

[0036] In the optimized RFM model, the preset linear constraints are updated as follows: min Δx,e ||Δx||1+||e||1;

[0037] Under the updated linear constraints, according to the geodetic coordinates of the control points, the equation system is constructed as follows: y=Ax initial +AΔx+e;

[0038] Solve the system of equations: Δx=N -1 A T y-Ax initial

[0039] Among them, Δx is the optimal solution of the offset transformation;

[0040] At the initial value x initial On the basis of Δx, we can get x=x initial +Δx, where x is the repaired RPC parameter.

[0041] Specifically, step S2 includes:

[0042] S21 determines a frame of image in the same space as a reference image based on the remote sensing image, and performs image enhancement on the two images so that the two images have interactive characteristics in grayscale;

[0043] S22 performs block matching on the two images using the constraints of the spatial geographic information of the remote sensing image, divides the remote sensing image into a plurality of sub-area image blocks to be matched, and obtains the corresponding reference sub-area image block according to the projection coordinates of each image block;

[0044] S23 uses the Shi-Tomasi algorithm to move the regional search window in all directions, observes the average value of the window's energy expansion, and then determines the central pixel point where the window energy change is greater than a specified threshold as a corner point;

[0045] S24 uses the SIFT descriptor with scale and rotation invariance to establish a feature descriptor to describe the features of the detected feature points;

[0046] After constructing the SIFT descriptor, S25 performs rough feature matching based on the KD tree nearest neighbor query algorithm, uses the KD tree for shortest distance search, and eliminates cluttered feature points that have no homonyms.

[0047] S26 uses a bidirectional matching strategy to initially match homonymous points, and then uses the global relaxation method to eliminate the remaining homonymous points whose reliability cannot be determined;

[0048] S27 uses random sampling consistency algorithm and least squares iteration method to eliminate false matching pairs in the candidate matching point set;

[0049] S28 uses a greedy algorithm to perform homogenization processing on the feature points to obtain a final uniformly distributed matching point pair, and obtains the first pixel coordinates corresponding to each matching point in the remote sensing image.

[0050] Specifically, the error distance in step S4 is calculated using any one of Euclidean distance, similarity, and angle cosine.

[0051] The beneficial effects of the present invention include at least:

[0052] (1) The present invention uses the L1 regularization method to optimize the RFM model. According to the preset constraints, the complexity of the RFM model is controlled to reduce overfitting, eliminate the problem of reduced geometric positioning accuracy caused by model errors in the RFM, and improve the geometric positioning accuracy of high-resolution satellite images.

[0053] (2) The present invention selects different RPC parameter repair strategies based on the error distance of the original RPC of the calculated remote sensing image, corrects or estimates the original RPC parameters, and further improves the geometric positioning accuracy of high-resolution satellite images.

[0054] (3) The present invention adopts an automatic control point extraction method based on geographic constraint block processing, Shi_Tomasi algorithm corner point detection, triangulated network model, small facet correction and SIFT algorithm. Based on the geographic constraint block strategy, the method improves the matching accuracy by combining multiple matching strategies to obtain evenly distributed control points. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0056] FIG1 is a schematic diagram of a method according to an embodiment of the present invention;

[0057] FIG2 is a technical flow chart of an embodiment of the present invention. DETAILED DESCRIPTION

[0058] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. It should be noted that as long as no conflict arises, the various embodiments and features of each embodiment of the present invention can be combined with each other, and the resulting technical solutions are all within the scope of protection of the present invention.

[0059] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or combinations thereof.

[0060] Referring to FIG. 1 and FIG. 2 , in an embodiment of the present invention, a remote sensing image RPC repair method is provided, which is characterized by comprising the following steps:

[0061] S1 obtains the image to be restored, reference image, geographic information data and DEM data;

[0062] A DEM (Digital Elevation Model) is a dataset of plane coordinates (X, Y) and elevations (Z) of regularly gridded points within a specific area. It primarily describes the spatial distribution of regional landforms. It is generated by collecting data (including sampling and measurement) using contour lines or similar stereo models, followed by data interpolation. A DEM is a virtual representation of landforms, from which information such as contour lines and slope maps can be derived.

[0063] S2 uses the SIFT algorithm to extract and match features between the image to be repaired and the reference image to obtain matching point pairs and the corresponding first pixel coordinates;

[0064] In an embodiment of the present invention, step S2 includes:

[0065] S21 determines a frame of image in the same space as a reference image based on the remote sensing image, and performs image enhancement on the two images so that the two images have interactive characteristics in grayscale;

[0066] In the embodiment of the present invention, the two images are enhanced by linearly stretching the enhancement coefficient.

[0067] S22 performs block matching on the two images using the constraints of the spatial geographic information of the remote sensing image, divides the remote sensing image into a plurality of sub-area image blocks to be matched, and obtains the corresponding reference sub-area image block according to the projection coordinates of each image block;

[0068] In an embodiment of the present invention, the remote sensing image is subjected to standard block segmentation, each image is divided into M×N pixels, the four corners of the sub-area image are respectively denoted as A1, A2, A3 and A4, and the image coordinates of the upper left boundary point A1 and the lower right boundary point A3 are respectively calculated and recorded as (I1, J1) and (I3, J3);

[0069] Convert the image coordinates into the corresponding geographic coordinates, expressed as follows:

[0070] Among them, (X i ,Y i ) represents the geographic coordinates, (X0, Y0) represents the geographic coordinates of the upper left corner of the input image, (I i ,J i ) represents the image coordinates of the i-th point, and a, b, c, d represent the transformation model parameters.

[0071] Through the mapping relationship between image coordinates and geographic coordinates, the geographic coordinates (X1, Y1) of the upper left boundary point A1 and the geographic coordinates (X2, Y2) of the lower right boundary point A3 of the sub-area image block to be matched are obtained;

[0072] The upper left boundary point B1, lower right boundary point B3 and the corresponding image coordinates (I'1, J'1) and (I'3, J'3) of the reference sub-area image block are obtained by inverse calculation of the geographic coordinates. According to the image coordinates of the upper left boundary point B1 and the lower right boundary point B3, the image coordinates of the upper right boundary point B2 and the lower left boundary point B4 of the reference sub-area image block are obtained.

[0073] Using the constraints of spatial geographic information of remote sensing images for block matching can effectively reduce redundant control points and improve matching accuracy.

[0074] Based on the corresponding sub-area image blocks to be matched and the reference sub-area image blocks, the Shi_Tomasi algorithm is used to extract feature points from each image block. First, the Shi_Tomasi algorithm is used to move the regional search window in various directions, and the energy expansion average of the window is observed. Then, the central pixel point whose window energy change is greater than the specified threshold is determined as the corner point, which is corner point detection.

[0075] S23 uses the Shi-Tomasi algorithm to move the regional search window in all directions, observes the average value of the window's energy expansion, and then determines the central pixel point where the window energy change is greater than a specified threshold as a corner point;

[0076] In an embodiment of the present invention, step S23 is specifically as follows:

[0077] Among them, (x, y) represents the coordinates of the window center, (u, v) represents the window translation, I represents the image grayscale, w(x, y) represents the window, and the window w(x, y) is represented by a Gaussian function, which is expressed as follows:

[0078] Here, σ represents the standard deviation.

[0079] For a small translation amount, the Taylor expansion of the grayscale change after the window translation can be approximately expressed as a formula containing the autocorrelation matrix, which is expressed as follows:

[0080] Where M represents the autocorrelation matrix, which is expressed as follows:

[0081] Among them, I x ,I y Represents the gradient value of the image grayscale in the x and y directions respectively.

[0082] Calculate the eigenvalues ​​λ1 and λ2 of the autocorrelation matrix M and compare the minimum value to the set threshold T. If min(λ1,λ2)>T, retain the feature point and it is considered a Shi-Tomasi corner point. This algorithm extracts relatively stable feature points and can effectively extract points even when the image is rotated, the illumination angle changes, or there is noise in the image. Furthermore, this algorithm can reduce the aggregation of control points, avoiding redundant control points and making the distribution of feature points more reasonable and uniform.

[0083] S24 uses the SIFT descriptor with scale and rotation invariance to establish a feature descriptor to describe the features of the detected feature points;

[0084] In an embodiment of the present invention, step S24 is specifically as follows:

[0085] Utilize the distribution characteristics of the gradient direction to calculate each direction parameter specified by the pixels in the neighborhood of the extreme point. Calculate the gradient value m(x, y) and the gradient direction θ(x, y) of the feature point in the scale image L(x, y) according to the following two formulas respectively, which are expressed as follows:

[0086] Determine the main direction of the feature point within the 8×8 region centered on the feature point by using the gradient value and the gradient direction. Rotate the coordinate axes of the image, then determine a 4×4 region centered on each key point, and set 8 directions in each region. By calculating the gradient histogram of each sub-region, form the best SIFT descriptor, that is, 4×4×8 feature vectors. Finally, normalize the feature vectors to remove the influence of illumination changes.

[0087] After constructing the SIFT descriptor, perform coarse feature matching based on the K-D tree nearest neighbor query algorithm. Use the K-D tree to search for the shortest distance, and eliminate the cluttered feature points without homologous points.

[0088] In the embodiment of the present invention, it is assumed that there are N data points in the K-dimensional space. At the i-th dimension where the data changes the most, the data space is equally divided into two parts. The truncated data m is the mean of the k / 2-th and (k + 1) / 2-th data in the i-th dimension. According to the construction order of the original k-d tree, compare the i-th dimension of the search data with the truncated data m of the i-th dimension of the k-d tree, so that the number of data falling at both ends of m differs by no more than 1. Divide the search points into child nodes according to size until there is only one data participating in the comparison. However, some of the shortest distance points are divided by adjacent other data points. After the k-d tree search is completed, it is still necessary to search the neighboring nodes M times (M < N) to determine the position where the minimum distance is located.

[0089] Use the two-way matching strategy to perform initial matching of homologous points, and then use the global relaxation method to eliminate the remaining homologous points whose reliability cannot be determined.

[0090] In the embodiment of the present invention, the steps of the two-way matching strategy include:

[0091] Take a feature point (x r , y r ) in the sub-region to be matched, and find the two points with the closest Euclidean distance corresponding to the reference sub-region. Among these two feature points, when the closest distance DN and the second closest distance DSN satisfy, then this point is the homologous point (xm, ym) of the image to be matched; preferably in this embodiment, the Threshold threshold T is set to 0.8; take the homologous point (xm, ym) in the sub-region to be matched as the target point, and obtain the homologous point (x' r , y' r ) on the remote sensing image; Compare (xr ,y r ) and (x' r ,y' r ) is less than one pixel; if it is greater, it is eliminated. After checking the bidirectional matching, there are still some points with the same name whose reliability cannot be determined. These points are screened using the global relaxation matching method to combine the matching points with the points in their neighborhood.

[0092] In an embodiment of the present invention, a global relaxation method is used to eliminate the remaining points of the same name whose reliability cannot be determined, and the context information in the neighborhood is utilized. The constraints and consistency between objects are considered, and the most consistent and compatible results are finally obtained through iterative calculation.

[0093] The specific steps include: determining the conjugate point i of a certain point j in the remote sensing image in the reference image, treating the target point j as a category, and treating the conjugate candidate point as the target A. i , consider image matching as classification, that is, overall fuzzy classification, and the category set is: C = {C1, C2, ..., C m ,C0}

[0094] Among them, Cj represents the image point to be matched, j = 1, 2, ..., m, C0 is not other point types. Using the basic matching algorithm, the target set determines r points i1, i2, ..., i for each image point j to be matched. r As the conjugate candidate point of point j, the target set is: A={A 11 ,A 12 ,....,A 21 ,A 22 ,.....,A i1 ,A i2 ,...,A m1 ,A m2}, where target A i1 ,A i2 ,....,A ir The correlation coefficients with remote sensing image point j are ρ1, ρ2, ....ρ r , then A ik ∈C j The probability is:

[0095] For each image point ja to be matched and another point jb within a certain range, as well as their corresponding conjugate candidate points h and k in the reference image, the others can be deduced by analogy, and Phk can also be determined. Pij can be modified iteratively and gradually optimized:

[0096] S27 uses random sampling consistency algorithm and least squares iteration method to eliminate false matching pairs in the candidate matching point set;

[0097] In an embodiment of the present invention, the RANSAC algorithm is used to initially estimate the initial homography matrix H between the input remote sensing image and the reference image, and then the transformed coordinates of each to-be-matched point after rough matching are estimated based on H, which is expressed as follows:

[0098] Among them, (x w ,y w ) represents the coordinates of the remote sensing image feature points, (X, Y) represents (x w ,y w )’s transformed coordinates.

[0099] Calculate the root mean square error of the i-th matching point, which is expressed as follows:

[0100] When the residual and RMS error of a matching point pair is greater than the set threshold, the matching point pair is removed. The residual and RMS error are iteratively calculated using the remaining matching point pairs until the residual and RMS error of all matching point pairs are less than the desired threshold.

[0101] S28 uses a greedy algorithm to perform homogenization processing on the feature points to obtain a final uniformly distributed matching point pair, and obtains the first pixel coordinates corresponding to each matching point in the remote sensing image.

[0102] Refit the obtained matching point pairs, calculate the residual and RMSE of each pair of control points, and sort the RMSE of the control points in ascending order. Select the control point with i=1 in the sequence as the first node and number it. Use Euclidean distance to calculate the distance between this point and all other control points. If the Euclidean distance is less than the threshold, it is considered that the two control points are too densely distributed and the control point is removed. Set i=i+1. Preferably, the distance threshold for feature point homogenization is 1000. Repeat the above steps until all control points are traversed and stop. Output the homogenized matching point pairs to obtain the first pixel coordinate (x1, y1) corresponding to the remote sensing image.

[0103] The first pixel coordinates of the matching points on the remote sensing image are converted into longitude and latitude coordinates (B, L) according to the geographic parameter information.

[0104] In an embodiment of the present invention, the orthophoto coordinate conversion parameters are obtained based on the orthophoto geographic information. The parameters are generally affine transformation matrices that can realize the conversion between the pixel coordinates and the latitude and longitude coordinates of the image, as shown in the following formula: B = Ax1 + By1 + CL = Dx1 + Ey1 + F

[0105] Calculate the longitude and latitude coordinates corresponding to the matching point according to the above formula. In the formula, (x1, y1) is the first pixel coordinate, (B, L) is the longitude and latitude coordinates, Then, the elevation corresponding to the feature points is calculated by interpolating the DEM using raster data interpolation methods such as bilinear interpolation and bicubic convolution interpolation.

[0106] S3 performs geographic coordinate conversion on the first pixel coordinate of the matching point pair according to the geographic information data to obtain the latitude and longitude coordinates of the matching point pair;

[0107] S4 interpolates the longitude and latitude coordinates of the matching point pairs according to the DEM data to obtain geodetic coordinates, and uses the matching point pairs corresponding to the geodetic coordinates as control points;

[0108] In an embodiment of the present invention, the elevation corresponding to the point is obtained based on the acquired latitude and longitude coordinates and DEM data using a bilinear polynomial interpolation method. Based on the four nearest data points, a bilinear polynomial can be determined:

[0109] Where Z is the elevation of the data point, a 00 、a 01 、a 10 、a 11 The linear coefficients between the elevation H and the longitude and latitude coordinates (B, L) are calculated using the four known data points. The elevation of the point is then interpolated from the longitude and latitude coordinates (B, L) and the calculated coefficients. The geodetic coordinates (B, L, H) of the control point are obtained through the above conversion steps.

[0110] The RFM model is a universal imaging geometry model for satellite remote sensing imagery. It is applicable to a wide range of sensors, including the latest aerospace and space sensor models. RFM-based imaging models do not require understanding the actual sensor characteristics or imaging process. They are simple, generalized imaging models that achieve approximately the same accuracy as strict imaging models.

[0111] S5 builds an RFM model, performs RPC forward calculation on the original RPC parameters of the image to be repaired and the geodetic coordinates, obtains the second pixel coordinates corresponding to the geodetic coordinates on the remote sensing image, and calculates the mean of the error distances of the corresponding first pixel coordinates and the second pixel coordinates respectively;

[0112] In an embodiment of the present invention, specifically, S5 includes the following steps:

[0113] S51 builds the RFM model:

[0114] Among them, P i(i=1,2,...,4) is a third-order polynomial, (l,s) is the normalized second pixel coordinate, (X,Y,Z) is the normalized geodetic coordinate of the control point, and all values ​​are between [-1,1];

[0115] S52 expands (l, s) into two linear equations according to the Taylor series: P1(X, Y, Z)-lP2(X, Y, Z)=0 P3(X, Y, Z)-sP4(X, Y, Z)=0;

[0116] S53 uses the geodetic coordinates of the control points to construct the two linear equations into a system of equations and solves the equations: y = Ax0 + e

[0117] Where y∈R m×1 is the observation vector, A∈R m×n is the coefficient matrix and e∈R m×1 is the residual vector, x0∈R n×1 is the RPCs vector; m = 2k, where k is the number of geodetic coordinates of the control points, and n is 78;

[0118] Solve the system of equations: x0=N -1 A T y

[0119] Where N = A T A is a normal matrix, x0∈R n×1 is the RPCs vector;

[0120] Based on the RFM model, S54 performs RPC forward calculation on the original RPC parameters of the remote sensing image and the geodetic coordinates of the control point to obtain the second pixel coordinates of the geodetic coordinates of the control point on the remote sensing image.

[0121] S55 calculates the average of the error distances between the first pixel coordinate and the second pixel coordinate.

[0122] In an embodiment of the present invention, the geodetic coordinates of the control point are changed into projection surface coordinates according to the projection transformation formula, and the coordinate conversion is performed based on the original RPC parameters of the remote sensing image based on the RFM model, and the geodetic coordinates of the control point corresponding to the second pixel coordinates on the remote sensing image are calculated from the projection surface coordinates;

[0123] The first-level image rigorous mathematical model refers to the principle of satellite imaging. A point on the ground is formed on the image plane of the remote sensing image through the satellite sensor. The positional relationship between these two points can be expressed by the formula: (x, y) = T-1(Dlat, Dlon, h). That is, given a geodetic longitude and latitude coordinate, its pixel coordinates on the remote sensing image can be uniquely determined, realizing the calculation process of this step.

[0124] In this embodiment of the present invention, the RFM model replaces the rigorous mathematical model of the first-level image. This means that the calculation process in this step is converted to an inverse transformation formula based on the RFM model. Simply read in the corresponding RPC parameters and bring them directly into the RFM model for solution.

[0125] Adjust the height h in the pixel coordinates (x, y, h) to h base , the pixel coordinates (x, y, h) are transformed by RFM forward transformation formula base ) is transformed into geodetic coordinates (D lat1 ,D lon1 ,h base );

[0126] The forward transformation formula of the RFM model can be expressed as: (D lat ,D lon )=T(x,y,h). The calculation process in the embodiment is to read the RPC parameters and then linearize the RFM equation:

[0127] According to the projection transformation formula, the geodetic coordinates (D lat1 ,D lon1 ,h base ) changes to the projection surface coordinate (D east ,D north );

[0128] This process is a projection transformation process, which converts the latitude and longitude coordinates into projection coordinates under the projection surface (i.e. the earth ellipsoid). The projection transformation formula can be expressed as:

[0129] According to the image information of the remote sensing image, coordinate conversion is performed, and the geodetic coordinates of the control point are calculated from the projection surface coordinates (Deast Dnorth) corresponding to the second pixel coordinates (x2, y2) of the remote sensing image; x2 = (D east -a) / d x y2=(D north -b) / d y

[0130] Among them, a, b, d x , d y Both can be obtained from the image information of remote sensing images. a, b are the geodetic coordinates of the upper left point of the remote sensing image, d x , d y The resolution in the X and Y directions. x , d y , so these geometric parameters can be directly used for coordinate transformation. According to this formula, the image plane coordinates can be obtained from the projection coordinates.

[0131] In the embodiment of the present invention, the first pixel coordinates (x 1i ,y 1i ) and each second pixel coordinate (x 2i ,y 2i ) is the mean of the error distances

[0132] S6 uses the L1 regularization optimization method based on the least squares adjustment principle and introduces the set linear constraints to obtain the optimized RFM model.

[0133] In an embodiment of the present invention, we propose a new L1 regularization method based on the robust 1-norm estimation principle. This estimation principle is applicable to both RPC and residual, generating a sparse model that can be used for feature selection. To a certain extent, L1 can also prevent overfitting. The preset linear constraint conditions in step S6 are:

[0134] Where ||x||1 is the sum of the absolute values ​​of the vector x, and γ is the regularization parameter.

[0135] By introducing the set linear constraints, the RFM model is optimized.

[0136] Get the optimal solution.

[0137] S7 makes a judgment based on the mean value of the error distance and the preset threshold, selects different repair strategies, and recalculates or corrects the RPC parameters.

[0138] Specifically, the different repair strategies selected in step S7 are:

[0139] When the mean value of the error distance meets the first judgment condition, the first restoration strategy is selected, and the RPC is recalculated according to the original RPC parameters of the image to be restored and the optimized RFM model to obtain the restored RPC parameters;

[0140] When the mean value of the error distance meets the second judgment condition, the second restoration strategy is selected, and RPC correction is performed according to the geodetic coordinates, the original RPC parameters of the image to be restored, and the optimized RFM model to obtain the restored RPC parameters.

[0141] In the embodiment of the present invention, in step S7, the first judgment condition is that the mean of the error distance is ≥ the preset threshold; the preset threshold is 30; when the mean of the error distance is Select the first repair strategy, use the optimized RFM model, and recalculate the RPC parameters according to the geodetic coordinates of the control points to obtain the repaired RPC parameters.

[0142] The corresponding first repair strategy is:

[0143] Using the optimized RFM model, the RPC parameters are recalculated according to the geodetic coordinates of the control points to obtain the repaired RPC parameters;

[0144] The first repair strategy includes:

[0145] In the optimized RFM model, the preset linear constraints are updated as follows: min x,e ||x||1+||e||1

[0146] Under the updated linear constraints, according to the geodetic coordinates of the control points, the equation system is constructed as follows: y=Ax+e

[0147] Solve the system of equations: x = N -1 A T y

[0148] Where x is the repaired RPC parameter.

[0149] The second judgment condition is that the mean of the error distance is less than the preset threshold;

[0150] When the mean error distance The second repair strategy is selected, and the RPC parameters are recalculated according to the geodetic coordinates of the control points using the optimized RFM model to obtain the repaired RPC parameters.

[0151] The corresponding second repair strategy is:

[0152] Using the optimized RFM model, the optimal solution of the offset transformation is estimated according to the geodetic coordinates of the control points, and the original RPC parameters are corrected to obtain the repaired RPC parameters.

[0153] The second repair strategy includes:

[0154] Use the original RPC parameter as the initial value x initial ;

[0155] In the optimized RFM model, the preset linear constraints are updated as follows: min Δx,e ||Δx||1+||e||1;

[0156] Under the updated linear constraints, according to the geodetic coordinates of the control points, the equation system is constructed as follows: y=Ax initial +AΔx+e;

[0157] Solve the system of equations: Δx=N -1 A T y-Ax initial

[0158] where Δx is the optimal solution of the offset transformation amount;

[0159] Based on the initial value x initial , adding Δx, that is, x = x initial +Δx, where x is the repaired RPC parameter.

[0160] The above repair strategy can be applied to three possible cases between m and n, namely m > n, m < n, and m = n.

[0161] For example, in the case of a limited number of GCPs (m < n), we would expect to select more RPCs rather than residuals, and compensate by increasing the cost coefficient of residuals relative to RPCs.

[0162] When m > n, the system of equations is overfitted, and a better reformulated form of the above linear constraint condition is: min x,e ||x||1 + c||e||1

[0163] Constructed as a system of equations: y = Ax + e

[0164] where c = n / m is the cost coefficient of e relative to x, and the linear constraint condition can be rewritten as: min x,e ||[x T ce T ||1

[0165] Constructed as a system of equations: y = Ax + c -1 ce

[0166] where c -1 = m / n, using v = ce (the scaled vector of residuals), the linear constraint condition is transformed into:

[0167] Constructed as a system of equations:

[0168] Introducing z = [x T , v T 2 , B = [A, c -1 I m , the following optimization problem is obtained: min z ||z||1 y = Bz

[0169] To transform the above optimization problem into a linear form, we use the method proposed by Emil-Simkuë. The vector z, generally speaking, has positive and negative values and can be obtained from two slack vectors ɑ and β: z i = α i - β i

[0170] ​i=1,...m+n,α i , β i ≥0

[0171] From the above formula, we can infer that |z i |=α i +β i ,

[0172] Where u=[1,.....,1] T ∈R(m+n)×1, the linear version of the optimization problem is as follows:

[0173] After determining α and β through the above linear optimization problem, z = α - β can be calculated. According to the definition of z, the first parameter n = 78 is the repaired RPC parameter.

[0174] When the mean error distance The second repair strategy is selected. The optimized RFM model is used to estimate the optimal solution of the offset transformation according to the geodetic coordinates of the control points. The original RPC parameters are corrected to obtain the repaired RPC parameters.

[0175] In an embodiment of the present invention, the second repair strategy in step S7 includes:

[0176] Use the original RPC parameter as the initial value x initial ;

[0177] In the optimized RFM model, the preset linear constraints are updated as follows: min Δx,e ||Δx||1+||e||1;

[0178] Under the updated linear constraints, according to the geodetic coordinates of the control points, the equation system is constructed as follows: y=Ax initial +AΔx+e;

[0179] Solve the system of equations: Δx=N -1 A T y-Ax initial

[0180] Among them, Δx is the optimal solution of the offset transformation;

[0181] After determining α and β through the above linear optimization problem, z = α - β can be calculated. According to the definition of z, the first parameter n = 78 is Δx. At the initial value x initial On the basis of Δx, we can get x=x initial +Δx, where x is the repaired RPC parameter.

[0182] The above embodiments are only used to illustrate the present invention, and are not intended to limit the present invention. Ordinary technicians in the relevant technical field may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions also fall within the scope of the present invention. The scope of patent protection of the present invention should be defined by the claims.

Claims

1. A remote sensing image RPC repair method, characterized in that: The following steps are involved: S1 obtains the image to be restored, reference image, geographic information data and DEM data; S2 uses SIFT algorithm to extract and match features between the image to be repaired and the reference image to obtain matching point pairs and corresponding first pixel coordinates; S3 performs geographic coordinate conversion on the first pixel coordinate of the matching point pair according to the geographic information data to obtain the latitude and longitude coordinates of the matching point pair; S4 interpolates the longitude and latitude coordinates of the matching point pairs according to the DEM data to obtain geodetic coordinates, and uses the matching point pairs corresponding to the geodetic coordinates as control points; S5 constructs an RFM model, performs RPC positive calculation on the original RPC parameters of the image to be repaired and the geodetic coordinates, obtains the second pixel coordinates corresponding to the geodetic coordinates on the remote sensing image, and respectively calculates the mean values ​​of the error distances of the corresponding first pixel coordinates and the second pixel coordinates; S6 uses the L1 regularization optimization method based on the least squares adjustment principle, introduces the set linear constraints, and obtains the optimized RFM model; S7 determines based on the mean value of the error distance and the preset threshold, selects different repair strategies, and recalculates or corrects the RPC parameters; Step S5 includes: S51 builds the RFM model: Among them, P i (i=1,2,...,4) is a third-order polynomial, (l,s) is the normalized coordinate of the second pixel, (X,Y,Z) is the normalized coordinate of the control point, and the values ​​are all between [-1,1]; S52 expands (l,s) into two linear equations according to the Taylor series: P1(X,Y,Z)-lP2(X,Y,Z)=0 P3(X,Y,Z)-sP4(X,Y,Z)=0; S53 uses the geodetic coordinates of the control points to construct the two linear equations into a system of equations and solves the equations: y=Ax0+e Where y∈R m×1 is the observation vector, A∈R m×n is the coefficient matrix and e∈R m×1 is the residual vector, x0∈R n×1 is the RPCs vector; m = 2k, where k is the number of geodetic coordinates of the control points, and n is 78; Solve the system of equations: x0=N -1 A T the Where N = A T A is a normal matrix, x0∈R n×1 is the RPCs vector; S54 performs RPC positive calculation on the original RPC parameters of the remote sensing image and the geodetic coordinates of the control point based on the RFM model to obtain the second pixel coordinates of the geodetic coordinates of the control point on the remote sensing image; S55 calculates the average of the error distances between the first pixel coordinate and the second pixel coordinate.

2. A remote sensing image RPC repair method according to claim 1, characterized in that: The preset linear constraint condition in step S6 is: Among them, ||x||1 is the sum of the absolute values ​​of the vector x, and γ is the regularization parameter.

3. A remote sensing image RPC repair method according to claim 1, characterized in that: The different repair strategies selected in step S7 are: When the mean value of the error distance meets the first judgment condition, the first restoration strategy is selected, and the RPC is recalculated according to the original RPC parameters of the image to be restored and the optimized RFM model to obtain the restored RPC parameters; When the mean of the error distances meets the second judgment condition, the second restoration strategy is selected, and RPC correction is performed according to the geodetic coordinates, the original RPC parameters of the image to be restored, and the optimized RFM model to obtain the restored RPC parameters.

4. A remote sensing image RPC repair method according to claim 3, characterized in that: In step S7, the first determination condition is that the mean of the error distances ≥ a preset threshold; The corresponding first repair strategy is: Using the optimized RFM model, the RPC parameters are recalculated according to the geodetic coordinates of the control points to obtain the repaired RPC parameters; The second judgment condition is that the mean of the error distance is less than the preset threshold; The corresponding second repair strategy is: Using the optimized RFM model, the optimal solution of the offset transformation is estimated according to the geodetic coordinates of the control points, and the original RPC parameters are corrected to obtain the repaired RPC parameters.

5. A remote sensing image RPC repair method according to claim 3, characterized in that: The first repair strategy in step S7 includes: In the optimized RFM model, the preset linear constraints are updated as follows: Minutes z,e ||z||1+||e||1 Under the updated linear constraints, the system of equations is constructed according to the geodetic coordinates of the control points: y=Az+e Solve the system of equations: z=N -1 A T y Among them, z is the RPC parameter after repair.

6. A remote sensing image RPC repair method according to claim 3, characterized in that: The second repair strategy in step S7 includes: Use the original RPC parameter as the initial value z initial ; In the optimized RFM model, the preset linear constraints are updated as follows: min Δz,e ||Δz||1+||e||1; Under the updated linear constraints, the system of equations is constructed according to the geodetic coordinates of the control points: y=Az initial +AΔz+e; Solve the system of equations: Δz=N -1 A T y-Az initial Among them, Δz is the optimal solution of the offset transformation; At the initial value z initial On the basis of , add Δz, that is, z = z initial +Δz, where z is the repaired RPC parameter.

7. A remote sensing image RPC repair method according to claim 1, characterized in that: Step S2 includes: S21 determines a frame of image in the same space as a reference image based on the remote sensing image, and performs image enhancement on the two images so that the two images have interactive characteristics in grayscale; S22 uses the constraints of the geographic parameter information of the remote sensing image to perform block matching on the two images, divides the remote sensing image into a plurality of sub-area image blocks to be matched, and obtains a corresponding reference sub-area image block according to the projection coordinates of each image block; S23 uses Shi_Tomasi algorithm to move the regional search window in all directions, observes the average value of the energy expansion of the window, and then determines the central pixel point whose window energy change is greater than the specified threshold as the corner point; S24 uses SIFT descriptors with scale and rotation invariance to establish feature descriptors to describe the features of the detected feature points; After constructing the SIFT descriptor, S25 performs rough feature matching based on the KD tree nearest neighbor query algorithm, uses the KD tree for shortest distance search, and removes cluttered feature points that have no homonymous points; S26 uses a bidirectional matching strategy to initially match homonymous points, and then uses the overall relaxation method to eliminate the remaining homonymous points whose reliability cannot be determined; S27 uses random sampling consistency algorithm and least squares iteration method to eliminate mismatched pairs in the candidate matching point set; S28 uses a greedy algorithm to perform homogenization processing on the feature points to obtain a final uniformly distributed matching point pair, and obtains the first pixel coordinates corresponding to each matching point in the remote sensing image; The error distance in step S4 is calculated by any one of Euclidean distance, similarity, and angle cosine.

Images