A method for block-wise multiscale neighborhood terrain analysis of digital elevation models

Abstract

This invention discloses a block-based multi-scale neighborhood terrain analysis method for digital elevation models, comprising the following steps: determining the global maximum neighborhood range; dividing the digital elevation model to be processed into at least one sub-block; reading all sub-blocks in blocks; performing neighborhood terrain calculations within the read sub-block data range, and outputting only the calculation results of the core region as the processing result of the sub-block; and splicing the calculation results of the core regions of all sub-blocks as the global result output. This invention improves the consistency and stability of multi-scale terrain analysis results by uniformly determining the maximum neighborhood range, enabling all neighborhood statistical calculations to be performed under the condition of complete neighborhood data support.

Description

Technical Field

[0001] This invention belongs to the fields of Geographic Information System (GIS), digital terrain analysis and large-scale raster data processing technology. Specifically, it relates to a block-based multi-scale neighborhood terrain analysis method for digital elevation models, which is applicable to neighborhood terrain statistical analysis and terrain feature identification under large-scale DEM data conditions. Background Technology

[0002] With the widespread application of high-resolution digital elevation models (DEMs) in engineering route selection, disaster assessment, and environmental simulation, the objects of terrain analysis have gradually expanded from local areas to large-scale continuous areas, resulting in an exponential increase in data scale. In such scenarios, the extraction of terrain features typically relies on neighborhood statistical operations. These operations are essentially sliding window spatial calculations, which place high demands on data access continuity and neighborhood integrity.

[0003] When the DEM data size is small, traditional whole-data loading methods can directly complete neighborhood statistics and classification calculations. However, under high-resolution, large-scale DEM conditions, existing technologies have the following shortcomings: 1. Lack of a block-based computing mechanism for large-scale DEMs Existing terrain analysis methods mostly rely on neighborhood calculations based on overall raster data, with the default data size allowing for one-time loading and processing. However, as the DEM resolution increases or the analysis range expands, this approach leads to a sharp increase in memory consumption, limiting the widespread application of these algorithms in practical engineering.

[0004] 2. The problem of neighbor information loss during block processing has not been resolved by the system. While some engineering practices employ block-based reading to reduce memory pressure, they fail to consider the dependency of neighborhood computation on block boundary information. This leads to missing neighborhood information in boundary regions, resulting in statistical biases or computational discontinuities. Boundary errors are further amplified, especially under multi-scale neighborhood overlay conditions.

[0005] 3. Lack of a unified control mechanism for multi-scale neighborhood ranges When performing neighborhood calculations at different spatial scales simultaneously, the radii of the neighborhoods at each scale are often different. If there is a lack of a unified neighborhood control range, it is difficult to ensure that all neighborhood calculations have complete data support under block processing conditions, thus affecting the consistency of the calculation results.

[0006] 4. Insufficient efficiency in calculating global statistical parameters Traditional global statistical methods typically rely on traversing the entire map once or repeatedly scanning the data, which is computationally inefficient under large-scale DEM conditions and is not suitable for working in conjunction with block processing mechanisms.

[0007] 5. Lack of scalable computing frameworks for engineering applications Existing methods focus more on the construction of terrain factors and classification rules themselves, while paying insufficient attention to the scalability and stability of the underlying computing framework, making it difficult to support the batch automatic processing needs of high-resolution DEMs in practical engineering. Summary of the Invention

[0008] To overcome the aforementioned shortcomings of existing technologies, this invention addresses the problems of high overall memory pressure, significant boundary errors in block processing, difficulty in unified management of multi-scale neighborhood ranges, and insufficient global statistical efficiency in existing large-scale digital elevation models (DEM) neighborhood statistical calculations. It proposes a block-based neighborhood calculation framework for large-scale DEMs.

[0009] This invention establishes an adaptive maximum neighborhood range determination mechanism and a dynamic four-way boundary expansion mechanism to ensure consistency between neighborhood calculation results and overall calculation results under block-based reading conditions. Combined with an online global statistical update method, it achieves efficient collaborative calculation of multi-scale terrain features under large-scale DEM conditions. Through these technical solutions, stable, continuous, and scalable processing of high-resolution, large-scale DEM data can be achieved under limited computing resources, improving the engineering applicability of terrain analysis and providing a unified and efficient computational framework for applications such as terrain identification, disaster assessment, and engineering site selection.

