A collaborative orthorectification method for panchromatic and multispectral images

Abstract

The invention discloses a method for collaborative orthorectification of panchromatic multispectral images in the field of satellite remote sensing image processing technology. The specific steps are as follows: (1) performing rough orthorectification on the panchromatic multispectral image to be corrected to preliminarily obtain the resolution and range of the image; (2) sorting the search order of the image homonymous points according to the resolution; (3) obtaining the data of the intersection area between the images according to the sorting order, generating a grayscale image by using LAB color space conversion and linear stretching and other enhancement methods, and automatically selecting homonymous points by grid division using an artificial intelligence algorithm; (4) screening homonymous points based on a homography matrix; (5) using the homonymous points to solve affine transformation parameters; (6) applying the affine transformation parameters to the image to be corrected and then correcting it. This solves the problem of mismatch between panchromatic and multispectral images and mismatch between the projection information of panchromatic and multispectral images and DEM, and improves the performance of orthorectification.

Description

Technical Field

[0001] The invention belongs to the technical field of satellite remote sensing image processing, and specifically is a method for collaborative orthorectification of panchromatic multispectral images. Background Art

[0002] Panchromatic multispectral image orthorectification is a key technology in remote sensing image processing, which aims to project the acquired remote sensing image into the coordinate system of the earth's surface, eliminate the terrain effect, and achieve accurate geographic information extraction. This technology is mainly used in aerial and satellite remote sensing systems, where panchromatic images provide high-resolution grayscale information and multispectral images provide color information in different bands. In GIS applications, users often need to integrate remote sensing images with geographic data for spatial analysis and mapping. Orthorectification is a necessary step to accurately map remote sensing images to the earth's surface to ensure consistency and accuracy between different data sources.

[0003] Orthorectification usually relies on digital elevation models to obtain terrain information. DEM (digital elevation model) provides elevation data of the surface, allowing the correction process to take into account the three-dimensional shape of the surface, thereby more accurately projecting the image onto the earth's surface. Panchromatic images have a strong terrain effect due to their high resolution, while multispectral images contain information from multiple bands, which helps to perform more classification and analysis of land objects.

[0004] However, the orthorectification process still has the following problems: (1) Panchromatic and multispectral images do not match. Due to different sensor performance or errors in the data acquisition process, there may be slight differences in the geometric characteristics of panchromatic and multispectral images. This difference will affect the effect of subsequent fusion processing.

[0005] (2) The projection information of panchromatic and multispectral images does not match the DEM. The rational polynomial parameters (rpc parameters) of panchromatic and multispectral images are provided by the satellite, while the DEM used for orthorectification is specified by the user. The mismatch between the two will cause the image to be projected on the wrong terrain, resulting in large deformation and distortion.

[0006] In the face of the above problems in the orthorectification process, the traditional solution is to use control points, manually puncture the same-name points on the panchromatic multispectral image, combine with a certain machine point finding algorithm, and recalculate the RPC parameters through adjustment.

[0007] However, the reference image of the control point may be quite different from the panchromatic and multispectral images to be corrected, which makes it difficult to select the same-name points. However, since the panchromatic and multispectral images are imaged at the same time, it is easier for them to find the same-name points. In addition, resolving the rational polynomial parameters requires a large number of control points. Generally, it is necessary to first establish the affine transformation of the image to be corrected and the control points to generate virtual control points for supplementation. Summary of the invention

[0008] In order to solve the problems existing in the above-mentioned orthorectification process, the purpose of the present invention is to provide a method for collaborative orthorectification of panchromatic and multispectral images based on reference images, sort the correction order according to the resolution of each image, use the panchromatic and multispectral image serial sorting method to search for image homonymous points, and improve the performance of orthorectification.

[0009] In order to achieve the above object, the technical solution of the present invention is as follows:

[0010] A method for collaborative orthorectification of panchromatic multispectral images comprises the following steps: Step 1: selecting a reference image, the reference image being a digital orthophoto image, and the reference image and the panchromatic multispectral image to be rectified have a geographically overlapping area;

