Control point extraction method and system suitable for geometric correction of remote sensing image

Abstract

The invention discloses a control point extraction method and system suitable for remote sensing image geometric correction. The method comprises the steps that a multi-scale Gaussian pyramid is constructed for a remote sensing image to be corrected and a reference image; calculating a feature response value of each layer of Gaussian pyramid image of the remote sensing image to be corrected and the reference image, and setting a dynamic threshold to screen feature points; performing feature point matching on the Gaussian pyramid image layer by layer through geometric constraint and radiation consistency to obtain a matching point set; eliminating gross errors in the matching point set by adopting a random sampling consistency algorithm, and selecting feature points on the to-be-corrected remote sensing image in matching point pairs meeting requirements as control points for geometric correction; and performing geometric correction on the remote sensing image to be corrected by using the screened control points and a quadratic polynomial model to obtain a corrected remote sensing image. The control point extraction method provided by the invention is efficient, robust and high in adaptability, and can meet the requirement of high-precision geometric correction.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of image feature extraction, and particularly relates to a control point extraction method and system suitable for geometric correction of remote sensing images. BACKGROUND

[0002] The geometric correction of remote sensing images is a key link in remote sensing data processing, and its purpose is to eliminate the geometric deformation of images caused by factors such as sensor attitude, terrain undulation, and earth curvature, so as to ensure the spatial position accuracy of the images. The core of geometric correction lies in the selection of control points, and the accuracy and distribution of control points directly affect the reliability of the correction results. The traditional control point extraction method mainly relies on manual selection, that is, the same name feature points are identified on the images and reference data (such as topographic maps, orthophotos, etc.) by manual visual interpretation. However, this method is inefficient and greatly influenced by the subjective experience of the operator, and it is difficult to meet the needs of large-scale remote sensing data processing.

[0003] With the development of computer vision and artificial intelligence technology, automatic or semi-automatic control point extraction methods have gradually become a research hotspot. Existing methods mainly include feature point extraction techniques based on image matching, such as Scale-invariant Feature Transform (SIFT), Speeded Up Robust Features (SURF), and feature matching algorithms based on deep learning. These methods have improved the efficiency of control point extraction to some extent, but still have high matching error rates and insufficient adaptability in complex scenes (such as cloud and fog coverage, dramatic changes in ground objects, and low contrast areas). In addition, traditional feature matching algorithms are sensitive to radiation differences (such as images of different time phases and different sensors), which can easily lead to mismatching and affect the accuracy of geometric correction.

[0004] Therefore, there is an urgent need for an efficient, robust, and adaptable control point extraction method and system that can automatically process multi-source remote sensing images, overcome matching difficulties in complex scenes, and optimize the distribution and accuracy of control points to meet the needs of high-precision geometric correction. SUMMARY

[0005] The present application provides a control point extraction method suitable for geometric correction of remote sensing images, which includes the following steps: Step 1: Constructing a multi-scale Gaussian pyramid for the to-be-corrected remote sensing image and the reference image, respectively; Step 2: Calculating the feature response value of each layer of the Gaussian pyramid image of the to-be-corrected remote sensing image and the reference image, and setting a dynamic threshold to screen the feature points; Step 3, feature point matching is performed on the Gaussian pyramid image layer by layer through geometric constraint and radiation consistency, and a matching point set is obtained; Step 4, the random sample consensus algorithm is used to eliminate gross errors in the matching point set, and the feature points on the remote sensing image to be corrected in the selected matching point pairs meeting the requirements are selected as control points for geometric correction; Step 5, using the control points selected in step 4, the remote sensing image to be corrected is geometrically corrected using a quadratic polynomial model to obtain a corrected remote sensing image.

[0006] Further, the number of layers of the Gaussian pyramid in step 1 is closely related to the size of the input image, and the calculation formula is as follows: (1) In the formula, L is the number of layers of the Gaussian pyramid, W is the width of the input image, H is the height of the input image, the input image is the remote sensing image to be corrected or the reference image, floor () represents the rounding operation towards negative infinity.

