Registration method for multi-temporal SAR remote sensing image

Abstract

The invention relates to the technical field of remote sensing image processing, and discloses a registration method for a multi-temporal SAR remote sensing image, which comprises the following steps: acquiring a to-be-processed multi-temporal SAR image, carrying out normalization processing on image pixels of the multi-temporal SAR image, and carrying out gamma correction on the multi-temporal SAR image; performing sub-region cutting on the corrected multi-temporal SAR image, and determining a multi-temporal scene SAR image of the same scene; the SAR-Harris detector is used for carrying out feature extraction on the SAR image of the multi-temporal scene to obtain feature points, and the SAR-Harris detector comprises an ROEWA operator and a multi-scale Harris corner detection operator; and carrying out feature matching on the feature points by adopting a rapid sampling consistency method so as to realize registration of the multi-temporal SAR image. According to the method, the registration precision and efficiency of the multi-temporal SAR remote sensing image can be effectively improved.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing image processing, and particularly to a registration method for multi-temporal SAR remote sensing images. Background Art

[0002] Synthetic Aperture Radar (SAR) has the advantages of two-dimensional high resolution, being unaffected by cloud and fog illumination, and operating all day and all weather, and can better adapt to harsh meteorological conditions and complex ground environments. In recent years, spaceborne SAR technology has developed rapidly, and SAR satellites have been able to provide various high-resolution, multi-temporal SAR data, with the highest resolution reaching 1 meter, providing various types of SAR datasets for urban areas. Multi-temporal SAR datasets have become an effective observation means in the fields of urban area planning, disaster prevention and control, land resource survey, etc.

[0003] In addition, the combination of multi-temporal SAR satellite data can effectively shorten the revisit period and improve the scene monitoring ability. When multiple SAR satellites form a constellation, a fixed area can be quickly revisited multiple times within a day, significantly enhancing the continuous observation ability of the SAR scene. In summary, the development of SAR satellites provides strong technical support for continuous tracking of scene changes.

[0004] In the actual application process, if dynamic monitoring of a scene is to be carried out, it is necessary to determine the position of the scene in multi-temporal SAR images and perform preprocessing means such as registration to facilitate subsequent change detection, target recognition and other steps. In terms of image registration, there are many classic feature extraction operators, such as the Scale Invariant Feature Transform (SIFT) method, the FAST (Features from Accelerated Segment Test) method, etc. These methods perform well in optical images, but are not suitable for multi-temporal SAR image scenes with speckle noise. The existing SAR-SIFT method replaces the area with speckle noise in the SAR image by detecting corner points in the image, and has achieved good results, but the registration efficiency of this method is generally average and often takes a long time in a large scene. Therefore, a fast and stable registration method still needs to be studied. The present application designs a registration method for multi-temporal SAR remote sensing images to achieve fast registration. Summary of the Invention

[0005] The purpose of the present invention is to provide a registration method for multi-temporal SAR remote sensing images, aiming to solve the above problems.

[0006] The present invention provides a registration method for multi-temporal SAR remote sensing images, including:

[0007] Obtain multi-temporal SAR images to be processed, perform normalization processing on the image pixels of the multi-temporal SAR images, and perform gamma correction on the multi-temporal SAR images;

[0008] Cut sub-regions from the corrected multi-temporal SAR images to determine multi-temporal scene SAR images of the same scene;

[0009] Use the SAR-Harris detector to extract features from the multi-temporal scene SAR images to obtain feature points. The SAR-Harris detector includes a ROEWA operator and a multi-scale Harris corner detection operator;

[0010] Adopt the random sample consensus method to perform feature matching on the feature points to achieve the registration of multi-temporal SAR images.

[0011] Preferably, before performing normalization processing on the image pixels of the multi-temporal SAR images, it includes:

[0012] Use the GDAL library to obtain the geotransformation information of the multi-temporal SAR images to be processed, and return a Dataset object. Read the pixel data and metadata of the multi-temporal SAR images. The pixel data and metadata include the position of the coordinate origin, the horizontal and vertical resolutions of the pixels; obtain the width, height, and number of bands of the multi-temporal SAR images;

[0013] Judge whether the origin in the geotransformation information is the default value (0, 0), and whether the horizontal resolution is the default value (1). If both the origin and the horizontal resolution are default values, do not process the multi-temporal SAR images.

[0014] Preferably, cutting sub-regions from the corrected multi-temporal SAR images to determine multi-temporal scene SAR images of the same scene includes:

[0015] Traverse the predefined scene information to detect whether the multi-temporal SAR images are in the defined dictionary;

[0016] Use the lonlat2geo and geo2imagexy functions to convert the longitude and latitude coordinates of the scene into pixel coordinates on the multi-temporal SAR images, and save the upper left vertex coordinates and upper right coordinates of the scene;

[0017] Process the pixel coordinates of the scene and convert them into integer values;

[0018] Determine whether the pixel coordinates of the scene are within the image range of the multi-temporal SAR images;

[0019] If the scene is within the image range, write the scene location information into a txt file, read the multi-temporal SAR images, convert the multi-temporal SAR images into 8-bit images, read the scene location information to crop the scene area of the multi-temporal SAR images, and save them to a new directory;

[0020] Add the information of the multi-temporal SAR images to the img_dict dictionary and execute the cropping program. If the image size of the multi-temporal SAR images meets the condition of exceeding the sliding window size, crop them according to the size of the sliding window and save the cropping results;

[0021] If the scene is not within the image range, continue to process the next scene. If no scene information meeting the target is found within the image range, print a prompt message, add the image file name to the img_dict["None"] list, and continue to process the next scene;

[0022] Until all processes are completed, print the img_dict dictionary. The img_dict dictionary contains a list of image file names for each scene, and print and output each scene and its corresponding number of images.

[0023] Preferably, the defined dictionary includes: area_info represents a dictionary containing scene information, the key represents the scene name, and the value represents a list containing longitude and latitude information.

[0024] Preferably, the formula model of the ROEWA operator is:

[0025] Vertical gradient:

[0026]

[0027] Horizontal gradient:

[0028]

[0029] Among them, f(x,y) = exp(-(|x| + |y|) / α) is a weighted exponential function, α is a parameter of the weighted exponential function, I(a + x, b + y) represents the image function, a and b are the center point coordinate positions, and x and y are the horizontal and vertical directions.

[0030] Preferably, when the regional weighting ratio R of direction i i,α is:

[0031]

[0032] Take the logarithm of R 1,a and R 3,a Then the new horizontal gradient G x,α and vertical gradient G of the ROEWA operatory,α It is defined as:

[0033] G x,α = log(R 1,α ), G y,α = log(R 3,α );

[0034] The gradient magnitude Mag α and the direction angle Ori α are calculated as follows:

[0035]

[0036] In the above formula, the exponentially weighted smoothing parameter α is the scale value. When calculating the image gradient, different α values are constructed to filter the image to establish the scale space.

[0037] Preferably, when using the multi-scale Harris corner detection operator, the scale space for feature point detection is constructed, and the exponentially weighted mean ratio is used as the calculation method of the gradient to obtain the SAR-Harris matrix C SH and the corner response function R SH as follows:

[0038]

[0039] R SH (x, y, α) = det(C SH (x, y, α)) - d · tr(C SH (x, y, α)) 2 ;

[0040] where G x,α , G y,α are the new horizontal gradient and vertical gradient of the ROEWA operator, is the Gaussian filtering function of the parameter ; according to the SAR-Harris matrix C SH the corner response function R SH can be calculated;

[0041] According to the SAR-Harris matrix C SH and the corner response function R SH the SAR-Harris function value of each pixel point is obtained, the extreme point detection map is obtained, and the scale space for extreme point detection at different scales is generated.

[0042] Preferably, in the scale space, the scale relationship between each layer of images is α i = α0 · k i , where α0 = 2 and k = 2 1 / 3, where \(i\) is the number of levels in the scale space, and \(i = 6\).

[0043] Preferably, the fast sample consensus method is used to perform feature matching on the feature points, including:

[0044] Adopting the nearest neighbor matching principle, based on two parameters with different precisions, screening the feature points for matching to obtain a set of rough matching points and a set of fine matching points;

[0045] Randomly select pairs of points that meet the model transformation requirements from the set of fine matching points, and calculate the affine transformation matrix;

[0046] Calculate the estimated coordinate error between the estimated point position and the actual matching point position in the set of rough matching points based on the affine transformation matrix;

[0047] Screen the pairs of points corresponding to the estimated coordinate error within the allowable error range;

[0048] After a set number of iterations, form a set with the screened pairs of points and calculate the root mean square error to complete the registration process.

[0049] Preferably, the set of rough matching points contains the set of fine matching points.