[0011] Step 2: Perform coarse orthorectification on the panchromatic and multispectral images to be corrected. Both the panchromatic and multispectral images to be corrected include rational polynomial model parameters. The coarse orthorectification step is to read the projection information of the reference image, perform coarse orthorectification on the panchromatic and multispectral images to be corrected according to the rational polynomial model parameters and the existing digital elevation model, and project them, so that the coarsely corrected panchromatic and multispectral images are consistent with the projection information of the reference image.

[0012] Step 3: Sort the panchromatic and multispectral image correction order, read the reference image and the resolution of the panchromatic and multispectral images after rough correction, calculate the resolution ratios between the reference image and the panchromatic and multispectral images, respectively, take the panchromatic or multispectral image with a smaller ratio as sorted image 1, and take the panchromatic or multispectral image with a larger ratio as sorted image 2;

[0013] Step 4: Search for the same-name points between the reference image and the sorted image 1 described in step 3, resample the reference image to make the resolution of the reference image consistent with the resolution of the sorted image 1 described in step 3, obtain the intersection area of ​​the reference image and the sorted image 1 and crop it, and establish a grid on the intersection area, loop each grid, convert the reference image and the sorted image 1 in the grid into grayscale images, and then use the lightglue algorithm to automatically search for the same-name points;

[0014] Step 5: Filter the points with the same name obtained in step 4 and use the RANSAC algorithm to filter the points with the same name;

[0015] Step 6: According to the same-name points selected in step 5, the affine transformation parameters of the sorted image 1 and the reference image are adjusted and solved. The adjustment method is the Levenberg-Marquardt least squares optimization algorithm.

[0016] Step 7: Perform affine transformation on the sorted image 1 based on the affine transformation parameters obtained in step 6, and perform step 2 again on the result after affine transformation to obtain the orthorectified result of the sorted image 1;

[0017] Step 8: Search for the same-name points between the orthorectified result of the sorted image 1 obtained in step 7 and the sorted image 2 obtained in step 3. First, resample the orthorectified result of the sorted image 1 to make the resolution of the sorted image 1 consistent with that of the sorted image 2 in step 3.

[0018] Substitute sorted image 1 and sorted image 2 into step 4, and search for the orthorectified result of sorted image 1 and the same-name points of sorted image 2;

[0019] Step 9: Filter the orthorectified results of sorted image 1 obtained in step 8 and the same-name points of sorted image 2, and use the RANSAC algorithm to filter the same-name points;

[0020] Step 10: According to the selected same-name points obtained in step 9, the affine transformation parameters of sorted image 2 and sorted image 1 are adjusted and solved, and the adjustment method is Levenberg-Marquardt least squares optimization algorithm;

[0021] Step 11: Perform an affine transformation on the sorted image 2 based on the affine transformation parameters obtained in step 10, and perform step 2 again on the result after the affine transformation to obtain the orthorectified result of the sorted image 2.

[0022] Further, in step 1, the image size of the geographical overlap area between the reference image and the panchromatic multispectral image to be corrected is less than or equal to 1024*1024 pixels, and the cloud coverage in the geographical overlap area is less than or equal to 20%.

[0023] Further, in step 3, the ratio of the resolution between the reference image and the panchromatic image and the multispectral image is calculated as follows:

[0024]

[0025]

[0026] where ratio ms , ratio pan Respectively represent the ratio of the resolution of the reference image to the panchromatic multispectral image, R ref ,R pan ,Rms Represent the resolution of the reference image, panchromatic image, and multispectral image, respectively.

[0027] Further, in step 4, the size of the grid is a square, and the side length of the square is a positive integer multiple of 1024.

[0028] Furthermore, in step 4, the method of converting the reference image and the sorted image 1 or the sorted image 2 in the grid into a grayscale image is to first determine the type of the sorted image:

[0029] If the reference image or sorted image is a panchromatic image, the image data in the grid is normalized to 0-1 using a 1% linear stretch. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image.

[0030] If the reference image or sorted image is a multispectral image, select the data of the first three bands, normalize the image data in the grid to the maximum value, convert the normalized data to the LAB color space, select the L band in the LAB color space and perform a 1% linear stretch to normalize it to 0-1. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image.

