Method for quickly matching cross-modal remote sensing data of environment along railway

Abstract

The invention discloses a quick matching method for cross-modal remote sensing data of an environment along a railway. The method comprises the following steps: S1, data acquisition; s2, generating an anisotropic weighted torque graph through phase consistency model diffusion; s3, performing phase orientation description matching to obtain initial matching point pairs; s4, estimating through an initial solution of the FSC to obtain an initial inner point pair set and an initial transformation matrix; s5, optimizing the initial transformation matrix by using a weighted average optimization algorithm to obtain an expanded optimal transformation matrix; s6, calculating a matching error according to the expanded optimal transformation matrix, and dynamically adjusting an error threshold value; and S7, checking convergence: carrying out iteration until convergence. According to the method, global optimization is considered, the convergence solution is prevented from falling into local optimum, the problems of high mismatching rate and incomplete interior point screening caused by spectral difference of different-source images are solved, and the registration precision is remarkably optimized.

Description

Technical Field

[0001] This invention relates to the field of cross-modal remote sensing image processing, and in particular to a method for rapid matching of cross-modal remote sensing data of railway environment. Background Technology

[0002] In railway defect detection and the elimination of safety hazards along railway lines, it is often necessary to combine multiple data sources for comprehensive analysis. For example, high-resolution satellite imagery can provide macroscopic information on the distribution of railways over a wide area, while UAV imagery can acquire more detailed information about railway tracks, roadbeds, and the surrounding environment. Furthermore, multi-source data such as ground monitoring data, laser point clouds, and thermal infrared imagery can also be used for railway facility condition monitoring and hazard identification. To achieve effective fusion of these multi-dimensional and multi-modal data, precise spatial matching and alignment are first required, enabling pixel-level or feature-level comparative analysis of data from different sources within the same coordinate system.

[0003] Matching UAVs and satellite imagery in railway environments offers a degree of localization due to linear structural features (such as parallel / intersecting straight lines of tracks and bridges). However, due to differences in scale and perspective, tracks in UAV imagery exhibit perspective curvature while appearing as regular straight lines in satellite imagery. Furthermore, facilities such as platforms and signal lights are prone to occlusion and distortion in UAV imagery. In addition, satellite imagery can only present macroscopic features, while UAV imagery can capture details such as sleepers and ballast, leading to scale mismatch issues in traditional matching methods. Coupled with dynamic interference such as seasonal changes and train occlusion, this further increases the complexity of matching and the difficulty in describing matching features.

[0004] In multimodal image matching within a railway environment, gross error removal is a crucial step in the registration process. Gross error removal methods can be broadly categorized into two types: traditional gross error removal methods and deep learning-based gross error removal methods. Traditional methods are mostly based on random sampling and consistency detection strategies, such as the random sampling consensus algorithm and its improved version, the fast sampling consensus algorithm. These methods are prone to getting trapped in local optima in high-noise and multimodal image conditions, leading to incomplete removal and affecting matching quality. In recent years, researchers have proposed some improvements to address the shortcomings of the random sampling consensus algorithm, such as the fast sampling consensus algorithm to accelerate the sampling process and MLESAC to enhance interior point selection accuracy. However, limitations still exist, particularly in approximating the global optimum.

[0005] Deep learning methods such as SuperGlue and LOFTR achieve more accurate same-point matching by globally modeling feature points through self-attention mechanisms. However, these methods require substantial computational resources and training data to achieve good results, and their generalization ability to multimodal imagery remains limited. Deep learning-based image feature matching methods, such as CNN feature fusion, demonstrate strong adaptability in gross error removal. Especially in satellite and UAV cross-modal images of railway perimeter environments, deep neural networks can extract more complex features, but these methods still heavily rely on large-scale computational resources.

[0006] In summary, although traditional gross error removal methods and deep learning methods each have their advantages and disadvantages, they still face challenges in matching multi-platform (satellite, UAV) and multi-modal (optical, infrared, depth map) remote sensing data in the environment surrounding railway lines. These challenges include the strong repetition of linear structures such as tracks, large differences in image perspective scale, difficulty in feature description due to large differences in nonlinear radiation between cross-modal data, severe interference from dynamic targets, and high requirements for matching computation efficiency. When optimizing these methods, there are still challenges such as poor global convergence and incomplete removal of mismatched points. Summary of the Invention