[0050] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0051] (1) The present invention provides a multi-temporal scene SAR image preprocessing scheme. By extracting the geographic transformation information in the SAR image, converting the longitude and latitude information into the projection coordinates in the image, thereby initially positioning the position of the scene in the multi-temporal SAR remote sensing image, and through normalization processing, correction, and sub-region segmentation, a multi-temporal SAR image of the same area is constructed;

[0052] (2) The present invention provides a registration method for multi-temporal SAR image scenes, designing a feature extraction and feature matching method based on invariant feature points, so as to quickly and accurately complete the registration of scene SAR images under the same scene. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the following drawings are only the embodiments of the present invention, and those of ordinary skill in the art can also obtain other drawings according to the provided drawings without creative efforts.

[0054] Figure 1 is a flowchart of a registration method for a multi-temporal SAR remote sensing image of the present invention;

[0055] Figure 2 It is a schematic flowchart of the multi-temporal SAR image processing of the present invention. Specific embodiments

[0056] Next, the technical solutions in the embodiments of the present application will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present application.

[0057] As Figure 1 shown, the present invention provides a registration method for multi-temporal SAR remote sensing images, including:

[0058] Obtain multi-temporal SAR images to be processed, normalize the image pixels of the multi-temporal SAR images, and perform gamma correction on the multi-temporal SAR images;

[0059] Cut sub-regions from the corrected multi-temporal SAR images to determine multi-temporal scene SAR images of the same scene;

[0060] Use the SAR-Harris detector to extract features from the multi-temporal scene SAR images to obtain feature points. The SAR-Harris detector includes the ROEWA operator and the multi-scale Harris corner detection operator;

[0061] Adopt the Random Sample Consensus (RANSAC) method to perform feature matching on the feature points to achieve the registration of multi-temporal SAR images.

[0062] The present invention is first a method for locating scene positions in multi-temporal synthetic aperture radar (SAR) images. This method utilizes the prior knowledge of the longitude and latitude coordinates of known scenes, extracts the geotransformation information in the SAR remote sensing images, and converts the longitude and latitude coordinates into projection coordinates in the images, thereby locating the approximate range of the scene. Read the multi-temporal SAR images to be processed, then normalize the image pixels, and then select appropriate gamma correction parameters for gamma correction according to the image mean to facilitate subsequent processing. Then, cut sub-regions according to the obtained projection coordinates, so as to locate and cut out the images of the same scene in multiple SAR remote sensing images taken at different times. After more refined alignment and correction, the data of the multi-temporal scene SAR images can be obtained for subsequent continuous observation of the same scene.

[0063] The registration method for multi-temporal scene SAR images is divided into two stages: feature extraction and feature matching. By using the points with invariant features in classical images, the affine transformation matrix between images is estimated through feature point extraction and matching, so as to complete the registration of other images to the reference image with one image as the standard.

[0064] The present invention can effectively improve the registration accuracy and efficiency of multi-temporal SAR remote sensing images. By normalizing and gamma correcting the multi-temporal SAR images, the brightness difference and noise influence between images of different time phases can be reduced, thus improving the accuracy of subsequent feature extraction and matching. At the same time, using the SAR-Harris detector for feature extraction can extract more stable and reliable feature points, further improving the stability of registration. In addition, adopting the fast sample consensus method for feature matching can quickly and accurately find the corresponding relationship between images of different time phases, thus realizing the fast registration of multi-temporal SAR images.

[0065] In some embodiments of the present application, before normalizing the image pixels of the multi-temporal SAR images, it includes: using the GDAL library to obtain the geotransformation information of the multi-temporal SAR images to be processed, and returning a Dataset object, reading the pixel data and metadata of the multi-temporal SAR images, where the pixel data and metadata include the position of the coordinate origin, the horizontal and vertical resolutions of the pixels; obtaining the width, height and number of bands of the multi-temporal SAR images; judging whether the origin in the geotransformation information is the default value (0, 0) and whether the horizontal resolution is the default value (1), if both the origin and the horizontal resolution are default values, then the multi-temporal SAR images are not processed.

[0066] It can be understood that after obtaining the geotransformation information and pixel data of the multi-temporal SAR images, by judging the rationality of the geotransformation information, the processing of invalid or incorrect image data can be effectively avoided, thus saving computing resources and improving processing efficiency. At the same time, the acquisition of the width, height and number of bands of the images provides a necessary information basis for subsequent processing. In addition, through the verification of the geotransformation information, the accuracy and reliability of the image data can be ensured, further improving the stability and accuracy of the registration of multi-temporal SAR images.

