An optical and SAR image registration method based on high-level ground feature compensation

Abstract

This invention discloses a method for registering optical and SAR images based on high-level ground feature compensation, comprising: Step 1, acquiring SAR image data and extracting potential dual-echo information from the SAR image to calculate the high-level ground feature compensation factor; Step 2, acquiring optical image data, resampling the optical image data, and calculating the high-level ground feature compensation factor using the morphological building index; Step 3, processing the SAR image data using the high-level ground feature compensation factor obtained in Step 1, and processing the optical image data using the high-level ground feature compensation factor obtained in Step 2, and extracting a threshold-adaptive feature point set; Step 4, performing feature point matching operation on the threshold-adaptive feature point set. This invention can effectively reduce the registration error caused by high-level ground features in local areas, optimize the registration effect of optical image data and SAR image data, and improve the accuracy of remote sensing images.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and in particular to an optical and SAR image registration method based on high-level ground feature compensation. Background Technology

[0002] In applications such as urban dynamic monitoring and disaster emergency rescue based on remote sensing technology, optical and synthetic aperture radar (SAR) imagery offer excellent data complementarity. Specifically, optical imagery possesses rich spectral, texture, and spatial detail information, which helps in the detailed characterization of geographic objects; SAR imagery, on the other hand, has all-weather, all-day imaging capabilities, and its diverse scattering characteristics can reflect information such as the height and spatial structure of geographic objects. Therefore, the joint processing of optical and SAR imagery has become a research hotspot in the field of remote sensing imagery.

[0003] The prerequisite for combining optical and SAR imagery is effective registration of heterogeneous images. Currently, commonly used remote sensing image registration methods can be divided into gray-scale-based and feature-based methods based on the different image information being registered. The former mainly involves statistically analyzing the gray-scale information of the images to determine the spatial transformation relationship between the two images. The main drawback of this method is its sensitivity to changes in gray-scale information, making it difficult to meet this requirement between heterogeneous remote sensing images with significant imaging differences. Feature-based methods, on the other hand, directly utilize features such as corner points, edges, and textures of the images for matching, and finally estimate the registration parameters. They are more robust to changes in image gray-scale and are therefore more suitable for registering heterogeneous remote sensing images. Currently, various methods proposed by the academic community do not consider the adaptive enhancement problem of control point pairs near tall features such as buildings in SAR images; they only consider global information and not local information. In fact, the difficulty in registration lies in the correction of tall features, which is also the main source of registration error. This is because tall features are more susceptible to projection differences due to the imaging angle, or differences in remote sensing data acquisition conditions may cause tall features to not be completely consistent in the images, resulting in positional differences. Therefore, it is necessary to design specific strategies for the local areas where high-rise features are located.

[0004] To enhance the ability to extract feature points from tall-rise features, it is first necessary to identify the potential local areas of these features. In optical imagery, the morphology building index (MBI), currently very popular, is based on establishing the relationship between the implicit features of buildings and morphological operator attributes, thus serving as a direct basis for identifying tall-rise features. However, due to the presence of speckle noise in SAR imagery, which results in a loss of shape information, and the fact that the spectral characteristics of SAR imagery are less rich than those of multispectral data, MBI is difficult to apply directly to SAR imagery. On the other hand, the unique dual-echo characteristic of SAR imagery reflects the secondary reflection phenomenon of radar beams after illuminating the ground and walls, providing a feasible approach for identifying features of a certain height. Nevertheless, an effective dual-echo detection method for SAR imagery is currently lacking. Summary of the Invention

[0005] To reduce the registration error caused by high-rise features in local areas and further optimize the registration effect of remote sensing images, this invention proposes an optical and SAR image registration method based on high-rise feature compensation.

[0006] An optical and SAR image registration method based on high-level ground feature compensation includes the following steps:

[0007] Step 1: Acquire SAR image data, extract potential double echo information from SAR images, and calculate the compensation factor for high-rise ground features in SAR images based on the potential double echo information.

[0008] Step 2: Acquire optical image data, resample the optical image data, and use the morphological building index to calculate the high-rise ground feature compensation factor of the optical image data.

[0009] Step 3: Process SAR image data using the SAR image high-level ground feature compensation factor obtained in Step 1, process optical image data using the optical image high-level ground feature compensation factor obtained in Step 2, and extract threshold adaptive feature point set.