[0007] For the first i layer Gaussian pyramid image, the standard deviation of the Gaussian function is The calculation method is as follows: (2) In the formula, is the initial standard deviation, f is the scale factor.

[0008] Further, the feature response value in step 2 is calculated using the self-similarity value of the pixel, and the specific calculation formula is as follows: (3) (4) In the formula, is the self-similarity value of the pixel point q ; W is a local window in the Gaussian pyramid image with the pixel point q as the center; represents the gray value of the pixel point q ; represents the average value of the pixel gray value in the local window W; represents the gray value of the pixel point q after the pixel point d moves q-d pixels; is a local window in the Gaussian pyramid image with the pixel point q-d as the center; represents the average value of the pixel gray value in the local window ; represents the standard deviation of the pixel gray value in the local window W; represents the standard deviation of pixel gray value in the local window ; represents the feature response value of the pixel point q ; k represents the number of pixel points participating in the calculation of the feature response value, the value is set according to actual needs; represents the self-similarity value of the i th pixel point and the current pixel point q .

[0009] The feature response value of each pixel point can be obtained through formula (4). After non-maximum suppression is performed on all feature response values, a plurality of local maximum values are obtained. A dynamic threshold T is set. If the local maximum value is greater than T, the pixel point corresponding to the local maximum value is a feature point.

[0010] The calculation formula of the dynamic threshold T is as follows: (5) In the formula, is the mean value of the feature response values of all pixel points before non-maximum suppression, is an adjustment factor, is the standard deviation of the feature response values of all pixel points before non-maximum suppression, and the calculation formula is as follows: (6) In the formula, N is the total number of feature response values before non-maximum suppression, is the i th feature response value.

[0011] Further, the step 3 adopts a bilateral matching method to obtain an initial matching point set of the Gaussian pyramid image of the to-be-corrected remote sensing image and the Gaussian pyramid image of the reference image , wherein and are coordinates of the i th feature point in the Gaussian pyramid image of the to-be-corrected remote sensing image and the Gaussian pyramid image of the reference image respectively, is the number of initial matching points; after the initial matching point set M is obtained, a fast iterative consistency algorithm is used to calculate a geometric transformation model H between the Gaussian pyramid image of the to-be-corrected remote sensing image and the Gaussian pyramid image of the reference image, and the geometric error of the feature points of the Gaussian pyramid image of the to-be-corrected remote sensing image and the Gaussian pyramid image of the reference image is calculated by using H. The specific calculation formula is as follows: (7) In the formula, represents the geometric error, represents a Euclidean distance calculation symbol.

[0012] For the initial matching point set M, the radiation difference is calculated using the normalized cross-correlation algorithm, and the formula is as follows: (8) In the formula, denotes the radiation difference, and denote the coordinates of the i-th feature point in the Gaussian pyramid image of the remote sensing image to be corrected and the Gaussian pyramid image of the reference image respectively, i and denote the standard deviations of the feature points in the Gaussian pyramid image of the remote sensing image to be corrected and the Gaussian pyramid image of the reference image respectively, denotes the covariance of and denotes the covariance of and and .

[0013] The feature point pair of is retained as the final matching result, and the value range of is .

[0014] Further, in step 4, the RANSAC algorithm is used to remove gross errors in the final matching point set obtained in step 3. The feature points on the remote sensing image to be matched after removing the gross errors are sorted in ascending order of geometric error, and the first feature points with an error less than pixels are selected as control points for subsequent geometric correction, and the values of and are set according to requirements.

[0015] Further, in step 5, the calculation formula of geometric correction is as follows: (9) In the formula, denotes the pixel coordinates of the remote sensing image to be corrected, denotes the pixel coordinates of the remote sensing image after geometric correction, is the geometric polynomial coefficient, which is obtained by least squares adjustment of the control points selected in step 4 and the corresponding feature points on the reference image.

[0016] The application also provides a control point extraction system suitable for remote sensing image geometric correction, which is used to realize the control point extraction method suitable for remote sensing image geometric correction as described above.