[0067] In some embodiments of the present application, for cutting sub-regions of the corrected multi-temporal SAR images to determine the multi-temporal scene SAR images of the same scene, it includes: traversing the pre-defined scene information to detect whether the multi-temporal SAR images are in the defined dictionary;

[0068] Use the lonlat2geo and geo2imagexy functions to convert the longitude and latitude coordinates of the scene into pixel coordinates on the multi-temporal SAR image, and save the coordinates of the upper left vertex and the upper right vertex of the scene;

[0069] Process the pixel coordinates of the scene and convert them to integer values;

[0070] Determine whether the pixel coordinates of the scene are within the image range of the multi-temporal SAR image;

[0071] If the scene is within the image range, write the scene location information to a txt file, read the multi-temporal SAR image, convert the multi-temporal SAR image into an 8-bit image, and read the scene location information to crop the scene area of the multi-temporal SAR image and save it to a new directory;

[0072] Add the information of the multi-temporal SAR image to the img_dict dictionary and execute the cropping program. If the image size of the multi-temporal SAR image meets the condition of exceeding the sliding window size, crop it according to the size of the sliding window and save the cropping result;

[0073] If the scene is not within the image range, continue to process the next scene. If no scene information that meets the target is found within the image range, print a prompt message, add the image file name to the img_dict["None"] list, and continue to process the next scene;

[0074] Until all processes are completed, print the img_dict dictionary. The img_dict dictionary contains a list of image file names for each scene, and print the output of each scene and its corresponding number of images.

[0075] The defined dictionary includes: area_info represents a dictionary containing scene information, where the keys represent scene names and the values represent a list containing longitude and latitude information.

[0076] It is understandable that by traversing the predefined scene information and detecting whether the multi-temporal SAR images are in the defined dictionary, the present invention can efficiently determine the scenes to be processed. Using the lonlat2geo and geo2imagexy functions to convert the longitude and latitude coordinates of the scene into pixel coordinates on the multi-temporal SAR images not only improves the accuracy of coordinate conversion but also provides accurate positioning information for subsequent cropping work. By performing integer value processing on the pixel coordinates of the scene, the feasibility of the cropping operation is ensured. At the same time, it is judged whether the pixel coordinates of the scene are within the image range of the multi-temporal SAR images, avoiding invalid cropping and further saving computing resources. After the scene is determined, the present invention writes the scene location information into a txt file and reads the multi-temporal SAR images for scene area cropping. This process is not only highly automated but also the cropping result is accurate. The cropped multi-temporal SAR images are saved to a new directory for subsequent processing and analysis. At the same time, the information of the multi-temporal SAR images is added to the img_dict dictionary, providing convenience for subsequent processing. If the image size of the multi-temporal SAR images meets the condition of exceeding the sliding window size, then cropping is performed according to the size of the sliding window. This step not only improves the processing efficiency but also ensures the uniformity and consistency of the cropping result.

[0077] If the scene is not within the image range or no scene information meeting the target is found within the image range, the present invention continues to process the next scene or prints a prompt message and adds the image file name to the img_dict["None"] list. This process not only improves the flexibility of processing but also ensures that all scenes are effectively processed. Until all processes are completed, the img_dict dictionary is printed. This dictionary contains a list of image file names for each scene, providing important reference information for subsequent analysis and processing. In addition, the number of images corresponding to each scene is printed and output, enabling the user to intuitively understand the data volume of each scene and further improving the practicality and reliability of the present invention.

[0078] In some embodiments of the present application, the formula model of the ROEWA operator is:

[0079] Vertical direction gradient:

[0080]

[0081] Horizontal direction gradient:

[0082]

[0083] Among them, \(f(x,y)=\exp(-(|x| + |y|) / \alpha)\) is a weighted exponential function, \(\alpha\) is the parameter of the weighted exponential function, \(I(a + x,b + y)\) represents the image function, \(a,b\) are the coordinates of the center point, and \(x,y\) are in the horizontal and vertical directions.

[0084] In some embodiments of the present application, when the regional weighting ratio \(R\) in the direction \(i\) i,α is:

[0085]

[0086] Taking the logarithm of \(R\) 1,a and \(R\) 3,a The new horizontal gradient \(G\) x,α and vertical gradient \(G\) y,α of the ROEWA operator are defined as:

[0087] \(G\) x,α =\(\log(R\) 1,α ), \(G\) y,α =\(\log(R\) 3,α );

[0088] The gradient magnitude \(Mag\) α and direction angle \(Ori\) α are calculated as:

[0089]

