An optical remote sensing image registration method based on multi-source feature fusion

Abstract

The application discloses a multi-source feature fusion-based optical remote sensing image registration method, and relates to the field of remote sensing image processing, and comprises the following steps: acquiring band gray images and selecting a reference band gray image; performing structure feature enhancement on the target and reference band gray images to obtain first and second structure feature images; performing mask operation on the first and second structure feature images to obtain third and fourth structure feature images; generating a regular grid control point to obtain a corner point; generating a hybrid feature point set; determining a reliable displacement vector set by using forward and reverse LK pyramid optical flow tracking; performing outlier elimination to obtain an inlier displacement vector set; performing linear interpolation on the expanded inlier set to construct an initial dense displacement field; obtaining a dense continuous elastic coordinate transformation field; and performing inverse mapping resampling on the target band gray image to obtain a registration image. The application solves the problem of feature matching failure caused by radiation difference in a cross-spectral band scene.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and specifically to an optical remote sensing image registration method based on multi-source feature fusion. Background Technology

[0002] In quantitative remote sensing applications using multispectral and hyperspectral methods, the accuracy of geometric registration between bands affects the reliability of these applications. For multispectral images acquired through multiple asynchronous exposures, spatial misalignment occurs between images of different bands due to spatial deviations in the detector arrays corresponding to each band on the focal plane. This misalignment, known as inter-band registration error, is caused by factors such as satellite platform flutter and optical system distortion. This error leads to inconsistent pixel positions of the same ground object in different bands, directly impacting the computational accuracy of quantitative remote sensing products and the reliability of quantitative applications such as spectral mixing analysis and change detection.

[0003] For band registration, existing methods are mainly divided into three categories: feature matching-based, geometric parameter calibration-based, and deep learning-based. Feature matching-based methods establish geometric correspondences between bands by extracting local feature descriptions of the image. However, these descriptions are based on local gray-level distributions and are highly sensitive to differences in radiometric response between bands. In cross-spectral registration scenarios such as visible light and near-infrared, the success rate of feature matching drops significantly due to the significant differences in radiometric response of the same ground object in different bands, making it difficult to guarantee registration accuracy. Geometric parameter calibration-based methods comprehensively utilize satellite orbital parameters, precise attitude quaternions, camera in-camera orientation elements, and focal plane detector arrangement parameters to construct a rigorous geometric model and calculate the precise correspondence between detectors of each band in the image-side coordinate system, thereby establishing geometric correspondences between bands. However, the registration accuracy of this method is highly dependent on the accuracy and timeliness of the on-orbit geometric calibration parameters. It is not adaptable to time-varying random residual errors caused by satellite flutter, thermal deformation, and complex on-orbit environments, and it is difficult to guarantee band registration accuracy under non-steady-state imaging conditions. Deep learning-based methods learn the geometric deformation field directly from image pairs through convolutional neural networks or optical flow estimation networks. However, the performance of such methods is highly dependent on large-scale training datasets, the generalization ability of the model is insufficient, there is a lack of sufficient adaptive training sample sets in cross-spectral scenarios, and the interpretability of the model inference process is poor.

[0004] Furthermore, in the existing registration process, the practice of resampling the image multiple times independently in each step leads to the gradual accumulation of radiometric interpolation errors, which destroys the radiometric fidelity of the image and affects the spectral consistency of subsequent quantitative analysis. Summary of the Invention

[0005] This invention aims to at least solve the aforementioned technical problems existing in the prior art. To this end, this invention provides an optical remote sensing image registration method based on multi-source feature fusion, the method comprising: Acquire multiple band grayscale images of optical remote sensing images, and select a reference band grayscale image from the multiple band grayscale images; Structural feature enhancement is performed on the target band grayscale image and the reference band grayscale image respectively to obtain a first structural feature image and a second structural feature image; An effective pixel mask is obtained based on the target band grayscale image and the reference band grayscale image. Masking operations are then performed on the first structural feature image and the second structural feature image based on the effective pixel mask to obtain the third structural feature image and the fourth structural feature image. Based on the third structural feature image, regular grid control points are generated, and corner points are detected on the fourth structural feature image to obtain corner points; A set of hybrid feature points is generated based on the control points of the regular grid and the corner points; Based on the aforementioned set of hybrid feature points, and using forward LK pyramid optical flow tracing and reverse LK pyramid optical flow tracing, a set of reliable displacement vectors is determined. Outlier values ​​are removed from the reliable displacement vector set to obtain the interior point displacement vector set; Anchor points are added at the four corners of the pixel coordinate grid to expand the interior points, resulting in an expanded set of interior points. Delaunay triangulation is used to perform linear interpolation on the expanded set of interior points to construct an initial dense displacement field. The initial dense displacement field is subjected to median filtering and Gaussian smoothing to obtain a dense continuous elastic coordinate transformation field; Using the dense continuous elastic coordinate transformation field as a spatial index, the original image of the target band is reverse-mapped and resampled to obtain a registered image.

[0006] Optionally, the step of performing structural feature enhancement on the target band grayscale image and the reference band grayscale image respectively to obtain a first structural feature image and a second structural feature image includes: The target band grayscale image and the reference band grayscale image are subjected to contrast-limited adaptive histogram equalization and Gaussian filtering to obtain a preprocessed image; the preprocessed image includes a first preprocessed image and a second preprocessed image; the first preprocessed image is the preprocessed image corresponding to the target band grayscale image; the second preprocessed image is the preprocessed image corresponding to the reference band grayscale image. The horizontal and vertical gradient responses of the preprocessed image are obtained using the Scharr operator. The gradient magnitude map of the preprocessed image is obtained based on the horizontal gradient response and the vertical gradient response; The gradient magnitude map is normalized to obtain a structural feature image.