[0017] Moreover, the system comprises a processor and a memory, the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the control point extraction method suitable for remote sensing image geometric correction as described above.

[0018] ​​​Alternatively, a readable storage medium is included, and the readable storage medium stores a computer program.

[0019] Compared with the prior art, the present application has the following advantages: 1) A multi-scale Gaussian pyramid is constructed for the to-be-corrected remote sensing image and the reference image, and the robustness of the feature detector to noise is enhanced; 2) The feature point matching is performed by using geometric constraints and radiometric consistency, and the accuracy and reliability of the matched point set are effectively improved, and the robustness of the algorithm in complex scenes is enhanced; 3) The secondary polynomial model and the high-quality control points after screening are used for geometric correction, and the corrected remote sensing image has high geometric accuracy, which can meet the needs of high-precision applications such as surveying and mapping and environmental monitoring. BRIEF DESCRIPTION OF DRAWINGS

[0020] In order to more clearly illustrate the technical solutions in the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or the prior art description. Obviously, the drawings in the following description are some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0021] Figure 1 The flowchart of the control point extraction method for remote sensing image geometric correction suitable for the embodiments of the present application.

[0022] Fig. 2(a)-Fig. 2(b) are the to-be-corrected high-resolution seven remote sensing image and the reference image of the embodiments of the present application, wherein Fig. 2(a) is the to-be-corrected high-resolution seven remote sensing image, and Fig. 2(b) is the reference image.

[0023] Fig. 3(a)-Fig. 3(d) are the chessboard diagrams of the local areas of the high-resolution seven remote sensing image after geometric correction and the reference image after superposition of the embodiments of the present application, wherein Fig. 3(a) is mainly vegetation and cultivated land, Fig. 3(b) is mainly cultivated land, Fig. 3(c) is mainly buildings and roads, and Fig. 3(d) is mainly farmland and water area. DETAILED DESCRIPTION

[0024] In order to make the purpose, technical solutions and advantages of the present application more clear, the technical solutions of the present application will be further described below in combination with the drawings and embodiments. Obviously, the described embodiments are some embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the protection scope of the present application.

[0025] Embodiment 1 As Figure 1As shown, the embodiment of the present application provides a control point extraction method suitable for geometric correction of remote sensing images, which comprises the following steps: Step 1, a multi-scale Gaussian pyramid is constructed for the to-be-corrected Gao Fen 7 remote sensing image and the reference image respectively.

[0026] The size of the to-be-corrected Gao Fen 7 remote sensing image and the reference image is the same, and the number of layers of the Gaussian pyramid is closely related to the size of the input image, and the calculation formula is as follows: (1) In the formula, L is the number of layers of the Gaussian pyramid, W is the width of the input image, H is the height of the input image, the input image is the to-be-corrected remote sensing image or the reference image, floor () represents the rounding operation towards negative infinity.

[0027] Fig. 2(a) is the to-be-corrected Gao Fen 7 remote sensing image of the embodiment, and Fig. 2(b) is the reference image of the embodiment.

[0028] For the first i layer Gaussian pyramid image, the standard deviation of the Gaussian function is The calculation method is as follows: (2) In the formula, is the initial standard deviation, f is the scale factor.

[0029] Step 2, the feature response value of each layer of the Gaussian pyramid image of the to-be-corrected Gao Fen 7 remote sensing image and the reference image is calculated, and a dynamic threshold is set to screen the feature points.

[0030] In each layer of the pyramid image, the self-similarity value of the pixel point q in the local window W centered on the pixel point q is calculated according to the following formula: (3) In the formula, is the self-similarity value of the pixel point q , the value range is ; W is the local window centered on the pixel point q in the Gaussian pyramid image, represents the gray value of the pixel point q , represents the average value of the pixel gray value in the local window W, represents the gray value of the pixel point q after the pixel point d moves q-d pixels, , is the gray value of the pixel pointq-d () as the center of the local window, Represents a local window The average grayscale value of the pixels within. This represents the standard deviation of pixel grayscale values ​​within a local window W. Represents a local window The standard deviation of pixel grayscale values ​​within the range.