[0090] In the above formula, the exponential weighted smoothing parameter \(\alpha\) is a scale value. When calculating the image gradient, different \(\alpha\) values are constructed to filter the image and establish a scale space.

[0091] It can be understood that by calculating the gradient of multi-temporal SAR images using the ROEWA operator, the edge features in the images can be effectively extracted, providing important feature information for subsequent registration work. The ROEWA operator introduces a weighted exponential function to weight the pixels around the center point, making the gradient calculation more accurate. At the same time, by adjusting the parameters of the weighted exponential function, the weighted range can be flexibly controlled to adapt to image features of different scales. When calculating the horizontal gradient and vertical gradient, logarithmic processing is performed on the parameters of the weighted exponential function and the coordinates of the center point, further improving the accuracy and stability of the gradient calculation. The calculation formulas for the gradient magnitude and direction angle are derived based on the new definition of the ROEWA operator, providing accurate gradient information for subsequent feature matching and registration work. In addition, by constructing different scale values to filter the image and establish a scale space, the ROEWA operator can adapt to image data of different resolutions, improving the applicability and accuracy of the registration method.

[0092] In some embodiments of the present application, referring to the second moment of the multi-scale Harris detection algorithm, the SAR-Harris detector used has the same original gradient calculation method, but uses fewer scale space levels to improve the matching efficiency. When using the multi-scale Harris corner detection operator, a scale space for feature point detection is constructed, and the exponentially weighted mean ratio is used as the gradient calculation method to obtain the SAR-Harris matrix and corner response function as follows:

[0093]

[0094] R SH (x,y,α) = det(C SH (x,y,α)) - d·tr(C SH (x,y,α)) 2 ;

[0095] where G x,α and G y,α are the new horizontal and vertical gradients of the ROEWA operator, is the Gaussian filtering function of the parameter ; according to the SAR-Harris matrix C SH the corner response function R SH can be calculated;

[0096] According to the SAR-Harris matrix C SH and the corner response function R SH the SAR-Harris function value of each pixel point is obtained, an extreme point detection map is obtained, and a scale space for extreme point detection at different scales is generated.

[0097] In some embodiments of the present application, in the scale space, the scale relationship between each layer of images is α i = α0·k i , where α0 = 2, k = 2 1 / 3 , i is the number of layers of the scale space, and i = 6. i is reduced from the original 8 - 10 to 6, thereby reducing redundant calculations and improving efficiency. Then, by setting a certain threshold t SH for the scale space function, edges and low-contrast points are filtered out, and feature points can be preliminarily extracted.

[0098] It is understandable that by adopting a multi-scale Harris corner detection operator, the corner features in SAR remote sensing images can be detected more effectively, improving the accuracy and stability of corner detection. At the same time, using the exponential weighted mean ratio as the calculation method for gradients can better adapt to the speckle noise characteristics of SAR images, reducing the impact of noise on gradient calculation, thereby improving the accuracy of corner detection. In addition, by constructing a scale space for detecting extreme points at different scales, the effective extraction and matching of different scale features in multi-temporal SAR remote sensing images can be achieved, further improving the accuracy and robustness of image registration.

[0099] In the feature matching stage, the main task is to complete the feature point detection work, and the specific implementation method is as follows.

[0100] 1. First, read the reference image img1 and the image to be registered img2, and perform image preprocessing, the calculation method of which is the same as that proposed in the invention content section. Convert the two images img1 and img2 into grayscale images, and then divide the pixel values in the grayscale images by 255 so that their value ranges become [0, 1], and save them to gray1 and gray2 respectively.

[0101] 2. Call the custom build_scale() function to create a SAR-Harris scale space, the calculation method of which is the same as that proposed in the invention content section. The specific steps are as follows

[0102] (1) First, initialize the variables. Arrays for storing the SAR corner function, gradients, and angles are created, and the size information of the input image is obtained.

[0103] (2) For each scale i, calculate the parameters at the current scale according to the given parameters, such as the blur radius and the convolution kernel size.

[0104] (3) Generate a Gaussian kernel weight matrix. A two-dimensional Gaussian weighted matrix W is generated according to the current scale for blurring the image.

[0105] (4) Calculate the SAR response function: Blur the image using the Gaussian kernel, and then calculate the SAR response function according to the SAR corner detection algorithm.

[0106] (5) Calculate the gradients: Calculate the gradients based on the blurred image and use logarithmic operations to avoid numerical instability.

[0107] (6) Calculate the angles: Calculate the angles based on the x and y components of the gradients and convert them to the range between 0 and 360 degrees.