[0007] Optionally, obtaining an effective pixel mask based on the target band grayscale image and the reference band grayscale image includes: The target band grayscale image and the reference band grayscale image are marked based on preset dark pixel thresholds and saturated pixel thresholds to generate a first invalid pixel mask; A morphological dilation operation is performed on the first invalid pixel mask to obtain a second invalid pixel mask; the second invalid pixel mask is an invalid pixel mask with expanded boundaries; Invert the second invalid pixel mask to obtain the valid pixel mask.

[0008] Optionally, determining the reliable displacement vector set based on the hybrid feature point set and employing forward LK pyramid optical flow tracing and reverse LK pyramid optical flow tracing includes: Starting from the set of mixed feature points, using the second structural feature image as the source image and the first structural feature image as the target image, forward LK pyramid optical flow tracing is performed to obtain the forward tracing target coordinates and the corresponding first point-by-point valid flags. Starting from the forward tracking target coordinates, and using the first structural feature image as the source image and the second structural feature image as the target image, reverse LK pyramid optical flow tracing is performed to obtain the reverse recovery coordinates and the corresponding second point-by-point valid flags. The forward and backward cyclic errors are obtained based on the coordinates of the mixed feature points and the reverse recovery coordinates. The tracking displacement amplitude is obtained based on the forward tracking target coordinates and the mixed feature point coordinates; The mixed feature points are filtered based on the first point-by-point valid flag, the second point-by-point valid flag, the forward and backward cyclic error, and the tracking displacement amplitude to obtain a reliable displacement vector set.

[0009] Optionally, the step of removing outliers from the reliable displacement vector set to obtain the interior point displacement vector set includes: A K-nearest neighbor spatial index is constructed for the control points corresponding to the reliable displacement vector, and the K nearest neighbors of each control point are queried. Based on the reliable displacement vectors corresponding to the K nearest neighbors of each control point, the horizontal and vertical displacements of the K nearest neighbors of each control point are obtained. The corresponding local median is obtained based on the horizontal and vertical displacements of the K nearest neighbors of the control point; Based on the reliable displacement vector of the control point and the local median of the horizontal displacement and the local median of the vertical displacement, a first median absolute deviation and a second median absolute deviation are obtained; the first median absolute deviation is the median absolute deviation in the horizontal direction; the second median absolute deviation is the median absolute deviation in the vertical direction. The control points are filtered based on a preset local deviation tolerance threshold, the first median absolute deviation, and the second median absolute deviation to obtain a set of interior point displacement vectors.

[0010] Optionally, after obtaining the set of interior point displacement vectors, the method further includes: Based on the set of interior point displacement vectors, determine the global median of the horizontal displacement and the global median of the vertical displacement. The third median absolute deviation is obtained based on the global median of the horizontal displacement. The fourth median absolute deviation is obtained based on the global median of the vertical displacement. A first global threshold is set based on the third median absolute deviation; the first global threshold is the threshold value of the absolute value of the difference between the interior point displacement vector and the global median value of the horizontal displacement. A second global threshold is set based on the fourth median absolute deviation; the second global threshold is the threshold value of the absolute value of the difference between the in-point displacement vector and the global median value of the vertical displacement. The set of interior point displacement vectors is filtered based on the first global threshold and the second global threshold to obtain the final set of interior point displacement vectors.

[0011] Optionally, the step of adding anchor points at the four corner points of the pixel coordinate grid to expand the interior points and obtain an expanded set of interior points includes: Add an anchor point at each of the four corner points of the pixel coordinate grid as a new interior point; For each anchor point, obtain the nearest first interior points to the anchor point and their corresponding displacements; The weight of the displacement of the first interior point is determined based on the Euclidean distance from the first interior point to the anchor point; The weighted median is obtained based on the displacement of the first interior point and the corresponding weight, and the weighted median is used as the displacement of the corresponding anchor point.

[0012] Optionally, the step of performing median filtering and Gaussian smoothing on the initial dense displacement field to obtain a dense continuous elastic coordinate transformation field includes: The initial dense displacement field is subjected to median filtering and Gaussian smoothing to obtain the processed displacement field; The processed displacement field is superimposed on the pixel coordinates of the target band grayscale image to obtain the coordinates of the pixels of the target band grayscale image in the reference band grayscale image, thereby obtaining the dense continuous elastic coordinate transformation field.

[0013] Optionally, the step of performing inverse mapping resampling on the original image of the target band using the dense continuous elastic coordinate transformation field as a spatial index to obtain a registered image includes: The pixel coordinates of the original target band image that should be sampled at each coordinate position in the registered image are obtained based on the dense continuous elastic coordinate transformation field. The pixel coordinates of the original image of the target band are rounded and bilinear interpolation is performed to obtain the registered image.

[0014] Optionally, the step of performing corner detection on the fourth structural feature image to obtain corner points includes: The corner points are obtained by performing corner detection on the fourth structural feature image using the Shi-Tomasi corner response function.

[0015] The present invention has at least the following beneficial effects: 1. This invention adopts a hybrid feature point strategy that combines regular grid and corner detection, and combines dynamic threshold bidirectional LK pyramid optical flow tracing and outlier removal to balance the spatial uniformity of feature point coverage and positioning accuracy. It significantly improves the retention rate of reliable displacement vectors in wide displacement distribution scenarios. 2. This invention extracts local geometric features that are independent of radiation intensity by combining contrast-adaptive histogram equalization with Scharr gradient magnitude, fundamentally solving the problem of feature matching failure caused by radiation differences in cross-spectral scenarios using traditional methods. 3. This invention constructs a dense, continuous elastic coordinate transformation field covering the entire image by combining Delaunay triangulation linear interpolation with elastic regularization, which can accurately model local elastic deformation of the image; and performs only one inverse mapping resampling on the original image of each band, which minimizes the superposition of radiometric interpolation errors caused by multiple resampling, and ensures the radiometric fidelity and spectral consistency of the image. Attached Figure Description