[0031] The feature response value is calculated using the self-similarity value of pixels, and the calculation formula is as follows: (4) In the formula, Represents pixels q The characteristic response value; k This indicates the number of pixels involved in calculating the feature response value; the value can be set according to actual needs. Indicates the first i The pixel and the current pixel q The self-similarity value.

[0032] The feature response value of each pixel can be obtained by formula (4). Non-maximum suppression is performed on all feature response values ​​to obtain multiple local maxima. A dynamic threshold T is set. If the local maximum value is greater than T, the corresponding pixel is a feature point.

[0033] The formula for calculating the dynamic threshold T is as follows: (5) In the formula, The mean value of the feature response values ​​of all pixels before non-maximum suppression; This is an adjustment factor, with a value range of [value range missing]. ; The standard deviation of the feature response values ​​of all pixels before non-maximum suppression is calculated using the following formula: (6) In the formula, N is the total number of characteristic response values ​​before nonmaximum suppression. For the first i Each characteristic response value.

[0034] Step 3: Perform feature point matching on the Gaussian pyramid image layer by layer through geometric constraints and radiation consistency to obtain the matching point set.

[0035] The initial set of matching points for the Gaussian pyramid images of the Gaofen-7 remote sensing image to be corrected and the reference image was obtained using a bilateral matching method. ,in and These are the Gaussian pyramid images of the Gaofen-7 remote sensing image to be corrected and the reference Gaussian pyramid image, respectively.i The coordinates of the feature points This represents the initial number of matched point pairs.

[0036] After obtaining the initial matching point set M, the Fast Sample Consensus (FSC) algorithm is used to calculate the geometric transformation model H between the Gaussian pyramid image of the Gaofen-7 remote sensing image to be corrected and the Gaussian pyramid image of the reference image. H is then used to calculate the geometric error of the feature points on the Gaussian pyramid images of the Gaofen-7 remote sensing image to be corrected and the reference image. The specific calculation formula is as follows: (7) In the formula, Indicates geometric error, This represents the symbol for Euclidean distance calculation.

[0037] For the initial set of matched points M, the radiation difference is calculated using the Normalized Cross-Correlation (NCC) algorithm, as shown in the following formula: (8) In the formula, Indicates differences in radiation. and These are the Gaussian pyramid images of the Gaofen-7 remote sensing image to be corrected and the reference Gaussian pyramid image, respectively. i The coordinates of the feature points and This represents the standard deviation of feature points in the Gaussian pyramid image of the Gaofen-7 remote sensing image to be corrected and the Gaussian pyramid image of the reference image. express and The covariance.

[0038] reserve The matching point pairs are used as the final matching result. The range of values ​​is .

[0039] Step 4: Use the Random Sample Consensus (RANSAC) algorithm to remove gross errors from the matching point set, and select the feature points on the Gaofen-7 remote sensing image to be corrected from the matching point pairs that meet the requirements as control points for geometric correction.

[0040] Since the geometric errors of the feature points extracted in step 3 are of different sizes and have a large variation range, if all of them are used as control points for geometric correction, the polynomial model obtained by fitting may deviate from the true position due to the control points with excessively large geometric errors, thereby reducing the geometric correction accuracy, and therefore, all of them are not suitable for being used as control points for geometric correction. The obtained control points are sorted in ascending order of geometric errors, and a certain number of control points are selected as the control points of the Gao Fen No. 7 remote sensing image.

[0041] In the present application, 4499 matching points are obtained by using normalized cross-correlation matching, 400 matching points are obtained by using RANSAC algorithm to remove gross errors in the final matching point set obtained in step 3, and the mean error is 2.01. The 400 matching points are sorted in ascending order of geometric errors, and the first 26 matching points with an error less than 1 pixel are selected as the feature points on the remote sensing image to be corrected for subsequent geometric correction.

[0042] In step 5, the control points selected in step 4 are used to perform geometric correction on the remote sensing image to be corrected by using a quadratic polynomial model.