[0108] (7) Calculate the Harris matrix: Calculate the Harris matrix required for corner detection based on the gradient information.

[0109] (8) Gaussian weighted smoothing Harris matrix: Smooth the Harris matrix using another Gaussian kernel to reduce the influence of noise.

[0110] (9) Calculate the SAR corner function: Calculate the SAR corner function based on the values of the Harris matrix, including the product of diagonal elements, the sum of diagonal elements, and the product of a constant term.

[0111] (10) Return the result: Return the calculated SAR-Harris corner function, gradient, and angle.

[0112] 3. Perform extreme value detection on the SAR-Harris function. Traverse each layer of the image pyramid, detect the feature points with higher response values in each layer, and save these feature points in an array.

[0113] (1) Initialize parameters: Set parameters such as the boundary width, SIFT direction feature peak ratio, etc., and initialize the number of key points and the key point array.

[0114] (2) Traverse the scale space: Traverse the SAR corner function at each scale.

[0115] (3) Traverse the image pixels: Traverse each pixel point, skipping the boundary region.

[0116] (4) Check for local maximum: For each pixel, check if it is a local maximum. If the condition is met, it is a candidate key point.

[0117] (5) Calculate the direction histogram: Accumulate the gradient information in the area around the candidate key point to obtain the direction histogram.

[0118] (6) Determine the main direction: Find the peak in the direction histogram to determine the main direction of the key point.

[0119] (7) Add key points: Store the position, scale, scale level, and main direction information of the determined key points in the key point array.

[0120] (8) Organize the key point array: Delete the blank initial rows and return the valid key point array.

[0121] 4. Finally, generate the feature descriptor by calculating the key point array and the gradient and angle information in the neighborhood. Call the custom function calc_descriptors() to calculate the descriptors of the feature points. This function will traverse each feature point, extract the gradient and angle information in the surrounding area centered on this point, and then calculate the descriptor of this feature point. The specific steps are as follows

[0122] (1) Initialize parameters: Set descriptor parameters, including the number of bins in the circular region and the number of bins in the descriptor histogram.

[0123] (2) Obtain the number of key points: Obtain the number M of key points from the key point array.

[0124] (3) Traverse key points: Traverse each key point.

[0125] (4) Extract key point information: Obtain the position, scale, scale level, and main direction information of the key points from the key point array.

[0126] (5) Obtain gradient and angle information at the current scale: Extract the gradient and angle at the current scale from the gradient and angle information according to the scale level of the key point.

[0127] (6) Calculate the descriptor: Calculate the descriptor of the region around the key point according to the descriptor formula.

[0128] (7) Store the descriptor: Store the calculated descriptor into the descriptor array.

[0129] (8) Return the result: Return the calculated descriptor array and the key point position information.

[0130] In some embodiments of the present application, the feature matching of the feature points is performed by using the Random Sample Consensus (RANSAC) method, including: adopting the nearest neighbor matching principle, screening the feature points based on two parameters with different precisions to obtain a set of rough matching points and a set of fine matching points; randomly selecting pairs of points that meet the model transformation requirements in the set of fine matching points, and calculating the affine transformation matrix; calculating the estimated coordinate error between the estimated point position and the actual matching point position in the set of rough matching points based on the affine transformation matrix; screening the pairs of points corresponding to the estimated coordinate error within the allowable error range; after a set number of iterations, forming a set with the screened pairs of points and calculating the root mean square error to complete the registration process.

[0131] The set of rough matching points contains the set of fine matching points.

[0132] It is understandable that by combining the strategies of rough matching and fine matching, the accuracy and efficiency of feature matching are improved. First, the nearest neighbor matching principle is used to screen feature points under different precision parameters, effectively reducing the possibility of incorrect matching. Second, the affine transformation matrix is calculated by randomly selecting the number of point pairs that meet the requirements of model transformation, providing a basis for subsequent coordinate error calculation. Then, the estimated coordinate error is calculated based on the affine transformation matrix, and the number of point pairs within the error tolerance range is screened, further improving the matching accuracy. Finally, through a set number of iterations, the matching result is gradually optimized until a satisfactory root mean square error is achieved, thus completing the registration process. This method not only improves the accuracy of multi-temporal SAR remote sensing image registration, but also enhances the stability and reliability of registration, providing strong support for subsequent image analysis and applications.

[0133] In the feature matching stage, the RANSAC method is used to complete the feature matching process.

[0134] First, the nearest neighbor matching principle is used to match and screen the feature point descriptors obtained in the previous step with two parameters of different precisions, obtaining a set of low-precision matching points C_big and a set of high-precision matching points C_small. The high-precision set C_small should be included in C_big.