[0010] To achieve the above objectives, the present invention provides the following technical solution: A block-based multi-scale neighborhood terrain analysis method for digital elevation models, characterized by the following steps: Determine the global maximum neighborhood range based on all neighborhood radii required for multi-scale neighborhood terrain analysis; The digital elevation model to be processed is divided into at least one sub-block. Each sub-block includes a core region and an extended region located outside the core region with a width not greater than the global maximum neighborhood range. Each sub-block is read in blocks, and based on the position of each sub-block in the digital elevation model, the reading range is dynamically determined and expanded outwards with the global maximum neighborhood range as a constraint, so as to read sub-block data including the core area and the complete extended area; Perform neighborhood terrain calculations within the range of the read sub-block data, and output only the calculation results of the core region as the processing result of the sub-block; The calculation results of the core regions of all sub-blocks are concatenated and output as the global result.

[0011] In a preferred embodiment, the global maximum neighborhood range is determined based on the maximum value among all neighborhood radii required for multi-scale neighborhood terrain analysis of the digital elevation model and the pixel radius corresponding to the relative undulation analysis scale.

[0012] In a preferred embodiment, the radius of the global maximum neighborhood range for:

[0013] in: This refers to the maximum neighborhood radius involved in multi-scale neighborhood statistical calculations. The pixel radius corresponds to the relative fluctuation analysis scale.

[0014] In a preferred embodiment, the multi-scale neighborhood radius set includes two neighborhood sets used to describe different spatial dimensions.

[0015] In a preferred embodiment, all sub-blocks are read sequentially one by one by traversal. The digital elevation model is divided into multiple sub-block regions in the row and column directions according to the preset block size. The actual size of the core region of each sub-block is dynamically determined based on the starting position of the current sub-block and its relationship with the data boundary. When the sub-block is located at the data boundary, its actual size is less than or equal to the preset block size.

[0016] In a preferred embodiment, the block-based multi-scale neighborhood terrain analysis method for a digital elevation model is characterized in that the expansion width of the core region in four directions is calculated as follows:

[0017]

[0018]

[0019]

[0020] in: The radius of the global maximum neighborhood. , These are the starting row index and starting column index of the core region of the current sub-block, respectively; , These are the height and width of the core region of the current sub-block, respectively; , These represent the total number of rows and columns in the digital elevation model, respectively. To expand the width upwards, To expand the width downwards, To expand the width to the left, To expand the width to the right.

[0021] In a preferred embodiment, the actual range of the digital elevation model data window that needs to be read is:

[0022]

[0023]

[0024]

[0025] The data window size is

[0026] in: To read the index of the top-left corner of the window, To read the index in the lower right corner of the window.

[0027] In a preferred embodiment, a fast neighborhood calculation method based on cumulative statistical structure is employed, including: Construct the cumulative statistics matrix for the current sub-block; Each pixel in the core region is traversed, and multi-scale neighborhood terrain calculations are performed during the traversal process, and global statistical parameters are updated online. The global statistical parameters include, but are not limited to, the number of global samples, the global mean, and the global sum of squared deviations. Write the calculation results of all pixels within the core area of ​​the current sub-block into the output raster file according to their original spatial positions.

[0028] In a preferred embodiment, the multi-scale neighborhood terrain calculation includes: Calculate the average elevation of the neighborhood surrounding the target pixel at different neighborhood scales; Calculate the multi-scale terrain feature index based on the difference between the elevation of the target pixel and the average elevation of its neighborhood at each scale. A comprehensive topographic feature index is obtained by weighted fusion of topographic feature indices at different scales. In a preferred embodiment, the aforementioned steps are repeated to perform multi-scale neighborhood statistical calculations on all sub-blocks and output the calculation results of the core region of the sub-block. The calculation results output by all sub-blocks are then stitched together to form complete output raster data.

