A multi-scale DEM generation method based on curvature wavelet transform

Abstract

The invention discloses a multi-scale DEM generation method based on curvature wavelet transform, comprising the following steps: obtaining the wavelet basis function and wavelet decomposition layer number of the curvature wavelet transform, as well as the characteristic information amount and output resolution of the DEM related parameters; obtaining the grid surface curvature of the DEM data, decomposing and reconstructing the grid surface curvature of the DEM data based on the wavelet basis function and the wavelet decomposition layer number, and obtaining the reconstruction matrix; based on the reconstruction matrix, extracting the terrain feature points in the DEM data with the characteristic information amount as the threshold; based on the terrain feature points in the DEM data, resampling with the output resolution as the resolution, and realizing DEM generation and synthesis. The invention expresses the terrain features in the DEM data based on the grid surface curvature, and adopts this feature as the basis for selecting the characteristic points of the DEM data, which can effectively identify the terrain and improve the multi-scale synthesis effect of the DEM.

Description

Technical Field

[0001] The invention belongs to the field of information science and technology, and in particular relates to a multi-scale DEM generation method based on curvature wavelet transform. Background Art

[0002] Digital Elevation Model (DEM) is a digital simulation of the actual terrain in the form of a set of ordered numerical arrays. DEM is one of the important data sources for the construction of the national spatial database. It is the terrain basic data for the construction of multi-scale Chinese real-life three-dimensional images. It has multi-scale and multi-level characteristics. In the fields of engineering construction, road planning, military operations, etc., multi-scale and multi-resolution DEM data plays an important role. One of the mainstream ways to construct multi-scale and multi-resolution DEM is to use the DEM synthesis method. The simplicity and efficiency of DEM synthesis contribute to the efficient construction of the national spatial database and promote the development of China's real-life three-dimensional images. It has important theoretical and practical significance.

[0003] At present, DEM integration methods are mainly divided into three categories: methods based on space / frequency filtering, methods based on structure selection, and methods based on the importance evaluation of DEM data points. Among them, the method based on space / frequency filtering mainly decomposes the information in DEM into high information load and low information load parts by filtering, and DEM integration can be achieved by weakening and eliminating the low information load part; the method based on structure selection extracts terrain features by identifying hidden information, and then combines the overall geographic information to evaluate the importance of different terrain features, and finally weakens or eliminates terrain features with lower importance to achieve DEM integration; the method based on the importance evaluation of DEM data points evaluates the importance of data points in DEM by setting certain rules for data of different types and regions, and then reconstructs the low-resolution DEM by selecting DEM data points with higher importance, which are generally the locations of key terrain, to achieve the purpose of DEM integration.

[0004] The method based on spatial / frequency filtering lacks the control of structural information, and its comprehensive results will lose the main terrain features to a certain extent; the method based on structure selection requires complex terrain structure information extraction technology, and the parameters are highly subjective, which makes it difficult to achieve practical purposes; the method based on DEM data point importance evaluation needs to extract important DEM data points first, but there is currently a lack of specific operable algorithms to extract the terrain represented as point elements in DEM. Therefore, it is urgent to propose a multi-scale DEM generation method based on curvature wavelet transform, which is suitable for the synthesis of large-scale DEM data. Summary of the invention

[0005] In order to solve the above technical problems, the present invention proposes a multi-scale DEM generation method based on curvature wavelet transform to solve the problems existing in the above prior art.

[0006] To achieve the above object, the present invention provides a multi-scale DEM generation method based on curvature wavelet transform, comprising the following steps:

[0007] Obtain the wavelet basis function and wavelet decomposition layer number of the curvature wavelet transform, as well as the characteristic information of DEM related parameters and output resolution;

[0008] Acquire the grid surface curvature of the DEM data, decompose and reconstruct the grid surface curvature of the DEM data based on the wavelet basis function and the wavelet decomposition layer number, and obtain a reconstruction matrix;

[0009] Based on the reconstruction matrix, the terrain feature points in the DEM data are extracted with the feature information amount as a threshold;

[0010] Based on the terrain feature points in the DEM data, resampling is performed with the output resolution as the resolution to achieve DEM generation and synthesis.

[0011] Optionally, the process of obtaining the grid surface curvature of DEM data includes: obtaining original DEM data, fitting the elevation values ​​of each pixel in the original DEM data and its eight adjacent pixels, and obtaining the grid surface curvature of the DEM data.