[0007] To address the problems existing in the prior art, this invention provides a fast matching method for cross-modal remote sensing data of railway environment. This method takes into account global optimization, which not only avoids the converged solution from getting trapped in local optima, but also solves the problems of high mismatch rate and incomplete interior point screening caused by spectral differences in heterogeneous images, thus significantly improving the registration accuracy.

[0008] Therefore, the present invention adopts the following technical solution:

[0009] A method for rapid matching of cross-modal remote sensing data along railway lines includes the following steps:

[0010] S1, Data Acquisition: Use existing historical optical satellite images or navigation map images with precise geolocation as reference images; use the acquired cross-modal remote sensing data of the railway line environment as the images to be matched;

[0011] S2, Phase Consistency Model Diffusion: Extract the structural features of the reference image and the image to be matched, perform nonlinear diffusion, calculate the phase consistency of the obtained nonlinearly diffused structural feature map to obtain a phase map; normalize according to the maximum and minimum moments of the phase map to generate an anisotropic weighted moment map.

[0012] S3, Phase-oriented description matching: Extend the phase feature direction of the anisotropic weighted moment graph, combine it with the log-polar coordinate description vector frame, perform iterative calculation to obtain the phase-oriented feature descriptor; perform feature description matching based on the phase-oriented feature descriptor to obtain the initial matching point pair;

[0013] S4, FSC Initial Solution Estimation: Randomly select the minimum number of identical point pairs from the initial matched point pairs, and perform initial estimation of the transformation matrix through FSC to obtain a coarse estimate of the transformation matrix. Perform interior point filtering to obtain an initial set of interior point pairs. and the initial transformation matrix ;

[0014] S5, Weighted average optimization: The initial set of interior pairs is obtained in S4. and the initial transformation matrix Based on this, the search space is expanded by introducing the Levy flight strategy, and the transformation matrix parameters are iteratively optimized by combining the weighted average position of the matching point pairs, further improving the accuracy of the transformation matrix and obtaining the expanded optimal transformation matrix. and the corresponding matching point pairs;

[0015] S6, the extended optimal transformation matrix obtained from S5. Calculate the matching error and dynamically adjust the error threshold to obtain more inlier matching pairs;

[0016] S7, Check convergence: Perform iteration, compare the residual value calculated in this iteration with the error threshold set before the iteration or determine whether the number of iterations has reached the maximum number of iterations limit, decide whether to stop the iteration, if not, return and execute S5-S7 until convergence.

[0017] In S4, for each pair of identical points... , satisfy Then, the least squares method is used to calculate the selected point pairs to obtain a rough estimate of the transformation matrix. :

[0018] , (3)

[0019] in, The minimum number of pairs of identical nodes;

[0020] The rough estimate of the transformation matrix It can be used to calculate the transformation error of matching point pairs.

[0021] The preferred option is the minimum number of identical pairs. =20.

[0022] The method for internal point filtering in S4 is as follows:

[0023] Calculate the transformation error for all matching point pairs. ,in yes After a rough estimate of the transformation matrix Calculated;

[0024] Set error threshold ,Will The points are taken as interior points to form an initial set of interior point pairs. ;

[0025] For the initial set of interior point pairs The initial transformation matrix is ​​obtained by performing a final least squares calculation on all interior points within the matrix. .

[0026] The preferred option is an error threshold. It is 3 pixels.

[0027] S5 includes the following steps:

[0028] S51, Calculate the weighted average position: Calculate the matching point pairs in the current population. Weighted average position Used to update the estimated transformation matrix, weighted average position. The calculation formula is:

[0029] , (4)

[0030] in:

[0031] For the first The weight of a matching pair is typically inversely proportional to its fitness value. During the initial calculation, The default value is 1;

[0032] The transformation parameters corresponding to this pair of matching points;

[0033] For matching pairs in the current population The number of interior points is determined during the first computation of the algorithm, using the initial set of interior point pairs obtained in S4. As a matching point pair in the current population ;

[0034] S52 updates the current transformation matrix using the Levy flight strategy. Matching points in the current population The extended optimal transformation matrix is ​​obtained. For the extended optimal transformation matrix Update accordingly. Recalculate and update the expanded optimal transformation matrix. The matching point pairs are located.

[0035] S52 includes the following steps:

[0036] S521 introduces the Levy flight strategy, resulting in the extended optimal transformation matrix. :

[0037] To expand the search space, the weighted average optimization algorithm introduces a Levy flight strategy, which generates the step size through the Levy distribution, enabling effective jumps globally and accelerating the search process. The step size is generated as follows:

[0038] , (5)

[0039] in, The initial transformation matrix obtained from S4 is used as the current transformation matrix during the first operation of the algorithm. As the current transformation matrix ; This is the adjustment coefficient for the Levy distribution; The random step size is generated by the Levy distribution; These are candidate solutions generated using the Levy flight strategy;

[0040] Candidate solutions With the current transformation matrix Compare, if candidate solutions If the fit is better, then replace the current transformation matrix. Meanwhile, using candidate solutions The corresponding set of intermediate interior points Replace matching pairs in the current population This serves as the starting point for subsequent iterations, ultimately yielding the expanded optimal transformation matrix. ;

[0041] S522, Update the extended optimal transformation matrix Matching point pair positions:

[0042] In the initialization phase of the weighted average optimization algorithm, the extended optimal transformation matrix obtained in S521 is... As the initial input, the expanded optimal transformation matrix Parameterized into a one-dimensional vector And generate a random population; optimize the objective function. As shown in the following formula:

[0043] (6)

[0044] in, It is the objective function of the weighted average optimization algorithm. This is the set of interior points currently being matched; For a one-dimensional vector A defined transformation matrix; The actual transformation matrix;

[0045] Weighted average optimization algorithm based on To optimize the calculation of the transformation parameters corresponding to the matching point pairs in the current matching interior point set. .

[0046] In each iteration, the transformation parameters corresponding to the matching point pairs in the population are... By weighted average position Update; the expanded optimal transformation matrix. The update formula is:

[0047] , (7)

[0048] in, The learning rate;

[0049] According to the updated Recalculate and update the extended optimal transformation matrix. The matching point pairs are located.

[0050] The preferred adjustment factor is the Levy distribution. The value range is 1.5 to 2.0.

[0051] The specific steps for S6 are as follows:

[0052] S61, Calculate the residual of the current matched point pair. The median or mean of the error is used as a representative statistic, as shown in the following formula:

[0053] (8)

[0054] in, Indicates the first The actual observed locations of each matching point pair; Indicates the first The optimized estimated positions of each matching point pair are the extended optimal transformation matrix obtained in S522. The positions of the matching points; It is the Euclidean norm;

[0055] S62, dynamically adjust the threshold based on the statistical characteristics of the residuals. To ensure that the size of the interior point set gradually increases, the dynamically adjusted threshold... The error threshold used to determine whether the current iteration has ended is calculated from the previous iteration; during the first iteration... For manual setting, threshold The adjustment method is as follows:

[0056] (9)

[0057] in, For adjustment coefficients; This is the median of the residuals for all current matching point pairs, reflecting the overall error level of the matching point pairs under the current transformation matrix; As an error threshold for determining whether the next iteration has ended;

[0058] After each iteration, the residual from the previous iteration is used as a weight for adjustment. By gradually attenuating the influence of outliers, the optimization process becomes more stable and converges faster. The weight adjustment calculation method is as follows:

[0059] (10)

[0060] in, This is the residual weight decay factor; As the weights for calculating the weighted average position in the next iteration.

[0061] The preferred one is the first iteration. It is 5 pixels; ; The value ranges from 0.1 to 0.5.

[0062] The fast matching method for cross-modal remote sensing data of railway environment of the present invention calculates the weighted average position of each matching point pair in the population by combining the advantages of FSC initial solution estimation and the global optimization characteristics of weighted average optimization algorithm, and performs global optimization by combining Levy flight strategy to improve the estimation accuracy of transformation matrix. The method of the present invention is used to iteratively adjust the extended optimal transformation matrix, dynamically change the error threshold and increase the number of correct corresponding point pairs.

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

[0064] 1. The method of this invention proposes a fusion framework based on weighted average optimization (WAA) and Levy flight strategy, which innovatively combines the fast interior point selection of the FSC algorithm with the global parameter optimization of WAA: by dynamically adjusting the transformation matrix parameters through weighted averaging and expanding the search space using the random step size generated by the Levy distribution, global optimization is taken into account, avoiding the convergence solution from getting stuck in local optima and ignoring the global transformation. This solves the convergence stagnation problem caused by local sampling bias in traditional methods (such as RANSAC and FSC), significantly improves the global convergence speed, and reduces the registration error in complex scenes such as depth images to less than 50% of that of traditional methods.