[0016] Figure 1 This is a flowchart of an optical remote sensing image registration method based on multi-source feature fusion provided in an embodiment of the present invention; Figure 2 This is the image to be registered in the embodiments of the present invention; Figure 3 This is the registered image in an embodiment of the present invention; Figure 4 This refers to local details of the image to be registered in an embodiment of the present invention. Figure 1; Figure 5 This refers to local details of the image to be registered in an embodiment of the present invention. Figure 2 ; Figure 6 This refers to local details of the registered image in an embodiment of the present invention. Figure 1 ; Figure 7 This refers to local details of the registered image in an embodiment of the present invention. Figure 2 . Detailed Implementation

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

[0018] This invention provides an optical remote sensing image registration method based on multi-source feature fusion, such as... Figure 1 As shown, it includes the following steps: Step 101: Obtain multiple band grayscale images of the optical remote sensing image, and select a reference band grayscale image from the multiple band grayscale images.

[0019] Specifically, the system reads multi-band raw image data acquired by an optical remote sensing camera. Each band image is W×H pixels in size, and the data type retains the original bit depth.

[0020] Percentile stretching and normalization are performed on the data of each band to convert it into an 8-bit grayscale image, and one of the bands is designated as the global reference band.

[0021] For the raw data of each band, extract all valid non-zero pixel values ​​to form a set of valid pixels. .like If the number of valid pixels is less than the preset minimum threshold, output an 8-bit image with all zeros and mark the band as invalid; otherwise, calculate the 2nd percentile on the set of valid pixels. With the 98th percentile Linear stretching and truncation of the entire image data are performed using the following expression:

[0022] In the formula, for The processed data of the pixels at that location. for Data of pixels before processing. The 2nd percentile in the set of valid pixels. The 98th percentile in the set of valid pixels. () indicates data truncation.

[0023] Among them, the denominator term To prevent division by zero anomalies when the effective dynamic range is extremely narrow, the calculation results are converted into 8-bit unsigned integers and used as input grayscale images for subsequent structural feature extraction and optical flow tracing.

[0024] Among all N bands, the band with index R is designated as the global reference band, and all subsequent target bands are registered based on this band. It should be noted that the global reference band, i.e., the R band, is determined according to the specific implementation scenario in this embodiment of the invention. It can be pre-specified according to the actual application scenario, for example, using a specific spectral band as the global reference band; it can also be automatically selected based on image quality evaluation indicators, such as selecting the band with the highest signal-to-noise ratio or the highest effective pixel ratio as the global reference band; or it can be an intermediate band as the global reference band.

[0025] Step 102: Perform structural feature enhancement on the target band grayscale image and the reference band grayscale image respectively to obtain the first structural feature image and the second structural feature image.

[0026] In one possible implementation, the step of performing structural feature enhancement on the target band grayscale image and the reference band grayscale image respectively to obtain a first structural feature image and a second structural feature image includes: The target band grayscale image and the reference band grayscale image are subjected to contrast-limited adaptive histogram equalization and Gaussian filtering to obtain a preprocessed image; the preprocessed image includes a first preprocessed image and a second preprocessed image; the first preprocessed image is the preprocessed image corresponding to the target band grayscale image; the second preprocessed image is the preprocessed image corresponding to the reference band grayscale image. The horizontal and vertical gradient responses of the preprocessed image are obtained using the Scharr operator. The gradient magnitude map of the preprocessed image is obtained based on the horizontal gradient response and the vertical gradient response; The gradient magnitude map is normalized to obtain a structural feature image.

[0027] Specifically, for any input 8-bit grayscale image I, the following three steps are performed to extract the structural feature image. The grayscale image here includes grayscale images of each band and a reference band grayscale image.

[0028] The first step is to apply contrast-constrained adaptive histogram equalization. The image is divided into several non-overlapping local blocks, and a gray-level histogram is independently constructed for each block, using contrast-constrained parameters. After cropping the histogram and evenly distributing the excess portion across all gray levels, histogram equalization is performed, followed by bilinear interpolation to eliminate block boundary artifacts, resulting in a contrast-enhanced image. .

[0029] The second step is... Apply Gaussian kernels are used for smoothing and noise reduction to obtain a smooth image. This suppresses the interference of high-frequency noise on subsequent gradient calculations.

[0030] The third step is to use the Scharr operator pair Calculate the gradient response in the horizontal direction separately. Gradient response in the vertical direction Scharr operator level kernel With vertical core Defined as:

[0031] Calculate gradient magnitude image :

[0032] right Perform Min-Max normalization as follows: Interval:

[0033] In the formula, This is a gradient magnitude image. for The minimum gradient magnitude in the figure, for The maximum gradient magnitude in the figure, This is the normalized gradient magnitude map, i.e., the structural feature map.

[0034] The output is then converted to an 8-bit unsigned integer as the final structural feature image. Compared to the Sobel operator, the Scharr operator has better rotation invariance and frequency domain response uniformity. The combination of contrast-limited adaptive histogram equalization and Scharr gradient magnitude processing retains only local edge and texture geometry information that is independent of radiation intensity, so that subsequent optical flow tracing can maintain stable matching performance even when there are significant differences in radiation response between bands.

[0035] Step 103: Obtain an effective pixel mask based on the target band grayscale image and the reference band grayscale image, and perform masking operations on the first structural feature image and the second structural feature image according to the effective pixel mask to obtain a third structural feature image and a fourth structural feature image.