[0012] Optionally, the process of decomposing the grid surface curvature of the DEM data includes: performing wavelet transform on the grid surface curvature of the DEM data based on the wavelet basis function and the wavelet decomposition layer number to obtain 1 low-frequency approximation matrix and 3 high-frequency detail matrices, wherein the high-frequency detail matrix includes a horizontal detail matrix, a vertical detail matrix and a diagonal detail matrix.

[0013] Optionally, the process of reconstructing the raster surface curvature of the DEM data includes: constructing three zero matrices of the same size as the high-frequency detail matrix, using the low-frequency approximation matrix and the three zero matrices to perform wavelet reconstruction to obtain a reconstructed matrix.

[0014] Optionally, the process of extracting terrain feature points in DEM data includes: obtaining positive sequence and negative sequence based on the reconstruction matrix; multiplying the positive sequence and negative sequence by the feature information amount to obtain the corresponding positive sequence feature quantity and negative sequence feature quantity; taking the positive sequence feature quantity and negative sequence feature quantity as thresholds, respectively, to extract elements in the positive sequence and negative sequence as feature point sets; traversing all feature point sets, mapping the row and column numbers of each feature point in the reconstruction matrix to the row and column numbers of the DEM data, and marking the terrain feature points. When the traversal is completed, all terrain feature points in the DEM data are obtained.

[0015] Optionally, the process of obtaining the positive value sequence is: arranging the positive values ​​in the reconstructed matrix from large to small to obtain the positive value sequence;

[0016] The process of obtaining the negative value sequence is as follows: the negative values ​​in the reconstructed matrix are arranged from small to large to obtain the negative value sequence.

[0017] Optionally, the process of realizing DEM generation and integration includes: performing cubic interpolation on terrain feature points in DEM data, and simultaneously resampling with the output resolution as the resolution to obtain a new resolution of DEM data, thereby realizing DEM generation and integration.

[0018] Compared with the prior art, the present invention has the following advantages and technical effects:

[0019] 1. The present invention realizes the expression of terrain features in DEM data based on the curvature of the grid surface. By performing curvature analysis on the DEM data, the features of the DEM data can be quickly and efficiently expressed. This feature is manifested as the horizontal and vertical change trend of the terrain. By using this feature as the basis for selecting feature points of DEM data, the terrain can be effectively identified and the multi-scale comprehensive effect of DEM can be improved.

[0020] 2. The present invention realizes multi-scale expression of terrain in DEM based on curvature wavelet transform, and uses the number of wavelet decomposition layers as different expression scales, which can realize the expression of DEM terrain features from multiple perspectives, and can efficiently realize the multi-scale expression of DEM terrain features, and also provide a high-quality data foundation for multi-scale DEM generation.

[0021] 3. The present invention selects terrain feature points by taking the feature information amount as a threshold, and generates DEM by combining the curvature wavelet transform results at different scales with the cubic interpolation method. It can realize the generation of multi-scale DEM on the basis of taking the feature information amount into consideration, and the generation effect is more in line with the spatial distribution characteristics of DEM data and human spatial cognition. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0023] Figure 1 A schematic diagram of a method flow of an embodiment of the present invention;

[0024] Figure 2 This is a schematic diagram of original DEM data with a resolution of 25 meters according to an embodiment of the present invention;

[0025] Figure 3 A curvature effect diagram of a grid surface calculated according to an embodiment of the present invention;

[0026] Figure 4 Result diagram of curvature wavelet transform when the parameter of the embodiment of the present invention is "bior3.7", 1; wherein (a) is a schematic diagram of a low-frequency approximation matrix, (b) is a schematic diagram of a horizontal detail matrix, (c) is a schematic diagram of a vertical detail matrix, and (d) is a schematic diagram of a diagonal detail matrix;

[0027] Figure 5 This is a diagram showing the effect of reconstructing the grid surface curvature after wavelet transformation according to an embodiment of the present invention;

[0028] Figure 6 This is a schematic diagram of terrain feature points extracted when the feature information amount is 30% according to an embodiment of the present invention;

[0029] Figure 7 This is a schematic diagram of the integrated DEM data with a resolution of 50 meters according to an embodiment of the present invention. DETAILED DESCRIPTION

[0030] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0031] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0032] Embodiment 1