[0029] The technical effects and advantages of this invention are as follows: 1. This application proposes a unified method for determining the global maximum neighborhood range. By uniformly analyzing the statistical radius of the multi-scale neighborhood and the analysis radius of local terrain undulations, the global maximum neighborhood range HALO is determined and used as a unified control scale for block reading and neighborhood expansion. By uniformly determining the maximum neighborhood range, this application enables all neighborhood statistical calculations to be performed under the condition of complete neighborhood data support, thereby improving the consistency and stability of multi-scale terrain analysis results.

[0030] 2. This application constructs a block-based traversal calculation method for large-scale digital elevation model data. By setting sub-block regions of fixed size, the entire DEM data is divided into multiple independently processable data sub-blocks, and read and calculated block by block in row and column order. Compared with the traditional whole data loading method, this application can complete terrain analysis processing without loading all DEM data at once, thereby significantly reducing memory resource consumption, improving the feasibility of the algorithm in ordinary computing environments, and supporting stable processing of large-scale, high-resolution DEM data.

[0031] 3. This application introduces a dynamic four-way boundary expansion strategy during the block reading process, and combines it with the data organization structure of the core area and the expansion area, so that any pixel in the core area of ​​the sub-block can obtain complete neighborhood data support when performing neighborhood statistical calculations. This effectively avoids the problem of missing neighborhood information that is common in traditional block processing methods, and makes the block calculation results consistent with the overall calculation results. At the same time, it avoids statistical errors at the block boundaries, thereby improving the accuracy and stability of large-scale DEM terrain analysis calculations. Attached Figure Description

[0032] Figure 1 This is a flowchart illustrating a block-based multi-scale neighborhood terrain analysis method for a digital elevation model according to this application.

[0033] Figure 2 This is a schematic diagram of the sub-block reading window in this application.

[0034] Figure 3 This is a schematic diagram of the original DEM (partial) in the embodiments of this application.

[0035] Figure 4 This is a schematic diagram of the canyon area extraction results (partial) in an embodiment of this application.

[0036] Figure 5 This is a schematic diagram of the extraction results (partial) of the high mountain area in an embodiment of this application. Detailed Implementation

[0037] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0038] like Figure 1 As shown, a block-based multi-scale neighborhood terrain analysis method based on a digital elevation model (DEM) is presented in this embodiment. The DEM uses DEM raster data as input data. DEM data is a regular grid structure, where each raster cell corresponds to the elevation value of a specific spatial location on the Earth's surface. DEM data typically includes the following basic information: total number of raster rows (H, height), total number of raster columns (W, width), cell spatial resolution (cell_sizecell, representing the actual length of a single pixel on the ground in meters), spatial reference information, and invalid value identifiers. The row and column dimensions are used to determine the spatial extent of the entire raster data and the subsequent traversal index range, while the cell resolution is used for the conversion between the neighborhood scale and the actual geographic distance. Before performing subsequent terrain feature analysis, the aforementioned basic metadata of the DEM is first read to determine the calculation range and related parameters.

[0039] Specifically, the following steps are included: The global maximum neighborhood range is determined based on all neighborhood radii required for multi-scale neighborhood terrain analysis. In other words, the maximum value among all neighborhood radii required for multi-scale neighborhood terrain analysis is taken as the global maximum neighborhood range. In this embodiment, the global maximum neighborhood range is determined based on the maximum value among all neighborhood radius sets required for multi-scale neighborhood terrain analysis of the digital elevation model and the pixel radius corresponding to the relative undulation analysis scale. Specifically, to achieve topographic statistical analysis at different spatial scales, multiple neighborhood radius parameters are pre-defined. The neighborhood radius is represented in pixels and is used to describe the neighborhood range extending outward from a target pixel. Let the neighborhood radius be r, then the corresponding neighborhood window size is:

[0040] That is, the rectangular window area formed by expanding outward by r pixels in both the row and column directions.

[0041] In this embodiment, two sets of neighborhood radius sets are set: one set is used to describe the neighborhood range at a smaller spatial scale, and the other set is used to describe the neighborhood range at a larger spatial scale.

[0042] Let the set of small-scale neighborhood radii be: ; The set of large-scale neighborhood radii is: ; in Represents the radius of the i-th small-scale neighborhood; Let represent the radius of the j-th large-scale neighborhood; n and m are the number of small-scale and large-scale radii, respectively. Different neighborhood radii correspond to different statistical analysis scales across spatial ranges. By performing statistical analysis at multiple neighborhood scales, we can simultaneously characterize local topographic relief features and regional topographic background changes over a larger area.

