Local geometric precise correction method suitable for high-resolution remote sensing image

Abstract

The application relates to a local geometric fine correction method suitable for high-resolution remote sensing images, relates to the field of remote sensing image geometric fine correction, and solves the problems that the existing method is calculated for the whole image, time and calculation resources are wasted, error matching is prone to occurring by using the existing matching method, geometric fine correction fails and the like, the rational polynomial RPC coefficient of the whole image is used for calculating the longitude and latitude range of the whole image; the longitude and latitude of a local area interested by a user are used for generating an image of the local area with rough geographic coordinates by using an indirect resampling method; a perceptual hashing algorithm is used for estimating the rough coordinate offset between the local area image and a reference base map, the part of the image corresponding to the local area image is cut from the reference base map, and the local area image to be processed is registered with the local reference base map by using Fourier transform blocking. The method has significant improvement in speed, accuracy and stability, and is suitable for various meter-level and sub-meter-level high-resolution remote sensing images.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of geometric rectification of remote sensing images, and in particular to a design of a local geometric rectification method suitable for high-resolution remote sensing images. BACKGROUND

[0002] Traditional geometric rectification methods of remote sensing images usually use ground control points, satellite orbit parameters, satellite attitude parameters and sensor optical parameters to establish a geometric rectification model, so as to correct the geometric deformation of remote sensing images caused by projection methods, sensor attitude errors, terrain undulations, atmospheric refraction and earth rotation, and eliminate the coordinate errors of image elements in the geographic coordinate system, so as to correct the image elements to the correct position. In the above method, ground control points are usually selected manually, which is low in efficiency and heavy in workload, and it is difficult to meet the requirements of automatic, fast and accurate processing under the current conditions of massive remote sensing images.

[0003] The current mainstream geometric rectification method based on image matching can meet the requirements of automation. This method uses a remote sensing image that has been processed and has accurate geographic coordinates as a reference base map, and matches the same name points between the image to be processed and the reference base map. Commonly used algorithms include Scale-Invariant Feature Transform (SIFT), Speeded Up Robust Feature (SURF), etc. The matching results are constrained by homography transformation or affine transformation, and the matching results with large errors are removed. The obtained same name point pairs are used as ground control points. The transformation relationship between the image to be processed and the geometric rectified image coordinate system can be described by a polynomial model. The coefficients of the polynomial can be solved by using the ground control points. The original image is resampled by using the indirect method, so as to realize the geometric rectification of the whole image, and then the image results of the local range interested by the user are obtained by cropping.

[0004] Since the RPC coefficients describe the transformation relationship between the pixel coordinates and the geographic coordinates of the whole image, the existing geometric rectification method based on image matching needs to calculate the whole image, and then crop the calculation results to obtain the local geometric rectified image interested by the user. In actual work, the user is often only interested in the local area of the image, and the calculation of the image outside the area interested by the user will waste time and computing resources. Therefore, the purpose of the present application is to only perform geometric rectification calculation on the local area of the image interested by the user, and to avoid the calculation of the whole image in the conventional method, so as to reduce the calculation amount and improve the efficiency.

[0005] Because the acquisition time between the image to be processed and the reference base map is often long, the ground features will change significantly; the difference between the sensor and solar radiation will cause differences in brightness and color between the images; in addition, the shooting area may have repetitive textures or sparse feature points. Commonly used corresponding point matching algorithms such as SIFT and SURF will have more mismatches due to the above factors, resulting in failure of geometric fine correction. Summary of the Invention

[0006] To address the problems of existing methods that waste time and computational resources by calculating the entire image and the tendency for incorrect matching to occur when using existing matching methods, leading to failure of geometric fine correction, this invention provides a local geometric fine correction method suitable for high-resolution remote sensing images.

[0007] A method for local geometric fine correction of high-resolution remote sensing images, which is implemented by the following steps:

[0008] Step 1: Select the entire image to be processed (I) In According to the entire image I In The RPC coefficients are used to calculate the I of the entire image. In The latitude and longitude range; then set the latitude and longitude range for a local area;