[0010] Step 4: Perform feature point matching operation on the threshold adaptive feature point set generated in Step 3;

[0011] Step 4 specifically includes the following steps;

[0012] Step 4.1: Use the nearest neighbor ratio algorithm to extract the initial set of matching point pairs;

[0013] Step 4.2, calculate the distance differences dx and dy between pairs of corresponding points;

[0014] Step 4.3, calculate the angle θ between the matching point pairs, θ = arctan(dx / dy);

[0015] Step 4.4: Use K-means clustering to analyze the angle θ in step 4.3, retain the largest cluster, and complete the registration.

[0016] Furthermore, the specific steps for extracting potential dual-echo information from SAR images in step 1 are as follows:

[0017] Step 1.1: Use the Otsu method to determine the threshold T, perform threshold segmentation, and retain pixels with gray values ​​greater than T;

[0018] Step 1.2: Edge detection is performed using the Canny operator;

[0019] Step 1.3: Extract straight line features using Hough transform and record the extraction results;

[0020] Line segments that satisfy the following two conditions are defined as constituting double echo information;

[0021] (1) The formal intersection point is the point where more than 10 curves intersect in the polar coordinate system;

[0022] (2) The maximum pixel interval of the line segment is 1.

[0023] Considering that building double echoes typically appear as prominent bright lines, step 1 first uses the Otsu method to determine the threshold T, retaining pixels with gray values ​​greater than T, which can effectively filter out weak edge information in the image. Furthermore, selecting only points where the formal intersections are more than 10 curves intersecting in the polar coordinate system, and line segments with a maximum pixel interval of 1, as double echo information can further exclude features with similar morphological characteristics, such as roads.

[0024] Furthermore, the formula for calculating the compensation factor for high-altitude features in SAR images in step 1 is as follows:

[0025]

[0026] Wherein, CFHGS represents the compensation factor for high-altitude ground features in SAR images, and X and Y represent the number of dual-echo pixels in the 4-neighborhood and 8-neighborhood, respectively.

[0027] At this point, the larger the CFHGS, the higher the probability that there are buildings in the local area. CFHGS provides an effective basis for judgment for subsequent image registration work.

[0028] Furthermore, the formula for calculating the morphological building index in step 2 is as follows:

[0029]

[0030] Where MBI represents the building index, DMP represents the differential morphological profile, D represents the number of directions, and S represents the scale of the linear structural element. s max s min These represent the maximum and minimum values ​​of the linear structuring element scale, respectively. d represents the direction of the linear structuring element, s represents the scale of the linear structuring element, and W-TH represents morphological white-hat reconstruction.

[0031] The Morphological Building Index (MBI) extracts buildings by establishing the relationship between implicit features and morphological operators. It is a mature technical index that provides a reliable foundation for the calculation of CFHGO.

[0032] Furthermore, the formula for calculating the optical image high-level ground feature compensation factor in step 2 is as follows:

[0033] CFHGO=1+(M+X′+Y′) / 3

[0034] Wherein, CFHGO represents the optical image high-level ground feature compensation factor, M represents the morphological index pixel, and X' and Y' are the mean values ​​of the MBI index in the 4-neighborhood and 8-neighborhood of M, respectively.

[0035] Similarly, the larger the CFHGO, the higher the probability of the presence of buildings in the local area. CFHGO provides an effective basis for judgment in subsequent image registration work.

[0036] Furthermore, the extraction step of the threshold adaptive feature point set in step 3 specifically includes:

[0037] Step 3.1: Based on the scale-invariant feature transform algorithm, the multiple sets of optical images resampled in Step 2 are convolved using a scale-variable Gaussian function G(x,y,σ) to construct a Gaussian scale space, expressed as follows:

[0038] L(x,y,σ)=G(x,y,σ)*I(x,y)

[0039] Where L represents the Gaussian convolution kernel, I represents the gray value, (x,y) represents the spatial coordinates, and σ represents the scale coordinates;

[0040] Step 3.2: Subtract the images of adjacent scales from the same group in Step 3.1 to obtain the Gaussian difference scale space, expressed as follows:

[0041] D(x,y,σ)=L(x,y,kσ)-L(x,y,σ)

[0042] Where D represents the Gaussian difference scale space, and k represents the scaling factor;

