A gangping depression terrain evaluation method based on multi-scale DEM fusion verification

Abstract

A kind of gao ping depression topography evaluation method based on multiscale DEM fusion verification.Method for automatically identifying and evaluating micro-terrain relates to the technical field, specifically relates to a kind of gao ping depression topography evaluation method based on multiscale DEM fusion verification.The method comprises the following steps: obtaining target county, measured DEM data is denoted as, the DEM data with precision of 12.5m is denoted as, three-level correction is carried out to the DEM data with 12.5m, weighted fusion is generated based on three-level correction result and, to generate fusion DEM;Error quality index is calculated, and topographic core factor and auxiliary factor are calculated based on fusion DEM;Factor confidence index is set, to determine whether topographic core factor and auxiliary factor meet confidence index;Gao ping depression decision tree is constructed, and the topography of target county is classified;Classification confidence is set, to determine whether classification result meets classification confidence.

Description

Technical Field

[0001] This invention relates to the field of automated identification and evaluation technology of micro-topography, specifically to a method for evaluating hilly, flat, and depression terrain based on multi-scale DEM fusion verification. Background Technology

[0002] Farmland topography assessment typically relies on DEM data. Topographic factors such as elevation, slope, aspect, relative elevation difference, and runoff characteristics are used to determine the relative position of plots within the regional topography, thus providing a basis for farmland quality assessment, farmland water and fertilizer management, drainage and irrigation, and land consolidation. In practical applications, available DEM data can be broadly categorized into two types: one is a national or regional scale DEM with relatively complete coverage, which is easy to acquire and suitable for continuous analysis at the county level and above; the other is a high-precision DEM generated through UAV aerial surveys, field measurements, or local mapping, which can reflect local plot undulations in greater detail, but has higher acquisition and processing costs and usually limited coverage. For the evaluation of hills, flatlands, and depressions at the farmland plot scale, using only low-resolution DEMs easily loses micro-topographical differences within and between plots, while relying solely on high-precision DEMs is insufficient to meet the cost and data coverage requirements for large-scale continuous applications.

[0003] Existing technologies include research on DEM error evaluation and correction. For example, the Chinese invention patent "A Pixel-Scale-Based Method for DEM Data Error Evaluation and Correction" (CN108038086A) uses original DEM data and elevation point data to obtain error point data. This error point data is then divided into modeling points and check points. A regression equation is established between the elevation error values ​​of the modeling points and the surface morphology feature values ​​to obtain the pixel-scale DEM error distribution. The check points are then used to evaluate the correction effect. The key to this method is using the statistical relationship between elevation point errors and surface morphology features to estimate and correct the accuracy of the DEM error, which can improve the overall accuracy of the original DEM. However, in farmland topography evaluation scenarios, the error of low-resolution DEMs is not only manifested as a uniform overall deviation but is also affected by landform type, sampling scale, local undulations, plot boundaries, and regional spatial location. If only a general pixel error regression method is used, the subsequent calculation of factors such as slope, runoff, and humidity index may still be affected by the propagation of local errors, thus affecting the classification results of hills, plains, and depressions.

[0004] In the area of ​​micro-topography classification, the Chinese invention patent application "A Method for Micro-Topography Classification in Semi-arid Regions Based on High-Precision Digital Elevation Model" (CN107330422A) uses UAVs to acquire high-precision images and extract high-precision DEMs. It further calculates factors such as slope aspect, slope gradient, inverse terrain DEM, positive and negative terrain distribution, and elevation, and then classifies the terrain into types such as ridges, shady slopes, sunny slopes, sandy plains, and sand dunes according to predefined classification rules. This method demonstrates that high-precision DEMs can be used for small-scale micro-topography identification, and that factors such as slope gradient, slope aspect, and positive and negative terrain distribution play a role in micro-topography classification. However, such methods typically rely on high-precision DEMs of specific regions and preset thresholds, and the classification objects are often targeted at specific landform types in semi-arid regions. When the evaluation object shifts to hills, plains, and depressions within county-level farmland, the plot scale is finer, the terrain undulation is relatively smaller, and the classification results are more sensitive to DEM errors. Directly applying fixed thresholds or specific landform classification rules can easily lead to inconsistent classification of the same type of terrain in different regions.

[0005] In summary, existing methods often treat DEM error correction, topographic factor calculation, and terrain classification as independent steps. Whether local errors still exist after DEM correction typically does not further affect the setting of classification thresholds or rules. Furthermore, when abnormal patches appear in the terrain classification results or do not match the actual landforms, there is a lack of mechanisms to reverse-check DEM errors or topographic factor calculation results. For county-level cultivated land evaluation, the characteristics of DEM errors and their impact on factors such as slope, runoff, and humidity index are not entirely the same across different townships, different geomorphic units, and different cultivated land types. Using uniform correction parameters and classification thresholds can easily lead to under-correction in some areas and over-correction in others, ultimately affecting the stability of the terrain evaluation results. Summary of the Invention

[0006] To address the aforementioned problems, the purpose of this invention is to propose a method for evaluating hilly, flat, and depression terrain based on multi-scale DEM fusion verification. By constructing a multi-scale DEM fusion verification, error type decomposition, and spatial heterogeneity adaptive correction mechanism, a fused DEM that balances coverage and local accuracy is generated. Furthermore, based on EQI adaptive thresholding and decision trees, automated quantitative classification of hills, flat areas, and depressions at the plot scale is achieved, thereby improving the accuracy, stability, and cross-regional applicability of farmland terrain evaluation.