[0031] Furthermore, in step 4, after the grid is established, the grid is thinned out and a subset of the grid is selected to search for points with the same name.

[0032] Further, in step 5 or step 9, the RANSAC algorithm screens the same-name points based on the homography matrix, the reprojection error of the algorithm is set to 1, the maximum number of iterations of the algorithm is set to 50,000, and the confidence of the algorithm is set to 0.999999.

[0033] Further, in step 6 or step 10, the loss function of the adjustment solution method is set to a linear loss function, the adjustment formula is an affine transformation formula, and the iterative initial values ​​of the six parameters in the affine transformation formula are all set to 0.

[0034] The following beneficial effects are achieved by adopting the above scheme: the present invention performs a rough orthorectification on the panchromatic multispectral image to be corrected, preliminarily obtains the resolution and range of the image, provides a basic reference for the accuracy of the subsequent orthorectification, and provides a basic reference for solving the mismatch of the panchromatic or multispectral image; sorts the search order of the image homonymous points according to the resolution, that is, it is convenient to pre-process the classification according to the sorting order, and strengthen the connection between the homonymous points to ensure the accuracy of the image projection;

[0035] Instead, they are processed sequentially according to the sorting order to obtain the data of the intersection areas between the images, and the grayscale image is generated using enhancement methods such as LAB color space conversion and linear stretching. By dividing the grid, the artificial intelligence algorithm is used to automatically select the points with the same name, and the difference between the projection information of the panchromatic or multispectral image and the DEM is reduced. The coordinated orthorectification of the panchromatic and multispectral images is realized, thereby improving the performance of the orthorectification. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 Flow chart of a method for collaborative orthorectification of panchromatic and multispectral images according to an embodiment of the present invention.

[0037] Figure 2 The reference image of the panchromatic multispectral image collaborative orthorectification method according to the embodiment of the present invention and the panchromatic multispectral image to be corrected have a geographically overlapping area map.

[0038] Figure 3 for Figure 2 Schematic diagram of grid division results.

[0039] Figure 4 for Figure 2 Grayscale image corresponding to the 10th grid in .

[0040] Figure 5 for Figure 2 The search results of the same-name points corresponding to the 20th grid in the figure.

[0041] Figure 6 for Figure 2 The orthorectification result map of sorted image 1.

[0042] Figure 7 This is the grayscale image corresponding to the 560th grid.

[0043] Figure 8 This is the search result map of the points with the same name corresponding to the 560th grid.

[0044] Fig. 9 for Figure 2 The orthorectification result of sorted image 2. DETAILED DESCRIPTION

[0045] The following is further described in detail through specific implementation methods:

[0046] The embodiment is basically as shown in the attached Figure 1 Figure 1 shows a method for collaborative orthorectification of panchromatic and multispectral images. The reference image used in this example is a Sentinel-2 image with a resolution of 10 meters. The image to be corrected comes from the BJ3A1 sensor, where the panchromatic image resolution is 0.5 meters and the multispectral image resolution is 2 meters.

[0047] Step 1: Select a reference image, which is a digital orthophoto image, and the reference image has a geographical overlap with the panchromatic multispectral image to be corrected;

[0048] For example, in this embodiment, a 10-meter resolution Sentinel-2 image is selected as the reference image. The reference image is a high-precision digital orthophoto with three bands and 8-bit depth encoding. The geographical overlap area between the reference image and the BJ3A1 image to be corrected is as follows: Figure 2 shown.

[0049] Step 2: Perform coarse orthorectification on the panchromatic and multispectral images to be corrected. Both the panchromatic and multispectral images to be corrected contain rational polynomial model parameters (rpc parameters); the coarse orthorectification step is to read the projection information of the reference image, perform coarse orthorectification and projection on the panchromatic and multispectral images to be corrected according to the rational polynomial model parameters and the existing digital elevation model (DEM), so that the coarsely corrected panchromatic and multispectral images are consistent with the projection information of the reference image, where the digital elevation model is NASA DEM with a plane resolution of 30 meters, and the coarsely corrected panchromatic and multispectral images are consistent with the projection information of the reference image.