[0065] 2. This invention solves the problems of high mismatch rate and incomplete interior point screening caused by spectral differences in heterogeneous images by dynamically adjusting the error threshold through iterative calculation of the median of the residuals. It increases the number of correct matching points and reduces RMSE to about 50% of the traditional method, significantly optimizing registration accuracy and robustness.

[0066] 3. The method of this invention introduces a residual weight decay factor to accelerate the iterative optimization process, which solves the problem of easily getting stuck in long-term loops in traditional parameter iterative optimization. It achieves a high registration success rate in data such as night light images, and the number of matching points is significantly increased, providing an efficient and reliable solution for multimodal image registration. Attached Figure Description

[0067] Figure 1 A flowchart of one embodiment of the fast matching method of the present invention;

[0068] Figure 2 This is a schematic diagram of the phase orientation feature descriptor construction in an embodiment of the present invention;

[0069] Figure 3 The multimodal remote sensing image test dataset in this embodiment of the invention includes: (a) multi-temporal images, (b) depth images and optical images, (c) navigation maps and optical images, and (d) night light images and optical images.

[0070] Figure 4 This is a mosaic image comparison diagram of the multimodal remote sensing image matching results in an embodiment of the present invention. Detailed Implementation

[0071] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0072] See Figure 1 The method for rapid matching of cross-modal remote sensing data of railway environment according to the present invention includes the following steps:

[0073] S1, Data Acquisition: Use existing historical optical satellite images or navigation map images with precise geolocation as reference images; use cross-modal remote sensing data of the railway line environment collected by optical sensors or depth sensors carried on UAVs or satellites as images to be matched.

[0074] S2, Phase Consistency Model Diffusion, includes the following steps:

[0075] S21. Structural features of the reference image and the image to be matched are extracted through anisotropic filtering, and nonlinear diffusion is performed to obtain a nonlinearly diffused structural feature map. Nonlinear diffusion can better preserve image edge information, facilitate the extraction of image structural features, and increase the richness of feature points to a certain extent.

[0076] S22. Phase consistency calculation is performed on the structural feature map obtained after nonlinear diffusion in S21 to obtain a phase map. Then, an anisotropic weighted moment map is generated by normalizing the maximum and minimum moments of the phase map, which is used for subsequent feature point extraction and construction of phase-oriented feature descriptors.

[0077] Phase consistency calculation is a frequency domain feature detection method that can effectively detect features such as corners and edges at the phase superposition of harmonic components in the Fourier transform of an image. It is insensitive to changes in signal amplitude and therefore exhibits strong illumination and contrast invariance. In this invention, the phase consistency model is extended to two-dimensional space, as shown in the extended formula below:

[0078] , (1)

[0079] in, The phase feature results calculated for the phase consistency model are used as image gradient magnitude in this invention; Weighting factors for frequency expansion; pixels on the image Amplitude at scale s and direction o in the Log-Gabor filter; To be a minimal constant that avoids the denominator being zero, it is generally taken as... =1×10 -10 ; This is the noise threshold, which is dynamically calculated based on the minimum scale response of the image and the filter characteristics to adapt to different noise levels. Its value range is generally 0.5~1.0. To prevent the value from being negative; It is a two-dimensional phase deviation function.

[0080] S3, Phase Orientation Description Matching, includes the following steps:

[0081] S31, the phase feature direction of the anisotropic weighted moment graph obtained in S22 is extended, and iteratively calculated using a log-polar coordinate descriptor vector frame to obtain the phase orientation feature descriptor. In one embodiment of the present invention, the phase feature direction of the anisotropic weighted moment graph is extended using an odd-symmetric filter of the Log-Gabor function.

[0082] Phase consistency calculation in S2 can only be used to calculate image feature values. Relying solely on the amplitude of phase features cannot achieve robust matching. Feature direction, on the other hand, reflects the direction in which image features change most significantly, which is crucial for the construction of phase orientation feature descriptors.

[0083] Odd-symmetric filters of the Log-Gabor function The convolution result represents the energy variation of the image in a specific direction. Therefore, this filter is used to construct phase-consistent directional features, as follows:

[0084] Since the calculation requires projecting the filtering results from multiple directions onto the horizontal (X) and vertical (Y) directions to obtain the horizontal energy A and the vertical energy C, it is necessary to extend the directional characteristics of the odd-symmetric filter of the Log-Gabor function. The extended formula is shown below:

[0085] , (2)

[0086] in, The directional characteristics are consistent with the phase. This is the convolution result of an odd-symmetric filter in the direction θ.