[0135] Next, randomly select the number of point pairs (such as 3 pairs) that meet the requirements of model transformation from the high-precision matching point set C_small, and use these point pairs to fit the affine transformation model to obtain the current transformation model H_i.

[0136] Use this transformation model H_i to calculate the coordinates of the matching points in the reference image img1 mapped to the image to be registered img2 in the C_big dataset, that is, estimate_points = H_i * img1_points.

[0137] Calculate the estimated coordinate error diff = |img2_points - estimate_points|, and count the number of point pairs m that do not exceed the specified error threshold (such as 1).

[0138] Repeat the selection of matching point pairs and calculate the error. After the loop reaches the specified number of times, output the matching point pair with the largest number, as well as the corresponding transformation model H, calculate the root mean square error, and complete the registration process.

[0139] Figure 2This is a demonstration case provided by the present invention. Two multi-temporal SAR remote sensing images of the same area are selected. First, the GDAL library of Python is used to obtain the geotransformation information of the images. The function of converting longitude and latitude coordinates to projected coordinates is used to obtain the projected coordinates of the scene on the image using the longitude and latitude information of the known scene. Then, image normalization processing, gamma correction, and positioning and cropping are performed to obtain the scene image. Finally, image slicing is completed through a sliding window with a size of 1024*1024, so as to obtain the preprocessing results of multi-temporal SAR images in the same area.

[0140] After that, the SAR-Harris detector is used for feature extraction, and the matching point pairs within the threshold range are obtained through the screening of random sample consensus, so as to complete the registration of multi-temporal slices. In the pictures of this group of demonstration experiments, a total of 45 matching points are detected, and the root mean square error is 0.732, indicating that the method provided by the present invention can quickly and stably complete the registration of multi-temporal scene SAR images.

[0141] Those skilled in the art should understand that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0142] The present application is described with reference to the flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each flow and / or block in the flowchart and / or block diagram can be implemented by computer program instructions, and the combination of the flows and / or blocks in the flowchart and / or block diagram can also be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate means for realizing the functions specified in Figure 1 one or more of the flows Figure 1 or blocks or multiple blocks.

[0143] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including instruction means, and the instruction means realizes the functions specified in Figure 1 one or more of the flows Figure 1 or blocks or multiple blocks.

[0144] These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus, so that a series of operation steps are performed on the computer or other programmable apparatus to generate a computer-implemented process, thereby providing instructions for implementing the steps of the function specified in one process or a plurality of processes and / or blocks Figure 1 in one block or a plurality of blocks Figure 1 in a block or a plurality of blocks.

[0145] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the above embodiments, those of ordinary skill in the art should understand that: modifications or equivalent replacements can still be made to the specific implementation manners of the present invention, and any modifications or equivalent replacements that do not depart from the spirit and scope of the present invention shall be covered by the protection scope of the claims of the present invention.

Claims

1. A registration method for multi-temporal SAR remote sensing images, characterized in that, Including: Obtain multi-temporal SAR images to be processed, normalize the image pixels of the multi-temporal SAR images, and perform gamma correction on the multi-temporal SAR images; Cut sub-regions from the corrected multi-temporal SAR images to determine multi-temporal scene SAR images of the same scene; Use the SAR-Harris detector to extract features from the multi-temporal scene SAR images to obtain feature points, and the SAR-Harris detector includes the ROEWA operator and the multi-scale Harris corner detection operator; Adopt the Random Sample Consensus (RANSAC) method to perform feature matching on the feature points to achieve the registration of multi-temporal SAR images.

2. The registration method for multi-temporal SAR remote sensing images according to claim 1, characterized in that, Before normalizing the image pixels of the multi-temporal SAR images, it includes: Use the GDAL library to obtain the geotransformation information of the multi-temporal SAR images to be processed, and return the Dataset object, read the pixel data and metadata of the multi-temporal SAR images, where the pixel data and metadata include the position of the coordinate origin, the horizontal and vertical resolutions of the pixels; obtain the width, height, and number of bands of the multi-temporal SAR images; Judge whether the origin in the geotransformation information is the default value (0, 0), and whether the horizontal resolution is the default value (1). If both the origin and the horizontal resolution are the default values, do not process the multi-temporal SAR images.