[0043] When it is necessary to convert the neighborhood radius to the actual ground distance, the spatial resolution of the DEM can be used for the conversion. The actual ground distance D corresponding to the neighborhood can be expressed as:

[0044] Where D represents the actual distance of the neighborhood radius on the ground (in meters); cell_size is the DEM cell spatial resolution.

[0045] In some terrain analysis processes, it is also necessary to statistically analyze local extreme values ​​within a certain spatial range to reflect the degree of terrain undulation in the region. Therefore, a pre-defined actual analysis distance R (in meters) is used to represent the spatial scale employed in local terrain undulation analysis. Since DEM calculations use pixels as the basic unit, this physical distance scale needs to be converted to the corresponding pixel radius. This conversion relationship can be expressed as:

[0046] in The calculated pixel radius, This indicates rounding operations; The corresponding neighborhood window size is:

[0047] This neighborhood range is used to describe the spatial scale involved in subsequent local fluctuation statistics. Through the above conversion, a consistent physical analysis range can be maintained under DEM conditions at different resolutions; To process large-scale DEM data under limited memory conditions, this invention employs a block-based processing approach, dividing the entire DEM data into multiple smaller sub-blocks and processing them block by block. Let the block height be... The block width is The units are all pixels. Each block corresponds to a core region, and its size is:

[0048] in This indicates the number of cells in each sub-block along the row direction; This indicates the number of cells in each sub-block along the column direction. The block size can be flexibly set according to the memory capacity and processing efficiency requirements of the computing device, thereby controlling memory usage while ensuring computing efficiency.

[0049] Through the above steps, this invention completes the basic preparation work for DEM input data and determines the main parameters involved in neighborhood analysis, including multi-scale neighborhood radius, local undulation analysis radius, and block processing size. These parameters provide a unified scale basis for subsequent global neighborhood range determination, DEM block traversal, and neighborhood statistical calculation.

[0050] Under multi-scale neighborhood statistical calculation conditions, different types of terrain analysis processes often correspond to different neighborhood scales. For example, multi-scale neighborhood statistical calculations involve multiple neighborhood windows with different radii, while local terrain undulation analysis corresponds to another set of neighborhood ranges. Therefore, when using a block-based processing method for DEM data calculation, it is necessary to first determine a unified neighborhood influence range to ensure that neighborhood data is not missing when performing neighborhood calculations under block conditions.

[0051] A unified analysis is performed on all radius parameters involved in neighborhood calculation to determine a global maximum neighborhood range, which is used to constrain the data expansion scale during subsequent block reading.

[0052] First, select the maximum value from the multi-scale neighborhood radius set. Let the small-scale neighborhood radius set be... The set of large-scale neighborhood radii is Then, the maximum neighborhood radius involved in multi-scale neighborhood statistical calculations can be expressed as:

[0053] Meanwhile, in the process of local fluctuation analysis, it is also necessary to statistically analyze local extreme values ​​within a certain spatial range, the corresponding pixel radius of which is... .

[0054] Since the two types of calculations mentioned above may involve different neighborhood scales, it is necessary to select a larger value from these two neighborhood radii as the unified neighborhood influence range. Let the global maximum neighborhood radius be... The calculation method is as follows:

[0055] This parameter represents the minimum expansion width that needs to be reserved when reading in blocks, in pixels. This parameter is used to ensure that, under block processing conditions, any pixel in the core area can obtain the required neighborhood data within the currently read data range when performing neighborhood statistical calculations.

[0056] By uniformly determining the maximum neighborhood radius, it can be ensured that each sub-block can provide sufficient neighborhood data support during subsequent DEM block traversal and data reading, thereby ensuring that the block calculation results are consistent with the overall calculation results and avoiding statistical errors caused by neighborhood window out-of-bounds.

[0057] In the subsequent DEM block processing, parameters It will serve as the basic expansion scale when reading sub-blocks, used to determine the neighborhood expansion range in the four directions of the sub-block, and provide complete data support for subsequent neighborhood statistical calculations.

[0058] The digital elevation model to be processed is divided into at least one sub-block. Each sub-block includes a core region and an extended region is set outside the core region. The width of the extended region is no greater than the global maximum neighborhood range.