[0087] However, since the convolution result is a convolution value in multiple directions, the phase direction feature values ​​have opposite directions, making it difficult to accurately describe the feature information of the image. Therefore, this invention expands the convolution result by taking the absolute value to reduce all the convolution result to the positive direction, and then obtains the absolute phase consistency direction feature through a non-negative constant term.

[0088] Figure 2This is a schematic diagram illustrating the construction of a phase-oriented feature descriptor in one embodiment of the present invention. The present invention selects a fixed circular neighborhood centered on the feature point and calculates its absolute phase consistency gradient and orientation features to generate a gradient histogram of the absolute phase consistency direction. This histogram is then evenly divided into 36 parts at 10° intervals. The absolute phase consistency gradient and orientation features of each part are statistically analyzed, and the peak direction of the histogram is selected as the principal direction of the feature point. After obtaining the principal direction of the feature point, the circular neighborhood is equally divided to form a grid of 41 sub-regions with approximately the same area. Finally, the 8-directional gradient and orientation features of the pixels in each sub-region are statistically analyzed to generate a phase-oriented feature descriptor containing a 328-dimensional feature vector.

[0089] S32, perform feature description matching based on phase orientation feature descriptors to obtain initial matching point pairs.

[0090] S4, FSC Initial Solution Estimation: The initial set of interior point pairs is obtained through the initial solution estimation of FSC (Fast Sample Consensus). and the initial transformation matrix The specific steps are as follows:

[0091] S41, randomly select the minimum number of identical point pairs from the initial matched point pairs obtained in S32, and perform an initial estimation of the transformation matrix to obtain a rough estimate of the transformation matrix. .

[0092] For each pair of identical points , satisfy Then, the least squares method is used to calculate the selected point pairs to obtain a rough estimate of the transformation matrix. :

[0093] , (3)

[0094] in, This represents the minimum number of pairs of identical nodes.

[0095] In one embodiment of the present invention, =20.

[0096] The coarsely estimated transformation matrix H can be used to calculate the transformation error of the matching point pair.

[0097] S42, perform interior point filtering to obtain the initial set of interior point pairs. :

[0098] Calculate the transformation error for all matching point pairs. , yes After a rough estimate of the transformation matrix Calculated.

[0099] Set error threshold ,Will The points are taken as interior points to form an initial set of interior point pairs. ;

[0100] For the initial set of interior point pairs The initial transformation matrix is ​​obtained by performing a final least squares calculation on all interior points within the matrix. .

[0101] In one embodiment of the present invention, It is 3 pixels.

[0102] S5, Weighted average optimization: Apply the weighted average optimization algorithm to the initial transformation matrix obtained in S4. Optimization was performed to further improve the accuracy of the transformation matrix.

[0103] The core objective of the weighted average algorithm (WAA) is to find the position of each matched pair in the population by weighting the average position of the matched pairs. By combining the Levy flight strategy, global optimization is performed to improve the estimation accuracy of the transformation matrix.

[0104] The specific steps for weighted average optimization are as follows:

[0105] S51, Calculate the weighted average position: Calculate the matching point pairs in the current population. Weighted average position Used to update the estimated transformation matrix, weighted average position. The calculation formula is:

[0106] , (4)

[0107] in:

[0108] For the first The weight of a matching pair is typically inversely proportional to its fitness value. During the initial calculation, The default value is 1;

[0109] The transformation parameters corresponding to this pair of matching points;

[0110] For matching pairs in the current population The number of interior points is determined during the first computation of the algorithm, using the initial set of interior point pairs obtained in S4. As a matching point pair in the current population .

[0111] S52 updates the current transformation matrix using the Levy flight strategy. Matching points in the current population It includes the following steps:

[0112] S521 introduces the Levy flight strategy, resulting in the extended optimal transformation matrix. :

[0113] To expand the search space, the weighted average optimization algorithm introduces a Levy flight strategy. This strategy generates the step size using the Levy distribution, enabling efficient jumps across the global scope and thus accelerating the search process. The step size is generated as follows:

[0114] , (5)

[0115] in, The initial transformation matrix obtained from S4 is used as the current transformation matrix during the first operation of the algorithm. As the current transformation matrix ; This is the adjustment coefficient for the Levy distribution, typically ranging from 1.5 to 2.0. is the random step size generated by the Levy distribution. These are candidate solutions generated through the Levy flight strategy, and their main function is to expand the search space.

