A multi-source remote sensing image registration method based on anisotropic diffusion description

Abstract

The present invention discloses a multi-source remote sensing image registration method based on anisotropic diffusion description, including: for the input heterogeneous image pairs, using different methods to perform gradient calculation; establishing anisotropic scale spaces for the heterogeneous images respectively, retaining the image edges while filtering out noise, and establishing Harris scale spaces on this basis to perform sub-region segmentation and feature point extraction within the sub-region; dividing the circular neighborhood of the feature into three parts in the radial direction and into eight equally spaced angles in the angular direction, resulting in a total of 24 sub-regions. The mapped diffusion function is used to describe the features within the sub-regions, and the feature descriptor is obtained after standardization, and the transformation parameter θ is then estimated. final This paper uses the diffusion function in the nonlinear diffusion equation to describe the feature points and their neighborhood to obtain the feature descriptor, which can reduce the radiation characteristic differences and noise interference in heterogeneous images and enhance the independence and robustness of the feature descriptor.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing image registration, and in particular relates to a multi-source remote sensing image registration method based on anisotropic diffusion description. Background Art

[0002] SAR (Synthetic Aperture Radar) images offer all-day, all-weather coverage, long range, and wide coverage, while optical images are more easily interpreted by human vision. Therefore, the fusion of optical and SAR images effectively combines the advantages of different image sensors. Image registration, a key step in the image fusion process, aligns two images of the same area captured at different times, angles, and weather conditions.

[0003] In the existing technology, image registration methods are mainly divided into two categories: region-based methods and feature-based methods. Specifically, region-based registration methods are mainly based on the similarity measurement of the same regions in the image, but are limited by the differences in imaging mechanisms. Feature-based registration methods mainly use point, line and region features to find the mapping relationship between two images to achieve registration.

[0004] The main processing steps of feature-based registration methods include feature extraction, feature representation, feature matching, and model building. Among them, feature representation can convert the geometric information of the feature neighborhood in the image, including gradient information, structural information, frequency domain information, etc., into a digital vector representation through mathematical statistics, which facilitates the subsequent matching between features with the same name. For the feature description in multi-source image registration, the following two properties need to be met: (1) independence. The descriptors between non-identical feature points should be highly separable so that the correct match can be obtained efficiently and accurately during feature matching; (2) robustness. For feature points with the same name in different source images, their feature descriptions need to be highly similar to avoid mismatching. Lowe et al. proposed the SIFT (Scale Invariant Feature Transform) descriptor, which uses the gradient modulus as the weight to calculate the gradient direction histogram within the square neighborhood of the feature point and obtain a one-dimensional vector descriptor of the feature point. Mikolajczyk et al. used the logarithmic polar coordinate system to replace the square area in SIFT to form GLOH (Gradient Location and Orientation Histogram), which obtained a more independent and robust feature descriptor.

[0005] It can be seen that the image registration methods in the existing technology all describe features based on image gradient information. However, due to differences in the radiation characteristics and noise interference of multi-source images, differences in image grayscale and gradient calculation results occur, which brings difficulties to the description of features and affects the independence and robustness of feature descriptors. Summary of the Invention

[0006] In order to solve the above problems existing in the prior art, the present invention provides a multi-source remote sensing image registration method based on anisotropic diffusion description. The technical problem to be solved by the present invention is achieved through the following technical solutions:

[0007] An embodiment of the present invention provides a multi-source remote sensing image registration method based on anisotropic diffusion description, comprising:

[0008] Acquire the optical image and SAR image to be registered;

[0009] After filtering the optical image and the SAR image respectively using a nonlinear diffusion equation, an anisotropic scale space ASS of the optical image is established based on the first gradient of the optical image. O , and establish an anisotropic scale space ASS based on the first gradient of the SAR image S ;

[0010] In the anisotropic scale space ASS O In the optical image, a first Hessian matrix is established point by point to generate a first Harris scale space, and in the anisotropic scale space ASS S In, establishing a second Hessian matrix for the SAR image point by point to generate a second Harris scale space;

[0011] Feature points are detected layer by layer in the first Harris scale space of the optical image and the second Harris scale space of the SAR image, and the main directions are assigned to the detected feature points, and then the feature vectors of the feature points are calculated;

[0012] Determine the reference image and the image to be registered from the optical image and the SAR image, and for the feature point P in the reference image O , calculate the Euclidean distance between it and each feature point in the image to be registered, and record the feature point pairs corresponding to the Euclidean distance that meet the first preset condition into the coarse matching point pair set CM;

[0013] After removing the erroneous points in the coarse matching point pair set CM, performing primary matching on the remaining feature points in the coarse matching point pair set CM to form a fine matching point pair set FM and determine the transformation parameter θ;

[0014] Generate a high confidence matching point pair set C based on the fine matching point pair set FM sample and low confidence matching point pair set C total , and for the high confidence matching point pair set C sample and the low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final ;

[0015] Based on the final matching point pair set FM final Estimate the transformation parameter θ between the optical image and the SAR image final , and according to the transformation parameter θ final Transform the image to be registered into the coordinate system of the reference image.

[0016] In one embodiment of the present invention, before the step of filtering the optical image and the SAR image using a nonlinear diffusion equation, the method further includes:

[0017] The first gradient of each target pixel in the optical image is calculated using the Sobel operator according to the following formula:

[0018]

[0019]

[0020] in, represent the magnitude and direction of the first gradient of each target pixel in the optical image, Respectively represent the first gradient value of each target pixel in the optical image in the horizontal direction and the vertical direction, is the intensity value image of the optical image after Gaussian smoothing, Respectively represent the templates of Sobel in the horizontal and vertical directions, β j is the scale of the optical image, j represents the scale layer in the scale space;