[0043] The calculation formula of geometric correction is as follows: (9) In the formula, X represents the pixel coordinates of the remote sensing image to be corrected, X represents the pixel coordinates of the remote sensing image to be corrected, is a geometric polynomial coefficient, which is obtained by using the 30 control points selected in step 4 and the corresponding feature points on the reference image to perform least squares adjustment.

[0044] Fig. 3(a)-Fig. 3(d) show the chessboard diagrams of the local areas of the corrected Gao Fen No. 7 remote sensing image and the reference image after superposition. Fig. 3(a) mainly shows vegetation and cultivated land, Fig. 3(b) mainly shows cultivated land, Fig. 3(c) mainly shows buildings and roads, and Fig. 3(d) mainly shows farmland and water area. It can be found from Fig. 3(a)-Fig. 3(d) that the edge features of the four local areas of the corrected Gao Fen No. 7 remote sensing image and the reference image are well aligned.

[0045] Example 2 Based on the same inventive concept, the present application also provides a control point extraction system suitable for geometric correction of remote sensing images, which comprises a processor and a memory, the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute a control point extraction method suitable for geometric correction of remote sensing images as described above.

[0046] Example 3 Based on the same inventive concept, the application further provides a control point extraction system suitable for geometric correction of remote sensing images, comprising a readable storage medium, wherein the readable storage medium stores a computer program, and the computer program, when executed, implements the control point extraction method suitable for geometric correction of remote sensing images.

[0047] In specific implementation, the method provided by the technical scheme of the application can be automatically run by a computer software technology, and the system device of the method, such as a computer readable storage medium storing the corresponding computer program of the technical scheme of the application and a computer device including the running of the corresponding computer program, should also be within the protection scope of the application.

[0048] The specific embodiments described herein are merely illustrative of the application. Various modifications or changes can be made to the described embodiments without departing from the spirit of the application. The modifications or changes are intended to be within the scope of the application as defined by the appended claims.

Claims

1. A method for extracting control points suitable for geometric correction of remote sensing images, characterized in that, Includes the following steps: Step 1: Construct multi-scale Gaussian pyramids for the remote sensing image to be calibrated and the reference image respectively; Step 2: Calculate the feature response values ​​of each layer of the Gaussian pyramid image of the remote sensing image to be corrected and the reference image, and set a dynamic threshold to filter feature points. Step 3: Perform feature point matching on the Gaussian pyramid image layer by layer through geometric constraints and radiation consistency to obtain the matching point set; Step 4: Use the random sampling consistency algorithm to remove gross errors in the matching point set, and select the feature points on the remote sensing image to be corrected that meet the requirements as control points for geometric correction. Step 5: Using the control points selected in Step 4, perform geometric correction on the remote sensing image to be corrected using a quadratic polynomial model to obtain the corrected remote sensing image.

2. The control point extraction method for geometric correction of remote sensing images as described in claim 1, characterized in that: In step 1, the number of Gaussian pyramid layers is closely related to the size of the input image, and the calculation formula is as follows: （1） In the formula, L represents the number of Gaussian pyramid levels, W represents the width of the input image, and H represents the height of the input image. The input image is either the remote sensing image to be corrected or a reference image. floor () indicates a rounding operation towards negative infinity; For the i Layered Gaussian pyramid image, standard deviation of Gaussian function The calculation method is as follows: （2） In the formula, The initial standard deviation, f is the scale factor.

3. The control point extraction method for geometric correction of remote sensing images as described in claim 1, characterized in that: In step 2, the feature response value is calculated using the self-similarity value of the pixel, and the calculation formula is as follows: （3） （4） In the formula, For pixels q The self-similarity value; W is the pixel value in the Gaussian pyramid image. q A local window centered on the user; Represents pixels q grayscale value; This represents the average grayscale value of pixels within a local window W; Represents pixels q move d The pixel after 1 pixel ( qd The grayscale value of ) It is a Gaussian pyramid image in pixels ( qd A local window centered on ) Represents a local window The average grayscale value of the pixels within the range; This represents the standard deviation of pixel grayscale values ​​within a local window W; Represents a local window Standard deviation of pixel grayscale values ​​within; Represents pixels q The characteristic response value; k This indicates the number of pixels involved in calculating the feature response value; the value can be set according to actual needs. Indicates the first i The pixel and the current pixel q The self-similarity value.