[0059] All sub-blocks are read in blocks, and based on the position of the sub-blocks in the digital elevation model, the reading range is dynamically determined and expanded in four directions with the global maximum neighborhood range as a constraint, so as to read the sub-block data containing the core area and the complete extended area; After preparing the DEM data and determining the global neighborhood range, this invention employs a block-based traversal method to process the entire DEM data block by block, in order to handle large-scale DEM data under limited computing resources. This method effectively reduces memory usage and improves processing efficiency by dividing the DEM raster data into several fixed-size sub-blocks and reading and calculating them sequentially block by block.

[0060] Let H be the total number of rows and W be the total number of columns in the input DEM data. Based on the block size parameters set in Section 4.1, the DEM data is divided into multiple sub-block regions in both the row and column directions. The size of each sub-block in the row direction is... The dimension in the column direction is The units for all of them are pixels.

[0061] During the traversal, the top-left corner pixel of the sub-block is used as the starting index position of the sub-block. Let the starting row index of the current sub-block be... The starting column index is By moving the starting position of sub-blocks sequentially in both row and column directions with fixed step sizes, a systematic traversal of the entire DEM data can be achieved.

[0062] Since sub-blocks at DEM boundaries may not completely cover the set block size, it is necessary to dynamically calculate their actual size based on their position within the DEM. Let the actual height of the current sub-block in the row direction be... The actual width in the column direction is The calculation method is as follows:

[0063]

[0064] in This is the starting line index of the current sub-block; This is the starting column index for the current sub-block. By using this method, we can ensure that sub-block regions can still be correctly constructed at DEM boundaries without index out-of-bounds issues.

[0065] In this invention, each sub-block corresponds to a core computing region, the size of which is:

[0066] This core region represents the actual processing range of the current sub-block within the DEM data. Subsequent neighborhood statistical calculations and result outputs are all based on this core region, while the area outside the sub-block is only used to provide the data support required for neighborhood calculations.

[0067] By using the above-mentioned block traversal and sub-block construction methods, the entire DEM data can be divided into multiple independently processable sub-block regions without changing the original DEM spatial structure, and a basic data organization structure can be provided for subsequent neighborhood expansion and statistical calculations.

[0068] After completing the DEM block traversal and determining the core region of each sub-block, this invention introduces a boundary expansion strategy during the reading process of each sub-block to ensure the integrity of the neighborhood data when performing neighborhood statistics calculations under block conditions. By reserving a data buffer area of ​​a certain width around the core region, it can be ensured that the neighborhood window of any cell within the core region will not exceed the currently read data range when performing neighborhood statistics.

[0069] This expansion width is determined by the global maximum neighborhood radius as described in the preceding sections. Control. By introducing this extended range during the sub-block reading phase, the neighborhood statistics under block processing conditions can be kept consistent with the overall DEM calculation results.

[0070] Let the index position of the top left corner of the current sub-block core region in the DEM be:

[0071] Its core area dimensions are:

[0072] Since the current sub-block's position in the DEM varies, the spatial extent it can expand in four directions also varies. Therefore, this invention dynamically calculates the expansion width in four directions based on the distance between the sub-block and the DEM boundary.

[0073] Let the expansion widths in the four directions be as follows: (Expand width upwards) (Expand width downwards) (Expand width to the left) (Expand width to the right); The calculation method is as follows:

[0074]

[0075]

[0076]

[0077] The above calculations ensure that when the sub-block is located within the DEM, it can be extended to the maximum neighborhood range HALO in all four directions; and when the sub-block is near the DEM boundary, the extension width will automatically shrink to avoid the reading range exceeding the original DEM data boundary.

[0078] After obtaining the expansion width in four directions, the actual range of the DEM data window to be read can be determined, such as... Figure 2 As shown.

[0079] Let the index of the top left corner of the reading window be... The index in the lower right corner of the read window is The calculation method is as follows:

[0080]

[0081]

[0082]

[0083] Therefore, the size of the reading window can be obtained as follows: This window corresponds to a continuous sub-region in the DEM data. Through the raster data reading interface, the data for this region can be loaded into memory, forming the complete data matrix of the current sub-block. Let the read sub-block data be: , where i and j represent the row and column indices within the sub-block, respectively. This data matrix contains the core computational region and the extended region surrounding the core region, thus providing complete data support for subsequent neighborhood statistical calculations.