[0021] The Adaptive ROEWA operator is used to calculate the local exponentially weighted average ratio of each target pixel in the SAR image in the horizontal and vertical directions. The first gradient of each target pixel in the SAR image is calculated according to the local exponentially weighted average ratio in the horizontal and vertical directions according to the following formula:

[0022]

[0023]

[0024] in, Respectively represent the amplitude and direction of the first gradient of each target pixel in the SAR image, Respectively represent the local exponential weighted average ratio of each target pixel in the SAR image in the horizontal and vertical directions, Respectively represent the first gradient value of each target pixel in the SAR image in the horizontal and vertical directions, α i is the scale of the scale layer where the target point is located in the SAR image.

[0025] In one embodiment of the present invention, in the anisotropic scale space ASS O In the optical image, a first Hessian matrix is established point by point to generate a first Harris scale space, and in the anisotropic scale space ASS S The step of establishing a second Hessian matrix point by point for the SAR image to generate a second Harris scale space comprises:

[0026] In the anisotropic scale space ASS O The second gradient of the SAR image after filtering by the nonlinear diffusion equation is calculated, and the first Hessian matrix is established point by point in the second gradient corresponding to each scale layer:

[0027]

[0028] in, They represent the second gradient values of each target pixel in the horizontal and vertical directions in the SAR image after filtering by the nonlinear diffusion equation, Denotes the variance as σ i Gaussian kernel, σ i =α i / β i ,× represents convolution operation;

[0029] According to the first Hessian matrix M S (α i ) Generate the first Harris scale space:

[0030] R S (α i )=det(M S (α i ))-d·tr(M S (α i )) 2

[0031] Among them, det represents the determinant of the calculation matrix, d represents the corner detection factor, tr represents the trace of the calculation matrix, R S (α i) represents the scale sequence {a1,...,α n}Generated Harris scale space;

[0032] In the anisotropic scale space ASS S The second gradient of the optical image after filtering by the nonlinear diffusion equation is calculated, and the second Hessian matrix is established point by point in the second gradient corresponding to each scale layer:

[0033]

[0034] in, They represent the second gradient values of each target pixel in the horizontal and vertical directions in the optical R image after filtering by the nonlinear diffusion equation, Denotes the variance as β i Gaussian kernel of

[0035] According to the second Hessian matrix M O (β i ) Generate the second Harris scale space:

[0036] R O (β i )=det(M O (β i ))-d·tr(M O (β i )) 2

[0037] Among them, R O (β i ) represents the scale sequence {β1,...,β n}Generated Harris scale space.

[0038] In one embodiment of the present invention, the feature vector of the feature point is calculated according to the following steps:

[0039] A circular neighborhood is established with the feature point as the center and the preset length as the radius, and a logarithmic coordinate system is established with the feature point as the pole;

[0040] Divide the circular neighborhood into three parts along the radial direction, and divide the inner circle and the two outer rings into 8 equal parts along the chord direction, forming 24 sub-areas;

[0041] After establishing a rectangular coordinate system with the feature point as the origin and the main direction of the feature point as the positive direction of the horizontal axis, the feature point in each sub-area is rotated to the corresponding rectangular coordinate system;

[0042] Divide 180° into 8 equal parts, and use the mapped diffusion function as the weight to sum the main directions of the feature points in each sub-region to obtain a one-dimensional feature vector with a length of 192 corresponding to each feature point.

[0043] In one embodiment of the present invention, a reference image and an image to be registered are determined from the optical image and the SAR image, and a feature point P in the reference image is O , respectively calculating the Euclidean distance between it and each feature point in the image to be registered, and recording the feature point pairs corresponding to the Euclidean distances that meet the first preset condition into the coarse matching point pair set CM, including:

[0044] One of the optical image and the SAR image is used as the reference image and the other as the image to be registered. A feature point P is selected from the feature point set of the reference image. O , and calculate the feature point P O The Euclidean distance between the feature vector of and the feature vector of each feature point in the image to be registered;

[0045] When the minimum value in the Euclidean distance satisfies the first preset condition, the feature point pair corresponding to the minimum value is recorded in the coarse matching point pair set CM.

[0046] In one embodiment of the present invention, the first preset condition is:

[0047]

[0048] Wherein, ED1 and ED2 represent the minimum value and the second minimum value in the Euclidean distance respectively, and Thres is a preset threshold.

[0049] In one embodiment of the present invention, after removing erroneous points from the coarse matching point pair set CM, performing primary matching on the remaining feature points in the coarse matching point pair set CM to form a fine matching point pair set FM and determining the transformation parameter θ includes:

[0050] Select any two feature point pairs [P k ,Q k ]、[P l ,Q l ], respectively feature point pairs [P k ,Q k ]The corresponding Euclidean distance and feature point pair [P l ,Q l ] The ratio of the Euclidean distances D kl ;

[0051] After traversing all feature point pairs, kl Perform histogram statistics;

[0052] Eliminate the feature point pairs corresponding to the minimum value in the histogram, and calculate the root mean square error of the Euclidean distance of the remaining feature point pairs;

[0053] Detect whether the root mean square error meets the second preset condition; if not, return to select any two feature point pairs from the coarse matching point pair set CM [P k ,Q k ]、[P l ,Q l ] step; if so, the point set consisting of the remaining feature point pairs is used as the SC-CM point set pair;

[0054] The cascaded sample consistency estimation algorithm is used to perform first-level matching on the SC-CM point set pairs to obtain the fine matching point set pairs FM and transformation parameters θ.