[0050] Step 3: Sort the panchromatic and multispectral image correction order. Read the resolution of the reference image and the roughly corrected panchromatic and multispectral image, where the resolution of the reference image is 10 meters, the resolution of the panchromatic image is 0.5 meters, and the resolution of the multispectral image is 2 meters. Then calculate the resolution ratio between the reference image and the panchromatic image and the multispectral image respectively, and take the panchromatic image or multispectral image corresponding to the smaller ratio as the sorted image 1, and take the panchromatic image or multispectral image corresponding to the larger ratio as the sorted image 2.

[0051] The calculation formula for the ratio of the resolution between the reference image and the panchromatic image and the multispectral image is:

[0052]

[0053]

[0054] where ratio ms , ratio pan Respectively represent the ratio of the resolution of the reference image to the panchromatic multispectral image, R ref ,R pan ,R ms Represent the resolution of the reference image, panchromatic image, and multispectral image, respectively.

[0055] Step 4: Search for the same-name points between the reference image and the sorted image 1 described in step 3. Resample the reference image to make the resolution of the reference image consistent with the resolution of the sorted image 1 described in step 3, obtain the intersection area of ​​the reference image and the sorted image 1 and crop it, and build a grid on the intersection area. Loop through each grid, convert the reference image and sorted image 1 in the grid into grayscale images, and use the lightglue algorithm to automatically search for the same-name points. The grid size is 1024*1024, and there are 49 grids in total, such as Figure 3 The schematic diagram of the grid division results is shown.

[0056] Among them, the reference image and sorted image 1 are both multispectral images, so only the data of the first three bands are selected, and the image data in the grid is normalized to the maximum value. The normalized data is converted to the LAB color space, and the L band in LAB is selected for 1% linear stretching to normalize it to 0-1. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image. Figure 4 As shown, the grayscale image corresponding to the 10th grid is displayed.

[0057] The lightglue algorithm described in step 4 is an algorithm proposed in the paper "LightGlue: Local Feature Matching at Light Speed", and its parameters are the input image and the resampling factor. In the present invention, the resampling factor is set to 1.

[0058] Step 5: Filter the same-name points filtered in step 4. Use the RANSAC algorithm to filter the same-name points. The RANSAC algorithm filters the same-name points based on the homography matrix, the reprojection error of the algorithm is set to 1, the maximum number of iterations of the algorithm is set to 50,000, and the confidence of the algorithm is set to 0.999999. Figure 5 The result of searching for points with the same name corresponding to the 20th grid is shown.

[0059] Step 6: According to the selected same-name points obtained in step 5, the affine transformation parameters of sorted image 1 and the reference image are adjusted. The adjustment method is the Levenberg-Marquardt least squares optimization algorithm, which uses a damped Gauss-Newton method. The loss function of this method is set to a linear loss function, the adjustment formula is the affine transformation formula, and the iterative initial values ​​of the six parameters in the affine transformation formula are all set to 0.

[0060] Step 7: Perform affine transformation on sorted image 1 based on the affine transformation parameters obtained in step 6, and perform step 2 again on the result after affine transformation to obtain the orthorectified result of sorted image 1. Figure 6The orthorectification results of sorted image 1 are shown. The resolution and range of the image are initially obtained to provide a basic reference for the accuracy of subsequent orthorectification and to solve the mismatch of panchromatic or multispectral images; the search order of the image homonymous points is sorted according to the resolution, which is convenient for early preprocessing of classification according to the sorting order, and the connection between homonymous points is strengthened to ensure the accuracy of image projection.

[0061] Step 8: Search for the same-name points between the orthorectified result of sorted image 1 obtained in step 7 and the sorted image 2 obtained in step 3. First, resample the orthorectified result of sorted image 1 to make the resolution of sorted image 1 consistent with that of sorted image 2 described in step 3. Find the intersection area of ​​the orthorectified result of sorted image 1 and sorted image 2 and crop it, and build a grid on the intersection area. Loop through each grid, convert the orthorectified result of sorted image 1 and sorted image 2 in the grid into grayscale images, and use the lightglue algorithm to automatically search for the same-name points. In this example, the grid size is 1024*1024, and a total of 720 grids are divided.