[0043] Step 3.3: Compare the point in the Gaussian difference scale space obtained in Step 3.2 with its eight neighborhood points and nine points in each of the two adjacent scales above and below. If the point is an extreme point, it is determined as a candidate feature point.

[0044] Step 3.4: Multiply the point in the Gaussian difference scale space obtained in step 3.2 by the corresponding compensation factor, and compare it with its eight neighborhood points and nine points in each of the two adjacent scales above and below. If the point is an extreme point, it is determined as a candidate feature point.

[0045] Step 3.5: Obtain the candidate feature points from Steps 3.3 and 3.4, which together form a candidate feature point set;

[0046] Step 3.6: Filter out unstable feature points from the candidate feature point set, expressed as follows:

[0047]

[0048]

[0049] In this context, CFHG represents CFHGS in SAR imagery and CFHGO in optical imagery. The index is set to 1.05 in SAR imagery and 1.5 in optical imagery. D(x,y,σ) represents the response value of the feature point in the Gaussian difference scale space. The threshold value represents the extreme point detection threshold.

[0050] Step 3.7: Calculate the principal curvature of each point in the candidate feature point set in Step 3.5, remove unstable edge key points, and obtain the final threshold adaptive feature point set.

[0051] Furthermore, the K value in the K-means clustering is set to 6.

[0052] Scale-invariant feature transform (SIFT) is commonly used in image processing. SIFT can detect key points in an image and is a local feature descriptor. It is based on points of interest in the local appearance of an object and is independent of the size and rotation of the image. It also has a high tolerance for changes in lighting, noise, and micro-viewpoints. Based on this, SIFT is highly significant, easy to obtain, and has minimal error, which can provide effective support for the registration of remote sensing images in this patent. In addition, due to the significant differences in imaging mechanisms, geometric features, and radiometric features between optical and SAR images, there are a large number of mismatches in the initial matching point pairs. This patent, based on the randomness of mismatches and the certain relationship between correct matching point pairs, calculates the distance difference between corresponding point pairs, and then calculates the angle between matching point pairs. K-means clustering is used to analyze the angles, retaining the largest cluster, which can minimize mismatches.

[0053] Beneficial effects:

[0054] This invention can effectively reduce the registration error caused by high-altitude features in local areas, further optimize the registration effect of optical image data and SAR image data, and significantly improve the accuracy of remote sensing images. Attached Figure Description

[0055] To more clearly illustrate the technical solution and effects of the present invention, the images used in the accompanying drawings retain the original colors of the optical images from the actual experimental data.

[0056] Figure 1 This is a flowchart of the image registration method of the present invention;

[0057] Figure 2 The optical images are from the first set of experimental data in the embodiment;

[0058] Figure 3 SAR images from the first set of experimental data in the embodiment;

[0059] Figure 4 Optical images from the second set of experimental data in the embodiment;

[0060] Figure 5 The SAR imagery is from the second set of experimental data in the example. Detailed Implementation

[0061] To enable those skilled in the art to better understand the present invention, the present invention will now be further described in conjunction with the accompanying drawings and specific embodiments.

[0062] Example

[0063] The basic experimental environment for this embodiment is a CPU-2.20GHz (Intel(R) Core(TM) i7-8750), RAM-8GB, and MATLAB version 2020b.

[0064] The experimental data used in this embodiment consists of two sets of remote sensing images of Sendai, Japan. Each set of experimental data includes one high-resolution optical image and one SAR image. The acquisition time for both is March 2011. The optical image is an optical image acquired by the IKONOS satellite, and the SAR image is a SAR image acquired by the TerraSAR-X satellite. In order to make the heterogeneous remote sensing images have the same resolution, the optical image was resampled in this experiment to obtain a set of heterogeneous images with the same resolution.

[0065] In the image registration process, three key parameters are mainly involved: image magnification, initial extreme point detection threshold, and Euclidean distance threshold. In this embodiment, the aforementioned three parameters of the two datasets are all set to (1, 0.03, 0.90).

[0066] like Figure 1 As shown, an optical and SAR image registration method based on high-level ground feature compensation includes the following steps.

[0067] Step 1: Acquire SAR image data collected by TerraSAR-X satellite, extract potential double echo information from SAR images, and calculate the high-altitude ground feature compensation factor of SAR images based on the potential double echo information.

[0068] The specific steps for extracting potential dual-echo information from SAR images in step 1 are as follows:

[0069] Step 1.1: Use the Otsu method to determine the threshold T, perform threshold segmentation, and retain pixels with gray values ​​greater than T;

[0070] Step 1.2: Edge detection is performed using the Canny operator;

[0071] Step 1.3: Extract straight line features using Hough transform and record the extraction results;

[0072] Line segments that satisfy the following two conditions are defined as constituting double echo information;

[0073] (1) The formal intersection point is the point where more than 10 curves intersect in the polar coordinate system;

[0074] (2) The maximum pixel interval of the line segment is 1.

[0075] The formula for calculating the compensation factor for high-altitude features in SAR images in step 1 is as follows:

[0076]

[0077] Wherein, CFHGS represents the compensation factor for high-altitude ground features in SAR images, and X and Y represent the number of dual-echo pixels in the 4-neighborhood and 8-neighborhood, respectively.

[0078] Step 2: Acquire optical image data collected by IKONOS satellite, resample the optical image data, and calculate the high-rise ground feature compensation factor of the optical image data using the morphological building index.

[0079] The formula for calculating the morphological building index in step 2 is as follows:

[0080]

[0081] Where MBI represents the building index, DMP represents the differential morphological profile, D represents the number of directions, and S represents the scale of the linear structural element. s max s min These represent the maximum and minimum values ​​of the linear structuring element scale, respectively. d represents the direction of the linear structuring element, s represents the scale of the linear structuring element, and W-TH represents morphological white-hat reconstruction.

[0082] The formula for calculating the optical image high-level ground feature compensation factor in step 2 is as follows:

[0083] CFHGO=1+(M+X′+Y′) / 3

[0084] Wherein, CFHGO represents the optical image high-level ground feature compensation factor, M represents the morphological index pixel, and X' and Y' are the mean values ​​of the MBI index in the 4-neighborhood and 8-neighborhood of M, respectively.

[0085] Step 3: Process SAR image data using the SAR image high-level ground feature compensation factor obtained in Step 1, process optical image data using the optical image high-level ground feature compensation factor obtained in Step 2, and extract threshold adaptive feature point set.

[0086] Step 3, the extraction of the threshold adaptive feature point set, specifically includes:

[0087] Step 3.1: Based on the scale-invariant feature transform algorithm, the multiple sets of optical images resampled in Step 2 are convolved using a scale-variable Gaussian function G(x,y,σ) to construct a Gaussian scale space, expressed as follows:

[0088] L(x,y,σ)=G(x,y,σ)*I(x,y)

[0089] Where L represents the Gaussian convolution kernel, I represents the gray value, (x,y) represents the spatial coordinates, and σ represents the scale coordinates;

[0090] Step 3.2: Subtract the images of adjacent scales from the same group in Step 3.1 to obtain the Gaussian difference scale space, expressed as follows:

[0091] D(x,y,σ)=L(x,y,kσ)-L(x,y,σ)

[0092] Where D represents the Gaussian difference scale space, and k represents the scaling factor;

[0093] Step 3.3: Compare the point in the Gaussian difference scale space obtained in Step 3.2 with its eight neighborhood points and nine points in each of the two adjacent scales above and below. If the point is an extreme point, it is determined as a candidate feature point.

[0094] Step 3.4: Multiply the point in the Gaussian difference scale space obtained in step 3.2 by the corresponding compensation factor, and compare it with its eight neighborhood points and nine points in each of the two adjacent scales above and below. If the point is an extreme point, it is determined as a candidate feature point.

[0095] Step 3.5: Obtain the candidate feature points from Steps 3.3 and 3.4, which together form a candidate feature point set;

[0096] Step 3.6: Filter out unstable feature points from the candidate feature point set, expressed as follows:

[0097]

[0098]

[0099] In this context, CFHG represents CFHGS in SAR imagery and CFHGO in optical imagery. The index is set to 1.05 in SAR imagery and 1.5 in optical imagery. D(x,y,σ) represents the response value of the feature point in the Gaussian difference scale space. The threshold value represents the extreme point detection threshold.

[0100] Step 3.7: Calculate the principal curvature of each point in the candidate feature point set in Step 3.5, remove unstable edge key points, and obtain the final threshold adaptive feature point set.

[0101] Step 4: Perform feature point matching operation on the threshold adaptive feature point set generated in Step 3.