3. The registration method for multi-temporal SAR remote sensing images according to claim 2, characterized in that, Cut sub-regions from the corrected multi-temporal SAR images to determine multi-temporal scene SAR images of the same scene, including: Traverse the predefined scene information to detect whether the multi-temporal SAR images are in the defined dictionary; Use the lonlat2geo and geo2imagexy functions to convert the longitude and latitude coordinates of the scene into pixel coordinates on the multi-temporal SAR images, and save the upper-left vertex coordinates and upper-right coordinates of the scene; Process the pixel coordinates of the scene and convert them into integer values; Determine whether the pixel coordinates of the scene are within the image range of the multi-temporal SAR images; If the scene is within the image range, write the scene position information into a txt file, read the multi-temporal SAR images, convert the multi-temporal SAR images into 8-bit images, and read the scene position information to crop the scene area of the multi-temporal SAR images and save them to a new directory; Add the information of the multi-temporal SAR images to the img_dict dictionary and execute the cropping program. If the image size of the multi-temporal SAR images meets the condition of exceeding the sliding window size, crop them according to the size of the sliding window and save the cropping results; If the scene is not within the image range, continue to process the next scene. If no target-compliant scene information is found within the image range, print a prompt message, add the image file name to the img_dict["None"] list, and continue to process the next scene; Until all processes are completed, print the img_dict dictionary. The img_dict dictionary contains a list of image file names for each scene, and print and output each scene and its corresponding number of images.

4. The registration method for multi-temporal SAR remote sensing images according to claim 3, characterized in that, The defined dictionary includes: area_info represents a dictionary containing scene information, where the keys represent scene names and the values represent a list containing longitude and latitude information.

5. The registration method for multi-temporal SAR remote sensing images according to claim 1, characterized in that The formula model of the ROEWA operator is as follows: Vertical direction gradient: Horizontal direction gradient: Among them, f(x, y) = exp(-(|x| + |y|) / α) is a weighted exponential function, α is a parameter of the weighted exponential function, I(a + x, b + y) represents the image function, a and b are the center point coordinate positions, and x and y are the horizontal and vertical directions.

6. The registration method for multi-temporal SAR remote sensing images according to claim 5, wherein When the regional weighting ratio R in the direction i i,α is as follows: For R 1,α and R 3,α take the logarithm, then the new horizontal gradient G x,α and vertical gradient G y,α are defined as: G x,α = log(R 1,α ), G y,α = log(R 3,α ); Gradient magnitude Mag α and direction angle Ori α The calculation formula is as follows: In the above formula, the exponential weighted smoothing parameter α is a scale value. When calculating the image gradient, different α values are constructed to filter the image to establish a scale space.

7. The registration method for multi-temporal SAR remote sensing images according to claim 6, characterized in that, When using the multi-scale Harris corner detection operator, a scale space for feature point detection is constructed, and the exponentially weighted mean ratio is used as the calculation method for the gradient to obtain the SAR-Harris matrix C SH and the corner response function R SH as follows: R SH (x, y, α) = det(C SH (x, y, α)) - d·tr(C SH (x, y, α)) 2 ; Among them, G x,α and G y,α are the new horizontal gradient and vertical gradient of the ROEWA operator, is the parameter of the Gaussian filtering function; According to the SAR-Harris matrix C SH the corner response function R SH can be calculated; According to the SAR-Harris matrix C SH and the corner response function R SH The SAR-Harris function value of each pixel is obtained, an extreme point detection map is obtained, and a scale space for extreme point detection at different scales is generated.

8. The registration method for multi-temporal SAR remote sensing images according to claim 7, characterized in that In the scale space, the scale relationship between each layer of images is α i = α0·k i , where α0 = 2 and k = 2 1 / 3 , i is the number of layers in the scale space, and i = 6.

9. The registration method for multi-temporal SAR remote sensing images according to claim 1, characterized in that, The fast sample consensus method is used to perform feature matching on the feature points, including: Adopting the nearest neighbor matching principle, based on two parameters with different precisions, the feature points are screened for matching to obtain a set of rough matching points and a set of fine matching points; Randomly select the number of point pairs that meet the model transformation requirements from the set of fine matching points, and calculate the affine transformation matrix; Based on the affine transformation matrix, calculate the estimated coordinate error between the estimated point position and the actual matching point position in the set of rough matching points; Screen the number of point pairs corresponding to the estimated coordinate error within the allowable error range; After a set number of iterations, the screened number of point pairs is combined to calculate the root mean square error, completing the registration process.

10. The registration method for multi-temporal SAR remote sensing images according to claim 9, wherein The set of rough matching points contains the set of fine matching points.

Images