[0116] Candidate solutions With the current transformation matrix Compare, if candidate solutions If the fit is better, then replace the current transformation matrix. Meanwhile, using candidate solutions The corresponding set of intermediate interior points Replace matching pairs in the current population This serves as the starting point for subsequent iterations, ultimately yielding the expanded optimal transformation matrix. .

[0117] The Levy flight strategy improves the optimization algorithm's ability to escape local optima by making jumps in the solution space.

[0118] S522, Update the extended optimal transformation matrix Matching point pair positions:

[0119] In the initialization phase of the weighted average optimization algorithm, the extended optimal transformation matrix obtained in S521 is... As the initial input, the expanded optimal transformation matrix Parameterized into a one-dimensional vector A random population is generated. The objective function is defined as minimizing the sum of matching errors for all interior points, as shown in the following equation:

[0120] (6)

[0121] in, It is the objective function of the weighted average optimization algorithm, where L is the set of currently matched interior points; For a one-dimensional vector A defined transformation matrix; This is the true transformation matrix.

[0122] Weighted average optimization algorithm based on To optimize the calculation of the transformation parameters corresponding to the matching point pairs in the current matching interior point set. .

[0123] In each iteration, the transformation parameters corresponding to the matching point pairs in the population are... By weighted average position Update the expanded optimal transformation matrix. The update formula is:

[0124] , (7)

[0125] in, This is the learning rate.

[0126] According to the updated Recalculate and update the extended optimal transformation matrix. The matching point pairs are located.

[0127] S6, Calculate the residual dynamic adjustment error threshold: To further improve matching accuracy, based on the extended optimal transformation matrix obtained in S522... The matching point pairs are located, the residuals are calculated, and the error threshold is dynamically adjusted to obtain more inlier point matching pairs. The specific steps are as follows:

[0128] S61, Calculate the residual of the current matched point pair. The median or mean of the error is used as a representative statistic, as shown in the following formula:

[0129] , (8)

[0130] in, Indicates the first The actual observed locations of each matching point pair; Indicates the first The optimized estimated position of each matching point pair is the extended optimal transformation matrix obtained in S522. The positions of the matching points; It is the Euclidean norm.

[0131] S62, dynamically adjust the threshold based on the statistical characteristics of the residuals. This is to ensure that the size of the interior point set gradually increases.

[0132] The error threshold used to determine whether the current iteration has ended is calculated from the previous iteration; the error threshold for the first iteration is... The value is manually set; in one embodiment of the invention, it is set to 5 pixels.

[0133] Dynamically adjust threshold The adjustment method is as follows:

[0134] (9)

[0135] in, To adjust the coefficient, ; This is the median of the residuals for all current matching point pairs, reflecting the overall error level of the matching point pairs under the current transformation matrix; This serves as an error threshold for determining whether the next iteration has ended.

[0136] After each iteration, the residual from the previous iteration is used as a weight for adjustment. By gradually attenuating the influence of outliers, the optimization process becomes more stable and converges faster. The weight adjustment calculation method is as follows:

[0137] (10)

[0138] in, This is the residual weight decay factor, typically ranging from 0.1 to 0.5; As the weights for calculating the weighted average position in the next iteration.

[0139] S7, Convergence Check: In each iteration, the transformation matrix estimated by the global optimization is used to re-select matching point pairs, and the newly obtained matching points are added to the existing matching point set. The convergence is checked by comparing the residual value calculated in this iteration with the error threshold set before the iteration (i.e., ...). Alternatively, determine if the maximum number of iterations has been reached, and decide whether to stop the iteration. If not, return and execute S5-S7 until convergence. After convergence, the final transformation matrix is ​​obtained. and the final set of pairs of points with the same name This is used for subsequent evaluation of the registration effect of heterogeneous images in a railway environment.

[0140] To compare the matching method of this invention with existing methods, the registration effect of heterogeneous images in a railway environment is evaluated using the final set of corresponding point pairs. Four sets of multimodal remote sensing images with significant registration difficulties are selected as test datasets. The test datasets are as follows: Figure 3 As shown, rotational differences, scale differences, and nonlinear radiometric differences are covered. For each image pair, quantitative verification is performed using the root-mean-square error (RMSE) and the number of matched corresponding points (MCN), where RMSE is measured in pixels.