[0084] Since the actual data sub-blocks read include extended regions, it is necessary to clearly identify the location of the core computation region within the sub-blocks so that only the core region is output after the neighborhood statistical calculations are completed.

[0085] Let the width of the upper boundary of the extended region in the row direction be... The width of the left boundary in the column direction is Then the index range of the core region in the sub-block data matrix is:

[0086]

[0087] The index range is defined as the core slice within the sub-block. In subsequent calculations, neighborhood statistics are performed across the entire sub-block data, and the final results are extracted and written only within the core region.

[0088] This method ensures that the neighborhood window of each core region cell can be fully acquired in the current sub-block data, while block-based computation does not cause neighborhood out-of-bounds issues, and the computation results between adjacent sub-blocks can maintain spatial continuity.

[0089] After completing the sub-block reading and core region localization, neighborhood statistical calculations and related terrain feature analyses can be performed on the current sub-block data. Once the calculations are complete, only the results corresponding to the core region are written to the output raster file.

[0090] Then, the memory occupied by the current sub-block data is released, and the next sub-block is read and processed in the same way. By using this block-by-block reading and processing method, continuous calculation of large-area DEM data can be achieved without loading the entire DEM data at once.

[0091] This method effectively reduces memory usage and improves the stability and scalability of large-scale DEM data processing while ensuring the integrity of neighborhood statistics.

[0092] Perform neighborhood terrain calculations within the range of the read sub-block data, and only output the calculation results of the core region as the processing result of the sub-block; A fast neighborhood calculation method based on cumulative statistical structure is adopted, and the calculation method includes the following steps: After the sub-block data is read, a two-dimensional cumulative statistical matrix is ​​first constructed for that sub-block data. Let the sub-block elevation data be... Here, i and j represent the row and column indices within the sub-block, respectively. By summing the elevation data along both rows and columns, a cumulative statistical matrix can be obtained. Its definition is:

[0093] This matrix represents the distance from the top left corner of the sub-block to the position. The sum of the elevations of all pixels.

[0094] Using this accumulation structure, the sum of elevations within an arbitrary rectangular neighborhood can be calculated in constant time, thus significantly reducing the complexity of neighborhood statistical computation. The elevation data is for the sub-block, where i and j represent the row index and column index within the sub-block, respectively. Each pixel in the core region is traversed, and multi-scale neighborhood terrain calculations are performed during the traversal process. Global statistical parameters are updated online, including the number of global samples, the global mean, and the global sum of squared deviations. The multi-scale neighborhood terrain calculation includes: Multi-scale neighborhood average elevation calculation; Let the target pixel position be If the neighborhood radius is r, then the corresponding neighborhood window size is:

[0095] In the process of multi-scale neighborhood statistics, it is necessary to calculate the average elevation of the neighborhood around the target pixel at different neighborhood radii.

[0096] By using a cumulative statistical matrix, the sum of elevations within the neighborhood can be quickly obtained, and the average elevation of the neighborhood can be further calculated. This applies to neighborhood radius sets at different scales. and The average elevation within the corresponding neighborhood range is calculated to obtain multi-scale neighborhood statistical results.

[0097] Multi-scale topographic index calculation; After obtaining the average elevation of the neighborhood at each scale, the topographic location features of the target pixel relative to the neighborhood background can be calculated.

[0098] Let the elevation of the target pixel be... The neighborhood average elevation at the k-th neighborhood scale is Then the difference in terrain location at this scale is:

[0099] By calculating the above differences at multiple scales, a multi-scale topographic location index set can be obtained.

[0100] In this invention, to integrate terrain information at different spatial scales, the terrain location indices at each scale are weighted and fused to obtain a comprehensive terrain location index. :

[0101] in This represents the weight coefficient corresponding to the k-th scale. This method allows for the simultaneous characterization of local terrain undulations and regional background terrain structure, thereby improving the stability of terrain recognition.

[0102] Write the calculation results of all pixels within the core area of ​​the current sub-block into the output raster file according to their original spatial positions.