[0007] The method includes the following steps: S1. Obtain the DEM data of the target county with an accuracy of 0.01m, denoted as S1. DEM data with an accuracy of 12.5m is denoted as and will Preprocessing is performed to obtain the same as Datasets with completely consistent spatial extent ; S2, targeting Perform three-level correction, based on the results of the three-level correction and Perform weighted fusion to generate a fused DEM; S3, Calculate the error quality index And based on the fused DEM, the core and auxiliary factors of the terrain are calculated; Set the factor credibility index. When the core topographic factor and auxiliary factor meet the credibility index, proceed to step S4. Otherwise, return to step S2; S4, based on Calculate the adaptive relative elevation threshold and utilize A decision tree for hilly, flat, and depression terrain was constructed using core and auxiliary topographic factors to classify the terrain of the target county. S5. Set the classification confidence level. When the classification result meets the classification confidence level, output the classification result. Otherwise, return to step S2.

[0008] Furthermore, the preprocessing includes: unifying resolution, cropping, and outlier removal.

[0009] Furthermore, the three-level correction includes: Township Correction: ,in, The spatial location of the point to be corrected. This represents the average error of the township to which the current point to be corrected belongs; Village-level correction: , This represents the average error of the village to which the current point to be corrected belongs; Spatial correlation error correction: , This represents the spatial correlation error value estimated by ordinary Kriging interpolation for the current point to be corrected.

[0010] Furthermore, the calculation formula for fused DEM is as follows: ,in, and These are the corresponding weighting coefficients.

[0011] Furthermore, the error quality index The formula for calculation is: ,in, The RMSE value is the 3×3 neighborhood window of the current spatial location of the point to be corrected. and They are respectively The maximum and minimum values.

[0012] Furthermore, the core terrain factors include: relative elevation. and slope ; The auxiliary factors include: plane curvature Pc, topographic humidity index TWI, and confluence accumulation TCA; The factor credibility index is set based on the coefficients of variation of the core terrain factor and auxiliary factors.

[0013] Furthermore, adaptive relative elevation threshold ,in, As the basic threshold for county-level elevation, Adjustment coefficient for landform type, This is the error quality adjustment factor. .

[0014] Furthermore, the decision tree uses slope as the root node and divides the target plot into a gentle terrain branch and a non-gentle terrain branch based on the comparison result between the slope and the preset slope threshold. In the gentle terrain branch and the non-gentle terrain branch, the comparison results of relative elevation and adaptive relative elevation threshold are used as secondary judgment conditions, and the secondary judgment results are verified by combining at least one auxiliary factor among plane curvature, topographic humidity index and runoff accumulation, so as to determine the classification result of the target county.

[0015] Furthermore, when the secondary determination result is inconsistent with the verification result obtained from the auxiliary factor, conflict arbitration is performed; The conflict arbitration determines the classification result of the target county according to the rule that core factors take precedence over auxiliary factors.

[0016] Furthermore, the formula for calculating the factor credibility index is: ,in, , and These are the weighting coefficients. For slope confidence level, The confidence level is the relative elevation. This represents the average confidence level of the auxiliary factors.

[0017] The beneficial effects of the method described in this invention are as follows: (1) Before using the 12.5m DEM for plot-scale analysis, this invention first introduces measured elevation data for verification, and calculates and corrects systematic deviations separately for spatial units such as townships and villages. After this processing, the DEM correction no longer uses a uniform correction amount for the whole county, but corrects it separately according to the actual error situation of each area, which can reduce the elevation deviation of the original 12.5m DEM and provide a more reliable data basis for the subsequent calculation of topographic factors such as slope, relative elevation, and runoff accumulation.

[0018] (2) This invention does not simply regard DEM error as a single overall error. Instead, it uses Moran's I, township-level error statistics, and experimental variogram analysis to determine whether the error has spatial clustering and different sources, and classifies it into systematic bias, spatially correlated error, and random error. For the above errors, zonal bias elimination, Kriging spatial interpolation, and weighted fusion are used respectively, which can avoid the problem that a single regression correction method is insufficient in handling different error sources, and make the correction results more consistent with the actual error distribution of different terrain areas within the county.

[0019] (3) This invention sets quality judgment indicators such as EQI, FCI, and CC between DEM correction, terrain factor calculation, and hill-flat-slope classification. When the quality of the preceding DEM is insufficient, the terrain factor fluctuates greatly, or the classification confidence is low, the system can adjust the fusion parameters, neighborhood window, or classification threshold according to the corresponding indicators, instead of directly outputting low-confidence results. This can reduce repeated manual checks and adjustments, and improve the stability and automation of the entire evaluation process.

[0020] (4) In determining the ridge, flat, and depression characteristics of a land parcel, this invention not only uses a single indicator such as elevation or slope, but also simultaneously introduces relative elevation, slope, plane curvature, topographic humidity index, and runoff accumulation. A decision tree is constructed by first determining the core factors and then verifying the auxiliary factors. Relative elevation and slope are used to determine the basic undulation and tilt of the land parcel, plane curvature is used to determine the local convex and concave morphology, and topographic humidity index and runoff accumulation are used to reflect the tendency of water accumulation and drainage. Therefore, it can more completely describe the ridge, flat, and depression characteristics of the land parcel, improving the accuracy and interpretability of the classification results.

[0021] (5) The relative elevation threshold and RMSE judgment threshold of this invention are not fixed, but can be adjusted according to the landform type, terrain complexity and DEM error quality. For different county farmland scenarios such as plains, hills and mountains, the relevant thresholds can be adjusted to achieve adaptation, so that the method is not only applicable to a single test area, but also easy to promote to practical applications such as high-standard farmland construction, cultivated land quality evaluation, drainage and irrigation and precision agricultural management. Attached Figure Description

[0022] Figure 1This is a schematic diagram of the classification and determination process for Gangpingwa as described in this invention; Figure 2 This is a schematic diagram of the Gangpingwa classification confusion matrix described in this invention; Figure 3 This invention provides a comparison of the spatial heterogeneity of RMSE in various townships. Figure 4 This is a schematic diagram of the D1-D2 elevation scatter plot and linear regression analysis described in this invention; Figure 5 This is a schematic diagram of the spatial distribution of the system deviations in each township as described in this invention; Figure 6 This is a schematic diagram of the cumulative distribution curve of absolute error described in this invention. Detailed Implementation

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