[0009] Step 2: Based on the latitude and longitude range of the local area in Step 1, generate local area image I using an indirect resampling method. Loc Traversing local area images I Loc All pixels, obtain I Loc Pixel values ​​of all pixels;

[0010] Step 3: Based on the latitude and longitude range of the local area obtained in Step 1, and the RPC coefficients, calculate the reference base map I for the local area. Ref In the reference base map I Ref Move the sliding window upwards to obtain the reference base image within the sliding window.

[0011] Step 4: Use the perceptual hashing algorithm to calculate the local area image I. Loc Reference base map in the sliding window pHash fingerprint, that is: the local area image I Loc Reference base map in the sliding window After downsampling to 100×100 pixels, DCT transformation is performed, and the pHash fingerprint of the upper left 10×10 area of ​​the DCT transformed image is extracted.

[0012] Step 5: Statistical analysis of local area images I Loc Reference basemap in all sliding windows The Hamming distance between the pHash fingerprints, and based on the calculated minimum reference basemap. Obtain the reference base image for the cut.

[0013] Step 6: Apply Fourier-Melin transform to I in a block-based manner. Loc and Perform a match, and then for I Loc Resampling enables geometrical correction of the image.

[0014] The beneficial effects of this invention are as follows: Starting from the area of ​​interest to the user, this invention generates an image with rough geographic coordinates for a local area by using indirect resampling based on the RPC coefficients and latitude and longitude range of the original whole image. Then, the local reference base map corresponding to the local area is obtained by using a perceptual hash algorithm, so that the subsequent Fourier Merlin matching calculation is only performed on the local area, which greatly reduces redundant calculations and effectively improves the speed of local geometric fine correction of high-resolution remote sensing images.

[0015] Through extensive experiments, this invention proposes a series of parameters and thresholds for geometrical fine correction of high-resolution remote sensing images in local areas using the Fourier Merlin matching algorithm. These parameters and thresholds can achieve more stable and accurate results than the classic SIFT and SURF matching algorithms, even when there are variations in ground features, radiometric and color differences between the local image to be processed and the reference base map.

[0016] This invention uses the RPC coefficients of the original whole image and the latitude and longitude range of the user's region of interest to generate images with rough geographic coordinates for local areas by indirect resampling. The calculation is performed only on the images of local areas, which improves efficiency and reduces the consumption of computing resources.

[0017] In this invention, the method uses a sliding window approach and the pHash algorithm to estimate the coarse coordinate offset between the local image to be processed and the reference base map, which greatly reduces redundant calculations in the matching process of corresponding points.

[0018] Based on extensive experiments, this invention has derived a threshold suitable for local geometric fine correction of high-resolution remote sensing images, such as the sliding window size (W). Loc H Loc ), the distance of each slide (0.1W) Loc 0.1H Loc The local image is downsampled to 100×100 pixels before DCT transformation, and the size of the upper left corner area of ​​the DCT transformation result is cropped (10×10 pixels).

[0019] The method of this invention uses Fourier-Melin transform in a block-based manner to perform corresponding point matching between the local image to be processed and the reference base map, which can achieve stable and accurate matching results even when there are differences in ground features, lighting and color between the images.

[0020] Based on the coarse coordinate offset between the image to be processed and the reference base map estimated by the perceptual hash algorithm, the image to be processed I is found. In Each block in the reference image The Fourier-Melin transform is calculated for the corresponding region. Based on extensive experiments, this invention has obtained thresholds suitable for high-resolution remote sensing images. For example, in the Fourier-Melin transform, the threshold for the scale coefficient of the first phase correlation operation is set to no more than 1.7, and the response threshold for the second phase correlation operation is set to no less than 0.05. For the matched points, the top 16 matched point pairs are selected as the final matching result, sorted from largest to smallest response value. If the number of matched points is less than 8, the matching is considered to have failed, and the geometric fine correction ends. Attached Figure Description

[0021] Figure 1 This is a flowchart of the local geometric fine correction method for high-resolution remote sensing images described in this invention;

[0022] Figure 2 To obtain the reference base image by using a sliding window method Schematic diagram;

[0023] Figure 3 For the original whole remote sensing image I In Image I of local area Loc Comparison chart;