[0055] In one embodiment of the present invention, the step of performing primary matching on the SC-CM point set pairs using the cascaded sample consistency estimation algorithm to obtain the fine matching point set pairs FM and the transformation parameters θ includes:

[0056] Randomly select three feature point pairs from the SC-CM point set and calculate the current transformation parameters;

[0057] transforming the feature points of the reference image in the coarse matching point pair set CM using the current transformation parameters, and calculating the error based on the feature points of the image to be registered in the coarse matching point pair set CM;

[0058] Counting the number num of feature point pairs that make the error smaller than a preset threshold;

[0059] Check whether the preset number of iterations has been reached; if not, return to the step of randomly selecting three feature point pairs from the SC-CM point set and calculating the current transformation parameters; if so, calculate the maximum number num of feature point pairs;

[0060] The feature point pairs corresponding to the maximum value of the feature point pair number num are obtained to form a fine matching point set pair FM, and the current transformation parameter is used as the transformation parameter θ.

[0061] In one embodiment of the present invention, the high confidence matching point pair set C sample and the low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final ; The steps include:

[0062] The feature points P of the reference image in the fine matching point pair set FM are I According to the transformation parameter θ, it is mapped to the image to be registered, and the nearest feature point Q is detected in the sub-area of the image to be registered. I ;

[0063] If so, then the feature point pair [P I ,Q I ] Recorded in the high confidence matching point set C sample and low confidence matching point pair set C total If not, then find the nearest feature point Q in the image to be registered J , and record the feature point pairs into the low confidence matching point pair set C total ;

[0064] The cascade sample consistency estimation algorithm is used to estimate the high confidence matching point pair set C sample and low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final .

[0065] Compared with the prior art, the present invention has the following beneficial effects:

[0066] An embodiment of the present invention provides a multi-source remote sensing image registration method based on anisotropic diffusion description, which uses the diffusion function in the nonlinear diffusion equation to describe feature points and their neighborhoods to obtain feature descriptors. It can reduce the radiation characteristic differences and noise interference in heterogeneous images, and enhance the independence and robustness of the feature descriptors.

[0067] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 This is a flow chart of a multi-source remote sensing image registration method based on anisotropic diffusion description provided by an embodiment of the present invention;

[0069] Figure 2 is another flow chart of a multi-source remote sensing image registration method based on anisotropic diffusion description provided by an embodiment of the present invention;

[0070] Figure 3a This is an example diagram of a SAR image provided by an embodiment of the present invention;

[0071] Figure 3b The embodiment of the present invention provides Figure 3a The gradient modulus of the SAR image shown;

[0072] Figure 3c The embodiment of the present invention provides Figure 3a Diffusion coefficient of the SAR image shown;

[0073] Figure 3d The embodiment of the present invention provides Figure 3a The diffusion function mapping result of the SAR image shown;

[0074] Figure 4 is a schematic diagram of calculating a characteristic vector provided by an embodiment of the present invention;

[0075] Figure 5a is an example diagram of an optical image provided by an embodiment of the present invention;

[0076] Figure 5b is another example diagram of a SAR image provided by an embodiment of the present invention;

[0077] Figure 6a The embodiment of the present invention provides Figure 5a The feature point detection results of the optical image shown;

[0078] Figure 6b The embodiment of the present invention provides Figure 5b The feature point detection results of the SAR image shown;

[0079] Figure 7 The embodiment of the present invention provides Figure 5a and Figure 5b Schematic diagram of feature matching results of optical-SAR image pairs shown;

[0080] Figure 8a is a schematic diagram of an image registration result provided by an embodiment of the present invention;

[0081] Figure 8b The embodiment of the present invention provides Figure 8a A local schematic diagram of the image registration results is shown. DETAILED DESCRIPTION

[0082] The present invention will be further described in detail below with reference to specific examples, but the embodiments of the present invention are not limited thereto.

[0083] Figure 1 This is a flow chart of a multi-source remote sensing image registration method based on anisotropic diffusion description provided by an embodiment of the present invention. Figure 1 As shown, an embodiment of the present invention provides a multi-source remote sensing image registration method based on anisotropic diffusion description, comprising:

[0084] S1. Acquire the optical image and SAR image to be registered;

[0085] S2. After filtering the optical image and SAR image respectively using the nonlinear diffusion equation, the anisotropic scale space ASS of the optical image is established based on the first gradient of the optical image. O , and establish anisotropic scale space ASS based on the first gradient of SAR image S ;

[0086] S3. In anisotropic scale space ASS O The first Hessian matrix is established point by point for the optical image to generate the first Harris scale space, and the first Harris scale space is generated in the anisotropic scale space ASS. S In the SAR image, the second Hessian matrix is established point by point to generate the second Harris scale space;

[0087] S4, performing feature point detection layer by layer in the first Harris scale space of the optical image and the second Harris scale space of the SAR image, assigning main directions to the detected feature points, and then calculating feature vectors of the feature points;

[0088] S5, determine the reference image and the image to be registered from the optical image and the SAR image, and for the feature point P in the reference image O , calculate the Euclidean distance between it and each feature point in the image to be registered, and record the feature point pairs corresponding to the Euclidean distance that meet the first preset condition into the coarse matching point pair set CM;

[0089] S6. After removing the erroneous points in the coarse matching point pair set CM, perform first-level matching on the remaining feature points in the coarse matching point pair set CM to form a fine matching point pair set FM and determine the transformation parameter θ;

[0090] S7. Generate a high confidence matching point pair set C based on the fine matching point pair set FM sample and low confidence matching point pair set C total , and match the high confidence point pair set C sample Match point pairs with low confidence set C total Perform secondary matching to form the final matching point pair set FM final ;