[0103] Repeat the above steps to complete multi-scale neighborhood statistical calculations for all sub-blocks and output the calculation results of the core region of the sub-block. Then, stitch together the calculation results of the core regions of all sub-blocks according to their original spatial locations to generate complete digital elevation model analysis results.

[0104] Example Taking the regional digital elevation model (DEM) data of Anhui Province as an example, this paper describes the block-based multi-scale neighborhood terrain feature calculation method proposed in this invention.

[0105] like Figure 3 As shown, the DEM data used in this embodiment is raster elevation data with a spatial resolution of approximately 12.5 m, in GeoTIFF format, with a raster size of 39538 × 42468 pixels, corresponding to a total number of pixels of 1.68 × 10⁻⁶. 9 Pixels. The uncompressed data size is approximately 6.26 GB. Due to the large data size, using the traditional whole-load method for neighborhood statistical calculation would result in significant memory consumption and reduced computational efficiency. Therefore, the block-based neighborhood calculation method proposed in this invention is used for processing.

[0106] First, the basic information of the DEM data is read, including the total number of raster rows, the total number of raster columns, and the pixel spatial resolution. Then, multi-scale neighborhood analysis parameters are set. In this embodiment, two sets of neighborhood scale parameters are used: small-scale neighborhood radius: 3, 5, 7 pixels; large-scale neighborhood radius: 15, 25 pixels.

[0107] Furthermore, to perform regional terrain undulation analysis, the local terrain analysis distance is set to 300 meters, which is approximately 24 pixels based on the DEM resolution. Considering the aforementioned neighborhood scale parameters, the global maximum neighborhood radius is determined in this embodiment as follows: HALO = 25 pixels This parameter controls the neighborhood expansion range during DEM block reading to ensure that the neighborhood statistics are consistent with the overall calculation results under block calculation conditions. Subsequently, the DEM data is processed through block traversal. In this embodiment, the block size is set as follows: = 512 = 512 That is, the core area of ​​each sub-block is 512 × 512 pixels. By moving the starting position of the sub-blocks with a fixed step size in the row and column directions, a systematic traversal of the entire DEM data can be achieved.

[0108] During each sub-block reading process, the core region is expanded in four directions according to the HALO parameters, thus forming a data sub-block containing a neighborhood buffer. The expansion width in four directions is dynamically calculated based on the DEM boundary position, forming the actual reading window. The read sub-block data includes the core region and its surrounding neighborhood expansion areas, ensuring that any cell within the core region can obtain complete neighborhood data when performing neighborhood statistical calculations.

[0109] After the sub-block data is read, a two-dimensional cumulative statistical matrix is ​​first constructed. Through the cumulative statistical structure, the statistical calculation of any rectangular neighborhood can be completed in constant time, thereby significantly reducing the computational complexity of multi-scale neighborhood operations.

[0110] Subsequently, the average elevation of the neighborhood is calculated at different neighborhood radii, and a multi-scale terrain location index is calculated based on the difference between the target pixel elevation and the average neighborhood elevation. Based on the calculated comprehensive terrain location index, terrain type identification can be performed on the DEM area. In this embodiment, the terrain is divided into two categories: (1) Canyon area like Figure 4 As shown, when the elevation of a target pixel is significantly lower than the average elevation of its neighborhood, i.e., the comprehensive terrain location index is less than a set threshold, the area is identified as a canyon or valley area.

[0111] (2) High mountain area like Figure 5 As shown, when the elevation of a target pixel is significantly higher than the average elevation of its neighborhood, i.e., the comprehensive terrain location index is greater than the set threshold, the area is identified as a mountain or ridge area.

[0112] Using the above method, mountain and canyon structures can be identified simultaneously under multi-scale neighborhood conditions, thereby enabling the automatic extraction of complex terrain structures.

[0113] After completing topographic statistical calculations and classification in the core area of ​​each sub-block, the calculation results are written to the output raster block by block, and the global statistical parameters are updated synchronously. Then, the memory occupied by the current sub-block data is released, and the next sub-block is read and processed in the same way. When all sub-blocks have been processed, the results of the core areas of each sub-block are stitched together according to their original spatial locations to generate a complete DEM topographic feature analysis result raster.