[0102] Step 4 specifically includes the following steps;

[0103] Step 4.1: Use the nearest neighbor ratio algorithm to extract the initial set of matching point pairs;

[0104] Step 4.2, calculate the distance differences dx and dy between pairs of corresponding points;

[0105] Step 4.3, calculate the angle θ between the matching point pairs, θ = arctan(dx / dy);

[0106] Step 4.4: Analyze the angle θ in Step 4.3 using K-means clustering. Set the K value in K-means clustering to 6 and retain the largest cluster to complete the registration.

[0107] To quantitatively evaluate the registration effect of this embodiment, this paper uses matching rate (MR), structural similarity (SSIM), and matching accuracy (CMR) for accuracy evaluation.

[0108]

[0109] Here, SMPP represents the filtered matching pairs, and MPP represents the initial matching pairs. The matching rate (MR) reflects the proportion of filtered matching pairs among the initial matching pairs.

[0110]

[0111] Where, μ x μ y and σ corresponds to the mean and variance of the two images, respectively. xy Let c1 and c2 represent the covariance, and c2 be constants. Structural similarity reflects the degree of image distortion in terms of brightness, contrast, and structure; the smaller the value, the greater the difference between the two images.

[0112]

[0113] Where, N o N represents the number of feature points on an optical image. s The CMR value represents the number of feature points on the SAR image. The matching accuracy reflects the reliability of the image matching algorithm. The higher the CMR value, the more accurate the transformation matrix is. When the SMPP is less than 4, the registration is generally considered to have failed.

[0114] This embodiment is compared with three existing registration strategies. The first existing strategy is SIFT (Scale-invariant feature transform), which uses Gaussian convolution kernels to construct initial extreme points in scale space, determines the location of feature points through stability judgment, and generates feature point descriptors using local gradient directions. The second existing strategy is SURF (Speeded Up Robust Features), which uses integral images for image convolution, uses an approximate Hessian matrix to detect extreme points, and constructs descriptors based on Haar wavelet features. The third existing strategy is proposed in the literature "A heterogeneous high-resolution remote sensing image registration method for high-rise features [J]. Electronic Measurement Technology, 2021, 44(16):12.", which uses shadow detection and image segmentation to screen high-rise features, extracts feature points based on phase consistency adaptive threshold, and introduces shadow area weighted feature vector distance and transformation error penalty factor to improve registration accuracy.

[0115] The table below compares the registration results of this embodiment with the three registration methods described above for the first set of experimental data.

[0116] Initial matching point pairs Matching points after filtering MR SSIM CMR This embodiment 51 20 39.2 0.0442 2×20 / (2084+1129) Strategy 1 51 14 27.5 0.0437 2×14 / (2096+1115) Strategy 2 30 3 10.0 / 2×3 / (1431+1129) Strategy 3 19 5 26.3 0.0437 2×5 / (754+575)

[0117] The table below compares the registration results of this embodiment with the three registration methods described above for the second set of experimental data.

[0118] Initial matching point pairs Matching points after filtering MR SSIM CMR This embodiment 47 18 38.3 0.0268 2×18 / (2220+1260) Strategy 1 51 15 29.4 0.0247 2×15 / (2153+1271) Strategy 2 30 3 10.0 / 2×3 / (1428+1143) Strategy 3 21 5 23.8 0.0248 2×5 / (786+648)

[0119] As shown in the table above, based on the extracted feature points, it can be seen that the matching point pairs obtained in the datasets of the two sets of experimental data in this embodiment have higher stability and better registration effect than the other three existing strategies.