4. The control point extraction method for geometric correction of remote sensing images as described in claim 3, characterized in that: In step 2, the feature response value of each pixel can be obtained through formula (4). Non-maximum suppression is performed on all feature response values ​​to obtain multiple local maxima. A dynamic threshold T is set. If the local maximum value is greater than T, the corresponding pixel is a feature point. The formula for calculating the dynamic threshold T is as follows: （5） In the formula, This is the mean of the feature response values ​​of all pixels before non-maximum suppression. As a regulating factor, The standard deviation of the feature response values ​​of all pixels before non-maximum suppression is calculated using the following formula: （6） In the formula, N is the total number of characteristic response values ​​before nonmaximum suppression. For the first i Each characteristic response value.

5. The control point extraction method for geometric correction of remote sensing images as described in claim 1, characterized in that: In step 3, a bilateral matching method is used to obtain the initial set of matching points between the Gaussian pyramid image of the remote sensing image to be corrected and the Gaussian pyramid image of the reference image. ,in and These are the Gaussian pyramid image to be corrected and the reference Gaussian pyramid image, respectively. i The coordinates of the feature points The initial number of matching point pairs is determined. After obtaining the initial matching point set M, the geometric transformation model H between the Gaussian pyramid image to be corrected and the reference Gaussian pyramid image is calculated using the fast sample consensus algorithm. H is then used to calculate the geometric error of the feature points between the Gaussian pyramid image to be corrected and the reference Gaussian pyramid image. The specific calculation formula is as follows: （7） In the formula, Indicates geometric error, This represents the symbol for Euclidean distance calculation.

6. The control point extraction method for geometric correction of remote sensing images as described in claim 5, characterized in that: In step 3, the radiation difference is calculated using the normalized cross-correlation algorithm for the initial matching point set M, as shown in the following formula: （8） In the formula, Indicates differences in radiation. and These are the Gaussian pyramid image to be corrected and the reference Gaussian pyramid image, respectively. i The coordinates of the feature points and This represents the standard deviation of feature points in the Gaussian pyramid image of the remote sensing image to be corrected and the Gaussian pyramid image of the reference image. express and covariance; reserve The feature point pairs are used as the final matching result. The range of values ​​is .

7. The control point extraction method for geometric correction of remote sensing images as described in claim 1, characterized in that: Step 4 uses the RANSAC algorithm to remove gross errors from the final set of matching points obtained in Step 3. The feature points on the remote sensing image to be corrected after removing gross errors are then sorted according to their geometric errors from smallest to largest, and the top [points] are selected. The error is less than Feature points of each pixel are used as control points for subsequent geometric correction. , The value can be set as needed.

8. The control point extraction method for geometric correction of remote sensing images as described in claim 1, characterized in that: The formula for calculating the geometric correction in step 5 is as follows: （9） In the formula, This represents the pixel coordinates of the remote sensing image to be corrected. Represents the pixel coordinates of the geometrically corrected remote sensing image. The coefficients are geometric polynomials, obtained by least-squares adjustment using the control points selected in step 4 and their corresponding feature points on the reference image.

9. A control point extraction system suitable for geometric correction of remote sensing images, characterized in that, It includes a processor and a memory, the memory being used to store program instructions, and the processor being used to call the program instructions in the memory to execute the control point extraction method for geometric correction of remote sensing images as described in any one of claims 1-8.

10. A control point extraction system suitable for geometric correction of remote sensing images, characterized in that, The method includes a readable storage medium on which a computer program is stored, which, when executed, implements a control point extraction method for geometric correction of remote sensing images as described in any one of claims 1-8.

Images