[0114] As can be seen from this embodiment, the block-based multi-scale neighborhood terrain feature calculation method proposed in this invention can achieve calculation of approximately 1.68 × 10⁻⁶ terrain features in a normal computer environment. 9 It enables stable processing of pixel-scale DEM data and can simultaneously identify both mountain and canyon terrain structures within the same computational framework, thereby significantly improving the efficiency and stability of large-scale DEM terrain analysis.

[0115] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0116] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A block-based multi-scale neighborhood terrain analysis method for digital elevation models, characterized in that, Includes the following steps: Determine the global maximum neighborhood range based on all neighborhood radii required for multi-scale neighborhood terrain analysis; The digital elevation model to be processed is divided into at least one sub-block. Each sub-block includes a core region and an extended region located outside the core region with a width not greater than the global maximum neighborhood range. Each sub-block is read in blocks, and based on the position of each sub-block in the digital elevation model, the reading range is dynamically determined and expanded outwards with the global maximum neighborhood range as a constraint, so as to read sub-block data including the core area and the complete extended area; Perform neighborhood terrain calculations within the range of the read sub-block data, and output only the calculation results of the core region as the processing result of the sub-block; The calculation results of the core regions of all sub-blocks are concatenated and output as the global result.

2. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 1, characterized in that, The global maximum neighborhood range is determined based on the maximum value among all neighborhood radii required for multi-scale neighborhood terrain analysis of the digital elevation model and the pixel radius corresponding to the relative undulation analysis scale.

3. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 2, characterized in that, The radius of the global maximum neighborhood range for: in: This refers to the maximum neighborhood radius involved in multi-scale neighborhood statistical calculations. The pixel radius corresponds to the relative fluctuation analysis scale.

4. A block-based multi-scale neighborhood terrain analysis method for a digital elevation model according to claim 2 or 3, characterized in that, The multi-scale neighborhood radius set includes two neighborhood sets used to describe different spatial sizes.

5. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 1, characterized in that, All sub-blocks are read sequentially by traversing through them. The digital elevation model is divided into multiple sub-block regions in the row and column directions according to the preset block size. The actual size of the core region of each sub-block is dynamically determined based on the starting position of the current sub-block and its relationship with the data boundary. When the sub-block is located at the data boundary, its actual size is less than or equal to the preset block size.

6. A block-based multi-scale neighborhood terrain analysis method for a digital elevation model according to claim 1 or 3, characterized in that, The calculation method for the expansion width of the core region in four directions is as follows: in: The radius of the global maximum neighborhood. , These are the starting row index and starting column index of the core region of the current sub-block, respectively; , These are the height and width of the core region of the current sub-block, respectively; , These represent the total number of rows and columns in the digital elevation model, respectively. To expand the width upwards, To expand the width downwards, To expand the width to the left, To expand the width to the right.

7. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 6, characterized in that, The actual range of the digital elevation model data window that needs to be read is: Data window size is in: To read the index of the top-left corner of the window, To read the index in the lower right corner of the window.

8. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 1, characterized in that, A fast neighborhood calculation method based on cumulative statistical structure is adopted, including: Construct the cumulative statistics matrix for the current sub-block; Each pixel in the core region is traversed, and multi-scale neighborhood terrain calculations are performed during the traversal process, and global statistical parameters are updated online. The global statistical parameters include, but are not limited to, the number of global samples, the global mean, and the global sum of squared deviations. Write the calculation results of all pixels within the core area of ​​the current sub-block into the output raster file according to their original spatial positions.

9. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 8, characterized in that, The multi-scale neighborhood terrain calculation includes: Calculate the average elevation of the neighborhood surrounding the target pixel at different neighborhood scales; Calculate the multi-scale terrain feature index based on the difference between the elevation of the target pixel and the average elevation of its neighborhood at each scale. A comprehensive topographic feature index is obtained by weighted fusion of topographic feature indices at different scales.

10. The method for block-based multi-scale neighborhood terrain analysis of a digital elevation model according to claim 1, characterized in that, Repeat the aforementioned steps to perform multi-scale neighborhood statistical calculations on all sub-blocks and output the calculation results of the core region of the sub-block. Then, stitch together the calculation results output by all sub-blocks into complete output raster data.

Images