[0062] Among them, sorted image 1 is a multispectral image, and sorted image 2 is a panchromatic image.

[0063] For sorted image 1, only the data of the first three bands are selected, and the image data in the grid is normalized to the maximum value. The normalized data is converted to the LAB color space, and the L band in LAB is selected for 1% linear stretching to normalize it to 0-1. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image.

[0064] For sorted image 2, a 1% linear stretch is used to normalize the image data in the grid to 0-1. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image. Figure 7 As shown, the grayscale image corresponding to the 560th grid is displayed.

[0065] Step 9: Filter the orthorectified results of sorted image 1 and the same-name points of sorted image 2 obtained in step 8. Use the RANSAC algorithm to filter the same-name points. The RANSAC algorithm filters the same-name points based on the homography matrix, the reprojection error of the algorithm is set to 1, the maximum number of iterations of the algorithm is set to 50,000, and the confidence of the algorithm is set to 0.999999. Figure 8 As shown, the search results of the same-name points corresponding to the 560th grid are displayed.

[0066] Step 10: According to the selected same-name points obtained in step 9, the affine transformation parameters of sorted image 2 and sorted image 1 are adjusted. The adjustment method is the Levenberg-Marquardt least squares optimization algorithm, which uses a damped Gauss-Newton method. The loss function of the method is set as a linear loss function, the adjustment formula is an affine transformation formula, and the iterative initial values ​​of the six parameters in the affine transformation formula are all set to 0.

[0067] Step 11: Perform affine transformation on sorted image 2 based on the affine transformation parameters obtained in step 10, and perform the processing of step 2 on the result after affine transformation to obtain the orthorectified result of sorted image 2. Fig. 9 As shown, the orthorectification result of sorted image 2 is displayed.

[0068] The images are processed sequentially according to the sorting order to obtain the data of the intersection area between the images, and the grayscale image is generated using enhancement methods such as LAB color space conversion and linear stretching. By dividing the grid, the artificial intelligence algorithm is used to automatically select the same-name points, reduce the difference between the projection information of the panchromatic or multispectral image and the DEM, and realize the collaborative orthorectification of the panchromatic and multispectral images, thereby improving the performance of the orthorectification and further improving the accuracy of GIS positioning.

[0069] The above is only an embodiment of the present invention, and the common knowledge such as the known specific structure and / or characteristics in the scheme is not described in detail here. It should be pointed out that for those skilled in the art, several deformations and improvements can be made without departing from the structure of the present invention, which should also be regarded as the protection scope of the present invention, and these will not affect the effect of the implementation of the present invention and the practicality of the patent. The scope of protection required by this application shall be based on the content of its claims, and the specific implementation methods and other records in the specification can be used to interpret the content of the claims.

Claims