[0024] Figure 4 For I Loc Image rendering after DCT transformation;

[0025] Figure 5 This is a schematic diagram of a pHash fingerprint.

[0026] Figure 6 Schematic diagram of Fourier Merlin transform and matching of corresponding points for block segmentation;

[0027] Figure 7 (7a) and (7b) are comparison images of the original full image and the reference base map;

[0028] Figure 8 Figures (8a) and (8b) show the comparison between the geometrically corrected local area image and the reference base map. Detailed Implementation

[0029] Combination Figures 1 to 8This embodiment describes a method for local geometric fine correction of high-resolution remote sensing images. The specific steps of this method are as follows:

[0030] (1) Let the entire image to be processed be I. In According to I In The RPC coefficients (Rational Polynomial Coefficients) are used to calculate I. In The maximum and minimum range of latitude and longitude coordinates, with the maximum value of the x-coordinate being... Minimum value is The maximum value of the ordinate is Minimum value is

[0031] (2) Let the maximum value of the latitude and longitude x-coordinates of the local area of ​​interest to the user be . Minimum value is Maximum value of y-axis Minimum value is Let the resolution of the entire image to be processed be r, then the image of the local region of interest to the user is I. Loc Width W Loc High H Loc They are respectively:

[0032]

[0033] The ceil function represents rounding up, W Loc H Loc The unit is pixels;

[0034] (3) Indirect resampling generates images of local areas with coarse geographic coordinates. Traversal I Loc Given all pixels, let the pixel coordinates of the current pixel be (i,j) (counting from 1). Combining the resolution r and the latitude and longitude range of the local area in step (2), the latitude and longitude X of the center of the pixel can be calculated. Loc (i,j), Y Loc (i,j), as shown in equation (2):

[0035]

[0036] Then, based on this latitude and longitude, its position in I is calculated using the RPC coefficients. In The corresponding pixel coordinates can be calculated using conventional bilinear interpolation to obtain I. Loc The pixel value of the current pixel. For I Loc Perform the above operation on all pixels to obtain I. Loc Pixel values ​​of all pixels;

[0037] (4) Calculate the required reference base map I for the local area based on the latitude and longitude range of the local area and the RPC positioning error. Ref The range, let the maximum value of the x-coordinate be... Minimum value is The maximum value of the ordinate is Minimum value is Read the reference basemap pixel values ​​within this range. For I Ref Take a sliding window, and slide it 0.1W distance in the east-west direction each time. Loc The north-south sliding distance is 0.1H. Loc The size of the sliding window is related to I Loc Same (e.g.) Figure 1 As shown), therefore, the maximum value of the horizontal coordinate of the sliding window is Minimum value is (where i represents the number of times the object slides from west to east), and the maximum value of the ordinate is... Minimum value is (Where j represents the number of times the slider moves from south to north). Obtain the reference base map within the sliding window, denoted as...

[0038] (5) Calculate I using the Perceptual Hash (pHash) algorithm. Loc and The specific implementation methods for each pHash fingerprint are as follows. I Loc and Downsample to 100×100 pixels. If the image contains multiple bands, convert it to grayscale image according to formula (3) (where Gray represents grayscale image, and R, G, and B represent red, green, and blue bands respectively). The purpose is to reduce the amount of computation.

[0039] Gray=0.299R+0.587G+0.114B (3)

[0040] downsampling and grayscale I Loc and A Discrete Cosine Transform (DCT) is performed to transform the image to the frequency domain. Original full-image remote sensing image I In Image I of the extracted local area Loc like Figure 2 As shown, I Loc The result after DCT transformation is as follows Figure 3 As shown;

[0041] Since most of the energy of the image is concentrated in the low frequency part, the upper left corner 10x10 pixel area of the image after DCT transformation is intercepted, which contains the low frequency information of the image. The mean value of the pixel value of the upper left corner 10x10 pixel area is calculated, and the area is binarized with the mean value as the threshold. If it is greater than the mean value, it is recorded as 1, and if it is less than the mean value, it is recorded as 0. The 100-bit pHash fingerprint can be obtained, as shown in Figure 4 I Loc The hash fingerprints of describe the low frequency information of the respective images, and binarization is performed with the mean value as the threshold, which can eliminate the influence of light and color differences;