[0024] Example 1 This embodiment provides a method for evaluating the terrain of Gangpingwa based on multi-scale DEM fusion verification. The method includes the following steps: S1. Obtain the DEM data of the target county with an accuracy of 0.01m, denoted as S1. DEM data with an accuracy of 12.5m is denoted as and will Preprocessing is performed to obtain the same as Datasets with completely consistent spatial extent ; The relevant operations in step S1 will be introduced with specific examples: Data Acquisition: In this embodiment, a high-precision DEM of Dongfeng County, Liaoyuan City, Jilin Province was acquired through field measurements, denoted as D1 (measured points n=104,265, covering 12 townships and 21 villages, with an accuracy of 0.01m). Based on the geospatial data cloud platform, a 12.5m DEM of the whole country was collected, denoted as D2, along with county-level land parcel vector data and county-level administrative boundary data.

[0025] Unified benchmarks: The coordinate system is CGCS2000, and the elevation benchmark is the 1985 National Elevation Benchmark; bilinear interpolation is used to resample D2 to the resolution of the measured D1 to ensure the consistency of the data spatial scale.

[0026] Range clipping and outlier handling: Using the administrative boundary of Dongfeng County as a mask, the resampled D2 was clipped to obtain a dataset D2' with the same spatial range as D1; ​​the three-standard-deviation method was also used to remove outliers such as null values, negative values, and outliers in the two types of DEMs (a total of 485 outliers were removed), while retaining valid data near the topographic structure lines.

[0027] S2, targeting Perform three-level correction, based on the results of the three-level correction and Perform weighted fusion to generate a fused DEM; The relevant operations in step S2 will be introduced with specific examples: Firstly, regarding and Multi-scale DEM fusion verification, error type decomposition, and spatial heterogeneity analysis were performed. The analysis process is as follows: (1) Sampling and error calculation: Taking Dongfeng County as the study area, 104,265 calibration points were evenly distributed, and the measured elevations of each point were extracted. ( ) and 12.5m DEM elevation ( ), The index value of the verification point. Calculate the following error index: Basic error index: Calculate the absolute error of each verification point ( ), Root Mean Square Error (RMSE), Mean Error (ME), Mean Absolute Error (MAE), Maximum Error and Error Threshold (e.g.) (The proportions); Calculated, the overall RMSE of Dongfeng County is 2.52m, ME is -0.095m, and MAE is 1.92m. It accounts for 50.87%. It accounts for 19.94%.

[0028] Linear Regression Analysis: A linear regression model D2' = a * D1 + b was established between D1 and D2', and the Pearson correlation coefficient and coefficient of determination R² were calculated. Experimental results show that the regression equation is D2' = 1.008753 * D1 + (-3.2046), Pearson r = 0.9957, and R² = 0.9915, indicating a very strong linear correlation between the two. However, the slope deviates from 1.0 and the intercept is not zero, indicating systematic proportional bias and offset bias.

[0029] Spatial structure index autocorrelation analysis: Moran's I index was used to quantify the spatial autocorrelation of DEM error. The calculation formula is as follows: ,in , The index value of the checkpoint, and Together they form a spatial point pair It is used to calculate the spatial correlation of errors between two verification points. For the first The verification point and the first The spatial weights between the verification points are used to represent the spatial correlation strength between the two points; in this embodiment, if the distance between the two points is... ,but ,like ,but , The sum of all spatial weights. , For the first The deviation between the point error value and the global average error. This indicates that the error at this point is higher than the average level. This indicates below average. The error sequence consists of the DEM error values ​​of all verification points. , This represents the mean of all verification points. After calculation, This indicates that the error exhibits moderate positive spatial autocorrelation and spatial clustering characteristics, necessitating spatial heterogeneity analysis and regional correction.

[0030] The spatial heterogeneity analysis of the error is as follows: Taking the 12 townships of Dongfeng County as the basic spatial units, the ME (systematic bias), RMSE, MAE, and error coefficient of variation (CV) of each township are calculated to analyze the spatial heterogeneity characteristics of the error. The verification results show that the RMSE range of each township is 1.64m (Yunding Town) to 2.91m (Liaoheyuan Town), the ME range is -1.47m (Sanhe Manchu and Korean Township) to +0.97m (Dongfeng County), and the error coefficient of variation (CV) ranges from 0.81 to 238.95, indicating that there is significant spatial heterogeneity in the DEM error, and a regional adaptive correction strategy is required.

[0031] Error type decomposition: Using the geostatistical experimental variogram analysis method, the DEM error is decomposed into three statistically independent components: Systematic bias variance ( ): Reflects the systematic differences in the mean error among different administrative spatial units (townships, villages), estimated by the variance of ME for each township. The measured results for Dongfeng County are as follows: This accounts for 8.3% of the total error variance.

[0032] Spatial correlation error variance ( ): Reflects the structural variation of error in space, estimated by the difference between the sill value and the nugget value of the variogram. The measured results in Dongfeng County are as follows: It accounts for 46.0% of the total error variance, making it the largest component of the error.