[0141] Based on the transformation relationship obtained by puncture points on the images of the test dataset as the true value and the error calculated, the fast matching method of cross-modal remote sensing data of railway environment of the present invention is compared with the Random Sample Consensus Algorithm (RANSAC) and the Fast Sample Consensus Algorithm (FSC). The comparison results are shown in Table 1.

[0142] Table 1 Comparison of Outlier Removal Methods

[0143]

[0144] In the matching scenarios of 12 groups of multimodal remote sensing images with rotation, scale and large nonlinear radiometric differences in the table, after removing erroneous matching points using the method of the present invention, the average number of matching point pairs for image registration is 158 pairs; while the results of FSC and RANSAC algorithms are only 140 pairs and 141 pairs, respectively, and the number of matching point pairs of the present invention is improved by 11%.

[0145] The average RMSE of the method of this invention is 2.68 pixels, while the RMSEs of the FSC and RANSAC algorithms are 5.57 and 4.91 pixels, respectively. The RMSE of the method of this invention is reduced by about 50% compared with the other two methods.

[0146] The results above demonstrate that, compared to the FSC and RANSAC algorithms, the method of this invention can estimate the global matching transformation, avoiding local optima in the matching result optimization that could lead to widespread misalignment when performing affine relation transformations on the overall image. Furthermore, while considering global optima, the method of this invention can also filter out more correct matching point pairs, exhibiting a significant advantage.

[0147] Comparison of the method of this invention with mosaic images optimized by FSC and RANSAC algorithms, for example Figure 4As shown in the figure, based on the transformation matrices calculated by different methods, the image to be matched is projected onto the reference image. The two images are displayed alternately at a certain distance, and the projection difference between the projected image and the reference image can be clearly seen through the alternating lines. It can be seen from the figure that the mosaicking effect of the image corrected by calculating the transformation matrix based on the matching point pairs optimized by the method of this invention is significantly better than that of the FSC and RANSAC methods, further verifying that the optimization focus of the method of this invention is meaningful for achieving global optimum.

Claims

1. A method for rapid matching of cross-modal remote sensing data of railway environment, characterized in that, Includes the following steps: S1, Data Acquisition: Use existing historical optical satellite images or navigation map images with precise geolocation as reference images; use the acquired cross-modal remote sensing data of the railway line environment as the images to be matched; S2, Phase Consistency Model Diffusion: Extract the structural features of the reference image and the image to be matched, perform nonlinear diffusion, calculate the phase consistency of the obtained nonlinearly diffused structural feature map to obtain a phase map; normalize according to the maximum and minimum moments of the phase map to generate an anisotropic weighted moment map. S3, Phase-oriented description matching: Extend the phase feature direction of the anisotropic weighted moment graph, and perform iterative calculations using the logarithmic polar coordinate description vector framework to obtain the phase-oriented feature descriptor; perform feature description matching based on the phase-oriented feature descriptor to obtain the initial matching point pair; S4, FSC Initial Solution Estimation: Randomly select the minimum number of identical point pairs from the initial matched point pairs, and perform initial estimation of the transformation matrix through FSC to obtain a coarse estimate of the transformation matrix. ; Perform interior point filtering to obtain an initial set of interior point pairs. and the initial transformation matrix ; S5, Weighted average optimization: The initial set of interior pairs is obtained in S4. and the initial transformation matrix Based on this, the search space is expanded by introducing the Levy flight strategy, and the transformation matrix parameters are iteratively optimized by combining the weighted average position of the matching point pairs, further improving the accuracy of the transformation matrix and obtaining the expanded optimal transformation matrix. and the corresponding matching point pairs; S6, the extended optimal transformation matrix obtained from S5. Calculate the matching error and dynamically adjust the error threshold to obtain more inlier matching pairs; S7, Check convergence: Perform iteration, compare the residual value calculated in this iteration with the error threshold set before the iteration or determine whether the number of iterations has reached the maximum number of iterations limit, decide whether to stop the iteration, if not, return and execute S5-S7 until convergence.

2. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 1, characterized in that, In S4, for each pair of corresponding points... , satisfy Then, the least squares method is used to calculate the selected point pairs to obtain a rough estimate of the transformation matrix. : ， （3） in, The minimum number of pairs of identical nodes; The rough estimate of the transformation matrix It can be used to calculate the transformation error of matching point pairs.

3. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 2, characterized in that: =20。 4. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 2, characterized in that, The method for internal point filtering in S4 is as follows: Calculate the transformation error for all matching point pairs. ,in yes After a rough estimate of the transformation matrix Calculated; Set error threshold ,Will The points are taken as interior points to form an initial set of interior point pairs. ; For the initial set of interior point pairs The initial transformation matrix is ​​obtained by performing a final least squares calculation on all interior points within the matrix. .

5. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 2, characterized in that: It is 3 pixels.

6. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 4, characterized in that, S5 includes the following steps: S51, Calculate the weighted average position: Calculate the matching point pairs in the current population. Weighted average position Used to update the estimated transformation matrix, weighted average position. The calculation formula is: ， （4） in: For the first The weight of a matching pair is typically inversely proportional to its fitness value; during the initial calculation, The default value is 1; The transformation parameters corresponding to this pair of matching points; For matching pairs in the current population The number of interior points is determined during the first computation of the algorithm, using the initial set of interior point pairs obtained in S4. As a matching point pair in the current population ; S52 updates the current transformation matrix using the Levy flight strategy. Matching points in the current population The extended optimal transformation matrix is ​​obtained. For the extended optimal transformation matrix Update accordingly. Recalculate and update the expanded optimal transformation matrix. The matching point pairs are located.

7. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 6, characterized in that, S52 includes the following steps: S521 introduces the Levy flight strategy, resulting in the extended optimal transformation matrix. : To expand the search space, the weighted average optimization algorithm introduces a Levy flight strategy, which generates the step size through the Levy distribution, enabling effective jumps globally and accelerating the search process. The step size is generated as follows: ， （5） in, The initial transformation matrix obtained from S4 is used as the current transformation matrix during the first operation of the algorithm. As the current transformation matrix ; This is the adjustment coefficient for the Levy distribution; The random step size is generated by the Levy distribution; These are candidate solutions generated using the Levy flight strategy; Candidate solutions With the current transformation matrix Compare, if candidate solutions If the fit is better, then replace the current transformation matrix. At the same time, it is used with candidate solutions The corresponding set of intermediate interior points Replace matching pairs in the current population This serves as the starting point for subsequent iterations, ultimately yielding the expanded optimal transformation matrix. ; S522, Update the extended optimal transformation matrix Matching point pair positions: In the initialization phase of the weighted average optimization algorithm, the extended optimal transformation matrix obtained in S521 is... As the initial input, the expanded optimal transformation matrix Parameterized into a one-dimensional vector And generate a random population; optimize the objective function. As shown in the following formula: （6） in, It is the objective function of the weighted average optimization algorithm. This is the set of interior points currently being matched; For a one-dimensional vector A defined transformation matrix; The actual transformation matrix; Weighted average optimization algorithm based on To optimize the calculation of the transformation parameters corresponding to the matching point pairs in the current matching interior point set. ; In each iteration, the transformation parameters corresponding to the matching point pairs in the population are... By weighted average position Update; the expanded optimal transformation matrix. The update formula is: ， （7） in, The learning rate; According to the updated Recalculate and update the extended optimal transformation matrix. The matching point pairs are located.

8. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 7, characterized in that: The value range is 1.5 to 2.

0.

9. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 7, characterized in that, The specific steps for S6 are as follows: S61, Calculate the residual of the current matched point pair. The median or mean of the error is used as a representative statistic, as shown in the following formula: （8） in, Indicates the first The actual observed locations of each matching point pair; Indicates the first The optimized estimated positions of each matching point pair are the extended optimal transformation matrix obtained in S522. The positions of the matching points; It is the Euclidean norm; S62, dynamically adjust the threshold based on the statistical characteristics of the residuals. To ensure that the size of the interior point set gradually increases, the dynamically adjusted threshold... The error threshold used to determine whether the current iteration has ended is calculated from the previous iteration; during the first iteration... For manual setting, threshold The adjustment method is as follows: （9） in, For adjustment coefficients; This is the median of the residuals for all current matching point pairs, reflecting the overall error level of the matching point pairs under the current transformation matrix; As an error threshold for determining whether the next iteration has ended; After each iteration, the residual from the previous iteration is used as a weight for adjustment. By gradually attenuating the influence of outliers, the optimization process becomes more stable and converges faster. The weight adjustment calculation method is as follows: （10） in, This is the residual weight decay factor; As the weights for calculating the weighted average position in the next iteration.

10. The method for rapid matching of cross-modal remote sensing data along railway lines according to claim 9, characterized in that: First iteration It is 5 pixels; ; The value ranges from 0.1 to 0.5.

Images