[0033] In this embodiment, a multi-scale DEM synthesis method based on curvature wavelet transform is provided. The method uses grid curvature analysis and image wavelet transform methods to respectively describe the terrain features hidden in the DEM data and extract the terrain feature points at multiple scales, and then uses the feature points as constraints to perform DEM data interpolation. It is essentially a multi-scale DEM synthesis method. First, it is necessary to set the hyperparameters of the curvature wavelet transform and DEM synthesis, and calculate the grid surface curvature of the DEM data; then, based on wavelet transform and reconstruction, the main curvature information hidden in the surface curvature is extracted, and the terrain feature points in the DEM data are extracted based on the feature information threshold; finally, cubic interpolation and resampling are used for heuristic DEM synthesis. This method can effectively improve the efficiency of DEM generation and perform multi-scale data synthesis.

[0034] To achieve the above purpose, this embodiment proposes a multi-scale DEM generation method based on curvature wavelet transform, comprising the following steps:

[0035] Step S1: Set the hyperparameters of curvature wavelet transform and DEM synthesis, including the wavelet basis function α of curvature wavelet transform, the number of wavelet decomposition layers β, and the characteristic information amount γ of DEM related parameters and the output resolution The dimensionless percentage measuring the amount of feature information ranges from 0 to 100%.

[0036] Step S2: Calculate the grid surface curvature of the DEM data; based on the original DEM data, the grid surface curvature of the DEM data is obtained by fitting the elevation values ​​of each pixel in the original DEM data with its eight adjacent pixels.

[0037] Step S3: extract the main information value of the DEM data grid surface curvature, and decompose and reconstruct the DEM surface curvature based on the given wavelet basis function α and wavelet decomposition layer number β;

[0038] In step S3, for the surface curvature C of the DEM data, the main information value in the surface curvature of the DEM data is extracted based on the curvature wavelet transform method, and the specific steps are as follows:

[0039] Step 3.1, take the wavelet basis function α and wavelet decomposition layer number β in step 1, perform wavelet transform on C, and the transformation result is a low-frequency approximation matrix CA and three high-frequency detail matrices, including the horizontal detail matrix CHD, the vertical detail matrix CVD, and the diagonal detail matrix CDD;

[0040] Step 3.2, construct three zero matrices ZM1, ZM2, and ZM3, the size of which is the same as that of the three high-frequency detail matrices, namely, the horizontal detail matrix CHD, the vertical detail matrix CVD, and the diagonal detail matrix CDD;

[0041] In step 3.3, the low-frequency approximate matrix CA in step 3.1 and the three zero matrices ZM1, ZM2, and ZM3 in step 3.2 are used to perform wavelet reconstruction to obtain the reconstruction matrix R, thereby compressing the secondary information in the surface curvature C of the DEM data and highlighting the main information.

[0042] Step S4: extracting terrain feature points in the DEM data based on the feature information amount γ as a threshold;

[0043] In step S4, for the reconstruction matrix R reconstructed by wavelet transform of DEM data, the terrain feature points in the DEM data are extracted based on the feature information amount γ as a threshold, and the specific steps are as follows:

[0044] Step 4.1, statistically analyze the matrix R reconstructed by wavelet transform of DEM data, arrange the positive values ​​in R from large to small, record them as positive value sequence PS, and arrange the negative values ​​from small to large, record them as negative value sequence NS;

[0045] Step 4.2, multiply the feature information amount γ by the number of elements in the positive sequence PS and the negative sequence NS to obtain the number of positive sequence features γPS and the number of negative sequence features γNS;

[0046] Step 4.3, taking the number of positive sequence features γPS and the number of negative sequence features γNS as thresholds, extract the first γPS elements in the positive sequence PS and the first γNS elements in the negative sequence NS as the feature point set FS;

[0047] Step 4.4, traverse the feature point set FS, map its row and column numbers in the wavelet transform reconstruction matrix R to the row and column numbers of the DEM data, and mark the position as a terrain feature point as 1. After the traversal of the feature point set FS is completed, mark the unmarked pixels in the DEM data as 0 to form a terrain feature point binary image FM;

[0048] Step S5: Based on the extracted DEM terrain feature points, a heuristic DEM synthesis is performed using a cubic interpolation method.

[0049] In step S5, a new resolution of the DEM data is obtained by performing cubic interpolation on the terrain feature points in the DEM data and resampling with the output resolution φ as the resolution, which is recorded as Dφ, to achieve DEM integration.

[0050] Embodiment 2

[0051] In this embodiment, the original DEM resolution is 25 meters. Figure 2 As shown; the target DEM resolution β is 50 meters. The basic tools required in this embodiment use the corresponding tools in the GIS software ArcGIS geoprocessing toolbox and the corresponding third-party toolkit in the programming language Python.

