A control point extraction method and system suitable for geometric correction of remote sensing images

Abstract

The application discloses a kind of control point extraction method and system suitable for remote sensing image geometric correction, including respectively to the construction of multiple scale Gaussian pyramid of the remote sensing image to be corrected and reference image;Calculate the feature response value of each layer Gaussian pyramid image of the remote sensing image to be corrected and reference image, set dynamic threshold to screen feature points;Through geometric constraint and radiation consistency, feature point matching is carried out on the Gaussian pyramid image layer by layer, and the matching point set is obtained;Using random sampling consistency algorithm, the gross error in the matching point set is eliminated, and the feature points on the remote sensing image to be corrected in the matching point pair meeting the requirements are selected as the control points for geometric correction;Using the screened control points, the remote sensing image to be corrected is geometrically corrected using the quadratic polynomial model, and the corrected remote sensing image is obtained.The control point extraction method proposed by the application is efficient, robust and highly adaptable, and can meet the demand of high-precision geometric correction.

Description

Technical Field

[0001] This invention belongs to the field of image feature extraction technology, and in particular relates to a method and system for extracting control points suitable for geometric correction of remote sensing images. Background Technology

[0002] Geometric correction of remote sensing imagery is a crucial step in remote sensing data processing. Its purpose is to eliminate geometric distortions caused by factors such as sensor attitude, terrain undulation, and Earth curvature, thereby ensuring the spatial accuracy of the imagery. The core of geometric correction lies in the selection of control points; the accuracy and distribution of these control points directly affect the reliability of the correction results. Traditional control point extraction methods primarily rely on manual selection, i.e., identifying corresponding feature points on the image and reference data (such as topographic maps and orthophotos) through visual interpretation. However, this method is inefficient and heavily influenced by the operator's subjective experience, making it difficult to meet the needs of large-scale remote sensing data processing.

[0003] With the development of computer vision and artificial intelligence technologies, automatic or semi-automatic control point extraction methods have gradually become a research hotspot. Existing methods mainly include image matching-based feature point extraction techniques, such as Scale-invariant Feature Transform (SIFT), Speeded-Up Robust Features (SURF), and deep learning-based feature matching algorithms. These methods have improved the efficiency of control point extraction to some extent, but they still suffer from high matching error rates and insufficient adaptability in complex scenes (such as cloud cover, drastic changes in ground features, and low-contrast areas). In addition, traditional feature matching algorithms are sensitive to radiometric differences in images (such as images from different time phases or different sensors), which can easily lead to mismatches and thus 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 challenges in complex scenarios, and optimize the distribution and accuracy of control points to meet the requirements of high-precision geometric correction. Summary of the Invention

[0005] This invention addresses the shortcomings of existing technologies by providing a control point extraction method suitable for geometric correction of remote sensing images, comprising the following steps:

[0006] Step 1: Construct multi-scale Gaussian pyramids for the remote sensing image to be calibrated and the reference image respectively;

[0007] 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.

[0008] 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;

[0009] 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.

[0010] 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.

[0011] Furthermore, 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:

[0012] (1)

[0013] 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.

[0014] For the i Layered Gaussian pyramid image, standard deviation of Gaussian function The calculation method is as follows:

[0015] (2)

[0016] In the formula, The initial standard deviation, f is the scale factor.

[0017] Furthermore, 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:

[0018] (3)

[0019] (4)

[0020] 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 ( q-dThe grayscale value of ) It is a Gaussian pyramid image in pixels ( q-d 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.

[0021] 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 maxima is greater than T, the corresponding pixel is a feature point.

[0022] The formula for calculating the dynamic threshold T is as follows:

[0023] (5)

[0024] 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:

[0025] (6)

[0026] In the formula, N is the total number of characteristic response values ​​before nonmaximum suppression. For the first i Each characteristic response value.

[0027] Furthermore, 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:

[0028] (7)

[0029] In the formula, Indicates geometric error, This represents the symbol for Euclidean distance calculation.

[0030] For the initial set of matched points M, the radiation difference is calculated using the normalized cross-correlation algorithm, as shown in the following formula:

[0031] (8)

[0032] 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 The covariance.

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

[0034] Furthermore, in step 4, the RANSAC algorithm is used 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 matched after removing gross errors are 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.

[0035] Furthermore, the calculation formula for geometric correction in step 5 is as follows:

[0036] (9)

[0037] 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.

[0038] The present invention also provides a control point extraction system suitable for geometric correction of remote sensing images, for implementing the control point extraction method for geometric correction of remote sensing images as described above.

[0039] Furthermore, it includes a processor and a memory, the memory being used to store program instructions, and the processor being used to call the stored instructions in the memory to execute a control point extraction method for geometric correction of remote sensing images as described above.

[0040] Alternatively, it may include a readable storage medium storing a computer program that, when executed, implements a control point extraction method for geometric correction of remote sensing images as described above.

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

[0042] 1) Construct multi-scale Gaussian pyramids for the remote sensing image to be calibrated and the reference image to enhance the robustness of the feature detector to noise; 2) Use geometric constraints and radiometric consistency for feature point matching to effectively improve the accuracy and reliability of the matching point set and enhance the robustness of the algorithm in complex scenarios; 3) Use a quadratic polynomial model and selected high-quality control points for geometric calibration to ensure that the calibrated remote sensing image has high geometric accuracy and can meet the needs of high-precision applications such as surveying and mapping and environmental monitoring. Attached Figure Description

[0043] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0044] Figure 1 This is a flowchart of a control point extraction method applicable to geometric correction of remote sensing images, as described in an embodiment of the present invention.