[0091] S8, based on the final matching point set FM final Estimate the transformation parameter θ between the optical image and the SAR image final , and according to the transformation parameter θ final Transform the image to be registered into the coordinate system of the reference image.

[0092] Figure 2 This is another flow chart of the multi-source remote sensing image registration method based on anisotropic diffusion description provided by an embodiment of the present invention. Figure 2 As shown, before the step of filtering the optical image and the SAR image respectively using the nonlinear diffusion equation, the method further includes:

[0093] The first gradient of each target pixel in the optical image is calculated using the Sobel operator according to the following formula:

[0094]

[0095]

[0096] in, represent the magnitude and direction of the first gradient of each target pixel in the optical image. To facilitate calculation, this embodiment simplifies the Sobel operator into a rectangular window. Respectively represent the first gradient value of each target pixel in the optical image in the horizontal direction and the vertical direction, is the intensity value image of the optical image after Gaussian smoothing, Respectively represent the templates of Sobel in the horizontal and vertical directions, β j is the scale of the optical image, j represents the scale layer in the scale space;

[0097] The Adaptive ROEWA operator is used to calculate the local exponentially weighted average ratio of each target pixel in the SAR image in the horizontal and vertical directions. The first gradient of each target pixel in the SAR image is calculated according to the local exponentially weighted average ratio in the horizontal and vertical directions according to the following formula:

[0098]

[0099]

[0100] in, Respectively represent the amplitude and direction of the first gradient of each target pixel in the SAR image, Respectively represent the local exponential weighted average ratio of each target pixel in the SAR image in the horizontal and vertical directions, Respectively represent the first gradient value of each target pixel in the SAR image in the horizontal and vertical directions, α i is the scale of the scale layer where the target point is located in the SAR image.

[0101] In the above formula, the local exponentially weighted average ratios of each target pixel in the SAR image in the horizontal and vertical directions are expressed as:

[0102]

[0103]

[0104] Among them, M and N represent the number of rows and columns of the calculation window respectively, and the window size is related to α i Proportional to, I(x,y) represents the image intensity at the coordinate (x,y) in the SAR image, Represents the local exponential weighted average of image intensity, loc = r, d, l, u, representing the four directions of up, down, left, and right respectively. Represents the adaptive constant values in the horizontal and vertical directions respectively.

[0105] It should be noted that, in this embodiment, Taking the logarithm to calculate the first gradient value of the target pixel in the SAR image in the horizontal and vertical directions can eliminate the interference of coherent speckle noise generated by the special imaging mechanism in the SAR imaging process on the subsequent image processing.

[0106] In addition, to ensure that the scale space of the optical image and the SAR image is consistent, when calculating the first gradient of the two, α i and β i Need to meet:

[0107] α i+1 / α i =k

[0108] β i+1 / β i =k

[0109] α1=β1

[0110] Among them, k represents the scale ratio between two adjacent scale layers in the scale space.

[0111] In the above step S2, the optical image and the SAR image are filtered using the nonlinear diffusion equation, as shown in the following equation:

[0112]

[0113] Among them, L represents the input image, i.e., optical image or SAR image, t represents the current scale, c(x, y, t) is the diffusion function, Δ, They represent the divergence and gradient calculation symbols respectively. div(·) represents the calculation of divergence. The diffusion function c(x, y, t) is a function of the image gradient modulus. Different calculation methods can be adopted depending on the focus of the image content, as follows:

[0114]

[0115] in, Represents the first gradient of the input image, and K is a constant.

[0116] For the nonlinear diffusion equation, the additive splitting operator can be used to solve it as follows:

[0117]

[0118] Finally, the anisotropic scale space ASS of the optical image is obtained O and the anisotropic scale space ASS of SAR images S .

[0119] In the above step S3, in the anisotropic scale space ASS O The first Hessian matrix is established point by point for the optical image to generate the first Harris scale space, and the first Harris scale space is generated in the anisotropic scale space ASS. S The steps of establishing the second Hessian matrix point by point for the SAR image to generate the second Harris scale space include:

[0120] S301, in anisotropic scale space ASS O The second gradient of the SAR image after filtering by the nonlinear diffusion equation is calculated, and the first Hessian matrix is established point by point in the second gradient corresponding to each scale layer:

[0121]

[0122] in, They represent the second gradient values of each target pixel in the horizontal and vertical directions in the SAR image after filtering by the nonlinear diffusion equation, Denotes the variance as σ i Gaussian kernel, σ i =α i / β i ,× represents convolution operation;

[0123] S302, according to the first Hessian matrix M S (α i ) Generate the first Harris scale space:

[0124] R S (α i )=det(M S (α i ))-d·tr(M S (α i )) 2

[0125] Among them, det represents the determinant of the calculation matrix, d represents the corner detection factor, which is generally 0.02 to 0.04, tr represents the trace of the calculation matrix, and R S (α i ) represents the scale sequence {a1,…,α n}Generated Harris scale space;

[0126] S303, in anisotropic scale space ASS SThe second gradient of the optical image after filtering by the nonlinear diffusion equation is calculated, and the second Hessian matrix is established point by point in the second gradient corresponding to each scale layer:

[0127]

[0128] in, They represent the second gradient values of each target pixel in the horizontal and vertical directions in the optical R image after filtering by the nonlinear diffusion equation, Denotes the variance as β i Gaussian kernel of

[0129] S304, according to the second Hessian matrix M O (β i ) Generate the second Harris scale space:

[0130] R O (β i )=det(M O (β i ))-d·tr(M O (β i )) 2