[0042] (6) The Hamming distance between the pHash fingerprints of Loc and all is calculated. The purpose of selecting the Hamming distance is to improve the speed through bit operation. The is expanded outward by 0.1W Loc and 0.1H Loc respectively, and the reference base map in this range is intercepted and recorded as The purpose of intercepting part of the reference base map is to reduce memory overhead. The center point geographic coordinates of I Loc are recorded as

[0043] (7) In a block manner, Fourier Mellin transform is used to match Loc and . This method can achieve stable and accurate matching results in the presence of differences in ground objects, light and color differences between images, and the purpose of blocking is to reduce memory overhead and enable parallel processing to improve running speed. The specific implementation steps are as follows:

[0044] Downsample Loc to the same resolution as to improve the accuracy of the matching result. Resample Loc into 6x6 blocks, and record the center point geographic coordinates of each block as where m=1,2,…6, k=1,2,…6;

[0045] The difference between the center point geographic coordinates of Loc and is considered as the approximate coordinate offset between the two images, and the geographic coordinates of the possible homonym on can be calculated​ The formula is as follows:

[0046]

[0047] by Centered on the same size as the individual blocks in step 51, in Capture a section of the image, such as Figure 5 As shown.

[0048] right and The two corresponding image blocks are subjected to Fourier-Melin transform. In the first phase correlation operation, when calculating the rotation angle and scale factor, the scale factor threshold is set to 1.7. If it is greater than this threshold, the match is considered to have failed. In the second phase correlation operation, when calculating the coordinate translation, the response threshold is set to 0.05. If it is less than this threshold, the match is considered to have failed.

[0049] Get the matching points with the same name, i.e. The corresponding Same point on If the number of identical point pairs is less than 8, the matching is considered to have failed and the geometric fine correction ends; if the number of identical point pairs exceeds 16, the first 16 matching point pairs are selected as the final matching result by sorting these identical point pairs in descending order of response value.

[0050] (8) Solve for the parameters of the second-order polynomial based on pairs of corresponding points, using the following formula:

[0051]

[0052] Then, using the conventional indirect method and bilinear interpolation to pair I... Loc Resampling is used to achieve geometrical fine correction of the image. The original image and the result of geometrical fine correction are shown below. Figure 6 , Figure 7 As shown, the color image represents the reference base map, and the grayscale image represents the image before and after geometric correction. The processing results demonstrate that this invention can effectively perform local geometric correction on high-resolution remote sensing images.

[0053] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0054] The above embodiments only express several implementation manners of the present application, and the description is more specific and detailed, but it should not be understood as a limitation on the scope of the patent. It should be noted that for ordinary skilled persons in the art, without departing from the concept of the present application, several modifications and improvements can be made, which are within the scope of protection of the present application. Therefore, the protection scope of the patent of the present application should be subject to the appended claims.

Claims