[0052] Reference Figure 1 This embodiment provides a multi-scale DEM generation method based on curvature wavelet transform, which specifically includes the following steps:

[0053] Step S1: Set the hyperparameters of curvature wavelet transform and DEM synthesis, including the wavelet basis function of curvature wavelet transform α = "bior3.7", the number of wavelet decomposition layers β = 1, the feature information amount γ = 30% and the DEM output resolution Meters; where it is a dimensionless percentage measuring the amount of feature information, ranging from 0-100%.

[0054] Step S2: Calculate the grid surface curvature of the DEM data;

[0055] Step S3: extract the main information value of the DEM data grid surface curvature, and decompose and reconstruct the DEM surface curvature based on the given wavelet basis function α and wavelet decomposition layer number β;

[0056] Step S4: extracting terrain feature points in the DEM data based on the feature information amount γ as a threshold;

[0057] Step S5: Based on the extracted DEM terrain feature points, a heuristic DEM synthesis is performed using a cubic interpolation method.

[0058] In the implementation of this example, the steps S1 and S2 first set the hyperparameters of the curvature wavelet transform and the DEM synthesis, including the wavelet basis function α="bior3.7" of the curvature wavelet transform, the wavelet decomposition layer number β=1 and the characteristic information amount γ=30%; then the raster surface curvature of the DEM data is calculated, based on the original DEM data, by fitting the elevation values ​​of each pixel in the DEM data and its eight adjacent pixels, the surface curvature of the DEM data is obtained; the raster surface curvature is calculated using the "ArcToolbox\Spatial Analysis\Surface Analysis\Curvature" tool in the GIS software ArcGIS. Figure 3 This is the calculated curvature effect of the grid surface.

[0059] In this embodiment, the step S3 extracts the main information value in the surface curvature C of the DEM data based on the curvature wavelet transform method, and the specific steps are as follows:

[0060] Step 3.1, take the wavelet basis function α = "bior3.7" in step 1 and the wavelet decomposition level β = 1, and perform wavelet transform on C. The transformation result is 1 low-frequency approximation matrix and 3 high-frequency detail matrices, including horizontal detail matrix, vertical detail matrix, and diagonal detail matrix, denoted as CA, CHD, CVD, and CDD; the wavelet transform uses the function in the third-party toolkit "pywt.wavedecn()" in the programming language Python to calculate according to the two input parameters. Figure 4 The curvature wavelet transform result calculated with parameters ("bior3.7",1).

[0061] Step 3.2, construct three zero matrices ZM1, ZM2, and ZM3. The size of the zero matrices is the same as that of the three high-frequency detail matrices. The zero matrices are constructed using the function in the third-party toolkit "numpy.zeros_like()" in the programming language Python.

[0062] In step 3.3, the low-frequency approximate matrix CA in step 3.1 and the three zero matrices ZM1, ZM2, and ZM3 in step 3.2 are used to perform wavelet reconstruction to obtain a reconstruction matrix R, thereby compressing the secondary information in the surface curvature C of the DEM data and highlighting the main information. The wavelet reconstruction uses the function in the third-party toolkit "pywt.waverecn()" in the programming language Python to reconstruct according to the input of a low-frequency approximate matrix CA and three zero matrices ZM1, ZM2, and ZM3. Figure 5 This is the effect of reconstructing the grid surface curvature after wavelet transform.

[0063] In this example, step S4, for the reconstructed matrix R of the DEM data after wavelet transformation, extracts the terrain feature points of the DEM elevation based on the characteristic information amount γ=30% as a threshold, and the specific steps are as follows:

[0064] Step 4.1, perform statistical analysis on the reconstructed matrix R of the DEM data after wavelet transformation, arrange the positive values ​​in R from large to small, recorded as the positive value sequence PS, and arrange the negative values ​​from small to large, recorded as the negative value sequence NS; the statistical analysis operation is performed using the function in the third-party toolkit "numpy.sort()" in the programming language Python.

[0065] Step 4.2: multiply the feature information γ by the number of elements in the positive sequence PS and the negative sequence NS to obtain the number of positive sequence features γ PS and the number of negative sequence features γ NS ; The number of positive and negative sequence features is calculated using the “*” operator in the programming language Python.