[0036] In one possible implementation, obtaining an effective pixel mask based on the target band grayscale image and the reference band grayscale image includes: The target band grayscale image and the reference band grayscale image are marked based on preset dark pixel thresholds and saturated pixel thresholds to generate a first invalid pixel mask; A morphological dilation operation is performed on the first invalid pixel mask to obtain a second invalid pixel mask; the second invalid pixel mask is an invalid pixel mask with expanded boundaries; Invert the second invalid pixel mask to obtain the valid pixel mask.

[0037] An effective pixel mask is generated based on the original grayscale image rather than the structural feature image to ensure that the judgment of overly dark and saturated pixels reflects the original radiation intensity information. The specific implementation is as follows: grayscale image of reference band Grayscale image of the target band Generate a joint invalid pixel mask according to the following expression. :

[0038] In the formula, The grayscale threshold for dark pixels. The grayscale threshold for saturated pixels. This is an indicator function.

[0039] The joint invalid pixel mask obtained here That is, the first invalid pixel mask.

[0040] Based on size All-unified structural element pairs The morphological dilation operation is performed as shown in the following expression:

[0041] In the formula, This is the invalid pixel mask after morphological dilation processing, i.e., the second invalid pixel mask. The size of the structural element for expansion. Let be the radius of the expansion, representing the... Round down to the nearest integer. In the first invalid pixel mask Pixel value at that location, Relative to The offset, Relative to The offset.

[0042] Invert the invalid mask after dilation to obtain the valid pixel mask. Only effective pixel mask The pixel position is used in subsequent feature point selection. Morphological dilation expands the boundary width of the invalid region outwards. Pixels are used to prevent feature points at the edges of the effective region from introducing incorrect matches due to invalid pixels in their optical flow tracing neighborhood.

[0043] Then, based on the effective pixel mask, masking operations are performed on the first structural feature image and the second structural feature image to obtain the third structural feature image and the fourth structural feature image.

[0044] Step 104: Generate regular grid control points based on the third structural feature image, and perform corner detection on the fourth structural feature image to obtain corner points.

[0045] In one possible implementation, the corner detection of the fourth structural feature image to obtain corner points includes: The corner points are obtained by performing corner detection on the fourth structural feature image using the Shi-Tomasi corner response function.

[0046] Specifically, with step size g in the effective internal region of the third structural feature image The following formula is used to generate rule grid control points:

[0047] In the formula, For the set of control points of a regular grid, The coordinates refer to the coordinates of the grid control points. Integer level, This indicates that the coordinates are integer multiples of g. Integer index, This indicates that the coordinates are integer multiples of g. Integer index.

[0048] Regular grids ensure that feature points are spatially uniformly covered within the effective area of ​​the image.

[0049] In effective pixel mask Within the constraints, corner detection is performed on the reference structural feature image using the Shi-Tomasi corner response function. Here, the structural feature image is the fourth structural feature image, obtained after masking. For each pixel position... Construct the gradient second-order moment matrix within its local neighborhood Ω centered at Ω. The expression is as follows:

[0050] In the formula, The coordinates are within the local neighborhood Ω. coordinates The horizontal gradient at point is, coordinates The vertical gradient at that point.

[0051] The corner response intensity is defined as the smallest eigenvalue of the second-order moment matrix of the gradient:

[0052] In the formula, The corner response intensity, This represents the smallest eigenvalue in the matrix.

[0053] Diagonal response intensity diagram Non-maximum suppression is performed, retaining response values ​​exceeding the maximum response value and quality factor of the entire graph. Find the extreme points of the product, and ensure that the distance between any two corner points is not less than g. Finally, output the set of corner points C, where the maximum number of corner points does not exceed the threshold. .

[0054] Regular grids ensure uniform spatial coverage but lack texture discrimination, while Shi-Tomasi corner point positioning has high accuracy but uneven spatial distribution. This invention combines the two to achieve both coverage uniformity and matching quality.

[0055] Step 105: Generate a set of mixed feature points based on the control points of the rule grid and the corner points.

[0056] Specifically, the set of grid control points G and the set of corner points C are merged, and duplicate coordinates are removed after rounding all coordinates. Then, out-of-bounds points exceeding the image boundaries are filtered out, along with their corresponding effective masks. Points with values ​​of zero will be output as the final set of blended feature points. .like If the number of feature points is less than the preset minimum feature point threshold, it degenerates into directly returning the identity coordinate transformation field.

[0057] Step 106: Based on the hybrid feature point set, and using forward LK pyramid optical flow tracing and reverse LK pyramid optical flow tracing, determine the reliable displacement vector set.

[0058] In one possible implementation, determining the reliable displacement vector set based on the hybrid feature point set and employing forward LK pyramid optical flow tracing and reverse LK pyramid optical flow tracing includes: Starting from the set of mixed feature points, using the second structural feature image as the source image and the first structural feature image as the target image, forward LK pyramid optical flow tracing is performed to obtain the forward tracing target coordinates and the corresponding first point-by-point valid flags. Starting from the forward tracking target coordinates, and using the first structural feature image as the source image and the second structural feature image as the target image, reverse LK pyramid optical flow tracing is performed to obtain the reverse recovery coordinates and the corresponding second point-by-point valid flags. The forward and backward cyclic errors are obtained based on the coordinates of the mixed feature points and the reverse recovery coordinates. The tracking displacement amplitude is obtained based on the forward tracking target coordinates and the mixed feature point coordinates; The mixed feature points are filtered based on the first point-by-point valid flag, the second point-by-point valid flag, the forward and backward cyclic error, and the tracking displacement amplitude to obtain a reliable displacement vector set.

[0059] Specifically, the LK optical flow algorithm is based on the assumption of brightness invariance and operates within a local window. Internal displacement Solve the following weighted least squares problem:

[0060] In the formula, It is a Gaussian weighted function. For local windows Pixel coordinates within, for Displacement in direction, for Displacement in direction, In the reference band grayscale image Pixels at that location, In the grayscale image of the target band The pixels at that location.

[0061] The pyramid multi-scale strategy recursively downsamples the image into an L-layer pyramid. After estimating the initial value of large displacement in the coarse-scale layer, it refines the image layer by layer to expand the effective displacement tracking range.

[0062] With a set of mixed feature points Starting from the reference structural feature image, i.e., the second structural feature image. The source image and the target structural feature image are the first structural feature image. For the target image, within the search window size Under the parameter setting of the number of pyramid levels L, forward LK pyramid optical flow tracing is performed to obtain the set of coordinates of the forward-traced target. and point-by-point valid mark .

[0063] The valid flag is a Boolean tracking success flag output by the LK pyramid optical flow algorithm for each input starting point: if a feature point fails to track due to reasons such as the search window going out of the image boundary, pyramid layer matching divergence, or excessive error, the point will be marked as invalid (invalid flag = 0); if the tracking converges and the result is reliable, it will be marked as valid (valid flag = 1).

[0064] by Starting from the target structural feature image, i.e., the first structural feature image. The source image and the reference structural feature image are the second structural feature image. For the target image, perform reverse LK pyramid optical flow tracing to obtain the reversed coordinate set. and point-by-point valid mark .

[0065] For each tracking point, calculate the forward and backward cycle errors. With tracking displacement amplitude :

[0066] In the formula, The coordinates of the mixed feature points, To track the target coordinates forward, This is for reverse coordinate recovery.

[0067] To track displacement amplitude Using as the independent variable, calculate the dynamic forward and backward consistency error threshold using the following formula. :

[0068] In the formula, This is the displacement proportionality coefficient. and These are the lower and upper bounds of the threshold, respectively. Tracking points that simultaneously satisfy the following four conditions are retained to form a reliable set of displacement vectors: ; ; ; ,in, This is the maximum allowed tracking displacement. If the number of valid points after filtering is less than the preset minimum, it degenerates into returning to the identity transformation field.

[0069] By adopting the dynamic threshold strategy described above, the tolerance for cyclic error is appropriately relaxed proportionally for tracking points with large displacement, while the upper limit of cyclic error is strictly controlled for tracking points with small displacement. Compared with a fixed threshold, this significantly improves the retention rate of reliable displacement vectors in scenarios with a wide range of displacement distribution.

[0070] Step 107: Remove outliers from the reliable displacement vector set to obtain the interior point displacement vector set.

[0071] In one possible implementation, the step of removing outliers from the reliable displacement vector set to obtain the interior point displacement vector set includes: A K-nearest neighbor spatial index is constructed for the control points corresponding to the reliable displacement vector, and the K nearest neighbors of each control point are queried. Based on the reliable displacement vectors corresponding to the K nearest neighbors of each control point, the horizontal and vertical displacements of the K nearest neighbors of each control point are obtained. The corresponding local median is obtained based on the horizontal and vertical displacements of the K nearest neighbors of the control point; Based on the reliable displacement vector of the control point and the local median of the horizontal displacement and the local median of the vertical displacement, a first median absolute deviation and a second median absolute deviation are obtained; the first median absolute deviation is the median absolute deviation in the horizontal direction; the second median absolute deviation is the median absolute deviation in the vertical direction. The control points are filtered based on a preset local deviation tolerance threshold, the first median absolute deviation, and the second median absolute deviation to obtain a set of interior point displacement vectors.

[0072] This invention employs a multi-level displacement vector outlier removal method, as detailed below: The first level involves K-nearest neighbor local outlier removal, constructing a K-dimensional tree spatial index for the set of effective control points, and for each control point... Query its K spatial nearest neighbor set Calculate the local median of horizontal and vertical displacements within the nearest neighbor set:

[0073]

[0074] In the formula, Control points The set of K spatial nearest neighbors The point in the middle, for Horizontal displacement, for The vertical displacement, Control points The local median of the horizontal displacement within the set of K spatial nearest neighbors. Control points The local median of the vertical displacement within the set of K spatial nearest neighbors.

[0075] Calculate the absolute deviation of the current control point displacement from the local median:

[0076] In the formula, Control points Horizontal displacement, Control points The vertical displacement, Control points The absolute deviation of the horizontal displacement from the local median, i.e., the absolute deviation of the first median. Control points The absolute deviation of the vertical displacement from the local median is called the second median absolute deviation.

[0077] Retain and satisfy simultaneously and The control points, among which This is the tolerance threshold for local deviations. If the total number of control points is less than K, then skip this stage of elimination.

[0078] The second level is the removal of outliers due to the absolute deviation of the global median.

[0079] In one possible implementation, after obtaining the set of interior point displacement vectors, the method further includes: Based on the set of interior point displacement vectors, determine the global median of the horizontal displacement and the global median of the vertical displacement. The third median absolute deviation is obtained based on the global median of the horizontal displacement. The fourth median absolute deviation is obtained based on the global median of the vertical displacement. A first global threshold is set based on the third median absolute deviation; the first global threshold is the threshold value of the absolute value of the difference between the interior point displacement vector and the global median value of the horizontal displacement. A second global threshold is set based on the fourth median absolute deviation; the second global threshold is the threshold value of the absolute value of the difference between the in-point displacement vector and the global median value of the vertical displacement. The set of interior point displacement vectors is filtered based on the first global threshold and the second global threshold to obtain the final set of interior point displacement vectors.

[0080] Specifically, for the set of interior points after the first-level screening, the global median of the horizontal and vertical displacements is calculated:

[0081] In the formula, This represents the global median of the horizontal displacement. The global median of the vertical displacement.

[0082] And the corresponding median absolute deviation:

[0083] In the formula, The third median absolute deviation, This represents the fourth median absolute deviation.

[0084] Set the global allowable range using a statistically adaptive approach:

[0085] In the formula, This is the amplification factor for the median absolute deviation. This is the minimum allowable margin. This represents the globally permissible range in the horizontal direction, i.e., the first global threshold. This represents the globally permissible range in the vertical direction, i.e., the second global threshold. This is a preset lower threshold for the global allowable range, used to prevent the erroneous removal of normal inliers when the local consistency of inliers is high and the statistical outlier is low, due to an excessively narrow global allowable range. Retaining inliers that meet the following conditions... and The interior point of .

[0086] It should be noted that interior points and control points represent the states of the same batch of spatial coordinate points at different stages in this embodiment of the invention. This embodiment of the invention filters interior points / control points using the above method, and the corresponding set of interior point displacement vectors can be obtained from the remaining interior points; that is, each point corresponds to a displacement vector.

[0087] The two-stage elimination process is executed in series, eliminating erroneous matching points from both the spatial neighborhood consistency and global statistical distribution levels, thereby significantly improving the input quality of subsequent dense displacement field interpolation.

[0088] Step 108: Add anchor points to the four corner points of the pixel coordinate grid to expand the interior points and obtain the expanded interior point set; use Delaunay triangulation to perform linear interpolation on the expanded interior point set to construct the initial dense displacement field.

[0089] In one possible implementation, the step of adding anchor points at the four corner points of the pixel coordinate grid to expand the interior points and obtain an expanded set of interior points includes: Add an anchor point at each of the four corner points of the pixel coordinate grid as a new interior point; For each anchor point, obtain the nearest first interior points to the anchor point and their corresponding displacements; The weight of the displacement of the first interior point is determined based on the Euclidean distance from the first interior point to the anchor point; The weighted median is obtained based on the displacement of the first interior point and the corresponding weight, and the weighted median is used as the displacement of the corresponding anchor point.

[0090] Specifically, to prevent extrapolation instability in sparse control point interpolation in the image boundary region, the coordinates of the four corner points of the pixel coordinate grid are... Add a virtual anchor point at each corner point or virtual anchor point, and for each corner point or virtual anchor point, search for the nearest point in space from the set of interior points. For each interior point, the weighted median of the horizontal and vertical displacements is calculated using the reciprocal of the Euclidean distance from each interior point to the corner point or virtual anchor point as the weight, and this median is used as the displacement assignment for that anchor point; if the number of valid interior points in the neighborhood of the corner point is insufficient... When this occurs, it degenerates into using the global median of all interior point displacements. and As anchor point displacement values, they constitute the expanded set of interior points. The set of interior points before expansion is the set of points corresponding to the set of interior point displacement vectors.