1. A local geometric precise correction method suitable for high-resolution remote sensing images, characterized in that: The correction method is realized by the following steps: Step one, selecting the whole image I to be processed In According to the RPC coefficient of the whole image I In , the latitude and longitude range of the whole image I In is calculated; then the local area latitude and longitude range is set; The maximum value of the local area latitude and longitude coordinate of interest of the user is set as The minimum value is The maximum value of the ordinate The minimum value is The resolution of the whole image I In to be processed is r, the width W Loc and the height H Loc of the image I Loc of the local area of interest of the user are respectively: In the formula, the ceil function represents rounding up; Step two, according to the latitude and longitude range of the local region in step one, generate the local region image I by indirect resampling Loc Step three, traverse all the pixels of the local region image I Loc , and obtain the pixel value of all the pixels of I Loc . Obtaining I Loc The specific process of the pixel value of all pixels is as follows: Set the pixel coordinate of the current pixel as (i, j), according to the resolution r of the whole image and the minimum value of the longitude and latitude horizontal and vertical coordinates of the local region in step two, calculate the longitude and latitude of the center of the current pixel, X Loc (i,j) and Y Loc (i,j) are the horizontal and vertical coordinates of the longitude and latitude of the center of the current pixel, respectively, which are expressed by the following formula: Then according to the longitude and latitude, the RPC coefficient is used to inverse calculate the corresponding pixel point coordinate in I In , and the pixel value of the current pixel in I Loc is calculated to obtain the pixel value of all pixels in I Loc . Loc ​ Step three, calculate the reference base map I of the local area according to the latitude and longitude range of the local area in step one and the RPC coefficient Ref , the reference base map I Ref is obtained by taking a sliding window on the reference base map I Reference base map I of local area Ref Take sliding window, each time slide 0.1W in east-west direction Loc , slide 0.1H in south-north direction Loc , the size of sliding window is same as I Loc , so the maximum of horizontal coordinate of sliding window is , the minimum is , the maximum of vertical coordinate of sliding window is , the minimum is Get reference base map in sliding window The xsteps is the number of sliding from west to east, ysteps is the number of sliding from south to north; Step four, using perceptual hashing algorithm, respectively calculate local area image I Loc and the reference base map in the sliding window pHash fingerprint, namely: local area image I Loc and the reference base map in the sliding window After down-sampling to 100x100 pixels, DCT transform is performed, and the pHash fingerprint of the 10x10 area at the upper left corner of the image after DCT transform is intercepted; Step five, counting the local area image I Loc Hamming distance between the pHash fingerprint of the reference base image in all sliding windows, and according to the minimum reference base image image I obtained by intercepting the reference base image Step 6: Apply Fourier-Melin transform to I in a block-based manner. Loc and Perform a match, and then for I Loc Resampling enables geometrical correction of the image; Step six, I Loc is down-sampled to the same resolution, I is re-sampled, and I Loc is divided into 6x6 blocks, and the center point of each block is marked as where m = 1, 2, … 6, k = 1, 2, … 6. Step six two, the difference between the center point geographic coordinates of I Loc , the difference between the center point geographic coordinates of I the corresponding geographic coordinates of the same point on I the corresponding geographic coordinates of the same point on I the corresponding geographic coordinates of the same point on I Step six three, a block of video is cut from the center of the video in step six one, with the same size as the single block in step six one. Step six four, the block of video in step six three is cut into two blocks of video, with the same size as the single block in step six one. Step six five, the block of video in step six four is cut into two blocks Step six four, to With The corresponding two images, Fourier-Mellin transform, get the same name point, that is Corresponding to the same name point On If the number of the same name point pair is less than 8, it is considered that the matching fails, and the geometric precision correction ends; if the number of the same name point pair exceeds 16, the same name point pairs are sorted according to the response value from large to small, and the first 16 matching point pairs are taken as the final matching result. Step six, according to the homonymic point pair, the parameters of the second order polynomial are solved, and the formula is as follows: In the formula, a0-a5, b0-b5 are the coefficients of the polynomial.

2. The local geometric precise correction method for high-resolution remote sensing images according to claim 1, characterized in that: In step four, if the image contains multiple bands, it is converted to a gray image according to the following formula: Gray = 0.299R + 0.587G + 0.114B In the formula, Gray represents the gray image, and R, G and B represent the red, green and blue bands respectively. I Loc With DCT transform is performed to transform the image into the frequency domain; a 10x10 pixel region at the upper left corner of the image after DCT transform is intercepted, the mean value of the pixel values in the 10x10 pixel region at the upper left corner is calculated, and the region is binarized using the mean value as the threshold, that is, if greater than the mean value, it is recorded as 1, and if less than the mean value, it is recorded as 0, to obtain a 100-bit pHash fingerprint.

3. The local geometric precise correction method for high-resolution remote sensing images according to claim 1, characterized in that: In step five, the minimum reference base map with the minimum Hamming distance is calculated. And in the minimum reference base map Each of the objects expands outward by 0.1W. Loc Expand outwards by 0.1H in both the north and south directions. Loc This range will be used as the reference base map for the cut. The extracted reference base map The geographical coordinates of the center point are I Loc The geographical coordinates of the center point are

Images