[0131] Among them, R O (β i ) represents the scale sequence {β1,...,β n}Generated Harris scale space.

[0132] It should be noted that the calculation process of the second gradient of the SAR image after filtering by the nonlinear diffusion equation is the same as the calculation method of the first gradient of the SAR image. The calculation process of the second gradient of the optical image after filtering by the nonlinear diffusion equation is the same as the calculation method of the first gradient of the optical image, so it is not repeated here.

[0133] In step S4 above, to ensure that the detected feature points are evenly distributed across the two images, this embodiment divides the two images into 5x5 regions of equal size before detection. Within each region, the target pixel with the largest corner response value is selected as a preliminary point. Next, maximum value detection is performed on each preliminary point, and the target pixel that passes the detection is considered a feature point. Furthermore, to ensure that feature points can be obtained even in flat areas, this embodiment also detects grid intersections.

[0134] Similarly, when calculating the principal direction of a feature point, the circular area surrounding the feature point is used as the neighborhood, and the radius of the circular area is proportional to the scale of the scale layer in which it is located. Due to the different imaging mechanisms of optical images and SAR images, the imaging results for the same target may show gradient direction flipping. Therefore, when calculating the principal direction of a feature point, the gradient directions of all feature points are normalized to [0°, 180°), and this interval is divided into 18 subintervals. After performing a gradient direction histogram statistics on the feature points in the neighborhood area, the resulting histogram is Gaussian smoothed and the maximum value of the interpolated histogram is obtained to obtain the principal direction of the feature point.

[0135] Specifically, in this embodiment, the feature vector of the feature point can be calculated according to the following steps:

[0136] A circular neighborhood is established with the feature point as the center and the preset length as the radius, and a logarithmic coordinate system is established with the feature point as the pole;

[0137] The circular neighborhood is divided into three parts along the radial direction, and the inner circle and the two outer rings are divided into 8 equal parts along the chord direction, forming 24 sub-areas;

[0138] After establishing a rectangular coordinate system with the feature point as the origin and the main direction of the feature point as the positive direction of the horizontal axis, the feature points in each sub-area are rotated to the corresponding rectangular coordinate system;

[0139] Divide 180° into 8 equal parts, and use the mapped diffusion function as the weight to sum the main directions of the feature points in each sub-region to obtain a one-dimensional feature vector with a length of 192 corresponding to each feature point.

[0140] Figure 3a is an example diagram of a SAR image provided by an embodiment of the present invention. Figure 3b The embodiment of the present invention provides Figure 3a The gradient modulus of the SAR image shown, Figure 3c The embodiment of the present invention provides Figure 3a The diffusion coefficient of the SAR image shown, Figure 3d The embodiment of the present invention provides Figure 3a The diffusion function mapping result of the SAR image shown in the figure is obtained by Figures 3a-3d It can be seen that compared with the gradient modulus, the diffusion function is clearer and more distinguishable in expressing the flat area and texture area of the SAR image.

[0141] Figure 4 is a schematic diagram of calculating the characteristic vector provided by the embodiment of the present invention. Figure 4 As shown, when calculating the eigenvector of a feature point, first take the feature point as the center and the length proportional to the current layer scale, such as 12σ iA circular neighborhood is established for the radius, and a logarithmic polar coordinate system is established with the feature point as the pole; then the circular neighborhood is divided into 24 sub-areas: specifically, along the radial direction according to 3σ i , 8σ i The circular neighborhood is divided into three parts, and the inner circle and the two outer rings are divided into 8 equal parts in terms of angle. In this way, the eigenvector calculation neighborhood divided into 24 sub-areas is obtained.

[0142] To ensure the invariance of each feature point to direction, a rectangular coordinate system is established with the main direction of the feature point as the positive direction of the horizontal axis of the rectangular coordinate system and the feature point as the origin. The coordinates of each pixel in the neighborhood of the feature point in the rectangular coordinate system are calculated. Gradient histogram statistics are then performed on the points in the sub-region: 180° is divided into 8 equal parts, and the mapped diffusion function is used as the weight and the gradient direction of each pixel is used as the index to sum them. Finally, for each feature point, a one-dimensional feature vector with a length of 192 is obtained. To reduce the sensitivity of the feature point to light intensity, the feature vector can be normalized after calculation. The large values in the feature vector caused by high intensity are limited and then normalized again to ensure the additive normalization of the feature vector.

[0143] In the above step S5, a reference image and an image to be registered are determined from the optical image and the SAR image, and the feature point P in the reference image is O , respectively calculating the Euclidean distance between it and each feature point in the image to be registered, and recording the feature point pairs corresponding to the Euclidean distances that meet the first preset condition into the coarse matching point pair set CM, including:

[0144] One of the optical image and the SAR image is used as the reference image and the other as the image to be registered. The feature point P is selected from the feature point set of the reference image. O , and calculate the feature point P O The Euclidean distance between the feature vector of and the feature vector of each feature point in the image to be registered;

[0145] When the minimum value in the Euclidean distance meets the first preset condition, the feature point pair corresponding to the minimum value is recorded in the coarse matching point pair set CM.

[0146] In this implementation, one of the optical image and the SAR image is set as the reference image and the other is set as the image to be registered. First, the feature points in the two images are roughly matched. Specifically, the feature point set in the reference image is denoted as O, and a feature point P is selected from it. O , then the feature point P OThe Euclidean distance between the feature vectors and the feature point set S of the image to be registered is calculated point by point. If the minimum and the second minimum values of the Euclidean distance meet the first preset condition shown below, the feature point pair corresponding to the minimum value is recorded in the coarse matching point pair set CM:

[0147]

[0148] Among them, ED1 and ED2 represent feature points P respectively. O The minimum Euclidean distance and the second minimum Euclidean distance between the feature vector of and the feature vectors of each feature point in the feature point set S, Thres is the preset threshold.

[0149] Furthermore, fine feature point matching is performed based on the coarse matching point pair set CM.

[0150] In the above step S6, after removing the erroneous points in the coarse matching point pair set CM, the steps of performing primary matching on the remaining feature points in the coarse matching point pair set CM to form a fine matching point pair set FM and determining the transformation parameter θ include:

[0151] S601, select any two feature point pairs from the coarse matching point pair set CM [P k ,Q k ]、[P l ,Q l ], respectively feature point pairs [P k ,Q k ]The corresponding Euclidean distance and feature point pair [P l ,Q l ] The ratio of the Euclidean distances D kl ;

[0152] S602, after traversing all feature point pairs, kl Perform histogram statistics;

[0153] S603, eliminating the feature point pairs corresponding to the minimum value in the histogram, and calculating the root mean square error of the Euclidean distance of the remaining feature point pairs;

[0154] S604, check whether the root mean square error meets the second preset condition; if not, return to select any two feature point pairs from the coarse matching point pair set CM [P k ,Q k ]、[P l ,Q l ] step; if so, the remaining feature point pairs are combined into a scale-constrained matching point pair set SC-CM;

[0155] S605 , performing first-level matching on the SC-CM point set pairs using the cascaded sample consistency estimation algorithm to obtain fine matching point set pairs FM and transformation parameters θ.