[0120] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. An optical and SAR image registration method based on high-layer feature compensation, characterized in that, The method comprises the following steps: Step 1, acquiring SAR image data, extracting potential double-bounce information of the SAR image, and calculating a high-rise ground feature compensation factor of the SAR image according to the potential double-bounce information of the SAR image; Step 2, acquiring optical image data, performing resampling processing on the optical image data, and calculating the resampling result of the optical image data by using a morphological building index to obtain a high-rise ground feature compensation factor of the optical image data; Step 3, processing the SAR image data by using the high-rise ground feature compensation factor of the SAR image obtained in step 1, processing the optical image data by using the high-rise ground feature compensation factor of the optical image obtained in step 2, and extracting a threshold adaptive feature point set; Step 4, performing a feature point matching operation on the threshold adaptive feature point set generated in step 3; The step 4 specifically comprises the following steps: Step 4.1, extracting an initial matching point pair set by using a nearest neighbor ratio algorithm; Step 4.2, calculating distance differences dx and dy of the same point pair; Step 4.3, calculating an angle θ between the matching point pairs, θ = arctan(dx / dy); Step 4.4, analyzing the angle θ in step 4.3 by using K-means clustering, retaining the largest class, and completing registration; The calculation formula of the high-rise ground feature compensation factor of the SAR image in the step 1 is CFHGS = 1 + (X + Y) / 2 Wherein, CFHGS represents the high-rise ground feature compensation factor of the SAR image, and X and Y represent the number of double-bounce pixels in the 4-neighborhood and 8-neighborhood of the double-bounce pixels respectively; The calculation formula of the high-rise ground feature compensation factor of the optical image in the step 2 is CFHGO = 1 + (M + X' + Y') / 3 2. The optical and SAR image registration method based on high layer object compensation according to claim 1, characterized in that, Wherein, CFHGO represents the high-rise ground feature compensation factor of the optical image, M represents a morphological index pixel, and X' and Y' are the 4-neighborhood and 8-neighborhood MBI index mean values of M respectively. The specific steps for extracting the potential double-bounce information of the SAR image in the step 1 are as follows: Step 1.1, determining a threshold T by using the OTSU method, performing threshold segmentation, and retaining pixels with a gray value greater than T; Step 1.2, performing edge detection by using a Canny operator; Step 1.3, extracting straight line features by using a Hough transform, and recording the extraction result; The line segments satisfying the following two conditions constitute double-bounce information; (1) The formal intersection point is a point where 10 or more curves intersect in a polar coordinate system; 3. The method according to claim 1, wherein the method is characterized by: (2) The maximum pixel interval of the line segment is 1. where MBI represents building index, DMP represents differential morphological profile, D represents the number of directions, S represents the scale of linear structuring element, s max , s min respectively represent the maximum and minimum of the scale of linear structuring element, d represents the direction of linear structuring element, s represents the scale of linear structuring element, and W-TH represents morphological white hat reconstruction.

4. The optical and SAR image registration method based on high layer object compensation according to claim 1, characterized in that, The calculation formula of the morphological building index in the step 2 is The extraction steps of the threshold adaptive feature point set in the step 3 specifically comprise: Step 3.1, based on a scale-invariant feature transform algorithm, performing convolution operation on the multiple groups of resampled optical images in step 2 by using a scale-variable Gaussian function G(x, y, σ) to construct a Gaussian scale space, and the expression is L(x, y, σ) = G(x, y, σ) * I(x, y) Wherein, L represents a Gaussian convolution kernel, I represents a gray value, (x, y) represents a spatial coordinate, and σ represents a scale coordinate; Step 3.2, difference the images of the same group of adjacent scales in step 3.1 to obtain a Gaussian difference scale space, expressed as D(x, y, σ) = L(x, y, kσ) - L(x, y, σ) Wherein, D represents the Gaussian difference scale space, k represents the scale factor; Step 3.3, compare the points in the Gaussian difference scale space obtained in step 3.2 with the eight neighborhood points and the nine points of the adjacent upper and lower scales, and if the point is an extreme point, it is determined as a candidate feature point; Step 3.4, multiply the points in the Gaussian difference scale space obtained in step 3.2 by the corresponding compensation factor, compare with the eight neighborhood points and the nine points of the adjacent upper and lower scales, and if the point is an extreme point, it is determined as a candidate feature point; Step 3.5, obtain the candidate feature points in step 3.3 and step 3.4 to form a candidate feature point set; Step 3.6, filter out unstable feature points in the candidate feature point set, expressed as Wherein, CFHG represents CFHGS in the SAR image, CFHGO in the optical image, index is set as 1.05 in the SAR image, index is set as 1.5 in the optical image, and D(x, y, sigma) represents the response value of the feature point in the Gaussian difference scale space, represents an extreme point detection threshold value; Step 3.7, calculate the principal curvatures of each point in the candidate feature point set in step 3.5, remove unstable edge key points, and obtain a final threshold adaptive feature point set.

5. The method according to claim 1, wherein the method is characterized by: The K value in the K-means clustering is set to 6.

Images