[0066] Step 4.3, respectively, the number of positive sequence features γ PS and the number of negative sequence features γ NS As the threshold, extract the first γ in the positive sequence PS PS elements and negative value sequence NS in the first γ NS elements as the feature point set FS; the feature point set is calculated using the index of the "list" object and the "*" operator in the programming language Python.

[0067] Step 4.4, traverse the feature point set FS, map its row and column numbers in the wavelet transform reconstruction matrix R to the row and column numbers of the DEM data, and mark the position as a terrain feature point as 1. After the traversal of the feature point set FS is completed, mark the unmarked pixels in the DEM data as 0 to form a terrain feature point binary image FM; wherein, the terrain feature point matrix is ​​constructed using comparison operators such as “|, >, <” in the programming language Python, and the terrain feature point binary image is generated using the function in the “PIL.save()” third-party toolkit in the programming language Python. Figure 6 This is the effect picture after selecting feature points.

[0068] In this embodiment, step S5 is to perform cubic interpolation on the terrain feature points in the DEM data and to Resample to resolution, generating resolution The DEM data is recorded as DEM synthesis is achieved. The cubic interpolation is performed using the function in the third-party toolkit "scipy.interpolate.griddata()" in the programming language Python, and the resampling is calculated using the "ArcToolbox\Data Management Tools\Raster\Raster Processing\Resampling" tool in the GIS software ArcGIS. Figure 7 This is the DEM data with a resolution of 50 meters after integration.

[0069] The above are only preferred specific implementations of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A multi-scale DEM generation method based on curvature wavelet transform, characterized in that: The following steps are involved: Obtaining wavelet basis functions and wavelet decomposition layer numbers of curvature wavelet transform, as well as DEM related parameters, the parameters including feature information amount and output resolution; Acquire the grid surface curvature of the DEM data, decompose and reconstruct the grid surface curvature of the DEM data based on the wavelet basis function and the wavelet decomposition layer number, and obtain a reconstruction matrix; Based on the reconstruction matrix, the terrain feature points in the DEM data are extracted with the feature information amount as a threshold; Based on the terrain feature points in the DEM data, resampling is performed at the output resolution to achieve DEM generation and synthesis; The process of extracting terrain feature points in DEM data includes: obtaining positive sequence and negative sequence based on reconstruction matrix; multiplying the feature information by positive sequence and negative sequence respectively to obtain the corresponding positive sequence feature quantity and negative sequence feature quantity; taking the positive sequence feature quantity and negative sequence feature quantity as thresholds respectively, extracting elements in positive sequence and negative sequence as feature point sets; traversing all feature point sets, mapping the row and column number of each feature point in the reconstruction matrix to the row and column number of DEM data, and marking the terrain feature point. After traversal is completed, all terrain feature points in DEM data are obtained; The process of obtaining the positive value sequence is as follows: the positive values ​​in the reconstructed matrix are arranged from large to small to obtain the positive value sequence; The process of obtaining the negative value sequence is as follows: the negative values ​​in the reconstructed matrix are arranged from small to large to obtain the negative value sequence.

2. The multi-scale DEM generation method based on curvature wavelet transform according to claim 1 is characterized in that: The process of obtaining the raster surface curvature of DEM data includes: obtaining the original DEM data, fitting the elevation values ​​of each pixel in the original DEM data and its eight adjacent pixels, and obtaining the raster surface curvature of the DEM data.

3. The multi-scale DEM generation method based on curvature wavelet transform according to claim 1 is characterized in that: The process of decomposing the grid surface curvature of DEM data includes: performing wavelet transform on the grid surface curvature of DEM data based on the wavelet basis function and the wavelet decomposition layer number, obtaining a low-frequency approximation matrix and three high-frequency detail matrices, wherein the high-frequency detail matrix includes a horizontal detail matrix, a vertical detail matrix and a diagonal detail matrix.

4. The multi-scale DEM generation method based on curvature wavelet transform according to claim 3 is characterized in that: The process of reconstructing the raster surface curvature of DEM data includes: constructing three zero matrices with the same size as the high-frequency detail matrix, using the low-frequency approximate matrix and the three zero matrices for wavelet reconstruction, and obtaining the reconstructed matrix.

5. The multi-scale DEM generation method based on curvature wavelet transform according to claim 1 is characterized in that: The process of realizing DEM generation and synthesis includes: performing cubic interpolation on the terrain feature points in the DEM data, and resampling at the output resolution to obtain a new resolution of the DEM data, thereby realizing DEM generation and synthesis.

Images