[0033] Random error variance ( (This reflects irregular, spatially uncorrelated random variation, estimated by the nugget value of the variation function. The measured results in Dongfeng County are...) This accounts for 42.0% of the total error variance.

[0034] The sum of the variances of the three components is 6.1066 m², which is 96.4% consistent with the measured total variance of 6.3371 m², verifying the rationality of the error type decomposition.

[0035] Error judgment criteria: Based on relevant literature, national standard (CH / T 9009-2010), and actual measurement results in Dongfeng County, multi-level error acceptance thresholds are set as shown in Table 1: Table 1

[0036] Note: The above threshold has been calibrated based on the verification results of 104,265 measured samples in Dongfeng County (overall RMSE=2.52m). In practical applications, it can be adaptively adjusted within ±20% according to the complexity of the county's terrain.

[0037] Based on the spatial heterogeneity analysis results above, it can be seen that D2' not only has an overall elevation deviation compared to D1, but also exhibits systematic differences with variations in township spatial units and residual errors with spatial clustering characteristics. If D2' is directly used for plot-scale terrain evaluation, the error will be further propagated to the calculation results of topographic factors such as slope, relative elevation, and runoff accumulation, affecting the accuracy of hill-flat-depression classification. Therefore, adaptive correction is needed in subsequent steps to improve the local elevation accuracy and reliability of terrain evaluation of the fused DEM.

[0038] This embodiment employs a correction strategy of "systematic bias elimination + spatial correlation error correction + random error weighted fusion": The first is a two-level hierarchical offset method consisting of "township-level offset + village-level fine correction". Township-level correction: The ME value for each township is calculated as the system bias, and township-by-town offset correction is performed on D2'. ,in, The spatial location of the point to be corrected. The average error of the township to which the current point to be corrected belongs is used to characterize the overall systematic offset of the 12.5m DEM relative to the measured elevation within that township. This correction eliminates the systematic offset of each township, making the ME of each township approach 0. After verification, the overall RMSE is reduced by approximately 1.8% after this correction.

[0039] Village-level correction: Based on the township-level correction, the remaining system bias is calculated using villages as spatial units, and a secondary offset correction is performed. , The average error of the village to which the current point to be corrected belongs is used to characterize the overall systematic shift of the 12.5m DEM within the village relative to the measured elevation. After verification, the overall RMSE was reduced by approximately 2.6% after the two-level correction.

[0040] Furthermore, based on the above two-level corrections, spatial correlation error correction is performed: Based on the experimental variograms of each township and the fitted spherical model parameters (null value, sill value, range), a spatial autocorrelation model is constructed.

[0041] Using ordinary kriging interpolation, with spatial coordinates As input (without using landform features), based on spatial autocorrelation structure pairs The spatial correlation error components in the data are subjected to optimal unbiased estimation to obtain spatially continuous spatial correlation error values, specifically: Verification point error after village-level correction Using the sample as an example, a spatial autocorrelation matrix is ​​constructed based on the fitted spherical variogram, and the optimal weights are obtained by solving the system of linear equations. For the current pixel position to be corrected The neighborhood modeling point errors are weighted and summed. This will give you the estimated spatial correlation error for that pixel.

[0042] Subtract the spatial correlation error of the Kriging estimate from D2': , This represents the spatial correlation error value estimated by ordinary Kriging interpolation for the current point to be corrected.

[0043] Furthermore, based on spatial correlation error correction, random error weighted fusion is performed: A weighted fusion strategy is used to generate a fused DEM by combining the corrected D2' with the measured D1. ); The formula for calculating the merged DEM is: ,in, and The corresponding weighting coefficients are dynamically allocated based on the magnitude of local errors. , , Current point to be calibrated The RMSE value of the 3×3 neighborhood window. In the evaluation of the correction effect, the RMSE of the fused DEM was reduced by 40%~60% compared with the original 12.5m DEM. for The maximum value.

[0044] S3, Calculate the error quality index And based on the fused DEM, the core and auxiliary factors of the terrain are calculated; Set the factor credibility index. When the core topographic factor and auxiliary factor meet the credibility index, proceed to step S4. Otherwise, return to step S2; The relevant operations in step S3 will be introduced with specific examples: Error Quality Index The formula for calculation is: ,in, The RMSE value is the 3×3 neighborhood window of the current spatial location of the point to be corrected. and They are respectively The maximum and minimum values. The value ranges from 0 to 1. A higher EQI indicates a better quality integrated DEM for the township.

[0045] Based on fused DEM ( Combining classic international methods for topographic analysis with domestic research findings, the following core topographic factors and auxiliary factors were calculated: Core factors: relative elevation The median elevation of the 100m×100m neighborhood (k=8~15 nearest neighbors) of the plot to be classified. Using this as a benchmark, the relative elevation of each point is calculated. This index eliminates the influence of the overall topographic relief of the region on the absolute elevation, and directly reflects the local topographic relief characteristics; The elevation of each pixel is extracted directly from the D-fusion raster elevation values, and the relative elevation is obtained by comparing it with the median elevation of the neighborhood. The calculation formula is as follows: ,in for Within the surrounding area The median elevation (m) of the nearest neighbor pixels.

[0046] slope :pass The elevation values ​​of the 3×3 neighborhood window raster were calculated using a third-order unweighted difference algorithm (Liu Xuejun et al. confirmed that this algorithm has the best slope calculation accuracy when there is a certain error in the DEM). The calculation formula is as follows: ,in,

[0047]

[0048] in, Only The elevation values ​​are for a 3x3 window grid. For grid resolution; the 3x3 windows are numbered as follows: Z1Z2Z3: Previous line (Northwest North Northeast); Z4Z5Z6: Current row (West center East); Z7Z8Z9: Next line (Southwest, Southeast, Southeast); Z5= : Elevation value of the center pixel; g: The raster resolution (m); cofactor: Plane curvature Pc: Plane curvature through The elevation values ​​of the 3×3 neighborhood window raster are calculated using second-order partial derivatives and used to determine the convexity or concavity of the terrain; convexity is positive (uplift characteristics), and concaveness is negative (depression characteristics). The calculation formula is as follows:

[0049] in: for The second partial derivative of the elevation surface along the x-direction is used to characterize the degree of curvature of the elevation in the east-west direction; if A larger absolute value indicates that the terrain exhibits more pronounced variations in elevation along the east-west direction.

[0050]

[0051] for The second partial derivative of the elevation surface along the y-direction is used to characterize the degree of curvature of the elevation in the north-south direction; if A larger absolute value indicates that the terrain exhibits more pronounced variations in elevation along the north-south direction.

[0052]

[0053] for The mixed second-order partial derivatives of the elevation surface with respect to the x and y directions are used to characterize the coupling relationship between east-west and north-south elevation changes, that is, the torsional changes of the terrain in the diagonal or oblique directions.

[0054]

[0055] Topographic Wetness Index (TWI): TWI does not directly from Instead of direct calculation, it is obtained indirectly through Slope and TCA, using the classic formula as follows: Where As represents a specific catchment area, and a high TWI value corresponds to a depression that is prone to water accumulation. in: AS: Specific catchment area (m² / m), AS = TCA / L; TCA: Convergence Accumulation (by...) (obtained through flow direction analysis calculations) L: Contour line unit width (taken as raster resolution g); Slope (in radians) .

[0056] Total Contributing Area (TCA): TCA through Global slope analysis determines the water flow direction, and then the number of upstream pixels is accumulated. The multi-flow direction algorithm (FD8) is used for calculation. A high TCA value reflects strong regional water collection capacity, providing support for depression identification. The calculation formula is as follows:

[0057] in: 1: The contribution of this pixel itself Flow direction All upstream neighboring pixel coordinates; Upstream pixels are allocated to The proportion of traffic; For the first The pixel number of the upstream water inflow (i.e., all slopes pointing to the current pixel). The first neighboring pixel indivual); The number representing one of the eight neighborhood directions ( These correspond to the eight directions: east, southeast, south, southwest, west, northwest, north, and northeast, respectively. First, iterate through the eight neighboring directions of the current pixel (…). Specifically, the relationship involves determining which directions the slope of neighboring pixels points towards the current pixel (i.e., the water flow originates from the neighboring pixels). (Flowing towards the center); these neighboring pixels that meet the conditions are the upstream water pixels, and then indexed. Each of their confluence contributions is accumulated individually.

[0058] Calculation method: For each 3×3 neighboring pixel, calculate its slope direction to the center pixel. The water flow direction points to the steepest descending direction. In the multi-flow direction algorithm, the water flow is proportionally distributed to multiple descending directions. ,in, For the first The slope of the upstream water pixel The slope represents the gradient of the neighboring pixels flowing towards the center pixel.

[0059] Initial condition: TCA of watershed boundary pixels = 1 This embodiment also sets a Factor Credibility Index (FCI). When the core topographic factor and auxiliary factor meet the credibility index, the process proceeds to step S4. Otherwise, return to step S2; The formula for calculating the factor credibility index is as follows: ,in For each factor ( The coefficient of variation of ΔH,Slope,Pc,TWI,TCA). This is an average value function. The FCI value ranges from 0 to 1; a higher FCI indicates more stable and reliable terrain factor calculations. The triggering conditions are as follows: When FCI ≥ 0.85, the terrain factor calculation is deemed stable and reliable, and the process proceeds directly to step 4. When 0.7 ≤ FCI < 0.85, adjust the neighborhood window parameter (increase the k value or window size) and recalculate the factor, then re-evaluate the FCI. When FCI < 0.7, the terrain factor calculation is determined to be unstable, and the process returns to step 2 to optimize the DEM fusion parameters (such as increasing the range parameter of the kriging interpolation).

[0060] S4, based on Calculate the adaptive relative elevation threshold and utilize A decision tree for hilly, flat, and depression terrain was constructed using core and auxiliary topographic factors to classify the terrain of the target county. The relevant operations in step S4 will be introduced with specific examples: Adaptive relative elevation threshold ,in, The basic elevation threshold for the county is 0.80m (0.80m for Dongfeng County, corresponding to the standard value for hilly areas). The landform type adjustment coefficient is determined based on the Roughness Index (Roughness Index = D1 standard deviation / D1 mean × 100%) for each township. This is the error quality adjustment factor. .

[0061] Landform type adjustment coefficient The range of values ​​for is shown in Table 2: Table 2

[0062] Design principles of adaptive threshold 1. The higher the topographic complexity, the more lenient the threshold: In areas with complex topography (such as Anshu Town, with a complexity of 8.56%), the surface has large natural undulations. If a strict threshold is used, a large number of natural undulations will be misclassified as hills. Therefore, the threshold needs to be relaxed to enhance the fault tolerance.

[0063] 2. The lower the EQI, the more lenient the threshold: In areas with large DEM errors (such as Liaoheyuan Town, EQI=0.00), the inaccuracy of the DEM itself will be transmitted to the ΔH calculation. Using a lenient threshold can avoid misjudging the DEM error as terrain features.

[0064] 3. The depression threshold is fixed at -0.5m: Depression identification focuses on the clear topographic feature of "below the neighborhood". Actual measurements show that plots that are more than 0.5m below the neighborhood have a clear tendency to accumulate water at the farmland scale. This threshold has good stability in different landform types and therefore does not change with townships.

[0065] like Figure 1 As shown, this embodiment adopts a three-layer decision-making architecture: priority determination of core factors, cross-validation of auxiliary factors, and conflict arbitration as a fallback, which includes a total of 7 decision nodes. The root node of the decision tree is the slope, because the slope is the most direct and stable terrain indicator to distinguish between hilly and non-hilly areas, and when there is a certain error in the DEM, the third-order unweighted difference algorithm has the best smoothing effect on the slope calculation.

[0066] Detailed explanation of the order of decision tree nodes (7 nodes): Node 1 (root node): Determines whether the slope is less than or equal to 2°. Distinguishes between "gentle terrain" and "steep terrain".

[0067] The judgment rules are as follows: When the slope ≤ 2°, the terrain is determined to be flat, possibly flat or low-lying, and the plot enters the left branch (node ​​2). When the slope is greater than 2°, the terrain is determined to be steep, possibly a hill or transition zone, and the plot enters the right branch (node ​​5). The reason for choosing slope as the root node is: 1. Slope is the most direct indicator for distinguishing between hilly terrain (Slope > 2°) and non-hilly terrain (Slope ≤ 2°). 2. The slope calculation uses a third-order unweighted difference algorithm, which has the best smoothing effect when there is a certain error in the DEM. 3. Slope is less sensitive to DEM errors than elevation difference, resulting in stronger calculation stability.

[0068] Node 2 (Left Branch): Relative Elevation Determination ? Applicable conditions: Plots whose root node is determined to be "yes" (Slope ≤ 2°); Purpose of the determination: In flat terrain, further distinguish between "hilly areas" and "non-hilly areas" by relative elevation.

[0069] The judgment rules are as follows: when When the relative elevation of a plot of land is significantly higher than that of its neighbors, it is classified as a hill (uplift type), and the result is output. when If the relative elevation of the plot is not significant, proceed to node 3. Note: To adapt the job location thresholds for each township, taking Dongfeng County as an example, each township... The range is 0.7050m (Nantunji Township) to 1.1233m (Anshu Town).

[0070] Node 3 (left branch): Depression determination ΔH < -0.5m? Applicable conditions: Plots where node 2 determines "No"; Purpose of determination: To identify depressions (micro-lowlands); Judgment rules: When ΔH < -0.5m, the relative elevation of the plot is determined to be significantly lower than that of the neighboring area, which may be a depression. The plot is then moved to node 4 for auxiliary factor verification. When ΔH ≥ -0.5m, the relative elevation of the plot is within the flat area, and the auxiliary factor verification is performed at node 4. The reason for fixing the depression threshold at -0.5m is that the key to depression identification is "below the neighborhood", and actual measurements show that plots that are more than 0.5m below the neighborhood have a clear tendency to accumulate water at the farmland scale. This threshold has good stability across different landform types.

[0071] Node 4 (left branch): Auxiliary factor verification Pc ≈ 0 and 4 ≤ TWI ≤ 8? Applicable conditions: Plots after node 3 (i.e., all flat plots where ΔH < ΔH_EQI); Objective: To improve classification reliability through cross-validation using plane curvature (Pc) and topographic moisture index (TWI).

[0072] Judgment rules: 1. Verification conditions on flat ground (must be met simultaneously): (1) |Pc| ≤ 0.5 (the curvature of the plane is close to 0 and the terrain has no obvious bumps or depressions); (2) 4 ≤ TWI ≤ 8 (the terrain humidity is at a moderate level, neither drought nor flood); If |Pc| ≤ 0.5 and 4 ≤ TWI ≤ 8, it is determined to be flat land, and the result is output.

[0073] 2. Depression verification conditions (must be met simultaneously): (1) Pc < 0 (concave slope); (2) TWI>8 (high humidity index, prone to water accumulation); (3) TCA>500 m² (large cumulative runoff volume, strong water collection capacity); If Pc < 0 and TWI > 8, it is determined to be a depression and the result is output.

[0074] If none of the above conditions are met, proceed to conflict arbitration (node ​​7). Node 5 (Right Branch): Relative Elevation Determination ? Applicable conditions: Plots where the root node is determined to be "No" (Slope>2°); Purpose of determination: To identify a position in steep terrain.

[0075] Judgment rules: (1) When At that time, the plot of land is both steep and higher than the surrounding area, meeting the conditions for a core hill site, and is therefore determined to be a hill site. The output result is then generated. (2) When At that time, although the plot was steep, the relative elevation was not significant, which may be a slope transition zone, entering node 6; Node 6 (Right Branch): Auxiliary factor verification Pc>0 and TWI<6? Applicable conditions: Plots where node 5 determines "No"; Purpose of the determination: To verify whether a site is a hilly area by using plane curvature and topographic humidity index. Judgment rules: 1. On-site verification conditions (must be met simultaneously): (1) Pc>0 (convex slope, topographic protrusion); (2) TWI < 6 (low humidity index, not easy to accumulate water, good drainage); If Pc>0 and TWI<6, it is determined to be a post (auxiliary verification type), and the result is output; If not satisfied, proceed to conflict arbitration (node ​​7).

[0076] Node 7 (Conflict Arbitration Node): Core Factor Priority Rule Applicable conditions: Plots whose auxiliary factor verification in node 4 or node 6 is "not satisfied"; The decision rules (executed in descending order of priority) are shown in Table 3: Table 3

[0077] Example of conflict resolution: Example 1: For a certain plot of land, the slope is 3° (>2°), and the ΔH is 0.6m ( Pc = -0.3 (<0) → Slope is called "hill", elevation is called "non-hill", curvature is called "concave" → According to P5 rule, Slope takes priority → hill; Example 2: A plot of land has a slope of 1.5° (≤2°), ΔH = -0.3m (> -0.5m), and TWI = 9.5 (>8). Slope indicates "flat land", elevation indicates "flat land", and humidity indicates "low-lying land". According to the P4 rule, the core factor takes priority. It is flat land (TWI anomalies are considered as local water accumulation).

[0078] In this embodiment, the principle of priority setting includes: Core factor priority principle: Slope and ΔH serve as the primary criteria for decision tree evaluation, while auxiliary factors (Pc, TWI, TCA) are only used for cross-validation and conflict arbitration. When the results of core factors and auxiliary factors conflict, the core factors take precedence.

[0079] Slope takes precedence over relative elevation: when Slope > 2° and In cases of conflict (i.e., the plot is relatively steep but the relative elevation is not significant), slope is used as the criterion, and the area is tended to be classified as a hill or transition zone. This is because slope is a direct reflection of terrain stability, while relative elevation is greatly influenced by the choice of neighboring benchmarks.

[0080] The "veto" power of auxiliary factors: Auxiliary factors do not have independent judgment power, but can play a "veto" role under specific conditions. For example, when Slope≤2° and ΔH is within the flat area, if TWI>12 (extremely high humidity), even if other factors all point to flat areas, conflict arbitration should be triggered to consider whether it is a hidden depression.

[0081] S5. Set the classification confidence level. When the classification result meets the classification confidence level, output the classification result. Otherwise, return to step S2.

[0082] The relevant operations in step S3 will be introduced with specific examples: This embodiment also sets classification confidence levels. The results used to assess the reliability of each land parcel classification are a key indicator for triggering reverse feedback.

[0083] The formula for calculating the factor credibility index is: ,in, , and These are the weighting coefficients. In this embodiment, , , ; For slope confidence level, ; (Fixed value, uniform across the county, based on farmland topography standards to distinguish between gentle / steep terrain); The confidence level is the relative elevation. ; (Adaptive value, varies by township, determined by "basic threshold 0.80m × (1 + landform coefficient + (1)"). The value is calculated as EQI)×0.25”, and the range is generally between 0.7050m and 1.1233m.

[0084] The average confidence level of the auxiliary factors. ,in, For plane curvature confidence, The confidence level for the topographic humidity index. The confidence level of the cumulative flow.

[0085] Based on the above formula, when CC≥0.9, the classification is considered highly reliable, and the classification result is directly output. When 0.8 ≤ CC < 0.9, the classification result is basically reliable, and reclassification is performed after fine-tuning the threshold; The fine-tuning threshold includes: Adjusting the adaptive threshold for job positions .

[0086] (1) And the proportion of positions in the workplace is greater than 40%.

[0087] (2) Furthermore, the proportion of positions held in the local area is less than 20%.

[0088] (3) ,and ,

[0089] (4) In other cases, two-way fine-tuning probes are conducted. Try each once The maximum fine-tuning range shall not exceed the original The maximum number of fine-tuning operations is 2, with a maximum range of ±10%. If CC is still <0.9 after two fine-tuning operations, reverse feedback is triggered (return to step 2). When CC < 0.8, the classification result is deemed unreliable, triggering a reverse feedback loop to return to step 2.

[0090] The software environment for this embodiment is: ArcGIS / GRASS GIS + Python (GDAL / NumPy / scipy / sklearn); The hardware environment for this embodiment is: a regular PC (with ≥8GB of memory); The parameters are set as follows: neighborhood window 100m, number of k nearest neighbors 8~15, step size of lag distance of variogram is 1 / 2 of the average point spacing, and threshold for hills / flatlands / depressions is adaptively adjusted according to the terrain complexity of the county.

[0091] Example 2 This embodiment further defines Embodiment 1. To verify the effectiveness of the classification scheme described in the embodiment, the following comparative experiment is conducted: Experiment 1: Using 1,776 measured samples from four villages in Dongfeng County (Shifeng Village, Lequn Village, Zhongxiang Village, and Jinshan Village) as the validation set, Gangpingwa was classified based on D1 (measured) and D2 (12.5m DEM), respectively. A confusion matrix was constructed, and the classification accuracy was statistically analyzed.

[0092] The test results are as follows Figure 2 As shown: Overall classification accuracy was 93.74%; recall rate for hilly areas was 98.8% (408 / 413), with a precision of 96.0%; recall rate for flat areas was 63.8% (30 / 47), with a precision of 71.4%; and recall rate for low-lying areas was 57.9% (11 / 19), with a precision of 91.7%.

[0093] Depend on Figure 2 It can be known that: The overall classification accuracy of 93.74% exceeds the technical requirement of 90%, verifying the effectiveness of the "decision tree + adaptive threshold" classification method. The hilly terrain was the best identified (98.8% recall and 96.0% precision), indicating that the hilly terrain features were the most identifiable in the 12.5m DEM, which is consistent with the significant terrain features of the hilly terrain, namely "large elevation differences and steep slopes". The recall rates for flat and low-lying areas were relatively low (63.8% and 57.9%, respectively), indicating that the resolution of the DEM is a major limitation in low-lying and flat areas. The classification accuracy can be significantly improved after DEM fusion correction.

[0094] Experiment 2: Calculate the RMSE and standard deviation (STD) of 12 townships, draw a horizontal bar chart (RMSE is the bar length, STD is the error bar), and mark the threshold of 2.5m for hilly areas.

[0095] The test results are as follows Figure 3 It can be seen that: Liaoheyuan Town has the highest RMSE (2.91m) and the largest STD; Yunding Town has the lowest RMSE (1.64m) and the smallest STD; RMSE and STD are positively correlated in all towns.

[0096] Depend on Figure 3 It can be known that: The positive correlation between RMSE and STD indicates that townships with large errors also have large error fluctuations, providing a dual-indicator verification for "error spatial heterogeneity". The significant differences in error bar lengths indicate that the degree of error dispersion varies across different townships, further supporting the necessity of regional adaptive correction. All townships had RMSE values ​​exceeding the threshold of 2.5m, verifying the overall conclusion that "12.5m DEM is insufficient in accuracy in Dongfeng County". Experiment 3: Using 104,265 calibration points as samples, a linear regression model D2 = a × D1 + b was established between D1 and D2. A scatter plot (5,000 sample points) was drawn and the regression line (red solid line) and 1:1 line (green dashed line) were overlaid, and the R² value was labeled.

[0097] The test results are as follows Figure 4 As shown: the regression equation is D2=1.008753×D1+(-3.2046); the Pearson correlation coefficient r=0.9957; the coefficient of determination R²=0.9915; the regression line basically coincides with the 1:1 line but there is a slight deviation.

[0098] Depend on Figure 4 It can be known that: R² = 0.9915 indicates that D1 and D2 have a very strong linear correlation, and the 12.5m DEM can explain 99.15% of the variation in measured elevation. The regression slope of 1.008753 deviates from 1.0, and the intercept of -3.2046 is not 0, indicating the existence of systematic proportional bias and offset bias, which need to be corrected through "systematic bias elimination". The slight deviation of the regression line from the 1:1 line is more pronounced at higher elevations, indicating that the higher the elevation, the more significant the proportional deviation of the DEM.

[0099] Experiment 4: Calculate the mean error ME (ME=mean(D2 - D1)) of the 12 townships, and draw a horizontal bar chart according to the ME value. Red indicates ME>0 (D2 is too high) and blue indicates ME<0 (D2 is too low).

[0100] The test results are as follows Figure 5 As shown: ME ranges from -1.47m (Sanhe Manchu and Korean Township) to +0.97m (Dongfeng County); ME < 0 in 7 townships (D2 is too low), and ME > 0 in 5 townships (D2 is too high); the maximum positive and negative deviation difference reaches 2.44m.

[0101] Depend on Figure 5 It can be known that: Significant differences in ME across townships (range 2.44m) quantitatively demonstrate the existence of a "systematic bias" component in DEM error, providing spatial evidence for the conclusion that "systematic bias variance accounts for 8.3%" in error type decomposition. The positive and negative distribution of ME shows no obvious pattern, indicating that the system bias is related to geographical location rather than terrain type, which further supports the correction principle of "spatial heterogeneity driven". The ME was highest in Dongfeng County (+0.97m) and lowest in Sanhe Manchu and Korean Township (-1.47m), providing specific offset parameters for township-level system deviation correction.

[0102] Experiment 5: Sort the absolute errors of 104,265 verification points, calculate the cumulative percentage, draw the cumulative distribution curve, and mark the RMSE (2.52m) position line.

[0103] The test results are as follows Figure 6 As shown: approximately 50% of samples have an absolute error <1.5m; approximately 70% of samples have an absolute error <2.5m (RMSE value); approximately 95% of samples have an absolute error <5m; the maximum absolute error is approximately 8m.

[0104] Depend on Figure 6 It can be known that: Approximately 50% of the samples had an absolute error of <1.5m, indicating that the DEM accuracy of half of the samples was acceptable. However, about 30% of the samples still had errors exceeding RMSE (2.52m), and about 5% of the samples had errors exceeding 5m, indicating that although the number of high-error samples was small, their impact was significant. The inflection point of the cumulative distribution curve is near RMSE, indicating that RMSE is an effective dividing point for distinguishing between "acceptable error" and "error requiring correction".

Claims

1. A method for evaluating the terrain of hilly and depression areas based on multi-scale DEM fusion verification, characterized in that, The method includes the following steps: S1. Obtain the DEM data of the target county with an accuracy of 0.01m, denoted as S1. DEM data with an accuracy of 12.5m is denoted as and will Preprocessing is performed to obtain the same as Datasets with completely consistent spatial extent ; S2, targeting Perform three-level correction, based on the results of the three-level correction and Perform weighted fusion to generate a fused DEM; S3, Calculate the error quality index And based on the fused DEM, the core and auxiliary factors of the terrain are calculated; Set the factor credibility index. When the core topographic factor and auxiliary factor meet the credibility index, proceed to step S4. Otherwise, return to step S2; S4, based on Calculate the adaptive relative elevation threshold and utilize A decision tree for hilly, flat, and depression terrain was constructed using core and auxiliary topographic factors to classify the terrain of the target county. S5. Set the classification confidence level. When the classification result meets the classification confidence level, output the classification result. Otherwise, return to step S2.

2. The method for evaluating the terrain of Gangpingwa based on multi-scale DEM fusion verification according to claim 1, characterized in that, The preprocessing includes: unifying resolution, cropping, and outlier removal.

3. The method for evaluating the terrain of Gangpingwa based on multi-scale DEM fusion verification according to claim 2, characterized in that, The three-level correction includes: Township Correction: ,in, The spatial location of the point to be corrected. This represents the average error of the township to which the current point to be corrected belongs; Village-level correction: , This represents the average error of the village to which the current point to be corrected belongs; Spatial correlation error correction: , This represents the spatial correlation error value estimated by ordinary Kriging interpolation for the current point to be corrected.

4. The method for evaluating the terrain of Gangpingwa based on multi-scale DEM fusion verification according to claim 3, characterized in that, The formula for calculating the merged DEM is: ,in, and These are the corresponding weighting coefficients.

5. The method for evaluating the terrain of Gangpingwa based on multi-scale DEM fusion verification according to claim 4, characterized in that, Error Quality Index The formula for calculation is: ,in, The RMSE value is the 3×3 neighborhood window of the current spatial location of the point to be corrected. and They are respectively The maximum and minimum values.

6. The method for evaluating the terrain of Gangpingwa based on multi-scale DEM fusion verification according to claim 5, characterized in that, The core terrain factors include: relative elevation. and slope ; The auxiliary factors include: plane curvature Pc, topographic humidity index TWI, and confluence accumulation TCA; The factor credibility index is set based on the coefficients of variation of the core terrain factor and auxiliary factors.

7. The method for evaluating hilly terrain based on multi-scale DEM fusion verification according to claim 6, characterized in that, Adaptive relative elevation threshold ,in, As the basic threshold for county-level elevation, Adjustment coefficient for landform type. This is the error quality adjustment factor. .

8. The method for evaluating hilly terrain based on multi-scale DEM fusion verification according to claim 7, characterized in that, The decision tree uses slope as the root node and divides the target plot into a gentle terrain branch and a non-gentle terrain branch based on the comparison result between the slope and the preset slope threshold. In the gentle terrain branch and the non-gentle terrain branch, the comparison results of relative elevation and adaptive relative elevation threshold are used as secondary judgment conditions, and the secondary judgment results are verified by combining at least one auxiliary factor among plane curvature, topographic humidity index and runoff accumulation, so as to determine the classification result of the target county.

9. A method for evaluating hilly terrain based on multi-scale DEM fusion verification according to claim 8, characterized in that, When the secondary determination result is inconsistent with the verification result obtained from the auxiliary factor, conflict arbitration shall be performed. The conflict arbitration determines the classification result of the target county according to the rule that core factors take precedence over auxiliary factors.

10. The method for evaluating hilly terrain based on multi-scale DEM fusion verification according to claim 9, characterized in that, The formula for calculating the factor credibility index is: ,in, , and These are the weighting coefficients. For slope confidence level, The confidence level is the relative elevation. This represents the average confidence level of the auxiliary factors.

Images