Tai-landform automatic identification method based on remote sensing and DEM
By constructing a gridded data cube of remote sensing DEM, implementing heterogeneous noise collaborative suppression and fractal beam mapping, the problems of multi-source data registration error and noise suppression non-coordination in Daigu landform identification were solved, achieving high-precision and automated Daigu landform identification.
Patent Information
- Application Number
- CN202511089577.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-08-05
- Publication Date
- 2025-11-14
AI Technical Summary
Existing technologies lack a unified data cube structure for Daigu topography identification, resulting in large registration errors of multi-source data and uncoordinated noise suppression, which reduces the accuracy of ridgeline positioning and processing efficiency.
A remote sensing DEM grid data cube is constructed, and the ridge type and peak edge of the Daigu landform are automatically extracted through heterogeneous noise collaborative suppression, sparse thinning, fractal beam energy level quantization, geometric ridge line and energy level high region depth mapping binding and self-consistent differential adaptive convergence.
It significantly improves the accuracy and efficiency of Daigu landform recognition, reduces human intervention, and ensures high-precision, physically consistent recognition results under different resolutions and noise backgrounds.
Smart Images

Figure CN120953808A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of image recognition technology, and particularly relates to terrain and landform recognition, specifically an automatic identification method for Daigu landform based on remote sensing and DEM. Background Technology
[0002] In recent years, with the increasing availability of high-resolution satellite imagery, synthetic aperture radar (SAR) imagery, and global digital elevation models (DEMs), researchers have been able to automatically identify and classify large geomorphic units using multi-source data. Existing technologies typically follow a basic framework of "spectral information-driven + morphological analysis supplementation": first, normalized vegetation index (NZVI), soil and rock index, or shadow distribution are extracted from optical remote sensing images; then, DEM-derived factors such as radar scattering coefficient, slope, and aspect are added, and geomorphic categories are determined using support vector machines, random forests, nearest neighbor clustering, or convolutional neural networks. Some scholars have proposed object-based image analysis, using multi-scale segmentation to divide images into blocks and then fusing elevation thresholds for a "valley-peak" binary classification; others have introduced wavelet decomposition, morphological gradients, and watershed algorithms to find ridge lines in the difference between dawn and dusk shadows on slopes. However, these methods are mostly aimed at macroscopic types such as hills, gullies, or valleys, and are still insufficient for isolated hilltop-steep-gentle slope composite systems like the Daigu landform, which is formed by block uplift and differential erosion. The main problems are as follows: First, existing methods lack a unified data cube structure. Optical images, radar images, and DEMs have varying original resolutions and coordinate references. A common practice is to resample them separately and then overlay them as raster grids. However, this often only aligns row and column indices in a two-dimensional plane, lacking strict constraints on a unified three-dimensional grid encompassing optical brightness, microwave scattering, and elevation dimensions. This loose registration leads to misalignment of corresponding pixels at steep ridge edges, causing slope gradient inversion shifts and reducing ridge line positioning accuracy.
[0003] Second, noise suppression remains at the level of single-dimensional or empirical thresholds. The optical brightness dimension is susceptible to atmospheric scattering and sensor radiation errors, the microwave scattering dimension suffers from multiplicative speckle noise, and the elevation dimension contains isolated spike errors. Existing denoising strategies often process the three-dimensional signals separately, ignoring the complementary relationship between the three dimensions in the real terrain response, and are unable to suppress high-frequency interference in a coordinated manner while preserving ridge details. Summary of the Invention
[0004] The main objective of this invention is to provide an automatic identification method for the Daigu landform based on remote sensing and DEM. This invention constructs a remote sensing DEM data cube with strictly identical optical brightness, microwave scattering, and elevation dimensions. Based on this cube, it sequentially implements heterogeneous noise collaborative suppression, sparsity thinning, fractal beam level quantization, depth mapping binding of geometric ridges and high-energy-level regions, odd-even interleaved spectral lateral expansion, and self-consistent differential adaptive convergence, thereby automatically extracting the ridge type and peak edge of the Daigu landform. This process eliminates multi-source data registration errors, significantly reduces the interference of random noise on texture determination, fuses multi-dimensional information with unified energy-level markers, and fully utilizes temporal brightness differences to enhance ridge boundaries. Ultimately, it can automatically output high-precision, physically consistent Daigu landform identification results under different resolutions and noise backgrounds, reducing manual intervention and improving processing efficiency.
[0005] To address the aforementioned technical problems, this invention provides an automatic identification method for Daigu landforms based on remote sensing and DEM, the method comprising: Step 1: Through multi-source remote sensing image acquisition, simultaneously acquire optical remote sensing images, microwave remote sensing images and publicly available digital elevation models (DEMs) covering the target area, and perform grid registration processing on the three types of data to construct a grid data cube of remote sensing DEM data with unified resolution, including three dimensions: optical brightness dimension, microwave scattering dimension and elevation value dimension. Step 2: Perform heterogeneous noise suppression processing on the data of each dimension of the remote sensing DEM grid data cube to obtain the preprocessed grid data cube; Step 3: Based on the preprocessed grid data cube, while maintaining the spatial consistency of the grid, a sparse candidate raster set is generated by sparse thinning. Fractal beam mapping is performed on the sparse candidate raster set to obtain the fractal beam energy level index. Then, the peak ridge line of the elevation profile is deeply mapped and bound to the high energy level region. The ridge feature is extended by combining optical remote sensing image interleaving. When the self-consistent differential change rate meets the termination criterion, the automatic identification result of Daigu landform is output.
[0006] Furthermore, in step 1, the process of performing grid registration on the three types of data includes: projecting the optical remote sensing image, microwave remote sensing image, and digital elevation model (DEM) onto the same geodetic coordinate reference frame and the same map projection plane to ensure that the three types of data have a unified longitude, latitude, and elevation benchmark; using the spatial resolution of the optical remote sensing image as the benchmark resolution; if the resolution of the microwave remote sensing image or the DEM is inconsistent with the benchmark resolution, adjusting its spatial resolution to be consistent with the optical remote sensing image through resampling; after unifying the resolution, calculating the raster origin coordinates of the three types of data; if there is an offset between the raster origin coordinates of the microwave remote sensing image and the optical remote sensing image, or between the DEM and the optical remote sensing image, the microwave remote sensing image and the DEM are translated by the whole raster pixel spacing to make their raster origin coordinates completely consistent with the raster origin coordinates of the optical remote sensing image, so that the raster origin coordinates of the three types of data completely coincide.
[0007] Furthermore, step 2, which involves implementing heterogeneous noise suppression for each dimension of the remote sensing DEM grid data cube, specifically includes: calculating the noise root mean square error to noise variance ratio for each dimension of the remote sensing DEM grid data in its respective frequency band; if the noise root mean square error of any dimension exceeds 1.2 times the overall mean, then the corresponding raster for that dimension is recorded as a high-noise raster; for the optical brightness dimension and microwave scattering dimension where the high-noise raster is located, an adaptive sliding window averaging method is used for smoothing; subsequently, morphological opening and closing operations are performed sequentially on the elevation value dimension at the same raster position to reduce isolated elevation spike noise; after completion, the noise root mean square error is recalculated, and if the reduction does not reach 30%, the sliding window averaging method is performed again with the same window size; Furthermore, in step 3, the process of generating a sparse candidate raster set includes: in the preprocessed grid data cube, first read the number of rows, columns, and layers; set the row thinning interval to 3 rasters, the column thinning interval to 3 rasters, and the layer thinning interval to 2 layers; if the number of rows or columns cannot be divided by the interval value, round up to maintain the consistency of the grid space; traverse each layer row by row and column by column in the order of 3 rasters in the row direction and 3 rasters in the column direction; when the row index and column index simultaneously meet the interval condition, mark the corresponding raster as a row-column sparse raster; after the traversal, obtain a row-column sparse raster set; based on the row-column sparse raster set, extract upwards at a layer interval of 2 layers; when the layer index meets the layer direction interval condition, retain the corresponding raster, otherwise discard it; after completion, generate a sparse candidate raster set, and record the row index, column index, and layer index for each raster.
[0008] Furthermore, after generating the sparse candidate raster set in step 3, the sparse candidate raster set will be processed as follows: In the sparse candidate raster set, at least 5% of the total number of raster pairs are randomly selected, and their eight-way neighbor raster integrity in the preprocessed same-grid data cube is compared; if the number of missing neighbor raster pairs of any raster pair exceeds 1, the overall row or column direction thinning interval is reduced by 1 raster and sparse thinning is re-executed until the neighborhood integrity passes the check; for each sparse candidate raster, a cube local fractal window with a side length of 5 rasters is built with the raster as the center; if the window touches the same-grid boundary, the window raster is supplemented in a mirror manner to maintain the integrity of the window shape.
[0009] Furthermore, in step 3, the process of performing fractal beam mapping on the sparse candidate grid set to obtain the fractal beam level index includes: within the local fractal window, statistically analyzing the gray level distribution of the optical brightness dimension, microwave scattering dimension, and elevation value dimension respectively; calculating the step difference frequency of the gray level distribution of each dimension, and concatenating the step difference frequencies of the three dimensions in the row direction, column direction, and layer direction to form a one-dimensional fractal difference sequence. The one-dimensional fractal difference sequence is matched step by step with the preset eight-level fractal beam energy level thresholds: if the sequence peak falls into the energy level threshold range 1, the corresponding grid beam energy level is marked as 1; if the sequence peak falls into the energy level threshold range 2, the corresponding grid beam energy level is marked as 2; and so on, until energy level 8. After mapping, a unique fractal beam energy level mark is generated for each sparse candidate grid. The fractal beam energy level marks of the sparse candidate grids are written into the index table in ascending order of row index, column index, and layer index. The index table contains three columns of position index and one column of energy level mark, for a total of four columns. After writing, the fractal beam energy level index is obtained.
[0010] Furthermore, in step 3, the process of deeply mapping and binding the peak ridges of the elevation profile with the high-energy-level regions includes: performing a differential scan of the digital elevation model (DEM) along the row and column directions, recording grids whose elevation values are greater than the elevation values in the four adjacent directions and whose elevation difference exceeds 3 meters as local peak grids; aggregating all local peak grids according to eight-directional connectivity to obtain continuous ridge segments; if the segment length is less than 5 grids, it is discarded; otherwise, it is retained and the peak ridge set is output; reading the fractal beam energy level index obtained in the previous operation; if the fractal beam energy level is greater than or equal to 6, the corresponding grid is marked as a high-energy-level grid; using the high-energy-level grid as a seed, according to eight-directional connectivity... Expand to form high-energy regions; if the area of a high-energy region is less than 9 grids, it is discarded; otherwise, it is retained and the set of high-energy regions is output; calculate the center coordinates of each peak ridge segment; search for high-energy regions in the set whose Euclidean distance from the center coordinates is no more than 5 grids, pair them up to generate pre-binding pairs; calculate the absolute value of the elevation difference for each pre-binding pair; if the absolute value of the elevation difference is no more than 8 meters, confirm the binding and record the ridge segment-high-energy region binding pair; otherwise, abandon the binding; if a single peak ridge segment and multiple high-energy regions meet the binding conditions, the high-energy region with the largest area is selected first; after binding is completed, output the depth-mapped binding list.
[0011] Furthermore, in step 3, the process of combining optical remote sensing images with interlaced ridge features includes: numbering the optical remote sensing images in ascending order by acquisition date, defining odd-numbered data as the first time-phase sequence and even-numbered data as the second time-phase sequence; freezing the optical brightness dimension data of the second time-phase sequence in the first time-phase sequence, and freezing the optical brightness dimension data of the first time-phase sequence in the second time-phase sequence, alternately synthesizing them into an interlaced spectral sequence; calculating the absolute value of the inter-phase optical brightness difference for each grid in the interlaced spectral sequence; if the absolute value exceeds 20 levels, the grid is marked as a high-gradient grid; otherwise, it is marked as a low-gradient grid. Gradient grid; Read each peak ridge segment in the depth mapping binding list and use its centerline direction as the reference direction; In the interlaced spectral sequence, starting from the endpoint of the peak ridge segment, search for 2 grids on each side of the reference direction; If a grid in the search path meets the following conditions, it is added to the ridge grid set: The grid is a high gradient grid; The absolute value of the difference between the optical brightness of the grid and the optical brightness of the endpoint grid is between 10 and 30 orders; The absolute value of the fractal beam energy level difference of the grid does not exceed 1; If the search path encounters a high energy level region boundary or a low gradient grid at any point, it immediately stops extending in that direction.
[0012] Furthermore, in step 3, the process of outputting the automatic identification result of the Daigu landform when the self-consistent difference rate of change meets the termination criterion includes: calculating the connected components of the ridge-type raster set; if the increase rate of the connected component area is less than 5%, the record expansion tends to saturate; when the increase rate of the connected component area after two consecutive extensions is less than 5%, the ridge feature direction consistency extension operation is terminated, and the interwoven extended ridge-type raster group is output; within the coverage area of the interwoven extended ridge-type raster group, the absolute difference of the raster values of the current iteration and the previous iteration is calculated for the optical brightness dimension, microwave scattering dimension, and elevation value dimension, respectively, and the raster values of the three dimensions are... The arithmetic mean of the absolute differences is used to obtain the self-consistent difference value of each grid cell. The arithmetic mean of the self-consistent difference values of all grid cells is used to obtain the average self-consistent difference of all grid cells in this iteration. The absolute difference between the average self-consistent difference of all grid cells in the previous iteration and the average self-consistent difference of all grid cells in the current iteration is recorded as the average self-consistent difference change rate. If the average self-consistent difference change rate is less than 0.02, it is considered a successful convergence. Otherwise, the interleaved spectral sequence construction operation is resumed and the iteration is restarted. When the average self-consistent difference change rate is less than 0.02 for two consecutive iterations, the termination criterion is met, and the entire iteration process is stopped.
[0013] This invention provides an automatic identification method for Daigu landforms based on remote sensing and DEM, which has the following beneficial effects: First, grid registration completely eliminates systematic errors in projection, resolution, and raster origin of multi-source images, enabling subsequent algorithms to perform calculations based on a one-to-one correspondence at the pixel level, fundamentally ensuring the consistency of spatial positioning. Second, heterogeneous noise collaborative suppression utilizes complementary features of optical, microwave, and elevation dimensions to significantly reduce random high-frequency interference while preserving true ridge details, providing clean input for subsequent texture analysis. Third, fractal beam energy level mapping transforms complex multi-dimensional texture energy into single-scale discrete markers, solving the feature imbalance problem caused by dimensional inconsistencies, and accurately locking the energy peaks of the ridge top and ridge edges through high-threshold screening. Subsequently, depth mapping pairing pairs geometric elevation ridge lines with high-energy-level regions, ensuring that the ridge skeleton simultaneously satisfies the dual physical meaning of spatial convex ridges and energy accumulation. Through interlaced spectral extension, the algorithm fully utilizes multi-temporal brightness differences to radiate and enhance the ridge flanks, significantly improving robustness under seasonal vegetation changes. Finally, an adaptive termination criterion using dual thresholds—the increase in the area of the connected domain and the average self-consistent difference rate—avoids reliance on empirical parameter tuning and ensures automatic convergence under different resolutions and noise backgrounds. The overall process of this invention combines spatial rigor, spectral synergy, and temporal consistency, significantly improving the positioning accuracy of the Daigu landform ridge, reducing the cost of manual intervention, and providing highly reliable and efficient data support for the research and resource assessment of related isolated ridge-like landforms. Attached Figure Description
[0014] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.
[0015] Figure 1 A schematic diagram of the method flow for an automatic identification method of Daigu landform based on remote sensing and DEM provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the cube structure of remote sensing DEM grid data provided in an embodiment of the present invention; Figure 3 A schematic diagram of a partial fractal window with a side length of 5 grids provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of the eight-level energy level classification system based on sequence peaks provided in the embodiments of the present invention in the coordinate system. Detailed Implementation
[0016] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.
[0017] refer to Figure 1 An automatic identification method for Daigu landform based on remote sensing and DEM, the method includes: Step 1: Through multi-source remote sensing image acquisition, simultaneously acquire optical remote sensing images, microwave remote sensing images and publicly available digital elevation models (DEMs) covering the target area, and perform grid registration processing on the three types of data to construct a grid data cube of remote sensing DEM data with unified resolution, including three dimensions: optical brightness dimension, microwave scattering dimension and elevation value dimension. In the automatic identification method for Daigu landform based on remote sensing and DEM described in this invention, step 1 first focuses on the synchronous acquisition and grid registration processing of three types of data: optical remote sensing imagery, microwave remote sensing imagery, and digital elevation model (DEM). The core objective is to provide a remote sensing DEM grid data cube with complete three-dimensional vectors, consistent spatial indexes, and uniform resolution for subsequent processing. Specifically, during the data acquisition phase, multiple constellation sensors with consistent orbital coverage are scheduled within the same time window to acquire optical and microwave remote sensing imagery of the target area in parallel. Simultaneously, DEM data for the corresponding spatial range is downloaded from a publicly available DEM database, keeping the time synchronization error within 3 hours to reduce the interference of transient surface changes on subsequent processing. After synchronous acquisition, the grid registration process immediately begins. First, all three types of data are projected onto the same geodetic coordinate reference frame and a unified map projection plane is used to ensure consistency in longitude, latitude, and elevation datum. Based on this, using the spatial resolution of optical remote sensing imagery as the baseline resolution, resampling is performed on microwave remote sensing imagery and digital elevation model (DEM) to ensure their pixel sizes are strictly consistent with those of the optical remote sensing imagery. After resolution normalization, the raster origin coordinates of each of the three types of data are calculated, and their row and column index systems are compared to detect relative displacement. If an integer pixel-level offset occurs between microwave remote sensing imagery and optical remote sensing imagery, or between DEM and optical remote sensing imagery, integer pixel-level translation is performed on microwave remote sensing imagery and DEM according to integer pixel distances to ensure that the raster origin coordinates of the three are completely coincident.
[0018] To avoid holes caused by inconsistencies in the edge ranges of different data sources, cropping was performed based on the largest overlapping area covered by the three types of data, ensuring that the remote sensing DEM grid data cubes were completely consistent in the number of rows, columns, and spatial extent. After unifying the resolution and origin position, for zero-value pixels caused by translation or cropping, grayscale values were reconstructed using bilinear interpolation on the optical brightness and microwave scattering dimensions, and elevation values were reconstructed using linear interpolation on the elevation dimension, ensuring that all raster positions had valid values. Subsequently, using the optical remote sensing image as a reference, the residual translation error of the microwave remote sensing image and the digital elevation model (DEM) relative to the reference image at the sub-pixel scale was monitored using image block matching methods. Based on the calculation results, sub-pixel-level translation corrections were performed on both, further improving the accuracy of the raster correspondence. Finally, at least 10% of all raster pairs were randomly sampled at the three-dimensional raster level to verify the consistency of row and column indices and spatial coordinates in the optical brightness, microwave scattering, and elevation dimensions. If the consistency was successful, the three types of data were formally merged into a remote sensing DEM grid data cube containing data in the optical brightness, microwave scattering, and elevation dimensions. This grid data cube not only ensures the consistency of spatial resolution and raster alignment but also lays a reliable data foundation for subsequent complex processing such as heterogeneous noise suppression, sparse thinning, and fractal beam mapping. This ensures that the Daigu landform features can be analyzed without differences across dimensions, thereby improving the accuracy and stability of the automatic identification results.
[0019] Step 2: Perform heterogeneous noise suppression processing on the data of each dimension of the remote sensing DEM grid data cube to obtain the preprocessed grid data cube; In the automatic identification method for Daigu landform based on remote sensing and DEM of this invention, the heterogeneous noise suppression processing in step 2 is not just a simple image denoising, but a multi-source consistency reconstruction approach oriented towards maintaining the three-dimensional collaborative features of the remote sensing DEM grid data cube. Its core principle lies in utilizing the complementary relationship between the optical brightness dimension, microwave scattering dimension, and elevation dimension in the landform response mechanism. Through spatial consistency constraints and spectral-shape coupling constraints, the noise energy is controllably redistributed among the three dimensions, and random fluctuations unrelated to the real landform pattern are dynamically weakened during the iteration process. First, at the theoretical level, the optical brightness dimension is mainly affected by solar incidence and surface reflection, and random noise manifests as high-frequency point-like brightness fluctuations; the microwave scattering dimension is affected by the difference in dielectric properties between the sensor coherence speckle and the surface, and its typical noise is multiplicative speckle; the elevation dimension produces isolated spikes due to data acquisition and interpolation calculation. Although the noise mechanisms differ across the three dimensions, for a target like the Daigu landform, which is characterized by ridge lines and peak morphology, there is a natural synergistic variation in the spatial location of the ridge lines. Therefore, this invention introduces cross-dimensional spatial correlation discrimination in the algorithm design: if a grid cell exhibits abnormally high-frequency fluctuations in both the optical brightness dimension and the microwave scattering dimension, but shows no obvious gradient change in the elevation dimension, it is considered pure noise; conversely, if the three dimensions exhibit synchronous gradients at similar locations, they are retained as genuine ridge line signals.
[0020] Based on this discrimination logic, the implementation process first calculates the local noise root mean square error for each of the three dimensions, and uses 1.2 times the overall mean as a dynamic threshold to mark high-noise gratings. Then, an adaptive sliding window averaging method is used for the optical brightness dimension and the microwave scattering dimension, but the window expansion and contraction rules are controlled by the elevation gradient: if the elevation gradient is significant, the window immediately contracts to preserve ridge details; if the elevation gradient is gentle, the window expands to enhance denoising, thus achieving a balance between spectral and morphological coupling. For isolated peaks in the elevation dimension, this invention employs serial morphological opening and closing operations. Opening operations remove isolated high points, and closing operations fill negative depressions. The size of the structural element is adaptively updated according to the local texture scale of the optical brightness dimension, ensuring consistent scale across the three dimensions. To prevent multiple iterations from weakening the true fine lines, this invention introduces an energy conservation compensation mechanism. After each round of processing, the total energy of the weakened pixels is statistically analyzed and injected back into the low-frequency components according to the contribution ratio of the three dimensions, suppressing high-frequency noise while maintaining the overall brightness and scattering intensity distribution. When the iterative convergence criterion is triggered—that is, the reduction in the three-dimensional noise mean square error is less than 30% for two consecutive rounds—the method will perform a spectral-shape consistency test on the randomly sampled graticule, comparing the local correlation coefficients before and after processing. If the improvement exceeds a predetermined threshold and the ridge connectivity in the elevation dimension remains stable, then the heterogeneous noise suppression processing is considered successful, and a preprocessed homogeneous data cube is output. This preprocessed homogeneous data cube achieves uniformity in spatial resolution, graticule alignment, and noise level, providing data support for subsequent sparse thinning and fractal beam mapping that balances texture clarity and numerical accuracy. This allows the fractal features of the Daigu landform ridgeline and the summit to be seamlessly mapped across dimensions, ultimately significantly improving the reliability and anti-interference capability of automatic identification.
[0021] Step 3: Based on the preprocessed grid data cube, while maintaining the spatial consistency of the grid, a sparse candidate raster set is generated by sparse thinning. Fractal beam mapping is performed on the sparse candidate raster set to obtain the fractal beam energy level index. Then, the peak ridge line of the elevation profile is deeply mapped and bound to the high energy level region. The ridge feature is extended by combining optical remote sensing image interleaving. When the self-consistent differential change rate meets the termination criterion, the automatic identification result of Daigu landform is output.
[0022] The theoretical basis for sparse thinning is the complementarity of Nyquist sampling and fractal dimension. When the known target terrain has strong self-similarity features in the dimensions of optical brightness, microwave scattering, and elevation, the full grid does not provide additional information density. Instead, it dilutes the feature variance due to high redundancy. Therefore, the algorithm first performs reversible thinning on the preprocessed grid data cube with a row spacing of 3 grids, a column spacing of 3 grids, and a layer spacing of 2 layers to obtain a sparse candidate grid set. The thinned set carries a complete spatial topological skeleton but significantly reduces the computational complexity of subsequent fractal feature analysis. Subsequently, fractal beam mapping is introduced, treating the cube window with a side length of 5 grids around each sparse candidate grid as a local fractal environment. The step difference frequencies of the optical brightness, microwave scattering, and elevation dimensions in the row, column, and layer directions are counted respectively. A one-dimensional fractal difference sequence is constructed by cascading, and then matched using an eight-level fractal beam energy level threshold. This transforms the originally incomparable multidimensional texture differences into single-scale energy level labels. The theoretical support for fractal beam mapping comes from the energy level partitioning law of fractal geometry: texture roughness at different scales can be described by finite discrete levels, and the Daigu landform ridge naturally corresponds to the high energy level region because it simultaneously presents high brightness gradient, high scattering gradient and high elevation gradient in three-dimensional space.
[0023] To extract potential ridges from the three-dimensional energy level field, the algorithm finds the peak ridges in the elevation profile through differential scanning in the elevation dimension. Then, it uses fractal beam energy level indexing to mark high-energy-level regions. The two are then deeply mapped and bound together under the conditions that the center-to-center distance does not exceed 5 grids and the absolute value of the elevation difference does not exceed 8 meters. This binding process embodies the dual consistency principle of "geometric constraints + energy level constraints": geometric constraints ensure that the ridges and high-energy-level regions are spatially aligned, while energy level constraints guarantee concentrated texture energy. Next, optical remote sensing image interleaving is introduced to extend the ridge features. An interleaved spectral sequence is constructed by alternating odd and even phase freezing, generating a high-gradient indicator field between phases. Using the centerline direction of the peak ridge segment as a reference, two grids are searched on both sides. Grids that simultaneously satisfy the conditions of high gradient, high brightness difference, and low energy level difference are incorporated into the ridge grid set. The theoretical significance of interlacing extension lies in enhancing ridge contours using multi-temporal radiation differences. Its fundamental assumption is that the morphological stability of the Daigu landform during seasonal changes leads to the true ridge lines exhibiting brightness difference boundaries across different temporal phases, while noise exhibits temporal randomness. As the ridge grid set continuously expands, the algorithm monitors the increase rate of connected region area in real time. When the increase rate is less than 5% for two consecutive rounds, the lateral gain is considered saturated. At this point, a self-consistent difference mechanism is used to evaluate the consistency of the 3D data in the current ridge mask across iterations.
[0024] The self-consistent difference method unifies the absolute differences of the optical brightness dimension, microwave scattering dimension, and elevation dimension into an average self-consistent difference, emphasizing cross-dimensional cooperative stability. The threshold for the average self-consistent difference change rate is set to 0.02, meaning that if the relative change between two rounds of differences is less than 2%, it is considered a single convergence; two consecutive rounds of convergence indicate that the overall system has reached a steady state. This threshold is not an empirical parameter but is derived from the statistical distribution of noise convergence rate based on Markov chain Monte Carlo simulations, which can significantly suppress high-frequency noise while ensuring ridge integrity. After the termination criterion is triggered, the algorithm freezes the ridge grid set, the high-energy-level region set, and the depth mapping binding list, and outputs the union of these three as the automatic identification result layer for the Daigu landform. This result not only preserves the spatial continuity of the ridge lines but also provides target weights for subsequent quality screening and sample training through fractal energy level attributes.
[0025] Furthermore, in the automatic identification method for Daigu landform based on remote sensing and DEM of this invention, grid registration processing is the starting point for all subsequent spatial-spectral-morphological fusion operations. Its principle lies in converting data from three sources—optical remote sensing images, microwave remote sensing images, and digital elevation models (DEMs)—with completely different physical quantity dimensions into rasterized three-dimensional data blocks with unified longitude, latitude, and elevation references through strict spatial benchmark unification and pixel-level alignment. Firstly, optical and microwave remote sensing images differ in their respective orbital inertial coordinate systems, attitude calculation models, and imaging scanning methods during acquisition. Digital elevation models (DEMs) typically originate from post-processing results of laser altimetry, stereo image pairs, or radar interferometry, and their native coordinate methods are often inconsistent. If these three types of data are directly input into the subsequent fractal beam mapping and depth mapping binding process, any pixel-level misalignment will be amplified during the iteration process, leading to positional drift or jagged edges in the ridge identification results. To avoid this problem, the grid registration process first performs coordinate transformation on the three types of data separately, projecting them onto the same geodetic coordinate reference frame and using the same map projection plane. At this point, all pixels are on the same reference ellipsoid in longitude and latitude coordinates. However, due to differences in sensor resolution, the pixel sizes are not consistent. Therefore, the spatial resolution of the optical remote sensing image is used as the reference resolution, and the spatial resolution of the microwave remote sensing image and the digital elevation model (DEM) is adjusted to have completely consistent pixel sizes through bilinear or cubic convolution resampling methods.
[0026] Optical remote sensing imagery serves as the baseline resolution for two main reasons: firstly, it boasts high radiometric resolution and rich detail; secondly, the Daigu landform, characterized by ridgelines and peaks as key features of moderate topographic relief, provides the most intuitive depiction of landform boundaries through its light and shadow distribution within the visible light range. Using other dimensions as the baseline would result in the loss of detail in the optical brightness dimension during resampling, thus reducing the energy level resolution of subsequent fractal beam mapping. After resolution normalization, grid registration further calculates the raster origin coordinates for each of the three data types—the absolute positions of pixels with row and column indices of zero on the map projection plane. Since different data sources may use different map sheet starting points during preprocessing and map clipping, when there are integer-pixel differences in the raster origin coordinates, the pixel grid will exhibit methodological translational misalignment. Therefore, if any integer-level offset is detected between microwave remote sensing imagery and optical remote sensing imagery, or between digital elevation model (DEM) and optical remote sensing imagery, the grid registration process immediately performs an integer translation according to the whole grid pixel spacing, translating the microwave remote sensing imagery and DEM as a whole to ensure that their grid origin coordinates are completely consistent with those of the optical remote sensing imagery, thus ensuring strict alignment of the three types of data in the row and column indexing system. Based on this, the grid registration process also needs to handle sub-pixel-level residual biases.
[0027] Although whole-pixel translation achieves pixel-level overlap, due to sensor imaging geometry, attitude drift, and resampling interpolation errors, microwave remote sensing images and digital elevation models (DEMs) may still have zero-point geometric differences at the sub-pixel scale compared to optical remote sensing images. To address this issue, the same-grid registration process employs block matching and a multi-resolution pyramid strategy. A feature window with uniform density distribution is selected as a reference in the optical remote sensing image. The maximum cross-correlation coefficient is searched in both the microwave remote sensing image and the DEM to estimate the sub-pixel displacement. Subsequently, a sub-pixel interpolation scheme is used to fine-tune the positions of both, ensuring that the ridge vertices of the same feature fall at the same grid intersection point in three dimensions. After completing the above process, at least 10% of all grid pairs are randomly sampled for double verification of the row and column index consistency of the three types of data and the actual coordinates of the map projection plane. If any random pair shows a longitude, latitude, or elevation deviation exceeding half a pixel, the process returns to the sub-pixel displacement estimation stage for iterative correction. Although the entire grid registration process is time-consuming, its significance lies in providing a spatially accurate and traceable remote sensing DEM grid data cube. This allows subsequent processes such as heterogeneous noise suppression, sparsity thinning, fractal beam mapping, and depth mapping to be carried out under the same row and column index, ensuring a one-to-one correspondence between the optical brightness dimension, microwave scattering dimension, and elevation value dimension at each pixel location. Only in this way can the brightness gradient, scattering gradient, and elevation gradient of the Daigu landform ridge achieve a true coordinated response in the fractal beam energy level field, further providing distortion-free input conditions for three-domain synchronous resonance tracking, thereby significantly improving the positioning accuracy and morphological integrity of the Daigu landform automatic identification.
[0028] Furthermore, in the automatic identification method for Daigu landform based on remote sensing and DEM of this invention, heterogeneous noise suppression processing is based on the three-dimensional synergistic characteristics of the same-grid data cube of remote sensing DEM. It achieves graded reduction of random noise energy in the optical brightness dimension, microwave scattering dimension, and elevation dimension through a combination of statistical discrimination, local smoothing, and morphological constraints. First, the method calculates the noise root mean square error and noise variance ratio for each of the three dimensions in their respective frequency bands and compares them with the global noise root mean square error. When the noise root mean square error of any dimension exceeds 1.2 times the overall mean, the corresponding raster is considered a high-noise raster, and its row index, column index, and layer index are recorded for subsequent directional processing. High-noise rasteres often appear as isolated bright spots or scattering spots in the optical brightness and microwave scattering dimensions. If they directly participate in subsequent fractal beam mapping, it will cause the energy level threshold calculation to shift upwards and mask the true ridge texture. To address this, the algorithm employs an adaptive sliding window averaging method to smooth the high-noise grid in both the optical brightness and microwave scattering dimensions. The initial window size is consistent with the baseline resolution of the remote sensing DEM grid data cube, and then automatically expands or shrinks based on the local brightness or scattering gradient: when the grayscale difference between the central grid and its neighbors exceeds one standard deviation of the background mean, the window expands to three times the pixel size to improve noise reduction; when the central grid is located in a significant gradient zone, the window shrinks to two pixels to protect the high-frequency information of the Daigu landform ridgeline edge.
[0029] After smoothing the optical brightness and microwave scattering dimensions, the method performs sequential morphological opening and closing operations on the elevation dimension at the same grid location to reduce isolated elevation spikes caused by stereo measurement errors or radar interference residuals. The opening operation first uses square structuring elements to remove isolated high points, and the closing operation then uses structuring elements of the same size to fill negative depressions, achieving overall smoothness of the elevation surface without sacrificing the continuity of the Daigu landform ridge. After joint suppression of the three dimensions, the algorithm recalculates the noise root mean square error. If the reduction in any dimension does not reach 30%, it indicates that residual noise still significantly interferes with the local energy level distribution. The method then maintains the same window size and performs the sliding window averaging method again, repeating the morphological operations until the noise root mean square error reduction meets the 30% condition or the number of iterations reaches a preset upper limit. At the end of each iteration, the method randomly samples at least 10% of all grids to compare the kurtosis and skewness of the optical brightness, microwave scattering, and elevation distributions before and after processing, to confirm that the main texture remains unchanged and to promptly identify the risk of oversmoothing. Furthermore, to avoid brightness drift caused by window iteration, energy conservation compensation is introduced. The weakened high-frequency energy is redistributed back to the low-frequency components according to the dimensional contribution ratio, thus ensuring the true consistency of the entire image in radiometrics. Through the above process, the preprocessed gridded data cube achieves synergistic reduction in noise levels across the three dimensions. Random high-frequency textures in the optical brightness and microwave scattering dimensions are significantly suppressed, isolated spikes in the elevation dimension are effectively reduced, and the mid-to-low-frequency structure of the Daigu landform ridges and peak morphology is fully preserved. This processing result provides a reliable, smooth, and realistically characteristic input for subsequent sparse thinning, fractal beam mapping, and depth mapping binding, laying a data foundation with controllable noise and sufficient information for the entire Daigu landform automatic identification process.
[0030] Furthermore, the preprocessed gridded data cube maintains the same row and column indices in the optical brightness, microwave scattering, and elevation dimensions, resulting in a number of complete grid cells often reaching tens of millions. Directly participating in fractal difference analysis not only incurs enormous computational costs but also leads to average dilution of texture energy due to high redundancy. To reduce computational complexity without disrupting the ridge topology of the Daigu landform, the method first reads the number of rows, columns, and layers of the preprocessed gridded data cube. The sampling period is set to 3 grid cells in the row direction, 3 grid cells in the column direction, and 2 layers in the layer direction. This period value is derived from the statistical results of typical ridge distances in the Daigu landform: the average half-power width of the three-dimensional ridge line in the row and column plane is approximately 7 grid cells, while the maximum abrupt change interval of the elevation ridge top in the layer direction does not exceed 3 layers. Therefore, using odd row and column intervals and even layer intervals can cover the ridge peak core while avoiding sampling aliasing. When the 3D dimension is divided by the interval and a remainder occurs, the method uniformly adopts an upward rounding and padding strategy, adding an equivalent virtual grid at the boundary to maintain the consistency of the same grid space and ensure the alignment of the traversal pointers.
[0031] Subsequently, a row-first traversal is performed at each layer. When both the row and column indices satisfy the modulo operation condition of 3 grid cells in both the row and column directions, the method marks the corresponding grid as a sparse grid. This mark is stored in the memory bitmap and also written to an independent index queue in real time, preparing for subsequent fast random access. After the traversal, a set of sparse grids is obtained, with a density approximately one-ninth that of the original planar grid. Based on this set, the method performs further selection by layer index from low to high, retaining grids that meet the conditions at a period of 2 layers in the layer direction, and discarding the rest, thus obtaining the final sparse candidate grid set. To ensure instruction-level parallelism for subsequent multi-threaded fractal beam computation, the method writes three-dimensional coordinate marks of the row, column, and layer indices to each retained grid in real time, and maintains row and column slice pointers in main memory, allowing different threads to access their respective slices without locks. This sparse candidate raster set achieves three key objectives: First, by using 332 sampling to simultaneously reduce the number of rasteres in the row-column plane and layer direction, the original data volume is compressed to about one-sixth, significantly shortening the execution time of fractal difference operations and eight-level energy matching; Second, by using odd intervals to form a stable checkerboard distribution on the plane, it ensures that at least one raster is sampled for any adjacent ridge top or slope shoulder, maintaining the geometric continuity of the Daigu landform ridge line; Third, by rounding up and filling with virtual rasteres, it supports the mirror expansion of subsequent local fractal windows at the boundaries, avoiding the impact of sampling holes on window integrity. The final sparse candidate grid set serves as the input skeleton for fractal beam mapping. It not only faithfully reconstructs the network backbone of the Daigu landform ridge in terms of spatial coverage, but also provides controllable and efficient data slices for multi-dimensional synchronous processing at the algorithm level. This lays a foundation for high accuracy and high performance in subsequent key steps such as energy level index generation, peak ridge line extraction from elevation profiles, and interlaced extended ridge features. It ensures that the automatic recognition of Daigu landform can maintain real-time processing speed and sub-pixel-level positioning accuracy even in large-scale data scenarios.
[0032] Furthermore, after the sparse candidate grid set is constructed, the method must verify the fidelity of its spatial sampling to the topology of the Daigu landform ridge. Therefore, a random sampling neighborhood integrity check mechanism is introduced as a dynamic feedback loop. The algorithm first uses a true random number generator to extract no less than 5% of the total number of grid pairs in the sparse candidate grid set. For each grid pair, eight-directional neighboring grids are retrieved in the preprocessed same-grid data cube. The number of missing neighboring grids is counted and the maximum missing value is recorded. The integrity of the eight-directional neighboring grids is a direct indicator of whether the sampling step size destroys the continuity of the local morphology. Because the elevation ridge top of the Daigu landform ridge line usually transitions between the row and column planes and the micro-slope shoulders and ridge top slopes in a single-pixel gradient manner, if any grid pair has two or more holes in the eight-directional neighborhood, it means that the current row direction thinning interval or column direction thinning interval has destroyed the local Nyquist sampling conditions of the area, which will cause the subsequent fractal difference sequence to have jump peaks, thereby distorting the fractal beam energy level threshold matching results. Therefore, if the number of missing neighboring rasters exceeds one, the method immediately reduces the row or column thinning interval by one raster and re-executes the sparse thinning process to obtain a new set of sparse candidate rasters and re-triggers random sampling inspection. This process can be repeated multiple times until all sampled raster pairs meet the threshold requirement that the number of missing neighbors does not exceed one.
[0033] The dynamic step size adjustment embodies the adaptive balance principle between "sampling density and computational complexity": increasing sampling density in densely curved ridge regions to preserve details, and maintaining the original sparsity in smooth regions to save computational resources. After neighborhood integrity is verified, the algorithm establishes a cubic local fractal window for each grid in the sparse candidate grid set. The window side length is fixed at five grids, ensuring that the window contains neighborhood information extending two grids in each of the row, column, and layer directions. Setting the five-grid window is an empirical value derived from the statistical average half-power width of the Daigu landform ridge, which can provide sufficient texture samples for fractal beam mapping while maintaining the sharpness of the ridge crest and the smoothness of the slope shoulder. If the cubic window touches the boundary of the same grid, the method uses a mirroring method to fill the window grid, that is, mirroring the extension position with the real grid value near the edge of the window to ensure the integrity of the window shape without introducing artificial discontinuities caused by zero padding. Mirror completion and random sampling verification together form a dual constraint of "boundary integrity - internal consistency," ensuring that each sparse candidate raster is in a local fractal environment with complete data and uniform scale. This guarantees the comparability of subsequent step difference frequency statistics in the dimensions of optical brightness, microwave scattering, and elevation. Through this series of self-checking, backtracking, and mirror expansion operations, the sparse candidate raster set not only maintains the coherence of the Daigu landform ridge network but also statistically meets the strict requirements of fractal analysis for local window integrity. Ultimately, this lays a reliable foundation for the accurate calibration of fractal beam level indexes and the geometric matching of depth mapping, ensuring that subsequent automatic identification of Daigu landforms can robustly extract the true landform skeleton against a background of multi-source data noise.
[0034] Furthermore, for the sparse candidate raster set that has passed the neighborhood integrity check, the method first constructs a cubic local fractal window with a side length of 5 grates around each raster. Inside this window, the gray-level distributions in the optical brightness dimension, microwave scattering dimension, and elevation dimension are statistically analyzed, and a discrete histogram is formed with the pixel gray-level value or elevation value on the horizontal axis and the occurrence frequency on the vertical axis. The histogram itself is only a display of probability density. To reveal the multi-scale jump characteristics of local texture, the method performs step difference on the histogram, counting the frequency difference between adjacent gray-levels in the form of "positive and negative one-steps", thereby obtaining the step difference frequency vector for each dimension. The step difference frequency vector mathematically corresponds to the first discrete derivative of local fluctuations, which can amplify the weak contribution of ridge or shoulder regions to the gray-level distribution. To ensure the comparability of three-dimensional difference information within the same framework, the algorithm follows the spatial order of row, column, and layer directions, sequentially concatenating the step difference frequencies of optical brightness, microwave scattering, and elevation values to obtain a one-dimensional fractal difference sequence. This sequence maintains the correspondence between rows, columns, and layers in terms of location and preserves the coupling of three-dimensional amplitudes in terms of value domain. Therefore, it not only preserves the comprehensive slope aspect information of the Daigu landform ridgeline in three-dimensional space but also avoids the bias caused by scale inconsistencies between dimensions.
[0035] Next, the method prepares a set of eight fractal beam energy level thresholds. These eight thresholds are obtained by discretizing the fractal dimension distribution of a large number of Daigu landform samples through statistical regression, representing continuous energy segments from background noise to strong ridge texture. After amplitude normalization, the one-dimensional fractal difference sequence is matched stepwise with the threshold table in the beam space: if the sequence peak falls into energy level threshold interval 1, the method marks the sparse candidate grid as fractal beam energy level 1, indicating that the three-dimensional texture fluctuation of the grid just crosses the lower limit of the lowest energy level; if the sequence peak falls into energy level threshold interval 2, it is marked as fractal beam energy level 2; and so on, until the peak falls into the highest threshold interval 8 and is marked as fractal beam energy level 8. Here, the "peak" is not simply the maximum value of the sequence, but rather the first local maximum point on the sequence energy density curve that satisfies the threshold interval closure condition and has the gradient signs reversed on both sides, so that the matching result can reflect the true dominant scale of the local texture. In this way, each sparse candidate grid is assigned a unique fractal beam energy level label. A larger label value indicates a more complex multidimensional texture and more dramatic slope aspect changes within the grid's window, making it more likely to correspond to the ridge crest or edge of the Daigu landform. To facilitate the subsequent depth mapping between peak ridge lines in the elevation profile and high-energy-level regions, the method writes the fractal beam energy level labels of the sparse candidate grids into an index table in ascending order of row index, column index, and layer index. This index table strictly maintains a four-column structure: three columns of position indices and one column of energy level labels, without any redundant fields, ensuring a one-to-one mapping between the fractal beam energy level domain and the physical space domain. Once generated, the index table constitutes the fractal beam energy level index, which serves as an important basis for detecting high-energy-level regions and also provides initial energy distribution values for subsequent probabilistic refraction feedback. Overall, fractal beam mapping, by statistically analyzing the gray-level distribution within a local fractal window, calculating the step difference frequency and concatenating them into a one-dimensional fractal difference sequence, and then matching it step by step with the eight-level fractal beam energy level threshold, endows sparse candidate grid sets with measurable energy level characteristics. This achieves cross-domain conversion from three-dimensional multi-source textures to single-scale energy markers, laying a unified and physically meaningful energy foundation for subsequent processes of automatic identification of Daigu landforms.
[0036] Furthermore, in the automatic identification method for Daigu landform based on remote sensing and DEM of this invention, the depth mapping between the peak ridge line of the elevation profile and the high-energy region is driven by the morphological constraints of the digital elevation model (DEM) and the energy constraints of the fractal beam energy level index. The core idea is to first locate the real elevation peak ridge with strict geometric difference logic, then use the fractal beam energy level field obtained by multi-dimensional texture fusion to lock the high-energy distribution area, and finally achieve one-to-one pairing through spatial proximity and elevation difference threshold, thereby establishing a ridge-shaped skeleton with clear physical meaning in the three-dimensional grid. In specific implementation, the algorithm performs a differential scan of the digital elevation model (DEM) along the row and column directions. As long as the elevation value of a certain grid exceeds the elevation values of the four adjacent directions and the elevation difference reaches or exceeds 3 meters, it can be identified as a local peak grid. This threshold takes into account the average abrupt height of the Daigu landform ridge top relative to the slope shoulder, so that the differential scan neither excessively picks up noise cusps nor misses the real ridge top. Next, all local peak gratings are aggregated according to eight-directional connectivity to obtain continuous ridge segments. If the segment length is less than 5 gratings, it is determined that it is insufficient to form a stable ridge and is discarded; otherwise, it is retained and the peak ridge set is output. This step uses grating length rather than curve geometry as the criterion, which can automatically adapt to changes in grating resolution. Subsequently, the fractal beam energy level index generated in the previous operation is read. All gratings with fractal beam energy levels greater than or equal to 6 are marked as high-energy gratings. The energy level threshold of 6 corresponds to the upper quantization segment of the fractal beam energy level quantization curve, which can effectively filter out background noise and low-energy slope textures.
[0037] Using high-energy-level grids as seeds, the algorithm expands along eight directions of connectivity to form high-energy-level regions. If the expanded region area is less than 9 grids, it is considered that the energy aggregation is insufficient to support a complete ridge and is discarded; otherwise, it is retained and the set of high-energy-level regions is output. To accurately map geometric ridges to high-energy-level regions, the system first calculates the center coordinates of each peak ridge segment, and then searches for high-energy-level regions in the set whose Euclidean distance from the center coordinates is no more than 5 grids. These regions are then paired to form pre-bonded pairs. The distance threshold is related to the previous sparse thinning step size to ensure that the search radius covers all possible energy-level peaks. For each pre-bonded pair, the absolute difference between the center elevation of the ridge segment and the average elevation within the high-energy-level region is calculated. If the absolute value of the elevation difference is no more than 8 meters, the binding is confirmed and the ridge segment-high-energy-level-region binding pair is recorded; otherwise, it is discarded to avoid confusing different terrain layers. If a single peak ridge segment matches multiple high-energy regions simultaneously within the dual threshold range of distance and elevation difference, the algorithm prioritizes binding the high-energy region with the largest area. This strategy is based on the statistical fact that the ridge crests of the Daigu landform generally correspond to the areas with the largest energy accumulation. After the above process is completed, a depth-mapped binding list is output. Each record corresponds to a unique pair of ridge segments and high-energy regions. They are collinear in space and resonate in energy, possessing both geometric convex ridge features and high fractal beam energy level attributes. This provides a stable seed for subsequent optical remote sensing image interleaving to expand ridge features, and lays the foundation for vertex-level input for three-domain synchronous resonance tracking. This enables the automatic identification of the Daigu landform to have strong robustness and repeatability in both macroscopic geometric context and microscopic energy details.
[0038] Furthermore, after completing the deep mapping binding, the system needs to further expand the ridge grid set, which is still in its skeletal form, by leveraging the radiometric differences in the time-series images. The theoretical basis of this method is that the reflectance fluctuations of the Daigu landform ridgeline under different observation dates, due to changes in solar altitude angle, vegetation humidity, and atmospheric scattering, exhibit relatively stable directional consistency, while the brightness changes of the surrounding slopes and valleys within the same time scale show more obvious randomness and phase misalignment. To fully utilize this temporal difference, the algorithm groups the optical remote sensing images numbered in ascending order by acquisition date into odd-even groups, defining odd-numbered data as the first temporal sequence and even-numbered data as the second temporal sequence. Then, an alternating freezing strategy is adopted: in the first temporal sequence, the optical brightness dimension data at the corresponding position of the second temporal sequence is fixed as a static constant; conversely, in the second temporal sequence, the optical brightness dimension data of the first temporal sequence is frozen. The resulting interlaced spectral sequence is equivalent to stretching the radiometric difference of the same grid in two temporal phases into an observable instantaneous contrast field, thereby highlighting the brightness contrast peaks at the ridgeline position caused by geometric exposure differences at the pixel level. Next, the algorithm calculates the absolute value of the interphase optical brightness difference for each grid in the interlaced spectral sequence. If the value exceeds 20 levels, it is marked as a high-gradient grid; otherwise, it is marked as a low-gradient grid. This threshold is determined by statistical results from multiple images and can shield low-amplitude noise caused by fine atmospheric particles without excessively rejecting the true ridge edge contrast.
[0039] Subsequently, each peak ridge segment in the depth mapping binding list is read, and its centerline direction is used as the reference direction. Two endpoints of the ridge segment are selected as starting points in the interlaced spectral sequence, and two grids are searched on each side of the reference direction. The search process requires simultaneous fulfillment of three conditions: First, the visited grid must be a high-gradient grid to ensure it truly represents a sudden edge of temporal difference in the interlaced spectral sequence; second, the absolute value of the optical brightness difference between the grid and the endpoint grid must be between 10 and 30 orders. This range excludes false edges with excessively low contrast and prevents extreme bright spots that may originate from cloud shadows or man-made objects due to excessively high contrast; finally, the absolute value of the fractal beam energy level difference of the visited grid must not exceed 1 to ensure that the expanded grid maintains continuity with the original ridge segment in multidimensional texture intensity. If all three conditions are met, the grid is merged into the ridge grid set and becomes the new endpoint for the next iteration; if the search path encounters a high-energy-level region boundary or the first low-gradient grid along the way, the algorithm immediately terminates the search in that direction to prevent the ridge from spreading to energy valleys or noise spots. By alternately performing left and right side searches and recursively updating endpoints, the ridge grid set gradually expands towards both wings while maintaining stability in the central direction, achieving integrated reinforcement of the ridge boundary brightness difference and fractal beam energy level slope. Since the search depth is limited to 2 grids, the system can quickly capture effective information from the ridge flanks while avoiding the misclassification of excessively distant slope shoulders or valley bottom pixels into the ridge range. After parallel expansion of multiple peak ridge segments, the ridge grid set not only obtains continuous edge bands on the plane but also maintains a highly consistent energy distribution in both the fractal beam energy level domain and the optical brightness difference domain. This provides complete data support for subsequent self-consistent differential evaluation and greatly improves the robustness of Daigu landform ridge identification in multi-temporal spectral backgrounds. It enables the automatic identification results to stably adapt to seasonal changes, vegetation disturbance, and atmospheric condition differences, while ensuring that the geometric coherence between the ridge crest and the valley edge is not disrupted.
[0040] Furthermore, after completing the interwoven and extended ridge features, the system needs to determine whether the current ridge grid set contains sufficiently complete Daigu landform information and can generate the final result accordingly. To this end, this invention first performs connected component analysis on the ridge grid set within a three-dimensional grid framework, aggregating all ridge grids based on eight-directional connectivity and calculating the connected component area. If the increase in connected component area is less than 5% compared to the previous iteration, it indicates that the lateral expansion's gain on the ridge skeleton has entered a marginal weakening zone, and the system records the "expansion approaching saturation" state accordingly. When the increase in connected component area is less than 5% for two consecutive iterations, it can be determined that the consistent extension of the ridge feature direction can no longer significantly expand the ridge coverage. Therefore, the current interwoven and extended ridge grid set is immediately frozen, and any search operations based on the reference direction are stopped. This processing logic originates from the spatial self-similarity characteristics of the ridge lines in the Daigu landform: at the grid scale, once the true ridge edge is completely covered, subsequent outward expansion only brings information about slopes and valleys, which weakens the energy concentration of the ridge line. Therefore, double stopping by the area and slope of the connected domain can prevent the ridge grid set from being disorderly expanded by noise. After entering the self-consistent evaluation stage, the algorithm only calculates the absolute difference between the corresponding grid values of the current iteration and the previous iteration for the spatial range covered by the interwoven and expanded ridge grid group, in terms of optical brightness, microwave scattering, and elevation. Then, the arithmetic mean of the absolute differences in the three dimensions is taken to obtain the self-consistent difference value of each grid, thus simultaneously measuring the overall consistency of radiation change, scattering change, and elevation change. The arithmetic mean of the self-consistent difference values of all ridge grids is taken again to obtain the average self-consistent difference of all grids in this iteration. The system compares this value with the corresponding index of the previous iteration, and the absolute difference is the average self-consistent difference rate of change.
[0041] When the average self-consistent difference rate of change is below 0.02, it indicates that the comprehensive difference amplitude of the three dimensions within the ridge mask has stabilized, and the system is considered to have successfully converged once. If the rate of change is still above the threshold, it indicates that some grids are not yet stable in terms of spectrum, scattering, or elevation. The algorithm will return to the interlaced spectral sequence construction stage and restart the interlacing freezing, high gradient labeling, and lateral search processes to further correct the residual anomaly regions in the ridge grid set. The design of this iterative mechanism aims to compress the residual differences between multi-source radiation and three-dimensional terrain to the sub-threshold level, so as to ensure that the final ridge grid is not only geometrically complete but also highly self-consistent in the multi-dimensional information field. If the average self-consistent difference rate of change is below 0.02 for two consecutive iterations, the termination criterion is met, and the iteration process is terminated immediately. At this time, the interlaced extended ridge grid set has dual characteristics of spatial steady state and radiation steady state and can be directly output as the automatic identification result of Daigu landform. The system then writes the grid group, along with the previously recorded row, column, and layer indices, as well as the fractal beam energy levels, into the results list, forming a complete output that integrates spatial coordinates, energy level attributes, and multidimensional radiation consistency. This provides a traceable data foundation for subsequent quality diagnosis, lateral comparison, and secondary training of machine learning. By using dual-threshold stopping judgment based on connected region area increment monitoring and average self-consistent differential rate of change, this invention effectively avoids the risks of overfitting noise or premature stopping that may miss ridge lines, ensuring that the automatic identification results of the Daigu landform achieve the optimal balance between accuracy and stability.
[0042] The following example demonstrates a complete and repeatable implementation of an automatic topographic identification method for Daigu based on remote sensing and DEM. The study area's center coordinates are set at 35.425°N, 117.825°E. The projected coordinate system is EPSG4490 with equal latitude and longitude Mercator. The planar area is 5000m × 5000m, corresponding to a pixel side length of 10m for optical remote sensing images, 20m for microwave remote sensing images, and 30m for the digital elevation model (DEM). Four optical images are used, and the microwave images were acquired synchronously within 3 hours. The DEM is sourced from NASA SRTM version 3 and supplemented with LiDAR.
[0043] First, the three types of data are projected onto a unified projection and resampled to a 10m reference pixel. The number of raster rows after projection is then determined. Column number Microwave image column offset was detected during grid origin alignment. Pixels, DEM row offset Pixels, therefore the microwave image was shifted 10m to the left and the DEM was shifted 20m down, making the three-dimensional cube Fully aligned, , , ,in For optical brightness dimension, For microwave scattering dimensions, Elevation is the dimension.
[0044] The heterogeneous noise suppression step first calculates the noise mean square error for each dimension. The calculation result is , , Overall mean Both the optical and microwave dimensions exceed the threshold. The corresponding pixel is marked as a set of high-noise raster cells. .exist Initial side length of the application window Adaptive moving average When the gray level difference between the center pixel and the mean gray level of its neighborhood is... Exceeding the standard deviation of the dimension At that time, Zoom in After processing the two dimensions, morphological opening and closing operations are performed on the elevation dimension. A square structuring element is used. Perform opening operations first. Then perform the closing operation. ; obtain updated elevation layers After one iteration, the root mean square error of the three-dimensional noise decreased to [value missing]. , , The decreases all exceeded 30%, and the cube was updated to .
[0045] Sparse candidate raster sampling with row thinning interval Column thinning interval thinning interval .Will Complete , Complete Construct a set ; can be obtained Grid. Random sampling For each grid pair, if any pair of eight neighboring cells has more than one hole, the value is reduced. Resample; calculate the maximum missing value as 1, and maintain the original sampling parameters. Establish a new sampling method at the center of each sparse raster. window histogram of optical brightness dimension In grayscale Keys are used, frequencies are used as values. Step difference Similarly, calculate , Serialized in row, column, and layer order to obtain .
[0046] right Normalize to [0, 100] and record the peak value. Let the energy level threshold sequence be... ;like Then the empowerment level Example sparse raster Calculated Therefore The index table entries are formatted as follows: .
[0047] exist Layer-to-elevation four-neighbor difference like Rice Store peak set Connectivity clustering yielded 214 ridge segments, of which 197 were retained with a length ≥ 5 pixels. (Energy level index table) The raster forms a high-energy-level seed set, which is continuously expanded to obtain 93 high-energy-level regions, retaining 80 regions with an area ≥ 9 pixels. For the center of each ridge segment... and the high-regional mass centers of each energy level calculate Simultaneously check The bindings were generated based on the satisfied data, resulting in a total of 164 depth mapping binding records. Optical images were sorted by date and grouped into odd and even groups, with odd groups... Even array Alternating freezing yields interlaced brightness. .calculate like mark ,otherwise For each bound ridge segment, take the centerline direction vector. At the endpoint along Each side moves forward 2 steps with a step size of 1 pixel. If the detected pixel satisfies... Then add it to the ridge grid set. Continue iterating until both the left and right sides encounter... Or the boundary of a high-energy-level region.
[0048] At the end of each iteration, the area of the ridge-shaped connected region is calculated. For example, the results of the first three rounds. Pixel Pixel Pixel, then The second round reaches saturation, and the third round again falls below 5%, triggering a halt to lateral expansion. This applies to each pixel within the ridge mask. calculate Summary average ; like Convergence is counted once. Actual measurement... , , ,correspond The third round of convergence was successful, and the fourth round continued iterating. , The iteration terminates after two consecutive convergences. The final size of the ridged grid set is... Pixels, along with energy level indices Write the 3D coordinates into the result file: .
[0049] in , For row and column indexes (unit cell). A value of 0 indicates that the output surface is located in the optical brightness layer. This file contains fractal beam energy level labels. It represents the automatic identification results of the Daigu landform and can be directly overlaid onto a GIS platform for visualization or used as training labels for deep learning models.
[0050] like Figure 2As shown, the remote sensing DEM gridded data cube of this invention is a unified resolution three-dimensional data structure constructed by performing gridded registration processing on three types of data: optical remote sensing imagery, microwave remote sensing imagery, and digital elevation model (DEM). This data cube includes three dimensions: optical brightness, microwave scattering, and elevation. During the construction of the data cube, the optical remote sensing imagery, microwave remote sensing imagery, and DEM are first projected onto the same geodetic coordinate reference frame and the same map projection plane to ensure that the three types of data have a unified longitude, latitude, and elevation reference. The spatial resolution of the optical remote sensing imagery is used as the reference resolution. If the resolution of the microwave remote sensing imagery or the DEM is inconsistent with the reference resolution, its spatial resolution is adjusted to be consistent with the optical remote sensing imagery through resampling. After unifying the resolution, the grid origin coordinates of the three types of data are calculated. If there is a misalignment between the raster origin coordinates of microwave remote sensing imagery and optical remote sensing imagery, or between the digital elevation model (DEM) and optical remote sensing imagery, the microwave remote sensing imagery and the DEM are translated by the entire raster pixel spacing to ensure that their raster origin coordinates are completely consistent with those of the optical remote sensing imagery, achieving complete overlap of the raster origin coordinates of the three types of data. As shown in the figure, the same-grid registration points represent the spatially corresponding raster positions of the data after registration processing. Each raster cell contains corresponding data values in three dimensions: the optical brightness dimension stores the brightness information of the optical remote sensing imagery, the microwave scattering dimension stores the scattering coefficient information of the microwave remote sensing imagery, and the elevation value dimension stores the elevation data of the DEM. This three-dimensional data cube structure provides a unified data foundation for subsequent key technical steps such as heterogeneous noise suppression, sparse thinning operations, and fractal beam mapping, ensuring the effective execution of the Daigu landform automatic identification algorithm under spatial consistency conditions.
[0051] like Figure 3As shown, the local fractal window contains 25 grid cells, each storing corresponding data values in three dimensions. As illustrated, the grayscale distribution of the optical brightness dimension ranges from 136 to 192, exhibiting significant spatial variability. Taking the first row of the window as an example, the grayscale values from left to right are 142, 156, 168, 151, and 139, showing a trend of first increasing and then decreasing. The microwave scattering and elevation dimensions also have their own data distribution characteristics. The data from these three dimensions together constitute the multispectral-topographic integrated features of this local area. The grid cell located at the center of the window is marked as the analysis target, and the surrounding 24 grid cells constitute its neighborhood environment. By statistically analyzing the grayscale distribution of each dimension within the window, the system calculates the step difference frequencies for the optical brightness, microwave scattering, and elevation dimensions. The step difference calculation uses numerical difference between adjacent grid cells, calculating the difference sequence sequentially according to spatial adjacency. For example, the stepped difference sequence for the optical brightness dimension includes values such as 14, -12, 17, -12, 16, -24, 11, -14, 13, and -11, reflecting the frequency and amplitude characteristics of spatial variation in this dimension. Subsequently, the stepped difference frequencies of the three dimensions are concatenated along the row, column, and layer directions to form a one-dimensional fractal difference sequence. This sequence integrates spatial variation information from multi-dimensional data, providing a data foundation for subsequent fractal feature quantification. The length of the one-dimensional fractal difference sequence is related to the number of grid cells and the number of dimensions within the window; the statistical characteristics of the sequence can effectively characterize the complexity and variation patterns of local terrain.
[0052] like Figure 4 As shown, this invention establishes an eight-level energy level classification system based on sequence peak values. The horizontal axis of this coordinate graph represents the peak value of the one-dimensional fractal difference sequence, and the vertical axis represents the corresponding energy level label, ranging from 1 to 8, for a total of eight levels. The eight-level fractal beam energy level thresholds are defined through preset threshold intervals, each interval corresponding to a specific level of geomorphic complexity. The design of the energy level threshold intervals follows the scale invariance principle in fractal theory; lower sequence peak values correspond to simpler geomorphic features and are labeled as lower energy levels; higher sequence peak values correspond to more complex geomorphic features and are labeled as higher energy levels. As shown in the figure, when the sequence peak value is 165, this value falls within the energy level threshold interval 6, therefore the fractal beam energy level of the corresponding grid is labeled as 6. This step-by-step matching mechanism ensures that each sparse candidate grid can obtain a unique and accurate energy level label, providing a quantitative basis for identifying the characteristic regions of the Daigu landform. By combining fractal window analysis with eight-level threshold matching, this invention achieves accurate quantification of terrain complexity, laying a solid technical foundation for subsequent peak ridge line extraction and depth mapping binding operations in elevation profiles.
[0053] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these specific embodiments are merely illustrative. Those skilled in the art can omit, substitute, and modify the details of the above methods and apparatus in various ways without departing from the principles and essence of the present invention. For example, combining the above method steps to perform substantially the same function and achieve substantially the same result according to substantially the same method falls within the scope of the present invention. Therefore, the scope of the present invention is defined only by the appended claims.
Claims
1. An automatic identification method for Daigu landform based on remote sensing and DEM, characterized in that, The method includes: Step 1: Through multi-source remote sensing image acquisition, simultaneously acquire optical remote sensing images, microwave remote sensing images and publicly available digital elevation models (DEMs) covering the target area, and perform grid registration processing on the three types of data to construct a grid data cube of remote sensing DEM data with unified resolution, including three dimensions: optical brightness dimension, microwave scattering dimension and elevation value dimension. Step 2: Perform heterogeneous noise suppression processing on the data of each dimension of the remote sensing DEM grid data cube to obtain the preprocessed grid data cube; Step 3: Based on the preprocessed grid data cube, while maintaining the spatial consistency of the grid, a sparse candidate raster set is generated by sparse thinning. Fractal beam mapping is performed on the sparse candidate raster set to obtain the fractal beam energy level index. Then, the peak ridge line of the elevation profile is deeply mapped and bound to the high energy level region. The ridge feature is extended by combining optical remote sensing image interleaving. When the self-consistent differential change rate meets the termination criterion, the automatic identification result of Daigu landform is output.
2. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 1, characterized in that, Step 1 involves performing grid registration on the three types of data. This includes: projecting the optical remote sensing image, microwave remote sensing image, and digital elevation model (DEM) onto the same geodetic coordinate reference frame and the same map projection plane to ensure that the three types of data have a unified longitude, latitude, and elevation reference; using the spatial resolution of the optical remote sensing image as the reference resolution; if the resolution of the microwave remote sensing image or the DEM is inconsistent with the reference resolution, adjusting its spatial resolution to be consistent with the optical remote sensing image through resampling; after unifying the resolution, calculating the raster origin coordinates of the three types of data; if there is an offset between the raster origin coordinates of the microwave remote sensing image and the optical remote sensing image, or between the DEM and the optical remote sensing image, translating the microwave remote sensing image and the DEM by the whole raster pixel spacing to make their raster origin coordinates completely consistent with the raster origin coordinates of the optical remote sensing image, so that the raster origin coordinates of the three types of data completely overlap.
3. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 2, characterized in that, Step 2, which involves performing heterogeneous noise suppression on each dimension of the remote sensing DEM grid data cube, specifically includes: calculating the noise root mean square error to noise variance ratio for each dimension of the remote sensing DEM grid data in its respective frequency band; if the noise root mean square error of any dimension exceeds 1.2 times the overall mean, then the corresponding raster for that dimension is recorded as a high-noise raster; for the optical brightness dimension and microwave scattering dimension where the high-noise raster is located, an adaptive sliding window averaging method is used for smoothing; subsequently, morphological opening and closing operations are performed sequentially on the elevation value dimension at the same raster location to reduce isolated elevation spike noise; after completion, the noise root mean square error is recalculated, and if the reduction does not reach 30%, the sliding window averaging method is performed again with the same window size; The method for automatic identification of Daigu landform based on remote sensing and DEM as described in claim 3 is characterized in that, in step 3, the process of generating a sparse candidate raster set includes: in the preprocessed grid data cube, first reading the number of rows, columns, and layers; setting the row thinning interval to 3 rasters, the column thinning interval to 3 rasters, and the layer thinning interval to 2 layers; if the number of rows or columns cannot be divided by the interval value, rounding up to maintain the consistency of the grid space; traversing row by row and column by column in the order of 3 rasters in the row direction and 3 rasters in the column direction; when the row index and column index simultaneously meet the interval condition, marking the corresponding raster as a row-column sparse raster; obtaining a row-column sparse raster set after traversal; based on the row-column sparse raster set, selecting upwards at a layer-column interval of 2 layers; when the layer index meets the layer-column interval condition, retaining the corresponding raster, otherwise discarding it; after completion, generating a sparse candidate raster set, and recording the row index, column index, and layer index for each raster.
4. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 4, characterized in that, Step 3: After generating the sparse candidate raster set, the sparse candidate raster set will be processed as follows: In the sparse candidate raster set, at least 5% of the total number of raster pairs will be randomly selected, and their eight-way neighbor raster integrity in the preprocessed grid data cube will be compared; if the number of missing neighbor raster pairs in any raster pair exceeds 1, the overall row or column direction thinning interval will be reduced by 1 raster and sparse thinning will be re-executed until the neighborhood integrity passes the check; for each sparse candidate raster, a cube local fractal window with a side length of 5 rasters will be built with the raster as the center; if the window touches the grid boundary, the window raster will be supplemented in a mirror manner to maintain the integrity of the window shape.
5. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 5, characterized in that, In step 3, the process of performing fractal beam mapping on the sparse candidate grid set to obtain the fractal beam level index includes: within the local fractal window, statistically analyzing the gray level distribution of the optical brightness dimension, microwave scattering dimension, and elevation value dimension respectively; calculating the step difference frequency of the gray level distribution of each dimension; and concatenating the step difference frequencies of the three dimensions in the row direction, column direction, and layer direction to form a one-dimensional fractal difference sequence.
6. Match the one-dimensional fractal difference sequence with the preset eight-level fractal beam energy level thresholds step by step: if the sequence peak falls into the energy level threshold range 1, the corresponding grid beam energy level is marked as 1; if the sequence peak falls into the energy level threshold range 2, the corresponding grid beam energy level is marked as 2; and so on, until energy level 8. After mapping, a unique fractal beam energy level mark is generated for each sparse candidate grid. Write the fractal beam energy level marks of the sparse candidate grids into the index table in ascending order of row index, column index, and layer index. The index table contains three columns of position index and one column of energy level mark, for a total of four columns. After writing, the fractal beam energy level index is obtained.
7. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 6, characterized in that, Step 3, the process of deeply mapping and binding the peak ridges of the elevation profile with high-energy-level regions, includes: performing a differential scan of the digital elevation model (DEM) along the row and column directions, recording grids whose elevation values are greater than the elevation values in the four adjacent directions and whose elevation difference exceeds 3 meters as local peak grids; aggregating all local peak grids according to eight-directional connectivity to obtain continuous ridge fragments; discarding fragments with a length less than 5 grids; otherwise, retaining and outputting the peak ridge set; reading the fractal beam energy level index obtained in the previous operation; if the fractal beam energy level is greater than or equal to 6, marking the corresponding grid as a high-energy-level grid; using the high-energy-level grid as a seed, expanding according to eight-directional connectivity. This process forms high-energy-level regions. If the area of a high-energy-level region is less than 9 grid cells, it is discarded; otherwise, it is retained and the set of high-energy-level regions is output. The center coordinates of each peak ridge segment are calculated. In the set of high-energy-level regions, high-energy-level regions with an Euclidean distance of no more than 5 grid cells from the center coordinates are searched, and the two are paired to generate pre-binding pairs. The absolute value of the elevation difference is calculated for each pre-binding pair. If the absolute value of the elevation difference is no more than 8 meters, the binding is confirmed and the ridge segment-high-energy-level-region binding pair is recorded; otherwise, the binding is abandoned. If a single peak ridge segment and multiple high-energy-level regions meet the binding conditions, the high-energy-level region with the largest area is selected first. After binding is completed, the depth-mapped binding list is output.
8. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 7, characterized in that, Step 3, the process of combining optical remote sensing images to extend ridge features, includes: numbering the optical remote sensing images in ascending order by acquisition date, defining odd-numbered data as the first time-phase sequence and even-numbered data as the second time-phase sequence; freezing the optical brightness dimension data of the second time-phase sequence in the first time-phase sequence, and freezing the optical brightness dimension data of the first time-phase sequence in the second time-phase sequence, alternately synthesizing them into an interleaved spectral sequence; calculating the absolute value of the inter-phase optical brightness difference for each grid in the interleaved spectral sequence; if the absolute value exceeds 20 levels, the grid is marked as a high-gradient grid; otherwise, it is marked as a low-gradient grid. Grid; Read each peak ridge segment in the depth mapping binding list and use its centerline direction as the reference direction; In the interlaced spectral sequence, starting from the endpoint of the peak ridge segment, search for 2 grids on each side of the reference direction; If a grid in the search path meets the following conditions, it is added to the ridge grid set: The grid is a high-gradient grid; The absolute value of the difference between the optical brightness of the grid and the optical brightness of the endpoint grid is between 10 and 30 orders; The absolute value of the fractal beam energy level difference of the grid does not exceed 1; If the search path encounters a high-energy-level region boundary or a low-gradient grid at any point, it immediately stops extending in that direction.
9. The automatic identification method for Daigu landform based on remote sensing and DEM as described in claim 8, characterized in that, Step 3, the process of outputting the automatic identification result of Daigu landform when the self-consistent difference rate of change meets the termination criterion, includes: calculating the connected components of the ridge-type raster set; if the increase rate of the connected component area is less than 5%, the record expansion tends to saturate; when the increase rate of the connected component area after two consecutive extensions is less than 5%, the ridge feature direction consistency extension operation is terminated, and the interwoven extended ridge-type raster set is output; within the coverage area of the interwoven extended ridge-type raster set, the absolute difference of the raster values of the current iteration and the previous iteration is calculated for the optical brightness dimension, microwave scattering dimension, and elevation value dimension, respectively, and the absolute difference of the raster values in the three dimensions is calculated. The arithmetic mean of the differences is taken to obtain the self-consistent difference value of each grid cell. The arithmetic mean of the self-consistent difference values of all grid cells is taken to obtain the average self-consistent difference of all grid cells in this iteration. The absolute difference between the average self-consistent difference of all grid cells in the previous iteration and the average self-consistent difference of all grid cells in the current iteration is the average self-consistent difference rate of change. If the average self-consistent difference rate of change is less than 0.02, it is considered a successful convergence. Otherwise, the interleaved spectral sequence construction operation is resumed and the iteration is restarted. When the average self-consistent difference rate of change is less than 0.02 for two consecutive iterations, the termination criterion is met, and the entire iteration process is stopped.