1. A method for collaborative orthorectification of panchromatic and multispectral images, characterized in that: The method comprises the following steps: Step 1: selecting a reference image, which is a digital orthophoto image, and the reference image has a geographically overlapping area with the panchromatic multispectral image to be corrected; Step 2: Perform coarse orthorectification on the panchromatic and multispectral images to be corrected. Both the panchromatic and multispectral images to be corrected include rational polynomial model parameters. The coarse orthorectification step is to read the projection information of the reference image, perform coarse orthorectification on the panchromatic and multispectral images to be corrected according to the rational polynomial model parameters and the existing digital elevation model, and project them, so that the coarsely corrected panchromatic and multispectral images are consistent with the projection information of the reference image. Step 3: Sort the panchromatic and multispectral image correction order, read the resolution of the reference image and the roughly corrected panchromatic and multispectral images, calculate the resolution ratios between the reference image and the panchromatic and multispectral images, respectively, take the panchromatic or multispectral image with the smaller ratio as sorted image 1, and take the panchromatic or multispectral image with the larger ratio as sorted image 2; Step 4: Search for the same-name points between the reference image and the sorted image 1 described in step 3, resample the reference image to make the resolution of the reference image consistent with the resolution of the sorted image 1 described in step 3, obtain the intersection area of ​​the reference image and the sorted image 1 and crop it, and establish a grid on the intersection area, loop each grid, convert the reference image and sorted image 1 in the grid into grayscale images, and then use the lightglue algorithm to automatically search for the same-name points; Step 5: Filter the points with the same name obtained in step 4 and use the RANSAC algorithm to filter the points with the same name; Step 6: According to the same-name points selected in step 5, the affine transformation parameters of the sorted image 1 and the reference image are adjusted and solved. The adjustment method is the Levenberg-Marquardt least squares optimization algorithm. Step 7: Perform affine transformation on the sorted image 1 based on the affine transformation parameters obtained in step 6, and perform step 2 again on the result after affine transformation to obtain the orthorectified result of the sorted image 1; Step 8: Search for the same-name points between the orthorectified result of the sorted image 1 obtained in step 7 and the sorted image 2 obtained in step 3. First, resample the orthorectified result of the sorted image 1 to make the resolution of the sorted image 1 consistent with that of the sorted image 2 in step 3. Find the intersection area of ​​the orthorectified result of sorted image 1 and sorted image 2 and crop it, and build a grid on the intersection area. Loop through each grid, convert the orthorectified result of sorted image 1 and sorted image 2 in the grid into grayscale images, and use the lightglue algorithm to automatically search for points with the same name; Step 9: Filter the orthorectified results of sorted image 1 obtained in step 8 and the same-name points of sorted image 2, and use the RANSAC algorithm to filter the same-name points; Step 10: According to the selected same-name points obtained in step 9, the affine transformation parameters of sorted image 2 and sorted image 1 are adjusted and solved, and the adjustment method is Levenberg-Marquardt least squares optimization algorithm; Step 11: Perform an affine transformation on the sorted image 2 based on the affine transformation parameters obtained in step 10, and perform step 2 again on the result after the affine transformation to obtain the orthorectified result of the sorted image 2.

2. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 1, the image size of the geographical overlap area between the reference image and the panchromatic multispectral image to be corrected is less than or equal to 1024*1024 pixels, and the cloud coverage in the geographical overlap area is less than or equal to 20%.

3. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 3, the ratio of the resolution between the reference image and the panchromatic image and the multispectral image is calculated as: ， ， where ratio ms , ratio pan Respectively represent the ratio of the reference image to the multispectral image and the panchromatic image resolution, R ref ,R pan ,R ms Represent the resolution of the reference image, panchromatic image, and multispectral image, respectively.

4. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 4, the grid is a square, and the side length of the square is a positive integer multiple of 1024.

5. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 4, the method of converting the reference image and the sorted image 1 in the grid into a grayscale image is to first determine the type of the sorted image: If the reference image or sorted image is a panchromatic image, the image data in the grid is normalized to 0-1 using a 1% linear stretch. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image. If the reference image or sorted image is a multispectral image, select the data of the first three bands, normalize the image data in the grid to the maximum value, convert the normalized data to the LAB color space, select the L band in the LAB color space and perform a 1% linear stretch to normalize it to 0-1. After stretching, values ​​less than 0 are set to 0, and values ​​greater than 1 are set to 1. The normalized result is used as a grayscale image.

6. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 4, after the grid is established, the grid is thinned out and some subsets of the grid are selected to search for points with the same name.

7. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 5 or step 9, the RANSAC algorithm selects the same-name points based on the homography matrix, the reprojection error of the algorithm is set to 1, the maximum number of iterations of the algorithm is set to 50,000, and the confidence of the algorithm is set to 0.999999.

8. The method for collaborative orthorectification of panchromatic and multispectral images according to claim 1, characterized in that: In step 6 or step 10, the loss function of the adjustment solution method is set to a linear loss function, the adjustment formula is an affine transformation formula, and the initial iteration values ​​of the six parameters in the affine transformation formula are all set to 0.

Images