[0091] Perform Delaunay triangulation on the extended interior point set P+ to divide the image region into several non-overlapping triangles, based on the full-resolution pixel coordinate grid. Each pixel on Locate the triangle containing it, and let the three vertices of the triangle be (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 ... Calculate the centroid coordinates It satisfies the following expression:

[0092] Linear interpolation of the displacements of the three vertices using the centroid coordinates as weights:

[0093]

[0094] In the formula, After linear interpolation Horizontal displacement at that location, After linear interpolation Vertical displacement at that point As vertices Horizontal displacement, As vertices Horizontal displacement, As vertices Horizontal displacement, As vertices The vertical displacement, As vertices The vertical displacement, As vertices The vertical displacement.

[0095] For pixels located outside the convex hull of the inner point that still have empty values, the nearest neighbor interpolation method is used to fill in the missing values, ensuring that the displacement field is fully covered without holes.

[0096] The initial dense displacement field can be obtained through the above steps.

[0097] Step 109: Perform median filtering and Gaussian smoothing on the initial dense displacement field to obtain a dense continuous elastic coordinate transformation field.

[0098] In one possible implementation, the step of performing median filtering and Gaussian smoothing on the initial dense displacement field to obtain a dense continuous elastic coordinate transformation field includes: The initial dense displacement field is subjected to median filtering and Gaussian smoothing to obtain the processed displacement field; The processed displacement field is superimposed on the pixel coordinates of the target band grayscale image to obtain the coordinates of the pixels of the target band grayscale image in the reference band grayscale image, thereby obtaining the dense continuous elastic coordinate transformation field.

[0099] Specifically, for the horizontal component of the initial dense displacement field With vertical component Apply sizes of Median filtering:

[0100] In the formula, The horizontal value is the result of median filtering. for pixel coordinates in This is the horizontal offset. This is the vertical offset. For the initial dense displacement field The horizontal component at that point is also after linear interpolation. The horizontal displacement at that location.

[0101] The vertical value after median filtering According to and The same calculation method is used, and will not be elaborated further here. Median filtering has the best suppression effect on outliers in the L1 norm sense, effectively eliminating displacement spike noise introduced by a small number of residual mismatches, while preserving the local edge transition structure of the displacement field.

[0102] Apply a standard deviation to the median filtering result. Spatial smoothing using Gaussian kernels:

[0103]

[0104] In the formula, For Gaussian kernel function, This is the displacement field after Gaussian smoothing. This is the result of median filtering.

[0105] Gaussian smoothing approximates the discretization effect of applying Tikhonov regularization smoothing constraints to the displacement field, and is particularly effective when the interior points are densely distributed. The larger the value, the closer the displacement field approaches a global rigid transformation; The smaller the value, the better the details of local elastic deformation are preserved; it can be flexibly adjusted according to the degree of local distortion in the actual image. and Values.

[0106] Coordinate transformation field generation: The Gaussian-smoothed displacement field is superimposed onto the target band pixel coordinates to generate the following dense elastic coordinate transformation field:

[0107]

[0108] in, Represents the target band coordinates The corresponding coordinates of the pixel at that location in the reference band are used for subsequent inverse mapping and resampling. This represents the horizontal displacement in the displacement field after Gaussian smoothing. This represents the vertical displacement in the displacement field after Gaussian smoothing.

[0109] Step 110: Using the dense continuous elastic coordinate transformation field as the spatial index, perform inverse mapping resampling on the original image of the target band to obtain the registered image.

[0110] In one possible implementation, the step of performing inverse mapping resampling on the original target band image using the dense continuous elastic coordinate transformation field as a spatial index to obtain a registered image includes: The pixel coordinates of the original target band image that should be sampled at each coordinate position in the registered image are obtained based on the dense continuous elastic coordinate transformation field. The pixel coordinates of the original image of the target band are rounded and bilinear interpolation is performed to obtain the registered image.

[0111] Specifically, the field is transformed using dense elastic coordinates. For spatial indexing, inverse mapping resampling is performed on the original images for each band. For each pixel location in the registered image... Its corresponding reference coordinates Generally, the coordinates are non-integer coordinates, and their radiation intensity values ​​are calculated using bilinear interpolation. Let... , That is to , Round down to get , The decimal offsets in the horizontal and vertical directions are respectively and The bilinear interpolation formula is:

[0112] in, For the original image of the corresponding band, a grayscale image is generated without percentile stretching and normalization. The boundary region uses a boundary pixel copying and filling strategy to process pixels outside the effective imaging area. To register the image The value at that location, This is the decimal offset in the horizontal direction. This is the decimal offset in the vertical direction. For the original image of the corresponding band The value at that location, For the original image of the corresponding band The value at that location, For the original image of the corresponding band The value at that location, For the original image of the corresponding band The value at that point. The resampled output is rounded to the nearest integer and truncated to the range of valid original image data. This preserves the original image bit depth, meaning it does not change the data type, and outputs the final registered image. Among these... This represents the maximum effective data range of the original image.

[0113] The registration effect using the method of the embodiments of the present invention can be referred to Figure 2 , Figure 3 , Figure 4 , Figure 5 , Figure 6 , Figure 7 , Figure 2 The image to be registered, also known as the original image. Figure 3 For the registered image, i.e. Figure 2 The registered image of the original image; Figure 4 This is a local detail image of the image to be registered. Figure 5 This is another detailed view of the image to be registered. Figure 6 A detailed view of the registered image, i.e. Figure 4 Detailed view of the corresponding registration image. Figure 7 Another local detail image for the registration image, namely Figure 5 Local detail image of the corresponding registered image.

[0114] In summary, in this embodiment of the invention, a hybrid feature point strategy that combines regular grid and corner detection is adopted, along with dynamic threshold bidirectional LK pyramid optical flow tracing and outlier removal, which balances the spatial uniformity of feature point coverage and positioning accuracy, and significantly improves the retention rate of reliable displacement vectors in scenarios with a wide range of displacement distribution.

[0115] This invention proposes a cross-spectral structural feature enhancement scheme. By combining contrast-restricted adaptive histogram equalization with Scharr gradient magnitude, it extracts local geometric features that are independent of radiation intensity, fundamentally solving the problem of feature matching failure caused by radiation differences in cross-spectral scenarios using traditional methods. The embodiments of the present invention construct a dense and continuous elastic coordinate transformation field covering the entire image by combining Delaunay triangulation linear interpolation with elastic regularization, which can accurately model the local elastic deformation of the image; and perform only one inverse mapping resampling on the original image of each band, which minimizes the superposition of radiometric interpolation errors caused by multiple resampling, and ensures the radiometric fidelity and spectral consistency of the image.

[0116] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0117] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A method for optical remote sensing image registration based on multi-source feature fusion, characterized in that, include: Acquire multiple band grayscale images of optical remote sensing images, and select a reference band grayscale image from the multiple band grayscale images; Structural feature enhancement is performed on the target band grayscale image and the reference band grayscale image respectively to obtain a first structural feature image and a second structural feature image; The step of performing structural feature enhancement on the target band grayscale image and the reference band grayscale image respectively to obtain a first structural feature image and a second structural feature image includes: The target band grayscale image and the reference band grayscale image are subjected to contrast-limited adaptive histogram equalization and Gaussian filtering to obtain a preprocessed image; the preprocessed image includes a first preprocessed image and a second preprocessed image; the first preprocessed image is the preprocessed image corresponding to the target band grayscale image; the second preprocessed image is the preprocessed image corresponding to the reference band grayscale image. The horizontal and vertical gradient responses of the preprocessed image are obtained using the Scharr operator. The gradient magnitude map of the preprocessed image is obtained based on the horizontal gradient response and the vertical gradient response; The gradient magnitude map is normalized to obtain a structural feature image; An effective pixel mask is obtained based on the target band grayscale image and the reference band grayscale image. Masking operations are then performed on the first structural feature image and the second structural feature image based on the effective pixel mask to obtain the third structural feature image and the fourth structural feature image. Based on the third structural feature image, regular grid control points are generated, and corner points are detected on the fourth structural feature image to obtain corner points; A set of hybrid feature points is generated based on the control points of the regular grid and the corner points; Based on the aforementioned set of hybrid feature points, and using forward LK pyramid optical flow tracing and reverse LK pyramid optical flow tracing, a set of reliable displacement vectors is determined. Outlier values ​​are removed from the reliable displacement vector set to obtain the interior point displacement vector set; Anchor points are added at the four corners of the pixel coordinate grid to expand the interior points, resulting in an expanded set of interior points. Delaunay triangulation is used to perform linear interpolation on the expanded set of interior points to construct an initial dense displacement field. The initial dense displacement field is subjected to median filtering and Gaussian smoothing to obtain a dense continuous elastic coordinate transformation field; Using the dense continuous elastic coordinate transformation field as a spatial index, the original image of the target band is reverse-mapped and resampled to obtain a registered image.