[0156] First, the scale constraint is applied to the coarse matching point pair set CM. Specifically, the scale constraint of any two feature point pairs in the coarse matching point pair set CM is calculated. k ,Q k ]、[P l ,Q l The ratio of the Euclidean distance between kl , where P k 、P l represents the feature points in the reference image, Q k , Q l Represents the feature points in the registered image. That is:

[0157]

[0158] After traversing all possible feature point pairs, ij Perform histogram statistics, the horizontal axis of the histogram is D ij The vertical axis represents the number of point pairs accumulated in each interval. Remove the feature point pair corresponding to the minimum value in the histogram and calculate the root mean square error (RMSE). Repeat the above steps until the RMSE changes within the preset range or until three feature point pairs remain.

[0159] After the scale constraint processing, the scale-constrained matching point pair set SC-CM is obtained. Then, the Fast Sample Consensus (FSC) algorithm is used for the first-level matching. The specific steps are as follows:

[0160] Randomly select three feature point pairs from the SC-CM point set and calculate the current transformation parameters;

[0161] The feature points of the reference image in the coarse matching point pair set CM are transformed using the current transformation parameters, and the error is calculated based on the feature points of the image to be registered in the coarse matching point pair set CM;

[0162] Count the number of feature point pairs num that make the error less than the preset threshold;

[0163] Check whether the preset number of iterations has been reached; if not, return to the step of randomly selecting three feature point pairs from the SC-CM point set and calculating the current transformation parameters; if so, calculate the maximum number of feature point pairs num;

[0164] The feature point pairs corresponding to the maximum value of the feature point pair number num are obtained to form a fine matching point set pair FM, and the current transformation parameter is used as the transformation parameter θ.

[0165] Furthermore, for the high confidence matching point pair set C sample Match point pairs with low confidence set Ctotal Perform secondary matching to form the final matching point pair set FM final ; The steps include:

[0166] The feature points P of the reference image in the fine matching point pair set FM are I According to the transformation parameter θ, it is mapped to the image to be registered, and the nearest feature point Q is detected in the sub-area of the image to be registered. I ;

[0167] If so, then the feature point pair [P I ,Q I ] Recorded in the high confidence matching point set C sample and low confidence matching point pair set C total ; If not, find the nearest feature point Q in the image to be registered J , and record the feature point pairs into the low confidence matching point pair set C total ;

[0168] The cascade sample consistency estimation algorithm is used to estimate the high confidence matching point pair set C sample and low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final .

[0169] According to the final matching point set FM final , combined with different image transformation models such as similarity transformation, affine transformation, projective transformation, etc., the transformation parameter θ between the optical image and the SAR image can be estimated final , according to the transformation parameter θ final The registration between the two images is completed by transforming the image to be registered into the coordinate system of the reference image through similarity transformation, projection transformation, and affine transformation. Generally, to examine and test the accuracy of the registration, the two registered images can be plotted together in the form of a checkerboard grid, and the registration effect can be reflected by comparing edges, regions, etc.

[0170] Figure 5a is an example diagram of an optical image provided by an embodiment of the present invention. Figure 5b This is another example of a SAR image provided by an embodiment of the present invention. In order to verify the performance of the descriptor used in the present invention, optical-SAR image pairs under different scenes were selected for registration experiments, and the SIFT and OS-SIFT algorithms were compared. Taking a single set of image pairs as an example, Figure 5a-5b As shown, the optical images used in this embodiment are collected from Google Earth, and the SAR images are collected from TerraSAR-X, both with a resolution of 3m.

[0171] Figure 6aThe embodiment of the present invention provides Figure 5a The feature point detection results of the optical image shown are Figure 6b The embodiment of the present invention provides Figure 5b The feature point detection results of the SAR image are shown in Figure 2. Figure 6a-6b As shown, since the present invention adopts a block extraction method, the number of features of the two images is relatively close. Figure 7 The embodiment of the present invention provides Figure 5a and Figure 5b The feature matching results for an optical-SAR image pair are shown in Figure 1. The descriptors used in this invention have improved independence and robustness, resulting in a higher number of matches and higher matching accuracy during feature matching. Table 1 compares the results of the three image registration methods, where CMR is the exact match rate, reflecting the stability of the algorithm, and RMSE is the root mean square error, reflecting the accuracy of the algorithm. Figure 8a is a schematic diagram of the image registration result provided by an embodiment of the present invention, Figure 8b The embodiment of the present invention provides Figure 8a The local schematic diagram of the image registration result is shown in . Figures 8a-8b As can be seen from Table 1, the registration method proposed in the present invention is superior to the other two existing algorithms in terms of algorithm stability and accuracy.

[0172] Table 1 Experimental results (* indicates registration failure)

[0173]

[0174] It can be seen from the above embodiments that the beneficial effects of the present invention are:

[0175] An embodiment of the present invention provides a multi-source remote sensing image registration method based on anisotropic diffusion description, which uses the diffusion function in the nonlinear diffusion equation to describe feature points and their neighborhoods to obtain feature descriptors. It can reduce the radiation characteristic differences and noise interference in heterogeneous images, and enhance the independence and robustness of the feature descriptors.

[0176] In the description of the present invention, the terms "first" and "second" are used for descriptive purposes only and should not be understood to indicate or imply relative importance or implicitly specify the number of the technical features indicated. Therefore, a feature specified as "first" or "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, "plurality" means two or more, unless otherwise specifically defined.

[0177] In the description of this specification, the reference terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" mean that the specific features, structures, materials, or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in any appropriate manner in any one or more embodiments or examples. In addition, those skilled in the art can combine and combine different embodiments or examples described in this specification.

[0178] Although the present application is described herein in conjunction with various embodiments, in the process of implementing the claimed application, those skilled in the art can understand and implement other variations of the disclosed embodiments by reviewing the drawings, the disclosure, and the appended claims.

[0179] The above is a further detailed description of the present invention in conjunction with specific preferred embodiments, and the specific implementation of the present invention should not be considered to be limited to these descriptions. For those skilled in the art of the present invention, without departing from the concept of the present invention, several simple deductions or substitutions can be made, which should be considered to fall within the scope of protection of the present invention.

Claims

1. A multi-source remote sensing image registration method based on anisotropic diffusion description, characterized in that: include: Acquire the optical image and SAR image to be registered; After filtering the optical image and the SAR image respectively using a nonlinear diffusion equation, an anisotropic scale space ASS of the optical image is established based on the first gradient of the optical image. O , and establish an anisotropic scale space ASS based on the first gradient of the SAR image S ; In the anisotropic scale space ASS O In the optical image, a first Hessian matrix is established point by point to generate a first Harris scale space, and in the anisotropic scale space ASS S In, establishing a second Hessian matrix for the SAR image point by point to generate a second Harris scale space; Feature points are detected layer by layer in the first Harris scale space of the optical image and the second Harris scale space of the SAR image, and the main directions are assigned to the detected feature points, and then the feature vectors of the feature points are calculated; Determine the reference image and the image to be registered from the optical image and the SAR image, and for the feature point P in the reference image O , calculate the Euclidean distance between it and each feature point in the image to be registered, and record the feature point pairs corresponding to the Euclidean distance that meet the first preset condition into the coarse matching point pair set CM; After removing the erroneous points in the coarse matching point pair set CM, performing primary matching on the remaining feature points in the coarse matching point pair set CM to form a fine matching point pair set FM and determine the transformation parameter θ; Generate a high confidence matching point pair set C based on the fine matching point pair set FM sample and low confidence matching point pair set C total , and for the high confidence matching point pair set C sample and the low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final ; Based on the final matching point pair set FM final Estimate the transformation parameter θ between the optical image and the SAR image final , and according to the transformation parameter θ final Transform the image to be registered into the coordinate system of the reference image.

2. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 1, characterized in that: Before the step of filtering the optical image and the SAR image respectively using a nonlinear diffusion equation, the method further includes: The first gradient of each target pixel in the optical image is calculated using the Sobel operator according to the following formula: in, represent the magnitude and direction of the first gradient of each target pixel in the optical image, Respectively represent the first gradient value of each target pixel in the optical image in the horizontal direction and the vertical direction, is the intensity value image of the optical image after Gaussian smoothing, Respectively represent the templates of Sobel in the horizontal and vertical directions, β j is the scale of the optical image, j represents the scale layer in the scale space; The Adaptive ROEWA operator is used to calculate the local exponentially weighted average ratio of each target pixel in the SAR image in the horizontal and vertical directions. The first gradient of each target pixel in the SAR image is calculated according to the local exponentially weighted average ratio in the horizontal and vertical directions according to the following formula: in, Respectively represent the amplitude and direction of the first gradient of each target pixel in the SAR image, Respectively represent the local exponential weighted average ratio of each target pixel in the SAR image in the horizontal and vertical directions, Respectively represent the first gradient value of each target pixel in the SAR image in the horizontal and vertical directions, α i is the scale of the scale layer where the target point is located in the SAR image.

3. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 2, characterized in that: In the anisotropic scale space ASS O In the optical image, a first Hessian matrix is established point by point to generate a first Harris scale space, and in the anisotropic scale space ASS S The step of establishing a second Hessian matrix point by point for the SAR image to generate a second Harris scale space comprises: In the anisotropic scale space ASS O The second gradient of the SAR image after filtering by the nonlinear diffusion equation is calculated, and the first Hessian matrix is established point by point in the second gradient corresponding to each scale layer: in, They represent the second gradient values of each target pixel in the horizontal and vertical directions in the SAR image after filtering by the nonlinear diffusion equation, Denotes the variance as σ i Gaussian kernel, σ i =α i / β i ,× represents convolution operation; According to the first Hessian matrix M S (α i ) Generate the first Harris scale space: R S (α i )=det(M S (α i ))-d·tr(M S (α i )) 2 Among them, det represents the determinant of the calculation matrix, d represents the corner detection factor, tr represents the trace of the calculation matrix, R S (α i ) represents the scale sequence {a1,...,α n }Generated Harris scale space; In the anisotropic scale space ASS S The second gradient of the optical image after filtering by the nonlinear diffusion equation is calculated, and the second Hessian matrix is established point by point in the second gradient corresponding to each scale layer: in, They represent the second gradient values of each target pixel in the horizontal and vertical directions in the optical R image after filtering by the nonlinear diffusion equation, Denotes the variance as β i Gaussian kernel of According to the second Hessian matrix M O (β i ) Generate the second Harris scale space: R O (β i )=det(M O (β i ))-d·tr(M O (β i )) 2 Among them, R O (β i ) represents the scale sequence {β1,...,β n }Generated Harris scale space.

4. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 1, characterized in that: Calculate the feature vector of the feature point as follows: A circular neighborhood is established with the feature point as the center and the preset length as the radius, and a logarithmic coordinate system is established with the feature point as the pole; Divide the circular neighborhood into three parts along the radial direction, and divide the inner circle and the two outer rings into 8 equal parts along the chord direction, forming 24 sub-areas; After establishing a rectangular coordinate system with the feature point as the origin and the main direction of the feature point as the positive direction of the horizontal axis, the feature point in each sub-area is rotated to the corresponding rectangular coordinate system; Divide 180° into 8 equal parts, and use the mapped diffusion function as the weight to sum the main directions of the feature points in each sub-region to obtain a one-dimensional feature vector with a length of 192 corresponding to each feature point.

5. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 1, characterized in that: Determine the reference image and the image to be registered from the optical image and the SAR image, and for the feature point P in the reference image O , respectively calculating the Euclidean distance between it and each feature point in the image to be registered, and recording the feature point pairs corresponding to the Euclidean distances that meet the first preset condition into the coarse matching point pair set CM, including: One of the optical image and the SAR image is used as the reference image and the other as the image to be registered. A feature point P is selected from the feature point set of the reference image. O , and calculate the feature point P O The Euclidean distance between the feature vector of and the feature vector of each feature point in the image to be registered; When the minimum value in the Euclidean distance satisfies the first preset condition, the feature point pair corresponding to the minimum value is recorded in the coarse matching point pair set CM.

6. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 5, characterized in that: The first preset condition is: Wherein, ED1 and ED2 represent the minimum value and the second minimum value in the Euclidean distance respectively, and Thres is a preset threshold.

7. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 5, characterized in that: After removing the erroneous points in the coarse matching point pair set CM, performing primary matching on the remaining feature points in the coarse matching point pair set CM to form a fine matching point pair set FM and determining the transformation parameter θ, the steps include: Select any two feature point pairs [P k ,Q k ]、[P l ,Q l ], respectively feature point pairs [P k ,Q k ]The corresponding Euclidean distance and feature point pair [P l ,Q l ] The ratio of the Euclidean distances D kl ; After traversing all feature point pairs, kl Perform histogram statistics; Eliminate the feature point pairs corresponding to the minimum value in the histogram, and calculate the root mean square error of the Euclidean distance of the remaining feature point pairs; Detect whether the root mean square error meets the second preset condition; if not, return to select any two feature point pairs from the coarse matching point pair set CM [P k ,Q k ]、[P l ,Q l ] step; if so, the point set consisting of the remaining feature point pairs is used as the SC-CM point set pair; The cascaded sample consistency estimation algorithm is used to perform first-level matching on the SC-CM point set pairs to obtain the fine matching point set pairs FM and transformation parameters θ.

8. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 7, characterized in that: The steps of performing first-level matching on the SC-CM point set pairs using the cascaded sample consistency estimation algorithm to obtain the fine matching point set pairs FM and the transformation parameters θ include: Randomly select three feature point pairs from the SC-CM point set and calculate the current transformation parameters; transforming the feature points of the reference image in the coarse matching point pair set CM using the current transformation parameters, and calculating the error based on the feature points of the image to be registered in the coarse matching point pair set CM; Counting the number num of feature point pairs that make the error smaller than a preset threshold; Check whether the preset number of iterations has been reached; if not, return to the step of randomly selecting three feature point pairs from the SC-CM point set and calculating the current transformation parameters; if so, calculate the maximum number num of feature point pairs; The feature point pairs corresponding to the maximum value of the feature point pair number num are obtained to form a fine matching point set pair FM, and the current transformation parameter is used as the transformation parameter θ.

9. The multi-source remote sensing image registration method based on anisotropic diffusion description according to claim 8, characterized in that: And for the high confidence matching point pair set C sample and the low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final ; The steps include: The feature points P of the reference image in the fine matching point pair set FM are I According to the transformation parameter θ, it is mapped to the image to be registered, and the nearest feature point Q is detected in the sub-area of the image to be registered. I ; If so, then the feature point pair [P I ,Q I ] Recorded in the high confidence matching point set C sample and low confidence matching point pair set C total If not, then find the nearest feature point Q in the image to be registered J , and record the feature point pairs into the low confidence matching point pair set C total ; The cascade sample consistency estimation algorithm is used to estimate the high confidence matching point pair set C sample and low confidence matching point pair set C total Perform secondary matching to form the final matching point pair set FM final .

Images