[0045] Figures 2(a) and 2(b) show the Gaofen-7 remote sensing image to be corrected and the reference image in the embodiment of the present invention, wherein Figure 2(a) is the Gaofen-7 remote sensing image to be corrected and Figure 2(b) is the reference image.

[0046] Figures 3(a)-3(d) are checkerboard maps of local areas after geometric correction of Gaofen-7 remote sensing images and superposition with reference images according to embodiments of the present invention. In Figure 3(a), the main features are vegetation and cultivated land; in Figure 3(b), the main features are cultivated land; in Figure 3(c), the main features are buildings and roads; and in Figure 3(d), the main features are farmland and water areas. Detailed Implementation

[0047] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be further described below in conjunction with the accompanying drawings and embodiments. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0048] Example 1

[0049] like Figure 1 As shown, this embodiment of the invention provides a control point extraction method suitable for geometric correction of remote sensing images, comprising the following steps:

[0050] Step 1: Construct multi-scale Gaussian pyramids for the Gaofen-7 remote sensing image to be calibrated and the reference image respectively.

[0051] The Gaofen-7 remote sensing image to be corrected is the same size as the reference image. The number of Gaussian pyramid layers is closely related to the size of the input image, and the calculation formula is as follows:

[0052] (1)

[0053] 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.

[0054] Figure 2(a) is the Gaofen-7 remote sensing image to be corrected in this embodiment, and Figure 2(b) is the reference image in this embodiment.

[0055] For the i Layered Gaussian pyramid image, standard deviation of Gaussian function The calculation method is as follows:

[0056] (2)

[0057] In the formula, The initial standard deviation, f is the scale factor.

[0058] Step 2: Calculate the feature response values ​​of each layer of the Gaussian pyramid image of the Gaofen-7 remote sensing image to be corrected and the reference image, and set a dynamic threshold to filter feature points.

[0059] In each layer of the pyramid image, pixels q In the local window W centered on the pixel point q The formula for calculating the self-similarity value is as follows:

[0060] (3)

[0061] In the formula, For pixels q The self-similarity value, with a range of values ​​being W represents the pixel representation of the Gaussian pyramid image. q A local window centered on the center. 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 ( q-d The grayscale value of ) in this embodiment , It is a Gaussian pyramid image in pixels ( q-d A local window centered on ) 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.

[0062] The feature response value is calculated using the self-similarity value of pixels, and the calculation formula is as follows:

[0063] (4)

[0064] 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.

[0065] 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.

[0066] The formula for calculating the dynamic threshold T is as follows:

[0067] (5)

[0068] 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:

[0069] (6)

[0070] In the formula, N is the total number of characteristic response values ​​before nonmaximum suppression. For the first i Each characteristic response value.

[0071] 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.

[0072] 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.

[0073] 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:

[0074] (7)

[0075] In the formula, Indicates geometric error, This represents the symbol for Euclidean distance calculation.

[0076] 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:

[0077] (8)

[0078] 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.

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

[0080] 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.

[0081] Because the geometric errors of the feature points extracted in step 3 vary considerably, using all of them as control points for geometric correction might cause the fitted polynomial model to deviate from its true position due to some control points with excessively large geometric errors, resulting in decreased geometric correction accuracy. Therefore, it is not suitable to use all of them as control points for geometric correction. The acquired control points are sorted according to their geometric errors from smallest to largest, and a certain number of control points are selected as control points for the Gaofen-7 remote sensing image.

[0082] This invention uses normalized cross-correlation matching to obtain 4499 matching points. After removing outliers from the final matching point set obtained in step 3 using the RANSAC algorithm, 400 matching points remain, with a mean error of 2.01. The 400 matching points are sorted according to their geometric errors from smallest to largest, and the first 26 matching points with errors less than 1 pixel are selected to align with feature points on the Gaofen-7 remote sensing image to be corrected as control points for subsequent geometric correction.

[0083] Step 5: Using the control points selected in Step 4, perform geometric correction on the Gaofen-7 remote sensing image to be corrected using a quadratic polynomial model.

[0084] The formula for calculating geometric correction is as follows:

[0085] (9)

[0086] In the formula, This represents the pixel coordinates of the Gaofen-7 remote sensing image to be corrected. Represents the pixel coordinates of the geometrically corrected Gaofen-7 remote sensing image. The coefficients are geometric polynomials, obtained by least-squares adjustment using the 30 control points selected in step 4 and their corresponding feature points on the reference image.

[0087] Figures 3(a)-3(d) show the checkerboard pattern of local areas after the Gaofen-7 remote sensing image and the reference image are overlaid. The main features in Figure 3(a) are vegetation and farmland; in Figure 3(b) are farmland; in Figure 3(c) are buildings and roads; and in Figure 3(d) are farmland and water features. From Figures 3(a)-3(d), it can be seen that the edge features of the four local areas in the corrected Gaofen-7 remote sensing image and the reference image are well aligned.

[0088] Example 2

[0089] Based on the same inventive concept, the present invention also provides a control point extraction system suitable for geometric correction of remote sensing images, including 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 the control point extraction method for geometric correction of remote sensing images as described above.

[0090] Example 3

[0091] Based on the same inventive concept, the present invention also provides a control point extraction system suitable for geometric correction of remote sensing images, including a readable storage medium on which a computer program is stored. When the computer program is executed, it implements the control point extraction method for geometric correction of remote sensing images as described above.

[0092] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0093] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to replace them, without departing from the spirit of the invention or exceeding the scope 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. 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 ( q d The grayscale value of ) It is a Gaussian pyramid image in pixels ( q d 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; 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. 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; 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 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.

4. The control point extraction method for geometric correction of remote sensing images as described in claim 3, 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 .

5. 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.

6. 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.

7. 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-6.

8. 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-6.

Images