2. The optical remote sensing image registration method based on multi-source feature fusion according to claim 1, characterized in that, Obtaining an effective pixel mask based on the target band grayscale image and the reference band grayscale image includes: The target band grayscale image and the reference band grayscale image are marked based on preset dark pixel thresholds and saturated pixel thresholds to generate a first invalid pixel mask; A morphological dilation operation is performed on the first invalid pixel mask to obtain a second invalid pixel mask; the second invalid pixel mask is an invalid pixel mask with expanded boundaries; Invert the second invalid pixel mask to obtain the valid pixel mask.

3. The optical remote sensing image registration method based on multi-source feature fusion according to claim 1, characterized in that, The determination of a reliable displacement vector set based on the hybrid feature point set, using forward LK pyramid optical flow tracing and reverse LK pyramid optical flow tracing, includes: Starting from the set of mixed feature points, using the second structural feature image as the source image and the first structural feature image as the target image, forward LK pyramid optical flow tracing is performed to obtain the forward tracing target coordinates and the corresponding first point-by-point valid flags. Starting from the forward tracking target coordinates, and using the first structural feature image as the source image and the second structural feature image as the target image, reverse LK pyramid optical flow tracing is performed to obtain the reverse recovery coordinates and the corresponding second point-by-point valid flags. The forward and backward cyclic errors are obtained based on the coordinates of the mixed feature points and the reverse recovery coordinates. The tracking displacement amplitude is obtained based on the forward tracking target coordinates and the mixed feature point coordinates; The mixed feature points are filtered based on the first point-by-point valid flag, the second point-by-point valid flag, the forward and backward cyclic error, and the tracking displacement amplitude to obtain a reliable displacement vector set.

4. The optical remote sensing image registration method based on multi-source feature fusion according to claim 1, characterized in that, The step of removing outliers from the reliable displacement vector set to obtain the interior point displacement vector set includes: A K-nearest neighbor spatial index is constructed for the control points corresponding to the reliable displacement vector, and the K nearest neighbors of each control point are queried. Based on the reliable displacement vectors corresponding to the K nearest neighbors of each control point, the horizontal and vertical displacements of the K nearest neighbors of each control point are obtained. The corresponding local median is obtained based on the horizontal and vertical displacements of the K nearest neighbors of the control point; Based on the reliable displacement vector of the control point and the local median of the horizontal displacement and the local median of the vertical displacement, a first median absolute deviation and a second median absolute deviation are obtained; the first median absolute deviation is the median absolute deviation in the horizontal direction; the second median absolute deviation is the median absolute deviation in the vertical direction. The control points are filtered based on a preset local deviation tolerance threshold, the first median absolute deviation, and the second median absolute deviation to obtain a set of interior point displacement vectors.

5. The optical remote sensing image registration method based on multi-source feature fusion according to claim 4, characterized in that, After obtaining the set of interior point displacement vectors, the process further includes: Based on the set of interior point displacement vectors, determine the global median of the horizontal displacement and the global median of the vertical displacement. The third median absolute deviation is obtained based on the global median of the horizontal displacement. The fourth median absolute deviation is obtained based on the global median of the vertical displacement. A first global threshold is set based on the third median absolute deviation; the first global threshold is the threshold value of the absolute value of the difference between the interior point displacement vector and the global median value of the horizontal displacement. A second global threshold is set based on the fourth median absolute deviation; the second global threshold is the threshold value of the absolute value of the difference between the in-point displacement vector and the global median value of the vertical displacement. The set of interior point displacement vectors is filtered based on the first global threshold and the second global threshold to obtain the final set of interior point displacement vectors.

6. The optical remote sensing image registration method based on multi-source feature fusion according to claim 1, characterized in that, The process of adding anchor points at the four corners of the pixel coordinate grid to expand the interior points, resulting in an expanded set of interior points, includes: Add an anchor point at each of the four corner points of the pixel coordinate grid as a new interior point; For each anchor point, obtain the nearest first interior points to the anchor point and their corresponding displacements; The weight of the displacement of the first interior point is determined based on the Euclidean distance from the first interior point to the anchor point; The weighted median is obtained based on the displacement of the first interior point and the corresponding weight, and the weighted median is used as the displacement of the corresponding anchor point.

7. The optical remote sensing image registration method based on multi-source feature fusion according to claim 1, characterized in that, The step of performing median filtering and Gaussian smoothing on the initial dense displacement field to obtain a dense continuous elastic coordinate transformation field includes: The initial dense displacement field is subjected to median filtering and Gaussian smoothing to obtain the processed displacement field; The processed displacement field is superimposed on the pixel coordinates of the target band grayscale image to obtain the coordinates of the pixels of the target band grayscale image in the reference band grayscale image, thereby obtaining the dense continuous elastic coordinate transformation field.

8. The optical remote sensing image registration method based on multi-source feature fusion according to claim 7, characterized in that, The step of performing inverse mapping resampling on the original image of the target band using the dense continuous elastic coordinate transformation field as a spatial index to obtain a registered image includes: The pixel coordinates of the original target band image that should be sampled at each coordinate position in the registered image are obtained based on the dense continuous elastic coordinate transformation field. The pixel coordinates of the original image of the target band are rounded and bilinear interpolation is performed to obtain the registered image.

9. The optical remote sensing image registration method based on multi-source feature fusion according to claim 1, characterized in that, The corner detection process for the fourth structural feature image, to obtain corner points, includes: The corner points are obtained by performing corner detection on the fourth structural feature image using the Shi-Tomasi corner response function.

Images